Dear list, 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 I understand the general principle of creating box constraints by creating pairs of positive and negative values within the constraint matrix ( http://tolstoy.newcastle.edu.au/R/help/05/10/14251.html), but when I pass this expanded constraint matrix to solve.QP it complains that Amat and dvec are incompatible. How should I expand dvec (and consequently Dmat) to accomodate the larger Amat? Moreover, I am unclear how to apply the meq equality constraint across more than one cell (i.e. rows summing to one) although I have attempted a guess below. Any help warmly received. Selwyn. ################ #examples below ################ library(quadprog) #pairs of population values over time mat = cbind( cal = c(12988,13408,13834,14267,14707,15155) , rest = c(152082,154201,156352,158536,160754,163006) ) (X = kronecker(diag(1, ncol(mat)), mat[-nrow(mat),])) (y = c(mat[-1,])) XX = (t(X) %*% X) Xy = t(X) %*% y # a working example of simply constrained solution # i.e. constrain 1st and 4th values to be <1, 2nd and 3rd >0 (a = diag(c(-1,1,1,-1))) (b = c(1,0,0,1)) # this is correct according to these constraints, but I need 0 >= a <= 1 for all a and each row to sum to 1 solve.QP(Dmat = XX,dvec=Xy,Amat = a,bvec=b) #my guess at the constraint matrix and b_vec... (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 kronecker(diag(rep(1,4)),c(1,-1)) #create positive and negative constraint pairs ) ) (b2 = c(1,1,rep(0:1,times=4))) solve.QP(Dmat = XX,dvec=Xy,Amat = a2,bvec=b2,meq=2) #complains that Amat and dvec are incompatible [[alternative HTML version deleted]]
[reposting as a plain text - apologies for the double posting] Dear list, 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 I understand the general principle of creating box constraints by creating pairs of positive and negative values within the constraint matrix (http://tolstoy.newcastle.edu.au/R/help/05/10/14251.html), but when I pass this expanded constraint matrix to solve.QP it complains that Amat and dvec are incompatible. How should I expand dvec (and consequently Dmat) to accomodate the larger Amat? Moreover, I am unclear how to apply the meq equality constraint across more than one cell (i.e. rows summing to one) although I have attempted a guess below. Any help warmly received. Selwyn. ################ #examples below ################ library(quadprog) #pairs of population values over time mat = cbind( cal = c(12988,13408,13834,14267, 14707,15155) , rest = c(152082,154201,156352,158536,160754,163006) ) (X = kronecker(diag(1, ncol(mat)), mat[-nrow(mat),])) (y = c(mat[-1,])) XX = (t(X) %*% X) Xy = t(X) %*% y # a working example of simply constrained solution # i.e. constrain 1st and 4th values to be <1, 2nd and 3rd >0 (a = diag(c(-1,1,1,-1))) (b = c(1,0,0,1)) # this is correct according to these constraints, but I need 0 >= a <1 for all a and each row to sum to 1 solve.QP(Dmat = XX,dvec=Xy,Amat = a,bvec=b) #my guess at the constraint matrix and b_vec... (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 kronecker(diag(rep(1,4)),c(1,-1)) #create positive and negative constraint pairs ) ) (b2 = c(1,1,rep(0:1,times=4))) solve.QP(Dmat = XX,dvec=Xy,Amat = a2,bvec=b2,meq=2) #complains that Amat and dvec are incompatible
G'day Selwyn, On Mon, 16 Feb 2009 10:11:20 +0000 Selwyn McCracken <selwyn.mccracken at 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 <= 1Note 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 incompatibleThe 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: statba at nus.edu.sg Singapore 117546 http://www.stat.nus.edu.sg/~statba