# [R] Trouble with optim

Damon Wischik djw1005 at cam.ac.uk
Sat Feb 1 02:56:02 CET 2003

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

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