状態空間モデル：(2)推定（概要とフィルタリングについて）

はじめに

前回、状態空間モデルの概念を簡単に紹介した。今回は、得られたデータに基づいて状態パラメータの「推定」をしたり、将来の観測値の「予測」をしたりする方法を説明する。状態空間モデルの状態の推定について、「フィルタリング」「状態予測」「平滑化」という概念がある。これらは、どの時点までのデータを使ってどの時点の状態を推定するかによる区別である。
フィルタリングは、新しい観測値 $y_t$ が得られたときに、その情報によって現在の状態 $\theta_t$ の推定を更新することをいう。時間順にデータが得られて、その度に現在の状態を更新したいときにフィルタリングが必要とされる。
状態予測は、現在の時点 $t$ までの観測値 $y_{1:t}$ から将来の状態 $\theta_{t+k} (k > 1)$ を予測することである。これは、将来の観測値 $y_{t+k}$ を予測したいときのステップとして用いられる。
平滑化は、現在の時点 $t$ までの観測値 $y_{1:t}$ に基づいて、過去の状態 $\theta_{1:t}$ を推定することを指す。まとまったデータを得たあとで、過去にさかのぼって状態を理解するような場合に用いられる。
この記事では、フィルタリングの説明をしていきたい。また、動的線形モデルの場合のアルゴリズムとして知られている「カルマンフィルタ」もあわせて述べる。状態予測や平滑化は別の記事にしたいと思う。

フィルタリング

フィルタリングは、前の時点の状態の確率分布 $p(\theta_{t-1}|y_{1:t-1})$ がすでにわかっていることを前提として、新しい観測値 $y_t$ を得たときに以下の手順で状態の確率分布 $p(\theta_t|y_{1:t})$ をアップデートする。なお、 $t = 1$ の場合には、 $p(\theta_{t-1}|y_{1:t-1})$ には初期状態 $\theta_0$ の（条件なし）確率分布を用いればよい。これは状態空間モデルを設定する際にすでに定義しているので、観測値を得るごとに $t = 1, 2, ...$ の順で状態を更新していくことができる。
f:id:mstour:20201223205538j:plain

(1) 時点 $t-1$ までの観測値に基づき、状態 $\theta_t$ の予測分布を計算する。

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

右辺の $p(\theta_{t-1}|y_{1:t-1})$ は前の時点までに得られたデータをもとに更新した状態の分布であり、 $p(\theta_t|\theta_{t-1})$ は状態方程式で定義しているので、（理論的には）計算可能である。

(2) 時点 $t-1$ までの観測値に基づき、観測値 $y_t$ の予測分布を計算する。

$\displaystyle p(y_t|y_{1:t-1}) = \int p(y_t|\theta_t) p(\theta_t|y_{1:t-1}) d \theta_t$

右辺の $p(\theta_t|y_{1:t-1})$ は(1)で求めた分布、また $p(y_t|\theta_t)$ は観測方程式で定義されているので、これも原理的に計算可能である。

(3) 新しい観測値 $y_t$ に基づき、状態の更新された分布「フィルタリング分布」を計算する。

$\displaystyle p(\theta_t|y_t) = \frac{ p(y_t|\theta_t) p(\theta_t|y_{1:t-1}) }{ p(y_t|y_{1:t-1}) }$

分子には(1)の予測分布と観測方程式、分母には(2)の予測分布を用いる。

これらの導出では、状態空間モデルの「 $\theta_t$ を与えたもとで $y_t$ と $y_{1:t-1}$ とが独立」という性質が重要な役割を果たしている。

カルマンフィルタ

前述のように一般的な状態空間モデルでのフィルタリングの推定式は具体的に書き下すことができるが、必ずしも計算は簡単ではない。しかしながら、動的線形モデルの場合には、同時分布・周辺分布・条件付き分布全てが正規分布となるため、上記の(1)〜(3)の分布は正規分布で具体的に特定できる。動的線形モデルの場合のフィルタリング解法はカルマンフィルタと呼ばれている。
いま、下記で表される動的線形モデル

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

を考える。前の時点でのフィルタリング分布を $p(\theta_{t-1}|y_{1:t-1}) = N(m_{t-1}, C_{t-1})$ とするとき（ $t=0$ の場合は、 $\theta_0$ の事前分布を用いる）、具体的に以下のように表現できる。

(1) 時点 $t-1$ までの観測値に基づく、状態 $\theta_t$ の予測分布は以下の正規分布になる。

$\displaystyle p(\theta_t|y_{1:t-1}) = N( G_{t-1} m_{t-1}, G_t C_{t-1} G_t^T + W_t )$

(2) 時点 $t-1$ までの観測値に基づく、観測値 $y_t$ の予測分布は以下の正規分布になる。

$\displaystyle p(y_t|y_{1:t-1}) = N( F_t a_t, F_t R_t F_t^T + V_t )$

ただし、 $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$ である。

(3) 新しい観測値 $y_t$ に基づく、状態のフィルタリング分布は以下の正規分布になる。

$\displaystyle p(\theta_t|y_t) = N( a_t + R_t F_t^T Q_t^{-1} (y_t - f_t), R_t - R_t F_t^T Q_t^{-1} F_t R_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$ である。

なおカルマンフィルタを用いる場合にも単純な計算では数値が不安定になることが知られており、様々な方法が提案されている。

まとめ

今回は状態空間モデルの推定についての概要と、推定方式の一つであるフィルタリングについて述べた。次回、将来の時点の予測や、過去の状態に関する平滑化について説明したい。

参考文献

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