15,919,132 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

## 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
 Congratulations for this piece of code. Right now is 2am in the morning and I'm doing a proyect for the university with a friend. The thing is that we couldn't figure out why our numbers turned so big when we put a lambda larger than 6. Thanks man!! Really appreciate it. - rodolfoRenomodified 23-Apr-20 1:41am.
 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
 Last Visit: 31-Dec-99 18:00     Last Update: 18-Jun-24 2:23 Refresh 1