From: Carlo Fezzi <c.fezzi_at_uea.ac.uk>

Date: Mon, 19 May 2008 17:22:48 +0100

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

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.
*