Abramowitz and Stegun approximation for cumulative normal distribution

Question

(Note: I know this looks like a programming question, but I'm OK with the programming part and just want to understand the mathematics.)

I found a bit of code to calculate the integral of the normal bell curve that I wrote 2 years ago, and I don't remember why or how it works. I just have a note that it's based on "Abramowitz & Stegun (1964)".

I want to verify that it's correct, and find out how accurate it is. I tried doing some internet searches but I didn't find anything that I could follow. (I haven't finished integral calculus yet.)

Basically I calculate standard deviation as $\mu$, define $p(\mu, x) = \frac{e^{-\frac{x^2}{2\mu^2}}}{\sqrt{2\pi\mu^2}} $

From there, the function $f(\mu, x)$, specifying the probability that a random sample is less than $x$, is calculated as:

$ t= \frac{1}{1+0.2316419 \frac{x}{\mu}} $

$ f(\mu, x) = 1 - p(\mu, x) (0.319381530t -0.356563782t^2 + 1.781477937t^3 -1.821255978t^4 + 1.330274429t^5) $

So, is this correct? Also, how accurate is it?

For reference, the original code is:

static double probability_density(double stddev, double x) {
    return exp(-x*x/(2*stddev*stddev)) / (sqrt(2.0*Pi*stddev*stddev));
}

static double normal_prob_lower(double stddev, double x) {
    //Cumulative Distribution Function
    //find probability that, in a normal distribution, a random value is lower than x
    //Abramowitz & Stegun (1964) approximation
    x/=stddev;
    stddev=1;
    bool rev=0;
    if (x<0) {
        x=-x;
        rev=1;
    }
    const double b[]= {0.2316419,0.319381530,-0.356563782,1.781477937,-1.821255978,1.330274429};
    double t=1.0/(1.0+b[0]*x);
    double sum=0;
    double tm=1;
    for (int i=1;i<=5;i++) {
        tm*=t;
        sum+=b[i]*tm;
    }
    double prob = 1-sum*probability_density(stddev, x);
    if (rev) {
        prob=1-prob;
    }
    return prob;
}

Answer 1

Yes, your function $f$ is correct, and the constants contain no typos. I pasted them in a own code and compared the function output to the CDF of the normal distribution from Maple.

You don't need to search the internet for the error term because it is given in Abramowitz/Stegun 26.2.17 and the value is $$\left|\epsilon\left(\frac{x}{\mu}\right)\right| < 7.5\times 10^{-8}$$

Here a picture of the actual absolute error $f(1,x) - P(x)\;$ confirming the given bound:

A few notes:

• Your notation is unusual because normally $\mu$ does not denote the standard deviation but the mean, $\sigma$ is mostly used for the std. dev. related parameter.

• Another often seen practice is to implement the standard normal CDF with $\mu=0, \sigma=1,$ and use this for the general Gaussian CDF.

• If you have access to a standard math library you will most probably find the error function erf and maybe the complementary function $\mathrm{erfc}(x) = 1-\mathrm{erf}(x)$ which nowadays are used to compute the normal CDF, i.e. the A&S function $P(x)$

$$P(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x} e^{-t^2/2} \mathrm{dt} = \int_{-\infty}^{x} Z(t) \mathrm{dt} = \frac{1}{2} \left( 1 + \mathrm{erf}\left( \frac{x}{\sqrt{2}}\right) \right) = \frac{1}{2} \mathrm{erfc}\left( -\,\frac{x}{\sqrt{2}}\right)$$