C# Poisson Cumulative Distribution for large Lambdas

• Download source - 8.7 KB

Introduction

Poisson Distribution in C# with Ramanujan’s factorial approximation

Background

If you need the Poisson cumulative distribution function for very large Lambdas, this code is an alternative. It’s evolved from the code in this article.

Using the Code

The formula for the Poisson probability mass function is:

In C#, this can be calculated as:

public double Cdf(long k)
    var e = Math.Pow(Math.E, -_lambda);
    long i = 0;
    var sum = 0.0;
    while (i <= k)
    {
        sum += e * Math.Pow(_lambda, i) / Factorial(i);
        i++;
    }
    return sum;
}

The problem is that both the nominator and denominator become very large as k and _lambda (_lambda = \(\lambda )\) increases, and the program crashes. This can be solved by using logarithms. For simplicity, the natural logarithm is used: \(log_{e}=ln\). The following rules are used in the calculation:

• ln(x)=n    => eⁿ=x
• ln(e^y) = y
• \(e^{ln(x)} = x    => e^{ln(x)}=y    => x=y,    for    x>= 0\)
• ln(x/y) = ln(x) - ln(y)
• ln(x*y) = ln(x) + ln(y)

The calculation uses logarithms and Ramanujan’s factorial approximation which says:

Please read more about it here.

Setting Math.Pow(_lambda, i) / Factorial(i) = (A) = (\frac{\lambda ^{i}}{i!}) gives

By using Ramanujan’s factorial approximation, we get:

Using C# notation and employing the fact that \(e^{ln\left ( x \right )}=x\), we get:

A = Math.Pow(e, i*ln(ℷ) -i*ln(i) + i - ln(i*(1 + 4*i*(1 + 2*i)))/6 - ln(π)/2))

This expression is used in the C# calculation as can be seen in the code below:

var log6ThTail = Math.Log(i * (1 + 4 * i * (1 + 2 * i)))/6;
var lnN = i * Math.Log(_lambda) - (i * Math.Log(i) - i + log6ThTail + logPiDivTwo);
n = Math.Pow(Math.E, lnN - _lambda);

Here's the code:

using System;

namespace PoissonEvaluator
{
	public class PoissonDistribution
	{
		private readonly double _lambda;

		public PoissonDistribution(double lambda = 1.0)
		{
			_lambda = lambda;
		}

		public double Pmf(long k)
		{
			if (k > 170 || double.IsInfinity(Math.Pow(_lambda, k)))
			{
				var logLambda = k * Math.Log(_lambda) - 
                                        _lambda - (k * Math.Log(k) - k +
					Math.Log(k * (1 + 4 * k * (1 + 2 * k))) / 6 + 
                                        Math.Log(Math.PI) / 2);
				return Math.Pow(Math.E, logLambda);
			}
			return Math.Pow(Math.E, -_lambda) * Math.Pow(_lambda, k) / Factorial(k);
		}

		public double Cdf(long k)
		{
			long i = 0;
			var sum = 0.0;
			var infinityIsFound = false;
			var eLambda = Math.Pow(Math.E, -_lambda);
			var logPiDivTwo = Math.Log(Math.PI) / 2;
			while (i <= k)
			{
				double n;
				if (infinityIsFound)
				{
					var log6ThTail = Math.Log(i * (1 + 4 * i * (1 + 2 * i))) / 6;
					var lnN = i * Math.Log(_lambda) - (i * Math.Log(i) - i +
						log6ThTail + logPiDivTwo);
					n = Math.Pow(Math.E, lnN - _lambda);
				}
				else
				{
					if (i > 170 || double.IsInfinity(Math.Pow(_lambda, i)))
					{
						infinityIsFound = true;
						var log6ThTail = Math.Log
                                                    (i * (1 + 4 * i * (1 + 2 * i))) / 6;
						var lnN = i * Math.Log(_lambda) - 
                                                            (i * Math.Log(i) - i +
							log6ThTail + logPiDivTwo);
						n = Math.Pow(Math.E, lnN - _lambda);
					}
					else
					{
						n = eLambda * Math.Pow(_lambda, i) / Factorial(i);
					}
				}

				sum += n;
				i++;
			}
			return (sum > 1) ? 1.0 : sum;
		}


		public double Factorial(long k)
		{
			long count = k;
			double factorial = 1;
			while (count >= 1)
			{
				factorial = factorial * count;
				count--;
			}

			return factorial;
		}
	}