SHOW:
|
|
- or go back to the newest paste.
| 1 | sourceCpp(code = ' | |
| 2 | // [[Rcpp::depends(RcppArmadillo)]] | |
| 3 | #include <RcppArmadillo.h> | |
| 4 | using namespace Rcpp; | |
| 5 | // [[Rcpp::export]] | |
| 6 | List fastLm(NumericVector yr, NumericMatrix Xr) {
| |
| 7 | int n = Xr.nrow(), k = Xr.ncol(); | |
| 8 | arma::mat X(Xr.begin(), n, k, false); | |
| 9 | arma::colvec y(yr.begin(), yr.size(), false); | |
| 10 | arma::colvec coef = arma::solve(X, y); | |
| 11 | arma::colvec resid = y - X*coef; | |
| 12 | double sig2 = arma::as_scalar(arma::trans(resid)*resid/(n- | |
| 13 | k)); | |
| 14 | arma::colvec stderrest = arma::sqrt( | |
| 15 | sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) ); | |
| 16 | return List::create(Named("coefficients") = coef,
| |
| 17 | Named("stderr") = stderrest);
| |
| 18 | }') | |
| 19 | X = matrix(rnorm(100*2), ncol = 2) | |
| 20 | y = 0.5 + X %*% c(0.3, -0.7) + rnorm(100) | |
| 21 | fastLm(y, X) | |
| 22 | ||
| 23 | sourceCpp(code = '// [[Rcpp::depends(RcppEigen)]] | |
| 24 | #include <RcppEigen.h> | |
| 25 | using namespace Rcpp; | |
| 26 | using Eigen::Map; | |
| 27 | using Eigen::MatrixXd; | |
| 28 | using Eigen::VectorXd; | |
| 29 | // [[Rcpp::export]] | |
| 30 | List fastLm2(NumericVector yr, NumericMatrix Xr) {
| |
| 31 | const Map<MatrixXd> X(as<Map<MatrixXd> >(Xr)); | |
| 32 | const Map<VectorXd> y(as<Map<VectorXd> >(yr)); | |
| 33 | int n = Xr.nrow(), k = Xr.ncol(); | |
| 34 | VectorXd coef = (X.transpose() * X).llt().solve(X.transpose() * y.col(0)); | |
| 35 | VectorXd resid = y - X*coef; | |
| 36 | double sig2 = resid.squaredNorm() / (n - k); | |
| 37 | VectorXd stderrest = (sig2 * ((X.transpose() * X).inverse()).diagonal()).array().sqrt(); | |
| 38 | return List::create(Named("coefficients") = coef,
| |
| 39 | Named("stderr") = stderrest);
| |
| 40 | }') | |
| 41 | fastLm2(y, X) |
