Hi, I've spent quite some time on MATH-753 [1], and I think I now have a satisfactory solution. The problem was to overcome the overflows which arise when computing the density of the Gamma distribution for large values of the argument and/or the scale parameter.
As I initially feared, what was proposed in the JIRA ticket leads to high floating point errors. I adapted a method proposed in BOOST [2] 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! 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 getLanczos(double): returns the value of the Lanczos sum. If you do not like this option, I can copy/paste the Lanczos approximation in the GammaDistribution class. I'm adverse to the latter option, as it leads to code duplication. What do you think? Best regards, Sébastien [1] https://issues.apache.org/jira/browse/MATH-753 [2] 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) [3] http://en.wikipedia.org/wiki/Lanczos_approximation, formula below "[...] the sum is recast into the following form:" --------------------------------------------------------------------- To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org For additional commands, e-mail: dev-h...@commons.apache.org
