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: rvarad...@jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h tml ---------------------------------------------------------------------------- -------- -----Original Message----- From: r-devel-boun...@r-project.org [mailto:r-devel-boun...@r-project.org] On Behalf Of epur...@stat.berkeley.edu Sent: Wednesday, November 18, 2009 1:10 PM To: r-de...@stat.math.ethz.ch Cc: r-b...@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 epur...@stat.berkeley.edu (510) 642-6154 (office) (510) 642-7892 (fax) ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
constrOptim2 <- function (theta, f, grad, ui, ci, mu = 1e-04, control = list(), method = if (is.null(grad)) "Nelder-Mead" else "BFGS", outer.iterations = 100, outer.eps = 1e-05, hessian=FALSE, ...) { if (!is.null(control$fnscale) && control$fnscale < 0) mu <- -mu R <- function(theta, theta.old, ...) { ui.theta <- ui %*% theta gi <- ui.theta - ci if (any(gi < 0)) return(NaN) gi.old <- ui %*% theta.old - ci bar <- sum(gi.old * log(gi) - ui.theta) if (!is.finite(bar)) bar <- -Inf f(theta, ...) - mu * bar } dR <- function(theta, theta.old, ...) { ui.theta <- ui %*% theta gi <- drop(ui.theta - ci) gi.old <- drop(ui %*% theta.old - ci) dbar <- colSums(ui * gi.old/gi - ui) grad(theta, ...) - mu * dbar } if (any(ui %*% theta - ci <= 0)) stop("initial value not feasible") obj <- f(theta, ...) r <- R(theta, theta, ...) for (i in 1L:outer.iterations) { obj.old <- obj r.old <- r theta.old <- theta fun <- function(theta, ...) { R(theta, theta.old, ...) } gradient <- function(theta, ...) { dR(theta, theta.old, ...) } a <- optim(theta.old, fun, gradient, control = control, method = method, hessian=hessian, ...) 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, ...) # if (obj > obj.old) # this is a bug if (obj > obj.old * sign(mu)) # this is the correct one break } if (i == outer.iterations) { a$convergence <- 7 a$message <- "Barrier algorithm ran out of iterations and did not converge" } if (mu > 0 && obj > obj.old) { a$convergence <- 11 a$message <- paste("Objective function increased at outer iteration", i) } if (mu < 0 && obj < obj.old) { a$convergence <- 11 a$message <- paste("Objective function decreased at outer iteration", i) } a$outer.iterations <- i a$barrier.value <- a$value a$value <- f(a$par, ...) a$barrier.value <- a$barrier.value - a$value a }
______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel