> 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