[R] Trouble with optim

[ripley@stats.ox.ac.uk]

I thinks it has found a *local* minimum.  It is not a global optimizer.

You appear to have made no attempt to scale the problem as the help page 
suggests you need to.  Please look into this, as from what little 
information you give I think the initial step has overshot.

On Sat, 1 Feb 2003, Damon Wischik wrote:

> 
> 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
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> http://www.stat.math.ethz.ch/mailman/listinfo/r-help
> 

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595