6.20 计算移动平均数

6.20.1 问题

你想要计算移动平均数。

6.20.2 解决方案

假设你的数据是带缺失值的噪声正弦波。

  1. set.seed(993)
  2. x <- 1:300
  3. y <- sin(x/20) + rnorm(300, sd = 0.1)
  4. y[251:255] <- NA

filter() 函数可以用来计算移动平均数。

  1. # 绘制未平滑的数据（灰色）
  2. plot(x, y, type = "l", col = grey(0.5))
  3. # 绘制网格
  4. grid()
  6. # 延迟平滑 当前样本和之前 19 个样本的平均数（红色）
  7. f20 <- rep(1/20, 20)
  8. f20
  9. #>  [1] 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05
  10. #> [11] 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05
  11. y_lag <- filter(y, f20, sides = 1)
  12. lines(x, y_lag, col = "red")
  15. # 对称性平滑：
  16. # 计算当前样本，时间上将来样本和过去10个样本的平均数（蓝色）
  17. f21 <- rep(1/21, 21)
  18. f21
  19. #>  [1] 0.04762 0.04762 0.04762 0.04762 0.04762 0.04762
  20. #>  [7] 0.04762 0.04762 0.04762 0.04762 0.04762 0.04762
  21. #> [13] 0.04762 0.04762 0.04762 0.04762 0.04762 0.04762
  22. #> [19] 0.04762 0.04762 0.04762
  23. y_sym <- filter(y, f21, sides = 2)
  24. lines(x, y_sym, col = "blue")

6.20 计算移动平均数 - 图1

filter() 会在遭遇缺失值时留下空缺，就像上面图中显示的一样。

一种处理缺失值的方法是简单地忽略它，不把它包含在平均数的计算中。这个功能可以由下面的函数实现：

  1. # x: the vector n: the number of samples centered: if
  2. # FALSE, then average current sample and previous (n-1)
  3. # samples if TRUE, then average symmetrically in past
  4. # and future. (If n is even, use one more sample from
  5. # future.)
  6. movingAverage <- function(x, n = 1, centered = FALSE) {
  8.   if (centered) {
  9.     before <- floor((n - 1)/2)
  10.     after <- ceiling((n - 1)/2)
  11.   } else {
  12.     before <- n - 1
  13.     after <- 0
  14.   }
  16.   # Track the sum and count of number of non-NA items
  17.   s <- rep(0, length(x))
  18.   count <- rep(0, length(x))
  20.   # Add the centered data
  21.   new <- x
  22.   # Add to count list wherever there isn't a
  23.   count <- count + !is.na(new)
  24.   # Now replace NA_s with 0_s and add to total
  25.   new[is.na(new)] <- 0
  26.   s <- s + new
  28.   # Add the data from before
  29.   i <- 1
  30.   while (i <= before) {
  31.     # This is the vector with offset values to add
  32.     new <- c(rep(NA, i), x[1:(length(x) - i)])
  34.     count <- count + !is.na(new)
  35.     new[is.na(new)] <- 0
  36.     s <- s + new
  38.     i <- i + 1
  39.   }
  41.   # Add the data from after
  42.   i <- 1
  43.   while (i <= after) {
  44.     # This is the vector with offset values to add
  45.     new <- c(x[(i + 1):length(x)], rep(NA, i))
  47.     count <- count + !is.na(new)
  48.     new[is.na(new)] <- 0
  49.     s <- s + new
  51.     i <- i + 1
  52.   }
  54.   # return sum divided by count
  55.   s/count
  56. }
  59. # 用比较厚的线条绘制和之前一样的图
  60. plot(x, y, type = "l", col = grey(0.5))
  61. grid()
  62. y_lag <- filter(y, rep(1/20, 20), sides = 1)
  63. lines(x, y_lag, col = "red", lwd = 4)  # 用红色表示延迟平均
  64. y_sym <- filter(y, rep(1/21, 21), sides = 2)
  65. lines(x, y_sym, col = "blue", lwd = 4)  # 用蓝色表示对称平均
  67. # 用上面定义的函数计算延迟平均
  68. y_lag_na.rm <- movingAverage(y, 20)
  69. lines(x, y_lag_na.rm, col = "green", lwd = 2)
  71. # 用上面定义的函数计算对称性平均
  72. y_sym_na.rm <- movingAverage(y, 21, TRUE)
  73. lines(x, y_sym_na.rm, col = "green", lwd = 2)

6.20 计算移动平均数 - 图2