Calculando el promedio móvil

186

Estoy tratando de usar R para calcular el promedio móvil sobre una serie de valores en una matriz. Sin embargo, la búsqueda normal de la lista de correo R no ha sido muy útil. Parece que no hay una función incorporada en R que me permita calcular promedios móviles. ¿Algún paquete proporciona uno? ¿O necesito escribir el mío?

Jared
fuente

Respuestas:

141
  • Rolling Means / Maximums / Medians en el paquete del zoológico (rollmean)
  • Promedios móviles en TTR
  • ma en previsión
f3lix
fuente
1
¿Cuál es el promedio móvil en R que no contiene valores futuros de una marca de tiempo dada? Lo revisé forecast::may contiene todo el vecindario, no está bien.
hhh
214

O simplemente puede calcularlo usando filtro, aquí está la función que uso:

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

Si usa dplyr, tenga cuidado de especificar stats::filteren la función anterior.

Matti Pastell
fuente
49
Debo señalar que "lados = 2" puede ser una opción importante en los casos de uso de muchas personas que no quieren pasar por alto. Si solo desea información final en su promedio móvil, debe usar lados = 1.
evanrsparks
36
Algunos años más tarde, pero ahora dplyr tiene una función de filtro, si usted tiene este paquete cargado usostats::filter
blmoore
sides = 2es equivalente a align = "center" para el zoo :: rollmean o RcppRoll :: roll_mean. sides = 1es equivalente a la alineación "correcta". ¿No veo una manera de hacer la alineación "izquierda" o calcular con datos "parciales" (2 o más valores)?
Matt L.
29

El uso cumsumdebe ser suficiente y eficiente. Suponiendo que tiene un vector x y desea una suma de n números

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

Como se señaló en los comentarios de @mzuther, esto supone que no hay NA en los datos. lidiar con ellos requeriría dividir cada ventana por el número de valores que no son NA. Aquí hay una forma de hacerlo, incorporando el comentario de @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

Esto todavía tiene el problema de que si todos los valores en la ventana son NA, habrá un error de división por cero.

aguja
fuente
8
Una desventaja de esta solución es que no puede manejar las faltas:cumsum(c(1:3,NA,1:3))
Jthorpe
Puedes hacerlo fácilmente manejando NAs haciendo cx <- c(0, cumsum(ifelse(is.na(x), 0, x))).
Ricardo Cruz
@ Ricardo Cruz: podría ser mejor eliminar las NA y ajustar la longitud del vector en consecuencia. Piense en un vector con una gran cantidad de NA: los ceros llevarán el promedio hacia cero, mientras que la eliminación de los NA dejará el promedio tal como está. Todo depende de sus datos y la pregunta que desea responder, por supuesto. :)
mzuther
@mzuther, actualicé la respuesta después de tus comentarios. Gracias por el aporte. Creo que la forma correcta de tratar con los datos faltantes no es extender la ventana (eliminando los valores de NA), sino promediando cada ventana por el denominador correcto.
pez pipa
1
rn <- cn [(n + 1): longitud (cx)] - cx [1: (longitud (cx) - n)] en realidad debería ser rn <- cn [(n + 1): longitud (cx)] - cn [1: (longitud (cx) - n)]
adrianmcmenamin
22

En data.table 1.12.0 nueva frollmeanfunción se ha añadido para calcular rápida y exacta rodar media cuidadosamente el manejo NA, NaNy +Inf, -Infvalores.

Como no hay un ejemplo reproducible en la pregunta, no hay mucho más que abordar aquí.

Puede encontrar más información sobre ?frollmeanen el manual, también disponible en línea en ?frollmean.

Ejemplos del manual a continuación:

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
fuente
10

El caToolspaquete tiene una velocidad de desplazamiento muy rápida / min / max / sd y algunas otras funciones. Solo he trabajado con runmeany runsdson los más rápidos de los otros paquetes mencionados hasta la fecha.

eddi
fuente
1
¡Esto es asombroso! Es la única función que hace esto de una manera agradable y simple. Y es 2018 ahora ...
Felipe Gerard
9

Puede usar RcppRollpara promedios móviles muy rápidos escritos en C ++. Solo llama a la roll_meanfunción. Los documentos se pueden encontrar aquí .

De lo contrario, este bucle for (más lento) debería ser el truco:

ma <- function(arr, n=15){
  res = arr
  for(i in n:length(arr)){
    res[i] = mean(arr[(i-n):i])
  }
  res
}
Cantdutchthis
fuente
3
¿Pueden explicarme en detalle cómo funciona este algoritmo? Porque no puedo entender la idea
Daniel Yefimov
Primero inicializa un vector de la misma longitud con res = arr. Luego hay un ciclo que itera comenzando en no, el decimoquinto elemento, hasta el final de la matriz. eso significa que el primer subconjunto del que toma la media es el arr[1:15]que llena el lugar res[15]. Ahora, prefiero establecer en res = rep(NA, length(arr))lugar de res = arrque cada elemento sea res[1:14]igual a NA en lugar de un número, donde no podríamos tomar un promedio completo de 15 elementos.
Evan Friedland
7

De hecho RcppRolles muy bueno.

El código publicado por cantdutchthis debe corregirse en la cuarta línea para que se fije la ventana:

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

Otra forma, que se ocupa de missings, se da aquí .

Una tercera forma, mejorar el código cantdutchthis para calcular promedios parciales o no, sigue:

  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
fuente
5

Para complementar la respuesta de cantdutchthis y 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
fuente
2

Aquí hay un código de ejemplo que muestra cómo calcular un promedio móvil centrado y un promedio móvil final utilizando la rollmeanfunción del paquete del zoológico .

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
Me gusta codificar
fuente
1

Aunque es un poco lento, también puede usar zoo :: rollapply para realizar cálculos en matrices.

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

donde x es el conjunto de datos, FUN = mean es la función; También puede cambiarlo a min, max, sd, etc. y el ancho es la ventana móvil.

Garima gulati
fuente
2
No es lento; Comparándolo con la base R, es mucho más rápido. set.seed(123); x <- rnorm(1000); system.time(apply(embed(x, 5), 1, mean)); library(zoo); system.time(rollapply(x, 5, mean)) En mi máquina es tan rápido que devuelve un tiempo de 0 segundos.
G. Grothendieck
1

Se puede usar el runnerpaquete para mover funciones. En este caso la mean_runfunción. El problema cummeanes que no maneja NAvalores, pero mean_runsí. runnerEl paquete también admite series temporales irregulares y las ventanas pueden depender de la fecha:

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

También se pueden especificar otras opciones como lag, y rodar solo atíndices específicos. Más en la documentación del paquete y la función .

GoGonzo
fuente
1

El paquete deslizante se puede utilizar para esto. Tiene una interfaz que ha sido diseñada específicamente para sentirse similar a ronronear. Acepta cualquier función arbitraria y puede devolver cualquier tipo de salida. Los marcos de datos incluso se repiten en fila. El sitio pkgdown está aquí .

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

La sobrecarga de los controles deslizantes y data.table frollapply()debe ser bastante baja (mucho más rápido que el zoológico). frollapply()Parece ser un poco más rápido para este ejemplo simple aquí, pero tenga en cuenta que solo se necesita una entrada numérica, y la salida debe ser un valor numérico escalar. Las funciones del control deslizante son completamente genéricas y puede devolver cualquier tipo de datos.

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
fuente
0
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))
}
Mohamed Galia
fuente
2
Por favor agregue una descripción para más detalles.
Farbod Ahmadian
Relacione su respuesta con la pregunta e incluya algún resultado que muestre que la pregunta ha sido respondida. Consulte Cómo responder para obtener orientación sobre cómo dar una buena respuesta.
Peter