Problem in optim(method="L-BFGS-B") (PR#1717)
ripley@stats.ox.ac.uk
ripley@stats.ox.ac.uk
Thu, 8 Aug 2002 10:25:11 +0100 (BST)
I can't reproduce this. I get (on Windows R 1.5.1)
iter 0 value 121.713156
iter 1 value 98.361657
iter 2 value 66.561121
iter 3 value 57.491023
iter 4 value 48.193824
iter 5 value 42.298821
iter 6 value 37.554272
iter 7 value 33.741241
iter 8 value 31.592504
iter 9 value 28.856846
iter 10 value 27.450574
iter 11 value 26.807136
iter 12 value 26.525653
iter 13 value 25.553618
iter 14 value 25.424473
iter 15 value 25.312266
iter 16 value 25.310750
iter 17 value 25.309939
iter 18 value 25.307804
iter 19 value 25.307030
iter 20 value 25.302144
iter 21 value 25.294243
iter 22 value 25.283499
iter 23 value 25.276034
iter 24 value 25.274281
iter 25 value 25.273023
iter 26 value 25.272970
iter 27 value 25.272845
iter 28 value 25.272737
iter 29 value 25.272729
iter 30 value 25.272662
iter 31 value 25.272652
iter 32 value 25.272652
iter 33 value 25.272652
final value 25.272652
converged
However, Linux gives similar problems. The underlying cause is that you
are using finite-difference approximations and not scaling your variables
as described on the help page.
On Fri, 28 Jun 2002 polzehl@wias-berlin.de wrote:
> Full_Name: Jörg Polzehl
> Version: 1.5.1
> OS: Windows 2000
> Submission from: (NULL) (193.175.148.198)
>
>
> When calculating MLE's in a variance component model using constrained
> optimization, i.e. optim(...,method="L-BFGS-B",...) I observed an inproper
> behaviour in cases where
> the likelihood function was evalueted at the constraint. Parameters and value of
> the
> function at the constraint seem to be returned in this case.
> I've tried to reproduce the effect in a very simple example:
>
> ###############################################################################
> fn <- function(par,x,y){
> ind <- 1:(length(x)/2)
> alpha <- par[1]
> beta <- par[2]
> sigmaq1 <- par[3]
> sigmaq2 <- par[4]
> sum((y[ind]-alpha-beta*x[ind])^2/sigmaq1+log(sigmaq1)+(y[-ind]-alpha-beta*x[-ind])^2/sigmaq2+log(sigmaq2))
> }
>
> set.seed(3)
> n <- 5
> x1 <- runif(n)
> x2 <- runif(n)
> y <- c(rnorm(x1,1+x1,1),rnorm(x2,1+x2,5))
> x <- c(x1,x2)
> par <- c(0,0,1,1)
> z <- optim(par,fn,method="L-BFGS-B",lower=c(-Inf,-Inf,1e-6,1e-6),x=x,y=y,control=list(trace=6,REPORT=1))
> ####################################################################################
>
> Note that the effect vanishes in case of other constraints for par[3:4].
>
> here is what happens with trace=1 :
>
> z <- optim(par,fn,method="L-BFGS-B",lower=c(-Inf,-Inf,1e-6,1e-6),x=x,y=y,control=list(trace=1,REPORT=1))
> iter 0 value 121.713156
> iter 1 value 98.361657
> iter 2 value 66.561121
> iter 3 value 57.491023
> iter 4 value 48.193824
> iter 5 value 42.298821
> iter 6 value 37.554272
> iter 7 value 33.741241
> iter 8 value 31.592504
> iter 9 value 28.856846
> iter 10 value 27.450574
> iter 11 value 4118799.069149
> final value 4118799.069149
> converged
> ----------------------------------------------------------------------------------
> and the essential part with trace = 6 :
>
> 4 variables are free at GCP on iteration 11
> LINE SEARCH 1 times; norm of step = 1.19245
> X = -3.34057 6.34639 1.18546 15.2136
> G = -1.08654 -1.14056 1.42527 -0.307525
> Iteration 11
>
> ---------------- CAUCHY entered-------------------
>
> There are 1 breakpoints
>
> Piece 1 f1, f2 at start point -4.6074e+000 3.1209e+001
> Distance to the next break point = 8.3175e-001
> Distance to the stationary point = 1.4763e-001
>
> GCP found in this segment
> Piece 1 f1, f2 at start point -4.6074e+000 3.1209e+001
> Distance to the stationary point = 1.4763e-001
> Cauchy X = -3.18017 6.51477 0.97505 15.259
>
> ---------------- exit CAUCHY----------------------
>
> 4 variables are free at GCP on iteration 12 #
> !!!!!!!!!!!!!!!!!!!!!!
> LINE SEARCH 0 times; norm of step = 2.30637
> X = -3.90499 7.26709 1e-006 16.8712
> G = 2.86794e+006 2.02055e+006 -4.1147e+009 -0.191468
>
> iterations 12
> function evaluations 15
> segments explored during Cauchy searches 12
> BFGS updates skipped 0
> active bounds at final generalized Cauchy point 0
> norm of the final projected gradient 4.1147e+009
> final function value 4.1188e+006
>
> X = -3.90499 7.26709 1e-006 16.8712
> F = 4.1188e+006
> final value 4118799.069149
> converged
>
>
> See especially the line marked with !!!!!!!!!!
>
> Regards,
>
> Jörg Polzehl
>
>
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>
--
Brian D. Ripley, ripley@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 272860 (secr)
Oxford OX1 3TG, UK Fax: +44 1865 272595
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._