zoo::rollmean
是一个有用的函数,它返回时间序列的滚动均值;对于长度为 n
、窗口大小为 k
的向量 x
,它返回向量 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倍?
RcppRoll
软件包提供了zoo::roll
函数更快的实现方式。 - ExperimenteR