Well, sure. But the [random] poisson_distribution function doesn't return the normal at high means. It returns 0 (slowly!). That's just wrong. It should either be fixed or throw an exception. And as I posted above, here's a fix: RealType product = RealType(0); for(result_type m = 0; ; ++m) { product += log(eng()); if(product <= -_mean) return m; } That's not NR code, that's just my having gotten rid of the exponentials in Knuth's algorithm by using logarithms. It actually simplifies things, as you can see. We can get rid of the init function and _exp_mean member variable in the code. Again, this still uses Knuth's algorithm, so it's tremendously slow. There is a rejection method algorithm on page 511 of the chapter_ten.pdf that you linked. Joel Eidsath On Dec 26, 2007 7:50 AM, Topher Cooper <topher@topherc.net> wrote:
Anyone attempting to generate non-uniform random numbers should consult "Non-Uniform Random Variate Generation" by Luc Devroye, which is available in full, for free use at:
http://cg.scs.carleton.ca/~luc/rnbookindex.html
Poisson generation is covered, specifically, in chapter 10:
http://cg.scs.carleton.ca/~luc/chapter_ten.pdf
starting on page 501.
How far off is the Poisson distribution from the normal at a mean of 750? How large a sample would you need to have, say, a one in a thousand chance of detecting the difference? (I don't know the answer, but if there is a problem at that extreme, its a question that should be answered before getting too fancy with "precise" solutions).
Topher