[R] constrOptim converging not to the optimal values

From: Carlo Fezzi <c.fezzi_at_uea.ac.uk>
Date: Mon, 19 May 2008 17:22:48 +0100


Dear helpers,

I am using constrOptim to minimize a function subject to inequality constraint. It works well when the number of parameters to optimize is low (e.g. 4) but when they are more (e.g. 10) it does not produce the expected results.

The function is minus the determinant of a binomial logit model with one explanatory variable.

Calling “xeta” the vector of explanatory variables the function to be optimized is:

negdet<-function(xeta)
{

	pi <- exp(a + b* xeta)/(1+exp(a + b* xeta))
	-( sum(pi*(1-pi))*sum(pi*(1-pi)*xeta^2) - sum(pi*(1-pi)*xeta)^2 ) 
}

With “a” intercept and “b” parameter of xeta in a standard binomial logit, i.e. p(i) = exp(a + b * xeta(i)) / (1 + exp(a + b*xeta(i) ) . I assume each element of xeta bounded between 0 and 10 and, in order to obtain an unique solution, I give an order to the elements of xeta (i.e. xeta(1) < xeta(2) < ... xeta(N)).

The solution should be to fix half of the elements of xeta to the values that make the probability pi = 0.18 and half of them to the values that make pi = 0.82. This maximizes the determinant. For xeta constituted by only few elements (e.g. 2 or 4) constrOptim works well, but when the number of elements increases (e.g. 10) the procedure does not converge to the optimal values.

Here is the code, with sample size = 4 and a = 5 and b= -2 the results are correct:

$par

[1] 3.272623 3.270943 1.727875 1.727867

$convergence

[1] 0

But with sample size = 10, the optimization have problems:

$par

 [1] 9.999999 9.999998 7.396835 3.281640 3.278376 3.277607 3.276871 1.737043  [9] 1.736261 1.732932

$convergence

[1] 0

I would be extremely grateful if anybody could point out a reason for this and in particular if you could suggest any other procedure I could use to solve this optimization problem.

Thanks a lot for your kindness,

Carlo

The code is below.

#### SAMPLE SIZE N <- 4

####

a<- 5
b<- -2

negdet<-function(xeta)
{

	pi <- exp(a + b* xeta)/(1+exp(a + b* xeta))
	-( sum(pi*(1-pi))*sum(pi*(1-pi)*xeta^2) - sum(pi*(1-pi)*xeta)^2 ) 
}

### BOUNDS low <- 0
high <- 10

###

### MATRIX REPRESENTAION OF THE CONSTRAINTS order<-c(rep(c(1,-1,rep(0,times=N-1)),N-2),1,-1) u1 <- matrix(order,ncol=N, nrow=N-1, byrow=TRUE)

bounds <- c(rep(c(1,-1,rep(0,times=N*2)),N-1),1,-1) u2 <- matrix(bounds, ncol=N,nrow=N*2)

u.s <-rbind(u1,u2)
c.s <-c(rep(0,N-1),rep(c(low,-high),N))

####

#### OPTIMIZATION a<-constrOptim(outer.iteration = 500, control = list(maxit=10000), theta=(high-low)/(N+1)* N:1, f=negdet, ui = u.s, ci=c.s, method="Nelder-Mead")

####



Carlo Fezzi
Senior Research Associate
Centre for Social Research on the Global Environment (CSERGE) Department of Environmental Sciences
University of East Anglia
Norwich (UK) NR2 7TJ
Telephone: +44(0)1603 591408
Fax: +44(0)1603 593739

R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Mon 19 May 2008 - 17:28:42 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 19 May 2008 - 17:30:38 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive