A Simple Poisson Random Number Generator

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):
	         Let L ← exp(−λ), k ← 0 and p ← 1.
	         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 . 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:

	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.

	QuantilePlot[RandomReal[GammaDistribution[22,10], 1000], data, 
	 Frame -> True]


As a gentle reminder, do not expect good performance from this two direct algorithms.