[jira] [Updated] (MATH-753) Better implementation for the gamma distribution density function

["Sébastien Brisard (JIRA)" <ji...@apache.org>]

     [ https://issues.apache.org/jira/browse/MATH-753?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]

Sébastien Brisard updated MATH-753:
-----------------------------------

    Attachment: lanczos.patch

{{lanczos.patch}} is the patch discussed in the the {{dev}} mailing-list:
{quote}
As I initially feared, what was proposed in the JIRA ticket leads to high floating point errors. I adapted a method proposed in [BOOST|http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_gamma/igamma.html] (formula (15))
with acceptable errors. Meanwhile, I've also managed to improve the accuracy of the computation of the density for the range of parameters where the previous implementation was already working: in this range, the accuracy _was_ about 300 ulps, and is now 1-2 ulps! {color:red}Note: I might have been too optimistic, here. There is a significant improvement, though{color}. I think this improvement is worth implementing.

The downside is that I need to expose the Lanczos implementation internally used by {{o.a.c.m3.special.Gamma}}. This approximation is so standard that I don't see it as a problem. I don't think that it reveals too much of the Gamma internals, since the javadoc of {{Gamma.logGamma}} states that it uses this approximation. So what I
propose is the addition of two methods in {{Gamma}}:
* {{double getLanczosG()}}: returns the {{g}} constant
* {{double lanczos(double)}}: returns the value of the Lanczos sum.
{quote}
                
> Better implementation for the gamma distribution density function
> -----------------------------------------------------------------
>
>                 Key: MATH-753
>                 URL: https://issues.apache.org/jira/browse/MATH-753
>             Project: Commons Math
>          Issue Type: Improvement
>    Affects Versions: 2.2, 3.0, 3.1
>            Reporter: Francesco Strino
>            Assignee: Sébastien Brisard
>            Priority: Minor
>              Labels: improvement, stability
>             Fix For: 3.1
>
>         Attachments: lanczos.patch
>
>
> 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.

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira