Ich begann zu spielen mit Rcpp
und möchte die fastLm
Funktion als ein Beispiel verwenden (auch weil es für mögliche spätere Arbeit nützlich ist). Ich weiß, dass fastLm
Teil des RcppArmadillo
Pakets ist, aber ich möchte es unter Verwendung sourceCpp
kompilieren. Der Code kann here gefunden werden und ist auch unten. Das erste Problem, auf das ich stoße, ist, dass ich sourceCpp("fastLm.cpp")
in R nach der Installation und dem Laden Rcpp
und RcppArmadillo
nicht einfach ausführen kann. Ich bekomme diesen Fehler error: RcppArmadillo.h: No such file or directory
und dann alle möglichen Dinge, die ich daraus ergebe.Mit `sourceCpp` kompilieren` fastLm`
Das zweite Problem ist, dass ich denke, ich muss einige Sachen in der fastLm.cpp
ändern. Meine Änderungen sind auch unten, aber ich bin sicher, dass etwas fehlt oder falsch ist. Ich integrierte #include <Rcpp.h>
und using namespace Rcpp;
und // [[Rcpp::export]]
, um die Funktion zu R zu exportieren, und ich änderte die Argumente von zu NumericVector
und NumericMatrix
. Ich sehe nicht, warum das nicht funktionieren sollte und eine ähnliche Anpassung ist wahrscheinlich für den Rückgabewert möglich?
fastLm.cpp
#include <RcppArmadillo.h>
extern "C" SEXP fastLm(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 = arma::solve(X, y); // fit model y ~ X
arma::colvec resid = y - X*coef; // residuals
double sig2 = arma::as_scalar(arma::trans(resid)*resid/(n-k));
// std.error of estimate
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
) ;
}
fastLm.cpp geändert
#include <Rcpp.h>
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::export]]
extern "C" SEXP fastLm(NumericVector yr, NumericMatrix Xr) {
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 = arma::solve(X, y); // fit model y ~ X
arma::colvec resid = y - X*coef; // residuals
double sig2 = arma::as_scalar(arma::trans(resid)*resid/(n-k));
// std.error of estimate
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
) ;
}
Schön. Das einzige, was jetzt fehlt, ist der try/catch-Block, den das Beispiel des OP ebenfalls weggelassen hat (aber das ist in der src). –