G'day Thomas,
On Mon, 3 Sep 2007 10:08:08 +0200
<[EMAIL PROTECTED]> wrote:
> What's wrong with my code?
>
> Require(quadprog)
library(quadprog) ? :)
> Dmat<-diag(1,7,7)
> # muss als quadratische Matrix eingegeben werden
> Dmat
> dvec<-matrix(0,7,1) # muss als Spaltenvektor eingegeben werden
> dvec
> mu<-0 # (in Mio. €)
> bvec<-c(1,mu,matrix(0,7,1)) # muss als Spaltenvektor eingegeben werden
> bvec
> mu_r<-c(19.7,33.0,0.0,49.7, 82.5, 39.0,11.8)
> Amat<-matrix(c(matrix(1,1,7),7*mu_r,diag(1,7,7)),9,7,byrow=T)
> # muss als Matrix angegeben werden, wie sie wirklich ist
> Amat
> meq<-2
> loesung<-solve.QP(Dmat,dvec,Amat=t(Amat),bvec=bvec,meq=2)
With the commands above, I get on my system:
> loesung<-solve.QP(Dmat,dvec,Amat=t(Amat),bvec=bvec,meq=2)
Error in solve.QP(Dmat, dvec, Amat = t(Amat), bvec = bvec, meq = 2) :
constraints are inconsistent, no solution!
Which is a bit strange, since Dmat is the identity matrix, so there
should be little room for numerical problems.
OTOH, it is known that the Goldfarb-Idnani algorithm can have problem
in rare occasions if the problem is "ill-scaled"; see the work of
Powell.
Changing the definition of Amat to
> Amat<-matrix(c(matrix(1,1,7),mu_r,diag(1,7,7)),9,7,byrow=T)
leads to a successful call to solve.QP. And since mu=0, I really
wonder why you scaled up the vector mu_r by a factor of 7 when putting
it into Amat.... :)
> loesung<-solve.QP(Dmat,dvec,Amat=t(Amat),bvec=bvec,meq=2)
Now, with the solution I obtained the rest of your commands show:
> loesung$solution %*% mu_r
[,1]
[1,] 6.068172e-15
> sum(loesung$solution)
[1] 1
> for (i in 1:7){
a<-loesung$solution[i]>=0
print(a)
}
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] FALSE
[1] FALSE
[1] TRUE
But:
> for (i in 1:7){
a<-loesung$solution[i]>=- .Machine$double.eps*1000
print(a)
}
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
> for (i in 1:7) print(loesung$solution[i])
[1] 0
[1] 0
[1] 1
[1] 0
[1] -4.669169e-17
[1] -1.436586e-17
[1] 8.881784e-16
This is a consequence of finite precision arithmetic. Just as you
should not compare to numeric numbers directly for equality but rather
that the absolute value of their difference is smaller than an
appropriate chosen threshold, you should not check whether a number is
bigger or equal to zero, but whether it is bigger or equal to an
appropriately chosen negative threshold. More information are given in
FAQ 7.31.
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: [EMAIL PROTECTED]
Singapore 117546 http://www.stat.nus.edu.sg/~statba
______________________________________________
R-help@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.
