'erf' function in C

[Jesper Lund]

It is possible to compress the code further at the expense of
clarity (warning: this is not for the faint of heart):

double Phi(double x)
{long double s=x,t=0,b=x,x2=x*x,i=1;
while(s!=t) s=(t=s)+(b*=x2/(i+=2));
return .5+s*exp(-.5*x2-.91893853320467274178L);
}

This avoids the unnecessary computation of i+1 followed by the
unnecessary conversion of the result from int to long double,
eliminates the unnecessary variable pwr, and precomputes x*x.

Have you checked this algorithm in the tail? For, say, x = -50, I either
get inf/nan or an incorrect answer (0.5), depending on which compiler I
use. The compiler in my tests is gcc, specifically DJGPP or MingW. The
problem is overflow in s.

Furthermore, Phi(x) is not increasing in x for values between -13 and -7,
see below (DJGPP using long double).

x Phi(x)
-13 -2.94464e-15
-12 -2.89807e-15
-11 6.27672e-16
-10 6.49085e-16
-9 6.65294e-16
-8 1.28093e-15
-7 1.28042e-12

Best regards,
Jesper Lund

 

[Deane Yang]

Yes, the paper cited below is the one I used to implement
the inverse cumulative normal function. In that paper
you recommended that such a function be used to generate
random normally distributed numbers, instead of the Box-Muller
method.

1) Would you still recommend this?
2) The approach suggested below seems to be aiming more
for accuracy, rather than speed. On the other hand,
the published algorithm is machine-dependent and therefore
a bit of a headache to maintain.

Is there another fast but portable algorithm with sufficient
accuracy available?

And what about the inverse cumulative normal distribution?
I used to implement an algorithm you published,
but do not like it because it is machine-dependent.
What in your opinion is the best algorithm?

With the improvement on my method for
evaluating cPhi(x) as R(x)*phi(x),
when the Taylor series for R(x) is about zero,
(suggested by Bill Daly), the simple C function

double Phi(double x)
{ long double t=0,b=1,s=x,pwr=x;
int i;
for(i=2;s!=t;i+=2)
{ b/=(i+1); pwr=pwr*x*x; t=s; s+=pwr*b;}
return .5+s*exp(-.5*x*x-.91893853320467274178L);
}

should serve very well for solving the equation
Phi(X)=U
for positive X, given U>1/2, by Newton's method:
Start with an initial estimate

x=sqrt(-1.6*ln(1-U^2));

then repeat

x=x-(Phi(x)-U)/phi(x);

until you get no further changes in x.


The paper that you refer to must be where I want
to generate a normal X by means of solving
Phi(X)=U, given a random U, uniform in [0,1),
by getting an integer j, formed from certain
bits of the exponent part of the floating point
representation of U. Then use j to access tabled
values A[j] and B[j] to initalize the Taylor series
for Phi inverse, using the fractional part of U
as the argument in the series. It is designed for
speed in generating a normal variate directly
as Phi^(-1)(U), but only provides accuracy to 6/7 digits,
(As I recall, the Taylor series was easy and fast only for the first few
terms.)

It is described in
``Normal (Gaussian) random variables in supercomputers'', (1991)
Journal of Supercomputing}, v 5, 49--55,
where the method is particularly suited for parallel computation.

George Marsaglia

 

[George Marsaglia]

Have you checked this algorithm in the tail? For, say, x = -50, I either
get inf/nan or an incorrect answer (0.5), depending on which compiler I
use. The compiler in my tests is gcc, specifically DJGPP or MingW. The
problem is overflow in s.

Furthermore, Phi(x) is not increasing in x for values between -13 and -7,
see below (DJGPP using long double).

x Phi(x)
-13 -2.94464e-15
-12 -2.89807e-15
-11 6.27672e-16
-10 6.49085e-16
-9 6.65294e-16
-8 1.28093e-15
-7 1.28042e-12

Best regards,
Jesper Lund

The attraction of the method for evaluating
cPhi(x)=1-Phi(x)
through the intermediary
R(x)=cPhi(x)/phi(x)
is that R(x) has an easily developed Taylor
series. Thus an easy way to evaluate R(z+h),
given a known value a=R(z), is through

b=z*a-1;
pwr=1;
s=a+h*b;
for(i=2;s!=t;i+=2)
{ a=(a+z*b)/i;
b=(b+z*a)/(i+1);
pwr=pwr*h*h;
t=s;
s=s+pwr*(a+h*b);
}

It happens that even using the Taylor series
for z=0 provides better accuracy for Phi(x)
than most users are likely to need for
applications in probability and statistics.
But at the expense of speed and restriction
to practical values of x , say, -7<x<7.

Those needing greater accuracy in the tails,
or, more generally, relative rather than
absolute accuracy, should use the above Taylor
expansion for various values of z, where
15-16 digit relative accuracy is attainable.
(The series about zero will provide only a
few digits of relative accuracy for large |x|.)

Some of my early implementations of the method
(one of which on newsgroups several years ago),
kept a table of R(z) for 128 or 64 values z,
and the series needed only a few steps to terminate.

Subsequent exchanges with interested parties led
to consideration of 8 then 4 and finally one value,
z=0, for which evaluation of the odd terms in the
Taylor series can be avoided, thanks to the
observation of Bill Daly.

So, those content with the accuracies indicated
in earlier parts of this thread may wish to use
the simpler Taylor series about zero. Those
seeking faster routines or greater relative
accuracy should consider the not-quite-as-easy
Taylor series for R(z+h), using as many tabled
values of R(z) as their speed and accuracy
requirements call for.

George Marsaglia

 

[Axel Vogt]

Have you checked this algorithm in the tail? For, say, x = -50, I either
get inf/nan or an incorrect answer (0.5), depending on which compiler I
use. The compiler in my tests is gcc, specifically DJGPP or MingW. The
problem is overflow in s.

For large arguments one might want to switch to asymptotics.
I think, that is what Genz makes in his approx for cdf normal.

 

[Jesper Lund]

For large arguments one might want to switch to asymptotics.
I think, that is what Genz makes in his approx for cdf normal.

I found the following algorithm (using Google)
http://tigger.smu.edu.sg/software/mnp-stuff/stat.ubc.ca/pnorms2.c
developed by W.D. Cody (there is a full reference in the link). It works
for the entire domain, and produces Phi(-37.5) = 4.60535e-308 with
double precision. I don't have access to the paper by Cody, but
according to the comments in the text, the algorithm is accurate to 18
significant decimal places.

 

[Jesper Lund]

The attraction of the method for evaluating
cPhi(x)=1-Phi(x)
through the intermediary
R(x)=cPhi(x)/phi(x)
is that R(x) has an easily developed Taylor
series. Thus an easy way to evaluate R(z+h),
given a known value a=R(z), is through

[stuff deleted]

Some of my early implementations of the method
(one of which on newsgroups several years ago),
kept a table of R(z) for 128 or 64 values z,
and the series needed only a few steps to terminate.

Subsequent exchanges with interested parties led
to consideration of 8 then 4 and finally one value,
z=0, for which evaluation of the odd terms in the
Taylor series can be avoided, thanks to the
observation of Bill Daly.

Thanks for explaining the finer details of your algorithm here. I have
already found the earlier versions of your algorithm [specifically, the
ones with 8 and 128 pre-calculated values for R(z)] using Google Groups.

Are you familiar with the algorithm developed by W.D. Cody? See this
link for a C implementation with a stated accuracy of 18 significant
decimal digits
http://tigger.smu.edu.sg/software/mnp-stuff/stat.ubc.ca/pnorms2.c

 

[Rob Johnson]

The attraction of the method for evaluating
cPhi(x)=1-Phi(x)
through the intermediary
R(x)=cPhi(x)/phi(x)
is that R(x) has an easily developed Taylor
series. Thus an easy way to evaluate R(z+h),
given a known value a=R(z), is through

b=z*a-1;
pwr=1;
s=a+h*b;
for(i=2;s!=t;i+=2)
{ a=(a+z*b)/i;
b=(b+z*a)/(i+1);
pwr=pwr*h*h;
t=s;
s=s+pwr*(a+h*b);
}

It happens that even using the Taylor series
for z=0 provides better accuracy for Phi(x)
than most users are likely to need for
applications in probability and statistics.
But at the expense of speed and restriction
to practical values of x , say, -7<x<7.

Those needing greater accuracy in the tails,
or, more generally, relative rather than
absolute accuracy, should use the above Taylor
expansion for various values of z, where
15-16 digit relative accuracy is attainable.
(The series about zero will provide only a
few digits of relative accuracy for large |x|.)

Some of my early implementations of the method
(one of which on newsgroups several years ago),
kept a table of R(z) for 128 or 64 values z,
and the series needed only a few steps to terminate.

Subsequent exchanges with interested parties led
to consideration of 8 then 4 and finally one value,
z=0, for which evaluation of the odd terms in the
Taylor series can be avoided, thanks to the
observation of Bill Daly.

So, those content with the accuracies indicated
in earlier parts of this thread may wish to use
the simpler Taylor series about zero. Those
seeking faster routines or greater relative
accuracy should consider the not-quite-as-easy
Taylor series for R(z+h), using as many tabled
values of R(z) as their speed and accuracy
requirements call for.

Another method for approximating erf(x) for large values of x is through
an asymptotic series. Formula [8] in

http://www.whim.org/nebula/math/cumnormdist.html

is an asymptotic expansion for the tail of the normal distribution,
which is easily converted to erf via the relation

erf(x) = 2 CND(x sqrt(2)) - 1

Rob Johnson <[email protected]>
take out the trash before replying