計算龐大數字的二項式概率

Question

我想計算Python上的二項式概率。我嘗試應用公式：

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

我得到的一些概率是無限的。我檢查了一些p = inf的值。其中之一，n = 45萬和k = 17。該值必須大於1e302，這是浮點處理的最大值。

然後我試圖用sum(np.random.binomial(n,p,numberOfTrials)==valueOfInterest)/numberOfTrials

這吸引numberOfTrials樣本和計算的時間值valueOfInterest繪製的平均數量。

這不會產生任何無限的價值。但是，這是一種有效的方式嗎？爲什麼這種方式不會提高任何無限的價值，而計算概率呢？

Answer 1

在日誌域中工作以計算組合和指數函數，然後將它們提升爲指數。

事情是這樣的：

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時發生下溢。但是，您仍然可以使用日誌概率。

Answer 2

我的事情，你應該做你的計算採用對數：

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

Answer 3

因爲你使用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]) 

Answer 4

爲了避免像無窮大這樣的多重性使用逐步乘法這樣。

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