Approximating the
Error Function
and the
Cumulative Normal Distribution
Daniel D. Warner
June 6, 1999
The Error Function
The error function is defined by
, and the complementary error function is defined by
.
The following approximation is a translation of a FORTRAN program by W. J. Cody, Argonne National Laboratory, NETLIB/SPECFUN, March 19, 1990. The main computation evaluates the near-minimax approximations presented in "Rational Chebyshev approximations for the error function" by W. J. Cody, Math. Comp., 1969, PP. 631-638.
Explanation of machine-dependent constants:
XMIN = the smallest positive floating-point number.
XINF = the largest positive finite floating-point number.
XNEG = the largest negative argument acceptable to ERFCX;
the negative of the solution to the equation
2*exp(x*x) = XINF.
XSMALL = argument below which erf(x) may be represented by
2*x/sqrt(pi) and above which x*x will not underflow.
A conservative value is the largest machine number X
such that 1.0 + X = 1.0 to machine precision.
XBIG = largest argument acceptable to ERFC; solution to
the equation: W(x) * (1-0.5/x**2) = XMIN, where
W(x) = exp(-x*x)/[x*sqrt(pi)].
XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to
machine precision. A conservative value is
1/[2*sqrt(XSMALL)]
XMAX = largest acceptable argument to ERFCX; the minimum
of XINF and 1/[sqrt(pi)*XMIN].
C------------------------------------------------------------------
C Machine-dependent constants for the CASIO FX2.0
C------------------------------------------------------------------
DATA XINF,XNEG,XSMALL/9.99999999999999D+99, -15.1514145, 1.0D-15/,
1 XBIG,XHUGE,XMAX/14.98912960D0, 1.581138834D+7, 5.642D+98/
>
Digits := 22;
XMIN := 1.00000000000000*10^(-99);
XINF := 9.99999999999999*10^99;
XNEG := -sqrt(log(XINF/2));
XSMALL := 1.0*10^(-15);
> W := x -> exp(-x^2)/(x*sqrt(Pi));
>
XBIG := fsolve(W(x)*(1-0.5/x**2) = XMIN, x=15);
XHUGE := 1/(2*sqrt(XSMALL));
XMAX := evalf(min(XINF,1/(sqrt(Pi)*XMIN)));
>
Erfcore := proc( xx::algebraic, k::integer ) local a, b, c, d, p, q, j, y, x, xn, xd, z;
a := [3.16112374387056560, 113.864154151050156, 377.485237685302021,
3209.37758913846947, 0.185777706184603153];
b := [23.6012909523441209, 244.024637934444173,
1282.61652607737228, 2844.23683343917062];
c := [0.564188496988670089, 8.88314979438837594, 66.1191906371416295,
298.635138197400131, 881.952221241769090, 1712.04761263407058,
2051.07837782607147, 1230.33935479799725, 0.0000000215311535474403846];
d := [15.7449261107098347, 117.693950891312499, 537.181101862009858,
1621.38957456669019, 3290.79923573345963, 4362.61909014324716,
3439.36767414372164, 1230.33935480374942];
p := [0.305326634961232344, 0.360344899949804439, 0.125781726111229246,
0.0160837851487422766, 0.000658749161529837803, 0.0163153871373020978];
q := [2.56852019228982242, 1.87295284992346047, 0.527905102951428412,
0.0605183413124413191, 0.00233520497626869185];
x := evalf(xx);
if (abs(x) < (1/2)) then
y := x^2;
xn := a[5]*y;
xd := y;
for j from 1 to 3 do
xn := (xn + a[j])*y;
xd := (xd + b[j])*y;
od;
z := evalf(x*(xn + a[4])/(xd + b[4]));
if k = 0 then RETURN( z ) else RETURN( 1-z ) fi;
elif abs(x) < 4 then
y := abs(x);
xn := c[9]*y;
xd := y;
for j from 1 to 7 do
xn := (xn + c[j])*y;
xd := (xd + d[j])*y;
od;
z := evalf(exp(-y^2)*(xn + c[8])/(xd + d[8]));
if (k = 0) then
if (x < 0) then RETURN( z-1 ) else RETURN( 1-z ) fi;
else
if (x < 0) then RETURN( 2-z ) else RETURN( z ) fi;
fi;
else
y := 1/x^2;
xn := p[6]*y;
xd := y;
for j from 1 to 4 do
xn := (xn + p[j])*y;
xd := (xd + q[j])*y;
od;
z := evalf(exp(-x^2)*(1/sqrt(Pi) - y*(xn + p[5])/(xd + q[5]))/abs(x));
if (k = 0) then
if (x < 0) then RETURN( z-1 ) else RETURN( 1-z ) fi;
else
if (x < 0) then RETURN( 2-z ) else RETURN( z ) fi;
fi;
fi;
end:
> Erf1 := x-> Erfcore(x,0):
> plot(Erf1, -5..5, title="Plot of the Approximated Error Function");
> plot((erf - Erf1), -5..5, numpoints=2048, title="Detailed Plot of the Errors");
> Erfc1 := x-> Erfcore(x,1):
> plot(Erfc1, -5..5, title="Plot of the Approximated Complementary Error Function");
> plot((erfc - Erfc1), -5..5, title="Detailed Plot of the Errors");
Normal Distribution
The Normal Probability Density
The normal probability density is given by
,
where
.
>
sigma := `sigma`: mu := `mu`:
Npd := (x,sigma,mu) -> (1/(sigma*sqrt(2*Pi)))*exp(-((x-mu)/(sigma*sqrt(2)))^2);
> plot(Npd(x,1,0),x=-3..3);
The Cumulative Normal Probability Distribution
The Normal Probability Distribution,
, is given by
,
where
is the lower boundary and
is the upper boundary.
> Ncd0 := (a,b,sigma,mu) -> int(Npd(x,sigma,mu),x=a..b);
>
Ncd1 := proc(a,b,sigma,mu) local u, v, s, t, z;
u := evalf((b-mu)/(sigma*sqrt(2.0)));
v := evalf((a-mu)/(sigma*sqrt(2.0)));
if (v <= u) then
s := 1;
else
s := -1;
t := u;
u := v;
v := t;
fi;
if (u < (-0.5)) or (0.5 < v) then
z := evalf(s*0.5*(Erfcore(abs(u),1) - Erfcore(abs(v),1)));
else
z := evalf(s*0.5*(Erfcore(u,0) - Erfcore(v,0)));
fi;
RETURN( z );
end;
> Ntest0 := x -> Ncd0(-infinity, x, 1.0, 0.0);
> Ntest1 := x -> Ncd1(-10000, x, 1.0, 0.0);
> plot((Ntest1), -3..3, title="Plot of the Approximated Cumulative Normal Distribution");
> plot((Ntest0 - Ntest1), -3..3, title="Detailed Plot of the Errors");
> plot(((Ntest0 - Ntest1)/Ntest0), -3..3, title="Detailed Plot of the Relative Errors");
Inverse Error Function
The inverse error function is defined by
returns the value
such that
.
The following approximation is based on the Matlab program by Cleve Moler. It uses a simple rational funciton approximation to the inverse error function followed by two steps of the Newton-Raphson iteration. The rational approximation is accurate to about 6 digits.
>
InvErf := proc( yy::algebraic ) local a, b, c, d, x, y, z;
a := [ 0.886226899, -1.645349621, 0.914624893, -0.140543331];
b := [-2.118377725, 1.442710462, -0.329097515, 0.012229801];
c := [-1.970840454, -1.624906493, 3.429567803, 1.641345311];
d := [ 3.543889200, 1.637067800];
y := evalf(yy);
if abs(y) < 0.7 then
z := y^2;
x := y*(((a[4]*z+a[3])*z+a[2])*z+a[1]) / ((((b[4]*z+b[3])*z+b[2])*z+b[1])*z+1);
elif (0 < y) then
z := sqrt(-log((1-y)/2));
x := (((c[4]*z+c[3])*z+c[2])*z+c[1]) / ((d[2]*z+d[1])*z+1);
else
z := sqrt(-log((1+y)/2));
x := -(((c[4]*z+c[3])*z+c[2])*z+c[1]) / ((d[2]*z+d[1])*z+1);
fi;
x := evalf(x - (Erfcore(x,0) - y) / (2/sqrt(Pi) * exp(-x^2)));
x := evalf(x - (Erfcore(x,0) - y) / (2/sqrt(Pi) * exp(-x^2)));
RETURN( x );
end:
> x0 := InvErf(-0.8);
> Erf1(x0);
>
Id0 := x -> x;
Id1 := x -> Erf1(InvErf(x));
> plot((Id0 - Id1), -0.99..0.99);
The Inverse Cumulative Normal Distribution
The Inverse Cumulative Normal Probability Distribution finds
given
,
> Ncd2 := (a,b,sigma,mu) -> (1/2)*(erf((b-mu)/(sigma*sqrt(2))) - erf((a-mu)/(sigma*sqrt(2))));
> InvN0 := (p,sigma,mu) -> fsolve(Ncd2(-100,x,sigma,mu)=p,x);
> InvN1 := (p,sigma,mu) -> evalf(mu + sqrt(2)*sigma*InvErf(2*p-1));
>
Ntest2 := x -> InvN0(x,1,0);
Ntest3 := x -> InvN1(x,1,0);
> plot((Ntest2 - Ntest3), 0..1, numpoints=1001, title="Detailed Plot of the Errors");
> plot(((Ntest2 - Ntest3)/Ntest2), 0..1, numpoints=1001, title="Detailed Plot of the Relative Errors");
> Ntest2(0.7);
> Ntest3( 0.999999999 );
> Erfc1( 6 );
>