[
https://issues.apache.org/jira/browse/MATH-585?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13052035#comment-13052035
]
Mikkel Meyer Andersen commented on MATH-585:
--------------------------------------------
Thanks for the comments.
While working with the constants, I fell over a new and extremely interesting
paper [1] "A Simple Method for Generating Gamma Variables" from 2001 by
Marsaglia and Tsang (I'm tempted to ask: who else? :-)). It relies on a fast
Gaussian random generator (which I believe we have). It is equivalent to the
speed of the Gaussian random generator. Furthermore, it is really simple to
implement. So in short, it looks very interesting. Actually, it also seems
faster than Ahrens and Dieter (1982). A few samples (using JDKRandomGenerator):
{quote}
Generating 1000000 random gamma(1.0, 1.1) with mean 1.1 and var
1.2100000000000002
Marsaglia and Tsang (2001) took 269.0 ms. with mean 1.0997977261198317 and var
1.2113772336053434
Ahrens and Dieter (1982) took 358.0 ms. with mean 1.1003522016214415 and var
1.212513390495578
Generating 1000000 random gamma(1.5, 2.0) with mean 3.0 and var 6.0
Marsaglia and Tsang (2001) took 279.0 ms. with mean 2.998893072520475 and var
5.978837507498615
Ahrens and Dieter (1982) took 357.0 ms. with mean 2.9975071864627463 and var
5.999613004242815
Generating 1000000 random gamma(15.0, 1.0) with mean 15.0 and var 15.0
Marsaglia and Tsang (2001) took 259.0 ms. with mean 14.997425134380348 and var
15.006604469206364
Ahrens and Dieter (1982) took 278.0 ms. with mean 14.993230125478366 and var
14.984980114735736
{quote}
It is implemented as follows:
{code}
final double d = shape - 0.333333333333333333;
final double c = 1.0 / (3*FastMath.sqrt(d));
do {
final double x = generator.nextGaussian();
final double v = (1+c*x)*(1+c*x)*(1+c*x);
if (v <= 0) {
continue;
}
final double xx = x*x;
final double u = this.nextUniform(0, 1);
// Squeeze
if (u < 1 - 0.0331*xx*xx) {
return scale*d*v;
}
if (FastMath.log(u) < 0.5*xx + d*(1 - v + FastMath.log(v))) {
return scale*d*v;
}
} while (true);
{code}
The algorithm is also usable for shape < 1, but Ahrens and Dieter (1974) seems
better in that case.
What do you think about using Marsaglia and Tsang (2001) instead of Ahrens and
Dieter (1982)? I now prefer Marsaglia and Tsang (2001) :-).
Phil, in regards to freely available online reference, I agree that it would be
great, but I also think that it is second to choosing a good algorithm.
[1]: http://dx.doi.org/10.1145/358407.358414
> Very slow generation of gamma random variates
> ---------------------------------------------
>
> Key: MATH-585
> URL: https://issues.apache.org/jira/browse/MATH-585
> Project: Commons Math
> Issue Type: Improvement
> Affects Versions: 2.2, 3.0
> Environment: All
> Reporter: Darren Wilkinson
> Assignee: Mikkel Meyer Andersen
> Labels: Gamma, Random
> Attachments: MATH585-1.patch, MATH585-4-gamma.patch
>
> Original Estimate: 6h
> Remaining Estimate: 6h
>
> The current implementation of gamma random variate generation works, but uses
> an inversion method. This is well-known to be a bad idea. Usually a carefully
> constructed rejection procedure is used. To give an idea of the magnitude of
> the problem, the Gamma variate generation in Parallel COLT is roughly 50
> times faster than in Commons Math.
--
This message is automatically generated by JIRA.
For more information on JIRA, see: http://www.atlassian.com/software/jira
