Вычисление биномиальной вероятности для огромных чисел

Я хочу вычислить биномиальные вероятности на питоне. Я попытался применить формулу:

probability = scipy.misc.comb(n,k)*(p**k)*((1-p)**(n-k))

Некоторые из вероятностей, которые я получаю, бесконечны. Я проверил некоторые значения, для которых p=inf. Для одного из них n=450 000 и k=17. Это значение должно быть больше 1e302, что является максимальным значением, обрабатываемым числами с плавающей запятой.

Затем я попытался использовать sum(np.random.binomial(n,p,numberOfTrials)==valueOfInterest)/numberOfTrials

Это отрисовывает выборки numberOfTrials и вычисляет среднее количество раз, когда значение valueOfInterest отрисовывается.

Это не увеличивает бесконечное значение. Однако является ли это допустимым способом продолжения? И почему этот способ не даст никакого бесконечного значения, в то время как вычисление вероятностей дает?


python numpy algorithm probability binomial-coefficients
person bigTree    schedule 05.03.2014    source источник

Ответы (4)


arrow_upward
7
arrow_downward

Работайте в логарифмической области, чтобы вычислить функции комбинации и возведения в степень, а затем возведите их в степень.

Что-то вроде этого:

combination_num = range(k+1, n+1)
combination_den = range(1, n-k+1)
combination_log = np.log(combination_num).sum() - np.log(combination_den).sum()
p_k_log = k * np.log(p)
neg_p_K_log = (n - k) * np.log(1 - p)
p_log = combination_log + p_k_log + neg_p_K_log
probability = np.exp(p_log)

Избавляет от числового недополнения/переполнения из-за больших чисел. В вашем примере с n=450000 и p = 0.5, k = 17 он возвращает p_log = -311728.4, т.е. т. е. журнал окончательной вероятности довольно мал, и, следовательно, при взятии np.exp происходит недополнение. Однако вы все еще можете работать с логарифмической вероятностью.

person Sudeep Juvekar    schedule 05.03.2014

arrow_upward
9
arrow_downward

Поскольку вы используете scipy, я подумал, что упомяну, что в scipy уже реализованы статистические распределения. Также обратите внимание, что при таком большом n биномиальное распределение хорошо аппроксимируется нормальным распределением (или пуассоновским, если p очень мало).

n = 450000
p = .5
k = np.array([17., 225000, 226000])

b = scipy.stats.binom(n, p)
print b.pmf(k)
# array([  0.00000000e+00,   1.18941527e-03,   1.39679862e-05])
n = scipy.stats.norm(n*p, np.sqrt(n*p*(1-p)))
print n.pdf(k)
# array([  0.00000000e+00,   1.18941608e-03,   1.39680605e-05])

print b.pmf(k) - n.pdf(k)
# array([  0.00000000e+00,  -8.10313274e-10,  -7.43085142e-11])
person Bi Rico    schedule 05.03.2014

arrow_upward
7
arrow_downward

Я думаю, вы должны делать все свои вычисления, используя логарифмы:

from scipy import special, exp, log
lgam = special.gammaln

def binomial(n, k, p):
    return exp(lgam(n+1) - lgam(n-k+1) - lgam(k+1) + k*log(p) + (n-k)*log(1.-p))
person hivert    schedule 05.03.2014
comment
Обратите также внимание на scipy.special функцию xlogy, которая более стабильна, чем, например, k*log(p). - person Ian Hincks; 18.08.2017

arrow_upward
1
arrow_downward

Чтобы избежать кратности, такой как ноль на бесконечность, используйте пошаговое умножение, как это.

def Pbinom(N,p,k):
    q=1-p
    lt1=[q]*(N-k)
    gt1=list(map(lambda x: p*(N-k+x)/x, range(1,k+1)))
    Pb=1.0
    while (len(lt1) + len(gt1)) > 0:
        if Pb>1:
            if len(lt1)>0:
                Pb*=lt1.pop()
            else:
                if len(gt1)>0:
                    Pb*=gt1.pop()
        else:
            if len(gt1)>0:
                Pb*=gt1.pop()
            else:
                if len(lt1)>0:
                    Pb*=lt1.pop()
    return Pb
person Алексей Патрашов    schedule 08.02.2015