Opened 11 years ago

Closed 8 years ago

#1540 closed Bugs (fixed)

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

Reported by: John Maddock Owned by: Steven Watanabe
Milestone: Boost 1.36.0 Component: random
Version: Boost 1.34.1 Severity: Problem
Keywords: Cc: jeidsath@…, pbristow@…

Description

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:

Change History (5)

comment:1 Changed 11 years ago by 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.

comment:2 Changed 11 years ago by John Maddock

Cc: pbristow@… added

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.

comment:3 Changed 9 years ago by Steven Watanabe

Component: Nonerandom

comment:4 Changed 9 years ago by Steven Watanabe

Owner: changed from jsiek to Steven Watanabe
Status: newassigned

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

comment:5 Changed 8 years ago by Steven Watanabe

Resolution: fixed
Status: assignedclosed

Fixed in [63126].

Note: See TracTickets for help on using tickets.