3

ExcelのNORMDIST（累積）関数のC ++実装を探しているときに、これをWebサイトで見つけました。

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

誰かがそれらがどこから来たのか、そしてそれらをどのように導き出すことができるのかを提案できますか？

4

1 に答える 1

10

数値レシピ、第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を一度だけ見つけることができます（数値レシピ、セクション5.8）。テイラー展開を使用しません。テイラー展開では保証できない、実数直線全体で近似が適切である必要があります。チェビシェフ近似は、L ^ 2ノルムでの最良の多項式近似であり、これは、見つけるのが非常に難しい最小多項式（= supノルムでの最良の多項式近似）の良い代替手段です。

ここのバージョンは少し異なります。代わりに、

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 3rd editionでは、6.2.2章の最後に、このタイプの非常に正確なC実装を提供します。t*exp(-z^2 + c0 + c1*t+ c2t^2 + c3*t^3 + ... + c9t^9)

于 2011-02-08T14:30:52.060 に答える