畳庵〜tatamiya practice〜

機械学習・統計モデリング練習帳。読書・source code読みの記録や，実装など。

R経由で学ぶPythonの状態空間時系列モデリング#1: 逐次推定法とカルマンフィルタ

時系列分析 状態空間モデル

はじめに

概要

この記事は、状態空間時系列モデリングを解説する全5本の記事の第2弾です。

前回の記事では、状態空間モデルの概要について紹介しました。

状態空間モデルでは、トレンドや季節性といった時系列の特徴をマルコフ過程に従う「状態」として仮定し、観測データをもとに状態の予測・推定を行います。

tatamiya-practice.hatenablog.com

今回は状態推定方法のうち逐次推定法と呼ばれるものの概要から入り、カルマンフィルタについて解説を行います。

状態推定・予測の方法には、大きく分けて逐次推定法とMarkov Chain Monte Carlo法(MCMC)の2種類があります。

このうち前者の逐次推定法については、モデルが線型かつシステムノイズ・観測ノイズがガウス分布に従う場合、カルマンフィルタといって近似やサンプリングに頼らず厳密に推定値・予測値の分布を求めることができます。

線型ガウスのモデルは状態空間時系列モデリングでも基本となり応用範囲も広いため、 カルマンフィルタはKFASやstatsmodelsといった多くのライブラリで実装されています。

そこで今回は逐次推定法とカルマンフィルタの仕組みを理解することで、次回以降ライブラリを用いてモデリングを行う際に、出力された値がどのような過程を経て算出されたか、処理の裏でどのような演算が行われているかイメージできるようになることを目指します。

f:id:tatamiyatamiyatatatamiya:20210110150254p:plain:w500

当初の予定では今回の記事でKFASを用いたRによる実行方法も解説する予定でしたが、理論面で解説する内容が膨らんだため、次回に回すことにしました。

なお、カルマンフィルタの漸化式などの詳細な導出は今回は行いません。

忙しい人のために

  • 状態空間モデリングで重要になる推定量は、観測値Y _ {1:T}=\lbrace y _ t\rbrace _ {t=1}^{\rm T}をもとにした状態\alpha _ tの事後分布P(\alpha _ t|Y _ {1:T})と（長期）予測分布P(\alpha _ {T+i}|Y _ {1:T})である。
  • 逐次推定手法では、以下のステップで事後分布P(\alpha _ t|Y _ {1:T})と予測分布P(\alpha _ {T+i}|Y _ {1:T})を推定する
    • 一期先予測分布P(\alpha _ t|Y _ {1:t-1})とフィルタリング分布P(\alpha _ t|Y _ {1:t})を交互に求めながら、時間順方向に１ステップずつ前進する。
    • P(\alpha _ T|Y _ {1:T})まで求めたら時間逆方向に１ステップずつ遡り事後分布P(\alpha _ t|Y _ {1:T})を推定する（平滑化）
    • P(\alpha _ {T}|Y _ {1:T})から時間順方向に１ステップずつ長期予測P(\alpha _ {T+i}|Y _ {1:T})を行う

f:id:tatamiyatamiyatatatamiya:20210110150615p:plain

  • 線型ガウスのモデルの場合、サンプリングや近似に頼らず推定値・予測値の分布を求められる（カルマンフィルタ）
    • 共分散行列と平均値の漸化式を解いていく

問題設定

系列データY_{1:T} = \lbrace y _ t \rbrace _ {t=1}^{\rm T}が観測値として与えられているとします。

状態空間モデルでは、各観測値の背後にマルコフ過程に従って状態遷移する状態\lbrace \alpha _ t\rbrace_{t=1}^{\rm T}が存在していると仮定します。

{ \begin{aligned} y_t &\sim P(y_t|\alpha_t)\\ \alpha_{t+1} &\sim P(\alpha_{t+1}|\alpha_t) \end{aligned} }

このような仮定のもと、全観測値Y _ {1:T}から状態\lbrace \alpha _ t\rbrace _ {t=1}^{\rm T}の推定と、未来の状態\alpha _ {t+i}, \, i>0の予測を行います。

求めるのは、下記の条件付き確率分布です：

状態推定：

{ \begin{aligned} P(\alpha_t|Y_{1:T}), \quad 0\le t \le T \end{aligned} }

予測：

{ \begin{aligned} P(\alpha_{T+i}|Y_{1:T}), \quad i>0 \end{aligned} }

f:id:tatamiyatamiyatatatamiya:20210110150813p:plain

なお、上記の考え方による状態推定は、\lbrace y _ t \rbrace _ {t=1}^{\rm T}から観測に伴うノイズを除去したスムーズな系列\lbrace \alpha _ t \rbrace _ {t=1}^{\rm T}を求めることにも使われるため、「平滑化」と表現されます。

逐次推定法の考え方

ここでは、事後分布P(\alpha _ t|Y _ {1:T})と予測分布P(\alpha _ {T+i}|Y _ {1:T})を、逐次的に推定する方法の概要を解説します。

逐次推定・予測では、状態の時間発展がひとつ前の時刻の状態にのみに依存するという設定を利用し、時間方向に1ステップずつ推定・予測を行っていきます。

全体の処理は、以下の3つのまとまりに分けることができます：

  • 一期先予測+フィルタリング
    • P(\alpha _ {t-1}|Y _ {1:t-1})からP(\alpha _ {t}|Y _ {1:t-1})を求める（一期先予測）
    • P(\alpha _ {t}|Y _ {1:t-1})からP(\alpha _ {t}|Y _ {1:t})を求める（フィルタリング）
    • t=1から順に上記を交互に繰り返し、P(\alpha _ T|Y _ {1:T})を求める
  • 平滑化
    • P(\alpha _ T|Y _ {1:T})からP(\alpha _ {T-1}|Y _ {1:T})を求める
    • 同様にして時間逆方向に遡り、全t=T,T-1, ..., 3, 2, 1についてP(\alpha _ t|Y _ {1:T})を求める
  • 長期予測
    • P(\alpha _ T|Y _ {1:T})からP(\alpha _ {T+1}|Y _ {1:T})を求める
    • 同様にして一ステップずつ進み、P(\alpha _ {T+i-1}|Y _ {1:T})からP(\alpha _ {T+i}|Y _ {1:T})を求めていく

上記のうち、一期先予測・フィルタリング〜平滑化の流れを図にして表すと以下のようになります。

f:id:tatamiyatamiyatatatamiya:20210110151047p:plain

一期先予測とフィルタリング

時刻tまでの観測値Y _ {1:t} = \lbrace y _ i \rbrace _ {i=1}^{\rm T}が与えられているもとで、次の時刻の状態\alpha _ {t+1}を予測することを考えます。

これは、時刻tまでの観測値Y _ {1:t}をもとに時刻tにおける状態の推定分布P(\alpha _ t|Y _ {1:t})が得られていれば、状態の発展式\alpha _ {t+1} \sim P(\alpha _ {t+1}|\alpha _ t)を用いて以下のように表せます：

一期先予測：

{ \begin{aligned} P(\alpha_{t+1}|Y_{1:t}) &= \int P(\alpha_{t+1}|\alpha_{t}) P(\alpha_{t}|Y_{1:t})d\alpha_{t} \end{aligned}\tag{1}\label{eq:one_step_prediction} }

では、P(\alpha _ t | Y _ {1:t})をどう求めれば良いかですが、以下のように書きなおします：

フィルタリング：

{ \begin{aligned} P(\alpha_t|Y_{1:t})&=P(\alpha_t|y_t, Y_{1:t-1})= \frac{P(\alpha_t, y_t|Y_{1:t-1})}{P(y_t|Y_{1:t-1})}\\ &=\frac{P(y_t|\alpha_t)P(\alpha_t|Y_{1:t-1})}{\int P(y_t|\alpha_t)P(\alpha_t|Y_{1:t-1}) d\alpha_t} \end{aligned}\tag{2}\label{eq:filtering} }

これは、前時刻までのデータY _ {1:t-1}=\lbrace y _ i \rbrace _ {i=1}^{t-1}から予測した\alpha _ tの分布P(\alpha _ t|Y _ {1:t-1})とP(y _ t|\alpha _ t)が分かっていれば求めることができます。

P(y _ t|\alpha _ t)については観測方程式y _ t \sim P(y _ t|\alpha _ t)を用いるのですが、この時y _ tには観測値が入ります。

そして、P(\alpha _ t|Y _ {1:t-1})を求めるのには、再度一期先予測の式を用います。

以上から、時刻t=1から始めてフィルタリングと一期先予測を1ステップずつ繰り返していくことで、最終的に全時刻t=1, 2, ..., TでP(\alpha _ t|Y _ {1:t-1}), P(\alpha _ t|Y _ {1:t})が求まります（下図）。

f:id:tatamiyatamiyatatatamiya:20210110151311p:plain

なお、初手t=1については一期先予測P(\alpha _ 1|Y _ {1:0})が存在しないため、適当な分布P(\alpha _ 1)を仮定して用います。1

平滑化

上記の一期先予測でP(\alpha _ T|Y _ {1:T})まで求まったら、今度は時間逆方向に折り返してP(\alpha _ t|Y _ {1:T})を求めていきます。

状態発展式を逆向きにしたものP(\alpha _ {T-1}|\alpha _ T, Y _ {1:T})がわかれば、P(\alpha _ T|Y _ {1:T})からP(\alpha _ {T-1}|Y _ {1:T})を求めることができます：

{ \begin{aligned} P(\alpha_{T-1}|Y_{1:T}) = \int P(\alpha_{T-1}|\alpha_T, Y_{1:T})P(\alpha_T|Y_{1:T})d\alpha_T \end{aligned} }

P(\alpha _ {T-1}|\alpha _ T, Y _ {1:T})=P(\alpha _ {T-1}|\alpha _ T, Y _ {1:T-1})が成り立つので2、これはベイズの定理を使って以下のように求めることができます：

{ \begin{aligned} P(\alpha_{T-1}|\alpha_T, Y_{1:T})&= \frac{P(\alpha_{T-1}, \alpha_T|Y_{1:T-1})}{P(\alpha_T|Y_{1:T-1})}\\ &= \frac{P(\alpha_T|\alpha_{T-1})P(\alpha_{T-1}|Y_{1:T-1})}{P(\alpha_T|Y_{1:T-1})} \end{aligned} }

P(\alpha _ T|\alpha _ {T-1})は状態発展の式から、P(\alpha _ {T-1}|Y _ {1:T-1})とP(\alpha _ T|Y _ {1:T-1})はそれぞれフィルタリングと一期先予測で求めています。

これにより時刻Tにおける推定状態分布P(\alpha _ {T-1}|Y _ {1:T})が求まれば、同様のステップを踏むことでP(\alpha _ {T-2}|Y _ {1:T}), P(\alpha _ {T-3}|Y _ {1:T}), ..., P(\alpha _ {2}|Y _ {1:T}), P(\alpha _ 1|Y _ {1:T})も求めることができます。

以上をまとめると以下の計算を繰り返していけばよいことがわかります：

{ \begin{aligned} P(\alpha_t | Y_{1:T}) &= \int P(\alpha_t|\alpha_{t+1}, Y_{1:t}) P(\alpha_{t+1}|Y_{1:T})d\alpha_{t+1}\\ &= \int \frac{P(\alpha_{t+1}|\alpha_t)P(\alpha_t|Y_{1:t})}{P(\alpha_{t+1}|Y_{1:t})} P(\alpha_{t+1}|Y_{1:T})d\alpha_{t+1}\\ \end{aligned}\tag{3}\label{eq:smoothing} }

f:id:tatamiyatamiyatatatamiya:20210110151632p:plain

長期予測

長期予測についても、時刻Tまでの状態推定が行えていれば簡単に計算できます。

時刻T+1については以下のようになります：

{ \begin{aligned} P(\alpha_{T+1}|Y_{1:T}) = \int P(\alpha_{T+1}|\alpha_T)P(\alpha_T|Y_{1:T})d\alpha_T \end{aligned} }

同様にt=T+2, T+3, ...と繰り返していくことは簡単にできます：

{ \begin{aligned} P(\alpha_{T+i}|Y_{1:T})=\int P(\alpha_{T+i}|\alpha_{T+i-1})P(\alpha_{T+i-1}|Y_{1:T})d\alpha_{T+i-1} \end{aligned}\tag{4}\label{eq:long_prediction} }

f:id:tatamiyatamiyatatatamiya:20210110151740p:plain

なお、観測値の予測分布を求めたい場合は、さらに観測式y _ {T+i} \sim P(y _ {T+i}|\alpha _ {T+i})を使って以下を計算します：

{ \begin{aligned} P(y_{T+i}|Y_{1:T}) = \int P(y_{T+i}|\alpha_{T+i})P(\alpha_{T+i}|Y_{1:T})d\alpha_{T+i} \end{aligned} }

カルマンフィルタ

以上の逐次推定法では条件付き確率分布ベースで議論を行いましたが、途中に出てきた積分は一般には容易に実行できません。 その場合、MCMCによりサンプリングを行うか、または近似的に分布を求める必要があります。

しかし、モデルが線型かつシステムノイズ・観測ノイズがガウス分布に従う場合、サンプリングや近侍に頼らず厳密に事後分布・長期予測分布を求めることができます。

線型ガウスの状態空間モデル

以下のように書けるモデルを考えます。

{ \begin{aligned} y_t &= Z\alpha_t + \varepsilon_t, \quad \varepsilon_t \sim N(0, H)\\ \alpha_{t+1} &= T\alpha_t + R\eta_t, \quad \eta_t \sim N(0, Q) \end{aligned}\tag{5}\label{eq:dlm} }

N(x|\mu, \Sigma)は期待値が\mu, 共分散行列が\Sigmaの正規分布を表します。 観測値y_tと状態\alpha_tはともにベクトルで表現しており、Z, H, T, Rはいずれも行列です。

これをy_t, \alpha_{t+1}の確率分布の形で書き直すと以下のようになります：

{ \begin{aligned} y_t &\sim N(y_t|Z\alpha_t, H)\\ \alpha_{t+1} &\sim N(\alpha_{t+1}|T\alpha_t, RQR^{\rm T}). \end{aligned}\tag{6}\label{eq:dlm_prob} }

f:id:tatamiyatamiyatatatamiya:20210110151857p:plain

前章までの表記とは、

{ \begin{aligned} P(y_t|\alpha_t) &= N(y_t|Z\alpha_t, H)\\ P(\alpha_{t+1}|\alpha_t) &= N(\alpha_{t+1}|Z\alpha_t, RQR^{\rm T}) \end{aligned}\tag{7}\label{eq:obs_state_prob} }

のように対応がつきます。

このモデルをみると、観測・状態発展ともにガウス分布で表されるほか、観測値y_tの期待値は状態\alpha_tと線型な関係にあり、また、\alpha_{t+1}の時間発展も平均的には前の時刻における値と線型な関係にあります。

このようなモデルを、線型ガウス状態空間モデルと呼ぶことにします。

なお、前回の記事で紹介した、典型的な3種類のモデル設定は、全て線型ガウス状態空間モデルに含まれます。

一期先予測における平均ベクトルと共分散行列の漸化式

このような線型ガウスのモデルについて、前章の要領で事後分布P(\alpha_t|Y _ {1:T})と長期予測分布P(\alpha _ {T+i}|Y _ {1:T})を求めていきます。

そのためには式(\ref{eq:one_step_prediction})〜(\ref{eq:long_prediction})に先ほどの式(\ref{eq:obs_state_prob})を代入していけば良いのですが、初期分布P(\alpha_1)にもガウス分布を仮定すると、途中全てガウス分布同士の積・積分なので、得られる事後分布、長期予測分布ほか、途中で求める一期先予測分布、フィルタリング分布も全てガウス分布になります。

これにより、平均ベクトル、共分散行列の漸化式を導くことができます。

例えば、フィルタリング分布が以下のように書けるとします：

{ \begin{aligned} P(\alpha_t|Y_{1:t}) &= N(\alpha_t|a_{t|t}, V_{t|t})\\ \end{aligned} }

すると、式(\ref{eq:one_step_prediction})と式(\ref{eq:obs_state_prob})から一期先予測分布は、

{ \begin{aligned} P(\alpha_{t+1}|Y_{1:t}) &= \int P(\alpha_{t+1}|\alpha_{t}) P(\alpha_{t}|Y_{1:t})d\alpha_{t}\\ &=\int N(\alpha_{t+1}|Z\alpha_t, RQR^{\rm T})N(\alpha_t|a_{t|t}, V_{t|t})d\alpha_t\\ &=N(\alpha_{t+1}|Za_{t|t}, RQR^{\rm T} + ZV_{t|t}Z^{\rm T}) \end{aligned} }

となります。

途中、

{ \begin{aligned} \int N(y|Ax+b, \Lambda_y)N(x|\mu, \Lambda_x)dx=N(y|A\mu+b, \Lambda_y + A\Lambda_xA^{\rm T}) \end{aligned} }

という関係を使いました（PRML 2章参照）。

ここで、一期先予測分布が

{ \begin{aligned} P(\alpha_{t+1}|Y_{1:t}) &= N(\alpha_{t+1}|a_{t+1|t}, V_{t+1|t})\\ \end{aligned} }

のように書けたとすると、

{ \begin{aligned} \begin{cases} a_{t+1|t}&=Ta_{t|t}\\ V_{t+1|t} &= TV_{t|t}T^{\rm T}+RQR^{\rm T}\\ \end{cases} \end{aligned} }

という関係が成り立つことがわかります。

平均ベクトルと共分散行列の漸化式（結果のみ）

上記と同様にして、

事後分布：

{ \begin{aligned} P(\alpha_t|Y_{1:T}) &= N(\alpha_t|\hat{a}_t, \hat{V}_t)\\ \end{aligned} }

長期予測分布：

{ \begin{aligned} P(\alpha_{T+i}|Y_{1:T}) = N(\alpha_{T+i}|\hat{\alpha}_{T+i}, \hat{V}_{T+i}) \end{aligned} }

のようにおきます。

そうすると、式(\ref{eq:one_step_prediction})〜(\ref{eq:long_prediction})から、以下のような平均ベクトルと共分散行列の漸化式が得られます：

一期先予測＋フィルタリング

{ \begin{aligned} \begin{cases} a_{t+1|t}&=Ta_{t|t}\\ V_{t+1|t} &= TV_{t|t}T^{\rm T}+RQR^{\rm T}\\ \end{cases} \end{aligned} }
{ \begin{aligned} \begin{cases} a_{t|t} &= a_{t|t-1}+V_{t|t-1}Z^{\rm T}F^{-1}(y_t - Za_{t|t-1})\\ V_{t|t} &= V_{t|t-1}-V_{t|t-1}Z^{\rm T}F^{-1}ZV_{t|t-1}\\ \end{cases} \end{aligned} }

ただし、

{ \begin{aligned} F \equiv H + ZV_{t|t-1}Z^{\rm T} \end{aligned} }

平滑化

{ \begin{aligned} \begin{cases} \hat{\alpha}_t &= a_{t|t} + A_t(\hat{\alpha}_{t+1} - a_{t+1|t})\\ \hat{V}_t &= V_{t|t} + A_t(\hat{V}_{t+1} - V_{t+1|t})A_t^{\rm T} \end{cases} \end{aligned} }

ただし、

{ \begin{aligned} A_t = V_{t|t}T^{\rm T}V_{t+1|t}^{-1} \end{aligned} }

長期予測

{ \begin{aligned} \begin{cases} \hat{\alpha}_{T+i}&=T\hat{\alpha}_{T+i-1}\\ \hat{V}_{T+i} &= T\hat{V}_{T+i-1}T^{\rm T}+RQR^{\rm T}\\ \end{cases} \end{aligned} }

まとめ

今回の記事では、逐次推定法の考え方から出発して、カルマンフィルタについて紹介しま た。

逐次推定法は

  • 一期先予測+フィルタリング
  • 平滑化
  • 長期予測

の3種類の操作から構成されます。

特に線型ガウス状態空間モデルであれば、カルマンフィルタといって、推定・予測分布をサンプリングや近似に頼らず解析的に求めることができました。

今回はカルマンフィルタの漸化式の導出は一部除いて行いませんでしたが、こちらについてはまた別の機会に記事にできたらと考えています。

次回の記事では今回解説した理論をベースに、RとPythonのライブラリを利用したモデリングを典型的なモデル設定を題材に行っていきたいと思います。

参考文献

  1. 初期状態分布P(\alpha _ 1)については、事前知識を仮定しない場合は標準偏差が無限大の裾野の広い分布を使うのが一般的です（散漫初期化）。詳細は野村2016などを参照してください。↩

  2. \alpha _ Tで条件づけられているため、\alpha _ {T-1}とy _ Tが独立であることに着目しました。この理屈を理解するためには、PRML 8章のグラフィカルモデリングの章を参照してください。↩