[MATH-753] Better implementation for the gamma distribution density function

• Type: Improvement
• Status: Closed
• Priority: Minor
• Resolution: Fixed
• Affects Version/s: 2.2, 3.0, 3.1
• Fix Version/s: 3.1
• Labels: improvement stability
• Assignee: Sebastien Brisard
• Reporter: Francesco Strino
• Created: 24/Feb/12 21:38
• Updated: 04/Mar/13 18:53
• Resolved: 16/May/12 05:39

Description

The way the density of the gamma distribution function is estimated can be improved.

It's much more stable to calculate first the log of the density and then exponentiate, otherwise the function returns NaN for high values of the parameters alpha and beta.

It would be sufficient to change the public double density(double x) function at line 204 in the file org.apache.commons.math.distribution.GammaDistributionImpl as follows:

return Math.exp(Math.log( x )*(alpha-1) - Math.log(beta)*alpha - x/beta - Gamma.logGamma(alpha));

In order to improve performance, log(beta) and Gamma.logGamma(alpha) could also be precomputed and stored during initialization.