Как быстро вычислять процентиль в C++/Rcpp

Question

У меня есть большой вектор, содержащий кучу двойных элементов. Учитывая массив вектора процентиля, например 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; 
} 

Я также попытаться вызвать функцию R quantile(x, c(.90, .91, .92, .93, .94, .95)) в Rcpp с помощью:

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 

I ожидайте вычисления скоростной процентной доли, но я полагаю, std::sort(tmp.begin(), tmp.end()) замедляет скорость. Есть ли лучший способ получить быстрый результат с использованием C++, RCpp/RcppAramdillo? Благодарю.

Answer 1

Ветвление в петлю может быть точно оптимизировано. Используйте вызовы std :: min/max с помощью int.

Я бы решить Percent расчета индексов массива таким образом:

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

Только эта линия в середине цикла выше:

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

Answer 2

В зависимости от того, сколько процентили вы должны рассчитать и насколько велики ваши векторы, вы можете сделать намного лучше (только O (N)), чем сортировать весь вектор (в лучшем случае O (N * log (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; 
} 

Другой альтернативой является IQAgent Алгоритм из численного рецепты 3. Издание Первоначально он предназначался для потоков данных, но вы можете обмануть его, разделив ваш большой datavector на более мелкие куски (например, 10K элементов) и вычислить процентили для каждого из блоков (где используется сортировка на 10K кусках). Если вы обрабатываете блоки по одному, каждый последующий блок немного изменит значения процентилей, пока вы не получите довольно хорошее приближение в конце. Алгоритм дал хорошие результаты (до 3-го или 4-го десятичного числа), но все еще медленнее, чем реализация n-го элемента.