Расчет скользящей средней
Я пытаюсь использовать R для вычисления скользящего среднего по ряду значений в матрице. Обычный поиск по списку рассылки R не очень помог. Кажется, в R нет встроенной функции, которая позволяла бы мне вычислять скользящие средние. Какие-нибудь пакеты предоставляют? Или мне нужно написать свое?
18 ответов
Или вы можете просто рассчитать его с помощью фильтра, вот функция, которую я использую:
ma <- function(x,n=5){filter(x,rep(1/n,n), sides=2)}
С помощью cumsum должно быть достаточно и эффективно. Предполагая, что у вас есть вектор x, и вы хотите получить текущую сумму из n чисел
cx <- c(0,cumsum(x))
rsum <- (cx[(n+1):length(cx)] - cx[1:(length(cx) - n)]) / n
Как отмечено в комментариях @mzuther, это предполагает, что в данных нет NA. чтобы справиться с этим, потребуется разделить каждое окно на количество значений, отличных от NA. Вот один из способов сделать это, включив комментарий @Ricardo Cruz:
cx <- c(0, cumsum(ifelse(is.na(x), 0, x)))
cn <- c(0, cumsum(ifelse(is.na(x), 0, 1)))
rx <- cx[(n+1):length(cx)] - cx[1:(length(cx) - n)]
rn <- cn[(n+1):length(cx)] - cx[1:(length(cx) - n)]
rsum <- rx / rn
Это все еще имеет проблему, что, если все значения в окне являются NA, тогда будет ошибка деления на ноль.
В data.table 1.12.0 new frollmean добавлена функция для быстрого и точного вычисления среднего значения NA, NaN а также +Inf, -Inf ценности.
Так как в этом вопросе нет воспроизводимого примера, здесь не так много вопросов.
Вы можете найти больше информации о ?frollmean в ручном режиме, также доступны онлайн на ?frollmean,
Примеры из руководства ниже:
library(data.table)
d = as.data.table(list(1:6/2, 3:8/4))
# rollmean of single vector and single window
frollmean(d[, V1], 3)
# multiple columns at once
frollmean(d, 3)
# multiple windows at once
frollmean(d[, .(V1)], c(3, 4))
# multiple columns and multiple windows at once
frollmean(d, c(3, 4))
## three above are embarrassingly parallel using openmp
caTools Пакет имеет очень быстрое скользящее среднее / мин / макс / сд и несколько других функций. Я работал только с runmean а также runsd и они являются самыми быстрыми из любых других пакетов, упомянутых на сегодняшний день.
Вот пример кода, показывающий, как вычислить центрированную скользящую среднюю и скользящую среднюю скользящую с помощьюrollmeanфункция из пакета зоопарка.
library(tidyverse)
library(zoo)
some_data = tibble(day = 1:10)
# cma = centered moving average
# tma = trailing moving average
some_data = some_data %>%
mutate(cma = rollmean(day, k = 3, fill = NA)) %>%
mutate(tma = rollmean(day, k = 3, fill = NA, align = "right"))
some_data
#> # A tibble: 10 x 3
#> day cma tma
#> <int> <dbl> <dbl>
#> 1 1 NA NA
#> 2 2 2 NA
#> 3 3 3 2
#> 4 4 4 3
#> 5 5 5 4
#> 6 6 6 5
#> 7 7 7 6
#> 8 8 8 7
#> 9 9 9 8
#> 10 10 NA 9
Вы могли бы использовать RcppRoll для очень быстрых скользящих средних, написанных на C++. Просто позвони roll_mean функция. Документы можно найти здесь.
В противном случае этот (более медленный) цикл for должен сработать:
ma <- function(arr, n=15){
res = arr
for(i in n:length(arr)){
res[i] = mean(arr[(i-n):i])
}
res
}
По факту RcppRoll очень хорошо.
Код, размещенный cantdutchthis, должен быть исправлен в четвертой строке, чтобы окно было исправлено:
ma <- function(arr, n=15){
res = arr
for(i in n:length(arr)){
res[i] = mean(arr[(i-n+1):i])
}
res
}
Другой способ, который обрабатывает пропуски, дан здесь.
Третий способ - улучшить cantdutchthis код для вычисления частичных средних значений или нет - следующим образом:
ma <- function(x, n=2,parcial=TRUE){
res = x #set the first values
if (parcial==TRUE){
for(i in 1:length(x)){
t<-max(i-n+1,1)
res[i] = mean(x[t:i])
}
res
}else{
for(i in 1:length(x)){
t<-max(i-n+1,1)
res[i] = mean(x[t:i])
}
res[-c(seq(1,n-1,1))] #remove the n-1 first,i.e., res[c(-3,-4,...)]
}
}
Вы можете рассчитать скользящее среднее вектора
x с шириной окна
k по:
apply(embed(x, k), 1, mean)
Для того, чтобы дополнить ответ cantdutchthis и Rodrigo Remedio;
moving_fun <- function(x, w, FUN, ...) {
# x: a double vector
# w: the length of the window, i.e., the section of the vector selected to apply FUN
# FUN: a function that takes a vector and return a summarize value, e.g., mean, sum, etc.
# Given a double type vector apply a FUN over a moving window from left to the right,
# when a window boundary is not a legal section, i.e. lower_bound and i (upper bound)
# are not contained in the length of the vector, return a NA_real_
if (w < 1) {
stop("The length of the window 'w' must be greater than 0")
}
output <- x
for (i in 1:length(x)) {
# plus 1 because the index is inclusive with the upper_bound 'i'
lower_bound <- i - w + 1
if (lower_bound < 1) {
output[i] <- NA_real_
} else {
output[i] <- FUN(x[lower_bound:i, ...])
}
}
output
}
# example
v <- seq(1:10)
# compute a MA(2)
moving_fun(v, 2, mean)
# compute moving sum of two periods
moving_fun(v, 2, sum)
Для этого можно использовать пакет слайдера. У него есть интерфейс, который был специально разработан, чтобы чувствовать себя похожим на мурлыканье. Он принимает любую произвольную функцию и может возвращать любой тип вывода. Кадры данных даже повторяются по строкам. Сайт pkgdown находится здесь.
library(slider)
x <- 1:3
# Mean of the current value + 1 value before it
# returned as a double vector
slide_dbl(x, ~mean(.x, na.rm = TRUE), .before = 1)
#> [1] 1.0 1.5 2.5
df <- data.frame(x = x, y = x)
# Slide row wise over data frames
slide(df, ~.x, .before = 1)
#> [[1]]
#> x y
#> 1 1 1
#>
#> [[2]]
#> x y
#> 1 1 1
#> 2 2 2
#>
#> [[3]]
#> x y
#> 1 2 2
#> 2 3 3
Накладные расходы как на слайдер, так и на data.table frollapply() должно быть довольно низким (намного быстрее, чем зоопарк). frollapply()выглядит немного быстрее для этого простого примера здесь, но обратите внимание, что он принимает только числовой ввод, а вывод должен быть скалярным числовым значением. функции слайдера полностью универсальны, и вы можете возвращать любой тип данных.
library(slider)
library(zoo)
library(data.table)
x <- 1:50000 + 0L
bench::mark(
slider = slide_int(x, function(x) 1L, .before = 5, .complete = TRUE),
zoo = rollapplyr(x, FUN = function(x) 1L, width = 6, fill = NA),
datatable = frollapply(x, n = 6, FUN = function(x) 1L),
iterations = 200
)
#> # A tibble: 3 x 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 slider 19.82ms 26.4ms 38.4 829.8KB 19.0
#> 2 zoo 177.92ms 211.1ms 4.71 17.9MB 24.8
#> 3 datatable 7.78ms 10.9ms 87.9 807.1KB 38.7
Для людей, просто желающих вычислить это самостоятельно, это не более чем:
# x = vector with numeric data
# w = window length
y <- numeric(length = length(x))
for (i in seq_len(length(x))) {
ind <- c((i - floor(w / 2)):(i + floor(w / 2)))
ind <- ind[ind %in% seq_len(length(x))]
y[i] <- mean(x[ind])
}
y
Но становится забавно делать его независимым от mean(), так что вы можете вычислить любую "движущуюся" функцию!
# our working horse:
moving_fn <- function(x, w, fun, ...) {
# x = vector with numeric data
# w = window length
# fun = function to apply
# ... = parameters passed on to 'fun'
y <- numeric(length(x))
for (i in seq_len(length(x))) {
ind <- c((i - floor(w / 2)):(i + floor(w / 2)))
ind <- ind[ind %in% seq_len(length(x))]
y[i] <- fun(x[ind], ...)
}
y
}
# and now any variation you can think of!
moving_average <- function(x, w = 5, na.rm = FALSE) {
moving_fn(x = x, w = w, fun = mean, na.rm = na.rm)
}
moving_sum <- function(x, w = 5, na.rm = FALSE) {
moving_fn(x = x, w = w, fun = sum, na.rm = na.rm)
}
moving_maximum <- function(x, w = 5, na.rm = FALSE) {
moving_fn(x = x, w = w, fun = max, na.rm = na.rm)
}
moving_median <- function(x, w = 5, na.rm = FALSE) {
moving_fn(x = x, w = w, fun = median, na.rm = na.rm)
}
moving_Q1 <- function(x, w = 5, na.rm = FALSE) {
moving_fn(x = x, w = w, fun = quantile, na.rm = na.rm, 0.25)
}
moving_Q3 <- function(x, w = 5, na.rm = FALSE) {
moving_fn(x = x, w = w, fun = quantile, na.rm = na.rm, 0.75)
}
Можно использовать runner пакет для движущихся функций. В этом случаеmean_runфункция. Проблема сcummean в том, что он не справляется NA ценности, но mean_run делает. runner пакет также поддерживает нерегулярные временные ряды, и окна могут зависеть от даты:
library(runner)
set.seed(11)
x1 <- rnorm(15)
x2 <- sample(c(rep(NA,5), rnorm(15)), 15, replace = TRUE)
date <- Sys.Date() + cumsum(sample(1:3, 15, replace = TRUE))
mean_run(x1)
#> [1] -0.5910311 -0.2822184 -0.6936633 -0.8609108 -0.4530308 -0.5332176
#> [7] -0.2679571 -0.1563477 -0.1440561 -0.2300625 -0.2844599 -0.2897842
#> [13] -0.3858234 -0.3765192 -0.4280809
mean_run(x2, na_rm = TRUE)
#> [1] -0.18760011 -0.09022066 -0.06543317 0.03906450 -0.12188853 -0.13873536
#> [7] -0.13873536 -0.14571604 -0.12596067 -0.11116961 -0.09881996 -0.08871569
#> [13] -0.05194292 -0.04699909 -0.05704202
mean_run(x2, na_rm = FALSE )
#> [1] -0.18760011 -0.09022066 -0.06543317 0.03906450 -0.12188853 -0.13873536
#> [7] NA NA NA NA NA NA
#> [13] NA NA NA
mean_run(x2, na_rm = TRUE, k = 4)
#> [1] -0.18760011 -0.09022066 -0.06543317 0.03906450 -0.10546063 -0.16299272
#> [7] -0.21203756 -0.39209010 -0.13274756 -0.05603811 -0.03894684 0.01103493
#> [13] 0.09609256 0.09738460 0.04740283
mean_run(x2, na_rm = TRUE, k = 4, idx = date)
#> [1] -0.187600111 -0.090220655 -0.004349696 0.168349653 -0.206571573 -0.494335093
#> [7] -0.222969541 -0.187600111 -0.087636571 0.009742884 0.009742884 0.012326968
#> [13] 0.182442234 0.125737145 0.059094786
Также можно указать другие параметры, например lag, и катиться только atконкретные индексы. Подробнее в документации по пакетам и функциям.
Вот простая функция с filter демонстрация одного из способов позаботиться о начале и завершении NA с заполнением и вычисление средневзвешенного значения (поддерживается filter) с использованием нестандартных весов:
wma <- function(x) {
wts <- c(seq(0.5, 4, 0.5), seq(3.5, 0.5, -0.5))
nside <- (length(wts)-1)/2
# pad x with begin and end values for filter to avoid NAs
xp <- c(rep(first(x), nside), x, rep(last(x), nside))
z <- stats::filter(xp, wts/sum(wts), sides = 2) %>% as.vector
z[(nside+1):(nside+length(x))]
}
Хотя и немного медленно, но вы также можете использовать zoo::rollapply для выполнения вычислений на матрицах.
reqd_ma <- rollapply(x, FUN = mean, width = n)
где x - набор данных, FUN = mean - функция; Вы также можете изменить его на min, max, sd и т. д., а width - это скользящее окно.
Еще одна полезная функция, если вы хотите, чтобы два конца ряда были не NA, а представляли собой рекурсивно рассчитанные скользящие средние:
smoothing = function(x, k=1) {
sapply(seq_along(x), function(i) {
i.min = max(i-k, 1)
i.max = min(i+k, length(x))
mean(x[i.min:i.max], na.rm=TRUE)
})
}
Пример:
x = 1:10/2
[1] 0,5 1,0 1,5 2,0 2,5 3,0 3,5 4,0 4,5 5,0
smoothing(x, 2)
[1] 1,00 1,25 1,50 2,00 2,50 3,00 3,50 4,00 4,25 4,50
vector_avg <- function(x){
sum_x = 0
for(i in 1:length(x)){
if(!is.na(x[i]))
sum_x = sum_x + x[i]
}
return(sum_x/length(x))
}
Я использую агрегат вместе с вектором, созданным rep(). Преимущество этого заключается в использовании cbind() для одновременной агрегации более 1 столбца в вашем фрейме данных. Ниже приведен пример скользящего среднего значения 60 для вектора (v) длины 1000:
v=1:1000*0.002+rnorm(1000)
mrng=rep(1:round(length(v)/60+0.5), length.out=length(v), each=60)
aggregate(v~mrng, FUN=mean, na.rm=T)
Обратите внимание, что первый аргумент в rep состоит в том, чтобы просто получить достаточное количество уникальных значений для скользящего диапазона, исходя из длины вектора и величины, подлежащей усреднению; второй аргумент сохраняет длину равной длине вектора, а последний повторяет значения первого аргумента столько раз, сколько период усреднения.
В совокупности вы можете использовать несколько функций (медиана, максимум, минимум) - например, показано среднее значение. Опять же, можно было бы использовать формулу с cbind, чтобы сделать это для более чем одного (или всех) столбцов в фрейме данных.