Hessian matrix with optim or numderiv package?

Question

I do maximum likelihood estimation for a loglikelihood function of a poisson distribution. After the estimation i want to compute the standard errors of the coeffients. For that i need the hessian matrix. Now i dont know which function i should use to get the hessian matrix , optim() or hessian() from the numderiv package.

Both function give me different solution. And if i try to compute Standard errors from the hessian that i get from optim i get one NaN entry in the result. Whats the difference between these two functions for the compution of the hessian matrix?

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 )))

Hessian matrix with optim or numderiv package?

0 Answers0