Re:ゼロから始めるML生活

どちらかといえばエミリア派です

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

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

tsunotsuno.hatenablog.com

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

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

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

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

カルマンフィルタ

カルマンフィルタとは

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

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

とします。

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


x_t = G_t x_{t-1} + w_t


y_t = F_t x_t + v_t

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

  1. 時点t-1でのフィルタリング分布:m_{t-1}, C_{t-1}
  2. 時点tでの更新手続き
    • 一期先予測分布
      • (平均) a_t \gets G_tm_{t-1}
      • (分散) R_t \gets G_tC_{t-1}G_{t}^{\mathrm{T}} + W_{t}
    • 一期先予測尤度
      • (平均) f_t \gets F_ta_{t}
      • (分散) Q_t \gets F_tR_{t}F_{t}^{\mathrm{T}} + V_{t}
    • カルマン利得 K_t \gets R_tF_{t}^{\mathrm{T}}Q_t^{-1}
    • 状態の更新
      • (平均) m_t \gets a_t + K_{t} [ y_t - f_t ]
      • (分散) C_t \gets [I-K_tF_t ]  R_t
  3. 時点tでのフィルタリング分布:m_{t}, C_{t}

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

qiita.com

カルマンフィルター - ORWiki

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

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

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

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

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

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

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

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

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

カルマンフィルタを使ってみる

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

  • 観測値 : y(Nile)
  • 状態 : m

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

フィルタリング

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

予測

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

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

  1. 時点t+(k-1)でのk-1期先予測分布:a_t(k-1), R_t(k-1)
  2. 時点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}
  3. 時点tでのフィルタリング分布:a_{t}(k), R_{t}(k)

平滑化

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

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

  1. 時点t-1でのフィルタリング分布:s_{t+1}, S_{t+1}
  2. 時点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}}
  3. 時点tでのフィルタリング分布:s_{t}, S_{t}

ローカルレベルモデルでやってみる

なんとな~く状態の推定ができそうだと言うことがわかったところで、今度はちゃんと分析してみます。

対象データ

観測値はナイルを使用していました。図だけ出すとこんな感じです。 データに関する事前の準備については、以前もやったのですがもう一度やってみます。

なんとなく見ていくと、5数要約としてはこんな感じですね。

  • 最小値:456.0
  • 中央値:893.5
  • 平均値:919.4
  • 最大値:1370.0

感覚的には外れ値もなさそうです。 自己相関(acf)としては特定の相関は見られない気がします。 スペクトルを見てみても、相関が見られないような感じですね。

ローカルレベルモデルの定義


x_t = x_{t-1} + w_t


y_t = x_t + v_t

いつもとほとんど同じ式ですね。 GとFが常に1なだけです。 Fが1なので、状態と観測値の間にホワイトノイズしか乗ってきていないことになります。 つまり、ノイズの分の誤差はあれど、それ以外は状態と観測値は一致することを意味しています。

分析

さっきとほとんど同じなんですが、一応分析の全体像を示します。

f:id:nogawanogawa:20180917223221j:plain:w500

上では、著者お手製のカルマンフィルタを使って推定しましたが、今回はRのdlmライブラリを駆使して推定するため、 コードの流れが変わりますが、全体像としてやることは同じです。 ライブラリを使うためのお作法による違いがあるだけです。

Nileのデータがローカルレベルモデルに従うと仮定して、モデルを構築します。 その後、フィルタリング・予測・平滑化によって状態を推定します。

モデルの構築

何はともあれモデルを作ってみます。 この辺はもう書いてあるとおりですね。。。

フィルタリング

ローカルレベルモデルになると若干モデルの支配方程式が変わるので、結果も自ずと変わってきます。 まずはフィルタリングをやってみます。

赤い実線がフィルタリングによる状態の推定の平均、緑の実践が95%信頼区間の境界を示しています。 代替がNileのグラフを覆っているので、そこそこの制度で見積もれているように見受けられます。

予測

次に予測です。

これは、未来の状態を求めているので、まあこうなりますよね。 上下の破線は95%信頼区間を表しています。 時間の経過とともに正確に推定できず、範囲が広がっていく様子がわかります。

平滑化

最後に平滑化です。

一番最後に一つおまけでグラフを付けてみました。 95%信頼区間の推移を色別に示しています。

フィルタリングに比べて平滑化のほうが気持ち範囲が狭まっていることがわかりますでしょうか? 平滑化によって、より詳細に状態が推定されていることがわかります。

結果の吟味

さて、ちゃんと出てきた結果のについて評価します。

観測値yに対して、仮定したモデルが適合しているかどうかを判断することになります。 実際には、様々なモデルや前提に対して状態を推定してそれぞれを相対的に比較することで、現象がどのモデルに従っているかを判断します。

今回はローカルレベルモデルしかやっていないので、基準となる指標を出してみます。

まず尤度について出しています。

次にイノベーション(予測誤差・残差)を出しています。 ここで言うイノベーションeとは観測値yと一期先予測尤度fに対して


e_t = y_t + f_t

で表されます。

これを一期先予測尤度の分散で割ったe_t / \sqrt{Q_t}を規格化イノベーションと言います。

適切にモデルが設定されていれば、規格化イノベーションはノイズの影響による正規分布に従うはずです。 今回はそちらも確認します。

最後に、QQプロットを出しています。 規格化イノベーションeのプロットがy=xの式に近づけば近づくほど正規分布に従っていることを表します。 まあ大体正規分布かなというレベルですね。

そんなわけで、参考値程度ではありますがなんとなく適切にモデルを適用できているように見受けられます。

感想

いやー、難しいです。 次はいろいろ分析をやってみます。