Магические числа в C ++ реализации функции Excel NORMDIST - PullRequest
Новая
       43

Магические числа в C ++ реализации функции Excel NORMDIST

3 голосов
/ 08 февраля 2011

При поиске C ++ реализации Excel NORMDIST (накопительный) функция, которую я нашел это на веб-сайте : 

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

Кто-нибудь может подсказать, откуда они берутся и как их можно получить?

1 Ответ

9 голосов
/ 08 февраля 2011

Ознакомьтесь с числовыми рецептами, глава 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))

для положительного х, где 

t = 2 / (2 + x)

и поскольку t находится между 0 и 1, вы можете найти P по чебышевскому приближению раз и навсегда (Числовые рецепты, раздел 5.8). Вы не используете расширение Тейлора: вы хотите, чтобы аппроксимация была хорошей во всей реальной линии, что расширение Тейлора не может гарантировать. Чебышевское приближение является наилучшим полиномиальным приближением в L ^ 2 норме , что является хорошей заменой очень сложному нахождению минимаксного полинома (= наилучшее полиномиальное приближение в верхней норме).

Версия здесь немного отличается. Вместо этого каждый пишет 

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

но процедура аналогична, и normCdf вычисляется напрямую, а не erf.

В частности и очень похожая «реализация», которую вы используете, несколько отличается от той, что обрабатывается в тексте, потому что она имеет вид b*exp(-a*z^2)*y(t), но это также чевишевский ок. к функции erfc (x), как вы можете видеть в этой статье Schonfelder (1978) [http://www.ams.org/journals/mcom/1978-32-144/S0025-5718-1978-0494846-8/S0025-5718-1978-0494846-8.pdf]

Также в «Числовых рецептах», 3-е издание, в конце главы 6.2.2 они обеспечивают очень точную реализацию на языке C t*exp(-z^2 + c0 + c1*t+ c2t^2 + c3*t^3 + ... + c9t^9)

...