Calculating a moving average

Problem

You want to calculate a moving average.

Solution

Suppose your data is a noisy sine wave with some missing values:

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

The filter() function can be used to calculate a moving average.

  1. # Plot the unsmoothed data (gray)
  2. plot(x, y, type="l", col=grey(.5))
  3. # Draw gridlines
  4. grid()
  6. # Smoothed with lag:
  7. # average of current sample and 19 previous samples (red)
  8. f20 <- rep(1/20, 20)
  9. f20
  10. #>  [1] 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05
  11. #> [18] 0.05 0.05 0.05
  12. y_lag <- filter(y, f20, sides=1)
  13. lines(x, y_lag, col="red")
  15. # Smoothed symmetrically:
  16. # average of current sample, 10 future samples, and 10 past samples (blue)
  17. f21 <- rep(1/21,21)
  18. f21
  19. #>  [1] 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905
  20. #>  [8] 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905
  21. #> [15] 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905 0.04761905
  22. y_sym <- filter(y, f21, sides=2)
  23. lines(x, y_sym, col="blue")

plot of chunk unnamed-chunk-2

filter() will leave holes wherever it encounters missing values, as shown in the graph above.

A different way to handle missing data is to simply ignore it, and not include it in the average. The function defined here will do that.

  1. # x: the vector
  2. # n: the number of samples
  3. # centered: if FALSE, then average current sample and previous (n-1) samples
  4. #           if TRUE, then average symmetrically in past and future. (If n is even, use one more sample from future.)
  5. movingAverage <- function(x, n=1, centered=FALSE) {
  7.     if (centered) {
  8.         before <- floor  ((n-1)/2)
  9.         after  <- ceiling((n-1)/2)
  10.     } else {
  11.         before <- n-1
  12.         after  <- 0
  13.     }
  15.     # Track the sum and count of number of non-NA items
  16.     s     <- rep(0, length(x))
  17.     count <- rep(0, length(x))
  19.     # Add the centered data 
  20.     new <- x
  21.     # Add to count list wherever there isn't a 
  22.     count <- count + !is.na(new)
  23.     # Now replace NA_s with 0_s and add to total
  24.     new[is.na(new)] <- 0
  25.     s <- s + new
  27.     # Add the data from before
  28.     i <- 1
  29.     while (i <= before) {
  30.         # This is the vector with offset values to add
  31.         new   <- c(rep(NA, i), x[1:(length(x)-i)])
  33.         count <- count + !is.na(new)
  34.         new[is.na(new)] <- 0
  35.         s <- s + new
  37.         i <- i+1
  38.     }
  40.     # Add the data from after
  41.     i <- 1
  42.     while (i <= after) {
  43.         # This is the vector with offset values to add
  44.         new   <- c(x[(i+1):length(x)], rep(NA, i))
  46.         count <- count + !is.na(new)
  47.         new[is.na(new)] <- 0
  48.         s <- s + new
  50.         i <- i+1
  51.     }
  53.     # return sum divided by count
  54.     s/count
  55. }
  58. # Make same plots from before, with thicker lines
  59. plot(x, y, type="l", col=grey(.5))
  60. grid()
  61. y_lag <- filter(y, rep(1/20, 20), sides=1)
  62. lines(x, y_lag, col="red", lwd=4)         # Lagged average in red
  63. y_sym <- filter(y, rep(1/21,21), sides=2)
  64. lines(x, y_sym, col="blue", lwd=4)        # Symmetrical average in blue
  66. # Calculate lagged moving average with new method and overplot with green
  67. y_lag_na.rm <- movingAverage(y, 20)
  68. lines(x, y_lag_na.rm, col="green", lwd=2)
  70. # Calculate symmetrical moving average  with new method and overplot with green
  71. y_sym_na.rm <- movingAverage(y, 21, TRUE)
  72. lines(x, y_sym_na.rm, col="green", lwd=2)

plot of chunk unnamed-chunk-3