15,921,295 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
{
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;
}
}

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.

## Comments and Discussions

 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
 It was a pleasure to read and to consider something I haven't thought about in a long while. (Fortunately, in statistics, there is often a nice denominator to cancel out most all of the terms in some larger factorials .. which very seldom, for me at least, get as big as you are discussing!) Running your program, I noticed (what I'd forgotten) about how big in fact the errors or differences actually become with such approximations. The Wiki article you reference talked entirely about Stirling's, not Ramanujan's approximation. But I found a more interesting discussion here: Is Ramanujan's approximation for the factorial optimal, or can it be tweaked? (answer below) - Mathematics Stack Exchange[^]
 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: 21-Jun-24 20:53 Refresh 1

General    News    Suggestion    Question    Bug    Answer    Joke    Praise    Rant    Admin

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