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

【问题标题】：Magic numbers in C++ implementation of Excel NORMDIST functionExcel NORMDIST 函数的 C++ 实现中的幻数
【发布时间】：2011-06-23 11:35:15
【问题描述】：

在寻找 Excel 的 NORMDIST 的 C++ 实现时（累积）我找到的函数this on a website:

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;
}

我有限的数学知识让我想到了Taylor series，但我无法确定这些数字的来源：

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

任何人都可以建议它们来自哪里，以及它们是如何得到的？

【问题讨论】：

嗯，3.14159265358979323846 非常接近 PI... ;)
关于 Math.StackEchange 的这个问题似乎是相关的：Abramowitz and Stegun approximation for cumulative normal distribution

标签： c++ excel function math magic-numbers

【解决方案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 之间，您可以通过Chebyshev approximation 一次性找到 P（数字食谱，第 5.8 节）。您不使用泰勒展开：您希望近似在整个实线上都很好，这是泰勒展开无法保证的。 Chebyshev 逼近是L^2 norm 中最好的多项式逼近，可以很好地替代很难找到的minimax polynomial（= sup norm 中的最佳多项式逼近）。

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

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

但过程类似，直接计算 normCdf，而不是 erf。

特别是非常相似，您使用的“实现”与文本中处理的“实现”有所不同，因为它的形式为b*exp(-a*z^2)*y(t)，但它也是 Chevishev 近似值。如您在 Schonfelder(1978)[http://www.ams.org/journals/mcom/1978-32-144/S0025-5718-1978-0494846-8/S0025-5718-1978-0494846-8.pdf] 的这篇论文中看到的 erfc(x) 函数

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

【讨论】：