泊松计算(erlang C)
Poisson calculation (erlang C)
我之前 post 编辑过这个,用户告诉我 post 它在 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)
例外:
0.45269
我的方法怎么样?有没有更好的方法来计算 Poisson 和 erlangC?
一些信息:Excel 具有 POISSON 函数,并且在 Excel 上它完美地工作...是否有一种方法可以查看算法(代码)EXCEL 用于泊松?
(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;
}
为了处理超出 double
范围的值,重新编码以使用 log 值。缺点 - 一些精度损失。
可以通过改进代码重新获得精度,但这里至少可以解决范围问题。
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
您的大型递归 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;
}
以上是 "okay",但由于 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
或多精度算术。或者,想出 equation/algorithm 的 "analog"。
我之前 post 编辑过这个,用户告诉我 post 它在 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)
例外:
0.45269
我的方法怎么样?有没有更好的方法来计算 Poisson 和 erlangC?
一些信息:Excel 具有 POISSON 函数,并且在 Excel 上它完美地工作...是否有一种方法可以查看算法(代码)EXCEL 用于泊松?
(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;
}
为了处理超出 double
范围的值,重新编码以使用 log 值。缺点 - 一些精度损失。
可以通过改进代码重新获得精度,但这里至少可以解决范围问题。
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
您的大型递归 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;
}
以上是 "okay",但由于 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
或多精度算术。或者,想出 equation/algorithm 的 "analog"。