4

问题

我的很多编程都涉及 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了五个数量级？关于如何解决这个问题的建议？

4

1 回答 1

5

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

您可以通过渐近扩展来解决此问题（参见betaln implementation，它必须处理相同的问题）。也就是说，您可以使用 Gamma(a + b) 的扩展 - Gamma(a) 用于 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 = 符号（'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 ⎞
──────────── + ────────────── + ────── - b⋅log⎜─⎟ + O⎜──；a → ∞⎟
      3 2 ⎝a⎠ ⎜ 4 ⎟
     aa⎝a⎠

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

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

于 2014-01-20T10:27:17.973 回答