[Rd] behavior of L-BFGS-B with trivial function triggers bug in stats4::mle
Ben Bolker
bolker at zoo.ufl.edu
Mon Jul 30 14:05:30 CEST 2007
[sent this last night, may have bounced, resending]
With R 2.5.1 ...
"L-BFGS-B" behaves differently from all of the
other optim() methods, which 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 that was originally fitted with "L-BFGS-B".
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) ## BUG
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