Hello All, I would like to optimize a (log-)likelihood function subject to a number of linear constraints between parameters. These constraints are equality constraints of the form A%*%theta=c, ie (1,1) %*% 0.8,0.2)^t = 1 meaning that these parameters should sum to one. Moreover, there are bounds on the individual parameters, in most cases that I am considering parameters are bound between zero and one because they are probabilities. My problems/questions are the following:
1) constrOptim does not work in this case because it only fits inequality constraints, ie A%*%theta > = c --- As a result when providing starting values constrOptim exits with an error saying that the initial point is not feasible, which it isn't because it is not in the interior of the constraint space. Is there an alternative to constrOptim that can handle such strict (equality) linear constraints? 2) Another option of course would be to reparametrize the problem as follows. I will illustrate with an example: I have parameters: > p [,1] [1,] 0.8 [2,] 0.2 [3,] 0.2 [4,] 0.8 [5,] 0.6 [6,] 0.1 [7,] 0.3 [8,] 0.1 [9,] 0.3 [10,] 0.6 [11,] 0.5 [12,] 0.5 and the following constraints (all these constraint amount to certain probabilities summing to one, ie the first two parameters sum to one, the second pair of parameters sum to one etc): > A [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [1,] 1 1 0 0 0 0 0 0 0 0 0 0 [2,] 0 0 1 1 0 0 0 0 0 0 0 0 [3,] 0 0 0 0 1 1 1 0 0 0 0 0 [4,] 0 0 0 0 0 0 0 1 1 1 0 0 [5,] 0 0 0 0 0 0 0 0 0 0 1 1 and the bounds on the constraints are: > ci [,1] [1,] 1 [2,] 1 [3,] 1 [4,] 1 [5,] 1 The bounds on the parameters are all zero and one. In order to reparametrize I could use z=ginv(b)%*%(p-ginv(A)%*%cc) with b=Null(t(A)) and optimize z instead of p using p=ginv(a)%*%ci + b%*%z My question is however what the bounds on z are? In the above example z is: > z [,1] [1,] -0.23333333 [2,] -0.03333333 [3,] -0.57974349 [4,] -0.37974349 [5,] -0.07974349 [6,] -0.18856181 [7,] -0.18856181 which conforms to the constraints, so these values can be used as an intial estimate. If I knew the bounds on z I could use optim with method="L-BFGS". I hope this problem is sufficiently clear to elicit response, if not please let me know. ingmar ps: to reproduce above example use the following: p=matrix(c(0.8, 0.2, 0.20, 0.8, 0.6, 0.1, 0.3, 0.1, 0.3, 0.6, 0.5, 0.5),nc=1) A = matrix(c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1),nc=12,byrow=T) ci=matrix(rep(1,5),nc=1) b=Null(t(A)) z=ginv(b)%*%(p-ginv(A)%*%cc) pp=ginv(a)%*%ci + b%*%z ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html