require( inline ) require( Rcpp ) cxxflags <- capture.output( RcppArmadillo:::CxxFlags() ) ldflags <- capture.output( Rcpp:::LdFlags() ) code <- ' #include extern "C" { SEXP file10d63af1 ( SEXP ys, SEXP Xs ); } SEXP file10d63af1 ( SEXP ys, SEXP Xs ) { Rcpp::NumericVector yr(ys); // creates Rcpp vector from SEXP Rcpp::NumericMatrix Xr(Xs); // creates Rcpp matrix from SEXP int n = Xr.nrow(), k = Xr.ncol(); arma::mat X(Xr.begin(), n, k, false); // reuses memory and avoids extra copy arma::colvec y(yr.begin(), yr.size(), false); arma::colvec coef = solve(X, y); // fit model y ~ X arma::colvec resid = y - X*coef; // residuals double sig2 = arma::as_scalar( trans(resid)*resid/(n-k) ); // std.error of estimate arma::colvec stderrest = sqrt( sig2 * diagvec( arma::inv(arma::trans(X)*X)) ); Rcpp::Pairlist res(Rcpp::Named( "coefficients", coef), Rcpp::Named( "stderr", stderrest)); return res; } ' writeLines( code, "armaex.cpp" ) Sys.setenv( PKG_LIBS = ldflags ) Sys.setenv( PKG_CXXFLAGS = cxxflags ) unlink( "armaex.so" ) system( "R CMD SHLIB armaex.cpp" ) dyn.load( "armaex.so" ) fx <- function( ... ){ .Call( "file10d63af1", ... ) } fx( log(trees$Volume), cbind(rep(1,31), log(trees$Girth) ) )