Właśnie w celu pracy nad moim programowaniem w C++/Rcpp, zrobiłem ujęcie w implementacji (próbka) standardowej funkcji odchylenia :Wydajność R statystyk :: sd() vs. arma :: stddev() vs. implementacja 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));
}
ja też postanowiłem przetestować funkcję stddev
z biblioteki Armadillo, biorąc pod uwagę, że można go nazwać na wektorze:
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
// [[Rcpp::export]]
double armaSD(arma::colvec inVec)
{
return arma::stddev(inVec);
}
Potem przykładowym te dwie funkcje wobec funkcji baza R sd()
dla kilku wektorów o różnym rozmiarze:
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
nie byłem strasznie zaskoczony faktem, że moje C++ funkcja cppSD()
był szybszy niż stats::sd()
dla mniejszych próbek, ale wolniej dla większych rozmiarów wektorów od stats::sd()
jest wektorowy. Jednak nie spodziewałem się pogorszenia wydajności z funkcji arma::stddev()
, ponieważ wydaje się, że działa również w wektorze. Czy jest jakiś problem ze sposobem, w jaki używam arma::stdev()
, czy jest to po prostu, że został napisany stats::sd()
(w C
) w taki sposób, że może on obsługiwać większe wektory o wiele bardziej wydajnie? Dowolne wejście będzie docenione.
Aktualizacja:
Chociaż moje pytanie było pierwotnie o prawidłowym korzystaniu z arma::stddev
, i nie tyle o próbując znaleźć najbardziej efektywny sposób można obliczyć odchylenie standardowe próbki, to jest bardzo ciekawy aby przekonać się, że funkcja cukrowa Rcpp::sd
działa tak dobrze. Żeby było trochę bardziej interesująca, że porównywana z arma::stddev
i Rcpp::sd
funkcje poniżej przeciwko RcppParallel
wersji, że przystosowanym od dwóch JJ Allaire za RCPP Gallery postów - here i 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
To było na moim laptopie działa 64 -bitowy Linux z procesorem i5-4200U CPU @ 1.60GHz; ale domyślam się, że różnica między par_sd
i sugar_sd
byłaby mniej znacząca na komputerze z systemem Windows.
I kod dla wersji RcppParallel
(co jest znacznie dłuższa i wymaga C++ 11 zgodnego kompilatora dla wyrażenia lambda używane w przeciążony operator()
z InnerProduct
funktora):
#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));
}
Oprócz kwestii wskazanej przez Dirka Eddelbuettela poniżej, twoja funkcja cppSD() ma dalsze problemy. Nie jest to niezawodny sposób obliczania odchylenia standardowego, ponieważ łatwo podlega on przepełnieniom zmiennoprzecinkowym w zmiennej sumarycznej. std :: accumulate() nie jest przeznaczony do liczb. Zamiast tego Twój kod powinien korzystać z bieżących statystyk jako podstawowej metody obliczeniowej lub awaryjnego. Funkcja cppSD() jest również nieefektywna: niepotrzebnie kopiuje dane do nowego wektora, a następnie ponownie zapisuje wektor. Zamiast tego powinien używać oryginalnych danych w trybie tylko do odczytu. – mtall
Powtarzające się: używanie funkcji takich jak std :: accumulate() lub std :: inner_product() do obliczenia średniej i odchylenia standardowego jest naprawdę złym pomysłem: nie są odporne na przepełnienia lub niedomiary zmiennoprzecinkowe. Standardowa funkcja odchylenia standardu Armadillo jest o wiele bezpieczniejsza, ponieważ uwzględnia potencjalne problemy. Używając powyższej metody parowania, otrzymasz szybsze wyniki. – mtall