scipy stats binom cdf returns nan

user_2886575 0
Viewed 1079

If I understand correctly, the cdf for a scipy.stats discrete distribution should return the sum of the probabilities of the values up to the parameter given.

Thus, scipy.stats.binom(7000000000, 0.5).cdf(6999999999) should return something almost exactly 1, because in 7 billion trials, with a 50/50 chance, the probability of getting a success on 7 billion minus 1 of them or less is pretty much certain. Instead, I get np.nan. In fact, for any value provided to .cdf EXCEPT 7 billion itself (or more), I get back np.nan.

What is going on here? Is there some limit to the numbers that scipy.stats distributions can handle that is not in the docs?

scipy python
1 Answers

TL; DR

Lack of floating point precision during internal calculations. Though scipy is a Python library, it's core is written on C and uses C numeric types.

Let me show you an example:

import scipy.stats

for i in range (13):
    trials = 10 ** i
    print(f"i: {i}\tprobability: {scipy.stats.binom(trials, 0.5).cdf(trials - 1)}")

And output is:

i: 0    probability: 0.5
i: 1    probability: 0.9990234375
i: 2    probability: 0.9999999999999999
i: 3    probability: 0.9999999999999999
i: 4    probability: 0.9999999999999999
i: 5    probability: 0.9999999999999999
i: 6    probability: 0.9999999999999999
i: 7    probability: 0.9999999999999999
i: 8    probability: 0.9999999999999999
i: 9    probability: 0.9999999999999999
i: 10   probability: nan
i: 11   probability: nan
i: 12   probability: nan

The reason lies in the CDF formula for binomial distribution (I can't embed images so here is link to wiki: https://en.wikipedia.org/wiki/Binomial_distribution

Inside the scipy sources we would see a refence to this implementation: http://www.netlib.org/cephes/doubldoc.html#bdtr

Deeply inside it involves division by trials (incbet.c, line 375: ai = 1.0 / a; here it's called a, but nwm). And if your trials is too large, the result of this division is so small that when we add this little number to another, not so little number, it doesn't actually change because we lack floating point precision here (only 64 bits so far). Then, after some more arithmetics, we try get logarithm from a number, but it's equal to zero as it hasn't changes when it should. And log(0) is not defined, which equals to np.nan.

user_78058770
Related