C# code for generating Poisson random values

The following C# code will generate random values from a Poisson distribution. It depends on a uniform random number generator function GetUniform such as you can find here. It also requires a function LogFactorial that returns the natural logarithm of the factorial of its argument. More on LogFactorial below.

public static int GetPoisson(double lambda)
{
    return (lambda < 30.0) ? PoissonSmall(lambda) : PoissonLarge(lambda);
}

private static int PoissonSmall(double lambda)
{
    // Algorithm due to Donald Knuth, 1969.
    double p = 1.0, L = Math.Exp(-lambda);
    int k = 0;
    do
    {
        k++;
        p *= GetUniform();
    }
    while (p > L);
    return k - 1;
}

private static int PoissonLarge(double lambda)
{
    // "Rejection method PA" from "The Computer Generation of 
    // Poisson Random Variables" by A. C. Atkinson,
    // Journal of the Royal Statistical Society Series C 
    // (Applied Statistics) Vol. 28, No. 1. (1979)
    // The article is on pages 29-35. 
    // The algorithm given here is on page 32.

    double c = 0.767 - 3.36/lambda;
    double beta = Math.PI/Math.Sqrt(3.0*lambda);
    double alpha = beta*lambda;
    double k = Math.Log(c) - lambda - Math.Log(beta);

    for(;;)
    {
        double u = GetUniform();
        double x = (alpha - Math.Log((1.0 - u)/u))/beta;
        int n = (int) Math.Floor(x + 0.5);
        if (n < 0)
            continue;
        double v = GetUniform();
        double y = alpha - beta*x;
        double temp = 1.0 + Math.Exp(y);
        double lhs = y + Math.Log(v/(temp*temp));
        double rhs = k + n*Math.Log(lambda) - LogFactorial(n);
        if (lhs <= rhs)
            return n;
    }
}

The LogFactorial function takes some care to write well. It would be easy to write a naive implementation of factorial and take the logarithm of the output, but this would have two problems. First, the most obvious way to compute factorial has an execution time proportional to the size of the argument. Second, the factorial function can overflow for relatively small arguments.

If you have access to a function that returns the log of the gamma function, it would be best to use it. Otherwise, you could write a function that simply looks up the value of log( n! ) from an array for some large value of n and uses an asymptotic approximation for values of n beyond the range of the lookup table. Such an implementation is available here. For details, see How to compute log factorial.

This code is in the public domain. Do whatever you want with it, no strings attached.

Need help with randomization?