[Rd] behavior of L-BFGS-B with trivial function triggers bug in stats4::mle

Ben Bolker bolker at zoo.ufl.edu
Mon Jul 30 00:04:23 CEST 2007


   With the exception of "L-BFGS-B", all of the
other optim() methods return the value of the function
when they are given a trivial function (i.e., one with no
variable arguments) to optimize.  I don't think this
is a "bug" in L-BFGS-B (more like a response to
an undefined condition), but it leads to a bug in stats4::mle --
 a spurious error saying that a better fit
has been found during profiling if one tries to profile
a 1-parameter model.

  I haven't dug quite all the way to the bottom of this
yet, but the attached code will clearly show the problem.

  In the version of mle that I've built (which has
gotten some ugly bells and whistles added) I added
a check which would be more or less equivalent
to check if length(start)==0 and then setting

    oout <- list(par=start, value=f(start),
                 hessian = matrix(numeric(0),0,0)
 
 or something along those lines.

  Or one could change L-BFGS-B to behave the same
as the other methods.

 cheers
   Ben Bolker

---------------------
library(stats4)

## using example from ?mle
x <- 0:10
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
ll <- function(ymax=15, xhalf=6)
  -sum(stats::dpois(y, lambda=ymax/(1+x/xhalf), log=TRUE))

## fix one parameter to get 1D profile
fit2 <- mle(ll, fixed=list(xhalf=6))
profile(fit2)

## same again with method="L-BFGS-B"
fit3 <- mle(ll, fixed=list(xhalf=6),method="L-BFGS-B")
profile(fit3)

ll0 <- function(zzz) {
  ymax <- 15
  xhalf <- 6
  -sum(stats::dpois(y, lambda=ymax/(1+x/xhalf), log=TRUE))
}

## try mle() with all-fixed parameters with various methods ...
methods = eval(formals(optim)$method)
sapply(methods,
       function(m) {
         -logLik(mle(ll, start=list(ymax=15,xhalf=6),
                     fixed=list(ymax=15,xhalf=6),method=m))
       })
##   Nelder-Mead          BFGS            CG      L-BFGS-B          SANN
##  3.389441e+01  3.389441e+01  3.389441e+01 5.048277e-270  3.389441e+01



More information about the R-devel mailing list