I sent this in first on 30 July. Now that UseR! is over I'm trying again
(slightly extended version from last time).
With R 2.5.1 or R 2.6.0 (2007-08-04 r42421)
"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. This is not
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".
One possible fix is to check for length(start)==0 and
return a dummy optim() result in that case (see patch included below).
The patch below fixes the basic problem, although I haven't
tested extensively to see what its other implications are.
Or one could change L-BFGS-B to behave the same as the other methods.
If I don't hear anything in a few days would it be appropriate
to submit this as a bug report?
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
*** mle.R 2007-07-27 11:50:38.000000000 -0400
--- src/library/stats4/R/mle.R 2007-08-13 17:47:11.000000000 -0400
***************
*** 56,62 ****
l[n] <- fixed
do.call("minuslogl", l)
}
! oout <- optim(start, f, method=method, hessian=TRUE, ...)
coef <- oout$par
vcov <- if(length(coef)) solve(oout$hessian) else matrix(numeric(0),0,0)
min <- oout$value
--- 56,66 ----
l[n] <- fixed
do.call("minuslogl", l)
}
! if (length(start)==0) {
! oout <- list(par=numeric(0),value=f(start))
! } else {
! oout <- optim(start, f, method=method, hessian=TRUE, ...)
! }
coef <- oout$par
vcov <- if(length(coef)) solve(oout$hessian) else matrix(numeric(0),0,0)
min <- oout$value
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
