「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。目下、データ分析・語学に力点を置いています。今月(2022年10月)からは多忙につき、日々の投稿数を減らします。

MENU

時系列解析の基礎(22/XX)

 以下の書籍

を中心に時系列解析を勉強していきます。

13. 線形・Gauss型状態空間モデルの逐次解法

13.1. カルマンフィルタ

 線形・\mathrm{Gauss}型状態空間モデルにおいて、推定すべきデータの真の値と点推定値の間の平均二乗誤差を最小にする意味で最適な逐次推定法はカルマンフィルタと呼ばれる。カルマンフィルタは非定常過程でも扱うことができる。カルマンフィルタで行なわれる計算はすべて線形計算で、正規分布の再生性から関心の対象となる分布はすべて正規分布である。

 カルマンフィルタにおいて前提とする線形・\mathrm{Gauss}型状態空間モデルは


\begin{aligned}
\boldsymbol{x}_t&=\boldsymbol{G}_t\boldsymbol{x}_{t-1}+\boldsymbol{w}_t,&\boldsymbol{w}_t\sim\mathcal{N}\left(\boldsymbol{0},\boldsymbol{W}_t\right)\\
y_t&=\boldsymbol{F}_t\boldsymbol{x}_{t}+v_t,&v_t\sim N\left(0,V_t\right)
\end{aligned}

で、\boldsymbol{G}_tp次状態遷移行列、\boldsymbol{F}_t1\times p[観測行列、\boldsymbol{W}_tp次状態雑音の共分散行列、そしてV_tは観測雑音の分散である。また\boldsymbol{x}_t\sim\mathcal{N}\left(\boldsymbol{m}_0,\boldsymbol{C}_0\right)で、\boldsymbol{m}_0は事前分布におけるp次元平均ベクトルで、\boldsymbol{C}_0p次分散共分散行列である。したがってこのモデルにおける推定パラメータ\boldsymbol{\theta}


\begin{aligned}
\boldsymbol{\theta}=\left\{\boldsymbol{G}_t,\boldsymbol{F}_t,\boldsymbol{W}_t,V_t,\boldsymbol{m}_0,\boldsymbol{C}_0\right\}
\end{aligned}

である。

13.1.1 カルマンフィルタリング

 カルマンフィルタリングにおけるフィルタリング分布、一期先予測分布および一期先予測尤度を

  フィルタリング分布N\left(\boldsymbol{m}_t,\boldsymbol{C}_t\right): 平均ベクトル\boldsymbol{m}_t,分散共分散行列\boldsymbol{C}_t
  一期先予測分布N\left(\boldsymbol{a}_t,\boldsymbol{R}_t\right): 平均ベクトル\boldsymbol{a}_t,分散共分散行列\boldsymbol{R}_t
  一期先予測尤度N\left(f_t,Q_t\right): 平均f_t,分散Q_t

と書くものとする。
 時点t-1でのフィルタリング分布からその翌時点t時点におけるフィルタリング分布は、以下のアルゴリズムで推定する。

    0. 時点t-1におけるフィルタリング分布パラメータとして

平均ベクトル\boldsymbol{m}_t,分散共分散行列\boldsymbol{C}_tを設定する。
    1. 時点tにおける更新を以下で行う:

      一期先予測分布

 (平均)\boldsymbol{a}_t\leftarrow\boldsymbol{G}_t\boldsymbol{m}_{t-1}

 (分散)\boldsymbol{R}_t\leftarrow\boldsymbol{G}_t\boldsymbol{C}_{t-1}{}^{t}\boldsymbol{G}_t+\boldsymbol{W}_t

一期先予測尤度

 (平均)f_t\leftarrow\boldsymbol{F}_t\boldsymbol{a}_t

 (分散)Q_t\leftarrow\boldsymbol{F}_t\boldsymbol{R}_t{}^{t}\boldsymbol{F}_t+V_t

カルマン利得

 \boldsymbol{K}_t\leftarrow\boldsymbol{R}_t{}^{t}\boldsymbol{F}_tQ_t^{-1}

状態の更新

 (平均)  \boldsymbol{m}_t\leftarrow\boldsymbol{a}_t+\boldsymbol{K}_t\left[y_t-f_t\right]

 (分散)\boldsymbol{C}_t\leftarrow\left[\boldsymbol{I}-\boldsymbol{K}_t\boldsymbol{F}_t\right]\boldsymbol{R}_t
    2. 時点tでのフィルタリング分布のパラメータを\boldsymbol{m}_t,\boldsymbol{C}_tと設定する。

これをt=1から繰り返すことで全時点のフィルタリング分布を求める。
 なおカルマン利得は、一期先予測分布を補正してフィルタリング分布を求める際の補正の程度を表す量である。具体的には、観測値y_tと一期先予測尤度の平均f_tとの差を予測誤差e_t=y_t-f_tと置くときに、補正において予測誤差をどの程度反映するかのウェイトという意味合いを持つ。カルマン利得の値には\boldsymbol{W}_tV_tの比が支配的な影響を及ぼす。

############################
### 自作カルマンフィルタ ###
############################

library("ggplot2")

# データ
y <- Nile
t_max <- length(y)

# 1時点分のカルマンフィルタリングを行なう関数
Kalman_filtering <- function(m_t_minus_1, C_t_minus_1, t){
  # 一期先予測分布
  a_t <- G_t %*% m_t_minus_1
  R_t <- G_t %*% C_t_minus_1 %*% t(G_t) + W_t
  
  # 一期先予測尤度
  f_t <- F_t %*% a_t
  Q_t <- F_t %*% R_t %*% t(F_t) + V_t
  
  # カルマン利得
  K_t <- R_t %*% t(F_t) %*% solve(Q_t)
  
  # 状態の更新
  m_t <- a_t + K_t %*% (y[t]-f_t)
  C_t <- (diag(nrow(R_t)) - K_t %*% F_t) %*% R_t
  
  # フィルタリング分布
  list(m = m_t, C = C_t,
       a = a_t, R = R_t)
}

# 線形ガウス型状態空間のパラメータを設定
G_t <- matrix(1, ncol = 1, nrow = 1)
W_t <- matrix(exp(7.29), ncol = 1, nrow = 1)
F_t <- matrix(1, ncol = 1, nrow = 1)
V_t <- matrix(exp(9.62), ncol = 1, nrow = 1)
m0 <- matrix(0, ncol = 1, nrow = 1)
C0 <- matrix(1e+7, ncol = 1, nrow = 1)

# フィルタリング分布
# 状態の領域を確保
m <- rep(NA_real_, t_max)
C <- rep(NA_real_, t_max)
a <- rep(NA_real_, t_max)
R <- rep(NA_real_, t_max)

# 時点t=1
KF <- Kalman_filtering(m_t_minus_1 = m0,C_t_minus_1 = C0,t = 1)
m[1] <- KF$m
C[1] <- KF$C
a[1] <- KF$a
R[1] <- KF$R

#時点t=2-t_max
for (t in 2:t_max){
  KF <- Kalman_filtering(m_t_minus_1 = m[t-1],C_t_minus_1 = C[t-1],t = t)
  m[t] <- KF$m
  C[t] <- KF$C
  a[t] <- KF$a
  R[t] <- KF$R
}

# フィルタリングの結果
filtering_data <- data.frame(time = as.numeric(time(y)),
                             var = "観測値",
                             val = as.numeric(y))
filtering_data <- rbind(filtering_data,data.frame(time = as.numeric(time(y)),
                                                  var = "平均値",
                                                  val = m))
filtering_data <- rbind(filtering_data,data.frame(time = as.numeric(time(y)),
                                                  var = "下側95%区間",
                                                  val = mapply(FUN = function(m,v)qnorm(p = 0.95,mean = m,sd = sqrt(v)),m,C)))
filtering_data <- rbind(filtering_data,data.frame(time = as.numeric(time(y)),
                                                  var = "下側 5%区間",
                                                  val = mapply(FUN = function(m,v)qnorm(p = 0.05,mean = m,sd = sqrt(v)),m,C)))

filtering_data$var <- factor(filtering_data$var,levels = c("観測値","平均値","下側95%区間","下側 5%区間"))


# グラフ
g <- ggplot(data = filtering_data,mapping = aes(x = time, y = val,group = var, color = var, linetype = var))
g <- g + geom_line(size = 0.6) + theme_bw()
g <- g + labs(title = "ナイル川の流量データに対するカルマンフィルタ適用結果",
              x = "年",y = paste0("ナイル川の流量[",expression(10^8),expression(m^3),"]"),group = "",color = "",
              linetype = "")
g <- g + theme(axis.text= element_text(size = 12),
               axis.text.x=element_text(size = 12),
               legend.title = element_blank(),
               legend.text = element_text(size = 12),
               legend.position = "top",
               axis.title = element_text(size = 12),
               axis.title.y = element_text(angle = 90,size = 12),
               plot.title= element_text(size = 16,hjust = 0.5))
g <- g + scale_linetype_manual(values = c("solid","solid","dashed","dashed"))
g <- g + scale_color_manual(values = c("black","red","blue","blue"))
plot(g)


プライバシーポリシー お問い合わせ