計算非常低的p值的Python

我已經達到了我的代碼中的一個點，我反覆計算的某些條件下，p值：計算非常低的p值的Python

from scipy.stats import hypergeom 
pval = min(hypergeom.sf(k, M, n, N) + hypergeom.pmf(k, M, n, N), 1)

這種方法適用於「小」 N的（在彈出許多成功的要素）。我試圖高達500

後，我與n=5000嘗試，我得到一個精確的錯誤，因爲計算的P值非常低，四捨五入爲0

我怎樣才能克服在Python這些precision errors？

來源

2017-06-19 Jack

https://docs.python.org/2/tutorial/floatingpoint.html德在這裏查看浮點數的問題和侷限性。 –

在相關說明中，您是通過蟒蛇漂浮物還是numpy漂浮物？ –

用於'k'，'M'和'N'的典型值是什麼？ –

您要計算的值小於可以使用64位浮點值表示的值。你在評論中給出的一個例子是k = 5007, M = 45956, n = 18969, N = 5267。對於M，n和N那些值時，下溢PMF爲0時k參數是3478：

In [46]: k = 5007 

In [47]: M = 45956 

In [48]: n = 18969 

In [49]: N = 5267 

In [50]: hypergeom.pmf(3476, M, n, N) 
Out[50]: 9.8813129168249309e-324 

In [51]: hypergeom.pmf(3477, M, n, N) 
Out[51]: 4.9406564584124654e-324 

In [52]: hypergeom.pmf(3478, M, n, N) 
Out[52]: 0.0

的標準方法來解決這一問題是與概率的對數工作。該SciPy的離散分佈具備的功能logpmf和logsf此：

In [53]: hypergeom.logpmf(3476, M, n, N) 
Out[53]: -743.80749253381509 

In [54]: hypergeom.logpmf(3477, M, n, N) 
Out[54]: -744.95722489454783 

In [55]: hypergeom.logpmf(3478, M, n, N) 
Out[55]: -746.10790755529888 

In [56]: hypergeom.logpmf(5007, M, n, N) 
Out[56]: -3952.1782915849763

爲了計算hypergeom.sf(k, M, n, N) + hypergeom.pmf(k, M, n, N)，您可以使用numpy.logaddexp：

In [58]: np.logaddexp(hypergeom.logsf(k, M, n, N), hypergeom.logpmf(k, M, n, N)) 
Out[58]: -3952.1508002445375

唯一不方便的是，進一步的計算和比較，必須立足於概率的對數。如果這不適用於您，則必須切換到提供更高精度浮點計算的庫（例如mpmath）。例如，以下功能使用mpmath計算PMF和生存函數：

def hypergeom_pmf(k, M, n, N): 
    tot, good = M, n 
    bad = tot - good 
    pmf = (mpmath.beta(good+1, 1) * mpmath.beta(bad+1,1) * mpmath.beta(tot-N+1, N+1)/
      (mpmath.beta(k+1, good-k+1) * mpmath.beta(N-k+1,bad-N+k+1) * mpmath.beta(tot+1, 1))) 
    return pmf 

def hypergeom_sf(k, M, n, N): 
    sf = (mpmath.binomial(N, k+1) * mpmath.binomial(M-N, n - k - 1)/mpmath.binomial(M, n) * 
      mpmath.hyp3f2(1, k + 1 - n, k + 1 - N, k + 2, M + k + 2 - n - N, 1)) 
    return sf

（在hypergeom_pmf(k, M, n, N)使用的表達式scipy.stats.hypergeom._logpmf從SciPy的的實現採取hypergeom_sf使用對the wikipedia page on the hypergeometric distribution給出的CDF式它。不一定是生存功能的最佳實現）

例如：

In [107]: import mpmath 

In [108]: mpmath.mp.dps = 40 

In [109]: k, M, n, N 
Out[109]: (5007, 45956, 18969, 5267) 

In [110]: hypergeom_pmf(k, M, n, N) 
Out[110]: mpf('3.897413335837289136238051958307757561884655e-1717') 

In [111]: hypergeom_sf(k, M, n, N) 
Out[111]: mpf('1.086314878026431217760059547783856962636701e-1718')

來源

2017-06-19 20:41:03

如果問題不是太多，你能否提供一個'mpmath'的快速示例？ – Jack

@Jack：我添加了一個使用'mpmath'來計算PMF和生存函數的例子。 –

計算非常低的p值的Python

回答

相關問題