scipy.special.gammaln 的精度

Question

问题

我的很多编程都涉及到 scipy.stats 中的统计函数。一个新问题需要计算beta-二项式分布的 pmf 。因为它具有解析形式，但没有出现在 scipy.stats 中，所以我需要自己为其 pmf 定义一个函数。我正在使用 scipy 版本 0.12.0 和 numpy 版本 1.7.0。

import numpy
from scipy.special import gammaln, betaln

def beta_binomial_pmf(k, n, K, N):
    # compute natural log of pmf
    ln_pmf = ( gammaln(n+1) - gammaln(k+1) - gammaln(n-k+1) ) + \
        - betaln(K+1,N-K+1) + betaln(K+k+1,N-K+n-k+1)
    return numpy.exp(ln_pmf)

在统计问题中，我试图解决 n 和 k 的值通常在 0 到 100 之间的范围内，但 K 和 N 可以大到 1e9。我的问题是这个函数将为不同的输入返回相同的值。

例子

k = 0
n = 5
K = numpy.array([12, 10, 8])
N = 101677958
beta_binomial(k, n, L, N)

结果数组是

array([ 0.99999928,  0.99999905,  0.99999928])

鉴于 K 的每个值都不同，这很奇怪。更好地了解数组中第一个和第三个值之间的相似性

1 - beta_binomial(k, n, L, N)
array([  7.15255482e-07,   9.53673862e-07,   7.15255482e-07])

---

gammaln函数精度的一个非常简单的测试是 1-(Gamma(N+1)/Gamma(N))/N。它很有用，因为如果您在纸上计算代数，结果正好是 0。

N = numpy.logspace(0,10,11)
1-numpy.exp(gammaln(N+1)-gammaln(N))/N
array([  0.00000000e+00,  -1.11022302e-15,   1.90958360e-14,
    -9.94537785e-13,  -4.96402919e-12,   7.74684761e-11,
    -1.70086167e-13,   1.45905219e-08,   2.21033640e-07,
    -7.64616381e-07,   2.54126535e-06])

题

我认识到可以计算的精度是有限的，但是在 N=1e7 附近会发生什么，使精度变化gammaln五个数量级？有关如何解决此问题的建议？

Answer 1

您的问题与减法中浮点精度的损失有关。这实际上并不取决于 Scipy 的 gammaln 和 betaln 的精度。问题在于，对于大 N，gammaln(N+1) 与 gammaln(N) 的数量级相同，但比 gammaln(N+1)-gammaln(N) 大得多。因此，当您计算差异时，您会损失 ~ log10(gammaln(N)) 位精度。这是浮点的普遍问题。

您可以通过渐近展开来解决这个问题（参见betaln implementation，它必须处理相同的问题）。即，您可以对 Gamma(a + b) - Gamma(a) 使用扩展 >> |b|, 1。在 Sympy 中：

在 [44]: def lnstirling3(z): return (z - sympify('1/2')) * log(z) - z + log(sqrt(2*pi)) + 1/(12*z) - 1/（360*z*z*z）

在 [45] 中：a, b = symbols('a, b')

在 [46] 中：(lnstirling3(a + b) - lnstirling3(a)).series(a, oo, 4)

 4 3 2 3 2 2                              
bbbbbbbb                          
?? - ?? + ?? - ?? + ?? - ?? ?? - ?                          
12 6 12 6 4 12 2 2 ?1? ?1 ?
???????????? + ???????????????? + ?????? - 博客？？？+ O???; 一种 ？??
      3 2 a?a? ? 4 ?
     啊？

类似的渐近公式可以以类似的方式为您的 pmf 导出，并且当参数具有大值时，它们可以代替通常的表达式使用。

编辑：如果您感到懒惰，可以将原始公式与mpmath一起使用，并通过mpmath.mp.dps. mpmath.mpf但是，在对它们求和之前，请务必先将k, n, K, N 转换为第一个。