【R】時系列分析の覚書(基本、カルマンフィルタ)

前回はウィナーフィルタまでやって力尽きました。

そんなわけで、今回もこちらの本をやっていきます。

今回はカルマンフィルタです。一応位置づけは8章で、基本になっています。

f:id:nogawanogawa:20180715133812j:plain:w500

（これで基本とか、世の中難しすぎだろ…）

カルマンフィルタ

何をやっているのかを見失うと訳わかんなくなるので、目的を明確にします。カルマンフィルタの目的は、

時系列データが与えられたときに、観測値から状態を推定すること

とします。

カルマンフィルタでは、前回同様、線形・ガウス型状態空間モデルを想定します。式としてはこんな感じです。

$x_t = G_t x_{t-1} + w_t$

$y_t = F_t x_t + v_t$

このモデルに従う時系列データで観測値から状態を推定することを考える際にカルマンフィルタを適用することができます。カルマンフィルタリングのアルゴリズムを書くとこんな感じです。

・・・さっぱりわかりません！ってことで、ちょっと他の方のブログを参考にする反則技を使用します。

非常に参考になりました。ありがとうございます。

もとになる制御工学なんかで出てくるカルマンフィルタの図がこちら。

f:id:nogawanogawa:20180916164429j:plain:w600

そんでもって、今回のケースに合うようにいろいろ情報を修正した版がこちら。

f:id:nogawanogawa:20180916164509j:plain:w600

これだとなんかわかる感じがします。ちなみに図の上半分は、モデルとしてはあるけれど実際の点線の中身は観測できません。

観測できないxによって決定される観測値yの値を、カルマンフィルタ側（下半分）で推定してあげます。推定していくと、途中でxの推定値から次のyを推定できるようになります。このyの推定値と実際のyの観測値を比較して推定の誤差を小さくしていきます。

このときに出てくるのがカルマンゲインと言うやつです。カルマンゲインでは、本当のxに対してxの推定値との誤差を最も小さくすることを目標に算出します。すなわちそれは、xの平均二乗誤差を最も小さくなったときに達成され、平均二乗誤差を最小にすることでカルマンゲインが得られます。

ここまで簡単にしてもやっぱりむずいっすね。。。尤度とか分布とか確率モデルを考える前に、普通に実数で考えたほうが良さそうです。

何はともあれなんとなくやってみるのが一番はやいので、やってみます。この本では、サンプルとしてナイル川の流量のデータを使用した状態の推定をやってみます。

理屈は難しかったんですが、実装は式のとおりに書いていくだけなんで、そこは意外と簡単です。

フィルタリングなんで、現在までの観測値を使って現在の状態を推定します。

予測なんで、現在までの観測値を使って未来の状態を推定します。前提として、フィルタリングが完了している（標準とする時刻まで状態の推定が完了している）必要があります。

予測のアルゴリズムを書いておくとこんな感じです。

時点t+(k-1)でのk-1期先予測分布： $a_t(k-1), R_t(k-1)$
時点t+kでの更新手続き
- k期先予測分布
  - (平均) $a_t(k) \gets G_{t+k}a_{k}(k-1)$
  - (分散) $R_t(k) \gets G_{t+k}R_{t}(k-1)G_{t+k}^{\mathrm{T}} + W_{t+k}$
時点tでのフィルタリング分布： $a_{t}(k), R_{t}(k)$

平滑化なんで、現在までの観測値を使用して過去の状態を推定します。こちらもフィルタリングがしてある（標準とする時刻まで状態の推定が完了している）前提になります。

平滑化のアルゴリズムとしてはこんな感じです。

時点t-1でのフィルタリング分布： $s_{t+1}, S_{t+1}$
時点tでの更新手続き
- 平滑化利得 $A_t \gets C_tG_{t+1}^{\mathrm{T}}R_{t+1}^{-1}$
- 状態の更新
  - (平均) $s_t \gets m_t + A_{t} [ s_{t+1} - a_{t+1}$ ]
  - (分散) $S_t \gets C_t + A_{t}[S_{t+1}-R_{t+1}$ ] $A_t^{\mathrm{T}}$
時点tでのフィルタリング分布： $s_{t}, S_{t}$