Generating Poisson Random Variables

**telk02** · 07-24-2012, 03:13 PM

Does anyone know how to generate Poisson random variables when the mean is above 745? Using the code below work up to a mean of 745, after that point it just loops infinitely. I think it is related to the machine epsilon of Excel but I don't know of a work around to make it operate above 745.

Public Function pos(Mean)
' Generate and return a random Poisson variable
' with mean Mean
Dim Value As Double
Dim Prod As Double
Value = Exp(-Mean)
Prod = 1
Do Until (Prod < Value)
pos = pos + 1
Prod = Prod * Rnd()
Loop
pos = pos - 1
End Function

---------- Post added at 02:13 PM ---------- Previous post was at 11:10 AM ----------

I'm still struggling with this...

**shg** · 07-24-2012, 03:47 PM

EXP(-745) is the limit for double-precision floating point.

I wonder (don't know) if it's true that pos(n)+pos(n) has the same distribution as pos(2n). If so, that's an easy(ish) way out.

A Poisson distribution is infinitely divisible; I think that means that the foregoing is true, but not sure.

**telk02** · 07-24-2012, 04:44 PM

Yes, the sum of two independent Poisson random variables is itself Poisson distribution. So I can just split up the lambda into multiple random variables when going above 745.

Thanks for the help.

**shg** · 07-24-2012, 04:56 PM

Knuth writes about another (and faster) method for large u by Ahrens and Dieter that starts with a gamma-distributed variate; it runs in O(log(u)) time instead of O(u) (the code you posted). But I don't have a routine for gamma variates ...