Re: [boost] [random] poisson_distribution problem and fix

From: Joel Eidsath (jeidsath_at_[hidden])
Date: 2007-12-26 13:33:37

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_at_[hidden]> 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