Calculando el promedio móvil

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?

• Rolling Means / Maximums / Medians en el paquete del zoológico (rollmean)
• Promedios móviles en TTR
• ma en previsión

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.

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.

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

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.

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
}

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,...)]
  }
}

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)

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

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.

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 .

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

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