R help - Oct 2007 - Error optimizing Poisson log-likelihood with L-BFGS-B

Dear R users,
I have following code, estimating Poisson log-likelihood a number of times:


poisson.loglik <- function(mu, y){
    n <- NROW(y)
    logl <- sum(y)*log(mu)-n*mu
    return(-logl)
    }

estimates <- numeric(1e5)

for(i in seq_along(estimates)){
estimates[i] <- optim(par=1, poisson.loglik, method="L-BFGS-B",
lower=0, y=rpois(10,lambda=2))$par
}

I cannot run this because I always get an error message:

Error in optim(par = 1, poisson.loglik, lower = 0, method =
"L-BFGS-B",  :
  non-finite finite-difference value [1]

Could someone enlighten me what that means and why this happens?

THanks in advance,
Sergey

Sergey Goriatchev wrote:
> 
> 
> [snip]
> 
> I cannot run this because I always get an error message:
> 
> Error in optim(par = 1, poisson.loglik, lower = 0, method =
"L-BFGS-B",  :
>   non-finite finite-difference value [1]
> 
> Could someone enlighten me what that means and why this happens?
> 
> THanks in advance,
> Sergey
> 
> 
  It happens because "L-BFGS-B" is willing to try points that are
exactly on
the boundary.

## modify likelihood function to allow debug option
poisson.loglik <- function(mu, y,debug=FALSE){
  n <- NROW(y)
  logl <- sum(y)*log(mu)-n*mu
  v = -logl
  if (debug) cat(mu,v,"\n")
  return(v)
}

estimates <- numeric(1e5)

## try the loop, but save yval every time so we can replicate the problem
for (i in 1:1000) {
  yval = rpois(10,lambda=2)
  estimates[i] <- optim(par=1, poisson.loglik, method="L-BFGS-B",
        lower=0, y=yval)$par
}

## now see what crashed (stopped after rep 33 for me)
optim(par=1, poisson.loglik, method="L-BFGS-B",
        lower=0, y=yval)$par

## yes, it crashed. Try debugging.
optim(par=1, poisson.loglik, method="L-BFGS-B",
        lower=0, y=yval, debug=TRUE)$par

## increase lower bound value ...
optim(par=1, poisson.loglik, method="L-BFGS-B",
        lower=0.002, y=yval, debug=TRUE)$par


for (i in 1:10000) {
  yval = rpois(10,lambda=2)
  estimates[i] <- optim(par=1, poisson.loglik, method="L-BFGS-B",
        lower=0.002, y=yval)$par
}

 ## works (and minimum estimate is not on the boundary)
-- 
View this message in context:
http://www.nabble.com/Error-optimizing-Poisson-log-likelihood-with-L-BFGS-B-tf4726144.html#a13514904
Sent from the R help mailing list archive at Nabble.com.