winston@cubsvax.UUCP (Ken Winston) (03/02/86)
Here is a question for numerical analysis types:
The standard normal distribution function is
N(t) = 1/sqrt(2*pi) * integral(-infinity to t) [exp(-z^2/2)]*dz.
I have seen the following approximation a few different times now as a way
of finding N(t) on a computer. It seems to be accurate in the usual range, for
t<=0:
N(t) ~= x*exp(-t^2/2)*(.436184 - .120167*x + .937298*x^2)/sqrt(2*pi),
where x = 1/(1+.3327*|t|)
For t>0 use N(t) = 1-N(-t).
My ignorance of numerical analysis is vast (not just half-vast). Does anybody
know where this comes from? What is the general method that gives rise to this
series?
Ken Winston
...{cmcl2,rna,cubsvax}!wealth!kenjimc@ucla-cs.UUCP (03/04/86)
In article <443@cubsvax.UUCP> winston@cubsvax.UUCP (Ken Winston) writes: >The standard normal distribution function is > N(t) = 1/sqrt(2*pi) * integral(-infinity to t) [exp(-z^2/2)]*dz. >I have seen the following approximation... >...Does anybody know where this comes from? Abramowitz and Stegun, Handbook of Mathematical Functions, Dover, 1965 (repr. 1972), paperback. (Copy of the original published by National Bureau of Standards.) This book contains over 1000 pages of functions, series approx- imations, numerical methods, tables, etc.etc. and is well worth the not outrageous price. For an approximation to erf(x) similar to the one you gave, as well as many formulae elsewhere in the book, they credit: Hastings, Approximations for Digital Computers, Princeton Univ. Press (1955) -- James F. Carter (213) 206-1306 UCLA-SEASnet; 2567 Boelter Hall; 405 Hilgard Ave.; Los Angeles, CA 90024 UUCP:...!{ihnp4,ucbvax,{hao!cepu}}!ucla-cs!jimc ARPA:jimc@locus.UCLA.EDU
weemba@brahms.BERKELEY.EDU (Matthew P. Wiener) (03/05/86)
In article <443@cubsvax.UUCP> winston@cubsvax.UUCP (Ken Winston) writes: >The standard normal distribution function is > > N(t) = 1/sqrt(2*pi) * integral(-infinity to t) [exp(-z^2/2)]*dz. > >I have seen the following approximation a few different times now as a way >of finding N(t) on a computer. It seems to be accurate in the usual range, for >t<=0: > > N(t) ~= x*exp(-t^2/2)*(.436184 - .120167*x + .937298*x^2)/sqrt(2*pi), > where x = 1/(1+.3327*|t|) > >For t>0 use N(t) = 1-N(-t). > >Does anybody >know where this comes from? What is the general method that gives rise to this >series? This particular representation is due to C Hastings _Approximations for Digital Computers_. First, get the function to look as geometrically like a polynomial (or other easy to compute function as possible). Asymptotes should be factored out first, for example. Second, economize the approximation. This consists of fudging the coefficients so that the error is spread around more uniformly, usually by replacing the highest terms by known lower degree approximants. There is a whole theory on how to do this, starting with the Chebyshev polynomials. In the particular case, one starts with the asymptotic expansion 2 2 -t 1 1 1.3 1.3.5 N(t) = -------- e --- [ 1 - ----- + ------ - ----- + ... ] sqrt(pi) 2 t 2 2 2 2 3 (2 t ) (2 t ) (2 t ) easily derived by repeatedly integrating by parts. It diverges, but the error is less than the first omitted term. It can be used as is for x>5 or so. (Depends on how much accuracy you need.) A Taylor series, either for N(t) or N(t)/exp(-t^2) converges too slowly for t>.7 or so. Where Hastings came up with the change of variables x=1/(1 + c t) is beyond me. In terms of the geometry of the error function, the resulting change gives one a much more polynomic geometry, and therefore a much better chance at getting good approximation. One then has to economize the coefficients. Consider T_4(x)=8x^4-8x^2+1. Pretty innocuos looking, isn't it? Actually, it satisfies the important property that |T_4(x)|<=1 for |x|<=1. Thus x^4 is approximately x^2-1/8, with error at most 1/8, on the interval [-1,1]. In general, T_n(x) is a polynomial with leading term 2^(n-1)*x^n, and can be used to simulate x^n with a lower degree polynomial by committing an error of 2^(1-n). For other situations and intervals, similar methods have been found. Suppose for example, you want a good degree 7 approximation to sine(x) on some interval. One takes the Taylor series out to degree 11, commit a Taylor series error by stopping, and then commit two further errors by economizing out the x^9 and x^11 terms. These three errors together are much smaller than the degree 7 Taylor series error, and are more uniformly distributed. A very readable account is in F Acton _Numerical Methods That Work_, chapters 1 and 12. Not only is it a good book (although be warned that some of his *particular* comments about certain methods are no longer valid--old methods get revived a lot in numerical analyis), it is the only numerical analysis book that makes it clear that THINKING is the most important technique for solving numerical problems on computers. It is also the FUNNIEST math book ever written. "A favorite form of lunacy among engineers ..." "If your author were writing another typical book on Numerical Analysis ..." "reducing ... costly n-dimensional wanderings through foggy Numberland." "... you end up by getting angry and throwing the guy out of your office." ucbvax!brahms!weemba Matthew P Wiener/UCB Math Dept/Berkeley CA 94720
martin@entropy.UUCP (Don Martin) (03/05/86)
Ken Winston writes: > I have seen the following approximation a few different times now as a way ....... ( more of quote at end ) <> This looks like one of Hasting's approximations. It is not a series approximtion but rather an ( approximate ) Chebychev fit that minimizes the maximum error. Hastings published a strange little book in about 1956 on numerical approximations. It is strang in that it has almost no text. It seems to be slides from lectures. It is quite interesting. For a rundown on various approximations of the normal see: Statistical Computing by Kennedy and Gentile, pp.89-96. Marcel Dekker, Inc. For this methodology ( min-max approximations ) see : The Approximation of Elementary Functions ( ? may be inacurate ) Hart et al. John Wiley ( sponsered by SIAM) > > I have seen the following approximation a few different times now as a way > of finding N(t) on a computer. It seems to be accurate in the usual range, for > t<=0: > > N(t) ~= x*exp(-t^2/2)*(.436184 - .120167*x + .937298*x^2)/sqrt(2*pi), > where x = 1/(1+.3327*|t|) > > For t>0 use N(t) = 1-N(-t). > > My ignorance of numerical analysis is vast (not just half-vast). Does anybody > know where this comes from? What is the general method that gives rise to this > series? > > Ken Winston > ...{cmcl2,rna,cubsvax}!wealth!ken Donald C. Martin, phone (206) 543 1044 SC-32, Dept. of Biostatistics., U. of Wash, Seattle WA, 98195 {decvax,ihnp4,ucbvax!lbl-csam}!uw-beaver!entropy!martin