4

問題

私のプログラミングの多くは、scipy.stats の統計関数に関係しています。新しい問題では、ベータ二項分布の 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 番目と 3 番目の値の類似性をよりよく理解するには

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

関数の精度の非常に簡単なテストは、gammaln1-(ガンマ(N+1)/ガンマ(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])

質問

gammaln計算できる精度には限界があることは認識していますが、精度が 5 桁変化する N=1e7 付近で何が起こるのでしょうか? この問題を回避する方法についての提案はありますか?

4

1 に答える 1

5

あなたの問題は、減算で浮動小数点精度が失われることに関係しています。実際、これは Scipy の gammaln と betaln の精度には依存しません。問題は、N が大きい場合、gammaln(N+1) は gammaln(N) と同じ桁ですが、gammaln(N+1)-gammaln(N) よりもはるかに大きいことです。結果として、差を計算すると、~ log10(gammaln(N)) 桁の精度が失われます。これは、浮動小数点の一般的な問題です。

これは、漸近展開によって回避できます (同じ問題に対処する必要があるbetaln 実装を参照)。つまり、Gamma(a + b) - Gamma(a) for 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 = symbol('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⎜──; あ → ∞⎟
      3 2 ⎝a⎠ ⎜ 4 ⎟
     ああ⎝⎠

同様の方法で、pmf に対して同様の漸近式を導き出すことができ、パラメーターの値が大きい場合に通常の式の代わりに使用できます。

EDIT : 怠惰な場合は、元の式をmpmathと一緒に使用して、 を介してより高い精度を有効にすることができますmpmath.mp.dps。mpmath.mpfただし、それらを合計する前に、必ず最初に k、n、K、N をキャストしてください。

于 2014-01-20T10:27:17.973 に答える