如何近似泊松分布 $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));
    }

但我不是很满意。我在 "Numerical recipes" 上搜索了一种评估泊松分布的好方法，但我没有找到任何相关信息。有更好的方法吗？

编辑：明确地说，我正在寻找最精确的方法来近似大型事件的概率，给定泊松分布，而不计算笨拙的 (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 的限制，这要么是好的，要么不好，取决于您的具体需求。

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

Answer 2

因为泊松概率服从简单递归

P(k,lam) = lam/k * P(k-1,lam)

一种可能性是对多项式使用霍纳规则之类的东西。即：

Sum{k|C[k]*P(k,lam)} = exp(-lam)*(C[0]+(lam/1)*(C[1]+(lam/2)*(..))..)

或

P = c[cut]
for k=cut-1 .. 0
  P = P*lam/(k+1) + C[k]
P *= exp(-lam)

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

How to approximate a Poisson distribution $P(n,lambda)$

c

approximation

numerical-analysis