0

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!

Y. Ma
  • 31
  • 1
  • 1
  • 3
  • 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 Jan 09 '17 at 18:33
  • @Thales Thank you. – Y. Ma Jan 09 '17 at 18:44
  • 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 Jan 09 '17 at 19:05
  • @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 Jan 09 '17 at 19:18

1 Answers1

0

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.

Bhas
  • 1,844
  • 1
  • 11
  • 9