Simple algorithm for generating Poisson distribution

Question

I found a very simple algorithm that draws values from a Poisson distribution from this project.

The algorithm's code in Java is:

public final int poisson(double a) {
        double limit = Math.exp(-a), prod = nextDouble();
        int n;
        for (n = 0; prod >= limit; n++)
            prod *= nextDouble();
        return n;
        }

nextDouble() is a function from the Random package in Java that returns a uniformly distributed random double, for example 0.885598042879084.

I can't understand how this creates a Poisson distribution.

Can someone explain?

Answer 1

It is related to the Poisson process: suppose $a$ is fixed and $N$ is the number of independent Exponential (mean 1) RV's to be added until the sum exceeds $a$. In this case, $N \sim Poisson(a)$.

In the code above, everything is anti-logged. For example, instead of adding Exponentials, they multiply Uniforms, due to the relation $Y = - \ln(U)$ is an Exponential (mean 1) RV whenever $U$ is Uniform[0,1].

Answer 2

From wikipedia:

The Poisson distribution is a discrete probability distribution that expresses the probability of a given number of events occurring in a fixed interval of time and/or space if these events occur with a known average rate and independently of the time since the last event.

I think that you could get something similar by adding together random values until you hit a maximum, as the random function has a "known average rate" and is basically "independent" of the last function call. It looks to me like multiplying it like so is just a transformation to simplify the code.

Answer 3

Yes, the previous answer is correct. This is simply a version of Knuth's algorithm. It gets slower the larger the parameter lambda.