Obliczanie średniej kroczącej

185

Próbuję użyć R do obliczenia średniej ruchomej dla szeregu wartości w macierzy. Jednak normalne wyszukiwanie listy mailingowej R nie było zbyt pomocne. Wydaje się, że nie ma wbudowanej funkcji w R, która pozwoli mi obliczyć średnie ruchome. Czy jakieś pakiety to zapewniają? Czy też muszę pisać własne?

Jared
źródło

Odpowiedzi:

140
  • Rolling Means / Maximums / Medians in the zoo package (rollmean)
  • MovingAverages in TTR
  • ma w prognozie
f3lix
źródło
1
Jaka jest średnia ruchoma w R niezawierająca przyszłych wartości danego znacznika czasu? Sprawdziłem forecast::mai zawiera całą okolicę, nie w porządku.
hhh
213

Lub możesz po prostu obliczyć to za pomocą filtra, oto funkcja, której używam:

ma <- function(x, n = 5){filter(x, rep(1 / n, n), sides = 2)}

Jeśli używasz dplyr, pamiętaj, aby określić stats::filterw powyższej funkcji.

Matti Pastell
źródło
49
Powinienem zaznaczyć, że „side = 2” może być ważną opcją w wielu przypadkach użycia, których nie chcą przeoczyć. Jeśli chcesz tylko końcowych informacji w średniej ruchomej, powinieneś użyć side = 1.
evanrsparks
35
Kilka lat później, ale dplyr ma teraz funkcję filtrowania, jeśli masz ten pakiet załadowany użyjstats::filter
blmoore
sides = 2jest równoważne align = "center" dla zoo :: rollmean lub RcppRoll :: roll_mean. sides = 1jest równoważne wyrównaniu z „prawym”. Nie widzę sposobu, aby wykonać wyrównanie do lewej strony lub obliczyć za pomocą danych „częściowych” (2 lub więcej wartości)?
Matt L.
29

Używanie cumsumpowinno być wystarczające i wydajne. Zakładając, że masz wektor x i chcesz bieżącą sumę n liczb

cx <- c(0,cumsum(x))
rsum <- (cx[(n+1):length(cx)] - cx[1:(length(cx) - n)]) / n

Jak wskazano w komentarzach @mzuther, zakłada się, że w danych nie ma NA. radzenie sobie z nimi wymagałoby podzielenia każdego okna przez liczbę wartości innych niż NA. Oto jeden ze sposobów, aby to zrobić, uwzględniając komentarz @ 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)] - cn[1:(length(cx) - n)]
rsum <- rx / rn

Nadal występuje problem polegający na tym, że jeśli wszystkie wartości w oknie są NA, wówczas wystąpi błąd dzielenia przez zero.

iglicznia
źródło
8
Jednym minusem tego rozwiązania jest to, że nie cumsum(c(1:3,NA,1:3))
radzi
W ten sposób możesz łatwo poradzić sobie z NA cx <- c(0, cumsum(ifelse(is.na(x), 0, x))).
Ricardo Cruz
@ Ricardo Cruz: może być lepiej usunąć NA i odpowiednio dostosować długość wektora. Pomyśl o wektorze z dużą ilością NA - zera przyciągną średnią do zera, a usunięcie NA pozostawi średnią tak, jak jest. Wszystko zależy oczywiście od twoich danych i pytania, na które chcesz odpowiedzieć. :)
mzuther
@mzuther, zaktualizowałem odpowiedź po komentarzach. Dzięki za wkład. Myślę, że właściwym sposobem radzenia sobie z brakującymi danymi nie jest rozszerzenie okna (poprzez usunięcie wartości NA), ale uśrednienie każdego okna za pomocą odpowiedniego mianownika.
pipefish
1
rn <- cn [(n + 1): długość (cx)] - cx [1: (długość (cx) - n)] powinna faktycznie być rn <- cn [(n + 1): długość (cx)] - cn [1: (length (cx) - n)]
adrianmcmenamin
22

W data.table 1.12.0 nowa frollmeanfunkcja została dodana do obliczenia szybki i dokładny toczenia średnią starannie obsługi NA, NaNa +Inf, -Infwartości.

Ponieważ w pytaniu nie ma powtarzalnego przykładu, nie ma tu wiele do rozwiązania.

Więcej informacji na temat ?frollmeaninstrukcji można znaleźć również w Internecie pod adresem ?frollmean.

Przykłady z instrukcji poniżej:

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
jangorecki
źródło
10

caToolsPakiet bardzo szybko toczenia średni / min / max / SD oraz kilka innych funkcji. Pracowałem tylko z nimi runmeani runsdsą one najszybsze z innych wymienionych do tej pory pakietów.

eddi
źródło
1
To jest niesamowite! Jest to jedyna funkcja, która robi to w ładny, prosty sposób. I jest teraz 2018 ...
Felipe Gerard
9

Możesz użyć RcppRolldo bardzo szybkich średnich ruchomych napisanych w C ++. Wystarczy wywołać roll_meanfunkcję. Dokumenty można znaleźć tutaj .

W przeciwnym razie ta (wolniejsza) pętla powinna załatwić sprawę:

ma <- function(arr, n=15){
  res = arr
  for(i in n:length(arr)){
    res[i] = mean(arr[(i-n):i])
  }
  res
}
cantdutchthis
źródło
3
Czy możesz wyjaśnić mi szczegółowo, jak działa ten algorytm? Ponieważ nie rozumiem pomysłu
Daniel Jefimow
Najpierw inicjuje wektor o tej samej długości res = arr. Następnie jest pętla, która iteruje od początku nlub do 15-tego elementu do końca tablicy. oznacza to, że pierwszym podzbiorem, który bierze pod uwagę, jest to, arr[1:15]które zajmuje miejsce res[15]. Teraz wolę ustawiać res = rep(NA, length(arr))zamiast, res = arrwięc każdy element res[1:14]równa się NA niż liczba, w której nie moglibyśmy wziąć pełnej średniej 15 elementów.
Evan Friedland
7

W rzeczywistości RcppRolljest bardzo dobry.

Kod wysłany przez cantdutchthis musi zostać poprawiony w czwartej linii do okna, aby naprawić:

ma <- function(arr, n=15){
  res = arr
  for(i in n:length(arr)){
    res[i] = mean(arr[(i-n+1):i])
  }
  res
}

Podano tutaj inny sposób radzenia sobie z brakami .

Trzeci sposób udoskonalenia tego kodu do obliczania średnich cząstkowych lub nie:

  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,...)]
  }
}
Rodrigo Remedio
źródło
5

W celu uzupełnienia odpowiedzi cantdutchthis i 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)
Cristóbal Alcázar
źródło
2

Oto przykładowy kod pokazujący, jak obliczyć środkową średnią ruchomą i końcową średnią ruchomą za pomocą rollmeanfunkcji z pakietu zoo .

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
Lubię kodować
źródło
1

Jest to trochę powolne, ale można również użyć zoo :: rollapply do wykonywania obliczeń na macierzach.

reqd_ma <- rollapply(x, FUN = mean, width = n)

gdzie x to zbiór danych, FUN = średnia jest funkcją; możesz także zmienić to na min, max, sd itp., a szerokość to przesuwane okno.

Garima gulati
źródło
2
To nie jest powolne ;. Porównując go do podstawy R, jest znacznie szybszy. set.seed(123); x <- rnorm(1000); system.time(apply(embed(x, 5), 1, mean)); library(zoo); system.time(rollapply(x, 5, mean)) Na mojej maszynie jest tak szybki, że zwraca czas 0 sekund.
G. Grothendieck,
1

Można użyć runnerpakietu do przenoszenia funkcji. W tym przypadku mean_runfunkcja. Problem cummeanpolega na tym, że nie obsługuje NAwartości, ale je mean_runobsługuje. runnerpakiet obsługuje również nieregularne szeregi czasowe, a okna mogą zależeć od daty:

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

Można również określić inne opcje, takie jak lag, i rzucić tylko atokreślone indeksy. Więcej w dokumentacji pakietu i funkcji .

GoGonzo
źródło
1

Można do tego użyć pakietu suwaka. Ma interfejs, który został specjalnie zaprojektowany, aby przypominał mruczenie. Przyjmuje dowolną funkcję i może zwrócić dowolny typ wyniku. Ramki danych są nawet iterowane względem wiersza. Strona pkgdown jest tutaj .

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

Obciążenie zarówno suwaka, jak i tabeli danych frollapply()powinno być dość niskie (znacznie szybsze niż zoo). frollapply()wygląda na to, że jest nieco szybszy w tym prostym przykładzie tutaj, ale należy pamiętać, że wymaga on tylko wprowadzania liczbowego, a wynik musi być skalarną wartością liczbową. funkcje suwaka są całkowicie ogólne i można zwrócić dowolny typ danych.

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
Davis Vaughan
źródło