[SciPy-User] Loss of precision with powers and factorials?

[Chris Cameron]

Hi,

I tried to do some work with the binomial probability formula, but I’m not getting the answers I expect. My code:

############
import numpy as np
import scipy as sp

n = 36000
i = np.array(range(0, 250))

np.sum(
    (sp.misc.factorial(n) / \
    (sp.misc.factorial(n-i) * sp.misc.factorial(i))) * \
    np.power(0.00625, i) * np.power((1-0.00625), (n-i)))
############

When I execute this I get ‘nan’.

In Mathematica this same formula and values get me a seemingly reasonable answer:
Sum[36000!/((36000 - i)!*i!)*0.00625^i*(1 - 0.00625)^(36000 - i), {i, 0, 250}] = 0.95406

I’m running Python 2.7.11, Numpy 1.10.4, and SciPy 0.17 from inside iPython 4.1.2


From trying to execute this piece by piece, it seems like the “np.power()” function is just returning “0” when the result is super small. Could the problem be not enough precision is being held?


Thanks!

Chris