计算非常低的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

回答

相关问题