Description
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:
x^a exp(-x) / gamma(a) = exp(a log(x) - x - loggamma(a)
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
Attachments
Issue Links
- is related to
-
NUMBERS-171 Update the InverseErfc function to support the full range of x in [0, 2]
- Closed
-
NUMBERS-172 Update the Erf function to use the Boost rational function approximtion
- Closed
-
STATISTICS-39 Correct the switch of function evaluation in the GammaDistribution
- Closed
-
STATISTICS-38 Optimise use of the Gamma function in the PoissonDistribution
- Closed
- relates to
-
NUMBERS-175 Add continued fraction implementations using a generator of terms
- Closed
-
NUMBERS-167 RegularizedGamma.P with precomputed LogGamma value
- Open
- links to