何も解らない。。

データサイエンス関連

カルマンフィルター

【カルマンフィルター】

観測方程式: \begin{equation} \begin{split} Y_{t} = \beta_{t}X_{t} + e_{t} \end{split} \end{equation}

状態方程式: \begin{equation} \begin{split} \beta_{t} = T\beta_{t-1} + \epsilon_{t} \end{split} \end{equation}

この時、カルマンゲインを \begin{equation} \begin{split} K_{t} := \frac{X_{t}Σ_{t|t-1}}{X_{t}^2Σ_{t|t-1}+σ_{e}^2} \end{split} \end{equation} と表すと、状態変数の期待値、分散のフィルタリングは以下、

\begin{equation} \begin{split} \hat{\beta_{t|t}} = \hat{\beta_{t|t-1}} + K_{t}(y_{t} - \hat{Y_{t|t-1}}) \cr \hat{Σ_{t|t}} = (1 - X_{t}K_{t}) \hat{Σ_{t|t-1}} \end{split} \end{equation}

で表せる。

適当な紹介にも程があります。順を追ってみていきましょう。

今回は以下の本を参考にまとめます。

はじめに

本記事の構成ですが、まず

  1. 状態空間モデルについて簡単に説明
  2. フィルタリングの直感的理解
  3. フィルタリングの公式の導出

という流れになっています。

尚、上記のフィルタリングの式の導出では、簡単のため、幾つかの仮定*1をおいたものになっております。あしからず。

そもそも状態空間モデルってなんだっけ

観測した結果についての関係式”観測方程式”と、観測できない状態に関する関係式”状態方程式”からなる時系列分析モデル、です。

…と言われてもなかなかピンと来ませんよね。。2つの式を再掲します。

観測方程式:  Y_{t} = \beta_{t}X_{t} + e_{t}

状態方程式 \beta_{t} = T\beta_{t-1} + \epsilon_{t}

具体的に考えてみます。

例えば、消費関数*2を例に考えてみます。 Y_{t}が消費、 X_{t}が所得、 \beta_{t}を限界消費性向とします。よく見る式だと上の観測方程式オンリーで  \beta_{t}が定数  \betaになったものではないでしょうか。

 Y_{t} = \beta X_{t} + e_{t}

つまり状態方程式のないノーマルな消費関数の式*3では限界消費性向が時間で不変であると仮定してモデリングをしていることになります。

これを状態方程式を考えることで、観測できない状態(限界消費性向)を時間で変化するものとして扱うことができます。このように状態空間モデルを使うことで、より柔軟なモデリングが可能になるというわけです。

フィルタリングとは(直感的な理解)

カルマンフィルターの目的は、確率変動する状態変数 \beta_{t}の平均と分散を推定することにあります。

本記事ではそのような推定の一つである”フィルタリング”をメインに解説しますが、他にも”一期先予測”と”スムージング”という推定があります。これらの直感的な関係性を参考書籍に出てくる例で確認します。本の中ではこれらの推定を”ミステリー小説の犯人予想”に喩えて説明しています。

まず一期先予測ですが、これは分かりやすい推定で、例えば”99ページまで本を読んだ状態で100ページ目で誰が犯人か予測する”というものです。 次にスムージングは、本を全て読み終えた状態で100ページ目*4での犯人予測がなぜ間違いだったか(正解だったか)を振り返るような作業にあたります。最後に本記事のメインのフィルタリングは、100ページ目で得た情報をもとに一期先予測の修正を行うような作業、と考えればなんとなく理解ができると思います。

ベイジアンアプローチによるフィルタリング公式の導出

さて、登場人物の説明が概ね完了したので、最後に実際にどう推定(フィルタリング)するのか、というところを詰めていきましょう。

参考書籍では、正規性を仮定し導出、正規性を仮定せず平均分散から導出、ベイズの定理から導出、の3タイプの導出が掲載されています。*5 本記事では、この中から”ベイズの定理から導出”(ベイジアンアプローチ)を紹介することにします。

前提条件の確認

まず、前提条件について確認しましょう。観測、状態方程式の誤差項について以下を仮定します。

\begin{equation} \begin{split} e_{t} \propto N(0, σ_{e}^2), \epsilon_{t} \propto N(0, σ_{\epsilon}^2) \cr E[e_{t}\epsilon_{s}] = 0 \qquad \forall s, t \cr E[e_{t}e_{s}] = 0, E[\epsilon_{t}\epsilon_{s}] = 0 \quad \forall s \neq t \end{split} \end{equation}

誤差項は正規分布に従い、系列相関なし、観測、状態方程式のそれぞれの誤差同士も相関を持たないということですね。 この前提のもと、時点 t-1で時点 tの予測を行って、時点 tで得られた観測値Y_{t}をもとに予測値の更新を行う…という作業を逐次的に行っていくことになります。

ベイズの式とカルマンフィルターの関係

次に今回の記事の主役であるベイズの式とカルマンフィルターの関係性を考えると以下のようになります。

 Pr(\beta_{t}|Y_{t}=y_{t})\propto Pr(Y_{t}|\beta_{t})Pr(\beta_{t})

カルマンフィルタの文脈で言うと、事前に予測していた状態変数(事前分布)を新たに観測した情報(尤度)をもとに現時点での状態変数の予測(事後分布)に更新(フィルタリング)している式だと読み取れます。

それでは事前分布と尤度を具体的に見ていきます。

事前分布

事前分布 Pr(\beta_{t})は、 Y_{t} = yが観測される前の時点( t-1)での \beta_{t}の分布です。

期待値と分散は時点 t-1で活用できる情報を使って推定される値、 \hat{\beta_{t|t-1}}, \hat{Σ_{t|t-1}}で、仮定よりこの分布は正規分布に従うため、  \hat{\beta_{t}}|\Omega_{t-1} \sim N(\hat{\beta_{t|t-1}}, \hat{Σ_{t|t-1}})となります。

事前分布:  Pr(\beta_{t}) = \frac{1}{\sqrt{2\pi \hat{Σ_{t|t-1}}}} exp(-\frac{1}{2} \frac{ {(\beta_{t} - \hat{\beta_{t|t-1}})}^2 }{\hat{Σ_{t|t-1}} })

尤度

尤度 Pr(Y_{t}|\beta_{t})は、 \beta_{t}がわかった状態で観測変数 Y_{t}を考えます。

観測方程式:  Y_{t} = \beta_{t}X_{t} + e_{t}

であることを思い出すと、尤度は Y_{t}|\beta_{t} \sim N(\beta_{t}X_{t}, σ_{e}^2)となります。

尤度:  Pr(Y_{t}|\beta_{t}) = \frac{1}{\sqrt{2\pi σ_{e}^2}} exp(-\frac{1}{2} \frac{ {(Y_{t} - \beta_{t}X_{t})}^2 }{ σ_{e}^2 })

事後分布

最後に事後分布ですが、事後分布は2通りの方法で表現出来ます。

まずは今まで見てきた事前分布、尤度を使って以下のように表現出来ます。

\begin{equation} \begin{split} Pr(\beta_{t}|Y_{t}) \propto \frac{1}{\sqrt{2\pi σ_{e}^2}} exp(-\frac{1}{2} \frac{ {(Y_{t} - \beta_{t}X_{t})}^2 }{ σ_{e}^2 }) \times \frac{1}{\sqrt{2\pi \hat{Σ_{t|t-1}}}} exp(-\frac{1}{2} \frac{ {(\beta_{t} - \hat{\beta_{t|t-1}})}^2 }{\hat{Σ_{t|t-1}} }) \cr \propto exp(-\frac{1}{2}(\frac{X_{t}^2}{σ_{e}^2} + \frac{1}{Σ_{t|t-1}})\beta_{t}^2 + (\frac{X_{t}y_{t}}{σ_{e}^2} + \frac{\beta_{t|t-1}}{Σ_{t|t-1}})\beta_{t}) \end{split} \end{equation}

細かい計算はキャッツアイしますが、指数部分を \betaに関して整理していけばこのようになります。

もう一つの方法は、

条件付き期待値: \hat{\beta_{t|t}} \equiv E[\beta_{t}|Y_{t} = y_{t}]

分散: \hat{Σ_{t|t}} \equiv Var[\beta_{t}|Y_{t} = y_{t}]

を使って表現します。

\begin{equation} \begin{split} Pr(\beta_{t}|Y_{t}) = \frac{1}{\sqrt{2\pi \hat{Σ_{t|t}}}} exp(-\frac{1}{2} \frac{ {(\beta_{t} - \hat{\beta_{t|t}})}^2 }{ \hat{Σ_{t|t}^2} }) \cr \propto exp(-\frac{1}{2}\frac{\beta_{t}^2 -2\beta_{t}\hat{\beta_{t|t}}}{\hat{Σ_{t|t}}}) \cr = exp(-\frac{1}{2} \frac{1}{\hat{Σ_{t|t}}} \beta_{t}^2 + \frac{\hat{\beta_{t|t}}}{\hat{Σ_{t|t}}} \beta_{t}) \end{split} \end{equation}

こちらも粛々と計算するのみですが、ポイントとしては確率変数 \beta_{t}に注目して式変形していくところでしょうか。推定値である、 \hat{\beta_{t|t}} \hat{Σ_{t|t}}は無視します。

これで準備が整いました。あとは上記の2式を比較することで、フィルタリングの公式を導出しましょう。

フィルタリングの公式

上記の事後分布の2通りの表現について、 \beta_{t} \beta_{t}^2に注目して比較することで、以下の式を得ます。

\begin{equation} \begin{split} \frac{1}{\hat{Σ_{t|t}}} = \frac{X_{t}^2}{σ_{e}^2} + \frac{1}{\hat{Σ_{t|t-1}}} \cr \frac{\hat{\beta_{t|t}}}{\hat{Σ_{t|t}}} = \frac{X_{t}y_{t}}{σ_{e}^2} + \frac{\hat{\beta_{t|t-1}}}{\hat{Σ_{t|t-1}}} \end{split} \end{equation}

これを式変形することで冒頭のフィルタリングの式が導出できます。まず、分散のフィルタリングの式から導出しましょう。

上記の最初の式を変形することで以下を得ます。

\begin{equation} \begin{split} \hat{Σ_{t|t}} = \frac{σ_{e}^2}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} \hat{Σ_{t|t-1}} \cr = (1 - X_{t}\frac{X_{t}\hat{Σ_{t|t-1}}}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2}) \hat{Σ_{t|t-1}} \cr = (1 - X_{t}K_{t}) \hat{Σ_{t|t-1}} \end{split} \end{equation}

1行目から2行目はややテクニカルですが、1行目の分子に X_{t}^2 \hat{Σ_{t|t-1}} - X_{t}^2 \hat{Σ_{t|t-1}} = 0を足してみると分かりやすいかもしれません。

最後に期待値のフィルタリングの式を導出します。こちらは少し厄介ですが、頑張って計算します。

\begin{equation} \begin{split} \hat{\beta_{t|t}} = \frac{X_{t} \hat{Σ_{t|t}}}{σ_{e}^2} y_{t} + \frac{\hat{Σ_{t|t}}}{\hat{Σ_{t|t-1}}} \hat{\beta_{t|t-1}} \cr = \hat{Σ_{t|t}} (\frac{X_{t} y_{t} \hat{Σ_{t|t-1}} + \hat{\beta_{t|t-1}} σ_{e}^2}{σ_{e}^2 \hat{Σ_{t|t-1}}}) \cr = \frac{σ_{e}^2}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} \hat{Σ_{t|t-1}} (\frac{X_{t} y_{t} \hat{Σ_{t|t-1}} + \hat{\beta_{t|t-1}} σ_{e}^2}{σ_{e}^2 \hat{Σ_{t|t-1}}}) \cr = \frac{\hat{\beta_{t|t-1}} σ_{e}^2}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} + \frac{X_{t} \hat{Σ_{t|t-1}}}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} y_{t} \cr = \frac{\hat{\beta_{t|t-1}} σ_{e}^2}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} + K_{t}y_{t} + (K_{t}\hat{\beta_{t|t-1}}X_{t} - K_{t}\hat{\beta_{t|t-1}}X_{t}) \cr = \frac{σ_{e}^2}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} \hat{\beta_{t|t-1}} + \frac{X_{t} \hat{Σ_{t|t-1}}}{X_{t}^2 \hat{Σ_{t|t-1}} + σ_{e}^2} X_{t} \hat{\beta_{t|t-1}} + K_{t}(y_{t} - \hat{\beta_{t|t-1}} X_{t}) \cr = \hat{\beta_{t|t-1}} + K_{t}(y_{t} - \hat{Y_{t|t-1}}) \end{split} \end{equation}

こちらの式変形でも K_{t}\hat{\beta_{t|t-1}}X_{t} - K_{t}\hat{\beta_{t|t-1}}X_{t} = 0をうまく使って式変形しています。

以上で目的のフィルタリング式が導出できました。お疲れ様でした。

記念に式を再掲しておきましょう。

\begin{equation} \begin{split} \hat{\beta_{t|t}} = \hat{\beta_{t|t-1}} + K_{t}(y_{t} - \hat{Y_{t|t-1}}) \cr \hat{Σ_{t|t}} = (1 - X_{t}K_{t}) \hat{Σ_{t|t-1}} \end{split} \end{equation}

なんだか愛着が湧いてきましたね。カルマンフィルターは個人的には小難しくてなんだかとっつきにくい印象だったのですが、この参考書籍のおかげで少し仲良くなれた気がします。 具体例なども豊富に取り扱っていて、読むためのバーもそんなに高くないので、とてもおすすめできる逸品です。

それでは今回はこの辺で、ごきげんよう

*1:具体的には、状態変数、観測変数がそれぞれ1つで、正規分布に従うことを仮定しています。正規性を仮定しないベイジアンアプローチについては、谷崎(1993)の第1章、北川(2005)の第14章を参照、とのことです。

*2:上で紹介した本でも事例として出てきます。

*3:別に限界消費性向が定数なものがデファクトというわけではないでしょうけど。。

*4:特に100ページ目でなくても良いです。

*5:とても贅沢な内容になっています。