1

R の次の関数に問題があります。

test <- function(alpha, beta, n){
  result <- exp(lgamma(alpha) + lgamma(n + beta) - lgamma(alpha + beta + n) - (lgamma(alpha) + lgamma(beta) - lgamma(alpha + beta)))
  return(result)
}

次の値を挿入すると、次のようになります。

betabinom(-0.03292708, -0.3336882, 10)

失敗し、NaN. これは、Excel で正確な関数を実装すると、数値ではない結果が得られるためです。J32 は のセルalpha、K32betaと L32 は のセルであるため、Excel での実装は簡単ですN。結果のセルの実装を以下に示します。

=EXP(GAMMALN(J32)+GAMMALN(L32+K32)-GAMMALN(J32+K32+L32)-(GAMMALN(J32)+GAMMALN(K32)-GAMMALN(J32+K32)))

したがって、関数はゼロより大きいアルファとベータ、およびゼロ以上のnに対してのみ定義されているため、これは正しい答えを与えるようです。したがって、ここで何が起こっているのだろうか？パッケージ Rmpf を試して数値精度を上げましたが、何もしていないようです。

ありがとう

4

1 に答える 1

2

tl;dr log(gamma(x)) は、あなたが考えているよりも、または Excel が考えているよりも一般的に定義されています。alpha関数がandの負の値を受け入れないようにしたい場合beta、または を返すようにしたい場合はNaN、手動でテストして適切な値 ( if (alpha<0 || beta<0) return(NaN)) を返すだけです。

これは数値の精度の問題ではなく、定義の問題です。ガンマ関数は、負の実数値に対して定義されています?lgamma。

ガンマ関数は (Abramowitz and Stegun セクション 6.1.1、255 ページ) によって定義されます。

ガンマ(x) = 整数_0^Inf t^(x-1) exp(-t) dt

ゼロおよび負の整数を除くすべての実数 'x' ('NaN' が返される場合)。

さらに、参考までにlgamma・・・

...そしてガンマ関数の絶対値の自然対数...

(原文強調)

curve(lgamma(x),-1,1)

ここに画像の説明を入力

gamma(-0.1)          ## -10.68629
log(gamma(-0.1)+0i)  ## 2.368961+3.141593i
log(abs(gamma(-0.1)) ## 2.368961
lgamma(-0.1)         ## 2.368961

Wolfram Alpha は 2 番目の計算に同意します。

于 2015-08-06T12:49:55.000 に答える