2012-12-03 7 views
5

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 
    ) ; 

} 

Antwort

9

Sie müssen die Abhängigkeit von RcppArmadillo mit dem Rcpp::depends Pseudo-Attribut, um anzuzeigen. Diese kümmert sich um zu finden RcppArmadillo Kopf- und Link gegen blas, lapack etc ...

#include <RcppArmadillo.h> 
// [[Rcpp::depends(RcppArmadillo)]] 

using namespace Rcpp; 

// [[Rcpp::export]] 
List 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 
    ) ; 

} 

Auch ist es sehr wichtig, dass Sie verwenden #include <RcppArmadillo.h> und nicht #include <Rcpp.h>. RcppArmadillo.h kümmert sich um Rcpp.hzur richtigen Zeit, und Reihenfolge der Include-Dateien ist sehr wichtig hier.

Sie können auch eine List zurückgeben und die extern "C" löschen.

+0

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). –