Algorithm to generate deviate for a Poisson distribution

342 Views Asked by At

According to a previous post, the following algorithm generates deviates of a Poisson distribution:

static Random random = new Random();
public static final int generatePoissonDeviate(double a) {
    double limit = Math.exp(-a);
    double product = random.nextDouble();
    int n;
    for( n = 0; product >= limit; n++ ) product *= random.nextDouble();
    return n;
}

What is "a" in this algorithm? It does not appear to be a mean.

When I run this algorithm it does not generate any zeroes.

I want to generate deviates of a Poisson distribution with a mean value N that includes zeroes. For example, assume number of children per family has a Poisson distribution and the mean number of children per family is 2.5. How can I generate deviates for individual families? Obviously the deviates must include the possibility of a family having zero children.

2

There are 2 best solutions below

2
On

Here is a link to a really efficient implementation of a routine taking a mean, using a flat random generator (on (0,1) and delivering a Poisson deviate:

http://www.fnal.gov/docs/working-groups/fpcltf/Pkg/CLHEP/Random/RandPoisson.cc

I wrote this some years ago. It is in C++, but you should be able to discern the algorithm from this code even if you don't want to use C++. Lhe lnGamma routine is based on Numerical Recipes.

0
On

The "a" in your algorithm is the mean of the Poisson distribution, you are simulating from: The output is the realization of a Poisson variate of mean $a$.

The algorithm should generate zeroes, and it looks fine when I run it. I don't have Eclipse on this machine, so here's the C# version, which generates plenty of zeroes.

    static Random random = new Random();
    public static int generatePoissonDeviate(double a)
    {
        double limit = Math.Exp(-a);
        double product = random.NextDouble();
        int n;
        for (n = 0; product >= limit; n++)
        {
            product *= random.NextDouble();
        }
        return n;
    }

    static void Main(string[] args)
    {
        for (int i = 0; i < 100; i++)
        {
            Console.WriteLine(generatePoissonDeviate(2.5));
        }
        Console.ReadKey();
    }