[Rcpp-devel] [rcpp-devel] arma::mat constructor vs Rcpp::as and interaction with .insert_col

[Dirk Eddelbuettel]

Dale,

Here is a simple fix:  just _join_ the matrix and the col vector.
Result first (now with set.seed(42), so rngs different)

R> sourceCpp("/tmp/dale.cpp")

R> set.seed(42)  ## needed !!

R> mat <- cbind(rnorm(50), rnorm(50), rnorm(50))

R> y <- rnorm(50)

R> fastLm2(y, mat) 
$coefficients
           [,1]
[1,] -0.0227235
[2,]  0.0858994
[3,] -0.0457201
[4,] -0.0441346

$stderr
         [,1]
[1,] 0.115021
[2,] 0.141599
[3,] 0.136838


R> lm(y ~ mat) 

Call:
lm(formula = y ~ mat)

Coefficients:
(Intercept)         mat1         mat2         mat3  
    -0.0227       0.0859      -0.0457      -0.0441  


R> 


The key part is to use join_rows() and to actually use the new matrix:

  arma::colvec ones = arma::ones<arma::colvec>(X.n_rows);
  arma::mat X2 = arma::join_rows(ones, X);

  arma::colvec coef = arma::solve(X2, y);
  arma::colvec resid = y - X2*coef;


Full function below. 

Hth,  Dirk



// [[Rcpp::export]]
Rcpp::List fastLm2(Rcpp::NumericVector yr, Rcpp::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 ones = arma::ones<arma::colvec>(X.n_rows);
  arma::mat X2 = arma::join_rows(ones, X);

  arma::colvec coef = arma::solve(X2, y);
  arma::colvec resid = y - X2*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 Rcpp::List::create(Rcpp::Named("coefficients") = coef,
			    Rcpp::Named("stderr")       = stderrest);
}



-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com