Poisson distribution very slow for large mean (and may also overflow).

[John Maddock]

Joel Eidsath reports:

The poisson_distribution code will fail for any mean larger than 750 or so (on my x86_64 machine). The problem is a comparison with exp(-_mean) which evaluates as zero right around there.

The easiest fix is to operate with logs instead of exponents. I changed the operator() code to be something like this:

RealType product = RealType(0);

for(result_type m = 0; ; ++m) {

product += log(eng()); if(product <= -_mean)

return m;

}

This also makes it possible to get rid of the init() function and the _exp_mean private member variable.

A far better fix would be to use an algorithm other than Knuth's for this. Knuth's algorithm is simple, but O(n).

References to better implementation methods include:

• Numerical Recipes
• ​http://statistik.wu-wien.ac.at/unuran
• "Non-Uniform Random Variate Generation" by Luc Devroye. Full text available for free at ​http://cg.scs.carleton.ca/~luc/rnbookindex.html with Poisson covered in chapt 10 at ​http://cg.scs.carleton.ca/~luc/chapter_ten.pdf starting p501.

[John Maddock]

And another reference:

Ahrens, J.H. and Dieter, U. (1982). Computer generation of Poisson deviates from modified normal distributions. ACM Trans. Math. Software 8, 163-179.

[John Maddock]

Similar overflow issues occur in the Boost.Math Poisson distribution PDF function: now fixed in SVN Trunk by calling gamma_p_derivative to do the hard work.

[Steven Watanabe]

I'm working on new implementations of a couple distributions, but it's slow going.

[Steven Watanabe]

Fixed in [63126].