[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
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
More information about the R-help
mailing list