Hessian matrix in optim() in R

Question

I'm having some trouble using optim() in R to solve for a likelihood involving an integral and obtain the hessian matrix for the 2 parameters. The algorithm converged but I get an error when I use the hessian=TRUE option in optim(). The error is:

Error in integrate(integrand1, lower = s1[i] - 1, upper = s1[i]) : non-finite function value

Also had a warning message of NAs

Here's my code:

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)

Thanks for your help!

I believe you evaluate dgamma(0, shape=.1, scale=.1), which is Inf, when integrating for i == 2 from s1[2]-1 (0) to s1[2] (1). What happens if you exclude 1 from s1? — Thales
You were told yesterday to use lower=c(.001,.001). Why are you not using that for argument lower? Change lower to c(0.01,0.01) and you will see a Hessian. — Bhas
@Bhas Yes I did use lower=c(0.01, 0.01) per your suggestion (copied the wrong code in this post). Any idea why the rate parameter in the gamma always converges to the lower bound and whether the estimate is valid? Thanks. — Y. Ma

Bhas Bhas · Accepted Answer · 2017-01-09T19:53:53

optim minimizes the function. You could plot the likelihood function given the argument rate. It needs a bit of a fiddling to get a plot. Do it like this:

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)

You will see that the function is minimal for the lowest value for rate. Same if you change from to .1. I cannot judge if the estimate is valid.

Hessian matrix in optim() in R

1 Answers