1

統合によって定義される分布の尤度関数を作成しようとしています。integrate() 関数を使用していますが、関数の残りの部分でこれを使用しようとすると、エラーが発生します。

「B(alpha + i, beta + 6 - i)/B(alpha, beta) のエラー: 二項演算子への非数値引数」

積分の値は、たとえば、「絶対誤差 < 0.00078 で 9.501501」です。trunc() を使用しようとしましたが、これも役に立ちません。私はRに比較的慣れていないので、これに対する既知の解決策はありますか? どんな助けでも大歓迎です！

B <- function(a,b){ 
   integrand <- function(t){(t^(a-1))*((1-t)^(b-1))} 
   integrate(integrand,lower=0,upper=1) 
} 
betalik <- function(alpha,beta){ 
    likelihood <- 0 Z <- c(37,22,25,29,34,49) 
    for(i in 1:6) 
       likelihood <- likelihood + 
         Z[i]*log((B(alpha+i,beta+6-i))/B(alpha,beta)) 
    return(likelihood)

}

ドリアン、

4

2 に答える 2

4

上記のコメントと回答のフォローアップ...

結果$valueからコンポーネントを抽出することに関する@Fhnuzoagのコメントを統合して、より適切にフォーマットされた元の関数を次に示します。integrate()

B <- function(a,b){
    integrand <- function(t){(t^(a-1))*((1-t)^(b-1))}
    integrate(integrand,lower=0,upper=1)$value
}
betalik <- function(alpha,beta){
    likelihood <- 0
    Z <- c(37,22,25,29,34,49)
    for(i in 1:6) likelihood <- likelihood +
        Z[i]*log((B(alpha+i,beta+6-i))/B(alpha,beta))
    return(likelihood)
}

ここで、@ dasonのコメントをチェックして、あなたのB関数がRの組み込み関数と同等であるということを確認しbetaます（ただし、Rの関数は確かに高速で、おそらくより正確です）。

all.equal(B(1.1,2.7),beta(1.1,2.7))  ## TRUE

'data'を個別に指定することを好みます。

Z <- c(37,22,25,29,34,49)

lbeta組み込み（log-beta）関数を使用し、ベクトル化された尤度関数の新しいバージョン：

blik2 <- function(alpha,beta,Z) {
    index <- 1:6
    sum(Z*(lbeta(alpha+index,beta+6-index)-lbeta(alpha,beta)))
}

all.equal(blik2(1.1,2.7,Z),betalik(1.1,2.7))
## small difference, blik2 is *probably* more
##   accurate ... "Mean relative difference: 6.406495e-08"

（これをもう少し一般化6して、コード内のの値をlength(Z)...に置き換えるとよいでしょう。）

于 2012-05-24T18:02:03.823 に答える
4

?integrate....の統合の例を使用します。

integrate(dnorm, -1.96, 1.96)
# 0.9500042 with absolute error < 1e-11

ここで実際に何が起こっているのですか？さて、統合関数は「統合」クラス S3 オブジェクトを作成します。これは基本的に、いくつかのフィールドを持つリストであり、次にそれに適用さprint()れ、順番に実行されますprint.integrate()。関数の出力は、実際には文字列「絶対誤差 < 1e-11 で 0.9500042」ではなく、単にそのように表示されるだけです。

によって作成された R オブジェクトがintegrate()実際に何であるかを取得するには、次のようにします。

obj = integrate(dnorm, -1.96, 1.96)
str(obj)
# List of 5
# $ value       : num 0.95
# $ abs.error   : num 1.05e-11
# $ subdivisions: int 1
# $ message     : chr "OK"
# $ call        : language integrate(f = dnorm, lower = -1.96, upper = 1.96)
# - attr(*, "class")= chr "integrate"

したがって、積分の値だけが必要な場合はvalue、計算のために、その関数によって作成されたリストのフィールドを抽出する必要があります。例えば

10*integrate(dnorm, -1.96, 1.96)$value
# [1] 9.500042
于 2012-05-24T17:53:11.743 に答える