Menghitung moving average

185

Saya mencoba menggunakan R untuk menghitung rata-rata bergerak atas serangkaian nilai dalam sebuah matriks. Pencarian milis R normal belum sangat membantu. Sepertinya tidak ada fungsi bawaan di R yang akan memungkinkan saya menghitung rata-rata bergerak. Apakah ada paket yang menyediakannya? Atau apakah saya perlu menulis sendiri?

Jared
sumber

Jawaban:

140
f3lix
sumber
1
Berapa rata-rata bergerak dalam R yang tidak mengandung nilai timestamp yang akan datang? Saya memeriksa forecast::madan itu berisi semua lingkungan, tidak benar.
hhh
214

Atau Anda cukup menghitungnya menggunakan filter, inilah fungsi yang saya gunakan:

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

Jika Anda menggunakan dplyr, berhati-hatilah untuk menentukan stats::filterfungsi di atas.

Matti Pastell
sumber
49
Saya harus menunjukkan bahwa "sisi = 2" mungkin merupakan opsi penting dalam kasus penggunaan banyak orang yang tidak ingin mereka abaikan. Jika Anda hanya ingin mengekor informasi dalam moving average Anda, Anda harus menggunakan sisi = 1.
evanrsparks
36
Beberapa tahun kemudian tetapi dplyr sekarang memiliki fungsi filter, jika Anda menggunakan paket inistats::filter
blmoore
sides = 2setara dengan menyelaraskan = "pusat" untuk kebun binatang :: rollmean atau RcppRoll :: roll_mean. sides = 1setara dengan perataan "benar". Saya tidak melihat cara untuk melakukan perataan "kiri" atau menghitung dengan data "parsial" (2 nilai atau lebih)?
Matt L.
29

Penggunaan cumsumharus memadai dan efisien. Dengan asumsi Anda memiliki vektor x dan Anda ingin jumlah n angka yang berjalan

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

Sebagaimana ditunjukkan dalam komentar oleh @mzuther, ini mengasumsikan bahwa tidak ada NAS dalam data. untuk mengatasinya akan membutuhkan pembagian setiap jendela dengan jumlah nilai non-NA. Inilah satu cara untuk melakukan itu, memasukkan komentar dari @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

Ini masih memiliki masalah bahwa jika semua nilai di jendela adalah NA maka akan ada pembagian dengan kesalahan nol.

pipefish
sumber
8
Satu kelemahan dari solusi ini adalah tidak dapat menangani kesalahan:cumsum(c(1:3,NA,1:3))
Jthorpe
Anda dapat dengan mudah membuatnya menangani NAS dengan melakukan cx <- c(0, cumsum(ifelse(is.na(x), 0, x))).
Ricardo Cruz
@ Ricardo Cruz: mungkin lebih baik untuk menghapus NAs dan menyesuaikan panjang vektor yang sesuai. Pikirkan vektor dengan banyak NAs - nol akan menarik rata-rata ke nol, sementara menghapus NAs akan meninggalkan rata-rata seperti semula. Itu semua tergantung pada data Anda dan pertanyaan yang ingin Anda jawab, tentu saja. :)
mzuther
@muther, saya memperbarui jawaban setelah komentar Anda. Terima kasih atas masukannya. Saya pikir cara yang benar untuk menangani data yang hilang tidak memperluas jendela (dengan menghapus nilai-nilai NA), tetapi dengan rata-rata setiap jendela dengan penyebut yang benar.
pipefish
1
rn <- cn [(n + 1): panjang (cx)] - cx [1: (panjang (cx) - n)] harus benar-benar rn <- cn [(n + 1): panjang (cx)] - cn [1: (panjang (cx) - n)]
adrianmcmenamin
22

Dalam data.tabel 1.12.0frollmean fungsi baru telah ditambahkan untuk menghitung rolling yang cepat dan tepat berarti penanganan NA, NaNdan +Inf, -Infnilai-nilai secara cermat .

Karena tidak ada contoh yang dapat direproduksi dalam pertanyaan, maka tidak banyak lagi yang perlu dibahas di sini.

Anda dapat menemukan lebih banyak info tentang ?frollmeansecara manual, juga tersedia online di ?frollmean.

Contoh dari manual di bawah ini:

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

The caToolspaket telah sangat cepat bergulir berarti / min / max / sd dan beberapa fungsi lainnya. Saya hanya bekerja dengan runmeandan runsddan mereka adalah yang tercepat dari paket lain yang disebutkan sampai saat ini.

eddi
sumber
1
Ini luar biasa! Ini adalah satu-satunya fungsi yang melakukan ini dengan cara yang bagus dan sederhana. Dan sekarang tahun 2018 ...
Felipe Gerard
9

Anda bisa menggunakan RcppRollrata-rata bergerak sangat cepat yang ditulis dalam C ++. Panggil saja roll_meanfungsinya. Documents dapat ditemukan di sini .

Kalau tidak, ini (lebih lambat) untuk loop harus melakukan trik:

ma <- function(arr, n=15){
  res = arr
  for(i in n:length(arr)){
    res[i] = mean(arr[(i-n):i])
  }
  res
}
tidak bisa ini
sumber
3
Bisakah Anda jelaskan saya secara detail, bagaimana cara kerja algoritma ini? Karena saya tidak dapat memahami ide
Daniel Yefimov
Pertama dia menginisialisasi vektor dengan panjang yang sama res = arr. Lalu ada loop yang dimulai mulai dari natau, elemen ke-15, ke ujung array. itu berarti bagian pertama yang diambilnya adalah arr[1:15]yang mengisi titik res[15]. Sekarang, saya lebih suka pengaturan res = rep(NA, length(arr))daripada res = arrsetiap elemen res[1:14]sama dengan NA daripada angka, di mana kami tidak bisa mengambil rata-rata penuh dari 15 elemen.
Evan Friedland
7

Padahal RcppRollsangat bagus.

Kode yang diposting oleh cantdutchini harus diperbaiki pada baris keempat ke jendela diperbaiki:

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

Cara lain, yang menangani kerugian, diberikan di sini .

Cara ketiga, memperbaiki kode cantdutch ini untuk menghitung rata-rata parsial atau tidak, mengikuti:

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

Untuk melengkapi jawaban cantdutchthis dan 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
sumber
2

Berikut adalah contoh kode yang menunjukkan cara menghitung rata-rata bergerak terpusat dan rata-rata bergerak tertinggal menggunakan rollmeanfungsi dari paket kebun binatang .

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
Saya Suka Kode
sumber
1

Meskipun agak lambat tetapi Anda juga dapat menggunakan zoo :: rollapply untuk melakukan perhitungan pada matriks.

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

di mana x adalah kumpulan data, FUN = mean adalah fungsinya; Anda juga dapat mengubahnya ke min, maks, sd dll dan lebar adalah jendela bergulir.

Garima gulati
sumber
2
Itu tidak lambat; Membandingkannya dengan basis R, jauh lebih cepat. set.seed(123); x <- rnorm(1000); system.time(apply(embed(x, 5), 1, mean)); library(zoo); system.time(rollapply(x, 5, mean)) Di komputer saya, kecepatannya sangat cepat sehingga mengembalikan waktu 0 detik.
G. Grothendieck
1

Satu dapat menggunakan runnerpaket untuk memindahkan fungsi. Dalam hal ini mean_runberfungsi. Masalah dengan itu cummeanadalah bahwa ia tidak menangani NAnilai, tetapi mean_runtidak. runnerpaket juga mendukung seri waktu yang tidak teratur dan windows dapat bergantung pada tanggal:

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

Satu juga dapat menentukan opsi lain seperti lag, dan hanya menggulung atindeks tertentu. Lebih banyak dalam dokumentasi paket dan fungsi .

GoGonzo
sumber
1

Paket slider dapat digunakan untuk ini. Ini memiliki antarmuka yang telah dirancang khusus agar terasa mirip dengan purrr. Ini menerima fungsi sewenang-wenang, dan dapat mengembalikan segala jenis output. Frame data bahkan diulang lebih dari baris bijaksana. Situs pkgdown ada di sini .

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

Overhead slider dan data.table frollapply()harus cukup rendah (jauh lebih cepat daripada kebun binatang). frollapply()terlihat sedikit lebih cepat untuk contoh sederhana ini di sini, tetapi perhatikan bahwa itu hanya membutuhkan input numerik, dan output harus berupa nilai numerik skalar. fungsi slider sepenuhnya generik, dan Anda dapat mengembalikan tipe data apa pun.

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
sumber
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
sumber
2
Silakan tambahkan deskripsi untuk detail lebih lanjut.
Farbod Ahmadian
Harap kaitkan jawaban Anda dengan pertanyaan dan sertakan beberapa output yang menunjukkan bahwa pertanyaan telah dijawab. Lihat Cara Menjawab untuk panduan cara membuat jawaban yang baik.
Peter