6
$\begingroup$

In answer https://stackoverflow.com/a/23119456/2421256, an approximation of the complementary normal CDF (ie $\frac{1}{\sqrt{2 \pi}} \int_x^{+\infty} e^{-\frac{t^2}{2}} dt$) was given for $0 \leq x \leq 7$ using a rational function as follows : $$ \frac{1}{\sqrt{2 \pi}} e^{\frac{-x^2}{2}} \frac{N_0 + N_1 x + ... + N_6 x^6}{M_0 + M_1 x + ... + M_7 x^7} $$

How can one derive such an approximation,without relying on an already implemented approximation (ie no methods based on curve fitting), and a bound for the absolute error ?

EDIT : I added a bounty because I believe this question is not at all trivial.

EDIT 2 : I'm adding here what I've been able to do for now. Indeed, as Ian pointed out, for $0 \leq x \leq \sqrt{2}$ a good approximation of the normal CDF can be derived easily, by using a Taylor series expansion of $e^{-\frac{t^2}{2}}$ and integrating it term by term between $0$ and $x$.

For $\sqrt{2} \leq x \leq 7$, I proved, using Mill's inequality and the fact that $\frac{\int_x^{+\infty} e^{-\frac{t^2}{2}} dt}{x e^{-\frac{x^2}{2}}}$ is increasing in $x$, that $$\mathbb{P} \left( Z > x \right) \leq \frac{x e^{-\frac{x^2}{2}}}{98 \pi} $$, which yields an error smaller than $2.39 \times 10^{-3} $, but it is far from the precision yielded by the rational approximation above (around $10^{-16}$, checked numerically since I don't have a rigourous bound for it).

EDIT 3 : if it's of help, the target error is something smaller than $10^{-9} $.

$\endgroup$
7
  • $\begingroup$ For slightly larger $x$ (I think it turns out to be $x>\sqrt{2}$), you can derive the bounds $p(x) e^{-x^2/2} \leq F(x) \leq q(x) e^{-x^2/2}$, where $p,q$ are rational functions that can be computed explicitly. This means that $e^{x^2/2} F(x)$ is "asymptotically like" a rational function. (Try plotting it; past about $x=30$ you get floating point issues, but until then you can see what I mean.) You can then try to approximate this instead, and then stitch together the approximation close to $0$ and close to $\infty$ using Pade approximation or similar techniques. $\endgroup$
    – Ian
    May 8, 2016 at 18:30
  • 1
    $\begingroup$ @BS. I don't have access to the paper, but references in literature suggest that the following paper may be of interest, giving an approximation precisely for the interval you are interested in: P. A. P. Moran, "Calculation of the normal distribution function", Biometrika 67 (1980), pp. 675-677. The approximation seems to be an expansion in trig functions best I can tell. $\endgroup$
    – njuffa
    May 18, 2016 at 14:34
  • 1
    $\begingroup$ C implementation of Moran is double cdf (double z) { double pi = 3.14159265358979323; double sqrt2 = sqrt (2.0); double sum = 0.0; for (int k = 1; k <= 12; k++) { sum += (1.0 / k) * exp (-k * k / 9.0) * sin (k * z * sqrt2 / 3); } return 0.5 + (z / (3 * sqrt2) + sum) / pi; } , accurate to about nine decimal places. Literature references indicate Moran's work is based on Anthony J. Strecok, "On the Calculation of the Inverse of the Error Function", Math. Comp., v. 22, 1968, pp. 144-158. $\endgroup$
    – njuffa
    May 18, 2016 at 15:14
  • 1
    $\begingroup$ Sorry, I linked the wrong paper in my preceding comment. The Strecok paper I referenced can be found here $\endgroup$
    – njuffa
    May 18, 2016 at 15:31
  • 1
    $\begingroup$ @BS. Also of interest may be this paper: George Marsaglia, "Evaluating the Normal Distribution", Journal of Statistical Software, July 2004, Volume 11, Issue 4, pp. 1-11. The approximation for $\Phi$ is claimed to be accurate to almost double precision over your targeted interval using a very simple expansion. Online copy can be found here $\endgroup$
    – njuffa
    May 18, 2016 at 16:56

3 Answers 3

3
+100
$\begingroup$

This isn't a definitive answer, but hopefully points you in the right direction. The short of it is, I don't know how the coefficients were derived, but it is likely related to some Pade, Remez, interpolating, or least-squares approximant. Read on for details.

The function we're talking about is $$ Q(x) = \frac{1}{\sqrt{2\pi}}\int_{x}^\infty e^{-t^2/2}\;dt = \frac{1}{2}\operatorname{erfc}\left(\frac{x}{\sqrt{2}}\right), $$ for $x>0$, where $\operatorname{erfc}$ is the complimentary error function.

The link you provide gives code for an approximation of the following form, slightly different than how it was stated in the original question (note the lack of normalization by $\sqrt{2\pi}$ and a different form of the polynomials). $$ \tilde{Q}(x) = e^{\frac{-x^2}{2}} \frac{(((((N_6 z+N_5)z+N_4)z+N_3)z+N_2)z+N_1)z+N_0}{((((((M_7 z+M_6)z+M_5)z+M_4)z+M_3)z+M_2)z+M_1)Z+M_0} $$ Why do we pick this form? Well, asymptotically, the $Q$ function is well approximated by $\frac{e^{-x^2/2}}{x\sqrt{2\pi}}$ as $x\rightarrow\infty$. See here for details on that. This justifies factoring out the asymptotic behavior and looking for a rational approximation of what remains. If we factored out the asymptotic form, we could approximate what remains with something that goes to $1$ for large $x$, i.e., something with equal order in the numerator and denominator. If we leave that extra factor of $x$ in the denominator from the asymptotic piece, we can account for it by absorbing it into the denominator of the rational function, giving one higher degree in the denominator polynomial. The same goes for the factor of $\sqrt{2\pi}$, it gets absorbed in to the approximation.

Now, how on earth did Hart come up with those $N_i$ and the $M_j$ coefficients? I thought it was by requiring that the Taylor series of $\tilde{Q}(x)e^{x^2/2}$ matches that of $Q(x)e^{x^2/2}$ at some a point, up to a high enough order to uniquely determine all the coefficients. This is the Pade approximation technique, and you should be able to find more about the error bounds for approximations resulting from this technique. I used Mathematica to calculate these below.

I thought that surely, if you calculate that Pade approximant around $x=0$, and scale the coefficients to match this weird $M_0=440.413735824752$ (Pade approximants are usually normalized so that the constant term in the denominator is $1$) you'd get the formula. Instead, I get:

$$ \begin{array}{ccc} i & N_i & M_i\\ 0 & 220.207 & 440.414\\ 1 & 152.909 & 711.438\\ 2 & 95.568 & 525.241 \\ 3 & 23.168 & 234.066 \\ 4 & 4.6336 & 61.9813\\ 5 & 0.34752 & 10.2661 \\ 6 & 0.154453 & 1.04799 \\ 7 & {} & 0.048886 \end{array} $$

It's a similar pattern to Hart's numbers, and the values are all in the same ballpark as Hart's numbers, but it is not a match. Then I thought that there is this weird breakpoint in the code that appears to be $\frac{10}{\sqrt{2}}\approx 7.07$. This lead me down the path of looking to see if the Pade approximant near $\frac{10}{\sqrt{2}}$ or at the halfway point of the interval, $\frac{10}{2\sqrt{2}}$, appropriately scaled, gave Hart's original coefficients. They don't.

So here are some possibilities.

  1. Hart's coefficients are for a Pade approximant around SOME point on the interval $[0,10/\sqrt{2}]$, and I don't have the patience to find which point.
  2. Hart used some other method to generate these coefficients, and they don't correspond to a Pade approximant at all. One candidate is the Remez algorithm for rational functions, which should guarantee the error oscillates equally across the domain of approximation. Another candidate is that the coefficients are derived from an interpolation between some fixed points with known values of the domain. Yet another is a least-squares fit between the rational approximant a a bunch of known values of the function.

This has been a long winded answer, and boils down to "I don't know". But I hope the techniques I mentioned can be fruitful avenues for you to continue looking into this. I would finally suggest getting the original "Algorithm 5666 for the error function" and seeing if there is an explanation of the methodology in book in which it appears.

$\endgroup$
2
  • $\begingroup$ Excellent answer and very useful insights. Props also for the references ! $\endgroup$
    – BS.
    May 18, 2016 at 22:43
  • $\begingroup$ I have not checked the introductory chapters in Hart, et. al. in many years, but am reasonably certain that the authors used variants of the Remez algorithm to derive the coefficients for many (probably most) minimax approximations in the book, both the polynomial and the rational ones. I understood the present question as a request for alternative approximations not derived by such purely numerical methods as the Remez algorithm. $\endgroup$
    – njuffa
    May 18, 2016 at 23:20
1
$\begingroup$

HINT

Required rational approximation for $\operatorname{erfc}$ function can be constructed analytically, on the basis of known Maclaurin series $$\sqrt\pi ze^{z^2}\operatorname{erfc}z = v,$$ $$v=1-\frac{1\cdot3}{2z^2}+\frac{1\cdot3\cdot5}{(2z^2)^2}-\frac{1\cdot3\cdot5\cdot7}{(2z^2)^3}+\frac{1\cdot3\cdot5\cdot7\cdot9}{(2z^2)^4}-\frac{1\cdot3\cdot5\cdot7\cdot9\cdot11}{(2z^2)^5}+\cdots,$$ and $$(1+a_1x+a_2x^2+a_3x^3+a_4x^4+\dots)^{-1}=$$ $$=1-a_1x+(a_1^2-a_2)x^2+(-a_1^3+2a_1a_2-a_3)x^3+(a_1^4-3a_1^2a_2+2a_1a_3+a_2^2-a_4)x^4+\dots$$ Then $$v=1+\frac1{v_1},\quad v_1=\frac{1}{v-1} $$$$v_1=\left(-\frac{1\cdot3}{2z^2}+\frac{1\cdot3\cdot5}{(2z^2)^2}-\frac{1\cdot3\cdot5\cdot7}{(2z^2)^3}+\frac{1\cdot3\cdot5\cdot7\cdot9}{(2z^2)^4}-\frac{1\cdot3\cdot5\cdot7\cdot9\cdot11}{(2z^2)^5}+\dots\right)^{-1} $$$$ = -\frac{2z^2}3\left(1 - \frac{5}{2z^2} + \frac{5\cdot7}{(2z^2)^2} - \frac{5\cdot7\cdot9}{(2z^2)^3} + \frac{5\cdot7\cdot9\cdot11}{(2z^2)^4}+\dots\right)^{-1} $$$${= -}\frac{2z^2}3\left(1 + \frac{5}{2z^2} + \frac{5^2-5\cdot7}{(2z^2)^2} + \frac{-(-5)^3+2\cdot(-5)\cdot(5\cdot7)-(-5\cdot7\cdot9)}{(2z^2)^3}+ \frac{(-5)^4-3\cdot(-5)^2\cdot(5\cdot7)+2\cdot(-5)\cdot(-5\cdot7\cdot9)+(5\cdot7)^2-5\cdot7\cdot9\cdot11}{(2z^2)^4}+\dots\right)$$ $$= -\frac{2z^2}3 - \frac53 + \frac{10}{3\cdot2z^2} - \frac{30}{(2z^2)^2}+ \frac{1090}{3(2z^2)^3}+\dots,$$ $$v_1{= -}\frac{2z^2}3 - \frac53 + \frac1{v_2},\quad v_2=\dfrac1{v_1 + \dfrac{2z^2}3 + \dfrac53}$$ $$v_2 = \left(\frac{10}{3\cdot2z^2} - \frac{30}{(2z^2)^2} + \frac{1090}{3\cdot(2z^2)^3}+\dots\right)^{-1} = \frac{3\cdot2z^2}{10}\left(1 - \frac{9}{2z^2} + \frac{109}{(2z^2)^2}+\dots\right)^{-1}$$ etc.

The advantage of this approach - that we can indefinitely increase the accuracy of the basic series, and thus - the approximation accuracy.

$\endgroup$
0
$\begingroup$

In Maple, you can do it like this (which uses the Remez algorithm):

with(numapprox):
Digits:= 20: 
Fapp:= minimax(F,x=0..7,[7,6],1,'maxerr');

Fapp := (0.19114016230602883109e-2+(0.16638131476248234372e-2+(0.73397713688511034855e-3+(0.18372599982443885556e-3+(0.25987590927945377226e-4+(0.16534245454532854223e-5+(2.9989826198693857705*10^(-11)-5.4618378185344946027*10^(-13)*x)*x)*x)*x)*x)*x)*x)/(0.38228032461211988467e-2+(0.6377781984229194914e-2+(0.46452864301009549948e-2+(0.19016818756092056203e-2+(0.46504884699484812014e-3+(0.65109362703586754431e-4+0.41464699141737875111e-5*x)*x)*x)*x)*x)*x)

The (approximate) maximum error:

maxerr;

8.138333*10^(-14)

To verify, plot the difference:

plot(Fapp-F,x=0..7);

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.