状態空間モデル：(3)推定（平滑化と予測）

はじめに

今回は、以前書いた状態空間モデルの話の続きを書いていきたい。
mstour.hatenablog.com

平滑化(Smoothing)

状態空間モデルのフィルタリングとは、現在までに得られているデータをもとにして、（観測されない）現在の状態を推定するという行為であった。一方、平滑化(Smoothing)というのは、今までに得られたデータを使って過去の状態を推定し直すという作業にあたる。
数式で書くと、現在の時点を $T$ とするとき、現在までの観測値 $y_{1:T}$ を用いて過去の時点 $t < T$ における状態の分布 $p(\theta_t|y_{1:T})$ を推定する、というのが平滑化である。
f:id:mstour:20201029204252j:plain
平滑化についても、フィルタリングと同様にして時点を一つずつ遡る漸化式を利用することができる。状態空間モデルの性質（状態 $\theta_t$ を与えたもとでは次の時点の状態 $\theta_{t+1}$ と $t$ までの観測値 $y_{1:t}$ とは条件付き独立になる）がこのような漸化式の導出を可能にしている。
(1) まず、観測値 $y_{1:T}$ の条件付きでの状態の逆向きの遷移確率が以下のようになる。

$\displaystyle p(\theta_t|\theta_{t+1}, y_{1:T}) = \frac{p(\theta_{t+1}|\theta_t) p(\theta_t|y_{1:t})}{p(\theta_{t+1}|y_{1:t})}$

(2) (1)を $\theta_{t+1}$ について周辺化することで、過去の状態の分布が得られる。

$\displaystyle p(\theta_t|y_{1:T}) = p(\theta_t|y_{1:t}) \int \frac{p(\theta_{t+1}|\theta_t)}{p(\theta_{t+1}|y_{1:t})}p(\theta_{t+1}|y_{1:T})d\theta_{t+1}$

(2)を求めるために必要な分布のうち、フィルタリング分布 $p(\theta_t|y_{1:t})$ 、フィルタリングの計算で用いた一期先の状態の分布 $p(\theta_{t+1}|y_{1:t})$ 、状態方程式のモデル $p(\theta_{t+1}|\theta_t)$ はすでに利用可能であり、 $p(\theta_{t+1}|y_{1:T})$ は時点 $T$ から順次遡っていけば算出することができる。ただし、積分計算が必要とされるので一般的に計算が容易とは言えない。
動的線形モデルの場合には、平滑化分布も正規分布で表すことができる。これはカルマンスムーザー(Kalman smoother)とも呼ばれている。以下の動的線形モデル

$\displaystyle y_t = F_t \theta_t + v_t \hspace{15pt} (v_t \sim N(0, V_t))$

$\displaystyle \theta_t = G_t \theta_{t-1} + w_t \hspace{15pt} (w_t \sim N(0, W_t))$

を考え、フィルタリングによって次の条件付き期待値および分散が得られているものとする。

$\displaystyle m_{t-1} := E(\theta_{t-1}|y_{1:t-1}), \\ C_{t-1} := V(\theta_{t-1}|y_{1:t-1}), \\ a_t := E(\theta_t|y_{1:t-1}) = G_t m_{t-1}, \\ R_t := V(\theta_t|y_{1:t-1}) = G_t C_{t-1} G_t^{T} + W_t, \\ f_t := E(y_t|y_{1:t-1}) = F_t a_t, \\ Q_t := V(y_t|y_{1:t-1}) = F_t R_t F_t^{T} + V_t, \\ m_t := E(\theta_t|y_{1:t}) = a_t + R_t F_t^{T} Q_t^{-1} (y_t - f_t), \\ C_t := V(\theta_t|y_{1:t}) = R_t - R_t F_t^{T} Q_t^{-1} F_t R_t.$

このとき、平滑化分布 $p(\theta_t|y_{1:T})$ は次の平均・分散をもつ正規分布 $N(s_t, S_t)$ で表される：

$\displaystyle s_t = m_t + C_t G_{t+1}^{T} R_{t+1}^{-1} (s_{t+1} - a_{t+1}), \\ S_t = C_t - C_t G_{t+1}^{T} R_{t+1}^{-1} (R_{t+1} - S_{t+1}) R_{t+1}^{-1} G_{t+1} C_t.$

ただし、

$\displaystyle s_{t+1} = E(\theta_{t+1}|y_{1:T}), \\ S_{t+1} = V(\theta_{t+1}|y_{1:T})$

である。

予測(Prediction)

フィルタリングでは計算の過程で状態の一期先の予測分布を導出する。これを利用すれば、将来の任意の時点の予測分布を得ることができる。 $k > 0$ とすると、
(1) 状態の $k$ 期先予測分布は次の通り

$\displaystyle p(\theta_{t+k}|y_{1:t}) = \int p(\theta_{t+k}|\theta_{t+k-1}) p(\theta_{t+k-1}|y_{1:t}) d \theta_{t+k-1}$

右辺の1つめの分布は状態方程式からすぐに得られる。2つめの分布は $k-1$ 期先予測分布なので、一期先予測分布から順次計算をしていけば得られていることになる。
(2) 観測値の $k$ 期先予測分布は次の通り

$\displaystyle p(y_{t+k}|y_{1:t}) = \int p(y_{t+k}|\theta_{t+k}) p(\theta_{t+k}|y_{1:t}) d \theta_{t+k}$

右辺の1つめの分布は観測方程式からすぐに得られる。2つめの分布は(1)で求めた状態の $k$ 期先予測分布である。
予測分布についても、動的線形モデルの場合は具体的な正規分布で書きくだせる。
(1) 状態の $k$ 期先予測分布は、以下の平均と分散をもつ正規分布になる。

$\displaystyle E(\theta_{t+k}|y_{1:t}) = G_{t+k} E(\theta_{t+k-1}|y_{1:t}) \\ V(\theta_{t+k}|y_{1:t}) = G_{t+k} V(\theta_{t+k-1}|y_{1:t}) G_{t+k}^{T} + W_{t+k}$

(2) 観測値の $k$ 期先予測分布は、以下の平均と分散をもつ正規分布になる。

$\displaystyle E(y_{t+k}|y_{1:t}) = F_{t+k} E(\theta_{t+k}|y_{1:t}) \\ V(y_{t+k}|y_{1:t}) = F_{t+k} V(\theta_{t+k}|y_{1:t}) F_{t+k}^{T} + V_{t+k}$

まとめ

2回に分けて、状態空間モデルの推定に関する話題を紹介した。データを得たあとでの状態や観測値の分布の推定は、どのデータを使ってどの時点の推定を行うかによってフィルタリング、平滑化、予測に分類することができる。そして、いずれの場合も状態空間モデルの性質から導かれる漸化式が存在する。動的線形モデルの場合にはそれらはすべて具体的な正規分布で書き表すことができる。
ただし、ここまでに述べた内容は、状態方程式・観測方程式のいずれも構造が既知であり（つまり動的線形モデルでいえば $F_t$ や $G_t$ は未知パラメータでなく固定値である）、かつ誤差のパラメータも既知である（ $V_t$ と $W_t$ が固定値）と想定した上での結果である。データからその背後にある現象をモデリングする場面では多くの場合これらは推定すべきパラメータとして考えたほうがよいはずなので、推定はもっと複雑な問題になってくる。このあたりの話はもう少し勉強しておきたい。

参考文献

[1] Petris他(2013), "Rによるベイジアン動的線型モデル", 朝倉書店.