对Rcpp NumericMatrix的列进行排序以进行中值计算

Question

我一直在测试Rcpp和RcppArmadillo来计算大矩阵的汇总统计数据.这比基础R colMeans或犰狳在约400万行,45列上快得多(快5或10倍).

colMeansRcpp <- cxxfunction(signature(X_="integer"), 
                            plugin='Rcpp',
                            body='
                            Rcpp::IntegerMatrix X = X_;
                            int ncol = X.ncol(); int nrow = X.nrow();                      
                            Rcpp::NumericVector out(ncol);
                            for(int col = 0; col < ncol; col++){
                              out[col]=Rcpp::sum(X(_, col));
                            }                             
                            return wrap(out/nrow);
                          ')

我真的想计算绘图的中位数和其他分位数 - 因为它需要排序它更需要C++外包.犰狳似乎有点慢,所以我想做一个类似于上面的代码排序,但我不能正确的语法...这是我正在尝试..

# OK I'm aware this floor(nrow/2) is not **absolutely** correct 
# I'm simplifying here
    colMedianRcpp <- cxxfunction(signature(X_="integer"), 
                          plugin='Rcpp',
                          body='
                          Rcpp::IntegerMatrix X = clone(X_);
                          int ncol = X.ncol(); int nrow = X.nrow();                           
                          Rcpp::NumericVector out(ncol);
                          for(int col = 0; col < ncol; col++){
                          X(_,col)= std::sort((X_,col).begin, (X_,col).end));
                          out[col]=X(floor(nrow/2), col));
                          }
                        return wrap(out);
                        ')

基本上就是这条线

X(_,col)= std::sort((X_,col).begin, (X_,col).end));

我不知道如何使用Rcpp sugar和std C++的混合物来表达"对柱子进行分类".对不起,我可以看到我正在做的是错的,但正确的语法提示将是可爱的.

ps我是对的我需要做这个clone()所以我不改变R对象？

编辑 我添加RcppArmadillo代码和基准比较来解决下面的答案/评论.基准测试只有50k行才能快速回复,但我记得它与更多类似.我意识到你是Rcpp的作者..非常感谢你的时间!

我想这可能是因为我正在用RcppArmadillo代码做一些愚蠢的事情,使它比基本的colMeans或Rcpp版本运行得慢得多？

colMeansRcppArmadillo <- cxxfunction(signature(X_="integer"), 
                                     plugin="RcppArmadillo",
                                      body='
                                      arma::mat X = Rcpp::as<arma::mat > (X_);
                                      arma::rowvec MD= arma::mean(X, 0);
                                      return wrap(MD);
                                    ')

基准是......

(mb = microbenchmark(
+                     colMeans(fqSmallMatrix), 
+                     colMeansRcpp(fqSmallMatrix), 
+                     colMeansRcppArmadillo(fqSmallMatrix),
+                     times=50))
Unit: milliseconds
                                 expr       min       lq    median        uq        max neval
              colMeans(fqSmallMatrix) 10.620919 10.63289 10.640819 10.648882  10.907145    50
          colMeansRcpp(fqSmallMatrix)  2.649038  2.66832  2.676709  2.700839   2.841012    50
 colMeansRcppArmadillo(fqSmallMatrix) 25.687067 26.23488 33.168589 33.792489 113.832495    50

Answer 1

您可以将列复制到新的矢量中

NumericVector y = x(_,j);

完整的例子:

library(Rcpp)
cppFunction('
  NumericVector colMedianRcpp(NumericMatrix x) {
    int nrow = x.nrow();
    int ncol = x.ncol();
    int position = nrow / 2; // Euclidian division
    NumericVector out(ncol);
    for (int j = 0; j < ncol; j++) {
      NumericVector y = x(_,j); // Copy the column -- the original will not be modified
      std::nth_element(y.begin(), y.begin() + position, y.end());
      out[j] = y[position];
    }
    return out;
  }
')
x <- matrix( sample(1:12), 3, 4 )
x
colMedianRcpp(x)
x   # Unchanged