Thanks Douglas, I will study your code and Eigen carefully. Best wishes!
Honglang Wang Office C402 Wells Hall Department of Statistics and Probability Michigan State University 1579 I Spartan Village, East Lansing, MI 48823 wangh...@msu.edu On Thu, Dec 6, 2012 at 1:33 PM, Douglas Bates <ba...@stat.wisc.edu> wrote: > I don't know as much about Armadillo as I do about Eigen so I will cheat > and write using RcppEigen instead of RcppArmadillo. Page 6 of the Eigen > tutorial at http://eigen.tuxfamily.org/dox/ discusses decompositions and > solving linear systems of equations. One of the simplest ways of solving a > weighted least squares system, if you only want the solution itself, is > with a QR decomposition of which there are 3 in Eigen. The simplest is > HouseholderQR (named after A.S. Householder who first described an > elementary reflector). For a least squares solution (assuming the model > matrix has full column rank) you simply use the solve method. For a > weighted least squares solution you premultiply both the model matrix and > the response by the square roots of the weights. > > It can be collapsed to a one-liner, as in the enclosed. A test that it > works for unweighted least squares (i.e. using unit weights) is > > > sourceCpp("/tmp/wtls.cpp") > > coef(lm(optden ~ carb, Formaldehyde))(Intercept) carb > 0.005085714 0.876285714 > > with(Formaldehyde, wtls(model.matrix(~carb), optden, rep.int(1, > length(optden)))) > $betahat > [1] 0.005085714 0.876285714 > > The use of the other decompositions (Cholesky, ColPivHouseholderQR, > Symmetric Eigendecomposition, SVD) are described in the RcppEigen vignette. > > By the way, the sourceCpp, evalCpp and cppFunction capabilities developed > by the Rcpp authors is very close to magic. They are to be congratulated > on a remarkable accomplishment. > > > On Wed, Dec 5, 2012 at 1:16 PM, Tim Triche, Jr. <tim.tri...@gmail.com>wrote: > >> this part will always make your code crawl along: >> >> >> On Wed, Dec 5, 2012 at 11:09 AM, Honglang Wang < >> wanghonglang2...@gmail.com> wrote: >> >>> arma::vec betahat = arma::inv(Inv)*arma::trans(D)*W*y; >>> >> >> first time I wrote a GLM engine, I wrote it the way statistics books >> illustrate it (i.e. actually invert things to do IWLS) and it crawled. I >> took it to a physicist friend who went through alternate stages of disgust >> and laughter then told me never to invert anything I didn't have to. >> >> You should take Doug Bates' advice, it could save you a great deal of >> time. Never bit fiddle when you can use a better algorithm. >> >> >> -- >> *A model is a lie that helps you see the truth.* >> * >> * >> Howard Skipper<http://cancerres.aacrjournals.org/content/31/9/1173.full.pdf> >> >> >> _______________________________________________ >> Rcpp-devel mailing list >> Rcpp-devel@lists.r-forge.r-project.org >> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel >> > >
_______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
