R中的optim()中的Hessian矩阵

Question

我在使用 R 中的 optim（） 求解涉及积分的可能性并获得 2 个参数的黑森矩阵时遇到了一些问题。算法收敛了，但是当我在optim（）中使用hessian=TRUE选项时出现错误。错误是：

积分错误（积分1， 下限 = s1[i] - 1， 上部 = 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)

感谢您的帮助！

Answer 1

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 也是如此。我无法判断估计是否有效。