Berwin, Many thanks for your reply - your solution is exactly what I was looking for and the additional advice is also much appreciated.
Have a good day :-) Selwyn On Mon, Feb 16, 2009 at 3:13 PM, Berwin A Turlach <ber...@maths.uwa.edu.au> wrote: > G'day Selwyn, > > On Mon, 16 Feb 2009 10:11:20 +0000 > Selwyn McCracken <selwyn.mccrac...@gmail.com> wrote: > >> I am trying to follow an example that estimates a 2x2 markov >> transition matrix across several periods from aggregate data using >> restricted least squares. >> >> I seem to be making headway using solve.QP(quadprog) as the >> unrestricted solution matches the example I am following, and I can >> specify simple equality and inequality constraints. However, I cannot >> correctly specify a constraint matrix (Amat) with box constraints per >> cell and equality constraints that span multiple cells. Namely the >> solution matrix I am aiming for needs to respect the following >> conditions: >> - each row must sum to 1 >> - each cell must >= 0 >> - each cell must <= 1 > > Note that this set of constraints contains some redundancy. If each > row has to sum to one and the summands are all negative then, > necessarily, each summand is at most one. So the last set of > constraints is not necessary. > > [...] >> (b2 = c(1,1,rep(0:1,times=4))) > > This should be > (b2 <- c(1,1, rep(0:(-1), times=4))) > > Since x <= 1 iff -x >= -1 > >> solve.QP(Dmat = XX,dvec=Xy,Amat = a2,bvec=b2,meq=2) #complains that >> Amat and dvec are incompatible > > The constraints of the QP are A^T b >= b0 and b0 is passed to b2 while > A is passed to Amat. Your matrix a2 contains A^T, so you have to pass > t(a2) to solve.QP: > > R> solve.QP(Dmat = XX,dvec=Xy,Amat = t(a2),bvec=b2,meq=2) > Error in solve.QP(Dmat = XX, dvec = Xy, Amat = t(a2), bvec = b2, meq = 2) : > constraints are inconsistent, no solution! > > While this might not be very helpful, the problem now is that the > entries in your matrix XX (and Xy) are quite large and this can lead to > numerical problems in the FORTRAN code that quadprog relies on. But > this is easily fixed. > > R> XX <- XX*1e-9 > R> Xy <- Xy*1e-9 > R> solve.QP(Dmat = XX,dvec=Xy,Amat = t(a2),bvec=b2,meq=2) > $solution > [1] 0.9367934 0.0000000 0.0632066 1.0000000 > > $value > [1] -63.38694 > > $unconstrainted.solution > [1] 1.004181694 0.002401266 0.012673485 1.012848987 > > $iterations > [1] 4 0 > > $iact > [1] 10 1 2 > > But, as pointed out earlier, your constraints contain some > redundancies, so it would be shorter to code them as: > > R> a2 = rbind( > c(1,0,1,0), #specify the two cells in the 'row' I want to sum to one(???) > c(0,1,0,1), # likewise > diag(rep(1,4)) > ) > R> b2 <- c(1,1, rep(0, times=4)) > R> solve.QP(Dmat = XX,dvec=Xy,Amat = t(a2),bvec=b2,meq=2) > $solution > [1] 0.9367934 0.0000000 0.0632066 1.0000000 > > $value > [1] -63.38694 > > $unconstrainted.solution > [1] 1.004181694 0.002401266 0.012673485 1.012848987 > > $iterations > [1] 4 0 > > $iact > [1] 1 2 4 > > > HTH. > > Cheers, > > Berwin > > =========================== Full address ============================= > Berwin A Turlach Tel.: +65 6515 4416 (secr) > Dept of Statistics and Applied Probability +65 6515 6650 (self) > Faculty of Science FAX : +65 6872 3919 > National University of Singapore > 6 Science Drive 2, Blk S16, Level 7 e-mail: sta...@nus.edu.sg > Singapore 117546 http://www.stat.nus.edu.sg/~statba > ______________________________________________ R-help@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.
