r - 与简单的Rcpp实现相比，为什么zoo::rollmean慢？

zoo::rollmean是一个有用的函数，它返回时间序列的滚动平均值；对于长度为x和窗口大小为n的矢量k，它返回矢量c(mean(x[1:k]), mean(x[2:(k+1)]), ..., mean(x[(n-k+1):n]))。

我注意到对于正在开发的某些代码，它似乎运行缓慢，因此我使用Rcpp包和一个简单的for循环编写了自己的版本:

library(Rcpp)
cppFunction("NumericVector rmRcpp(NumericVector dat, const int window) {
  const int n = dat.size();
  NumericVector ret(n-window+1);
  double summed = 0.0;
  for (int i=0; i < window; ++i) {
    summed += dat[i];
  }
  ret[0] = summed / window;
  for (int i=window; i < n; ++i) {
    summed += dat[i] - dat[i-window];
    ret[i-window+1] = summed / window;
  }
  return ret;
}")

令我惊讶的是，此版本的函数比zoo::rollmean函数快得多:

# Time series with 1000 elements
set.seed(144)
y <- rnorm(1000)
x <- 1:1000
library(zoo)
zoo.dat <- zoo(y, x)

# Make sure our function works
all.equal(as.numeric(rollmean(zoo.dat, 3)), rmRcpp(y, 3))
# [1] TRUE

# Benchmark
library(microbenchmark)
microbenchmark(rollmean(zoo.dat, 3), rmRcpp(y, 3))
# Unit: microseconds
#                  expr     min       lq       mean    median        uq       max neval
#  rollmean(zoo.dat, 3) 685.494 904.7525 1776.88666 1229.2475 1744.0720 15724.321   100
#          rmRcpp(y, 3)   6.638  12.5865   46.41735   19.7245   27.4715  2418.709   100

加速甚至适用于更大的向量:

# Time series with 5 million elements
set.seed(144)
y <- rnorm(5000000)
x <- 1:5000000
library(zoo)
zoo.dat <- zoo(y, x)

# Make sure our function works
all.equal(as.numeric(rollmean(zoo.dat, 3)), rmRcpp(y, 3))
# [1] TRUE

# Benchmark
library(microbenchmark)
microbenchmark(rollmean(zoo.dat, 3), rmRcpp(y, 3), times=10)
# Unit: milliseconds
#                  expr        min         lq       mean     median         uq        max
#  rollmean(zoo.dat, 3) 2825.01622 3090.84353 3191.87945 3206.00357 3318.98129 3616.14047
#          rmRcpp(y, 3)   31.03014   39.13862   42.67216   41.55567   46.35191   53.01875

为什么简单的Rcpp实现比zoo::rollmean快100倍？

最佳答案

在zoo中四处逛逛，看来rollmean.*方法都是在R中实现的。

而您是用C++实现的。打包的R代码可能还会再执行一些检查，例如pp，所以打败它并不奇怪吗？

关于r - 与简单的Rcpp实现相比，为什么zoo::rollmean慢？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/30090336/