Computing the Normal Distribution Function

Every once in a while, I need to evaluate the normal distribution function $\Phi(x)$:

$$ \Phi(x) = \frac{1}{\sqrt{2 \pi}} \int_{-\infty}^x \! e^{-\frac{1}{2}t^2} \, dt $$

Unfortunately, it is not always available in the standard math libraries, and hence I have to implement a “good-enough” version myself. Here are some options.

I usually do not need great precision — an error of $10^{-4}$ or so is usually quite acceptable. On the other hand, I need an implementation that is simple and robust, that will work on all platforms, and does not require care and feeding or particular skills and knowledge in numerical computation.

The following old chestnut is straight from Abramowitz and Stegun, entry 26.2.17. Its accuracy (absolute deviation from true value) is better than $7.5 \cdot 10^{-8}$ — that is good enough for some semi-serious scientific computations.

def f1(x):
    u = 1.0/(1.0 + 0.2316419*math.fabs(x))
    y = u*(0.319381530
           + u*(-0.356563782
                + u*(1.781477937
                     + u*(-1.821255978
                          + u*1.330274429))))
    
    z = math.exp(-0.5*x*x)/math.sqrt(2.0*math.pi)
    
    if x >= 0.0:
        return 1.0 - z*y
    else:
        return z*y

The following is a newer approximation. It is nice, because the formula is simple and requires only two magic constants. Its total deviation is less than $4.5 \cdot 10^{-4}$.

def f2(x):
    a = 0.647 - 0.021*math.fabs(x)
    if x >= 0.0:
        return 0.5*(1.0 + math.sqrt(1.0 - math.exp(-a*x*x)))
    else:
        return 0.5*(1.0 - math.sqrt(1.0 - math.exp(-a*x*x)))

Here is another, similarly straightforward formula, but with significantly better accuracy. The total deviation is less than $6.25 \cdot 10^{-5}$.

def f3(x):
    x = x/math.sqrt(2.0*math.pi)
    return 0.5*(1.0 + math.tanh(19.5*x - 55.5*math.atan(35.0*x/111.0)))

In fact, we can incorporate the factor $1/\sqrt{2 \pi}$ in the previous expression into the numerical coefficients. After some additional polishing, we obtain a variant of the previous function, with an accuracy of $3.5 \cdot 10^{-5}$! As a true one-liner, with only three numerical constants, each given to at most 4 digits, this version may well be the “go-to” function for most applications.

def f3a(x):
    return 0.5*(1.0 + math.tanh(7.7784*x - 55.49*math.atan(0.1258*x)))

The last one is from Abramowitz and Stegun again, entry 26.2.18. It is interesting because it is a strictly rational approximation and does not require any special functions, only basic numerical operations, so that it also works where there is no math library available at all. Its total deviation is less than $2.5 \cdot 10^{-4}$.

def f4(x):
    u = x if x >= 0.0 else -x

    y = 1 + u*(0.196854 + u*(0.115194 + u*(0.000344 + u*0.019527)))
    y *= y  # y^2
    y *= y  # y^4
    
    if x >= 0.0:
        return 1.0 - 0.5/y
    else:
        return 0.5/y  # 1 - 1.0-0.5/y

Here is a plot of the deviations from the “true” value for all functions, except the first. (The error of the first is so small as to be indistinguishable from zero on the scale of this graph.)

Resources

• The sample implementation: cdf.py

• Abramowitz & Stegun: Handbook of Mathematical Functions
  A Google search will turn up many further locations for this standard reference.

• The function f2(x) is based on the paper:
  A Simple Approximation for Normal Distribution Function by Medhat Edous, Omar Eidous (Mathematics and Statistics 6(4): 47-49, 2018)

• The functions f3(x) and f3a(x) are inspired by the paper (beware of potential typos):
  High Accurate Simple Approximation of Normal Distribution Integral by Hector Vazquez-Leal, Roberto Castaneda-Sheissa, Uriel Filobello-Nino, Arturo Sarmiento-Reyes, and Jesus Sanchez Orea1 (Mathematical Problems in Engineering: Volume 2012, Article ID 124029)

Instead of providing weblinks, I suggest doing a Google search on the titles.