[R] linear constraint optim with bounds/reparametrization

[Ingmar Visser]

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