c++ - 如何在 C++/Rcpp 中进行快速百分位数计算

Question

我有一个包含一堆 double 元素的大 vector 。给定一个百分位数 vector 数组，例如 percentile_vec = c(0.90, 0.91, 0.92, 0.93, 0.94, 0.95)。我目前正在使用 Rcpp sort 函数对大 vector 进行排序，然后找到相应的百分位值。主要代码如下:

// [[Rcpp::export]]
NumericVector sort_rcpp(Rcpp::NumericVector& x)
{
  std::vector<double> tmp = Rcpp::as<std::vector<double>> (x);    // or NumericVector tmp = clone(x);
  std::sort(tmp.begin(), tmp.end());
  return wrap(tmp);
}

// [[Rcpp::export]]
NumericVector percentile_rcpp(Rcpp::NumericVector& x, Rcpp::NumericVector& percentile)
{
  NumericVector tmp_sort = sort_rcpp(x);
  int size_per = percentile.size();
  NumericVector percentile_vec = no_init(size_per);
  for (int ii = 0; ii < size_per; ii++)
  {
    double size_per = tmp_sort.size() * percentile[ii];
    double size_per_round;
    if (size_per < 1.0)
    {
      size_per_round = 1.0;
    }
    else
    {
      size_per_round = std::round(size_per);
    }
    percentile_vec[ii] = tmp_sort[size_per_round-1];  // For extreme case such as size_per_round == tmp_sort.size() to avoid overflow
  }
  return percentile_vec;
}

我还尝试在 Rcpp 中调用 R 函数 quantile(x, c(.90, .91, .92, .93, .94, .95)) 使用:

sub_percentile <- function (x)
{
  return (quantile(x, c(.90, .91, .92, .93, .94, .95)));
}  

source('C:/Users/~Call_R_function.R')

下面列出了 x=runif(1E6) 的测试休息:

microbenchmark(sub_percentile(x)->aa, percentile_rcpp(x, c(.90, .91, .92, .93, .94, .95))->bb)
#Unit: milliseconds
              expr      min       lq     mean   median       uq       max   neval
  sub_percentile(x) 99.00029 99.24160 99.35339 99.32162 99.41869 100.57160   100
 percentile_rcpp(~) 87.13393 87.30904 87.44847 87.40826 87.51547  88.41893   100

我期待一个快速的百分位数计算，但我假设 std::sort(tmp.begin(), tmp.end()) 会降低速度。有没有更好的方法可以使用 C++、RCpp/RcppAramdillo 快速获得结果？谢谢。

Answer 1

肯定可以优化循环中的分支。对 int 使用 std::min/max 调用。

我会用这种方式解决数组索引的百分比计算:

uint PerCentIndex( double pc, uint size )
{
    return 0.5 + ( double ) ( size - 1 ) * pc;
}

只有上面循环中间的这一行:

percentile_vec[ii] 
 = tmp_sort[ PerCentIndex( percentile[ii], tmp_sort.size() ) ];