Definite Integral of e^(-x^2) and e^((-x^2)/2))

Date: 08/18/99 at 10:58:10
From: Mark W. Miller
Subject: Integral of e^(-x^2) or e^((-x^2)/2))

Dear Dr. Math,

I have tried every other venue for learning the answer to my question, 
and now I'm turning to you. I wish to create the cumulative normal 
frequency distribution found in the back of virtually every 
introductory statistics book.

I know from your Web site and my own calculus text (Ellis and Gulick, 
1982) that the integral of e^(-x^2) from 0 to infinity is 0.5*sqrt(pi) 
and the integral of e^((-x^2)/2)) from -infinity to +infinity is 
sqrt(2*pi). However, I cannot seem to take the next step from there to 
integrate these functions from, for example, x = 0 to 1.96, which 
should be 0.4750 for the standard normal curve (Snedecor and Cochran, 
1989).

Stigler (_The History of Statistics_, 1986 p. 134) gives two different 
solutions for e^(-x^2); one for 0 to x, and one for x to +infinity. 
These solutions were derived by LaPlace and neither seems to create 
the table I want. I don't know what I am doing wrong.

I am not a mathematician, and am only a self-taught statistician with 
degrees in ecology. Any advice you can give would be appreciated, even 
if it is just a reference to a mathematical statistics textbook for 
upper level undergraduate statistics majors.

Sincerely,

Mark W. Miller

Date: 08/18/99 at 14:25:24
From: Doctor Anthony
Subject: Re: integral of e**(-x**2) or e**((-x**2)/2))

The standard normal integral is

      = 1/sqrt(2pi).INT[e^(-x^2/2).dx]

There is no closed form of the integral, except between limits 0 and 
infinity, so we need to find an approximate integral for the range of, 
say x = -1 to x = +1.

We can expand e^(-x^2/2) to a sufficient number of terms to get the 
accuracy we require.

     e^(-x^2)2 = 1 + (-x^2/2) + (-x^2/2)^2/2! + (-x^2/2)^3/3! + ...

               = 1 - x^2/2 + x^4/8 - x^6/48 + x^8/384 - ...

and integrating term by term we get

     x - x^3/6 + x^5/40 - x^7/336 + x^9/3456 - x^11/42240 + ....

Putting in the limits 0 and 1 (we double the result to get -1 to +1)

     = 1 - 1/6 + 1/40 - 1/336 + 1/3456 - 1/42240 + ...

     = 0.85562282

and doubling this gives  1.711245641.

Finally divide by sqrt(2pi) to get 0.682688238.

So the area under the standard normal curve from x = -1 to x = +1 is

     0.682688

Clearly you can use this method for integration between any two values 
of x.

- Doctor Anthony, The Math Forum
  http://mathforum.org/dr.math/   

Date: 01/10/2002 at 21:06:04
From: Nick de Peyster
Subject: Formula for Cumulative Standard Normal Distribution 

When I tested the above approach against the NORMSDIST function in 
Excel, the answers agreed pretty well until Z > 1 and Z <- 1. At 
points farther than 1 from Z = 0, the differences between Excel and 
the formula above became very large.

In order to improve the accuracy of the formula, do I need to extend 
the series to include more terms, or is there a simpler formula 
altogether?

Thanks,
Nick

Date: 01/10/2002 at 21:37:53
From: Doctor Schwa
Subject: Re: Formula for Cumulative Standard Normal Distribution 

Hi Nick,

A pretty simple approximation, using one polynomial and one 
exponentiation, is at

   Normal Distribution approximations - Bernt Arne Odegaard
   http://finance.bi.no/~bernt/gcc_prog/algoritms/algoritms/node59.html   

Scroll down the page to where you see "cumulative normal
distribution" and you get


// cumulative univariate normal distribution.
// This is a numerical approximation to the normal distribution.  
// See Abramowitz and Stegun: Handbook of Mathemathical functions
// for description.  The arguments to the functions are assumed 
// normalized to a (0,1 ) distribution. 

double N(double z) {
    double b1 =  0.31938153; 
    double b2 = -0.356563782; 
    double b3 =  1.781477937;
    double b4 = -1.821255978;
    double b5 =  1.330274429; 
    double p  =  0.2316419; 
    double c2 =  0.3989423; 

    if (z >  6.0) { return 1.0; }; // this guards against overflow 
    if (z < -6.0) { return 0.0; };
    double a=fabs(z); 
    double t = 1.0/(1.0+a*p); 
    double b = c2*exp((-z)*(z/2.0)); 
    double n = ((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t; 
    n = 1.0-b*n; 
    if ( z < 0.0 ) n = 1.0 - n; 
    return n; 
}

I just tried it out on Excel and it seemed to give the right answer to 
at least 7 decimal places for all values of z.

- Doctor Schwa, The Math Forum
  http://mathforum.org/dr.math/