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

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*^{n}=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:

$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:

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

I work as Senior Developer mainly in Microsoft environment and my strenghts are SQL, C#, Vue.js, Angular and ReactJS.