8

Pythonで二項確率を計算したい。私は式を適用しようとしました：

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

私が得る確率のいくつかは無限です。p=inf の値をいくつかチェックしました。そのうちの 1 つは、n=450,000 で k=17 です。この値は、float によって処理される最大値である 1e302 より大きくなければなりません。

私はそれから使用しようとしましたsum(np.random.binomial(n,p,numberOfTrials)==valueOfInterest)/numberOfTrials

これは numberOfTrials サンプルを描画し、値 valueOfInterest が描画される平均回数を計算します。

これは無限の値を上げません。しかし、これは続行する有効な方法ですか？そして、確率を計算するのに、なぜこの方法は無限の値を上げないのでしょうか?

4

4 に答える 4

7

対数を使用してすべての計算を行う必要があります。

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))
于 2014-03-05T15:46:27.227 に答える
7

対数ドメインで作業して、組み合わせ関数と累乗関数を計算し、それらを指数に上げます。

このようなもの：

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.4np.exp

于 2014-03-05T15:43:09.930 に答える
1

ゼロのような多重度を無限大のように回避するには、このように段階的な乗算を使用します。

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
于 2015-02-08T16:33:34.497 に答える