Opened 10 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:
- 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.
Attachments (0)
Change History (5)
comment:1 Changed 10 years ago by
comment:2 Changed 10 years ago by
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
Component: | None → random |
---|
comment:4 Changed 8 years ago by
Owner: | changed from jsiek to Steven Watanabe |
---|---|
Status: | new → assigned |
I'm working on new implementations of a couple distributions, but it's slow going.
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.