データ分析関連のまとめ

データ分析・機械学習周りのもくもく会LTやイベント参加をまとめていきます

状態空間モデルの基礎まとめ②(カルマンフィルタ)

技術キャッチアップ

状態空間モデルの概要を知り、使えるようにするため基礎からまとめていきます。
今回は線形ガウス空間モデルからカルマンフィルタの性質を調べる所までです
(間違っている点等がある場合、指摘いただければ幸いです)
※前回進めたのはここまで

yhiss.hatenablog.com

今回の目的

  • 下記を理解する事を目的としています
    • カルマンフィルタはそもそも何をするのか？
    • どのような仮定をおいているか？

線形・ガウス状態空間モデル

線形・ガウス状態空間モデルは下記のように表される。

{\omega}_t = F_t{\omega}_{t-1} + G_t u_t
y_t = H_t{\omega}_t + v_t
{\omega_t}:K次元ベクトル、u_t:L次元ベクトル、y_t・v_t:M次元ベクトル、F_t・G_t・H_T:行列

  • 上の式が状態の時間を表すシステムモデル、下の式が観測と予測の差を表す観測モデル

  • 上記の行列は現象に関係する法則や統計等に基づく関係式・仮定に基づく既知の行列となっている。
    今回は係数が時間変化しないとして以下のように簡単な形にする。

{\omega}_t = {\omega}_{t-1} + u_t
y_t = {x_t}^T{\omega}_t + v_t

導出の前に、カルマンフィルタでは以下の仮定を置いている

・u_t \sim N(0,U_t):平均0、分散U_tの正規分布
・v_t \sim N(0,D_t):平均0、分散D_tの正規分布
・u_tと{\omega}_t:独立 p(u_t|{\omega}_t)=p(u_t)
・v_tと{\omega}_t:独立 p(v_t|{\omega}_t)=p(v_t)
・マルコフ性
p({\omega}_t|{\omega}_{t-1},y_{1:t}) = p({\omega}_t|{\omega}_{t-1})
p(y_t|{\omega}_{1:t},y_{1:t-1}) = p(y_t|{\omega}_t)
x_{1:t}=x_1,x_2,...,x_t
・時刻t-1の\omegaのフィルタ分布は正規分布

まずシステムモデル、観測モデルにおけるノイズが正規分布、そして状態を表す変数ωはそれぞれのノイズに対して独立だと仮定されている。

カルマンフィルタにおける一期先予測

今回、t-1時点までの情報を用いてt時点での一期先予測を求めるためまずは以下の式

p({\omega}_t|{\omega}_{t-1})={\int}p({\omega}_t|{\omega}_{t-1},u_t)p(u_t)du_t
u_tが正規分布に従う仮定から
p({\omega}_t|{\omega}_{t-1})=(2{\pi})^{-\frac{K}{2}}|U_t|^{-\frac{1}{2}}exp(-\frac{1}{2}({\omega}_t-{\omega}_{t-1})^TU_t^{-1}({\omega}_t-{\omega}_{t-1}))
→{\omega}_t|{\omega}_{t-1} \sim N({\omega}_{t-1},U_t)

上記のように条件付確率はt-1時点での状態変数を平均とする正規分布になる(分散は変化なし)

そしてt-1時点でのフィルタの分布を下記のように仮定

{\omega}_{t-1}|y_{1:t-1} \sim N({\omega}_{t-1|t-1},Q_{t-1|t-1})
y_t = {x_t}^T{\omega}_t + v_t

上記の仮定から、一期先予測の分布は

平均E({\omega}_{t-1}|y_{1:t-1}) = {\omega}_{t|t-1}
分散共分散行列V({\omega}_{t-1}|y_{1:t-1}) = Q_{t|t-1} = Q_{t-1|t-1} +U_t = E(({\omega}_t - {\omega}_{t|t-1})({\omega}_t - {\omega}_{t|t-1})^T|y_{1:t-1})

一期先予測を行う際にノイズの分だけ分散が大きくなっていることがわかる。

カルマンフィルタ

先ほどの仮定からフィルタ分布は以下の式となる

{\omega}_t|y_{1:t} \sim N({\omega}_{t|t},Q_{t|t}):正規分布
{\omega}_{t|t} = {\omega}_{t|t-1}+K_tv_t
v_t = y_t - x_t^T{\omega}_{t|t-1}:イノベーション(予測値と観測値の差)
Q_{t|t} = Q_{t|t-1} - K_t x_t^T Q_{t|t-1}
K_t = Q_{t|t-1}x_t(x_t^T Q_{t|t-1} x_t + D_t)^{-1}
 K_t:カルマンゲイン(状態変数がt→t-1に更新されるときの変化量)
x_t^T Q_{t|t} x_t \leq x_t^T Q_{t|t-1} x_t
  • 時刻t-1のフィルタ分布:正規分布→時刻t-1の一期先予測分布:正規分布→時刻tのフィルタ分布:フィルタ分布:正規分布と表されている
    →或る時刻のフィルタ分布:正規分布ならそれ以降の時刻において一期先予測とフィルタの分布が正規分布となる
  • また、一期先予測では大きくなっていた分散だがフィルタ後に小さくなっている。

カルマンフィルタにおける流れ

  • カルマンフィルタの関係式は上述の一期先予測分布とフィルタ分布となる
  • また状態の推定は下記のような流れとなっている
    • ある時刻(仮にt-1)までの情報を基に一期先予測を行う(誤差分散増える)
    • 実際に値を観測する
    • 一期先予測の結果と観測値をつかって時刻tでの状態をフィルタで推定する(誤差分散減る)

留意事項(パラメータ調整等)

  • カルマンフィルタには仮定や係数の性質から下記のような留意事項が存在する
  • システムモデルでは難しいが、観測モデルにおいて下記を確認する。
    • 目的変数と説明変数が線形の関係にある必要がある
    • ノイズが正規分布である
  • パラメータの設定について
    • 観測ノイズの分散とシステムノイズの分散は係数の変動を決めるので重要なパラメータとなる
      • 観測ノイズ:過去のデータから算出した予測の平均二乗誤差から見積もる。予測が難しい場合には大きくするのが望ましい。
      • システムノイズ:対象の時刻とその一期前の係数の差の二乗和から見積もる
      • 最尤法に基づく手法やEMアルゴリズムに基づく方法もある
  • カルマンフィルタによる係数更新の適切さの確認
  • イノベーションの分布を時系列でプロットし、ばらつきが大きくなったりしていないか確認する

まとめ

  • カルマンフィルタ:一期先予測とフィルタを使って状態を推定する係数を逐次学習させている
  • ノイズに対して正規分布の仮定をおき、或る時点においてフィルタの分布が正規分布ならばそれ以降のフィルタ分布は全て正規分布となる
  • ノイズが正規分布であることや説明変数と目的変数が線形であること等性質を考慮したモデリングが必要となる

Next Action

  • 以下を行っていく予定です
    • カルマンフィルタを用いた実際のモデリング
    • 非線形なデータにも対応できる状態空間モデルの習得

参考文献

谷崎 久志:非線形・非正規状態空間モデルの推定について
http://www.lib.kobe-u.ac.jp/repository/00056064.pdf
気象庁予報部:ガイダンスの解説
https://www.jma.go.jp/jma/kishou/books/nwpreport/64/No64_all.pdf
佐藤整尚:状態空間モデルを用いた金融時系列分析
https://www.carf.e.u-tokyo.ac.jp/old/research/archive/-2014aug/20110926sato.pdf