with - rx axis




적응 롤링 윈도우 기능-R의 최고 성능 (2)

나는 R에서 롤링 / 슬라이딩 윈도우 함수의 관점에서 약간의 성능 향상을 기대하고있다. 이것은 어떤 순서화 된 관측 데이터 세트에서도 사용될 수있는 매우 일반적인 작업이다. 내 연구 결과 중 일부를 공유하고 싶습니다. 아마 누군가가 피드백을 제공하여 더 빨리 만들 수있을 것입니다.
중요한 것은 case align="right" 와 롤링 윈도우의 width 을 벡터 (우리의 관찰 벡터와 같은 길이)로 초점을 align="right" 입니다. 스칼라로서 width 가있는 경우에는 이미 zooTTR 패키지에서 아주 잘 개발 된 함수가 있습니다. 그 중 일부는 Fortran을 사용하기 때문에 매우 어려울 것입니다 (그러나 여전히 사용자 정의 FUN은 아래 언급 된 wapply 사용하여 더 빠를 수 있습니다) .
RcppRoll 패키지는 뛰어난 성능으로 인해 언급할만한 가치가 있지만, 지금까지는 그 질문에 대한 답변을 얻을 수있는 기능이 없습니다. 누군가가 질문에 대답 할 수 있다면 좋을 것입니다.

다음 데이터가 있다고 가정 해보십시오.

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

그리고 우리는 가변 롤링 윈도우 width 갖는 x 벡터에 롤링 함수를 적용하고자합니다.

set.seed(1)
width = sample(2:4,length(x),TRUE)

이 특별한 경우에 우리는 c(2,3,4) sample 에 적응하는 롤링 함수를 가질 것이다.
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")

적응 형 이동 평균의 다른 변형 또는 다른 함수로 width 인수를 생성하기 위해 모든 표시기를 사용할 수 있습니다.

최고 실적을 찾고 있습니다.


참고로, '회전'할 창 길이가 하나 RcppRoll 경우 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)
)

나에게 준다.

> 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

조금 더 빠릅니다.) 패키지는 사용자가 롤링 기능을 정의하고 사용할 수있을 정도로 유연합니다 (C ++ 사용). 장래에 여러 창 너비를 허용하도록 패키지를 확장 할 수는 있지만 올바른 결과를 얻으려면 까다로울 것입니다.

RcppRoll 을 정의하고 싶다면 그렇게 할 수 있습니다 - RcppRoll 은 자신 만의 C ++ 함수를 정의하여 전달하고 '롤링'기능을 생성 할 수있게합니다. rollit 은 다소 멋진 인터페이스를 제공하지만 rollit_rawRcpp::cppFunction 사용하는 것처럼 C ++ 함수를 직접 작성할 수 있습니다. 이 철학은 특정 창에서 수행하려는 계산 만 표현하면되며 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)
)

나에게 준다.

> 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

rollit 인터페이스 나 rollit_raw 통해 전달 된 C ++ 함수를 통해 특정 창에서 수행하려는 계산을 표현할 수있는 한 (인터페이스는 약간 딱딱하지만 여전히 기능적입니다), 실제로 좋은 모양.


C + +로 할 필요가없는 4 가지 솔루션을 선택 했으므로 쉽게 찾을 수 있습니다.

# 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)
}

아래는 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