Evaluate definite interval of normal distribution

Question

I know that an easy to handle formula for the CDF of a normal distribution is somewhat missing, due to the complicated error function in it.

However, I wonder if there is a a nice formula for $N(c_{-} \leq x < c_{+}| \mu, \sigma^2)$. Or what the "state of the art" approximation for this problem might be.

Answer 1

It depends on exactly what you are looking for. Below are some brief details and references.

Much of the literature for approximations centers around the function $$ Q(x) = \int_x^\infty \frac{1}{\sqrt{2\pi}} e^{-\frac{u^2}{2}} \, \mathrm{d}u $$

for $x > 0$. This is because the function you provided can be decomposed as a simple difference of the function above (possibly adjusted by a constant). This function is referred to by many names, including "upper-tail of the normal distribution", "right normal integral", and "Gaussian $Q$-function", to name a few. You'll also see approximations to Mills' ratio, which is $$ R(x) = \frac{Q(x)}{\varphi(x)} $$ where $\varphi(x) = (2\pi)^{-1/2} e^{-x^2 / 2}$ is the Gaussian pdf.

Here I list some references for various purposes that you might be interested in.

Computational

The de-facto standard for computing the $Q$-function or the related complementary error function is

W. J. Cody, Rational Chebyshev Approximations for the Error Function, Math. Comp., 1969, pp. 631--637.

Every (self-respecting) implementation uses this paper. (MATLAB, R, etc.)

"Simple" Approximations

Abramowitz and Stegun have one based on a polynomial expansion of a transformation of the input. Some people use it as a "high-precision" approximation. I don't like it for that purpose since it behaves badly around zero. For example, their approximation does not yield $\hat{Q}(0) = 1/2$, which I think is a big no-no. Sometimes bad things happen because of this.

Borjesson and Sundberg give a simple approximation which works pretty well for most applications where one only requires a few digits of precision. The absolute relative error is never worse than 1%, which is quite good considering its simplicity. The basic approximation is $$ \hat{Q}(x) = \frac{1}{(1-a) x + a \sqrt{x^2 + b}} \varphi(x) $$ and their preferred choices of the constants are $a = 0.339$ and $b = 5.51$. That reference is

P. O. Borjesson and C. E. Sundberg. Simple approximations of the error function Q(x) for communications applications. IEEE Trans. Commun., COM-27(3):639–643, March 1979.

Here is a plot of its absolute relative error.

The electrical-engineering literature is awash with various such approximations and seem to take an overly intense interest in them. Many of them are poor though or expand to very strange and convoluted expressions.

You might also look at

W. Bryc. A uniform approximation to the right normal integral. Applied Mathematics and Computation, 127(2-3):365–374, April 2002.

Laplace's continued fraction

Laplace has a beautiful continued fraction which yields successive upper and lower bounds for every value of $x > 0$. It is, in terms of Mills' ratio,

$$ R(x) = \frac{1}{x+}\frac{1}{x+}\frac{2}{x+}\frac{3}{x+}\cdots , $$

where the notation I've used is fairly standard for a continued fraction, i.e., $1/(x+1/(x+2/(x+3/(x+\cdots))))$. This expression doesn't converge very fast for small $x$, though, and it diverges at $x = 0$.

This continued fraction actually yields many of the "simple" bounds on $Q(x)$ that were "rediscovered" in the mid-to-late 1900s. It's easy to see that for a continued fraction in "standard" form (i.e., composed of positive integer coefficients), truncating the fraction at odd (even) terms gives an upper (lower) bound.

Hence, Laplace tells us immediately that $$ \frac{x}{x^2 + 1} < R(x) < \frac{1}{x} \>, $$ both of which are bounds that were "rediscovered" in the mid-1900's. In terms of the $Q$-function, this is equivalent to $$ \frac{x}{x^2 + 1} \varphi(x) < Q(x) < \frac{1}{x} \varphi(x) . $$ An alternative proof of this using simple integration by parts can be found in S. Resnick, Adventures in Stochastic Processes, Birkhauser, 1992, in Chapter 6 (Brownian motion). The absolute relative error of these bounds is no worse than $x^{-2}$, as shown in this related answer.

Notice, in particular, that the inequalities above immediately imply that $Q(x) \sim \varphi(x)/x$. This fact can be established using L'Hopital's rule as well. This also helps explain the choice of the functional form of the Borjesson-Sundberg approximation. Any choice of $a \in [0,1]$ maintains the asymptotic equivalence as $x \to \infty$. The parameter $b$ serves as a "continuity correction" near zero.

Here is a plot of the $Q$-function and the two Laplace bounds.

C-I. C. Lee has a paper from the early 1990's that does a "correction" for small values of $x$. See

C-I. C. Lee. On Laplace continued fraction for the normal integral. Ann. Inst. Statist. Math., 44(1):107–120, March 1992.

---

Durrett's Probability: Theory and Examples provides the classical upper and lower bounds on $Q(x)$ on pages 6–7 of the 3rd edition. They're meant for larger values of $x$ (say, $x > 3$) and are asymptotically tight.

Hopefully this will get you started. If you have a more specific interest, I might be able to point you somewhere.

Answer 2

I suppose I'm too late the hero, but I wanted to comment on cardinal's post, and this comment became too big for its intended box.

For this answer, I'm assuming $x >0$; appropriate reflection formulae can be used for negative $x$.

I'm more used to dealing with the error function $\mathrm{erf}(x)$ myself, but I'll try to recast what I know in terms of Mills's ratio $R(x)$ (as defined in cardinal's answer).

There are in fact alternative ways for computing the (complementary) error function apart from using Chebyshev approximations. Since the use of a Chebyshev approximation requires the storage of not a few coefficients, these methods might have an edge if array structures are a bit costly in your computing environment (you could inline the coefficients, but the resulting code would probably look like a baroque mess).

For "small" $|x|$, Abramowitz and Stegun give a nicely behaved series (at least better behaved than the usual Maclaurin series):

$$R(x)=\sqrt{\frac{\pi}{2}}\exp\left(\frac{x^2}{2}\right)-x\sum_{j=0}^\infty\frac{2^j j!}{(2j+1)!}x^{2j}$$ (adapted from formula 7.1.6)

Note that the coefficients of $x^{2j}$ in the series $c_j=\frac{2^j j!}{(2j+1)!}$ can be computed by starting with $c_0=1$ and then using the recursion formula $c_{j+1}=\frac{c_j}{2j+3}$. This is convenient when implementing the series as a summation loop.

---

cardinal gave the Laplacian continued fraction as a way to bound Mills's ratio for large $|x|$; what is not as well-known is that the continued fraction is also useful for numerical evaluation.

Lentz, Thompson and Barnett derived an algorithm for numerically evaluating a continued fraction as an infinite product, which is more efficient than the usual approach of computing a continued fraction "backwards". Instead of displaying the general algorithm, I'll show how it specializes to the computation of Mills's ratio:

$\displaystyle Y_0=x,\,C_0=Y_0,\,D_0=0$
$\text{repeat for }j=1,2,\dots$
$$D_j=\frac1{x+jD_{j-1}}$$ $$C_j=x+\frac{j}{C_{j-1}}$$ $$H_j=C_j D_j$$ $$Y_j=H_j Y_{j-1}$$ $\text{until }|H_j-1| < \text{tol}$
$\displaystyle R(x)=\frac1{Y_j}$

where $\text{tol}$ determines the accuracy.

The CF is useful where the previously mentioned series starts to converge slowly; you will have to experiment with determining the appropriate "break point" to switch from the series to the CF in your computing environment. There is also the alternative of using an asymptotic series instead of the Laplacian CF, but my experience is that the Laplacian CF is good enough for most applications.

---

Finally, if you don't need to compute the (complementary) error function very accurately (i.e., to only a few significant digits), there are compact approximations due to Serge Winitzki. Here is one of them:

$$R(x)\approx \frac{\sqrt{2\pi}+x(\pi-2)}{2+x\sqrt{2\pi}+x^2(\pi-2)}$$

This approximation has a maximum relative error of $1.84\times 10^{-2}$ and becomes more accurate as $x$ increases.

Answer 3

(This reply originally appeared in response to a similar question, subsequently closed as a duplicate. The O.P. only wanted "an" implementation of the Gaussian integral, not necessarily "state of the art." In his comments it became apparent that a relatively simple, short implementation would be preferred.)

---

As comments point out, you need to integrate the PDF. There are many ways to perform the integral. Long ago, when computations were slow and expensive, David Hill worked out an approximation using simple arithmetic (rational functions and an exponentiation). It has double precision accuracy for typical arguments (between $-8.5$ and $+8.5$, approximately). In 1973 he published a Fortran version in Applied Statistics called ALNORM.F. Over the years I have ported this to various environments which did not have a Normal (Gaussian) integral or which had suspect ones (such as Excel).

A MatLab version (with appropriate attributions) is available at http://people.sc.fsu.edu/~jburkardt/m_src/asa005/alnorm.m. A completely undocumented version of the original Fortran code appears on a "Koders Code Search" (sic) site.

Many years ago I ported this to AWK. This version may be more congenial for the modern developer to port due to its C-like (rather than Fortran) syntax and some additional comments I inserted when developing and testing it, because I needed to enhance its accuracy. It appears below.

For those without much experience porting scientific/math/stats code, some words of advice: one single typographical mistake can create serious errors that might not be easily detectable. (Trust me on this, I've made lots of them.) Always, always create a careful and exhaustive test. Because the normal integral/Gaussian integral/error function is available in so many tables and so much software, it's simple and fast to tabulate a huge number of values of your ported function and systematically compare (i.e., with the computer, not by eye) the values to correct ones. You can see such a test at the beginning of my code: it produces a table of values in -8.5:8.5 (by 0.1) which can be piped (via STDOUT) to another program for systematic checking.

Another testing approach--for those with enough numerical analysis background to know how to estimate expected errors--would be to numerically differentiate the values and compare them to the PDF (which is readily computed).

By the way: this code is only for the case with a mean of $0$ and unit standard deviation ("sigma"). But that's all one needs: to integrate from $-\infty$ to $x$ when the mean is $\mu$ and the SD is $\sigma$, just compute $z = (x-\mu)/\sigma$ and apply alnorm to it.

Edit

I tested a port of alnorm to Mathematica, which computes the values to arbitrary precision. To compare the results, here is a plot of the natural log of the ratios of upper tail values $1 - \Phi(z)$ with $z\ge 1$. (A positive relative error means alnorm is too large.)

The values are always accurate to $4 \times 10^{-11}$ relative to the vanishingly small tail probabilities. You can see where the calculation switches to an asymptotic formula (at $z=16$) and it is evident that this formula becomes extremely accurate as $z$ increases. The plot stops at $z=\sqrt(2\times 708) \approx 37.6$ because here is where double-precision exponentiation begins underflowing.

For example, alnorm[-6.0] returns $9.865\ 876\ 450\ 315E-10$ while the true value, equal to $\frac{1}{2}\text{erfc}(3\sqrt{2})$, is approximately $9.865\ 876\ 450\ 377E-10$, first differing in the twelfth decimal digit.

NB As part of this edit, I changed UPPER_TAIL_IS_ZERO from 15. to 16. in the code: it makes the result a tiny bit more accurate for $Z$ between $15$ and $16$. (End of edit.)

#----------------------------------------------------------------------#
#   ALNORM.AWK
#   Compute values of the cumulative normal probability function.
#   From G. Dallal's STAT-SAK (Fortran code).
#   Additional precision using asymptotic expression added 7/8/92.
#----------------------------------------------------------------------#
BEGIN {
    for (i=-85; i<=85; i++) {
        x = i/10
        p = alnorm(x, 0)
        printf("%3.1f %12.10f\n", x, p)
    }
    exit
}
function alnorm(z,up,    y,aln,w) {
#
#    ALGORITHM AS 66 APPL. STATIST. (1973) VOL.22, NO.3:
#    Hill,  I.D.  (1973).  Algorithm AS 66.  The normal  integral.
#                          Appl. Statist.,22,424-427.
#
#    Evaluates the tail area of the standard normal curve from
#    z to infinity if up, or from -infinity to z if not up.
#
#    LOWER_TAIL_IS_ONE, UPPER_TAIL_IS_ZERO, and EXP_MIN_ARG
#    must be set to suit this computer and compiler.

    LOWER_TAIL_IS_ONE = 8.5     # I.e., alnorm(8.5,0) = .999999999999+
    UPPER_TAIL_IS_ZERO = 16.0   # Changes to power series expression
    FORMULA_BREAK = 1.28        # Changes cont. fraction coefficients
    EXP_MIN_ARG = -708          # I.e., exp(-708) is essentially true 0

    if (z < 0.0) {
        up = !up
        z = -z
    }
    if ((z <= LOWER_TAIL_IS_ONE) || (up && z <= UPPER_TAIL_IS_ZERO)) {
        y = 0.5 * z * z
        if (z > FORMULA_BREAK) {
            if (-y > EXP_MIN_ARG) {
                aln = .398942280385 * exp(-y) / \
                  (z - 3.8052E-8 + 1.00000615302 / \
                  (z + 3.98064794E-4 + 1.98615381364 / \
                  (z - 0.151679116635 + 5.29330324926 / \
                  (z + 4.8385912808 - 15.1508972451 / \
                  (z + 0.742380924027 + 30.789933034 / \
                  (z + 3.99019417011))))))
            } else {
                aln = 0.0
            }
        } else {
            aln = 0.5 - z * (0.398942280444 - 0.399903438504 * y / \
              (y + 5.75885480458 - 29.8213557808 / \
              (y + 2.62433121679 + 48.6959930692 / \
              (y + 5.92885724438))))
        }
    } else {
        if (up) {   # 7/8/92
            # Uses asymptotic expansion for exp(-z*z/2)/alnorm(z)
            # Agrees with continued fraction to 11 s.f. when z >= 15
            # and coefficients through 706 are used.
            y = -0.5*z*z
            if (y > EXP_MIN_ARG) {
                w = -0.5/y  # 1/z^2
                aln = 0.3989422804014327*exp(y)/ \
                    (z*(1 + w*(1 + w*(-2 + w*(10 + w*(-74 + w*706))))))
                # Next coefficients would be -8162, 110410
            } else {
                aln = 0.0
            }
        } else {
            aln = 0.0
        }
    }
    return up ? aln : 1.0 - aln
}
### end of file ###