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 | e...@debian.org | http://dirk.eddelbuettel.com
_______________________________________________
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
