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