Siginificant figures calculating zprob

[Josiah Carlson]

> I want to calculate zprob(the area under the normal curve)
> with python and I managed to find a function from the
> internet.
> 
> But the problem is that the function calculates the result
> with only a few significant figures. If I want to get the 
> 20th number of the result(z-prob) what should I do? I want
> to get the confidence level of "6sigma" and all I get at the
> moment is "1".
> 
> I remember that Python's long type has unlimited number of
> significant figures as long as the memory allows. What about
> floating point numbers?

Python floats are normally IEEE 754 FP values, which generally have 53 
bits of precision.

It's been a while since I did this kind of stuff, but after doing a few 
tests, the summation using Simpson's rule for numerical integration 
converges fairly quickly for most znorm ranges.  I use 4096 divisions 
beacuse it seems to be fast and accurate to at least 10 decimal places 
for the ranges I've tested.

More divisions may or may not increase precision, but for most tasks, I 
would imagine the code given at the end would be sufficient.  There are 
various float casts, and floating point constants given in the code. 
With the tests that I did, it seems that some are necessary for higher 
precision.

Oh, I hope this is precise enough...
 >>> div = 1024
 >>> while div < 16*1024:
...     print div, repr(znorm_range(0, 1, div))
...     div *= 2
...
1024 0.3413447460685452
2048 0.3413447460685422
4096 0.34134474606854304
8192 0.34134474606854265
 >>>

  - Josiah



from math import sqrt, e, pi

def znorm_x(x,c=1.):
     return c/sqrt(2.*pi)*e**(-(x*x)/2.)

def znorm_range(low, high, divisions=4096):
     #divisions needs to be even, we'll increment it if necessary
     divisions += divisions%2
     #necessary if either low or high are integers
     low, high = float(low), float(high)
     inc = (high-low)/float(divisions)
     x = low+inc
     t = znorm_x(low) + znorm_x(high)
     c = 4.
     while x < high:
         t += znorm_x(x,c)
         x += inc
         c = 6.-c
     return (high-low)*t/3./divisions