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

Question

Хотя ищет C реализации ++ из NORMDIST (кумулятивного) функции в Excel Я нашел 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.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)) 

для положительных значений х, где

t = 2/(2 + x) 

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

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

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

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

В частности и очень аналогично «реализация», которую вы используете, несколько отличается от той, которая обрабатывается в тексте, потому что она имеет форму b*exp(-a*z^2)*y(t), но это также пример Чевишева. функции ERFC (х), как вы можете увидеть в этой статье Schönfelder (1978) [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)