# Approximation of a Normal Distribution function

I am reviewing and documenting a software application (part of a supply chain system) which implements an approximation of a Normal Distribution function; the original documentation mentions the same/similar formula quoted here

$$\phi(x) = {1\over \sqrt{2\pi}}\int_{-\infty}^x e^{-{1\over 2} x^2} \ dx$$

This is approximated with what looks like an asymptotic series like here however the expression is slightly different:

$$\phi(x) \simeq {1\over x\sqrt{2\pi}} \Bigl( 1 - {1\over x^2} + {3\over x^4} - {15\over x^6} + {105 \over x^8}\Bigr)$$

The original document quotes a "Cambridge Statistical Tables" book which I don't have; also, the software uses a precalculated table for values of x between [-8.2, 8,2] and uses the approximation function for values of x outside that interval.

I need to find some reference that explains why this particular formula was chosen, and how accurate the approximation really is. (The documentation claims that an error less then the 7th decimal is "not bad")

Also, the constant 1 ⁄ √ 2π does it have a name? All I was able to find was that "this expression ensures that the total area under the curve Φ(x) is equal to one"

Does it matter if the constant has only 10 decimals when the double type in Java (IEEE 754) has a precision of aproximate 16 decimals?

• In the first equation the exponent should have a minus. In the second one, you have an unclosed paren. The expansion is wrong. Fragments of it resemble the correct expansion for integral from x to $\infty$, but the series shouldn't be in a denominator, and you're missing an overall factor of $\exp(-x^2/2)$. – Mio Mar 26 '10 at 11:32
• See here mathworld.wolfram.com/Erf.html , there's a section about the asymptotic expansion and how it's obtained. Mind the factors of 1/2. – Mio Mar 26 '10 at 11:40
• Thanks Mio, I fixed the parts which were wrong. The Erf page describe asymptotic series where the x variables are nominators while the expression I found in my application has the x variables on the denominator; the original documentation mentions the formula as described above hence I am still puzzled. – Adrian Mar 26 '10 at 12:42
• I don't see how that asymptotic expansion makes any sense at all --- it is going to zero for large x while Phi(x) is going to one. PS You may want to look at Abramowitz & Stegun (Handbook of Mathematical Functions) which gives a lot of details on various functions, expansions etc. I don't have my copy here but I'll bet they have something useful. – Apollo Apr 12 '10 at 14:35
• Hope this doesn't sound too much like nitpicking but I'd rather see a different name than $x$ for the variable of integration. – Sylvain JULIEN Mar 6 '18 at 16:59

This is the (divergent) asymptotic development for $f(x)=\int_x^{\infty} e^{-{1\over 2}t^2} \ dt$ given by

$$f(x) \sim e^{-{x^2\over 2}}\ (\ {1\over x}\ + \ \Sigma_{k=1}^{\infty}\ {(-1)^k(2k-1)!\over 2^{k-1}(k-1)!}\ {1\over x^{2k+1}}\ )$$ or $$f(x) \sim e^{-{x^2\over 2}}\ (\ {1\over x}\ - {1\over x^3} + {3\over x^5} - {15\over x^7} + {105\over x^9} - {945\over x^{11}}...)$$

It is easily obtained from an integration by parts. The rest is given by an explicit integral. From its expression, one can check that f(x) is in fact squeezed between two consecutive sums of the series. As a result, we have the bound, for all x>0,

$$0 \geq f(x) - e^{-{x^2\over 2}}\ (\ {1\over x}\ - {1\over x^3} + {3\over x^5} - {15\over x^7} + {105\over x^9}) \geq - e^{-{x^2\over 2}}\ {945\over x^{11}}$$

Divergent series were standard tools at the beginning of the XXe century. I can't provide a reference, but the book of Hardy "Divergent series" may be a starting point for a bibliographic search.

The divergent asymptotics was chosen because it is "good" at infinity. When truncated to some power (say 9), it has the correct behavior when x goes to infinity, that is, it goes to zero, just like f(x). So it gives an interesting approximation of f for all sufficiently big value of x. This of course is not the case of a development obtained by truncating a converging series in positive powers of x.

The source code of R (http://cran.r-project.org/mirrors.html ) contains implementation (in C) of various distributions, and among them, of a normal one (see, for example, the file src/nmath/pnorm.c in the source tree of the latest version 2.10.1). The approximation they are using seems to be different (and more precise) than yours, but there are references in comments. Probably you may get a hint by looking there.

For the normal CDF:

$$\Phi(x) \approx \frac{1}{1+e^{-1.65451*x}}$$

worked very well for me.

• What is the downvote for? – Alt Aug 4 '17 at 1:50
• It certainly is not mine, but I am guessing the downvoter was frustrated to not know where this expression came from. – Igor Rivin Aug 9 '17 at 16:19
• I needed a fast way of computing the normal CDF, and the similarity in the shape of this function and the sigmoid function ($\sigma(t) = \frac{1}{1+e^{-t}}$) attracted my attentions. Then, I ran a computer simulation to tune the coefficient $a$ such that $\int_{x1}^{x2}\|\Phi(x)-\sigma(at) \|_2^2dx$ is minimized for some large magnitude negative $x_1$ and positive $x_2$. $a$ turned out to be $−1.65451$. – Alt Aug 9 '17 at 16:52
• The down voter was I (now reversed). Thanks for the explanation, which explains what your approximation is trying to optimize. Note that the OP asked about approximation for $x>8.2$, for which your approximation, while OK in terms of absolute difference between the actual value and the approximation, is bad in terms of relative difference. – Yoav Kallus Aug 9 '17 at 17:36