From epurdom at stat.berkeley.edu  Wed Nov 18 19:10:15 2009
From: epurdom at stat.berkeley.edu (epurdom at stat.berkeley.edu)
Date: Wed, 18 Nov 2009 19:10:15 +0100 (CET)
Subject: [Rd] bug in '...' of constrOptim (PR#14071)
Message-ID: <20091118181016.4F3362831FFC@mail.pubhealth.ku.dk>

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)



[Rd] bug in '...' of constrOptim (PR#14071)
epurdom at stat.berkeley.edu