Адаптивная скользящая средняя – верхняя производительность в R

Я ищу некоторые улучшения производительности с точки зрения функций качения / сдвига windows в R. Это довольно распространенная задача, которая может быть использована в любом упорядоченном наборе данных наблюдений. Я хотел бы поделиться некоторыми своими выводами, возможно, кто-то сможет предоставить обратную связь, чтобы сделать это еще быстрее.
Важно отметить, что я фокусируюсь на align="right" case align="right" и width виде вектора (такой же длины, как и наш вектор наблюдения). В случае, если у нас есть width как скаляр, есть уже очень хорошо развитые функции в пакетах zoo и TTR которые очень сложно обыграть, поскольку некоторые из них даже используют Fortran (но все же пользовательские FUN могут быть быстрее, используя упомянутый ниже wapply ) ,
Пакет RcppRoll заслуживает упоминания из-за его высокой производительности, но до сих пор нет функции, которая отвечает на этот вопрос. Было бы здорово, если бы кто-то мог продлить его, чтобы ответить на вопрос.

Рассмотрим, что у нас есть следующие данные:

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

участок куска make_x

И мы хотим применить функцию качения по вектору x с переменной width прокатки.

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

В этом конкретном случае у нас была бы функция качения, адаптивная к 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") 

участок куски make_results

Любой индикатор может использоваться для создания аргумента width виде различных вариантов адаптивных скользящих средних или любой другой функции.

Ищите лучшую производительность.

    Для справки, вы должны обязательно проверить RcppRoll если у вас есть только одна длина windows, чтобы «перевернуть»:

     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 ++). Я могу расширить пакет в будущем, чтобы позволить несколько оконных ширин, но я уверен, что будет сложно получить право.

    Если вы хотите определить prod самостоятельно, вы можете это сделать - RcppRoll позволяет вам определять свои собственные функции C ++ для передачи и генерации функции «прокатки», если хотите. rollit дает несколько более приятный интерфейс, в то время как rollit_raw просто позволяет вам писать C ++-функцию самостоятельно, как вы можете сделать с Rcpp::cppFunction . Философия заключается в том, что вам нужно будет только вычислить вычисления, которые вы хотите выполнить в определенном окне, и RcppRoll может позаботиться об итерации над windowsми некоторого размера.

     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 или с помощью функции C ++, переданной через rollit_raw (интерфейс которой немного жесткий, но все еще функциональный), вы находитесь в хорошая фигура.

    Я выбрал 4 доступных решения, которые не нужно делать с C ++, довольно легко найти или 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) } 

    Ниже приведены тайминги для функции 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 

    Как-то люди пропустили сверхбыструю скорость 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 
    Давайте будем гением компьютера.