[
https://issues.apache.org/jira/browse/NUMBERS-174?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17443993#comment-17443993
]
Alex Herbert commented on NUMBERS-174:
--------------------------------------
{quote}but I would not define any other such functions (i.e. Lower.value(...)
and Upper.value(...)
{quote}
So are you suggesting not exposing the computation of the incomplete gamma
function?
These functions do overflow with moderate arguments as the value of gamma(z)
used for normalisation grows more than exponentially. Matlab does not have the
standard functions. It has the incomplete gamma functions as the regularised
versions. Then it provides scaled versions. See [Matlab
gammainc|https://uk.mathworks.com/help/matlab/ref/gammainc.html].
{code:java}
p(a, z) scaled = p(a, z) * gamma(a+1) e^z / x^a
{code}
"This scaling cancels out the asymptotic behavior of the function near 0, which
avoids underflow with small arguments."
Leaving out the standard incomplete gamma functions would add only the
following to the API:
{code:java}
LogGamma.value(z, int[] sign);
RegularizedGamma.P.derivative(double a, double x);
RegularizedGamma.Q.derivative(double a, double x); /* Optional */
{code}
The derivative is all that is required for the Gamma function to be used in the
statistics project for the PDF and log PDF.
I am fine with adding just this to the API. Then looking at the Complementary
functions API. The usage of the Gamma functions in statistics requires the
regularized gamma function with fixed a (Gamma, Nakagami) or fixed z (Poisson).
Note that the RegularizedBeta class has a complement too which is not currently
in the API (see NUMBERS-169). There is also this old ticket from Math regarding
accuracy of the RegularizedBeta NUMBERS-52. So my next step was to look at
porting the Boost beta function. This could also make used of the Complementary
functions interface. The Boost beta functions have standard and regularised
functions and the derivative of the regularised function so porting them should
be able to use a similar process.
See:
[Boost incomplete
beta|https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_beta/ibeta_function.html]
[Boost beta
derivative|https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_beta/beta_derivative.html]
I also note that there is this ticket regarding the LanczosApproximation:
[NUMBERS-47]. Since the class is only public for use in statistics then it
could be marked as deprecated when statistics no longer depends on it and uses
the new gamma function to compute the derivative.
> Update Gamma functions using the Boost implementation
> -----------------------------------------------------
>
> Key: NUMBERS-174
> URL: https://issues.apache.org/jira/browse/NUMBERS-174
> Project: Commons Numbers
> Issue Type: Improvement
> Components: gamma
> Affects Versions: 1.0
> Reporter: Alex Herbert
> Assignee: Alex Herbert
> Priority: Major
> Fix For: 1.1
>
>
> The current regularised incomplete gamma functions for P and Q compute using
> a series representation for all input. This method is not robust to extreme
> arguments.
> The Gamma functions are used in Commons Statistics in the Gamma and Poisson
> distributions.
> Large values of the mean in the Poisson distribution have low
> precision for the CDF. The distribution also has a configurable
> threshold for convergence of the gamma function evaluation
> (STATISTICS-38). Ideally this implementation detail should be
> removed from the API.
> The Gamma distribution has a function switch based on a threshold to
> compute the PDF. This threshold results in incorrect values for
> certain parameters (STATISTICS-39).
> The function switch in the Gamma distribution is based on the
> documentation for the Boost gamma functions [1]. However it does not
> implement all the suggested optimisations detailed in the most recent
> Boost documentation. This includes avoiding using the converging
> series of the gamma function when the convergence is slow or unstable,
> i.e. addresses the need for a configurable convergence threshold in
> the Poisson distribution.
> One part of the evaluation of the incomplete gamma functions require
> accuracy in the leading term:
> {noformat}
> x^a exp(-x) / gamma(a) = exp(a log(x) - x - loggamma(a){noformat}
> When x and a are large then using logs to compute this leads to
> cancellation. Use of logs to compute this term is done in the current
> implementation in numbers for RegularizedGamma P and Q. It is the
> source of low precision for the CDF for the Poisson distribution when
> the mean is large.
> I propose to port the Boost gamma functions to numbers. This will be a
> process similar to the recent port of the error functions to the same
> package in numbers. These functions improved both accuracy and
> speed over the existing implementation. Once the gamma functions are
> ported a comparison can be made between the existing and new
> implementations.
> [1]
> [https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html]
--
This message was sent by Atlassian Jira
(v8.20.1#820001)
