One more bug in `constrOptim' that I have reported before and have also suggested a fix, but this has also not been corrected:
if (obj > obj.old) break This is obviously wrong when the user is trying to maximize by specifying control$fnscale = -1. The correct statement is the following: if (obj > obj.old * sign(mu)) break Ravi. -----Original Message----- From: Ravi Varadhan [mailto:rvarad...@jhmi.edu] Sent: Thursday, June 17, 2010 10:55 AM To: 'Ravi Varadhan'; 'Duncan Murdoch'; 'John Nolan' Cc: 'r-devel@r-project.org' Subject: RE: [Rd] constrOptim( ): conflict between help page and code In `constrOptim', there is also a mistake in the reporting of function and gradient counts. These counts, as reported, correspond to that of the vary last "inner" iteration. However, they should be cumulatively summed up over each outer iteration. I have fixed this and the convergence criterion issue that I mentioned before. Currently, constrOptim can only use the Nelder-Mead when the analytic gradient is not specified. I have added a numerical gradient option so that it can use BFGS or other optimization functions even when the analytic gradient is not specified. It would be desirable if Tom Lumley can commit these changes to constrOptim. In the meanwhile, I can send the function with these upgrades to anyone who is interested. Best, Ravi. -----Original Message----- From: r-devel-boun...@r-project.org [mailto:r-devel-boun...@r-project.org] On Behalf Of Ravi Varadhan Sent: Thursday, June 17, 2010 9:53 AM To: 'Duncan Murdoch'; 'John Nolan' Cc: r-devel@r-project.org Subject: Re: [Rd] constrOptim( ): conflict between help page and code Nolan, You are correct that there is inconsistency. The feasible region should be ui %*% theta - ci > 0, so that log(ui %*% theta - ci) is defined. There is a more serious problem in termination criterion for iterations: if (is.finite(r) && is.finite(r.old) && abs(r - r.old)/(outer.eps + abs(r - r.old)) < outer.eps) break Ideally convergence is achieved when |r - r.old| is small. But according to the above code, the logical condition inside the if(.) will be TRUE only when abs(r - r.old) < (outer.eps)^2 (for small outer.eps). For example, let |r - r.old| = outer.eps. The above condition becomes: if (0.5 < outer.eps) break, which will obviously never happen for values of outer.eps < 1/2. Now, if outer.eps = 1.e-05 (default) and we let |r - r.old| <= 1.e-10, then the above condition will be satisfied. In short, the termination criterion used is too stringent. Better termination criteria are: if (is.finite(r) && is.finite(r.old) && abs(r - r.old) < outer.eps) break or if (is.finite(r) && is.finite(r.old) && abs(r - r.old)/(1 + abs(r + r.old)/2) < outer.eps) break Best, Ravi. -----Original Message----- From: r-devel-boun...@r-project.org [mailto:r-devel-boun...@r-project.org] On Behalf Of Duncan Murdoch Sent: Thursday, June 17, 2010 4:38 AM To: John Nolan Cc: r-devel@r-project.org Subject: Re: [Rd] constrOptim( ): conflict between help page and code John Nolan wrote: > There is a contradiction between what the help page says and what constrOptim > actually > does with the constraints. The issue is what happens on the boundary. > I don't know if this was a recent change, but the current help says: "The starting value must be in the interior of the feasible region, but the minimum may be on the boundary." The boundary is not in the interior. Duncan Murdoch > The help page says > The feasible region is defined by ‘ui %*% theta - ci >= 0’, > but the R code for constrOptim reads > if (any(ui %*% theta - ci <= 0)) > stop("initial value not feasible") > The following example shows that when the initial point is on the boundary of > the > feasibility region, you get the above error message and execution stops. > > >> fn <- function(x) { return(sum(x)) } >> >> ui <- diag(rep(1,2)) >> ci <- matrix(0,nrow=2,ncol=1) >> constrOptim( c(0,0), fn, NULL, ui, ci ) >> > Error in optim(theta.old, fun, gradient, control = control, method = method, > : > function cannot be evaluated at initial parameters > >> version >> > platform i386-pc-mingw32 > arch i386 > os mingw32 > system i386, mingw32 > status > major 2 > minor 10.0 > year 2009 > month 10 > day 26 > svn rev 50208 > language R > version.string R version 2.10.0 (2009-10-26) > > In contrast, at a different place in constrOptim - the internal function R - > the boundary of the feasibility region is allowed: if (any(gi < 0)) > return(NaN), > and it seems to explicitly allow boundaries at another place: > allowing gi==0 and interpreting log(gi) as -Inf. > > > John > > ........................................................................... > > John P. Nolan > Math/Stat Department > 227 Gray Hall > American University > 4400 Massachusetts Avenue, NW > Washington, DC 20016-8050 > > jpno...@american.edu > 202.885.3140 voice > 202.885.3155 fax > http://academic2.american.edu/~jpnolan > ........................................................................... > ______________________________________________ > R-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel > ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
