c - 如何近似泊松分布 $P(n,lambda)$

Question

我必须实现一种算法来计算泊松函数的总和，每个泊松函数都带有乘法常数:

其中 C(k) 是正常数 <1，cut 是一个截止值，因为原则上总和应该取无限个 k，而 lambda 是一个在我的情况下可能从 20 到 100 变化的数字。我已经尝试过我的代码中的直接实现:

    #include<quadmath.h>
   ... //Definitions of lambda and c[k]...
    long double sum=0;
    for(int k=0;k<cutoff;++k)
    {
     sum=sum + c[k] powq(lambda,k)/tgamma(k+1)* (1.0/expq(lambda));
    } 

但是我不太满意。我在“数值食谱”上搜索了评估泊松分布的好方法，但我没有找到任何相关信息。有更好的方法吗？

<小时/>

编辑:要明确的是，我正在寻找最精确的方法来近似大型事件的概率，给定泊松分布，而无需计算尴尬的 (lambda^k)/k!因素!

Answer 1

嗯，一个简单的改进就是手动计算，并缓存 lambda^k 和 (k+1)!，因为它们在上一次迭代中的值可以用来快速计算当前迭代中的相应值， O(1) 计算。

另外，由于1.0/exp(lambda)是一个常量，所以你应该提前计算一次¹

#include<quadmath.h>
... //Definitions of lambda and c[k]...

const long double e_coeff   = 1.0 / expq(lambda);
long double inv_k_factorial = 1.0l;
long double lambda_pow_k    = 1.0l;
long double sum             = 0.0l;
for(int k=0; k < cutoff; ++k)
{
 lambda_pow_k *= lambda;
 inv_k_factorial /= (k + 1);
 sum += (c[k] * lambda_pow_k * inv_k_factorial);
} 
sum *= e_coeff;

现在这三个函数调用及其各自的开销已经完全从循环中消失了。

现在，我尝试使用与您在编写问题时相同的数据类型。由于您的评论表明 lambda 大于 1.0(我相信，由于 lambda_pow_k 快速减小，因此没有相对误差增长)这里丢失的任何意义都取决于 long double 的限制>，这是好还是坏，取决于您的具体需求。

<小时/>

1. 现在的编译器很聪明。所以它可以以任何方式进行优化，但我认为最好将不太明显的优化留给优化器。即使将代码交给非优化编译器，您的代码也不应降低性能。