# Simple algorithm for generating Poisson distribution

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?

• what is nextDouble()? May 7, 2014 at 16:18
• @PA6OTA Good question, I added the meaning and a link to the API. Sorry for leaving it out initially May 7, 2014 at 16:22

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].

• doesn't help too much... May 7, 2014 at 16:55
• And yet this is exactly the reason why the algorithm works.
– Did
Dec 1, 2014 at 9:28

Suppose that we have $k$ independent and identically distributed exponential random variables $X_1, \dotsc, X_k$ with parameter $\lambda$. If we define a counting process $\{N(t)\}_{t \geq 0}$ such that $S_k := X_1 + \dotsb + X_k$ is the occurring time of the $k$-th event, then this is a Poisson process with rate $\lambda$. This means that $\mathbb{P}\{N(t) \leq k\} = \mathbb{P}\{S_k \geq t\}$, or more informally that you can describe the distribution of $N(t)$ through the distribution of $S_k$.

Now consider that if $X$ is an exponential random variable with parameter $\lambda$, then $Y := F_X(X) = 1 - e^{-\lambda X}$ is a standard uniform random variable (here $F_X$ is the cumulative distribution function of $X$). Note also that $$F_X^{-1}(Y) = -\frac{1}{\lambda}\log(1 - Y)$$ is an exponential random variable with parameter $\lambda$. Since also $U := 1 - Y$ is a standard uniform random variable, the random variable $$X' = -\frac{1}{\lambda}\log U$$ is again an exponential random variable with parameter $\lambda$.

We want to count (at most) $k$ events in a given time interval $[0, t]$ (that we can suppose without loss of generality to be $[0, 1]$, i.e. we arbitrarily fix $t = 1$), that is we want $N(1) \leq k$. By what we observed at the beginning, this corresponds to the condition $S_k \geq 1$, which, in its turn, corresponds to the following: $$\sum_{i = 1}^k X_i \geq 1 \iff \sum_{i = 1}^k -\frac{1}{\lambda}\log U_i \geq 1 \iff -\frac{1}{\lambda}\log\bigg( \prod_{i = 1}^k U_i \bigg) \geq 1 \iff \prod_{i = 1}^k U_i \leq e^{-\lambda},$$ which is precisely what is described in Knuth's algorithm.

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.

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