Fast calculation of quantiles for each column of a matrix using Rcpp

colQuantiles.cpp

#include <Rcpp.h>
using namespace Rcpp;

// q is a vector of quantiles to calculate for each column of x
// [[Rcpp::export]]
NumericMatrix colQuantiles(NumericMatrix x, NumericVector q) {
  int c = x.ncol();
  int r = x.nrow();
  int n = q.length();
  NumericMatrix out(n, c);
  NumericVector y(r);
  IntegerVector z(n);
  
  // calculate positions for quantiles
  for (int i=0; i<n; i++){
    z[i] = std::floor(q(i)*(r-1) + 0.5);
  }

  
  
  // sort all elements
  for (int i=0; i<c; i++){
    y = x(_, i);
    std::sort(y.begin(), y.end());
    
    for (int j=0; j<n; j++){
      out(j, i) = y(z(j));
    }
  }
  return out;
}

/*** R
colQuantiles(matrix(1:12, 3, 4), c(0, .9, 1))
*/