scipy stats binom cdf возвращает nan

Question

Если я правильно понимаю, cdf для дискретного распределения scipy.stats должно возвращать сумму вероятностей значений вплоть до заданного параметра.

Таким образом, scipy.stats.binom(7000000000, 0.5).cdf(6999999999) должен возвращать что-то почти точно 1, потому что в 7 миллиардах испытаний с вероятностью 50/50 вероятность успеха на 7 миллиардах минус 1 из них или меньше почти наверняка. Вместо этого я получаю np.nan. Фактически, за любое значение, предоставленное .cdf ЗА ИСКЛЮЧЕНИЕМ 7 миллиардов (или более), я получаю обратно np.nan.

Что здесь происходит? Существуют ли ограничения на число, с которыми могут работать дистрибутивы scipy.stats, которых нет в документах?

Answer 1

TL; DR

Отсутствие точности с плавающей запятой при внутренних расчетах. Хотя scipy - это библиотека Python, ее ядро ​​написано на C и использует числовые типы C.

---

Позвольте мне показать вам пример:

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)}")

И вывод:

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

Причина заключается в формуле CDF для биномиального распределения (я не могу встраивать изображения, поэтому вот ссылка на вики: https://en.wikipedia.org/wiki/Binomial_distribution

Внутри scipy источников мы увидим ссылку на эту реализацию: http://www.netlib.org/cephes/doubldoc.html#bdtr

Глубоко внутри это включает деление на trials (incbet.c, line 375: ai = 1.0 / a; здесь это называется a, но nwm). И если ваш trials слишком велик, результат этого деления настолько мал, что когда мы добавляем это маленькое число к другому, а не к такому маленькому, оно на самом деле не меняется, потому что здесь нам не хватает точности с плавающей запятой (только 64 бита). до сих пор). Затем, после некоторой арифметики, мы пытаемся получить логарифм от числа, но оно равно нулю, поскольку оно не меняется, когда должно. И log(0) не определено, что равно np.nan.