(Note: I know this looks like a programming question, but I'm OK with the programming part and just want to understand the mathematics.)
I found a bit of code to calculate the integral of the normal bell curve that I wrote 2 years ago, and I don't remember why or how it works. I just have a note that it's based on "Abramowitz & Stegun (1964)".
I want to verify that it's correct, and find out how accurate it is. I tried doing some internet searches but I didn't find anything that I could follow. (I haven't finished integral calculus yet.)
Basically I calculate standard deviation as $\mu$, define $p(\mu, x) = \frac{e^{-\frac{x^2}{2\mu^2}}}{\sqrt{2\pi\mu^2}} $
From there, the function $f(\mu, x)$, specifying the probability that a random sample is less than $x$, is calculated as:
$ t= \frac{1}{1+0.2316419 \frac{x}{\mu}} $
$ f(\mu, x) = 1 - p(\mu, x) (0.319381530t -0.356563782t^2 + 1.781477937t^3 -1.821255978t^4 + 1.330274429t^5) $
So, is this correct? Also, how accurate is it?
For reference, the original code is:
static double probability_density(double stddev, double x) {
return exp(-x*x/(2*stddev*stddev)) / (sqrt(2.0*Pi*stddev*stddev));
}
static double normal_prob_lower(double stddev, double x) {
//Cumulative Distribution Function
//find probability that, in a normal distribution, a random value is lower than x
//Abramowitz & Stegun (1964) approximation
x/=stddev;
stddev=1;
bool rev=0;
if (x<0) {
x=-x;
rev=1;
}
const double b[]= {0.2316419,0.319381530,-0.356563782,1.781477937,-1.821255978,1.330274429};
double t=1.0/(1.0+b[0]*x);
double sum=0;
double tm=1;
for (int i=1;i<=5;i++) {
tm*=t;
sum+=b[i]*tm;
}
double prob = 1-sum*probability_density(stddev, x);
if (rev) {
prob=1-prob;
}
return prob;
}