精度scipy.special.gammaln公司

Question

问题

我的许多编程都涉及scipy.stats公司. 一个新的问题需要计算beta-binomial分布的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。我的问题是这个函数将为不同的输入返回相同的值。在

示例

^{pr2}$

结果数组是

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

你的问题与loss of floating point precision in subtractions有关。这实际上并不取决于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中：

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

In [45]: a, b = symbols('a, b')

In [46]: (lnstirling3(a + b) - lnstirling3(a)).series(a, oo, 4)

 4    3    2      3    2         2                              
b    b    b      b    b    b    b    b                          
── - ── + ──   - ── + ── - ──   ── - ─                          
12   6    12     6    4    12   2    2        ⎛1⎞    ⎛1        ⎞
──────────── + ────────────── + ────── - b⋅log⎜─⎟ + O⎜──; a → ∞⎟
      3               2           a           ⎝a⎠    ⎜ 4       ⎟
     a               a                               ⎝a        ⎠

类似的渐近公式可以用类似的方式推导出pmf，当参数值较大时，可以用它们代替通常的表达式。在

编辑：如果您感到懒惰，可以将原始公式与mpmath一起使用，并通过mpmath.mp.dps打开更高的精度。但是，在求和之前，一定要先将k，n，k，n转换成mpmath.mpf。在