15,920,053 members
Articles / General Programming / Algorithms
Tip/Trick

# C# Poisson Cumulative Distribution for large Lambdas

Rate me:
8 Apr 2020CPOL1 min read 25.5K   234   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}$

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
{

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

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

 First Prev Next
 Nice work man!! You saved us!! Member 1481063922-Apr-20 19:13 Member 14810639 22-Apr-20 19:13
 Large Lambdas Dan Buskirk9-Apr-20 8:28 Dan Buskirk 9-Apr-20 8:28
 Excellent! Chris Maunder28-Jun-18 4:54 Chris Maunder 28-Jun-18 4:54
 Enjoyed your article asiwel25-Nov-17 17:36 asiwel 25-Nov-17 17:36
 Re: Enjoyed your article Gunnar S11-Dec-17 3:27 Gunnar S 11-Dec-17 3:27
 Thanks. By using the tweaking described in the linked article the error reduction is apparent in the beginning, but then the zeroes start to disappear, until at 11! the error overtakes the decimal point. With the addition of the extra term described in the article this happens at 16! instead of at 11!. A small improvement in others words.modified 12-Dec-17 5:01am.
 Last Visit: 31-Dec-99 18:00     Last Update: 19-Jun-24 16:49 Refresh 1