Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: genpoisson & lambda>80

From (Roberto G. Gutierrez, StataCorp.)
Subject   Re: st: genpoisson & lambda>80
Date   Mon, 08 Sep 2003 10:48:46 -0500

Benoit Dulong <> asks:

> I searched the Stata listeserver archive and I have not seen any reported
> problem for -- genpoisson --.

> If X = Poisson distribution with parameter lambda, mean(X) = lambda

> --genpoisson-- seems to work for lambda < 80.
> --genpoisson-- seems to have problems for lambda > 80 ?

> Is this a known problem ?


to which Nick Cox <> replied:

> The problem looks like one of numerical analysis.

> I suspect that the algorithm used is being stretched beyond breaking point
> by your larger means. Clearly it never generates large enough high values
> for you to get the means you want. I suspect that also arises from the
> difficulty of holding exp(-mean) accurately when mean is large.

> I note also that a Poisson with mean 200 has skewness 1/sqrt(200) =
> .07071068 and kurtosis 3 + 1/200 = 3.005, so in practice you could not
> distinguish a sample from such a Poisson and one from a (discretised)
> Gaussian with mean 200 and sd sqrt(200).

It was indeed a numeric problem associated with the algorithm being used to
generate the Poisson deviates.  I have updated -genpoisson- to use a different
(and more robust) algorithm, and this new algorithm seems to hold up better
for larger values of lambda.  I tested values of lambda as high as 10,000 and
got statistically valid results.  

Benoit can get the new version -genpoisson- by reinstalling/updating the
-gendist- package located on my user site.

   . net from

Note that the algorithm (as well as the old one) loops over the possible
values of the Poisson deviate, and so for large lambda this can take some
time.  Perhaps, at some future date, I will replace the offending ado-loop
with a plugin.

*   For searches and help try:

© Copyright 1996–2022 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index