1 #include <RcppArmadillo.h>
   2 
   3 extern "C" SEXP fastLm(SEXP ys, SEXP Xs) {
   4     Rcpp::NumericVector yr(ys);         // creates Rcpp vector from SEXP
   5     Rcpp::NumericMatrix Xr(Xs);         // creates Rcpp matrix from SEXP
   6     int n = Xr.nrow(), k = Xr.ncol();
   7 
   8     arma::mat X(Xr.begin(), n, k, false);       // reuses memory and avoids extra copy
   9     arma::colvec y(yr.begin(), yr.size(), false);
  10 
  11     arma::colvec coef = arma::solve(X, y);      // fit model y ~ X
  12     arma::colvec res = y - X*coef;      // residuals
  13 
  14     double s2 = std::inner_product(res.begin(), res.end(), res.begin(), double())/(n - k);
  15                         // std.errors of coefficients
  16     arma::colvec stderr = arma::sqrt(s2 * arma::diagvec( arma::inv(arma::trans(X)*X) ));    
  17 
  18     return Rcpp::List::create(Rcpp::Named("coefficients") = coef,
  19                   Rcpp::Named("stderr")       = stderr,
  20                   Rcpp::Named("df")           = n - k
  21                   );
  22 }
  23