Numerical method for approximating the standard Normal distribution cdf with mean 0 and variance 1

Question

The standard Normal distribution probability density function is $$p(x) = \frac{1}{\sqrt{2\pi}}e^{-x^2/2},\int_{-\infty}^{\infty}p(t)\,dt = 1$$ i.e., mean 0 and variance 1. The cumulative distribution function is given by the improper integral $$P(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-t^2/2}\,dt$$ Describe a numerical method for approximating $P(x)$ given a value of $x$ to a prescribed absolute error $\tau$. Your solution should be as efficient as possible. Justify your answer.

I believe we need to consider the case where $x > 0$ and $x < 0$ and then use one a numerical approximating method that is the most accurate to use. That is the crux of the problem, there are many methods, does one need to try each and everyone to see which one is best. Is there way of determining which one is the best given this particular problem before actually performing any computations. Any suggestions is greatly appreciated.

Answer 1

Say we look at $x<0$ only. Then consider a change of variable of the form $u=\frac{t}{1-t}$ to have$$P(x)=\int_{-1}^{\frac{x}{1-x}}\frac{e^{-\frac{u^2}{2 (u+1)^2}}}{\sqrt{2 \pi } (u+1)^2}du$$ which is manageable stuff since $$\lim_{u\to -1}\frac{e^{-\frac{u^2}{2 (u+1)^2}}}{\sqrt{2 \pi } (u+1)^2}=0.$$

Answer 2

First note that cumulative distribution of Normal distribution $$ P(x) = \frac{1}{2}+\frac{1}{2}\mathrm{erf}\left(\frac{x}{\sqrt{2}}\right)$$ where $\mathrm{erf}$ is the error function.

Abramowitz and Stegun gave several approximations of the error function, check this and that. I hope that it was helpful.

Answer 3

Before embarking on crafting a custom implementation, it seems advisable to check whether the CDF of the standard normal distribution is supported as a built-in function in the programming environment of your choice. For example, MATLAB offers a function normcdf, as does CUDA.

If no implementation of normcdf is available, the next thing to check is whether the programming environment offers an implementation of the complementary error function, i.e. $\mathrm{erfc}$. This exists as ERFC in Fortran 2008, and as erfc() in C99 and C++11, for example. The CDF of the standard normal distribution is related to $\mathrm{erfc}$ as follows:

$$P(x) = \frac{1}{2} \mathrm{erfc} (-x \sqrt{\frac{1}{2}})$$

The reason for using $\mathrm{erfc}$ instead of $\mathrm{erf}$ is to avoid subtractive cancellation that leads to inaccuracy in the tails. Note that for $x < 0$, there is error magnification, that is, the relative error incurred in scaling the argument to $\mathrm{erfc}$ is magnified. Computing $x\sqrt{\frac{1}{2}}$ should be performed as accurately as possible, and if the final result $P(x)$ is subject to tight ulp error bounds, additional compensation may be needed. A technique for computing such a product with maximum accuracy with the help of a fused-multiply add (FMA) operation is presented in this paper:

Nicolas Brisebarre and Jean-Michel Muller, "Correctly Rounded Multiplication by Arbitrary Precision Constants." IEEE Transactions on Computers, Vol. 57, No. 2, February 2008, pp. 165-174 (draft online)

The FMA operation, which computes $ab+c$ with a single rounding at the end is available as the standard math function fma() in C99 and C++11.

If a given computational environment does not provide a built-in way to compute $\mathrm{erfc}$, you might want to look into using or porting robust code for the computation of error functions by W. J. Cody which can be found on Netlib. That Fortran code is based on the following paper:

W. J. Cody, "Rational Chebyshev approximations for the error function." Mathematics of Computation, Vol 23., No. 107 (1969), pp. 631-637. (online)

If Cody's code is not suitable for your work (e.g. due to licensing issues), a custom implementation of $\mathrm{erfc}$ could be created based on the following paper:

M. M. Shepherd and J. G. Laframboise, "Chebyshev approximation of $(1 + 2x)\exp(x^2)\operatorname{erfc} x$ in $0 \leqslant x < \infty$." Mathematics of Computation, Vol. 36, No. 153 (Jan., 1981), pp. 249-253 (online)

I used the algorithm from this paper as the basis for the implementation of erfc() for a shipping parallel computation platform, but used a freshly generated polynomial minimax approximation rather than the Chebyshev approximation from the paper. Tools like Maple or Mathematica offer functionality for generating such minimax approximations, or you could generate your own using the Remez algorithm.