泊松计算（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; 
} 

爱尔朗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！但它无法处理它。我认为原始数据类型太小而无法正确工作，或者你说什么？

顺便说一句：这是阶乘函数I用于第一个泊松：

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) 

如果我使用的值高于170，用于getErlangC的像的第一个参数（米）：

输入：

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

输出：

naN (incorrect) 

的例外：

0.45269 

怎么样我的做法？有没有更好的方法来计算泊松和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; 
} 

Answer 2

为了应付值超过double范围，重码使用日志值。下行 - 一些精确损失。

精度可以通过改进的代码重新获得，但这里至少可以解决范围问题。

OP代码的轻微变体如下：用于比较。

long double factorial(unsigned m) { 
    long double f = 1.0; 
    while (m > 0) { 
    f *= m; 
    m--; 
    } 
    return f; 
} 

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

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

建议修改

#ifdef M_PI 
    #define MY_PI M_PI 
#else 
    #define MY_PI 3.1415926535897932384626433832795 
#endif 

// log of n! 
// 
// Gosper Approximation of Stirling's Approximation 
// http://mathworld.wolfram.com/StirlingsApproximation.html 
// n! about= sqrt(pi*(2*n + 1/3.)) * pow(n,n) * exp(-n) 
static double ln_factorial(unsigned n) { 
    if (n <= 1) return 0.0; 
    double x = n; 
    return log(sqrt(MY_PI * (2 * x + 1/3.0))) + log(x) * x - x; 
} 


double getPoisson_2(unsigned m, double u, bool cumu) { 
    double ret = 0.0; 
    if (cumu) { 
    // Simplify term calculation. `mul` does not get too large nor small. 
    double mul = exp(-u); 
    for (unsigned i = 0; i < m; i++) { 
     ret += mul; 
     mul *= u/(i + 1); 
     // printf("ret:% 10e mul:% 10e\n", ret, mul); 
    } 
    } else { 
    // ret = (exp(-u) * pow(u, m))/(factorial(m)); 
    double ln_ret = -u + log(u) * m - ln_factorial(m); 
    return exp(ln_ret); 
    } 
    return ret; 
} 

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

测试代码

void ErTest(unsigned m, double u, double p, double expect) { 
    printf("m:%4u u:% 14e p:% 14e", m, u, p); 
    printf(" E0:% 14e", expect); 
    double y1 = getErlang(m, u, p); 
    printf(" E1:% 14e", y1); 
    double y2 = getErlang_2(m, u, p); 
    printf(" E2:% 14e", y2); 
    puts(""); 
} 

int main(void) { 
    ErTest(50, 48, 0.96, 0.694456); 
    ErTest(100, 96, 0.96, 0.5872811); 
    ErTest(500, 487, 0.974, 0.45269); 
} 

m: 50 u: 4.800000e+01 p: 9.600000e-01 E0: 6.944560e-01 E1: 6.944556e-01 E2: 6.944562e-01 
m: 100 u: 9.600000e+01 p: 9.600000e-01 E0: 5.872811e-01 E1: 5.872811e-01 E2: 5.872813e-01 
m: 500 u: 4.870000e+02 p: 9.740000e-01 E0: 4.526900e-01 E1:   nan E2: 4.464746e-01 

Answer 3

你大递归factorial是一个问题，因为它可能会产生一个堆栈溢出，以及一个值溢出。 pow也可能变大。

下面是逐步的东西结合的方式：以上

double 
getPoisson(double m, double u, bool cumu) 
{ 
    double sum = 0; 
    double facto = 1; 
    double u_i = 1; 
    double ehu = exp(-u); 
    double cur = ehu; 

    // u_i -- pow(u,i) 
    // cur -- current/last term in series 
    // sum -- sum of terms 

    for (int i = 0; i < m; i++) { 
     cur = (ehu * u_i)/facto; 
     sum += cur; 

     u_i *= u; 
     facto *= (i + 1); 
    } 

    return cumu ? sum : cur; 
} 

是“还行”，但仍可能会溢出，因为u_i和facto方面的一些价值观。

这里是一个将术语作为比率组合的替代方案。这是不太可能溢出：

double 
getPoisson(double m, double u, bool cumu) 
{ 
    double sum = 0; 
    double ehu = exp(-u); 
    double cur = ehu; 
    double ratio = 1; 

    // cur -- current/last term in series 
    // sum -- sum of terms 
    // ratio -- u^i/factorial(i) 

    for (int i = 0; i < m; i++) { 
     cur = ehu * ratio; 
     sum += cur; 

     ratio *= u; 
     ratio /= (i + 1); 
    } 

    return cumu ? sum : cur; 
} 

上述威力仍然产生一些大的值。如果是这样，您可能必须使用long double，quadmath或多精度算术。或者，拿出方程/算法的“模拟”。