泊松計算（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或多精度算術。或者，拿出方程/算法的「模擬」。