[R] Trouble with optim

[Damon Wischik]

I am having trouble with optim. It claims to have converged to a minimum,
yet it has in the course of the optimization visited many points which are
closer to optimal. I would be grateful for any explanation of this
behaviour.

I'm trying to estimate the parameters in the model
  X ~ Binomial(1,p) * NegBin(mu,theta).
So I define a log likelihood function, and invoke optim thus:
  o <- optim (c(.7,10.3,8.5), fn=loglikfun, method="L-BFGS-B",
              lower=c(.001,.001,.001), upper=c(1,Inf,Inf))
The initial parameters are an educated guess. I made the loglikfun print
out some information every time it's called. An (edited) output is like
this:

Eval fn at 0.7 10.3 8.5 --- Val = 42.70597
Eval fn at 0.701 10.3 8.5 --- Val = 42.70595
Eval fn at 0.699 10.3 8.5 --- Val = 42.70603
...
Eval fn at 0.7425713 21.12820 0.001 --- Val = 64.99
Eval fn at 0.7425713 21.12920 0.002 --- Val = 60.20449
Eval fn at 0.7425713 21.12920 0.001 --- Val = 64.99

> o$val
[1] 64.99
> o$convergence
[1] 0
> o$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

So optim thinks it has found an optimum (i.e. minimum). But my initial
guess is better than optim's answer; and optim has visited many points
which are better than its final answer.

If I choose a different initial guess, like c(.7,10.3,1), optim reaches
the answer I expect.

What is going on?

Damon Wischik.


I'm running R 1.6.2 on Windows 2000.

Sys.info() tells me:

sysname: Windows
version: (build 2195) Service Pack 2
machine: x86
release: NT 5.0

R.Version() tells me:

platform: i386-pc-mingw32
arch: i386
os: mingw32
system: i386, mingw32
status: 
R.version.string: R version 1.6.2, 2003-01-10


My full code is here:

dnegbinp <- function(x,prob,mu,size,log=FALSE) {
  ll <- ifelse(x==0,
          log(1-prob+prob*dnbinom(0*x,mu=mu,size=size)),
          log(prob)+dnbinom(x,mu=mu,size=size,log=TRUE));
  if (log) ll else exp(ll);
  }

x <- c(1,3,0,19,32,0,0,2,23,23)

loglikfun <- function(p) {
    cat("Eval fn at ")
    cat(p)
    s <- -sum(dnegbinp(x, prob=p[[1]],mu=p[[2]],size=p[[3]], log=TRUE))
    cat(" --- Val = ")
    cat(s)
    cat("\n")
    s }

o <- optim (c(.7,10.3,8.5), fn=loglikfun, method="L-BFGS-B",
       lower=c(.001,.001,.001), upper=c(1,Inf,Inf))
o$val
o$convergence
o$message