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

Posted 2023-02-16

【中文标题】如何在 C++/Rcpp 中进行快速百分位数计算

【英文标题】：How to do fast percentile calculation in C++/Rcpp

【发布时间】：2015-05-19 12:41:59

【问题描述】：

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

// [[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 获得快速结果？谢谢。

【问题讨论】：

您可能已经意识到这一点，但这些函数产生的结果略有不同。

好吧排序将是 O(n log(n)) 并且你不能比排序向量更好。之后您正在对向量进行线性搜索以找到相应的元素。你可能会从 binary search 中受益，因为你有一个排序的向量。

@nurssell 你说的完全正确，我也很好奇R是怎么做percentile计算的。我注意到对于runif(1E6)，两个结果略有不同，在我的容忍范围内。

@NathanOliver 感谢您的意见。会看看

@Alvin 我相信this 是base R 的quantile 函数的实现。

【参考方案1】：

循环中的分支肯定可以优化。使用带整数的 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() ) ];

【参考方案2】：

根据您必须计算的百分位数和向量的大小，您可以比对整个向量进行排序（最多 O(N*log(N))）做得更好（仅 O(N)）。

我必须计算 1 个百分位数的向量 (>=160K) 元素，所以我做了以下事情：

void prctile_stl(double* in, const dim_t &len, const double &percent, std::vector<double> &range) 
// Calculates "percent" percentile.
// Linear interpolation inspired by prctile.m from MATLAB.

double r = (percent / 100.) * len;

double lower = 0;
double upper = 0;
double* min_ptr = NULL;
dim_t k = 0;

if(r >= len / 2.)      // Second half is smaller
    dim_t idx_lo = max(r - 1, (double) 0.);
    nth_element(in, in + idx_lo, in + len);             // Complexity O(N)
    lower = in[idx_lo];
    if(idx_lo < len - 1) 
        min_ptr = min_element(&(in[idx_lo + 1]), in + len);
        upper = *min_ptr;
        
    else
        upper = lower;
    
else                   // First half is smaller
    double* max_ptr;
    dim_t idx_up = ceil(max(r - 1, (double) 0.));
    nth_element(in, in + idx_up, in + len);             // Complexity O(N)
    upper = in[idx_up];
    if(idx_up > 0) 
        max_ptr = max_element(in, in + idx_up);
        lower = *max_ptr;
        
    else
        lower = upper;
    

// Linear interpolation
k = r + 0.5;        // Implicit floor
r = r - k;
range[1] = (0.5 - r) * lower + (0.5 + r) * upper;

min_ptr = min_element(in, in + len);
range[0] = *min_ptr;

另一种选择是来自 Numerical Recepies 3rd 的 IQAgent 算法。埃德。 它最初是为数据流设计的，但你可以通过将大数据向量分成更小的块（例如 10K 元素）并计算每个块的百分位数（使用 10K 块的排序）来欺骗它。如果你一次处理一个块，每个连续的块都会稍微修改百分位数的值，直到最后你得到一个很好的近似值。该算法给出了很好的结果（最高到小数点后第 3 位或第 4 位），但仍然比第 n 个元素的实现要慢。