c - 泊松计算(erlang C)

Question

我之前发布过这个，用户告诉我把它发布到 codereview 上。我做到了，他们关闭了它……所以在这里又一次:(我删除了旧问题)

我有这些公式:

我需要 erlangC 公式的泊松公式:

我试图在 C 中重建公式:

double getPoisson(double m, double u, bool cumu)
{
    double ret = 0;
    if(!cumu)
    {
        ret = (exp(-u)*pow(u,m)) / (factorial(m));
    }
    else
    {
        double facto = 1;
        double ehu = exp(-u);
        for(int i = 0; i < m; i++)
        {
            ret = ret + (ehu * pow(u,i)) / facto;
            facto *= (i+1);
        }
     }
     return ret;
}

Erlang C 公式:

double getErlangC(double m, double u, double p)
{
    double numerator = getPoisson(m, u, false);
    double denominator = getPoisson(m, u, false) + (1-p) * getPoisson(m, u, true);
    return numerator/denominator;
}

主要问题是，getPoisson 中的m 参数值很大 (>170) 所以它要计算 >170!但它无法处理。我认为原始数据类型太小而无法正常工作，或者你怎么说？

顺便说一句:这是我用于第一个泊松的阶乘函数:

double factorial(double n)
{
    if(n >= 1)
        return n*factorial(n-1);
    else
        return 1;
}

一些示例:

输入:

double l = getErlangC(50, 48, 0.96);
printf("%g", l);

输出:

0.694456 (correct)

输入:

double l = getErlangC(100, 96, 0.96);
printf("%g", l);

输出:

0.5872811 (correct)

如果我对 getErlangC 的第一个参数 (m) 使用大于 170 的值，例如:

输入:

double l = getErlangC(500, 487, 0.974);
printf("%g", l);

输出:

naN (incorrect)

异常(exception):

0.45269

我的方法如何？有没有更好的方法来计算 Poisson 和 erlangC？

一些信息:Excel 具有 POISSON 函数，在 Excel 上它可以完美运行...是否有办法查看 EXCEL 用于 POISSON 的算法(代码)？

Answer 1

(pow(u, m)/factorial(m)) 可以表示为递归循环，每个元素显示为 u/n，其中每个 n 都是 m! 的一个元素。

double ratio(double u, int n)
{
    if(n > 0)
     {
        // Avoid the ratio overflow by calculating each ratio element
        double val;
        val = u/n;
        return val*ratio(u, n-1);
      }
    else
      {
         // Avoid division by 0 as power and factorial of 0 are 1
        return 1;
      }
}

注意，如果你想避免递归，你也可以把它作为一个循环来做

double ratio(double u, int n)
{
    int i;
    // Avoid the ratio overflow by calculating each ratio element
    // default the ratio to 1 for n == 0
    double val = 1;
    // calculate the next n-1 ratios and put them into the total
    for (i = 1; i<=n; i++)
      {
        // Put in the next element of the ratio 
        val *=  u/i;
      }
    // return the final value of the ratio
    return val;
}