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 . 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, 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`"] QuantilePlot[RandomReal[GammaDistribution[22,10], 1000], data, Frame -> True]
As a gentle reminder, do not expect good performance from this two direct algorithms.