[R] constrOptim converging not to the optimal values

[Carlo Fezzi]

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