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

[Ben Bolker]

  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