那些遇到过的问题

自适应移动平均 - R中的最佳性能

jan*_*cki 16 r zoo mapply rollapply data.table

我正在寻找R中滚动/滑动窗口函数方面的一些性能提升.这是一个非常常见的任务,可用于任何有序的观测数据集.我想分享一些我的发现,也许有人能够提供反馈,使其更快.
重要的是我专注于案例align="right"和自适应滚动窗口,因此width是一个向量(与我们的观察向量相同的长度).如果我们有width标量,那么已经有非常好的函数zoo和TTR包非常难以击败(4年后:它比我预期的要容易),因为其中一些甚至使用Fortran(但仍然是用户定义的)使用下面提到的FUN可以更快wapply.
RcppRoll由于其出色的性能,包值得值得一提,但到目前为止还没有能够回答这个问题的功能.如果有人可以扩展它以回答这个问题,那将会很棒.

考虑一下我们有以下数据:

x = c(120,105,118,140,142,141,135,152,154,138,125,132,131,120)
plot(x, type="l")
Run Code Online (Sandbox Code Playgroud)

大块make_x的情节

我们希望在x带有可变滚动窗口的矢量上应用滚动函数width. 

set.seed(1)
width = sample(2:4,length(x),TRUE)
Run Code Online (Sandbox Code Playgroud)

在这种特殊情况下,我们将不得不滚动功能适应sample的c(2,3,4).
我们将应用mean功能,预期结果:

r = f(x, width, FUN = mean)
print(r)
##  [1]       NA       NA 114.3333 120.7500 141.0000 135.2500 139.5000
##  [8] 142.6667 147.0000 146.0000 131.5000 128.5000 131.5000 127.6667
plot(x, type="l")
lines(r, col="red")
Run Code Online (Sandbox Code Playgroud)

大块make_results的情节

任何指标都可用于产生自width变量作为自适应移动平均线的不同变体或任何其他函数.

寻找最佳表现.

jan*_*cki 23

2018年12月更新

最近在data.table中已经有效地实现了自适应滚动功能 - 更多信息在froll手册中.另外,已经鉴定了使用碱R的有效替代解决方案(fastama下文).不幸的是,Kevin Ushey的回答没有解决这个问题,因此不包括在基准测试中.由于毫无意义地比较微秒,基准的规模已经增加.

set.seed(108)
x = rnorm(1e6)
width = rep(seq(from = 100, to = 500, by = 5), length.out=length(x))
microbenchmark(
  zoo=rollapplyr(x, width = width, FUN=mean, fill=NA),
  mapply=base_mapply(x, width=width, FUN=mean, na.rm=T),
  wmapply=wmapply(x, width=width, FUN=mean, na.rm=T),
  ama=ama(x, width, na.rm=T),
  fastama=fastama(x, width),
  frollmean=frollmean(x, width, na.rm=T, adaptive=TRUE),
  frollmean_exact=frollmean(x, width, na.rm=T, adaptive=TRUE, algo="exact"),
  times=1L
)
#Unit: milliseconds
#            expr          min           lq         mean       median           uq          max neval
#             zoo 32371.938248 32371.938248 32371.938248 32371.938248 32371.938248 32371.938248     1
#          mapply 13351.726032 13351.726032 13351.726032 13351.726032 13351.726032 13351.726032     1
#         wmapply 15114.774972 15114.774972 15114.774972 15114.774972 15114.774972 15114.774972     1
#             ama  9780.239091  9780.239091  9780.239091  9780.239091  9780.239091  9780.239091     1
#         fastama   351.618042   351.618042   351.618042   351.618042   351.618042   351.618042     1
#       frollmean     7.708054     7.708054     7.708054     7.708054     7.708054     7.708054     1
# frollmean_exact   194.115012   194.115012   194.115012   194.115012   194.115012   194.115012     1
Run Code Online (Sandbox Code Playgroud) 
ama = function(x, n, na.rm=FALSE, fill=NA, nf.rm=FALSE) {
  # more or less the same as previous forloopply
  stopifnot((nx<-length(x))==length(n))
  if (nf.rm) x[!is.finite(x)] = NA_real_
  ans = rep(NA_real_, nx)
  for (i in seq_along(x)) {
    ans[i] = if (i >= n[i])
      mean(x[(i-n[i]+1):i], na.rm=na.rm)
    else as.double(fill)
  }
  ans
}
fastama = function(x, n, na.rm, fill=NA) {
  if (!missing(na.rm)) stop("fast adaptive moving average implemented in R does not handle NAs, input having NAs will result in incorrect answer so not even try to compare to it")
  # fast implementation of adaptive moving average in R, in case of NAs incorrect answer
  stopifnot((nx<-length(x))==length(n))
  cs = cumsum(x)
  ans = rep(NA_real_, nx)
  for (i in seq_along(cs)) {
    ans[i] = if (i == n[i])
      cs[i]/n[i]
    else if (i > n[i])
      (cs[i]-cs[i-n[i]])/n[i]
    else as.double(fill)
  }
  ans
}
Run Code Online (Sandbox Code Playgroud)

老答案:

我选择了4个可用的解决方案,这些解决方案不需要C++,很容易找到或谷歌.

# 1. rollapply
library(zoo)
?rollapplyr
# 2. mapply
base_mapply <- function(x, width, FUN, ...){
  FUN <- match.fun(FUN)
  f <- function(i, width, data){
    if(i < width) return(NA_real_)
    return(FUN(data[(i-(width-1)):i], ...))
  }
  mapply(FUN = f, 
         seq_along(x), width,
         MoreArgs = list(data = x))
}
# 3. wmapply - modified version of wapply found: https://rmazing.wordpress.com/2013/04/23/wapply-a-faster-but-less-functional-rollapply-for-vector-setups/
wmapply <- function(x, width, FUN = NULL, ...){
  FUN <- match.fun(FUN)
  SEQ1 <- 1:length(x)
  SEQ1[SEQ1 <  width] <- NA_integer_
  SEQ2 <- lapply(SEQ1, function(i) if(!is.na(i)) (i - (width[i]-1)):i)
  OUT <- lapply(SEQ2, function(i) if(!is.null(i)) FUN(x[i], ...) else NA_real_)
  return(base:::simplify2array(OUT, higher = TRUE))
}
# 4. forloopply - simple loop solution
forloopply <- function(x, width, FUN = NULL, ...){
  FUN <- match.fun(FUN)
  OUT <- numeric()
  for(i in 1:length(x)) {
    if(i < width[i]) next
    OUT[i] <- FUN(x[(i-(width[i]-1)):i], ...)
  }
  return(OUT)
}
Run Code Online (Sandbox Code Playgroud)

以下是prod功能的时间.mean功能可能已在内部优化rollapplyr.所有结果相同.

library(microbenchmark)
# 1a. length(x) = 1000, window = 5-20
x <- runif(1000,0.5,1.5)
width <- rep(seq(from = 5, to = 20, by = 5), length(x)/4)
microbenchmark(
  rollapplyr(data = x, width = width, FUN = prod, fill = NA),
  base_mapply(x = x, width = width, FUN = prod, na.rm=T),
  wmapply(x = x, width = width, FUN = prod, na.rm=T),
  forloopply(x = x, width = width, FUN = prod, na.rm=T),
  times=100L
)
Unit: milliseconds
                                                       expr       min        lq    median       uq       max neval
 rollapplyr(data = x, width = width, FUN = prod, fill = NA) 59.690217 60.694364 61.979876 68.55698 153.60445   100
   base_mapply(x = x, width = width, FUN = prod, na.rm = T) 14.372537 14.694266 14.953234 16.00777  99.82199   100
       wmapply(x = x, width = width, FUN = prod, na.rm = T)  9.384938  9.755893  9.872079 10.09932  84.82886   100
    forloopply(x = x, width = width, FUN = prod, na.rm = T) 14.730428 15.062188 15.305059 15.76560 342.44173   100

# 1b. length(x) = 1000, window = 50-200
x <- runif(1000,0.5,1.5)
width <- rep(seq(from = 50, to = 200, by = 50), length(x)/4)
microbenchmark(
  rollapplyr(data = x, width = width, FUN = prod, fill = NA),
  base_mapply(x = x, width = width, FUN = prod, na.rm=T),
  wmapply(x = x, width = width, FUN = prod, na.rm=T),
  forloopply(x = x, width = width, FUN = prod, na.rm=T),
  times=100L
)
Unit: milliseconds
                                                       expr      min       lq   median       uq      max neval
 rollapplyr(data = x, width = width, FUN = prod, fill = NA) 71.99894 74.19434 75.44112 86.44893 281.6237   100
   base_mapply(x = x, width = width, FUN = prod, na.rm = T) 15.67158 16.10320 16.39249 17.20346 103.6211   100
       wmapply(x = x, width = width, FUN = prod, na.rm = T) 10.88882 11.54721 11.75229 12.19790 106.1170   100
    forloopply(x = x, width = width, FUN = prod, na.rm = T) 15.70704 16.06983 16.40393 17.14210 108.5005   100

# 2a. length(x) = 10000, window = 5-20
x <- runif(10000,0.5,1.5)
width <- rep(seq(from = 5, to = 20, by = 5), length(x)/4)
microbenchmark(
  rollapplyr(data = x, width = width, FUN = prod, fill = NA),
  base_mapply(x = x, width = width, FUN = prod, na.rm=T),
  wmapply(x = x, width = width, FUN = prod, na.rm=T),
  forloopply(x = x, width = width, FUN = prod, na.rm=T),
  times=100L
)
Unit: milliseconds
                                                       expr       min       lq   median       uq       max neval
 rollapplyr(data = x, width = width, FUN = prod, fill = NA) 753.87882 781.8789 809.7680 872.8405 1116.7021   100
   base_mapply(x = x, width = width, FUN = prod, na.rm = T) 148.54919 159.9986 231.5387 239.9183  339.7270   100
       wmapply(x = x, width = width, FUN = prod, na.rm = T)  98.42682 105.2641 117.4923 183.4472  245.4577   100
    forloopply(x = x, width = width, FUN = prod, na.rm = T) 533.95641 602.0652 646.7420 672.7483  922.3317   100

# 2b. length(x) = 10000, window = 50-200
x <- runif(10000,0.5,1.5)
width <- rep(seq(from = 50, to = 200, by = 50), length(x)/4)
microbenchmark(
  rollapplyr(data = x, width = width, FUN = prod, fill = NA),
  base_mapply(x = x, width = width, FUN = prod, na.rm=T),
  wmapply(x = x, width = width, FUN = prod, na.rm=T),
  forloopply(x = x, width = width, FUN = prod, na.rm=T),
  times=100L
)
Unit: milliseconds
                                                       expr      min       lq    median        uq       max neval
 rollapplyr(data = x, width = width, FUN = prod, fill = NA) 912.5829 946.2971 1024.7245 1071.5599 1431.5289   100
   base_mapply(x = x, width = width, FUN = prod, na.rm = T) 171.3189 180.6014  260.8817  269.5672  344.4500   100
       wmapply(x = x, width = width, FUN = prod, na.rm = T) 123.1964 131.1663  204.6064  221.1004  484.3636   100
    forloopply(x = x, width = width, FUN = prod, na.rm = T) 561.2993 696.5583  800.9197  959.6298 1273.5350   100
Run Code Online (Sandbox Code Playgroud)

  • 不同意,我认为这没有错.SO是知识库.我从一个人那里学到了一些东西. (20认同)
  • @Dirk,我不想只为一个共享创建一个博客.SO也是一个非常好的讨论板.关于滚动窗口FUN的问题在SO上很常见,但我没有找到任何好的基准.我仍然希望在没有去Cpp的情况下获得该领域的一些性能提升,所以我期待着更好地回答我的问题. (5认同)
  • @DirkEddelbuettel,SO确实鼓励用户提出并回答他们自己的问题,只要这是通过正确的礼仪完成的.见[here](http://meta.stackexchange.com/questions/17845/etiquette-for-answering-your-own-question). (5认同)
  • 让我换一种说法,如果其他人提出这个问题,并且MusX已经回答了问题,那么他的回答将被视为一种积极的贡献.而且这也不是一个无用的问题,所以我真的看不出这是多么糟糕. (5认同)
  • 对不起,没有耐心的meta.se肚脐凝视这里.曾经是一个很好的资源.现在它被那些不能/不会自己做任何研究,重复问题和其他滥用平台的人所淹没.将它用作个人CMS仍然不适合我. (4认同)
  • 根据Kevin RcppRoll的信息处理`width`只作为标量,所以默认情况下它不会解决问题.我不知道Cpp,但我相信我可能会做一些调整来处理这个问题,就像我从r-bloggers那里做的那样 - 它也只接受标量`width`.使用R-cpp接口可以胜过大多数任务.没有'width`作为标量,很难比较时间,即使这样,最好的非cpp解决方案与cpp解决方案时间比较也是好的. (4认同)
  • @DirkEddelbuettel,同意很多懒惰,但你必须认识到这里根本不是问题.显然,MusX花了很多时间将他的问题和答案放在一起.我已经阅读了很多标准的R,并且他描述的问题通常没有记录.还有更糟糕的行为要贬低/大喊大叫. (2认同)
  • 这个问题有一个很好的(基于C++的)解决方案---由Kevin作为CRAN包提供.如果没有提供任何细节,OP就会"无法测试",并采用大量非编译解决方案来浪费每个人的时间.这是一个贡献？ (2认同)

Kev*_*hey 22

作为参考,你一定要看看RcppRoll你是否只有一个窗口长度来"滚动":

library(RcppRoll) ## install.packages("RcppRoll")
library(microbenchmark)
x <- runif(1E5)
all.equal( rollapplyr(x, 10, FUN=prod), roll_prod(x, 10) )
microbenchmark( times=5,
  rollapplyr(x, 10, FUN=prod),
  roll_prod(x, 10)
)
Run Code Online (Sandbox Code Playgroud)

给我

> library(RcppRoll)
> library(microbenchmark)
> x <- runif(1E5)
> all.equal( rollapplyr(x, 10, FUN=prod), roll_prod(x, 10) )
[1] TRUE
> microbenchmark( times=5,
+   zoo=rollapplyr(x, 10, FUN=prod),
+   RcppRoll=roll_prod(x, 10)
+ )
Unit: milliseconds
     expr        min         lq     median         uq         max neval
      zoo 924.894069 968.467299 997.134932 1029.10883 1079.613569     5
 RcppRoll   1.509155   1.553062   1.760739    1.90061    1.944999     5
Run Code Online (Sandbox Code Playgroud)

它的速度要快一些;)并且软件包足够灵活,用户可以定义和使用自己的滚动功能(使用C++).我可以在将来扩展包以允许多个窗口宽度,但我相信要做到正确将是棘手的.

如果你想定义prod自己,你可以这样做 - RcppRoll允许你定义你自己的C++函数来传递并生成一个'滚动'函数,如果你愿意的话.rollit提供了一个更好的界面,而rollit_raw只是让你自己编写一个C++函数,有点像你可能会做的Rcpp::cppFunction.理念是,您应该只需要表达您希望在特定窗口上执行的计算,并且RcppRoll可以负责迭代某些大小的窗口.

library(RcppRoll)
library(microbenchmark)
x <- runif(1E5)
my_rolling_prod <- rollit(combine="*")
my_rolling_prod2 <- rollit_raw("
double output = 1;
for (int i=0; i < n; ++i) {
  output *= X(i);
}
return output;
")
all.equal( roll_prod(x, 10), my_rolling_prod(x, 10) )
all.equal( roll_prod(x, 10), my_rolling_prod2(x, 10) )
microbenchmark( times=5,
  rollapplyr(x, 10, FUN=prod),
  roll_prod(x, 10),
  my_rolling_prod(x, 10),
  my_rolling_prod2(x, 10)
)
Run Code Online (Sandbox Code Playgroud)

给我

> library(RcppRoll)
> library(microbenchmark)
> # 1a. length(x) = 1000, window = 5-20
> x <- runif(1E5)
> my_rolling_prod <- rollit(combine="*")
C++ source file written to /var/folders/m7/_xnnz_b53kjgggkb1drc1f8c0000gn/T//RtmpcFMJEV/file80263aa7cca2.cpp .
Compiling...
Done!
> my_rolling_prod2 <- rollit_raw("
+ double output = 1;
+ for (int i=0; i < n; ++i) {
+   output *= X(i);
+ }
+ return output;
+ ")
C++ source file written to /var/folders/m7/_xnnz_b53kjgggkb1drc1f8c0000gn/T//RtmpcFMJEV/file802673777da2.cpp .
Compiling...
Done!
> all.equal( roll_prod(x, 10), my_rolling_prod(x, 10) )
[1] TRUE
> all.equal( roll_prod(x, 10), my_rolling_prod2(x, 10) )
[1] TRUE
> microbenchmark(
+   rollapplyr(x, 10, FUN=prod),
+   roll_prod(x, 10),
+   my_rolling_prod(x, 10),
+   my_rolling_prod2(x, 10)
+ )

> microbenchmark( times=5,
+   rollapplyr(x, 10, FUN=prod),
+   roll_prod(x, 10),
+   my_rolling_prod(x, 10),
+   my_rolling_prod2(x, 10)
+ )
Unit: microseconds
                          expr        min          lq      median          uq         max neval
 rollapplyr(x, 10, FUN = prod) 979710.368 1115931.323 1117375.922 1120085.250 1149117.854     5
              roll_prod(x, 10)   1504.377    1635.749    1638.943    1815.344    2053.997     5
        my_rolling_prod(x, 10)   1507.687    1572.046    1648.031    2103.355    7192.493     5
       my_rolling_prod2(x, 10)    774.381     786.750     884.951    1052.508    1434.660     5
Run Code Online (Sandbox Code Playgroud)

所以,实际上,只要你能够通过rollit接口或通过C++函数表达你想在特定窗口中执行的计算rollit_raw(其界面有点僵硬,但仍然有用),你就可以了. .

  • 我更新了我的答案,举例说明了如何使用`RcppRoll`的用户定义函数; 特别是,在这里表达"prod"的两种方式. (4认同)

Lau*_*riK 5

不知怎的,人们已经错过了runmed()基地R(统计数据包)中的超快速.就我理解原始问题而言,它不具有自适应性,但对于滚动中位数,它是快速的!这里roll_median()与RcppRoll 比较.

> microbenchmark(
+   runmed(x = x, k = 3),
+   roll_median(x, 3),
+   times=1000L
+ )
Unit: microseconds
                 expr     min      lq      mean  median      uq     max neval
 runmed(x = x, k = 3)  41.053  44.854  47.60973  46.755  49.795 117.838  1000
    roll_median(x, 3) 101.872 105.293 108.72840 107.574 111.375 178.657  1000
Run Code Online (Sandbox Code Playgroud)

归档时间：

查看次数：

5555 次

最近记录：

6 年，8 月 前
R data.table与rollapply 5
更多相关链接
相关归档
Emacs ESS模式 - 评论区域的标签 49
为什么使用as.factor()而不仅仅是因子() 45
如何将`==`行为扩展到包含NA的向量？ 30
您是使用attach()还是通过名称或切片调用变量？ 26
如何在ggplot2图例中使用下标[R] 25
使用R将PDF文件转换为文本文件以进行文本挖掘 17
重新编码/重新定位具有不同级别的数据框架因子 10
R:将数据集拆分为四分位数/十分位数.什么是正确的方法？ 10
如何使用R的data.table包来对键值的否定进行子集化？ 5
优化R中的循环 2
难疑归档
使用Git将特定文件重置或还原到特定版本？ 4255
我应该在MySQL中使用日期时间或时间戳数据类型吗？ 2598
如何从YouTube API获取YouTube视频缩略图？ 2291
如何在Python中获取列表中的元素数量？ 1846
在Vim中复制整行 1555
如何用"喜欢"查询MongoDB？ 1331
JavaScript中是否有"null coalescing"运算符？ 1305
如何查看仅一个用户提交的git日志？ 1178
基于容器宽度的字体缩放 1083
如何在JavaScript中将浮点数转换为整数？ 1043