0

R で optim() を使用して積分を含む可能性を解き、2 つのパラメーターのヘッセ行列を取得するのに問題があります。アルゴリズムは収束しましたが、optim() で hessian=TRUE オプションを使用するとエラーが発生します。エラーは次のとおりです。

integrate(integrand1, lower = s1[i] - 1, upper = s1[i]) のエラー: 非有限関数値

NAの警告メッセージもありました

これが私のコードです：

s1=c(1384,1,1219,1597,2106,145,87,1535,290,1752,265,588,1188,160,745,237,479,39,99,56,1503,158,916,651,1064,166,635,19,553,51,79,155,85,1196,142,108,325  
 ,135,28,422,1032,1018,128,787,1704,307,854,6,896,902)

LLL=function (par) {

  integrand1 <- function(x){ (x-s1[i]+1)*dgamma(x, shape=par[1], rate=par[2]) }
  integrand2 <- function(x){ (-x+s1[i]+1)*dgamma(x, shape=par[1],rate=par[2]) }



  likelihood = vector() 

  for(i in 1:length(s1)) {likelihood[i] = 
    log( integrate(integrand1,lower=s1[i]-1,upper=s1[i])$value+ integrate(integrand2,lower=s1[i],upper=s1[i]+1)$value )  
  }

  like= -sum(likelihood)
  return(like)

}




optim(par=c(0.1,0.1),LLL,method="L-BFGS-B", lower=c(0,0))
optim(par=c(0.1,0.1),LLL,method="L-BFGS-B", lower=c(0,0), hessian=TRUE)

ご協力いただきありがとうございます！

4

1 に答える 1

0

optim関数を最小化します。引数 を指定して尤度関数をプロットできますrate。プロットを取得するには、少し手を加える必要があります。次のようにします。

z2 <- function(rate) {
    par <- numeric(2)
    par[1] <- .68
    par[2] <- rate
    y <- LLL(par)
    y
}

z1 <- Vectorize(z2,vectorize.args="rate")

curve(z1, from=.001,to=1)

の最小値に対して関数が最小であることがわかりますrate。fromに変更しても同じです.1。見積もりが正しいかどうか判断できません。

于 2017-01-09T19:53:53.990 に答える