Excel NORMDIST函数的C++实现中的幻数

Question

在寻找Excel NORMDIST(累积)函数的C++实现时,我在网站上找到了这个:

static double normdist(double x, double mean, double standard_dev)
{
    double res;
    double x=(x - mean) / standard_dev;
    if (x == 0)
    {
        res=0.5;
    }
    else
    {
        double oor2pi = 1/(sqrt(double(2) * 3.14159265358979323846));
        double t = 1 / (double(1) + 0.2316419 * fabs(x));
        t *= oor2pi * exp(-0.5 * x * x) 
             * (0.31938153   + t 
             * (-0.356563782 + t
             * (1.781477937  + t 
             * (-1.821255978 + t * 1.330274429))));
        if (x >= 0)
        {
            res = double(1) - t;
        }
        else
        {
            res = t;
        }
    }
    return res;
}

我有限的数学知识让我想到泰勒系列,但我无法确定这些数字的来源:

0.2316419, 0.31938153, -0.356563782, 1.781477937, -1.821255978, 1.330274429

任何人都可以建议他们来自哪里,以及如何衍生出来？

Answer 1

查看数字配方,第6.2.2章.近似值是标准的.回想起那个

NormCdf(x) = 0.5 * (1 + erf(x / sqrt(2)))
erf(x) = 2 / (sqrt(pi)) integral(e^(-t^2) dt, t = 0..x)

写erf作为

1 - erf x ~= t * exp(-x^2 + P(t))

对于正x,在哪里

t = 2 / (2 + x)

由于t在0和1之间,你可以一次性找到P by Chebyshev近似(Numerical Recipes,5.8节).你没有使用泰勒展开:你希望近似在整个实线上是好的,泰勒展开不能保证.切比雪夫近似是L ^ 2范数中的最佳多项式近似,这是非常难以找到的极小极大多项式的良好替代(=支持范数中的最佳多项式近似).

这里的版本略有不同.相反,一个人写道

1 - erf x = t * exp(-x^2) * P(t)

但过程类似,并且直接计算normCdf,而不是erf.

特别地,非常类似地,您使用的"实现"与文本中处理的实现有所不同,因为它具有形式,b*exp(-a*z^2)*y(t)但它也是Chevishev约.你可以在Schonfelder(1978)的论文中看到erfc(x)函数[ http://www.ams.org/journals/mcom/1978-32-144/S0025-5718-1978-0494846-8/ S0025-5718-1978-0494846-8.pdf ]

同样在Numerical Recipes第3版中,在6.2.2章的最后,它们提供了非常精确的类型的C实现 t*exp(-z^2 + c0 + c1*t+ c2t^2 + c3*t^3 + ... + c9t^9)