epurdom at stat.berkeley.edu
2009-Nov-18 18:10 UTC
[Rd] bug in '...' of constrOptim (PR#14071)
Dear all,
There appears to be a bug in how constrOptim handles ... arguments that
are suppose to be passed to optim, according to the documentation. This
means you can't get the hessian to be returned, for example (so this is
a real problem, and not just a question of mistaken documentation).
Looking at the code, it appears that a call to the user-defined f
includes the ..., when the ... should only be passed to the optim call.
I get around it by putting a dummy 'hessian=TRUE' argument in my f
function, but this is not how the function should work.
I show the relevant problem in the code, give an example, and also give
my sessionInfo() call below.
Thanks,
Elizabeth
##Copy of the relevant segment of the code of constrOptim and where I
think the problem might be:
for (i in 1L:outer.iterations) {
obj.old <- obj
r.old <- r
theta.old <- theta
fun <- function(theta, ...) { ##this one's okay
R(theta, theta.old, ...)##this one's okay
}
gradient <- function(theta, ...) {##this one's okay
dR(theta, theta.old, ...)##this one's okay
}
a <- optim(theta.old, fun, gradient, control = control,
method = method, ...)##this one's okay
r <- a$value
if (is.finite(r) && is.finite(r.old) && abs(r -
r.old)/(outer.eps +
abs(r - r.old)) < outer.eps)
break
theta <- a$par
obj <- f(theta, ...)##this one's NOT okay
if (obj > obj.old)
break
}
###Here is an example modified from the examples of the help page of
constrOptim:
> fr <- function(x) { ## Rosenbrock Banana function
+ x1 <- x[1]
+ x2 <- x[2]
+ 100 * (x2 - x1 * x1)^2 + (1 - x1)^2
+ }
> grr <- function(x) { ## Gradient of 'fr'
+ x1 <- x[1]
+ x2 <- x[2]
+ c(-400 * x1 * (x2 - x1 * x1) - 2 * (1 - x1),
+ 200 * (x2 - x1 * x1))
+ }
> constrOptim(c(-1.2,0.9), fr, grr, ui=rbind(c(-1,0),c(0,-1)), ci=c(-1,-1))
$par
[1] 0.9999761 0.9999522
$value
[1] 5.708626e-10
$counts
function gradient
6 1
$convergence
[1] 0
$message
NULL
$outer.iterations
[1] 14
$barrier.value
[1] -0.0001999198
> constrOptim(c(-1.2,0.9), fr, grr, ui=rbind(c(-1,0),c(0,-1)),
ci=c(-1,-1),hessian=TRUE)
Error in f(theta, ...) : unused argument(s) (hessian = TRUE)
#add in a dummy argument to my f function
> fr <- function(x,hessian=TRUE) { ## Rosenbrock Banana function
+ x1 <- x[1]
+ x2 <- x[2]
+ 100 * (x2 - x1 * x1)^2 + (1 - x1)^2
+ }
> constrOptim(c(-1.2,0.9), fr, grr, ui=rbind(c(-1,0),c(0,-1)),
ci=c(-1,-1),hessian=TRUE)
$par
[1] 0.9999761 0.9999522
$value
[1] 5.708626e-10
$counts
function gradient
6 1
$convergence
[1] 0
$message
NULL
$hessian
[,1] [,2]
[1,] 801.9599 -399.9905
[2,] -399.9905 199.9952
$outer.iterations
[1] 14
$barrier.value
[1] -0.0001999198
###My Session info
> sessionInfo()
R version 2.10.0 (2009-10-26)
x86_64-apple-darwin9.8.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] GEOquery_2.10.0 RCurl_1.2-1 bitops_1.0-4.1
Biobase_2.5.8 biomaRt_2.1.0 projectManager_1.0 lattice_0.17-26
RColorBrewer_1.0-2
[9] XML_2.6-0
loaded via a namespace (and not attached):
[1] grid_2.10.0 tools_2.10.0
--
Elizabeth Purdom
Assistant Professor
Department of Statistics
UC, Berkeley
Evans Hall, Rm 433
epurdom at stat.berkeley.edu
(510) 642-6154 (office)
(510) 642-7892 (fax)
Hi Elizabeth,
I have fixed this problem. See the attached function. There is also
another bug that can create a problem when control$fnscale = -1, i.e. when
the objective function is to be maximized. I have commented this in the
code and have also fixed this problem.
There is a word of caution. Hessian is problematic when the converged
parameter estimate is on the boundary of the feasible region, since the
derivatives are not defined there.
Also, note that the default method is Nelder-Mead. If you want to use
"BFGS", you have to specify the gradient function. I have written an
extension of constrOptim that can use numerical gradients, and can also
handle nonlinear constraints.
Hope this helps,
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvaradhan at jhmi.edu
Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-devel-bounces at r-project.org [mailto:r-devel-bounces at r-project.org]
On Behalf Of epurdom at stat.berkeley.edu
Sent: Wednesday, November 18, 2009 1:10 PM
To: r-devel at stat.math.ethz.ch
Cc: R-bugs at r-project.org
Subject: [Rd] bug in '...' of constrOptim (PR#14071)
Dear all,
There appears to be a bug in how constrOptim handles ... arguments that
are suppose to be passed to optim, according to the documentation. This
means you can't get the hessian to be returned, for example (so this is
a real problem, and not just a question of mistaken documentation).
Looking at the code, it appears that a call to the user-defined f
includes the ..., when the ... should only be passed to the optim call.
I get around it by putting a dummy 'hessian=TRUE' argument in my f
function, but this is not how the function should work.
I show the relevant problem in the code, give an example, and also give
my sessionInfo() call below.
Thanks,
Elizabeth
##Copy of the relevant segment of the code of constrOptim and where I
think the problem might be:
for (i in 1L:outer.iterations) {
obj.old <- obj
r.old <- r
theta.old <- theta
fun <- function(theta, ...) { ##this one's okay
R(theta, theta.old, ...)##this one's okay
}
gradient <- function(theta, ...) {##this one's okay
dR(theta, theta.old, ...)##this one's okay
}
a <- optim(theta.old, fun, gradient, control = control,
method = method, ...)##this one's okay
r <- a$value
if (is.finite(r) && is.finite(r.old) && abs(r -
r.old)/(outer.eps +
abs(r - r.old)) < outer.eps)
break
theta <- a$par
obj <- f(theta, ...)##this one's NOT okay
if (obj > obj.old)
break
}
###Here is an example modified from the examples of the help page of
constrOptim:
> fr <- function(x) { ## Rosenbrock Banana function
+ x1 <- x[1]
+ x2 <- x[2]
+ 100 * (x2 - x1 * x1)^2 + (1 - x1)^2
+ }
> grr <- function(x) { ## Gradient of 'fr'
+ x1 <- x[1]
+ x2 <- x[2]
+ c(-400 * x1 * (x2 - x1 * x1) - 2 * (1 - x1),
+ 200 * (x2 - x1 * x1))
+ }
> constrOptim(c(-1.2,0.9), fr, grr, ui=rbind(c(-1,0),c(0,-1)), ci=c(-1,-1))
$par
[1] 0.9999761 0.9999522
$value
[1] 5.708626e-10
$counts
function gradient
6 1
$convergence
[1] 0
$message
NULL
$outer.iterations
[1] 14
$barrier.value
[1] -0.0001999198
> constrOptim(c(-1.2,0.9), fr, grr, ui=rbind(c(-1,0),c(0,-1)),
ci=c(-1,-1),hessian=TRUE)
Error in f(theta, ...) : unused argument(s) (hessian = TRUE)
#add in a dummy argument to my f function
> fr <- function(x,hessian=TRUE) { ## Rosenbrock Banana function
+ x1 <- x[1]
+ x2 <- x[2]
+ 100 * (x2 - x1 * x1)^2 + (1 - x1)^2
+ }
> constrOptim(c(-1.2,0.9), fr, grr, ui=rbind(c(-1,0),c(0,-1)),
ci=c(-1,-1),hessian=TRUE)
$par
[1] 0.9999761 0.9999522
$value
[1] 5.708626e-10
$counts
function gradient
6 1
$convergence
[1] 0
$message
NULL
$hessian
[,1] [,2]
[1,] 801.9599 -399.9905
[2,] -399.9905 199.9952
$outer.iterations
[1] 14
$barrier.value
[1] -0.0001999198
###My Session info
> sessionInfo()
R version 2.10.0 (2009-10-26)
x86_64-apple-darwin9.8.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] GEOquery_2.10.0 RCurl_1.2-1 bitops_1.0-4.1
Biobase_2.5.8 biomaRt_2.1.0 projectManager_1.0 lattice_0.17-26
RColorBrewer_1.0-2
[9] XML_2.6-0
loaded via a namespace (and not attached):
[1] grid_2.10.0 tools_2.10.0
--
Elizabeth Purdom
Assistant Professor
Department of Statistics
UC, Berkeley
Evans Hall, Rm 433
epurdom at stat.berkeley.edu
(510) 642-6154 (office)
(510) 642-7892 (fax)
______________________________________________
R-devel at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: constrOptim2.txt
URL:
<https://stat.ethz.ch/pipermail/r-devel/attachments/20091118/e57287dc/attachment.txt>