The Poisson distribution describes
the probability of observing k events at a
given length of time if the events occur
independently at a constant rate λ.
It is often neccessary to simulate random variables
with a Poisson distribution, especially in physics
Monte Carlo simulations.
Sadly most
random number generators only give
uniformly distributed
random numbers, however using a uniform random number
pU~U(0,1)
one can calculate a Poisson distributed random
number
pPo~Po(λ).
The idea is simple, given λ, calculate the
cumulative probability S(n)
of observing n or less events for all
k>=0, and generate a uniform random number
pU.
The transformed random number is the first
n for which
pU
>=S(n).
This method is some times called the cumulant method and works for most probability distributions, but is most handy when calculating S(n) is easy.
Using the Poisson distribution function the sum can be
written as
S(n)=
Σk=0,n
e-λ
λn/k!
In an implementation this way of calculating
S(n) is slow, especially the
factorial and exponent parts of the expression.
To optimize one can use the iterative
form of the distribution:
P(n+1)=
P(n)λ/k
where
P(0)=
e-λ
It is also clever to put an upper limit on k
for the unlikely case that
pU=1
or else one might get stuck in an infinite loop.
If λ is constant for many calls it might be a good
idea to precalculate all S(n) and
store them in a lookup table.
Below is a
C
function that implements the idea.
Note that you have to provide your own
UniformRandomNumber().
(Code written by Inexpensive,
free to use as you like)
const int PoissonRandomNumber(const double lambda)
{
int k=0; //Counter
const int max_k = 1000; //k upper limit
double p = UniformRandomNumber(); //uniform random number
double P = exp(-lambda); //probability
double sum=P; //cumulant
if (sum>=p) return 0; //done allready
for (k=1; k<max_k; ++k) { //Loop over all k:s
P*=lambda/(double)k; //Calc next prob
sum+=P; //Increase cumulant
if (sum>=p) break; //Leave loop
}
return k; //return random number
}