Guest User

Examples_Rcpp_Attributes

a guest
Jun 30th, 2014
326
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.  
  2.  
  3. sourceCpp(code = '#include <Rcpp.h>
  4. using namespace Rcpp;
  5.  
  6. // [[Rcpp::export]]
  7. DataFrame mypnorm(NumericVector x){
  8.     int n = x.size();
  9.     NumericVector y1(n), y2(n), y3(n);
  10.    
  11.     for (int i=0; i<n; i++){
  12.         y1[i] = ::Rf_pnorm5(x[i], 0.0, 1.0, 1, 0);
  13.         y2[i] = R::pnorm(x[i], 0.0, 1.0, 1, 0);
  14.     }
  15.     y3 = pnorm(x);
  16.     return DataFrame::create(Named("R") = y1,
  17.                  Named("Rf_") = y2,
  18.                  Named("sugar") = y3);
  19. }')
  20. mypnorm(runif(10, -3, 3))
  21.  
  22. sourceCpp(code = '
  23. // [[Rcpp::depends(RcppArmadillo)]]
  24. #include <RcppArmadillo.h>
  25. using namespace Rcpp;
  26. // [[Rcpp::export]]
  27.  
  28. List fastLm(NumericVector yr, NumericMatrix Xr) {
  29.     int n = Xr.nrow(), k = Xr.ncol();
  30.     arma::mat X(Xr.begin(), n, k, false);
  31.     arma::colvec y(yr.begin(), yr.size(), false);
  32.     arma::colvec coef = arma::solve(X, y);
  33.     arma::colvec resid = y - X*coef;
  34.     double sig2 = arma::as_scalar(arma::trans(resid)*resid/(n-k));
  35.     arma::colvec stderrest = arma::sqrt(sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) );
  36.     return List::create(Named("coefficients") = coef,
  37.                            Named("stderr") = stderrest);
  38. }')
  39. X = matrix(rnorm(100*2), ncol = 2)
  40. y = 0.5 + X %*% c(0.3, -0.7) + rnorm(100)
  41. fastLm(y, X)
  42.  
  43. sourceCpp(code = '// [[Rcpp::depends(RcppEigen)]]
  44. #include <RcppEigen.h>
  45. using namespace Rcpp;
  46. using Eigen::Map;
  47. using Eigen::MatrixXd;
  48. using Eigen::VectorXd;
  49. // [[Rcpp::export]]
  50.  
  51. List fastLm2(NumericVector yr, NumericMatrix Xr) {
  52.     const Map<MatrixXd> X(as<Map<MatrixXd> >(Xr));
  53.     const Map<VectorXd> y(as<Map<VectorXd> >(yr));
  54.     int n = Xr.nrow(), k = Xr.ncol();
  55.     VectorXd coef = (X.transpose() * X).llt().solve(X.transpose() * y.col(0));
  56.     VectorXd resid = y - X*coef;
  57.     double sig2 = resid.squaredNorm() / (n - k);
  58.     VectorXd stderrest = (sig2 * ((X.transpose() * X).inverse()).diagonal()).array().sqrt();
  59.     return List::create(Named("coefficients") = coef,
  60.                 Named("stderr") = stderrest);
  61. }')
  62. fastLm2(y, X)
Advertisement
Add Comment
Please, Sign In to add comment