Generating random numbers with a Poisson distribution

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 p_U~U(0,1) one can calculate a Poisson distributed random number p_Po~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 p_U. The transformed random number is the first n for which p_U >=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^-λ λⁿ/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 p_U=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
}

The naive cumulative probability method, though the easiest method to think of, is not necessarily the simplest or fastest method, and the cumulative version of the coefficients has numerical stability issues related to repeated floating-point multiplication and division, which can lead to a distorted distribution. In fact, if generating random numbers is faster than floating point multiplication and division, the following method, developed by Knuth, will be faster:

const int PoissonRandomNumber(const double lambda)
{
  int k=0;
  const double target=exp(-lambda);
  double p=UniformRandomNumber();
  while (p>target)
  {
    p*=UniformRandomNumber();
    k+=1;
  }
  return k;
}

Although the code's operation is easy to understand, that it actually produces a Poisson distribution is somewhat harder. One explanation is to note that, if u is a uniformly-distributed random number, -ln(u) has an exponential distribution, which is the distribution that describes the time between events in a Poisson process; the Poisson distribution describes the number of events that occur in a specific amount of time.