Click here to Skip to main content
15,867,330 members
Articles / General Programming / Algorithms
Tip/Trick

C# Poisson Cumulative Distribution for large Lambdas

Rate me:
Please Sign up or sign in to vote.
5.00/5 (18 votes)
8 Apr 2020CPOL1 min read 24.8K   232   21   5
Implementation of the Poisson Cumulative Distribution function for large Lambdas
If you need the Poisson cumulative distribution function for very large Lambdas, this code is an alternative.

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:

$F \left ( k,\lambda \right ) = \sum_{i=0}^{k} \frac{e^{-\lambda }\lambda ^{i}}{i!}$

In C#, this can be calculated as:

C#
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    => en=x
  • ln(ey) = 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:

$ln \mathbf{\left (n!\right )}\approx n\mathbf{ln}\left (n\right )-n+\frac{ln\left ( n\left ( 1+4n\left ( 1+2n \right ) \right ) \right )}{6}+\frac{ln\left ( \pi \right )}{2}$

Please read more about it here.

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

$ln(A)=ln\left ( \frac{\lambda ^{i}}{i!} \right )=> ln(A)=ln(\lambda ^{i})-ln(i!)$

By using Ramanujan’s factorial approximation, we get:

$ln\left ( A \right )\approx \mathbf{i*ln\left ( \lambda \right )}-\mathbf{\left (i*ln\left ( i \right )-i+ln\left ( i*\left ( 1+4*i*\left ( 1+2*i \right ) \right ) \right )/6+\frac{ln\left ( \pi \right )}{2}\right )}$

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:

C#
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:

C#
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;
		}
	}

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


Written By
Software Developer (Senior)
Sweden Sweden
I work as Senior Developer mainly in Microsoft environment and my strenghts are SQL, C#, Vue.js, Angular and ReactJS.

Comments and Discussions

 
PraiseNice work man!! You saved us!! Pin
Member 1481063922-Apr-20 19:13
Member 1481063922-Apr-20 19:13 
QuestionLarge Lambdas Pin
Dan Buskirk9-Apr-20 8:28
Dan Buskirk9-Apr-20 8:28 
QuestionExcellent! Pin
Chris Maunder28-Jun-18 4:54
cofounderChris Maunder28-Jun-18 4:54 
PraiseEnjoyed your article Pin
asiwel25-Nov-17 17:36
professionalasiwel25-Nov-17 17:36 
GeneralRe: Enjoyed your article Pin
Gunnar S11-Dec-17 3:27
Gunnar S11-Dec-17 3:27 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.