Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- // [[Rcpp::depends(RcppArmadillo)]]
- #include <RcppArmadillo.h>
- // [[Rcpp::depends(RcppParallel)]]
- #include <RcppParallel.h>
- #include <string>
- using namespace arma;
- void chk_mat(const mat& x, const std::string& varName, const std::string& type){
- if (!is_finite(x))
- Rcpp::stop(varName + " must be numerical.\n");
- }
- arma::vec quantileCpp(const arma::vec& x, const arma::vec& probs){
- chk_mat(x, "x", "double");
- chk_mat(probs, "probs", "double");
- if (any(probs < 0) || any(probs > 1))
- Rcpp::stop("There is values in probs outside [0, 1].");
- vec x_sort = sort(x);
- vec index = (x.n_elem - 1) * probs;
- uvec lo = conv_to<uvec>::from(floor(index)), hi = conv_to<uvec>::from(ceil(index));
- vec h = (index - conv_to<vec>::from(lo));
- vec qs = (1 - h) % x_sort.elem(lo) + h % x_sort.elem(hi);
- return qs;
- }
- struct Worker_bs_stat: public RcppParallel::Worker {
- const vec& X;
- const vec& Y;
- const vec& probs;
- mat& bsStat;
- Worker_bs_stat(const vec& X, const vec& Y, const vec& probs, mat& bsStat):
- X(X), Y(Y), probs(probs), bsStat(bsStat) {}
- void operator()(std::size_t begin, std::size_t end)
- {
- for (uword i = begin; i < end; ++i)
- {
- uvec indx_x = randi<uvec>(X.n_elem, distr_param(0, X.n_elem-1));
- uvec indx_y = randi<uvec>(Y.n_elem, distr_param(0, Y.n_elem-1));
- bsStat.col(i) = quantileCpp(X.elem(indx_x), probs) - quantileCpp(Y.elem(indx_y), probs);
- }
- }
- };
- // [[Rcpp::export]]
- arma::vec dist_test(const arma::vec& X, const arma::vec& Y, const arma::vec& probs, const double& bsSize){
- // check data
- chk_mat(X, "X", "double");
- chk_mat(Y, "Y", "double");
- chk_mat(probs, "probs", "double");
- if (bsSize <= 0 || std::abs(bsSize - floor(bsSize)) > 1e-6)
- Rcpp::stop("bsSize must be a positive number.\n");
- mat bsDiff = zeros<mat>(probs.n_elem, bsSize);
- Worker_bs_stat bs_stat(X, Y, probs, bsDiff);
- RcppParallel::parallelFor(0, bsSize, bs_stat);
- vec stat = (quantileCpp(X, probs) - quantileCpp(Y, probs)) / stddev(bsDiff, 0, 1);
- return stat;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement