Función de ventana móvil adaptativa: máximo rendimiento en R

data.table zoo

Question

Estoy buscando algunas mejoras de rendimiento en términos de funciones de ventana rodante / deslizante en R. Es una tarea bastante común que se puede utilizar en cualquier conjunto de datos de observaciones ordenadas. Me gustaría compartir algunos de mis hallazgos, tal vez alguien podría proporcionar comentarios para hacerlo aún más rápido.
Una nota importante es que me concentro en el caso de align="right" y el width ventana móvil como vector (la misma longitud que nuestro vector de observación). En caso de que tengamos el width como escalar, ya existen funciones muy bien desarrolladas en zoo y paquetes TTR que serían muy difíciles de superar ya que algunos de ellos incluso usan Fortran (pero aún así las FUNs definidas por el usuario pueden ser más rápidas usando los siguientes wapply ) .
RcppRoll pena mencionar el paquete RcppRoll debido a su gran rendimiento, pero hasta ahora no hay ninguna función que responda a esa pregunta. Sería genial si alguien pudiera extenderlo para responder la pregunta.

Considere que tenemos una siguiente información:

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

Y queremos aplicar la función de balanceo sobre x vector con width variable de ventana móvil.

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

En este caso particular, tendríamos una función de desplazamiento adaptable a la sample de c(2,3,4) .
Aplicaremos la función mean , resultados esperados:

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")

Se puede emplear cualquier indicador para producir un argumento de width como diferentes variantes de promedios móviles adaptativos, o cualquier otra función.

Buscando un rendimiento superior.

Answer 1

Como referencia, definitivamente debe verificar RcppRoll si solo tiene una longitud de ventana única para ''pasar'':

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

me da

> 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

Es un poco más rápido;) y el paquete es lo suficientemente flexible para que los usuarios definan y usen sus propias funciones de desplazamiento (con C ++). Puedo extender el paquete en el futuro para permitir múltiples anchos de ventana, pero estoy seguro de que será complicado hacerlo bien.

Si quiere definir el prod usted mismo, puede hacerlo: RcppRoll permite definir sus propias funciones de C ++ para pasar y generar una función ''rolling'' si lo desea. rollit ofrece una interfaz algo más agradable, mientras que rollit_raw simplemente te permite escribir una función C ++, algo parecido a lo que harías con Rcpp::cppFunction . Siendo la filosofía, solo debe expresar el cálculo que desea realizar en una ventana en particular, y RcppRoll puede ocuparse de iterar en ventanas de algún tamaño.

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

me da

> 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

Entonces, en la medida en que sea capaz de expresar el cálculo que desea realizar en una ventana particular a través de la interfaz rollit o con una función C ++ pasada a través de rollit_raw (cuya interfaz es un poco rígida, pero sigue siendo funcional), usted está en buena forma.

Answer 2

Elegí 4 soluciones disponibles que no necesitan hacer para C ++, bastante fáciles de encontrar o google.

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

A continuación se muestran los tiempos para la función prod . mean función mean puede estar ya optimizada dentro de rollapplyr . Todos los resultados son iguales.

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

Answer 3

De alguna manera la gente se ha perdido el runmed() ultrarrápido runmed() en la base R (paquete de estadísticas). No es adaptativo, por lo que entiendo la pregunta original, pero para una mediana móvil, ¡es RÁPIDO! Comparando aquí con roll_median() desde 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