Come eseguire il calcolo rapido del percentile in C++/Rcpp

Question

Ho un vettore grande contenente un gruppo di elementi doppi. Data una matrice di vettore percentile, ad esempio percentile_vec = c(0.90, 0.91, 0.92, 0.93, 0.94, 0.95). Attualmente sto usando la funzione Rcpp sort per ordinare il vettore grande e quindi trovare il corrispondente valore percentile. Ecco i codici principali:

// [[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; 
} 

Inoltre provo a chiamare funzione R quantile(x, c(.90, .91, .92, .93, .94, .95)) in Rcpp utilizzando:

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

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

Il test riposa per x=runif(1E6) sono elencati di seguito:

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 si aspetta un calcolo rapido del percentile, ma suppongo che std::sort(tmp.begin(), tmp.end()) rallenti la velocità. C'è un modo migliore per ottenere un risultato veloce usando C++, RCpp/RcppAramdillo? Grazie.

Answer 1

La ramificazione in un ciclo potrebbe essere sicuramente ottimizzata. Utilizza le chiamate std :: min/max con ints.

vorrei risolvere cento calcolo degli indici degli array in questo modo:

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

Solo questa riga nel mezzo del ciclo precedente:

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

Answer 2

A seconda di quanti percentili si deve calcolare e quanto sono grandi i tuoi vettori, puoi fare molto meglio (solo O (N)) che ordinare l'intero vettore (nel migliore dei casi O (N * log (N))).

ho dovuto calcolare 1 percentile di vettori (> = 160K) elementi così quello che ho fatto è stata la seguente:

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; 
} 

Un'altra alternativa è l'algoritmo IQAgent da Numerica Ricette 3 °. Ed. Inizialmente era previsto per i flussi di dati, ma è possibile ingannarlo suddividendo il tuo grande datavector in blocchi più piccoli (ad esempio 10K elementi) e calcolare percentili per ciascuno dei blocchi (dove viene utilizzato un ordinamento sui blocchi 10K). Se elaborate i blocchi uno alla volta, ciascun blocco successivo modificherà leggermente i valori dei percentili, fino ad ottenere un'approssimazione abbastanza buona alla fine. L'algoritmo ha dato buoni risultati (fino al 3 ° o al 4 ° decimale), ma era ancora più lento rispetto all'implementazione dell'elemento n-esimo.