わカルマンを目指し、まずいちばん簡単な例を作ってみる。
※内容の保証はなく、また予告なく書き換える可能性があります。
セットアップ
まず離散かつ1次元で考えます。
\begin{align}
X_0 &= \mu_0 + w_{0} \\
X_{k+1} &= F X_{k} + w_{k+1} \\
Y_{k} &= X_{k} + v_k
\end{align}
とします。Fは適当な係数で、w_{k},v_{k}はすべて独立な正規分布とし、更に単純化のためにこれらの分散も一様にしましょう。
\begin{align}
w_k &\sim \mathcal{N}(0,\sigma_0^2)\\
v_k &\sim \mathcal{N}(0,\sigma_1^2)
\end{align}
状況としては、X_{k}がX_{k+1} = FX_{k}という時間発展をしたいところ、この時間発展自体がガウスノイズw_kの影響をうけており、それを観測した値がY_{k}である、としたいところ、この観測にも別のガウスノイズv_kの影響がある、という状況です。目標は、Y_0\dots Y_kまでの情報による、X_kの、ある意味で最良の推定値を求めることです。
プロトタイプ:2乗誤差最小と条件付き確率
ここで、ちょっと測度論チックな話をします。次の有名な事実があります。
ある確率空間(\Omega,\mathcal{F})の確率変数Xをとり、これが平均と分散を持つとします。この確率空間の加法族\Sigmaの部分加法族\mathcal{G}\subset \mathcal{F}をとると、\mathcal{G}への条件付き平均を考えることができます。
\mathbb{E} \left[X|\mathcal{G}\right]
これの構成方法としては、L^1の発想としてはG\in \mathcal{G} \mapsto \int_G X(\omega) dP(\omega)を(\Omega,\mathcal{G})上の測度だと思ってRadon-Nikodym微分を取るもの、L^2の発想としてはL^2(\Omega,\mathcal{G})\rightarrow L^2(\Omega,\mathcal{F})の共役をとるものがあると思います。がこれはここでは触れず、そういう事実だけ認めましょう。
この条件付き平均の性質として、実はこれがL^2の意味で射影である、というものがあります。どういうことかというと、\mathbb{E} \left[X|\mathcal{G}\right]は、Y \in L^2(\Omega, \mathcal{G})のうち \mathbb{E}\left[(X-Y)^2\right]を最小にするものです。
\mathcal{G}のような部分加法族は、確率変数の組から生成できます。なので、先の話に沿って言えば、Y_0 \dots Y_kを可測にする最小の加法族を\mathcal{Y}_kとすると、\mathbb{E} \left[X_k|\mathcal{Y}_k\right]は、Z \in L^2(\Omega, \mathcal{Y}_k)のうち \mathbb{E}\left[(X_k-Z)^2\right]を最小にするものです。つまり、2乗誤差を最小にする推定量は、参照してよい観測量の生成加法族の条件付き平均と思って良いということです。
そういうわけで、本稿で今後\mathbb{E} \left[X_k|\mathcal{Y}_k\right]をどのように求めるかを考えます。
測度論ベースではこの条件付き平均は再び\Omega上の関数ですが、\mathcal{G}を具体的な確率変数の生成加法族に取った時は、その確率変数が取る値の上の関数として書けることが知られています。(ドゥーブ=ディンキンの補題)。したがって、\mathbb{E} \left[X_k|\mathcal{Y}_k\right]を、たかだかy_0\dots y_kに依存する変数として求めることを考えます。
条件付き平均は、もし条件付き確率分布が存在するなら、その平均そのものです。なので、そもそものY_0,\dots Y_k条件つきX_k
X_k |_{Y_0=y_0,\dots,Y_k=y_k}
の分布を求めることができれば、その平均として与えられます。
プロトタイプ:2次元ガウス
突然ですが話を変えて、2変数ガウス分布を考えます。なぜかというと、実はここでの計算があとでの導出に本質的に効いてくるからです。
X \sim \mathcal{N}(\mu_1,\sigma_1^2)\\
Y \sim \mathcal{N}(\mu_2,\sigma_2^2)
とします。これは独立です。この上でZ = X + Y \sim \mathcal{N}(\mu_1 + \mu_2 , \sigma_1^2 + \sigma_2^2)を取ります。これはX,Yと独立ではありません。
ここで、Zの値がわかっているときに、Xの値の条件付き分布がどうなるか考えます。これには同時分布を求めて、変数変換して積分すればいいです。が、ガウス積分の変数変換は、はい、つらく単調です。退屈でミスを誘発する計算を飛ばして結果だけ見ると、Z=z条件付きXは、
X|_{Z=z} \sim \mathcal{N}\left(\frac{\sigma_2^2}{\sigma_1^2 + \sigma_2^2}\mu_1 + \frac{\sigma_1^2}{\sigma_1^2 + \sigma_2^2}(z-\mu_2),\frac{\sigma_1^2\sigma_2^2}{\sigma_1^2 + \sigma_2^2} \right)
つまり、再び正規分布になっています。そしてやや複雑ではありますが、その平均と分散はちゃんともとの分布から計算出来ています。特に平均ですが、\mu_2=0の場合を考えてみてください。\mu_2=0のときは、X,Yの分散による、z,\mu_1の内分点になっており、分散は平均に関係なく幾何平均です。
なぜこれが重要なのか、それは離散線形カルマンフィルタは、この条件付き平均を帰納的に回したものと思えるからです。
(最も簡単な)カルマンフィルタ
セットアップを思い出します。
\begin{align}
X_0 &= \mu_0+ w_{0} \\
X_{k+1} &= F X_{k} + w_{k+1} \\
Y_{k} &= X_{k} + v_k\\
w_k &\sim \mathcal{N}(0,\sigma_0^2)\\
v_k &\sim \mathcal{N}(0,\sigma_1^2)
\end{align}
X_k|_{Y_0=y_0,\dots,Y_k=y_k}を計算したいのでした。まずk=0を考えてみましょう。
\begin{align}
X_0 &= \mu_0 + w_0 \\
Y_0 &= \mu_0 + w_0 + v_0
\end{align}
なのだから、先のプロトタイプから明らかに
X_0 |_{Y_0 = y_0} \sim \mathcal{N}\left(\frac{\sigma_1^2}{\sigma_0^2 + \sigma_1^2}\mu_0 + \frac{\sigma_0^2}{\sigma_0^2 + \sigma_1^2}y_0,\frac{\sigma_0^2\sigma_1^2}{\sigma_0^2 + \sigma_1^2}\right)
と求まります。つまり、ガウスです。そもそもノイズの発生源がガウスであり、発展方程式も全部ガウス、ガウスの条件付分布もガウスなので、これからはガウスしか出てくる可能性はありません。そこで
X_k |_{Y_0 = y_0,\dots,Y_k=y_k} = \mathcal{N}(\alpha_k,\beta_k)
とおき、\alpha_k,\beta_kがどうなるかを見ていきましょう。
次はk=1です。
\begin{align}
X_1 &= FX_0 + w_1 \\
Y_1 &= X_1 + v_1
\end{align}
ですが、ここでY_0,Y_1の条件つけを考えます。ここで一気に条件付けを考えるのは難しいので、段階的にいきます。まず、X_1 = F X_0 + w_1について、Y_0条件付けを先に計算しておきます。Y_0はX_0としか相関しないのでそのまま条件付けができます。
\begin{align}
& X_1 |_{Y_0 = y_0,Y_1 = y_1}\\
= & X_1 |_{Y_0=y_0} |_{Y_1 =y_1}\\
= & (FX_0|_{Y_0 = y_0} + w_1) |_{Y_1 = y_1}
\end{align}
X_0|_{Y_0 = y_0}は\mathcal{N}(\alpha_0,\beta_0)だったので、これにw_1を加えてFX_0|_{Y_0 = y_0} + w_1 \sim \mathcal{N}(F\alpha_0,F^2\beta_0+\sigma_0^2)、これとv_1を加えたY_1との条件つき分布なので
\begin{align}
X_1 |_{Y_0 = y_0,Y_1 = y_1} &\sim N(\alpha_1,\beta_1) \\
\alpha_1 &= \frac{\sigma_1^2}{F^2\beta_0 + \sigma_0^2 + \sigma_1^2}F\alpha_0 + \frac{F^2\beta_0 +\sigma_0^2}{F^2\beta_0 + \sigma_0^2 + \sigma_1^2}y_1\\
\beta_1 &= \frac{(F^2\beta_0 +\sigma_0^2)\sigma_1^2}{F^2\beta_0 + \sigma_0^2 + \sigma_1^2}
\end{align}
この計算は帰納的に繰り返すことができ、条件付き分布の正規分布パラメータに関する漸化式が得られます。
\begin{align}
\alpha_{k+1} &= \frac{\sigma_1^2}{F^2\beta_k + \sigma_0^2 + \sigma_1^2}F\alpha_k + \frac{F^2\beta_k +\sigma_0^2}{F^2\beta_k + \sigma_0^2 + \sigma_1^2}y_{k+1}\\
\beta_{k+1} &= \frac{(F^2\beta_k +\sigma_0^2)\sigma_1^2}{F^2\beta_k + \sigma_0^2 + \sigma_1^2}
\end{align}
これはy_{k+1}が得られるたびに更新することができます。推定を行うことを考えると、\alpha_kは条件付き分布の期待値なので、ここでの求めたかった値そのものです。それとは別に\beta_kも計算可能であり、これによってその推定がどの程度信頼できるかを評価できます。
解釈
カルマンフィルタのお気持ちとして、
- 直前の予測からの、ノイズなしモデルによる時間発展(イノベーション)
- 現在の観測値がノイズなしだったと仮定したときの逆算
を、いい感じに合成する、という見方ができます。
実際\alpha_{k+1}は
\begin{align}
G_{k+1} &= \frac{F^2\beta_k +\sigma_0^2}{F^2\beta_k + \sigma_0^2 + \sigma_1^2}\\
\alpha_{k+1} &= (1-G_{k+1})F\alpha_k + G_{k+1}y_{k+1}
\end{align}
と書くことができ、F\alpha_kは、ノイズを無視した場合の過去の推定値からの発展、y_{k+1}は現在の観測値をノイズを無視して決め打ちしたときのx_{k+1}です。これの結合重みG_{k+1}が、いわゆるカルマンゲインと呼ばれるもので、過去の履歴や新たに加わるノイズ分散の大きさから、ちょうどよい結合重みになるように更新されます。
Discussion