Photo Credit: WeGraphics

A Simple Poisson Random Number Generator

I recently need to build a simulation relied upon poisson process. Well, I came across Knuth’s algorithm which seems to be famous:

algorithm poisson random number (Knuth):
init:
Let L ← exp(−λ), k ← 0 and p ← 1.
do:
k ← k + 1.
Generate uniform random number u in [0,1] and let p ← p × u.
while p > L.
return k − 1.

Nice and beautiful.

The underlying logic of this algorithm employs a fundamental property of Poisson process: the time interval between two consecutive events is exponential distributed. Suppose the Poisson distribution we want to simulate is:

Then the corresponding exponential distribution for events time interval and its CDF would be

Inverse transforming approach is applied here: we generate uniform distributed variable u which falls in [0, 1] and transform it according to:

then a is a random variable follows exponential distribution, representing time interval. Replacing 1-u with u, and sum these time intervals will give us:

That is what the fifth line of Knuth’s algorithm is doing.

Finally we just need to count how many time intervals totally contained in a unit time, that is the Poisson number.

As some people noted, this algorithm’s time cost is linearly dependent with $\lambda$ and there is potential risk inherited in the term $\exp(−\lambda)$. Luckily, my simulation lives well with small $\lambda$.

I have a simple attached for you to play with. Just stroke gcc poisson.cpp -o poisson and then ./poisson 12 20; first parameter is $\lambda$ and second is the number of random numbers you want. To test using q-q plot, dump the result to Mathematica:

Needs["StatisticalPlots"]
QuantilePlot[RandomInteger[PoissonDistribution[22], 1000], data,
Frame -> True]

We can also extend this approach to simulate a gamma distribution function. Since $\alpha$ is an exponential distributed variable, the sum of $\alpha$ follows gamma distribution.

This interesting property facilitates us to write code , which just adds n exponential distributed variables with parameter theta together. Similar mathematica test justifies it.

Needs["StatisticalPlots"]
Frame -> True]