はじめに
こんにちは。普段、時系列データの逐次予測を行うときは、状態空間モデルを構築して粒子フィルタを使っていますが、そもそも一番簡単な線形ガウスモデルの場合に有効なカルマンフィルタをしっかりと勉強したことがなかったので、今回まとめてみました。（本記事では、計算方法の確認程度にとどめていますので、より詳しく知りたい方は、他の良い記事や良書を適宜参照してください。）
状態空間モデル
時系列データをモデリングする枠組みとして、状態空間モデルという考え方があります。これは、観測値が観測される過程を、対象が時間変化する過程と観測データが観測される過程に分けてモデリングします。状態空間モデルの説明には、下図のようなグラフィカルモデルによる表現が用いられます。
$x_t$は観測対象の状態量、$y_t$は観測値を表しています。観測値が観測される裏側で観測対象が時間変化していく様子が表されています。これを確率分布を用いて表すと、以下のようになります。
\begin{align}
x_t &\sim p(x_t | x_{t-1}) \\
y_t &\sim p(y_t | x_t)
\end{align}
これらはそれぞれ、状態方程式（またはシステム方程式）と観測方程式と呼ばれます。ここで重要なのは、図中で示されている通り、状態ベクトルは一時刻前の状態ベクトルに依存するということで、これをマルコフ性と言います。
さて、時系列データを観測したときに我々がやりたいこととして、以下のようなことが挙げられます。
- 観測値$y_t$が観測されたのちに、次の時刻の状態$x_{t+1}$または観測値$y_{t+1}$を予測する
- 観測値$y_t$が観測された時の観測対象の状態$x_t$を推定する
これらは、それぞれ一期先予測やフィルタリングなどと呼ばれています。状態空間モデルを式（１）、（２）のように表した場合、一期先予測やフィルタリングは以下のように表現できます。
まず一期先予測ですが、時刻$t-1$のフィルタ分布$p(x_{t-1}|y_{t-1})$が与えられているとすると、以下のようにして一期先状態予測が計算できます（2行目から3行目は状態空間モデルのマルコフ性より）。
\begin{align}\begin{split}
p(x_t | y_{1:t-1}) &= \int{p(x_t, x_{t-1} | y_{1:t-1})}dx_{t-1}\\
&=\int{p(x_t | x_{t-1}, y_{1:t-1})p(x_{t-1} | y_{1:t-1})}dx_{t-1}\\
&=\int{p(x_t | x_{t-1})p(x_{t-1}|y_{1:t-1})}dx_{t-1}
\end{split}\end{align}
次に、フィルタ分布$p(x_t|y_t)$は、時刻$t$の観測値$y_t$が得られた後、以下のようにして求めます。
\begin{align}\begin{split}
p(x_t|y_{1:t}) &= \frac{p(x_t, y_t|y_{1:t-1})}{p(y_t|y_{1:t-1})} \\
&= \frac{p(y_t|x_t, y_{1:t-1}) p(x_t|y_{1:t-1})}{p(y_t|y_{1:t-1})} \\
&= \frac{p(y_t|x_t) p(x_t|y_{1:t-1})}{p(y_t|y_{1:t-1})} \\
&= \frac{p(y_t|x_t) p(x_t|y_{1:t-1})}{\int p(y_t, x_t|y_{1:t-1}) dx_t} \\
&= \frac{p(y_t|x_t) p(x_t|y_{1:t-1})}{\int p(y_t|x_t) p(x_t|y_{1:t-1}) dx_t}
\end{split}\end{align}
このようにして、一期先予測とフィルタを求めることができますが、式（３）、（４）には積分計算が含まれており、これでは容易に実装することができません。
線形ガウス型状態空間モデル
一期先予測とフィルタ推定の計算方法を確認しましたが、実は状態空間モデルが線形ガウス型である場合には、簡単な行列計算で求めることができます。線形ガウス型の状態空間モデルは、以下のように表されます。
\begin{align}
x_t &= F_t x_{t-1} + G_t v_t \\
y_t &= H_t x_t + w_t
\end{align}
ここで状態ベクトル$x_n$の次元を$p$、$F_t$は$p \times p$の行列、$G_t$は$p \times m$の行列、$v_t$は$m$次元の正規分布$N(0, Q)$に従う乱数、$H_t$は$q \times p$の行列、$w_t$は$q$次元の正規分布$N(0, R)$に従う乱数とします。つまり、全て線形の計算であり、なおかつ状態変動のノイズや観測ノイズがガウス分布であることを表しています。時刻$t=0$の状態ベクトルの平均と分散共分散行列をそれぞれ$x_{0|0}$、$V_{0|0}$とすると、一期先予測とフィルタ推定が以下のように計算できます（※それなりに計算が大変なのでここでは省略します）。
（1期先予測）
\begin{align}
x_{t|t-1} &= F_t x_{t-1|t-1} \\
V_{t|t-1} &= F_t V_{t-1|t-1} F_t^T + G_tQ_tG_t^T
\end{align}
（フィルタ）
\begin{align}
K_t &= V_{t|t-1} H_t^T (H_t V_{t|t-1} H_t^T + R_t)^{-1} \\
x_{t|t} &= x_{t|t-1} + K_t (y_t - H_t x_{t|t-1}) \\
V_{t|t} &= V_{t|t-1} - K_t H_t V_{t|t-1}
\end{align}
実行例
それでは、カルマンフィルタを使って状態推定する様子の簡単な例を示します。対象のデータは、ランダムウォークする状態$x_t$に標準偏差0.3の正規分布に従う観測ノイズが加わったローカルレベルモデルのデータを用意しました。今回使用した観測データを以下に示します。
この時系列を時刻$t=0$から逐次的に読み込んで状態の推定、一期先予測を繰り返しています。
using Random
using Plots
using LinearAlgebra
Random.seed!(123)
N = 1000
steps = randn(N)
random_walk = cumsum(steps)
noise = 0.3 .* randn(N)
observed_walk = random_walk .+ noise
F_t = 1.0
H_t = 1.0
Q_t = 1.0
R_t = 0.3^2
x_t_t_1 = zeros(1, N)
V_t_t_1 = zeros(1, 1, N)
K_t = zeros(1, 1, N)
x_t_t_1[1] = observed_walk[1]
V_t_t_1[1] = 1.0
for t in 2:N
# 一期先予測
x_t_t_1[t] = F_t * x_t_t_1[t-1]
V_t_t_1[:, :, t] = F_t * V_t_t_1[:, :, t-1] * F_t' + Q_t
# フィルタ
K_t[:, :, t] = V_t_t_1[:, :, t] * H_t' * inv(H_t * V_t_t_1[:, :, t] * H_t' + R_t)
x_t_t_1[t] = x_t_t_1[t] + K_t[:, :, t] * (observed_walk[t] - H_t * x_t_t_1[t])
V_t_t_1[:, :, t] = (1 - K_t[:, :, t] * H_t) * V_t_t_1[:, :, t]
end
状態を推定した結果を以下に示します。
黒色が真の状態で、赤色線がカルマンフィルタを使って推定したフィルタ推定量です。観測値にはノイズが加わっているにも関わらず、真の状態に近い推定量が得られていることがわかります。
扱わなかった内容
以下の内容は実用上は重要ですが、今回は取り上げませんでした。また別の記事で詳しく説明する予定です。
- 状態ベクトルの初期（$x_{0|0}$や$V_{0|0}$）の設定方法
- パラメータの決定方法
- 拡大状態ベクトル
- 非線形モデルの推定（アンサンブルカルマンフィルタなど）
参考文献
- データ同化入門、樋口 知之(編著)ほか、朝倉書店
- カルマンフィルタ、野村俊一、共立出版
- 経済・ファイナンスのためのカルマンフィルター入門、森平爽一郎、朝倉書店