次回列分析においては、ARIMAやSARIMAのように、過去の観測値をもとにモデル化する手法がある一方で、状態空間モデルは、観測値以外の見えない状態量があると仮定し、その状態により観測値が生起されるとしてモデル化する手法です。
この状態空間モデルにおいて、仮定する潜在変数がガウス分布に従うという仮定した（これは結構強い仮定である）システムを線形動的システム(Linear Dynamical System) という。このように、状態変数がガウス分布に従うと呼ばない仮定をおかない、非線形動的システム(NonLinear Dynamical System)も存在する。（拡張カルマンフィルタ などがこのモデルである）
ローカルモデルは、線形動的システムの1種である。今回は、ローカルモデルを用いて時系列を予測する手法について、数学的な背景についてわかりやすく解説した後で、実際のデータセットを用いて時系列予測と異常検知を行う。
連続状態空間モデルの数学的な背景
まず、ローカルレベルモデルを扱う前に、一般化された連続空間モデルの概要について解説する。連続状態空間モデルでは、見えない状態変数(隠れ変数)$\bm{z}$を導入します。この状態変数$\bm{z}$によって、観測値として、実際の時系列データである$y$が得られるということである。
グラフィカルモデルに表すと下記のようになる。
これを、もう少し定式化してみましょう。一般的に、状態空間モデルは、状態方程式と観測方程式でシステムを表現できます。
\begin{equation}
\begin{split}
\bm{z_t} = \bm{F_t}(\bm{z_{t-1}}) + \bm{G_t}(\bm{v_t})
\end{split}
\end{equation}
\begin{equation}
\begin{split}
\bm{y_t} = \bm{H_t}(\bm{z_{t}}) + \bm{w_t}
\end{split}
\end{equation}(1)式が状態方程式で、(2)式が観測方程式である。制御や機械学習など学問分野によって(1)や(2)式の呼び方はまちまちである。(1)はシステムモデル、(2)は観測モデル 等の呼ばれ方をすることもある。
大量の変数が出てしまうが、これらは1つ1つ丁寧に理解する必要性がある。
まず、$\bm{zt}$ であるが、これが時刻$t$における観測できない状態量である。状態変量はモデルの設計者が決められるが、多変量の状態を持つことができ、$k$次元のデータで表現することができる。続いて、$\bm{v_t}$は、状態ノイズと呼ばれている。上のグラフィカルモデルを見るとわかるが、内部の状態が時刻によって$z_{a} → z_{a+1}$ に変化するが、そのノイズである。
連続空間モデルの場合、このノイズがどのような分布に従うかは設計者が決められるが、今回は簡単のため、このノイズもガウス分布に従うと仮定する。
\begin{equation}
\begin{split}
\bm{v_t} \sim \mathcal{N} (0, \bm{Q_t})
\end{split}
\end{equation}すると、(3)式のように表現することができる。$\mathcal{N}$は、ガウス分布に従うことを示している。今回は、この$\bm{v_t}$はノイズと考えているため、平均は0である。またこのノイズの分散共分散行列を$\bm{Q_t}$としている。 このノイズの次元は任意であり、$m$次元である。
続いて(1)に登場する$\bm{F_t}$であるが、これは状態遷移行列と呼ばれる、状態量$\bm{z_{t-1}}$を$\bm{z_{t}}$に変更するための係数行列である。
この$\bm{F_t}$は予測モデルにおけるパラメータなのだが、観測できない状態を変更するためのパラメータである。例えば、重力や空気抵抗下での飛行機の位置を予測するモデルにおいては、この状態遷移行列$\bm{F_t}$は重力や空気抵抗、摩擦といった条件がパラメタとして反映されることになる。
続いて、(1)式の$\bm{G_t}$であるが、これは状態空間のおけるノイズの遷移ベクトルであり、$\bm{G_t} \in \mathcal{R}^{k \times m}$のベクトルである。状況によって、システムノイズの分布も刻一刻と変更することをモデル化しているための係数である。
続いて、(2)式の変数の説明に入る。
(2)式に登場する、\bm{y_t}が、実際に観測されるデータである。単変量でも多変量でも考えることができるか、多変量で考える場合は、その次元を$l$とする。
続いて、$\bm{H_t}$は、状態変数$\bm{z}$から$\bm{y}$に変換する行列である。状態空間から観測できる空間に写像するような行列となっており、観測行列(observation matrix)と呼ばれることもある。今回、観測される$\bm{y}$がl次元の場合は、$ l \times k$の行列になる。
続いて、$\bm{w_t}$は、観測時に発生するノイズであり、このノイズも$\bm{v_t}$と同じようにガウス分布に従うと仮定すると、平均は0、分散共分散行列を$\bm{R_t}$とすることができる。
線形ガウス型モデルのモデル
改めて、状態の変化が線形変化であり、ノイズもガウス分布に従うという線形ガウス型モデルを定式化する。
\begin{equation}
\begin{split}
\bm{z_t} = \bm{F_t}(\bm{z_{t-1}}) + \bm{G_t}\bm{v_t}, \quad \bm{v_t} \sim \mathcal{N}(0, \bm{Q_t}) \\
\bm{y_t} = \bm{H_t}\bm{z_t} + \bm{w_t}, \quad \bm{w_t} \sim \mathcal{N}(0, \bm{R_t})
\end{split}
\end{equation}ここで、(4)式の上側が状態方程式で、下側が観測方程式である。
線形ガウス型モデルの状態推定と予測
この時系列モデルにおいて、主要な問題は時系列が与えられたときに、その時の状態や将来の状態を予想することである。もし与えられた時系列データから任意の時点の状態を予想することができれば、(4)式の観測方程式で実際に観測できる値は計算することができるからである。
実際、時刻$t$までの観測値が与えられていた時、時刻$t$より前の状態を推定することを、平滑化(smoothing、スムージング)、現時刻$t$の状態を推定することを、フィルタリング、現時刻tより将来のことを予測することを、予測(prediction)と呼ぶ。
今回は時系列データの異常検知が題材なので、過去の状態を推定する平滑化は扱わず、こちらに関しては別の記事で扱うことにする。この記事では、現在の状態を推定するフィルタリングと予測について扱ってみる。
将来予測の定式化
今、時刻$t-1$までの時系列データ$\bm{Y_{1:t-1}}$が得られている時、時刻$t$の状態$\bm{x_t}$を予測するのは、
\begin{equation}
\begin{split}
p(\bm{x}_t | \bm{Y}_{1:t-1}) &= \int_{- \infin}^{\infin}
p(\bm{x}_t , \bm{x}_{t-1} | \bm{Y}_{1:t-1})
d \bm{x}_{t-1} \\
&=
\int_{- \infin}^{\infin}
p(\bm{x}_t | \bm{x}_{t-1}, \bm{Y}_{1:t-1}) p( \bm{x}_{t-1} | \bm{Y}_{1:t-1})
d \bm{x}_{t-1} \\
&=
\int_{- \infin}^{\infin}
p(\bm{x}_t | \bm{x}_{t-1}) p( \bm{x}_{t-1} | \bm{Y}_{1:t-1})
d \bm{x}_{t-1} \\
\end{split}
\end{equation}
という定式化をすることができる。これが予測分布の式である。(5)式は、状態空間モデルやカルマンフィルタ の参考書で至る所で登場するが、このような式変形にひよってはいけない。この式変形においては、単純に条件付き確率分布の連鎖律と、マルコフ連鎖による仮定しか用いていない。また、実務においても(5)式は覚える必要ない。一度式変形を理解して自分で納得すれば、実装もできるし、同様の手法を用いた論文も読めるようになる。
(5)式の解説の前に、現在の状態を推定するフィルタリングについても定式化する。
\begin{equation}
\begin{split}
p(\bm{x}_t | \bm{Y}_{1:t}) &=
p(\bm{x}_t | \bm{y}_t, \bm{Y}_{1:t-1}) \\
&=
\frac{p(\bm{x}_t, \bm{y}_t| \bm{Y}_{1:t-1})}{p(\bm{y}_t| \bm{Y}_{1:t-1})} \\
&=
\frac{p(\bm{y}_t| \bm{x}_t, \bm{Y}_{1:t-1}) p(\bm{x}_t | \bm{Y}_{1:t-1}) }{p(\bm{y}_t| \bm{Y}_{1:t-1})} \\
&=
\frac{p(\bm{y}_t| \bm{x}_t) p(\bm{x}_t | \bm{Y}_{1:t-1}) }{p(\bm{y}_t| \bm{Y}_{1:t-1})} \\
\end{split}
\end{equation}
(5)式、および、(6)式によって、将来の状態予測、および現在の状態予測の式が得られる。
ここまでで、線形ガウス型モデルにおける、山の30%を登ったところである。次は、この(5)式と(6)式を観察してみよう。
カルマンフィルタ のアルゴリズムの導入
ここで、カルマンフィルタ のアルゴリズムを考える。
カルマンフィルタ では、観測値と予測値から、逐次的に状態の予測とフィルタリングを続けていくが、この時、$t-1$までの時系列データ $\bm{Y}_{t-1}$が得られた時の、状態の予測分布 $p(\bm{x}_{t} | \bm{Y}_{t-1})$ の平均を$bm{x}_{t|t-1}$、分散共分散行列を$V_{t|t-1}$とすると、下記のような更新をすることができる。
1期先予測分布の平均$bm{x}_{t|t-1}$ と分散共分散行列の更新式 $V_{t|t-1}$
\begin{equation}
\begin{split}
\bm{x}_{t|t-1} &= \bm{F}_t \space \bm{x}_{t-1|t-1} \\
\bm{V}_{t|t-1} &= \bm{F}_t \space \bm{V}_{t-1|t-1} \bm{F'}_t + \bm{G}_t \bm{Q}_t \bm{G'}_t
\end{split}
\end{equation}Pythonで線形ガウス型モデルの実装
今、予測する値を$\bm{y_{t+1}}$としましょう。要は現在の時刻が$t$のとき、1秒後の時刻の値、$t+1$秒後の値、$\bm{y_{t+1}}$を予測する問題を考えます。ここで、現在の時刻$t$までのデータ$\bm{Y_{1:t}}$は既知のものであるとします。
要は、現時刻までのデータは存在している状況の中で、将来の予測をしていこうと考えているということです。
状態空間モデルのローカルモデルの定式化
それでは、まずローカルモデルによる定式化をしていきましょう。ローカルモデルでは、見えない状態変量(隠れ変数)$z$を導入します。
グラフィカルモデルのように考えるのであれば、上の図のようにあります。観測値は$y$ですが、この観測量は状態量$\bm{z}$ によって決まるというようなモデルを考えます。
その前に、解きたい問題というのは下記のようになります。
\begin{equation}
\begin{split}
p(y_{t+1} | \bm{Y}_{1:t})
\end{split}
\end{equation}(1)式の意味としては、時刻1~tまでの全データ$\bm{Y}_{1:t}$を利用して、$y_{t+1}$を予測するのが問題です。これが、1期先の時系列の予測分布となります。
しかし、実際には、上のグラフィカルモデルのように、見えない潜在変数$\bm{z}$を仮定してモデル化をします。
ここで、求めたい(1)式の分布がガウス分布になることを仮定すると、(1)式は平均ベクトルと共分散行列を用いて、下記のように表現することができます。
\begin{equation}
\begin{split}
p(y_{t+1} | \bm{Y}_{1:t}) = \mathcal{N} (y_{t+1} | \bm{\mu_t}, \bm{\Sigma_t} )
\end{split}
\end{equation}予測分布が正規分布になると仮定すると、時刻$t$における、ガウス分布の2つのパラメタ$\bm{\mu_t}, \bm{\Sigma_t}$を求める問題に帰着したことがわかりました。