ill-identified diary

所属組織の見解などとは一切関係なく小難しい話しかしません

トップ > 時系列分析 > [R] [教材]アニメーションで学ぶカルマンフィルタ

[R] [教材]アニメーションで学ぶカルマンフィルタ

時系列分析 統計学 R 教材 視覚化
この記事は最終更新日から3年以上が経過しています

概要

  1. すごく今更感があるが, カルマンフィルタのフィルタリングの話.
  2. アニメーションを作ってみたかっただけともいう.
  3. 簡単な説明なのでもっと具体的な話は他の文献で勉強して欲しい.

カルマンフィルタのアニメーションを作成している記事は既に

があるのだが, カルマンフィルタが逐次に情報を更新し, (長期の) 予測値を理想の値にどう近づけていくかをアニメーションで表そうと思ったので作ってみた.

今回はカルマンフィルタの中で一番基本的な線形・正規分布のモデルで説明する. あまり厳密な解説をするわけではないので, 詳しくは参考文献などを参照されたい.

簡単な説明

時系列モデルというとおそらく一番よく知られているのは ARIMA だろうが, 通常ARIMA は全てのデータを集めてからまとめて処理することで推定, あるいは機械学習の文脈でいうところのバッチ処理, によって推定するのが普通である. これに対してカルマンフィルタは, 1時点のデータが追加されるたびにパラメータを更新することができる*1. つまり機械学習でいうオンライン学習としての側面がある.

第2に, 欠測値のあるデータに対しても適用可能であるという点がある. ARIMA であっても欠測値をなんらかのルールに従いデータを補間することは可能だが, カルマンフィルタはその性質上, 欠測値があってもそのまま計算を続けることができる.

第3に, より多彩な形状のモデルに適用可能である. 詳しくは状態空間モデルについて調べてほしいが, 特筆すべきは, 時間によって値の変わる時変パラメータを含むモデルの推定も可能である, という点である.

以下が50期間ぶんのアニメーション. 赤点が観測点, 緑色が事前推定値で青点が事後推定値, 未来の範囲は灰色で表し, この範囲の予測値を青破線で表した. 事前推定・事後推定の意味はこの後に解説している. この画像の作り方やモデルの設定についても後のセクションを参照されたい. 時系列データの更新のたびに、予測線がそれらしくなっていくのがわかる.

2017/3/14 更新: 95%信頼区間も表示するようにした.

f:id:ill-identified:20170314195830g:plain

2017/3/14 更新: さらにt = 20~29 の10期間を欠測値とした場合の挙動もアニメ化した

f:id:ill-identified:20170314200226g:plain

狭義のカルマンフィルタは線形かつ誤差項が独立正規分布になるようなモデルだけを対象とするが, Unscend カルマンフィルタや拡張カルマンフィルタ, そして粒子フィルタ*2はさらに一般の時系列モデルに対して適応できる*3. カルマンフィルタは状態空間モデルで表現できる時系列モデルであれば適用できる. 今回はシンプルな線形ガウシアン (正規) 状態空間モデルを考える. 状態空間モデルとは一般に, 状態 (state) 方程式または遷移 (transition) 方程式 と呼ばれる差分方程式と, 時系列のラグを伴わない観測方程式で構成される(先に紹介した記事と記号の使い方が違うので注意):

\begin{align} \begin{cases} \boldsymbol{x}_{t+1} & =\mathbf{Z}\boldsymbol{x}_{t}+\mathbf{Q}\boldsymbol{\varepsilon}_{t}\\ \boldsymbol{y}_{t} & =\mathbf{H}\boldsymbol{x}_{t}+\mathbf{R}\boldsymbol{\eta}_{t} \end{cases}\end{align}

この\boldsymbol{\varepsilon}_{t} と \boldsymbol{\eta}_{t}は独立正規分布である. \boldsymbol{x}_{t} は \(t\) 時点の未知のパラメータで, 状態変数または単に状態と呼ばれる. つまり \boldsymbol{x}_{t} は時変パラメータということだが, 時間を通して一定なパラメータであっても可能である. 一方で \boldsymbol{y}_{t} は観測値であり, 実際に観測できなければならない. 係数行列 \mathbf{Z},\mathbf{Q},\mathbf{H},\mathbf{R} も時間によって変化してもよいが, ここでは簡単のため固定している.

数式ではなく直感的なイメージで言うなら, 状態方程式は時間に応じて変化する現象のメカニズム, たとえば株価の推移とか, 運動する物体とか, を表現するのに用いる. そしてこのメカニズムは, 各時点においてどうなっているか(状態) は全てを観測できなかったり, 間接的にしか観察できなかったり, あるいは観測時の誤差が混入していたりする. そこで各時点で実際に観察された値と, メカニズムを関連付けるのが観測方程式である. そのため, 状態方程式は t+1 期の状態と t 期の状態を関連付ける方程式だが, 観測方程式はある1時点 t だけを含む方程式でなければならない.

状態変数のパラメータはカルマンフィルタによって推定できるが, 誤差項の分散は推定できない. また係数行列もモデルの設定によっては自明になるが, 自明でない成分がある場合は, これらを最尤法などを併用して推定することになる.

フィルタリング

カルマンフィルタで逐次的に予測をする操作はフィルタリング*4と呼ばれる.

\mathcal{Y}_{t} を \{\boldsymbol{y}_{1},\cdots,\boldsymbol{y}_{t}\}, t 時点までに得られたデータとする. このとき, \(t-1\) 時点までの情報に基づく \boldsymbol{x}_{t} の期待値を事前状態推定値または1期先予測値という. これを \hat{\boldsymbol{x}}_{t\mid t-1} で表す.

\begin{align} \hat{\boldsymbol{x}}_{t\mid t-1}:= & \mathrm{E}\boldsymbol{x}_{t}\mid\mathcal{Y}_{t-1}\end{align}

さらに, t 時点の情報も含めて求めた \boldsymbol{x}_{t} の期待値を事後状態推定値またはフィルタ化推定値といい, \hat{\boldsymbol{x}}_{t\mid t} と表す.

\begin{align} \hat{\boldsymbol{x}}_{t\mid t}:= & \mathrm{E}\boldsymbol{x}_{t}\mid\mathcal{Y}_{t}\end{align}

また, 観測値 \boldsymbol{y}_{t} と予測値 \hat{\boldsymbol{y}}_{t} との差を\boldsymbol{v}_{t}:=\boldsymbol{y}_{t}-\hat{\boldsymbol{y}}_{t} で表す. これをイノベーションという.

フィルタリングは以下の手順を t=1,2,\cdots に対して繰り返し行うことで, \boldsymbol{x}_{t} の推定値を事前状態推定値から事後状態推定値へと更新しつつ時間 t を進めていく.

  1. t 時点のイノベーション \(\boldsymbol{v}_{t}:=\boldsymbol{y}_{t}-\hat{\boldsymbol{y}}_{t}\) と, \(\boldsymbol{v}_{t}\) の分散値を求める.
  2. 前回 (t-1時点) の予測値 \(\hat{\boldsymbol{x}}_{t\mid t-1}\)を, 新しく得られた \(\boldsymbol{y}_{t}\) の情報を加味して更新し, 事後状態推定値 \(\hat{\boldsymbol{x}}_{t\mid t}\) を計算する. さらにその分散値も求める.
  3. 事後状態推定値 \hat{\boldsymbol{x}}_{t\mid t} と状態方程式から, 次回 \(t+1\) の事前状態予測値 \(\hat{\boldsymbol{x}}_{t+1\mid t}\) を求める. さらにその分散値も求める.

事前推定値と実際の観測値の差であるイノベーションは予測値と実現値の誤差と解釈できる. そしてイノベーションが大きいほど, 事前状態推定値から事後状態推定値へ更新したときの変化の幅は大きくなる. 大雑把なイメージで言うなら, カルマンフィルタとは毎度毎度, 予測の誤差を計算し, 誤差が大きければこれを重く受け止め, 予測値の上方・下方修正を行うアルゴリズム, と捉えることができる

ここでは状態空間モデルの線型性と誤差項の正規分布を仮定しているため, これらの計算は全て行列演算だけで処理できる. 正式な式やなぜこれでうまく行くのかという説明は適当な参考文献, 例えば 足立 and 丸田 (2012) や 野村 (2016) などを.

長期予測

時系列モデルを構築する場合の多くは, 将来の長期にわたる予測を知りたいという要請から来る. カルマンフィルタの予測も1期先だけでなく, もっと未来の予測値をも求めることができる. 具体的には, 上記のフィルタリングの手順と同様に行う. ただし, 1期より先の予測の場合は \boldsymbol{y}_{t}が得られないため, 事前状態推定値 \(\hat{\boldsymbol{x}}_{t\mid t-1}\) から事後状態推定値 \(\hat{\boldsymbol{x}}_{t\mid t}\) への更新時には推定値 (とその分散値) は変化しない. 事後状態推定値 \(\hat{\boldsymbol{x}}_{t\mid t}\) から次期の事前状態推定値 \(\hat{\boldsymbol{x}}_{t+1\mid t}\) への計算は構造パラメータによって決まるので, 長期予測において, より未来のy_t が得られないからといって, 状態推定値を計算しても1期前から値が変わらない, ということにはならない. これは先ほどのアニメーションを見ても明らかである.

欠測値についても, これと同様に観測値が存在しないことから, 事前状態推定から事後状態推定への更新時に値を変えずに持ち越すことで計算できる.

カルマンフィルタの第3の機能として, 平滑化 またはスムージング (smoothing) と呼ばれるものがあるが*5, これは気が向いたら別の機会に説明する.

実演方法

今回は ARMA(1,1) + 線形トレンドのある乱数を生成し, これをカルマンフィルタでフィルタリングと予測をおこなう様子を表した. モデルを数式で表すと以下のようになる.

\begin{align} y_{t}= & \phi y_{t-1}+\theta\varepsilon_{t-1}+\delta t+\varepsilon_{t},\\ \varepsilon_{t}\sim & \mathcal{N}(0,1)\end{align}

状態空間モデルの状態方程式は必ず1階の差分方程式でなければならないので, 移動平均 \theta\varepsilon_{t-1} と 線形トレンド部分 \(\delta t\) をうまく扱う必要がある. 以下のように補助的な状態変数 \(v_{t}:=\theta\varepsilon_{t}\), \(\delta_{t}=\delta_{t-1}+\delta_{2,t}\), \(\delta_{2,t}=\delta_{2,t-1}\) を追加した4次元の状態空間モデルとして扱う必要がある. つまり,

\begin{align} \begin{bmatrix}y_{t}\\ v_{t}\\ \delta_{t}\\ \delta_{2,t} \end{bmatrix}= & \begin{bmatrix}\phi & 1 & 1 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 1 \end{bmatrix}\begin{bmatrix}y_{t-1}\\ v_{t-1}\\ \delta_{t-1}\\ \delta_{2,t-1} \end{bmatrix}+\begin{bmatrix}1 & 0 & 0 & 0\\ \theta & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \end{bmatrix}\begin{bmatrix}\varepsilon_{t}\\ \varepsilon_{t}\\ \varepsilon_{t}\\ \varepsilon_{t} \end{bmatrix},\\ y_{t}= & \begin{bmatrix}1 & 0 & 0 & 0\end{bmatrix}\begin{bmatrix}y_{t}\\ v_{t}\\ \delta_{t}\\ \delta_{2,t} \end{bmatrix}+\mathbf{O}\boldsymbol{\eta}_{t}\end{align}

という状態空間モデルで表現できる. 状態変数ベクトル \boldsymbol{x}_t の第一成分は y_t だが, 計算上は観測値とは別の変数として扱う. さらに, 予測したいのは y_{t} だが, この観測方程式からわかるように, 今回は状態変数ベクトルの第1成分は \(y_{t}\) と一致するため, 状態変数をそのまま予測値として扱うことができる. また, 厳密にはこれは線形トレンドではなく, いわゆるローカル線形トレンドモデルを ARMA に合成したものである. 実際にトレンドの大きさがわかっているのならば毎回の反復計算で定数を足せば良いだけなのだが, 状態空間モデルの標準形の範囲で表現するためにこうした.

プログラムは R を使う. Tusell (2011) によれば, R にはカルマンフィルタを利用できるパッケージとして, dse, dlm, sspir, FKF, KFAS があるが, sspir は作者音信不通により現在は CRAN から削除されている. 機能の多彩さでいうならば dlm と KFAS に長所があるが, 実行速度に関しては dlm は比較的遅く, KFAS は低次元のモデルに限って言えば速いが, 次元が増えると遅くなりやすいようだ. FKF は Fast Kalman Filter の略だけあって特に速いが, フィルタリング (とプロット) しかできないという無骨な仕様である. つまり, FKF を使う場合は先に述べたように状態変数いがいのパラメータを推定するプログラムを別途用意しなければならない.

今回はアニメーションの作成にフィルタリング時の事前・事後の状態推定値, さらに長期予測についても同様に2つの状態推定値が必要となる. しかし先に挙げたパッケージは事後推定値だけを出力するものが多いので, これらのプログラムを1から書いた. このプログラムはカルマンフィルタの挙動をわかりやすくするためのものなので, 計算の精度や速度は考慮していないので実用的ではない. 実用上は dlm か KFAS あたりを使うことが多いと思う. これらの使い方は他のサイトや文献を参照してほしい(なお, 野村 (2016) は KFAS パッケージの使い方にも言及がある.).

以下はプログラム


giste49d2bc38d5a0e6ad6f1c88a58250ddd


参考文献

カルマンフィルタの入門教科書になりそうなものを2種類挙げる. 2冊はそれぞれ切り口が違う. 足立・丸田本は制御理論の文脈で書かれており, 野村本は状態空間モデルと期待値 (と分散) の計算を意識して式を説明している. 今回の記事の範囲をより詳しく知りたいのならどちらでも良いので, 自分の持つ知識と照らしあわせて近い方を選ぶと良いと思う. 足立・丸田本は matlab のサンプルプログラムを掲載し, 野村本は R の KFAS パッケージを利用したサンプルプログラムを掲載している. また, 野村本には粒子フィルタの説明もあるが, 足立・丸田本にはない.

*1:実際には ARIMA に対しても逐次最小二乗法というものがあるが, 実質的に線形カルマンフィルタと同じである.

*2:逐次モンテカルロ法とも呼ばれる.

*3:他にもアンサンブルカルマンフィルタとか適応フィルタとか, 分野によっていろいろバリエーションがある.

*4:分野によっては「濾波」とも.

*5:厳密にはフィルタリングとスムージングは異なる概念なので, カルマンフィルタではなくカルマンスムーザと呼ぶべきだと主張する人もいるが, 平滑化の機能も含めて一連のテクニックをカルマンフィルタと呼ぶことが多い.

ill-identified 読者になる

広告を非表示にする
コメントを書く