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

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.

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);
}

• Great! As I explained on my other answer, the application is actually Java (due to needing a cross-platform GUI) and I edited the code to work in C because I suspect C code is of more interest to scientific computation community than Java code. I verified the C++ variant and it seems to be internally using a method very similar to mine. I know Java has uncommons math which contains Poisson RNG, but that is horribly slow if the mean is million or more. Using existing libraries is always preferable. I tested that both C++ and gsl are fast in generating values with a large mean. – juhist Jul 6 '17 at 13:39

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.

• Why do you reimplement lgamma? It's in the C Standard Library. – Henri Menke Jul 5 '17 at 22:28
• Also drand48 is UN*X specific and your code is not portable. – Henri Menke Jul 5 '17 at 22:32
• Actually, I implemented this for a Java application but wanted my answer to be C code because C is more widely used for scientific computation than Java. (Ok, Fortran would be even more widely used than C but it differs so much from Java that porting the code would be non-trivial.) Java has no lgamma. Yes, drand48 is Unix-specific but then again, somebody copypasting the code will notice it doesn't compile and then replace drand48 with whatever RNG they are using. – juhist Jul 6 '17 at 13:35