目录
介绍
用Ramanujan的阶乘近似值计算C#中的泊松分布
背景
如果需要用于很大型Lambda的泊松累积分布函数,则可以使用此代码。
使用代码
泊松概率质量函数的公式为:
在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;
}
问题在于,随着k和_lambda(_lambda = λ)增加,分母都变得非常大,并且程序崩溃。这可以通过使用对数来解决。为简单起见,自然对数是用于:。计算中使用以下规则:
计算使用对数和Ramanujan的阶乘近似表示:
请在此处详细了解。
设置Math.Pow(_lambda, i) / Factorial(i)= (A) = (\frac{\lambda ^{i}}{i!})给出了
通过使用Ramanujan的阶乘近似,我们得到:
使用C#表示法并利用,我们得到:
A = Math.Pow(e, i*ln(ℷ) -i*ln(i) + i - ln(i*(1 + 4*i*(1 + 2*i)))/6 - ln(π)/2))
可以在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);
这是代码:
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;
}
}