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

[josef.pktd at gmail.com]

On Thu, May 5, 2016 at 4:43 PM, Chris Cameron <chris at upnix.com> wrote:
> 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?


Floating precision does not work well enough for this.

The source of scipy.stats.distribution is now for most parts a good
place to see how to preserve numerical precision.

Essentially, don't do that, work in log space and use scipy special
function like lngamma instead of factorial and similar.

Josef


>
>
> Thanks!
>
> Chris
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> https://mail.scipy.org/mailman/listinfo/scipy-user