[R] Generating routine for Poisson random numbers

[Paul Meagher]

> You can generate Poisson random numbers from a Poisson process like this:
>
> rfishy<-function(lambda){
>     t <- 0
>     i <- -1
>     while(t<=lambda){
>         t<-t-log(runif(1))
>         i<-i+1
>     }
>     return(i)
> }

This is a nice compact algorithm for generating Poisson random numbers.  It
appears to work.  I am used to seeing a Poisson counting process implemented
using a "num frames" parameter, a success probability per frame, and
counting the unform numbers in the 0 to 1 range that fall below the success
probability over "num frames".  It is interesting to see an implementation
that bypasses having to use a num_frames parameter, just the lambda value,
and relies instead on incrementing a t value on each iteration and counting
the number of times you can do so (while t <= lambda).

> The name of this generator is descriptive, not a pub. It is very slow for
> large lambda, and incorrect for extremely large lambda (and possibly for
> extremely small lambda). If you only wanted a fairly small number of
> random variates with, say, 1e-6<lambda<100, then it's not too bad.

One could impose a lambda range check so that you can only invoke the
function using a lamba range where the Poisson RNG is expected to be
reasonably accurate.  The range you are giving is probably the most commonly
used range where a Poisson random number generator might be used?  Brian
Ripley also mentioned that the counting process based implementation would
not work well for large lambdas.  Do you encounter such large lambdas in
practice?   Can't you always, in theory, avoid such large lambdas by
changing the size of the time interval you want to consider?

> But why would anyone *want* to code their own Poisson random number
generator, except perhaps as an interesting student exercise?

Yes this is meant as an interesting exercise for someone who wants to
understand how to implement probability distributions in an object oriented
way (I am writing an article introducing people to probability modelling).
I am looking for a compact algorithm that I can easily explain to people how
it works and which will be a good enough rpois() approximation in many
cases.  I don't want to be blown out of the water for suggesting such an
algorithm to represent a Poisson RNG so if you think it is inappropriate to
learn about what how a Poisson RNG works using the above described
generating process,  then I would be interested in your views.

Thank you for your thoughts on this matter.

Regards,
Paul Meagher


>
>
> -thomas
>
>