只是为了锻炼我的C++ 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) );
}
我还决定测试一下 stddev 函数，因为它可以在向量上调用。
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
// [[Rcpp::export]]
double armaSD(arma::colvec inVec)
{
return arma::stddev(inVec);
}
然后，我将这两个函数与基本的R函数进行了基准测试 sd() 为几个大小不一的向量。
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
我并没有感到非常惊讶，因为我的C++函数是： cppSD() 快于 stats::sd() 对于较小的样本，但对于较大尺寸的向量，速度较慢，因为 stats::sd() 是矢量化的。然而，我没有想到会出现这样的性能下降。arma::stddev() 函数，因为它似乎也是以矢量化的方式运行。是否我使用的 arma::stdev()还是说 stats::sd() 已写好 C 我假设）这样的方式，它可以更有效地处理更大的向量？如果有任何意见，我们将不胜感激。
更新:
虽然我的问题原本是关于正确使用 "A "的问题。arma::stddev而不是试图找到最有效的方法来计算样本标准差，非常有趣的是，我们可以看到。Rcpp::sd 糖功能表现得如此出色。为了让事情变得更有趣，我对以下功能进行了基准测试 arma::stddev 和 Rcpp::sd 下面的功能对一个 RcppParallel 我从JJ Allaire的两篇Rcpp Gallery文章中改编的版本----。此处 和 此处:
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
这是在我的笔记本上，运行64位linux，CPU为i5-4200U @ 1.60GHz的处理器；但我猜测两者之间的区别在于 par_sd 和 sugar_sd 在Windows机器上将不那么可观。
而且 RcppParallel 的版本（该版本相当长，并且需要一个兼容C++11的编译器来编译重载的 operator() 的 InnerProduct functor）。)
#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) );
}
你在实例化Armadillo对象的方式上犯了一个微妙的错误--这会导致复制，从而降低性能。
使用一个 const arma::colvec & invec 取而代之，一切都好办了。
R> sourceCpp("/tmp/sd.cpp")
R> library(microbenchmark)
R> X <- rexp(500)
R> microbenchmark(armaSD(X), armaSD2(X), sd(X), cppSD(X))
Unit: microseconds
expr min lq median uq max neval
armaSD(X) 3.745 4.0280 4.2055 4.5510 19.375 100
armaSD2(X) 3.305 3.4925 3.6400 3.9525 5.154 100
sd(X) 22.463 23.6985 25.1525 26.0055 52.457 100
cppSD(X) 3.640 3.9495 4.2030 4.8620 13.609 100
R> X <- rexp(5000)
R> microbenchmark(armaSD(X), armaSD2(X), sd(X), cppSD(X))
Unit: microseconds
expr min lq median uq max neval
armaSD(X) 18.627 18.9120 19.3245 20.2150 34.684 100
armaSD2(X) 14.583 14.9020 15.1675 15.5775 22.527 100
sd(X) 54.507 58.8315 59.8615 60.4250 84.857 100
cppSD(X) 18.585 19.0290 19.3970 20.5160 22.174 100
R> X <- rexp(50000)
R> microbenchmark(armaSD(X), armaSD2(X), sd(X), cppSD(X))
Unit: microseconds
expr min lq median uq max neval
armaSD(X) 186.307 187.180 188.575 191.825 405.775 100
armaSD2(X) 142.447 142.793 143.207 144.233 155.770 100
sd(X) 382.857 384.704 385.223 386.075 405.713 100
cppSD(X) 181.601 181.895 182.279 183.350 194.588 100
R>
这是基于我的版本，你的代码，所有的东西都是一个文件，而且... ... armaSD2 正如我所建议的那样 -- -- 导致胜利的表现。
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
#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 = (*iter - mean)*(*iter - mean);
*iter = temp;
}
double sd = std::accumulate(inVec.begin(), inVec.end(), 0.0);
return std::sqrt( sd / (n-1) );
}
// [[Rcpp::export]]
double armaSD(arma::colvec inVec) {
return arma::stddev(inVec);
}
// [[Rcpp::export]]
double armaSD2(const arma::colvec & inVec) { return arma::stddev(inVec); }
/*** R
library(microbenchmark)
X <- rexp(500)
microbenchmark(armaSD(X), armaSD2(X), sd(X), cppSD(X))
X <- rexp(5000)
microbenchmark(armaSD(X), armaSD2(X), sd(X), cppSD(X))
X <- rexp(50000)
microbenchmark(armaSD(X), armaSD2(X), sd(X), cppSD(X))
*/
我认为。sd Rcpp糖中内置的函数效率更高。请看下面的代码。
#include <RcppArmadillo.h>
//[[Rcpp::depends(RcppArmadillo)]]
#include <vector>
#include <cmath>
#include <numeric>
using namespace Rcpp;
//[[Rcpp::export]]
double sd_cpp(NumericVector& xin){
std::vector<double> xres(xin.begin(),xin.end());
int n=xres.size();
double sum=std::accumulate(xres.begin(),xres.end(),0.0);
double mean=sum/n;
for(std::vector<double>::iterator iter=xres.begin();iter!=xres.end();++iter){
double tmp=(*iter-mean)*(*iter-mean);
*iter=tmp;
}
double sd=std::accumulate(xres.begin(),xres.end(),0.0);
return std::sqrt(sd/(n-1));
}
//[[Rcpp::export]]
double sd_arma(arma::colvec& xin){
return arma::stddev(xin);
}
//[[Rcpp::export]]
double sd_sugar(NumericVector& xin){
return sd(xin);
}
> sourcecpp("sd.cpp")
> microbenchmark(sd(X),sd_cpp(X),sd_arma(X),sd_sugar(X))
Unit: microseconds
expr min lq mean median uq max neval
sd(X) 47.655 49.4120 51.88204 50.5395 51.1950 113.643 100
sd_cpp(X) 28.145 28.4410 29.01541 28.6695 29.4570 37.118 100
sd_arma(X) 23.706 23.9615 24.65931 24.1955 24.9520 50.375 100
sd_sugar(X) 19.197 19.478 20.38872 20.0785 21.2015 28.664 100
关于Rcpp函数的2个问题。
1) 你有多大可能需要标准差而不需要平均值？如果不需要标准差而不需要平均值的情况不常见，为什么不同时返回这两个值，而不是要求R用户进行2次函数调用，这实际上是计算了两次平均值。
2）函数内部克隆向量有什么原因吗？
using namespace Rcpp;
// [[Rcpp::plugins(cpp14)]]
// [[Rcpp::export]]
NumericVector cppSD(NumericVector rinVec)
{
double n = (double)rinVec.size();
double sum = 0;
for (double& v : rinVec)
sum += v;
double mean = sum / n;
double varianceNumerator = 0;
for(double& v : rinVec)
varianceNumerator += (v - mean) * (v - mean);
return NumericVector::create(_["std.dev"] = sqrt(varianceNumerator/ (n - 1.0)),
_["mean"] = mean);
}