sourceCpp(code = '
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::export]]
List fastLm(NumericVector yr, NumericMatrix Xr) {
int n = Xr.nrow(), k = Xr.ncol();
arma::mat X(Xr.begin(), n, k, false);
arma::colvec y(yr.begin(), yr.size(), false);
arma::colvec coef = arma::solve(X, y);
arma::colvec resid = y - X*coef;
double sig2 = arma::as_scalar(arma::trans(resid)*resid/(n-
k));
arma::colvec stderrest = arma::sqrt(
sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) );
return List::create(Named("coefficients") = coef,
Named("stderr") = stderrest);
}')
X = matrix(rnorm(100*2), ncol = 2)
y = 0.5 + X %*% c(0.3, -0.7) + rnorm(100)
fastLm(y, X)

sourceCpp(code = '
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::export]]
List fastLm(NumericVector yr, NumericMatrix Xr) {
int n = Xr.nrow(), k = Xr.ncol();
arma::mat X(Xr.begin(), n, k, false);
arma::colvec y(yr.begin(), yr.size(), false);
arma::colvec coef = arma::solve(X, y);
arma::colvec resid = y - X*coef;
double sig2 = arma::as_scalar(arma::trans(resid)*resid/(n-
k));
arma::colvec stderrest = arma::sqrt(
sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) );
return List::create(Named("coefficients") = coef,
Named("stderr") = stderrest);
}')
X = matrix(rnorm(100*2), ncol = 2)
y = 0.5 + X %*% c(0.3, -0.7) + rnorm(100)
fastLm(y, X)

sourceCpp(code = '// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
using namespace Rcpp;
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::VectorXd;
// [[Rcpp::export]]
List fastLm2(NumericVector yr, NumericMatrix Xr) {
const Map<MatrixXd> X(as<Map<MatrixXd> >(Xr));
const Map<VectorXd> y(as<Map<VectorXd> >(yr));
int n = Xr.nrow(), k = Xr.ncol();
VectorXd coef = (X.transpose() * X).llt().solve(X.transpose() * y.col(0));
VectorXd resid = y - X*coef;
double sig2 = resid.squaredNorm() / (n - k);
VectorXd stderrest = (sig2 * ((X.transpose() * X).inverse()).diagonal()).array().sqrt();
return List::create(Named("coefficients") = coef,
Named("stderr") = stderrest);
}')
fastLm2(y, X)