How to generate Poisson-distributed random numbers quickly and accurately?

Question

I have attempted to create Poisson-distributed random numbers, seeing that it is not so easy as the simple multiplicative algorithm works accurately only if the mean is less than 500. Using logarithms is slightly slower, but it works at arbitrary mean values. However, it becomes extremely slow if the mean is very large.

My current algorithm is here in C-like language, assuming the functions nextUniformRandomDouble (returning uniform value between 0.0, inclusive, and 1.0, exclusive) and nextGaussianRandomDouble are defined:

int result = 0;
if (mean < 500.0)
{
    double L = exp(-mean);
    double p = 1;
    do
    {
        result++;
        p *= nextUniformRandomDouble();
    }
    while (p > L);
    result--;
}
else if (mean < 100000.0)
{
    double t = 0.0;
    while (t <= 1.0)
    {
        t += -log(1.0 - nextUniformRandomDouble())/mean;
        result++;
    }
    result--;
}
else
{
    result = (int)(nextGaussianRandomDouble() * sqrt(mean) + mean);
    if (result < 0)
    {
        // Probably never happens, but let's be safe
        result = 0;
    }
}

However, I am dissatisfied of this solution because of the low performance at values slightly below 100000 and the inaccuracy of the Gaussian approximation at values slightly above 100000.

Is there a better way of doing the same? What I'm looking for is extreme accuracy on one hand and extreme performance on the other hand.

Answer 1

If you can use C++11 you can use the built-in Poisson distribution

#include <iostream>
#include <random>

static std::random_device rd;
static std::mt19937 gen(rd());

int main()
{
  std::poisson_distribution<int> pd(5);

  for (int i = 0; i < 10000; ++i)
    std::cout << pd(gen) << '\n';
}

I mean, come on, it doesn't get any easier than this. And it's fast!

Here the output plotted with gnuplot

If you have to stick with plain C you could write a library which exposes the C++11 Mersenne Twister or use the GSL.

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int main()
{
  const gsl_rng_type * T;
  gsl_rng * r;

  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  for (int i = 0; i < 10000; ++i)
    printf("%d\n", gsl_ran_poisson(r, 5));

  gsl_rng_free(r);
}

Answer 2

I found the answer from the book Numerical Recipes in C, but unfortunately the code in that book is copyrighted. So, I re-derived a similar code that is slightly different but does the same thing.

The idea is that we select a random number from distribution $g(x)$ that has the property that $f(x) < c g(x)$ always with a finite constant $c$, where $f(x)$ is our Poisson distribution consisting of boxes of width 1.

Now, the selected $g(x)$ is $g(x) = \frac{1}{\pi \gamma}\frac{\gamma^2}{(x-x_0)^2+\gamma^2}$ which is the Cauchy distribution. The CDF is $y = F(x) = \frac{1}{\pi}\arctan(\frac{x-x_0}{\gamma}) + \frac{1}{2}$ and the inverse CDF is $x = x_0 + \gamma\tan(\pi(y-\frac{1}{2}))$. Thus, we generate $y$ randomly from uniform distribution between 0.0 and 1.0, and then calculate $x$. The parameters we select are $x_0 = \lambda$ and $\gamma = \sqrt{\lambda}$ where $\lambda$ is the parameter (mean value) of the Poisson distribution.

This distribution has negative values as well, so every time a negative value is obtained, the $y$ and $x$ need to be recalculated.

For this choice of distribution, it is true that $f(\lfloor x \rfloor)$ is always smaller than $2.4g(x)$ given any value of $\lambda$. So, we calculate the ratio $\frac{f(\lfloor x \rfloor)}{2.4g(x)}$ and reject if a uniformly distributed random variable $U$ between 0.0 and 1.0 is $U > \frac{f(\lfloor x \rfloor)}{2.4g(x)}$.

To calculate $f(\lfloor x \rfloor) = f(m) = \exp(-\lambda)\frac{\lambda^m}{m!}$, we calculate $\exp(m\log\lambda - \lambda - \textrm{lgamma}(m+1))$. The function $\textrm{lgamma}$ has a series representation.

Thus, the C code follows. I have tested that the C code gives me integers that have very accurate mean and also very accurate standard deviation. Thus, I believe the code is correct. The break-even point 60 was empirically determined so that the two algorithms given this point are equally fast.

The code is here:

double lgamma(double xx)
{
  double pi = 3.14159265358979;
  double xx2 = xx*xx;
  double xx3 = xx2*xx;
  double xx5 = xx3*xx2;
  double xx7 = xx5*xx2;
  double xx9 = xx7*xx2;
  double xx11 = xx9*xx2;
  return xx*log(xx) - xx - 0.5*log(xx/(2*pi)) +
         1/(12*xx) - 1/(360*xx3) + 1/(1260*xx5) - 1/(1680*xx7) +
         1/(1188*xx9) - 691/(360360*xx11);
}

int poissrnd_small(double mean)
{
  double L = exp(-mean);
  double p = 1;
  int result = 0;
  do {
    result++;
    p *= drand48();
  } while (p > L);
  result--;
  return result;
}

int poissrnd_large(double mean)
{
  double r;
  double x, m;
  double pi = 3.14159265358979;
  double sqrt_mean = sqrt(mean);
  double log_mean = log(mean);
  double g_x;
  double f_m;

  do {
    do {
      x = mean + sqrt_mean*tan(pi*(drand48()-1/2.0));
    } while (x < 0);
    g_x = sqrt_mean/(pi*((x-mean)*(x-mean) + mean));
    m = floor(x);
    f_m = exp(m*log_mean - mean - lgamma(m + 1));
    r = f_m / g_x / 2.4;
  } while (drand48() > r);
  return (int)m;
}

int poissrnd(double mean)
{
  if (mean < 60)
  {
    return poissrnd_small(mean);
  }
  else
  {
    return poissrnd_large(mean);
  }
}

Edit: I hereby place the code under public domain. Thus, you are free to use it however you want to, no need to attribute anything. The code has been entirely written by me.