[R] MAXIMIZATION WITH CONSTRAINTS

[Thomas Lumley]

On Fri, 8 Dec 2006, miri.dago at virgilio.it wrote:

> Dear R users, Iÿÿm a graduate students and in my master thesis I must
> obtain the values of the parameters x_i which maximize this Multinomial 
> logÿÿlikelihood function log(n!)-sum_{i=1]^4 log(n_i!)+sum_ {i=1}^4 n_i 
> log(x_i)

>under the following constraints:
>a) sum_i x_i=1, x_i>=0, 
>b) x_1<=x_2+x_3+x_4
>c)x_2<=x_3+x_4
>I have been using the 
>ÿÿConstrOptimÿÿ R-function with the instructions I report below, and I 
>have tried to implement them with different values of ÿÿnÿÿ. BUT I have 
>encountered 2 problems:

The problem is that the first constraint is not an inequality but an 
equality. Writing it as two inequalities results in the feasible region 
for the optimization being a very narrow slice of four-dimensional space, 
which makes the optimization difficult.

There are at least two ways to fix the problem.  The first is to note that 
the loglikelihood is monotone in each x, so that sum_i x_i <=1 is 
sufficient when maximizing.  The second is to reparametrize in terms of 
three parameters.

Minimization is more challenging, because the loglikelihood does not have 
a minimum. It is negative infinite when any x_i is zero and the 
corresponding n is non-zero.

 	-thomas


My R instructions

n=c(10,20,3,5)
n1=n[1]
n2=n
[2]
n3=n[3]
n4=n[4]

logfr=function(x) { ##function to maximize 
x1= x
[1] 
x2= x[2] 
x3= x[3]
x4= x[4]
log(factorial(sum(n)))-sum(log
(factorial(n)))+sum(n*log(x)) 
} 
grr.log <- function(x) { ## Gradient 
of 'log fr'
x1=x[1] 
x2=x[2] 
x3=x[3]
x4=x[4]
return(n/x)
} 
par.start= 
c(.19999999,.15,.4,.25)
constr.coeff = rbind(diag(1,4,4),c(-1,1,1,1),c
(0,-1,1,1),c(-1,-1,-1,-1), c(1,1,1,1))
constr.tn= c(0,0,0,0,0,0,-1,.
9999999)
min= constrOptim(par.start, logfr, grr.log, ui=constr.coeff, 
ci=constr.tn)
max=constrOptim(par.start, logfr, grr.log, ui=constr.
coeff, ci=constr.tn, control=list(fnscale=-1))

______________________________________________
R-help at stat.math.ethz.ch 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.

Thomas Lumley			Assoc. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle