Solo per lo scopo di lavorare sulla mia programmazione C++/Rcpp, ho preso una decisione sull'implementazione di una funzione di deviazione standard (campione) :Prestazioni di statistiche R :: sd() vs. arma :: stddev() rispetto all'implementazione di Rcpp
#include <Rcpp.h>
#include <vector>
#include <cmath>
#include <numeric>
// [[Rcpp::export]]
double cppSD(Rcpp::NumericVector rinVec)
{
std::vector<double> inVec(rinVec.begin(),rinVec.end());
int n = inVec.size();
double sum = std::accumulate(inVec.begin(), inVec.end(), 0.0);
double mean = sum/inVec.size();
for(std::vector<double>::iterator iter = inVec.begin();
iter != inVec.end(); ++iter){
double temp;
temp= (*iter - mean)*(*iter - mean);
*iter = temp;
}
double sd = std::accumulate(inVec.begin(), inVec.end(), 0.0);
return std::sqrt(sd/(n-1));
}
ho anche deciso di testare la funzione stddev dalla libreria Armadillo, considerando che può essere chiamato un vettore:
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
// [[Rcpp::export]]
double armaSD(arma::colvec inVec)
{
return arma::stddev(inVec);
}
Poi Benchmarked queste due funzioni contro la funzione di base R sd() per alcuni vettori di dimensioni variabili:
Rcpp::sourceCpp('G:/CPP/armaSD.cpp')
Rcpp::sourceCpp('G:/CPP/cppSD.cpp')
require(microbenchmark)
##
## sample size = 1,000: armaSD() < cppSD() < sd()
X <- rexp(1000)
microbenchmark(armaSD(X),sd(X), cppSD(X))
#Unit: microseconds
# expr min lq median uq max neval
# armaSD(X) 4.181 4.562 4.942 5.322 12.924 100
# sd(X) 17.865 19.766 20.526 21.287 86.285 100
# cppSD(X) 4.561 4.941 5.321 5.701 29.269 100
##
## sample size = 10,000: armaSD() < cppSD() < sd()
X <- rexp(10000)
microbenchmark(armaSD(X),sd(X), cppSD(X))
#Unit: microseconds
# expr min lq median uq max neval
# armaSD(X) 24.707 25.847 26.4175 29.6490 52.455 100
# sd(X) 51.315 54.356 55.8760 61.1980 100.730 100
# cppSD(X) 26.608 28.128 28.8885 31.7395 114.413 100
##
## sample size = 25,000: armaSD() < cppSD() < sd()
X <- rexp(25000)
microbenchmark(armaSD(X),sd(X), cppSD(X))
#Unit: microseconds
# expr min lq median uq max neval
# armaSD(X) 66.900 67.6600 68.040 76.403 155.845 100
# sd(X) 108.332 111.5625 122.016 125.817 169.910 100
# cppSD(X) 70.320 71.0805 74.692 80.203 102.250 100
##
## sample size = 50,000: cppSD() < sd() < armaSD()
X <- rexp(50000)
microbenchmark(armaSD(X),sd(X), cppSD(X))
#Unit: microseconds
# expr min lq median uq max neval
# armaSD(X) 249.733 267.4085 297.8175 337.729 642.388 100
# sd(X) 203.740 229.3975 240.2300 260.186 303.709 100
# cppSD(X) 162.308 185.1140 239.6600 256.575 290.405 100
##
## sample size = 75,000: sd() < cppSD() < armaSD()
X <- rexp(75000)
microbenchmark(armaSD(X),sd(X), cppSD(X))
#Unit: microseconds
# expr min lq median uq max neval
# armaSD(X) 445.110 479.8900 502.5070 520.5625 642.388 100
# sd(X) 310.931 334.8780 354.0735 379.7310 429.146 100
# cppSD(X) 346.661 380.8715 400.6370 424.0140 501.747 100
non ero terribilmente sorpreso per il fatto che il mio C++ funzione cppSD() era più veloce di stats::sd() per i campioni più piccoli, ma più lento per i vettori di dimensioni più grandi in quanto è stats::sd() Vettorializzare. Tuttavia, non mi aspettavo di vedere un tale degrado delle prestazioni dalla funzione arma::stddev() poiché sembra funzionare anche in modo vettorializzato. C'è un problema con il modo in cui sto usando arma::stdev(), o è semplicemente che stats::sd() è stato scritto (in C sto assumendo) in modo tale da poter gestire i vettori più grandi in modo molto più efficiente? Qualsiasi input sarebbe apprezzato.
Aggiornamento:
Anche se la mia domanda era in origine circa l'uso corretto dei arma::stddev, e non tanto di cercare trovare il modo più efficiente possibile calcolare la deviazione standard del campione, è molto interessante per vedere che la funzione zucchero Rcpp::sd si comporta così bene. Per rendere le cose un po 'più interessante, ho benchmark le arma::stddev e Rcpp::sd funzioni basso contro una versione RcppParallel che ho adattato da due dei messaggi Rcpp Galleria di JJ Allaire - here e here:
library(microbenchmark)
set.seed(123)
x <- rnorm(5.5e06)
##
Res <- microbenchmark(
armaSD(x),
par_sd(x),
sd_sugar(x),
times=500L,
control=list(warmup=25))
##
R> print(Res)
Unit: milliseconds
expr min lq mean median uq max neval
armaSD(x) 24.486943 24.960966 26.994684 25.255584 25.874139 123.55804 500
par_sd(x) 8.130751 8.322682 9.136323 8.429887 8.624072 22.77712 500
sd_sugar(x) 13.713366 13.984638 14.628911 14.156142 14.401138 32.81684 500
Questo è stato il mio portatile in esecuzione 64 -bit linux con processore i5-4200U a 1.60GHz; ma suppongo che la differenza tra par_sd e sugar_sd sia meno sostanziale su una macchina Windows.
E il codice per la versione RcppParallel (che è notevolmente più lungo, e richiede un compilatore compatibile C++ 11 per l'espressione lambda utilizzato in sovraccarico operator() del InnerProduct funtore):
#include <Rcpp.h>
#include <RcppParallel.h>
// [[Rcpp::depends(RcppParallel)]]
// [[Rcpp::plugins(cpp11)]]
/*
* based on: http://gallery.rcpp.org/articles/parallel-vector-sum/
*/
struct Sum : public RcppParallel::Worker {
const RcppParallel::RVector<double> input;
double value;
Sum(const Rcpp::NumericVector input)
: input(input), value(0) {}
Sum(const Sum& sum, RcppParallel::Split)
: input(sum.input), value(0) {}
void operator()(std::size_t begin, std::size_t end) {
value += std::accumulate(input.begin() + begin,
input.begin() + end,
0.0);
}
void join(const Sum& rhs) {
value += rhs.value;
}
};
/*
* based on: http://gallery.rcpp.org/articles/parallel-inner-product/
*/
struct InnerProduct : public RcppParallel::Worker {
const RcppParallel::RVector<double> x;
const RcppParallel::RVector<double> y;
double mean;
double product;
InnerProduct(const Rcpp::NumericVector x,
const Rcpp::NumericVector y,
const double mean)
: x(x), y(y), mean(mean), product(0) {}
InnerProduct(const InnerProduct& innerProduct,
RcppParallel::Split)
: x(innerProduct.x), y(innerProduct.y),
mean(innerProduct.mean), product(0) {}
void operator()(std::size_t begin, std::size_t end) {
product += std::inner_product(x.begin() + begin,
x.begin() + end,
y.begin() + begin,
0.0, std::plus<double>(),
[&](double lhs, double rhs)->double {
return ((lhs-mean)*(rhs-mean));
});
}
void join(const InnerProduct& rhs) {
product += rhs.product;
}
};
// [[Rcpp::export]]
double par_sd(const Rcpp::NumericVector& x_)
{
int N = x_.size();
Rcpp::NumericVector y_(x_);
Sum sum(x_);
RcppParallel::parallelReduce(0, x_.length(), sum);
double mean = sum.value/N;
InnerProduct innerProduct(x_, y_, mean);
RcppParallel::parallelReduce(0, x_.length(), innerProduct);
return std::sqrt(innerProduct.product/(N-1));
}
In aggiunta al problema indicato da Dirk Eddelbuettel di seguito, la funzione cppSD() ha ulteriori problemi. Non è un metodo robusto per calcolare la deviazione standard, poiché è facilmente influenzato da overflow in virgola mobile nella variabile sum. std :: accumulate() non è progettato per la numerics. Invece, il tuo codice dovrebbe utilizzare statistiche in esecuzione come metodo di calcolo primario o come fallback. La funzione cppSD() è anche inefficiente: copia inutilmente i dati su un vettore nuovo e quindi riscrive il vettore. Invece, dovrebbe usare i dati originali in modo di sola lettura – mtall
Questo vale ripetere: l'uso di funzioni come std :: accumulate() o std :: inner_product() per calcolare la media e la deviazione standard è una pessima idea: non sono robusti per overflow o underflow in virgola mobile. La funzione di deviazione standard di Armadillo è molto più sicura, in quanto tiene conto di questi potenziali problemi. Usando il metodo di parallelizzazione sopra riportato otterrai risultati cattivi più rapidamente. – mtall