KalmanFilter (univariate) with ggplot2 animation Pt2

KalmanFilter02.R

set.seed(1)

# 観測系列のサンプルサイズ
n <- 200

# 真の値が時間変化する
actual <- 3.0 + cumsum(rnorm(n, sd = 0.2))

# 観測される値 (誤差は標準偏差2の正規分布とする)
observed <- actual + rnorm(n, mean = 0, sd = 2)

# 結果保存用のvector
xhat <- rep(0, length.out = n)
P <- rep(0, length.out = n)
K <- rep(0, length.out = n)

# 誤差
Q = 1e-5
R = 0.1 ** 2

for (k in seq(2, n)) {
  # predict
  xhat.m <- xhat[k-1]
  P.m <- P[k-1] + Q

  # update
  S <- R + P.m
  K[k] <- P.m / S
  xhat[k] <- xhat.m + K[k] * (observed[k] - xhat.m)
  P[k] <- (1 - K[k]) * P.m
}

xhat

library(animation)
library(ggplot2)

d <- data.frame(x = seq(1, n), actual = actual,
                observed = observed, 
                fitted = xhat, P = P, K = K)

saveGIF({
  for (i in seq(2, n, by = 2)) {
    tmp <- head(d, i)
    p <- ggplot(tmp, aes(x = x)) +
      geom_point(aes(y = observed)) +
      geom_line(aes(y = fitted), colour = 'blue') +
      annotate(geom = 'text', x = i, y = tmp$fitted[i] - 0.5,
               label = 'Fitted value', colour = 'blue', hjust = 1) +
      geom_line(aes(y = actual), colour = 'red') +
      annotate(geom = 'text', x = i, y = tmp$actual[i] + 0.5,
               label = 'True value', colour = 'red', hjust = 1) +
      ylim(-1, 8) +
      scale_x_continuous(breaks = seq(0, n, by = 20)) +
      xlab('') + ylab('')
    print(p)
  }
}, interval = 0.2, movie.name = "kalmanfilter02_01.gif",
ani.width = 600, ani.height = 400)

library(tidyr)

saveGIF({
  for (i in seq(2, n, by = 2)) {
    tmp <- head(d, i)
    tmp <- gather(tmp, 'variable', 'value', c(P, K))
    p <- ggplot(tmp, aes(x = x, y = value, colour = variable)) +
      geom_line() + 
      facet_wrap(~ variable, ncol = 1, scale = 'free_y') +
      scale_x_continuous(breaks = seq(0, n, by = 20)) +
      xlab('') + ylab('')
    print(p)
  }
}, interval = 0.2, movie.name = "kalmanfilter02_02.gif",
ani.width = 600, ani.height = 400)