Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- sourceCpp(code = '#include <Rcpp.h>
- using namespace Rcpp;
- // [[Rcpp::export]]
- DataFrame mypnorm(NumericVector x){
- int n = x.size();
- NumericVector y1(n), y2(n), y3(n);
- for (int i=0; i<n; i++){
- y1[i] = ::Rf_pnorm5(x[i], 0.0, 1.0, 1, 0);
- y2[i] = R::pnorm(x[i], 0.0, 1.0, 1, 0);
- }
- y3 = pnorm(x);
- return DataFrame::create(Named("R") = y1,
- Named("Rf_") = y2,
- Named("sugar") = y3);
- }')
- mypnorm(runif(10, -3, 3))
- sourceCpp(code = '
- // [[Rcpp::depends(RcppArmadillo)]]
- #include <RcppArmadillo.h>
- 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);
- arma::colvec y(yr.begin(), yr.size(), false);
- arma::colvec coef = arma::solve(X, y);
- arma::colvec resid = y - X*coef;
- double sig2 = arma::as_scalar(arma::trans(resid)*resid/(n-k));
- arma::colvec stderrest = arma::sqrt(sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) );
- return List::create(Named("coefficients") = coef,
- Named("stderr") = stderrest);
- }')
- X = matrix(rnorm(100*2), ncol = 2)
- y = 0.5 + X %*% c(0.3, -0.7) + rnorm(100)
- fastLm(y, X)
- sourceCpp(code = '// [[Rcpp::depends(RcppEigen)]]
- #include <RcppEigen.h>
- using namespace Rcpp;
- using Eigen::Map;
- using Eigen::MatrixXd;
- using Eigen::VectorXd;
- // [[Rcpp::export]]
- List fastLm2(NumericVector yr, NumericMatrix Xr) {
- const Map<MatrixXd> X(as<Map<MatrixXd> >(Xr));
- const Map<VectorXd> y(as<Map<VectorXd> >(yr));
- int n = Xr.nrow(), k = Xr.ncol();
- VectorXd coef = (X.transpose() * X).llt().solve(X.transpose() * y.col(0));
- VectorXd resid = y - X*coef;
- double sig2 = resid.squaredNorm() / (n - k);
- VectorXd stderrest = (sig2 * ((X.transpose() * X).inverse()).diagonal()).array().sqrt();
- return List::create(Named("coefficients") = coef,
- Named("stderr") = stderrest);
- }')
- fastLm2(y, X)
Advertisement
Add Comment
Please, Sign In to add comment