雑記：カルマンフィルタとか何とかカルマンフィルタとか

参考文献：データ同化入門 (予測と発見の科学) | 樋口知之 |本 | 通販 | Amazon
文字の置き方と表式をだいたい上の本に準拠していますが、違うこともあります。導出の仕方は上の本とは異なります。間違っていたらご指摘ください。今回はカルマンフィルタの話しかないです。キャラクターに元ネタはないです。

次回：まだ

f:id:cookie-box:20180513082851p:plain:w60

カルマンフィルタとその仲間たちみたいな何とかフィルタって色々ありますよね。

f:id:cookie-box:20180513082908p:plain:w60

逐次ベイズ推定のアルゴリズムのこと？確かに色々あるね。拡張カルマンフィルタとか、粒子フィルタとか。

その辺ってなんかわかった気になってしまいませんか？線形でベイジアンだったらちゃんと解けばよくて、そうじゃなかったら現在の $x_t$ の周りで1次形近似したり、 $x_t$ の分布を粒子のアンサンブルで近似すればいいんだ、って感じで…。

なんで！？いうほどわかった気にならないよ？

というか、そもそも何をしたかったのかもよくわからなくなってしまって…。

目的は状態の推定だよね！？

よく ARIMA のような伝統的な時系列モデルから状態空間モデルへ、って説明されますけど、話が変わってる気がするんですよね。ARIMA では AR や MA の各次数の係数を求めて、状態空間モデルは状態の分布を推定していくんですよね。なんか話が逆になってませんか？モデル $f(x_t)$ に対して時系列データ $x_t$ を所与として $f$ の係数を求めたいのか、 $f$ の形は所与として $x_t$ を推定したいのか…。

うーん…状態って言葉が曖昧なのかな？観測値の真の値であるかのように説明されることも多いけど、というよりは状態って「トラッキングしたいもののベクトル」なんだよね。所与なのはあくまで観測値 $y_t$ で、ARIMA ではむしろ AR や MA の各次数の係数たちが $x_t$ で（ふつう時不変だから $x$ だけど）。状態空間モデルではモデルの形 $f$ は所与とするけどもしそのパラメータが未知ならそれは $x_t$ 側に含めると思うよ。なんていうか結局「確信しているものを所与として不確かなものを推定したい」だね。トートロジーだけど。

なるほど…。じゃあ、このブログの2018年4月23日の記事にある、「（状態空間モデルで）もし『状態』が観測される変数そのもので、『観測モデル』が恒等写像だったら（伝統的な時系列モデルと）一緒の構図になる」というのは対応がおかしかったんですね。より正しくは、「もし『状態』が ARIMA モデルの係数（及び必要なステップだけの過去から現在までのノイズや観測値）で、『観測モデル』がそれらの線形和（ARIMA モデルの形）だったら一緒の構図になる」ですね。観測モデルは観測といっても日常的な意味のように「観測時に入り混じる誤差を考慮したモデル」ではなく、「不確かな変数を知覚できる数値に変換するモデル」なのに、筆者はそこを踏み外して、思考停止でシステムモデルに ARIMA を代入していたんですね。

う、うん。たぶんそうだね…。あくまでそういう状態空間モデルにすれば ARIMA と「係数を求めたい」という目的が大雑把にそろうというだけで、実際には ARIMA モデルの方では逐次的に係数の分布を更新してはいかないと思うけどね。あと、もしそのように修正する前の元の文章の立場に立つなら、「この時系列データはこの係数のARIMAモデルで時間発展する」ということの方に絶対の確信があって、知覚された時系列データがそれに従っていないならそっちがあやしい、観測ノイズによって真の値とずれているんだろう、真の値が知りたい、という姿勢にはなるね。

不確かなものを推定したいという目的はわかったと思います。では早速各種フィルタについてちゃんと追って行きたいと思うんです。まず一般的に、一期先予測とフィルタリングは以下の式で表せます。
　　 $\displaystyle p(x_t | y_{1:t-1}) = \int p(x_t | x_{t-1}) p(x_{t-1} | y_{1:t-1}) dx_{t-1}$ 　（一期先予測分布）
　　 $\displaystyle p(x_t | y_{1:t}) = \frac{p(y_t | x_t) p(x_t | y_{1:t-1})}{p(y_t | y_{1:t-1})} = \frac{p(y_t | x_t) p(x_t | y_{1:t-1})}{\displaystyle \int p(y_t | x_t) p(x_t | y_{1:t-1}) dx_t}$ 　（フィルタ分布）

また端折ったね…。まず前提として、ここでは時刻 $t-1$ までの観測値 $y_{1:t-1}$ が得られたもとでの状態の分布 $p(x_{t-1} | y_{1:t-1})$ を所与とするんだね。かつ、時刻 $t-1$ の状態を $x_{t-1}$ と固定したもとでの次の時刻の状態の分布 $p(x_t | x_{t-1})$ も計算できるとする。これはシステムモデルを決めておきなさいってことだね。これらがわかれば、一期先予測の式は $x_{t-1}$ の分布にわたって次の時刻の $x_t$ がどうなるかの分布を積分するだけ。それで次にフィルタの式は、現時点までの観測値 $y_{1:t}$ が出そろったもとでの状態の分布 $p(x_t | y_{1:t})$ はどうなっているか、ということだけど…これは観測モデル $p(y_t | x_t)$ とさっきの一期先予測の式 $p(x_t | y_{1:t-1})$ の積からベイズの定理により求まるね。ただ、単にこれらの積だと確率分布が正規化されないから、 $p(y_t | y_{1:t-1})$ で割る必要がある。これは、 $p(y_t | y_{1:t-1})$ は $p(y_t | x_t) p(x_t | y_{1:t-1})$ を $x_t$ の分布にわたって積分すれば求まる。あと、ここまで暗に、下図のグラフィカルモデルで表される従属構造が仮定されていることに注意が必要だね。つまり、 $p(x_t | x_{1:t-1}, y_{1:t-1}) = p(x_t | x_{t-1})$ 及び $p(y_t | x_{1:t}, y_{1:t-1}) = p(y_t | x_t)$ が成り立つ。マルコフ性ともいうね。そうじゃなきゃ上の式にならない。

※ ここに状態空間モデルのグラフィカルモデルを貼る。

あ、はい、色々飛ばしてました…補足ありがとうございます。それでまずカルマンフィルタの導出です。ここからは、システムモデルと観測モデルがそれぞれ以下のようにかけると仮定しています。
　　 $\begin{cases} x_t = F_t x_{t-1} + G_t v_t \; \; \; \; & v_t \sim N_m(0, Q_t) \\ y_t = H_t x_t + w_t & w_t \sim N_l(0, R_t)\end{cases}$
上式に出てくる以下の文字は以下のような次元のベクトルや行列です。 $x_t \in \mathbb{R}^k, \; y_t \in \mathbb{R}^l, \; F_t \in \mathbb{R}^{k \times k}, \; G_t \in \mathbb{R}^{k \times m}, \; H_t \in \mathbb{R}^{l \times k}, \; Q_t \in \mathbb{R}^{m \times m}, \; R_t \in \mathbb{R}^{l \times l}$
このうち $Q_t$ と $R_t$ は正規ノイズの分散共分散行列ですね。それで、 $x_t, \; v_t, \; y_t, \; w_t$ が最初から決まっていない、確率変数です。目的は、 $x_t$ の分布をトラッキングすることです。

うん、観測値より状態の次元の方が大きいなら $H_t$ は横長の、逆なら縦長の行列になるね。

そうですね。それで、先に結論を書きます。
　　 $x_{t-1} \, | \, y_{1:t-1} \sim N_k (\mu_{t-1|t-1} , \, V_{t-1|t-1})$
を仮定すると、
　　 $x_{t} \, | \, y_{1:t-1} \sim N_k (\mu_{t|t-1} , \, V_{t|t-1}) = N_k (F_t \mu_{t-1|t-1} , \, F_t V_{t-1|t-1} {F_t}^{\top} + G_t Q_t {G_t}^{\top})$ 　（一期先予測分布）
　　 $x_{t} \, | \, y_{1:t} \sim N_k (\mu_{t|t} , \, V_{t|t}) = N_k (\mu_{t|t-1} + K_t e_t , \, V_{t|t-1} - K_t H_t V_{t|t-1})$ 　（フィルタ分布）
但し、
　　 $e_t = y_t - H_t \mu_{t|t-1} \; \in \mathbb{R}^l$
　　 $K_t = V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1} \; \in \mathbb{R}^{k \times l}$
が成り立ちます。時刻 $t-1$ までの観測値 $y_{1:t-1}$ が得られたもとで推定されている状態の分布がガウシアンならば、一期先予測もガウシアンであり、さらに時刻 $t$ のフィルタ分布もガウシアンであるということです。かつ、その一期先予測とフィルタ分布の平均ベクトルと分散共分散行列は上式のように1つ手前の平均ベクトルや分散共分散行列の式として陽に書き下せます。

有名な結論だね。

でもどうしてこういうことになるのか。

まず一期先予測には以下をつかおうか。つまり、多変量正規分布にしたがう確率ベクトルは線形変換しても多変量正規分布にしたがうし、多変量正規分布にしたがう確率ベクトルどうしを足しても多変量正規分布にしたがう。証明は、前者は $M_{Ax}(t) \equiv E(\exp(t^{\top} Ax))=M_x(A^{\top} t)$ を、後者は、互いに独立な確率ベクトル $x_1, \, x_2$ について $M_{x_1 + x_2}(t) = M_{x_1}(t) M_{x_2}(t)$ を利用すれば導出できるかな。あ、 $M_x(t)$ はモーメント母関数で、これは確率分布に対して一意に定まる関数で、多変量正規分布のモーメント母関数は $M_x(t) = \exp(t ^{\top} \mu + \frac{1}{2} t ^{\top} \Sigma t)$ だね。

$x \sim N_k (\mu_x, \, \Sigma_x)$ で $y = Ax$ のとき、 $y \sim N_l (A \mu_x, \, A \Sigma_x A^{\top})$

$x_1$ と $x_2$ が独立で、 $x_1 \sim N_k (\mu_{x_1}, \, \Sigma_{x_1}), \; x_2 \sim N_k (\mu_{x_2}, \, \Sigma_{x_2})$ のとき、 $x_1 + x_2 \sim N_k (\mu_{x_1} + \mu_{x_2}, \, \Sigma_{x_1} + \Sigma_{x_1})$

それらをつかう……？えっと、いま $x_{t-1} \, | \, y_{1:t-1}$ が多変量正規分布にしたがうと仮定しているので、その線形変換である $F_t x_{t-1} \, | \, y_{1:t-1}$ も多変量正規分布にしたがう。それとは独立にノイズ $G_t v_t$ も多変量正規分布にしたがうから… $x_{t} \, | \, y_{1:t-1} = F_t x_{t-1} \, | \, y_{1:t-1} + G_t v_t$ は多変量正規分布にしたがいますね。それで、その平均ベクトルと分散共分散行列は、
　　 $\mu_{t|t-1} \equiv {\rm E}(x_{t} \, | \, y_{1:t-1}) = {\rm E}(F_t x_{t-1} \, | \, y_{1:t-1}) + E(G_t v_t) = F_t \mu_{t-1|t-1}$
　　 $V_{t|t-1} \equiv {\rm Var}(x_{t} \, | \, y_{1:t-1}) = {\rm Var}(F_t x_{t-1} \, | \, y_{1:t-1}) + {\rm Var}(G_t v_t) = F_t V_{t-1|t-1} {F_t}^{\top} + G_t Q_t {G_t}^{\top}$

一期先予測分布が導出できたね。フィルタ分布はベイズの定理を使うとどうなる？さしあたり正規化因子は考えなくていいよ。

ベイズの定理をつかうと、
$\begin{split} p(x_{t} \, | \, y_{1:t}) &\propto p(y_t | x_t) p(x_t | y_{1:t-1}) \\ & = \exp \left( -\frac{1}{2} (y_t - H_t x_t)^{\top} {R_t}^{-1} (y_t - H_t x_t) \right) \exp \left( -\frac{1}{2} (x_t - \mu_{t|t-1} )^{\top} {V_{t|t-1}}^{-1} (x_t - \mu_{t|t-1} ) \right) \\ & = \exp \left( -\frac{1}{2} (y_t - H_t x_t)^{\top} {R_t}^{-1} (y_t - H_t x_t) -\frac{1}{2} (x_t - \mu_{t|t-1} )^{\top} {V_{t|t-1}}^{-1} (x_t - \mu_{t|t-1} ) \right) \end{split}$
ここからどうすれば…。

$\begin{split} (y_t - H_t x_t)^{\top} {R_t}^{-1} (y_t - H_t x_t) + (x_t - \mu_{t|t-1} )^{\top} {V_{t|t-1}}^{-1} (x_t - \mu_{t|t-1} ) \end{split}$
この式を $x_t$ があらわれない項は無視して展開して整理すると、

$\begin{split} & - {y_t}^{\top} {R_t}^{-1} H_t x_t - {(H_t x_t)}^{\top} {R_t}^{-1} y_t + {(H_t x_t)}^{\top} {R_t}^{-1} H_t x_t + {x_t}^{\top} {V_{t|t-1}}^{-1} x_t - {x_t}^{\top} {V_{t|t-1}}^{-1} \mu_{t|t-1} - {\mu_{t|t-1}}^{\top} {V_{t|t-1}}^{-1} x_t \\ &= {x_t}^{\top} ({H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1}) x_t - {x_t}^{\top} ({H_t}^{\top} {R_t}^{-1} y_t + {V_{t|t-1}}^{-1} \mu_{t|t-1}) - ({y_t}^{\top} {R_t}^{-1} H_t - {\mu_{t|t-1}}^{\top} {V_{t|t-1}}^{-1}) x_t \end{split}$

こうなるけど、 $x_t$ に挟まれている ${H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1}$ は、逆行列補題より逆行列 $V_{t|t-1} - ( V_{t|t-1} {H_t}^{\top} (H_t V_{t|t-1} {H_t}^{\top} + R_t)^{-1}) H_t V_{t|t-1}$ をもつ。加えて、 ${x_t}^{\top}$ に右側からかかっている行列と $x_t$ に左側からかかっている行列が対称だから、上の式は結局 $\begin{split} & ({x_t} - \alpha)^{\top} ({H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1}) (x_t -\alpha) \end{split}$ の形に整理することができる。これを展開して係数比較すれば $\alpha$ が求まるね。

※ ここに逆行列補題をかく。

そして、 $\begin{split} & ({x_t} - \alpha)^{\top} ({H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1}) (x_t -\alpha) \end{split}$ の形に整理できるということは、フィルタ分布 $p(x_{t} \, | \, y_{1:t})$ が多変量正規分布にしたがうということに他ならないよね。平均 $\alpha$ で、分散共分散行列 $({H_t}^{\top} {R_t}^{-1} H_t + {V_{t|t-1}}^{-1})^{-1}$ のね。

…ごめん寝てた。

なんで寝るの！？部長のために式展開したのに！

次のフィルタの話に移りたいです。ざっとモンテカルロ手法まで復習したく…というのも、そろそろ部活動の予算申請の時期なので、我らがベイズ統計部も膨大な計算資源の必要性を申請書にしたためなければ。

部活動予算でどんなマシン買う気なの！？運動部じゃないんだからそんな予算下りないよ？

（次回があれば）つづく