2

我对泊松分布的对数似然函数进行最大似然估计。在估计之后，我想计算系数的标准误差。为此，我需要粗麻布矩阵。现在我不知道应该使用哪个函数从 numderiv 包中获取 hessian 矩阵、optim() 或 hessian()。

这两个功能都给了我不同的解决方案。如果我尝试从我从 optim 得到的 hessian 计算标准错误，我会在结果中得到一个 NaN 条目。这两个函数用于计算 hessian 矩阵有什么区别？

logLikePois <- function(parameter,y, z) {

  betaKoef <- parameter

  lambda <- exp(betaKoef %*% t(z))

  logLikeliHood <- -(sum(-lambda+y*log(lambda)-log(factorial(y))))
  return(logLikeliHood)

}
grad <- function (y,z,parameter){
    betaKoef <- parameter

    # Lamba der Poissonregression
    lambda <- exp(betaKoef%*%t(z))
   gradient <-  -((y-lambda)%*%(z))
   return(gradient)
  }

data(discoveries)
disc <- data.frame(count=as.numeric(discoveries),
                   year=seq(0,(length(discoveries)-1),1))

yearSqr <- disc$year^2

formula <- count ~ year + yearSqr

form <- formula(formula)

model <- model.frame(formula, data = disc)

z <- model.matrix(formula, data = disc)

y <- model.response(model)

parFullModell <- rep(0,ncol(z))

optimierung <- optim(par = parFullModell,gr=grad, fn = logLikePois,
                     z = z, y = y, method = "BFGS" ,hessian = TRUE)

optimHessian <- optimierung$hessian

 numderivHessian <- hessian(func = logLikePois, x = optimierung$par, y=y,z=z)

sqrt(diag(solve(optimHessian)))
sqrt(diag(solve(numderivHessian )))
4

0 回答 0