基礎から分かる時系列分析 2.カルマンフィルタ (数理的な導出)

## 2.1カルマンフィルタの適用条件
シミュレーションと観測演算子が共に行列で与えられた状態空間モデル

```math
\begin{align}
{\bf x}_t=F_t{\bf x}_{t-1}+{\bf v}_t,\ \ \ \ \ {\bf v}_t\sim \mathcal{N}({\bf 0},Q_t)\\
{\bf y}_t=H_t{\bf x}_t+{\bf w}_t,\ \ \ \ \ \ {\bf w}_t\sim \mathcal{N}({\bf 0},R_t)
\end{align}
```

において、全ての確率変数が正規分布に従う場合(線形ガウス状態空間モデル)、状態を逐次推定するカルマンフィルタを適用することができる。

## 2.2一期先予測
昨日までのデータから今日の状態を推定することが一期先予測である。これは

```math
\begin{align}
p({\bf x}_t|{\bf y}_{1:t-1})&=\int p({\bf x}_t,{\bf x}_{t-1}|{\bf y}_{1:t-1})d{\bf x}_{t-1}\\
&=\int p({\bf x}_{t-1}|{\bf y}_{1:t-1})p({\bf x}_t|{\bf x}_{t-1})d{\bf x}_{t-1}\tag{1}
\end{align}
```

と書ける。(1)に対応するアルゴリズムを導出したい。今、k期までのデータが与えられた時のj期の状態が従う分布を

```math
p({\bf x}_j|{\bf y}_{1:k})\equiv\mathcal{N}({\bf x}_{j|k},V_{j|k})
```

と表記する。すると(1)のうち

```math
p({\bf x}_{t-1}|{\bf y}_{1:t-1})=\mathcal{N}({\bf x}_{t-1|t-1},V_{t-1|t-1})
```

は既知(一期前にフィルタ分布によって導出ずみ)であるから、

```math
p({\bf x}_t|{\bf　x}_{t-1})
```

を求めれば良い。これを${\bf v}_t$に関する周辺化積分形に展開すると、

```math
p({\bf x}_t|{\bf x}_{t-1})=\int p({\bf x}_t|{\bf x}_{t-1},{\bf v}_t)p({\bf v}_t)d{\bf v}_t
```

システムモデルより、x_{t-1}とv_tの実現値が与えられればx_tは唯一の確定値$F_t{\bf x}_{t-1}+{\bf v}_t$を取るので、

```math
p({\bf x}_t|{\bf x}_{t-1},{\bf v}_t)=\delta ({\bf x}_t-F_t{\bf x}_{t-1}-{\bf v}_t)
```

今、正規分布の分散共分散行列を$\epsilon^2 I$とすると、

```math
p({\bf x}_t|{\bf x}_{t-1},{\bf v}_t)=\mathcal{N}(F_t{\bf x}_{t-1}+{\bf v}_t,\epsilon^2I),\ \ \ \ \ as\ \ \ \epsilon\rightarrow 0
```

すなわち

```math
p({\bf x}_t-F_t{\bf x}_{t-1}|{\bf v}_t)=\mathcal{N}({\bf v}_t,\epsilon^2I)
```

これを周辺化することができれば、x_{t-1}を定数扱いして線形変換することでp(x_t|x_{t-1})も求まる。

(2)をよく見ると、条件付き分布の期待値が条件付けている変数${\bf v}_t$に依存している。そこで補題Aを用いれば

```math
p({\bf x}_t-F_t{\bf x}_{t-1})=\mathcal{N}({\bf 0},Q_t+\epsilon^2I)
```

すなわち

```math
p({\bf x}_t|{\bf x}_{t-1})=\mathcal{N}(F_t{\bf x}_{t-1},Q_t+\epsilon^2I)
```

さらに補題Aを

```math
\begin{align}
p({\bf x}_{t-1}|{\bf y}_{1:t-1})&=\mathcal{N}({\bf x}_{t-1|t-1},V_{t-1|t-1})\\
p({\bf x}_t|{\bf x}_{t-1})&=\mathcal{N}(F_t{\bf x}_{t-1},Q_t+\epsilon^2I)
\end{align}
```

に適用すれば、

```math
p({\bf x}_t|{\bf y}_{1:t-1})=\mathcal{N}(F_t{\bf x}_{t-1|t-1},F_tV_{t-1|t-1}\ ^TF_t+Q_t)
```

となる。

## 2.2 フィルタ

$t$期までのデータから$t$期の状態を推定するのがフィルタ分布である。これは

```math
\begin{align}
p({\bf x}_t|{\bf y}_{1:t})&=\frac{p({\bf x}_t,{\bf y}_{1:t})}{p({\bf y}_{1:t})}\\
&= \frac{p({\bf x}_t,{\bf y}_{1:t})}{\int p({\bf x}_t,{\bf y}_{1:t})d{\bf x}_t}\\
&=\frac{p({\bf y}_t|{\bf x}_t)p({\bf x}_t|{\bf y}_{1:t-1})}{\int p({\bf y}_t|{\bf x}_t)p({\bf x}_t|{\bf y}_{1:t-1})d{\bf x}_t}\tag{2}
\end{align}
```

と書ける。(2)に対応するアルゴリズムを導出したい。観測モデルより一時点尤度は

```math
\begin{align}
p({\bf y}_t|{\bf x}_t)=\mathcal{N}(H_t{\bf x}_t,R_t)
\end{align}
```

である。また

```math
p({\bf x}_t|{\bf y}_{1:t-1})=\mathcal{N}({\bf x}_{t|t-1},V_{t|t-1})
```

は一期先予測より既知なので、補題Bを適用すれば

```math
p({\bf x}_t|{\bf y}_{1:t})=\mathcal{N}({\bf x}_{t|t},V_{t|t})
```
が得られ、

```math
\begin{align}
{\bf x}_{t|t}&={\bf x}_{t|t-1}+K_t({\bf y}_t-H_t{\bf x}_{t|t-1})\nonumber\\
V_{t|t}&=V_{t|t-1}K_tH_tV_{t|t-1}\nonumber\\
K_t&=V_{t|t-1}\ ^TH_t(H_tV_{t|t-1}\ ^TH_t+R_t)^{-1}
\end{align}
```

となる。$K_t$をカルマンゲインと呼ぶ。

### 補題A,B

補題A：

```math
p({\bf z})=\mathcal{N}({\bf \mu}_{{\bf z}},\Sigma_{{\bf z}}),\ p({\bf x}|{\bf z})=\mathcal{N}(F{\bf z},\Sigma_{{\bf x}|{\bf z}})
```
ならば、

```math
p({\bf x})=\mathcal{N}(F{\bf \mu}_{{\bf z}},F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})
```

補題B：

```math
p({\bf z})=\mathcal{N}({\bf \mu}_{{\bf z}},\Sigma_{\bf {z}}),\ p({\bf x}|{\bf z})=\mathcal{N}(F{\bf z},\Sigma_{{\bf x}|{\bf z}})
```

ならば、

```math
\begin{align}
p({\bf z}|{\bf x})&=\mathcal{N}({\bf \mu}_{{\bf z}|{\bf x}},\Sigma_{{\bf z}|{\bf x}})\\
{\bf \mu}_{{\bf z}|{\bf x}}&={\bf \mu}_{{\bf z}}+\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}({\bf x}-F{\bf \mu})\\
\Sigma_{{\bf z}|{\bf x}}&=\Sigma_{{\bf z}}-\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}F\Sigma_{{\bf z}}\nonumber
\end{align}
```

(証明)

補題A、Bを証明する。

```math
\begin{eqnarray}
p({\bf z})=\mathcal{N}({\bf \mu}_{{\bf z}},\Lambda_{{\bf z}}^{-1}),\ p({\bf x}|{\bf z})=\mathcal{N}(F{\bf z},\Lambda_{{\bf x}|{\bf z}}^{-1})\ \ \ \ \ \Lambda_{{\bf z}}=\Sigma_{{\bf z}}^{-1},\ \Lambda_{{\bf x}|{\bf z}}=\Sigma_{{\bf x}|{\bf z}}^{-1}
\end{eqnarray}
```

なので、まずは

```math
{\bf y}=\begin{pmatrix}{\bf z}\\{\bf x}\end{pmatrix}
```
と置いて、z,xの同時分布$p({\bf y})$を求める。

```math
\begin{align}
\log p({\bf y})&=\ ^T({\bf z}-{\bf \mu}_{{\bf z}})\Lambda_{{\bf z}}({\bf z}-{\bf \mu}_{{\bf z}})+\ ^T({\bf x}-F{\bf z})\Lambda_{{\bf x}|{\bf z}}({\bf x}-F{\bf z})+const\\
&&\\
&=\ ^T{\bf z}\Lambda_{{\bf z}}{\bf z}+\ ^T{\bf z}\ ^TF\Lambda_{\bf {x}|{\bf z}}F{\bf z}-\ ^T{\bf z}\ ^TF\Lambda_{{\bf x}|{\bf z}}{\bf x}-\ ^T{\bf x}\Lambda_{\bf {x}|{\bf z}}F{\bf z}+\ ^T{\bf x}\Lambda_{{\bf x}|{\bf z}}{\bf x}\\
&\ \ -\ ^T{\bf z}\Lambda_{{\bf z}}{\bf \mu}_{{\bf z}}-\ ^T{\bf \mu}_{{\bf z}}\Lambda_{{\bf z}}{\bf z}+const\\
&&\\
&=\begin{pmatrix}{\bf z} & {\bf x}\end{pmatrix}\begin{pmatrix}
\Lambda_{{\bf z}}+\ ^TF\Lambda_{{\bf x}|{\bf z}}F & -\ ^TF\Lambda_{{\bf x}|{\bf z}}\\
-\Lambda_{{\bf x}|{\bf z}}F &\Lambda_{{\bf x}|{\bf z}}
\end{pmatrix}\begin{pmatrix}{\bf z} \\ {\bf x}\end{pmatrix}\\
&\ \ -\begin{pmatrix}{\bf z}&{\bf x}\end{pmatrix}\begin{pmatrix}\Lambda_{{\bf z}}{\bf \mu}_{{\bf z}}\\ {\bf 0}\end{pmatrix}-\begin{pmatrix}\ ^T{\bf \mu}_{{\bf z}}\Lambda &{\bf 0}\end{pmatrix}\begin{pmatrix}{\bf z}\\{\bf 0}\end{pmatrix}+const\\
&&\\
&=\ ^T{\bf y}\Lambda_{{\bf y}}{\bf y}-\ ^T{\bf y}A-\ ^TA{\bf y}+const\\
&&\\
&=\ ^T({\bf y}-\Lambda_{{\bf y}}^{-1}A)\Lambda_{{\bf y}}({\bf y}-\Lambda_{{\bf y}}A)+const
\end{align}
```

ただし、

```math
\begin{align}
\Lambda_{{\bf y}}=\begin{pmatrix}
\Lambda_{{\bf z}}+\ ^TF\Lambda_{{\bf x}|{\bf z}}F & -\ ^TF\Lambda_{{\bf x}|{\bf z}}\\
-\Lambda_{{\bf x}|{\bf z}}F &\Lambda_{{\bf x}|{\bf z}}
\end{pmatrix},\ \ A=\begin{pmatrix}\Lambda_{{\bf z}}{\bf \mu}_{{\bf z}}\\ {\bf 0}\end{pmatrix}
\end{align}
```

と置いた($\Lambda_{{\bf y}}$は同時分布の精度行列)。したがって、同時分布の分散共分散行列と平均ベクトルは、

```math
\begin{align}
\Sigma_{{\bf y}}&=\Lambda_{{\bf y}}^{-1}=\begin{pmatrix}
\Sigma_{{\bf z}} & \Sigma_{{\bf z}}\ ^TF\\
F\Sigma_{{\bf z}} & \Sigma_{{\bf x}|{\bf z}}+F\Sigma_{{\bf z}}\ ^TF
\end{pmatrix}\\
{\bf \mu}_{{\bf y}}&=\Lambda_{{\bf y}}^{-1}A=\begin{pmatrix}{\bf \mu}_{{\bf z}}\\ F{\bf \mu}_{{\bf z}}\end{pmatrix}
\end{align}
```

補題Cにより、${\bf x}$の周辺分布は

```math
p({\bf x})=\mathcal{N}({\bf \mu}_{{\bf z}},\Sigma_{{\bf x}|{\bf z}}+F\Sigma_{{\bf z}}\ ^TF)
```

なので、補題Aは示された。

また補題Cにより、xが与えられた時zは正規分布に従い、その分散は

```math
\begin{align}
\Sigma_{{\bf z}|{\bf x}}=(\Lambda_{{\bf z}}+\ ^TF\Lambda_{{\bf x}|{\bf z}}F)^{-1}=\Sigma_{{\bf z}}-\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}F\Sigma_{{\bf z}}
\end{align}
```

平均は

```math
\begin{align}
{\bf \mu}_{{\bf z}|{\bf x}}={\bf \mu}_{{\bf z}}+\{\Sigma_{{\bf z}}-\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}F\Sigma_{{\bf z}}\}\ ^TF\Sigma_{{\bf x}|{\bf z}}^{-1}({\bf x}-F{\bf \mu}_{{\bf z}})
\end{align}
```

ここで、

```math
\begin{align}
&\ \ \Sigma_{{\bf z}}-\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}F\Sigma_{{\bf z}}\ ^TF\Sigma_{{\bf x}|{\bf z}}^{-1}\\
&= \Sigma_{{\bf z}}\ ^TF\{\Sigma_{{\bf x}|{\bf z}}^{-1}+(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}F\Sigma_{{\bf z}}\ ^TF\Sigma_{{\bf x}|{\bf z}}\}\\
&= \Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}\{(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})\Sigma_{{\bf x}|{\bf z}}^{-1}-F\Sigma_{{\bf z}}\ ^TF\Sigma_{{\bf x}|{\bf z}}^{-1}\}\\
&=\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}
\end{align}
```

なので、

```math
\begin{align}
{\bf \mu}_{{\bf z}|{\bf x}}={\bf \mu}_{{\bf z}}+\Sigma_{{\bf z}}\ ^TF(F\Sigma_{{\bf z}}\ ^TF+\Sigma_{{\bf x}|{\bf z}})^{-1}({\bf x}-F{\bf \mu})
\end{align}
```

以上より、補題Bが示された。


### 補題C

x_a,x_bの同時分布が

```math
\begin{align}
p(\begin{pmatrix}{\bf x}_a\\ {\bf x}_b\end{pmatrix})&=\mathcal{N}({\bf \mu},\Sigma)\\
{\bf \mu}&=\begin{pmatrix}{\bf \mu}_a\\ {\bf \mu}_b\end{pmatrix}\\
\Sigma&=\begin{pmatrix}\Sigma_{aa} & \Sigma_{ab}\\ \Sigma_{ba} &\Sigma_{bb}\end{pmatrix}\\
\Lambda=\Sigma^{-1}&=\begin{pmatrix}\Lambda_{aa} & \Lambda_{ab}\\ \Lambda_{ba} &\Lambda_{bb}\end{pmatrix}
\end{align}
```

で与えられているとき、x_aの周辺分布は

```math
\begin{align}
p({\bf x}_a)=\mathcal{N}({\bf \mu}_a,\Sigma_{aa})
\end{align}
```

x_b$が与えられた時のx_aの条件付き分布は

```math
\begin{align}
p({\bf x}_a|{\bf x}_b)&=\mathcal{N}({\bf \mu}_{a|b},\Lambda_{aa}^{-1})\\
{\bf \mu}_{a|b}&={\bf \mu}_a-\Lambda_{aa}^{-1}\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)
\end{align}
```

(証明)

補題Cを証明する。まずは周辺分布p(x_a)を求める。

```math
\begin{align}
p({\bf x}_a)&=\int p({\bf x}_a,{\bf x}_b)d{\bf x}_b\\
&\propto\int\exp\{-\frac{1}{2}\ ^T({\bf x}-{\bf \mu})\Lambda({\bf x}-{\bf \mu})\}d{\bf x}_b
\end{align}
```

指数関数の中身をx_bで平方完成する。

```math
\begin{align}
&\ \ ^T({\bf x}-{\bf \mu})\Lambda({\bf x}-{\bf \mu})\\
&&\\
&=\ ^T({\bf x}_a-{\bf \mu}_a)\Lambda_{aa}({\bf x}_a-{\bf \mu}_a)+\ ^T({\bf x}_a-{\bf \mu}_b)\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)\\
&\ \ +\ ^T({\bf x}_b-{\bf \mu}_b)\Lambda_{ba}({\bf x}_a-{\bf \mu}_a)+\ ^T({\bf x}_b-{\bf \mu}_b)\Lambda_{bb}({\bf x}_b-{\bf \mu}_b)\\
&&\\
&=\ ^T{\bf x}_b\Lambda_{bb}{\bf x}_b+\ ^T{\bf x}_b\Lambda_{ba}({\bf x}_a-{\bf \mu}_a)\\
&\ -\ ^T{\bf x}_b\Lambda_{bb}{\bf \mu}_b-\ ^T{\bf \mu}_b\Lambda_{bb}{\bf x}_b+\ ^T({\bf x}_a-{\bf \mu}_a)\Lambda_{ab}{\bf x}_b\\
&\ \ +\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf \mu}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}_a\Lambda_{aa}{\bf \mu}_a-\ ^T{\bf x}_a\Lambda_{ab}{\bf \mu}_b-\ ^T{\bf \mu}_b\Lambda_{ba}{\bf x}_a+const\\
&&\\
&=\ ^T{\bf x}_b\Lambda_{bb}{\bf x}_b-\ ^T{\bf x}_b\{\Lambda_{bb}{\bf \mu}_b-\Lambda_{ba}({\bf x}_a-{\bf \mu}_a)\}-\{\ ^T{\bf \mu}_b\Lambda_{bb}-\ ^T({\bf x}_a-{\bf \mu}_a)\Lambda_{ab}\}{\bf x}_b\\
&\ \ +\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}(\Lambda_{aa}{\bf \mu}_a+\Lambda_{ab}{\bf \mu}_b)-(\ ^T{\bf \mu}_a\Lambda_{aa}+\ ^T{\bf \mu}_b\Lambda_{ba}){\bf x}_a+const\\
&&\\
&=\ \ ^T{\bf x}_bA^{-1}{\bf x}_b-\ ^T{\bf x}_bB-\ ^TB{\bf x}_b\\
&\ \ +\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}(\Lambda_{aa}{\bf \mu}_a+\Lambda_{ab}{\bf \mu}_b)-(\ ^T{\bf \mu}_a\Lambda_{aa}+\ ^T{\bf \mu}_b\Lambda_{ba}){\bf x}_a+const\\
&&\\
&=\ ^T({\bf x}_b-AB)A^{-1}({\bf x}-AB)-\ ^TB\ ^TAB\\
&\ \ +\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}(\Lambda_{aa}{\bf \mu}_a+\Lambda_{ab}{\bf \mu}_b)-(\ ^T{\bf \mu}_a\Lambda_{aa}+\ ^T{\bf \mu}_b\Lambda_{ba}){\bf x}_a+const
\end{align}
```

ただし、

```math
\begin{align}
A=\Lambda_{bb}^{-1},\ \ B=\Lambda_{bb}{\bf \mu}_b-\Lambda_{ba}({\bf x}_a-{\bf \mu}_a)
\end{align}
```

と置いた。ゆえ、

```math
\begin{align}
p({\bf x}_a)&\propto\exp[\\
&\ \ -\frac{1}{2}\{-\ ^TB\ ^TAB+\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}(\Lambda_{aa}{\bf \mu}_a+\Lambda_{ab}{\bf \mu}_b)-(\ ^T{\bf \mu}_a\Lambda_{aa}+\ ^T{\bf \mu}_b\Lambda_{ba}){\bf x}_a\}]
\end{align}
```

右辺の指数部分の中身を変形してまとめると、

```math
\begin{align}
&\ \ ^T{\bf x}_a(-\ ^T\Lambda_{ba}\Lambda_{bb}^{-1}\Lambda_{ba}+\Lambda_{aa}){\bf x}_a\\
&\ \ -\ ^T{\bf x}_a(-\ ^T\Lambda_{ba}{\bf \mu}_b-\ ^T\Lambda_{ba}\Lambda_{bb}^{-1}\Lambda_{ba}{\bf \mu}_a+\Lambda_{aa}{\bf \mu}_a+\Lambda_{ab}{\bf \mu}_b)\\
&\ \ -(\ ^T{\bf \mu}_b\Lambda_{ba}{\bf x}_a-\ ^T{\bf \mu}_a\ ^T\Lambda_{ba}\Lambda_{bb}^{-1}\Lambda_{ba}+\ ^T{\bf \mu}_a\Lambda_{aa}+\ ^T{\bf \mu}_b\Lambda_{ab}){\bf x}_a\\
&=\ ^T{\bf x}_aB^{-1}{\bf x}_a-\ ^T{\bf x}_aC-\ ^C{\bf x}_a\\
&=\ ^T({\bf x}_a-BC)B^{-1}({\bf x}_a-BC)+const
\end{align}
```

ただし、

```math
\begin{align}
B&=(\Lambda_{aa}-\Lambda_{ab}\Lambda_{bb}^{-1}\Lambda_{ba})^{-1}\\
C&=(\Lambda_{aa}-\Lambda_{ab}\Lambda_{bb}^{-1}\Lambda_{ba}){\bf \mu}_a
\end{align}
```

したがって、

```math
\begin{align}
{\bf \mu}_a&=BC={\bf \mu}_a\\
\Sigma_a&=B=\Sigma_{aa}
\end{align}
```

つぎに条件付き分布p(x_a|x_b)を求める。

```math
\begin{align}
^T({\bf x}-{\bf \mu})\Sigma^{-1}({\bf x}-{\bf \mu})&=
^T\begin{pmatrix} {\bf x}_a-{\bf \mu}_a \\ {\bf x}_b-{\bf \mu}_b\end{pmatrix}
\begin{pmatrix}\Lambda_{aa} & \Lambda_{ab}\\ \Lambda_{ba} & \Lambda_{bb}\end{pmatrix}
\begin{pmatrix}{\bf x}_a-{\bf \mu}_a\\ {\bf x}_b-{\bf \mu}_b\end{pmatrix}\\
&&\\
&=\ ^T({\bf x}_a-{\bf \mu}_a)\Lambda_{aa}({\bf x}_a-{\bf \mu}_a)+\ ^T({\bf x}_a-{\bf \mu}_b)\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)\\
&\ \ +\ ^T({\bf x}_b-{\bf \mu}_b)\Lambda_{ba}({\bf x}_a-{\bf \mu}_a)+\ ^T({\bf x}_b-{\bf \mu}_b)\Lambda_{bb}({\bf x}_b-{\bf \mu}_b)\\
&&\\
&=\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}_a\Lambda_{aa}{\bf \mu}_a+\ ^T{\bf x}_a\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)\\
&\ \ -\ ^T{\bf \mu}_a\Lambda_{aa}{\bf x}_a+\ ^T({\bf x}_b-{\bf \mu}_b)\Lambda_{ba}{\bf x}_a+const\\
&&\\
&=\ ^T{\bf x}_a\Lambda_{aa}{\bf x}_a-\ ^T{\bf x}\{\Lambda_{aa}{\bf \mu}_a-\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)\}\\
&\ \ -\{\ ^T{\bf \mu}_a\Lambda_{aa}-\ ^T({\bf x}_b-{\bf \mu}_b)\Lambda_{ba}\}{\bf x}_a+const\\
&&\\
&=\ ^T{\bf x}_aA^{-1}{\bf x}_a-\ ^T{\bf x}_aB-\ ^TB{\bf x}_a+const\\
&&\\
&=\ ^T({\bf x}_a-AB)A^{-1}({\bf x}_a-AB)+const
\end{align}
```

ただし、

```math
\begin{align}
A=\Lambda_{aa}^{-1},\ B=\Lambda_{aa}{\bf \mu}_a-\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)
\end{align}
```

と置いた。したがって、条件付き分布p(x_a|x_b)の平均と分散は、

```math
\begin{align}
{\bf \mu_{a|b}}&=&\Lambda_{aa}^{-1}\{\Lambda_{aa}{\bf \mu}_a-\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)\}\\
&={\bf \mu}_a-\Lambda_{aa}^{-1}\Lambda_{ab}({\bf x}_b-{\bf \mu}_b)\\
\Sigma_{a|b}&=\Lambda_{aa}^{-1}
\end{align}
```

となる。(証明終わり)

つぎに$\Lambda_{aa}^{-1},\Lambda_{ab}$を求める(このままの方が簡潔なので普通は求めない)。今、

```math
\begin{align}
\begin{pmatrix}\Sigma_{aa} &\Sigma_{ab}\\ \Sigma_{ba} & \Sigma_{bb}\end{pmatrix}^{-1}=\begin{pmatrix}\Lambda_{aa} &\Lambda_{ab}\\ \Lambda_{ba} & \Lambda_{bb}\end{pmatrix}
\end{align}
```

であるから、

```math
\begin{align}
\Lambda_{aa}&=(\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba})^{-1}\\
\Lambda_{ab}&=-(\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba})^{-1}\Sigma_{ab}\Sigma_{bb}^{-1}
\end{align}
```

より、

```math
\begin{align}
{\bf \mu}_{a|b}&={\bf \mu}_a+\Sigma_{ab}\Sigma_{bb}^{-1}({\bf x}_b-{\bf \mu}_b)\\
\Sigma_{a|b}&=\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}
\end{align}
```

となった。条件付き分布の平均はx_bの線形関数。


## 2.3 尤度

尤度$p({\bf y}_{1:T})$は

```math
\begin{align}
p({\bf y}_{1:T})=\Pi_{t=1}^Tp({\bf y}_{t}|{\bf y}_{1:t-1})
\end{align}
```

ここで、

```math
\begin{eqnarray}
p({\bf x}_t|{\bf y}_{1:t-1})&=\mathcal{N}({\bf x}_{t|t-1}|V_{t|t-1})\\
p({\bf y}_t|{\bf x}_t)&=\mathcal{N}(H_t{\bf x}_t,R_t)
\end{eqnarray}
```

を補題Aに適用して、

```math
\begin{align}
p({\bf y}_t|{\bf y}_{1:t-1})=\mathcal{N}(H_t{\bf x}_{t|t-1},H_tV_{t|t-1}\ ^TH_t+R_t)
\end{align}
```

となる。したがって、y_tの次元数をmと置くと対数尤度は

```math
\begin{align}
\log p({\bf y}_{1:T})&=\sum_{t=1}^T\log p({\bf y}_t|{\bf y}_{1:t-1})\\
&=-\frac{1}{2}\sum_{t=1}^T\log|D_{t|t-1}|-\frac{Tm}{2}\log(2\pi)\\
&\ \ -\frac{1}{2}\sum_{t=1}^T\ ^T({\bf y}_t-H_t{\bf x}_{t|t-1})D_{t|t-1}^{-1}({\bf y}_t-H_t{\bf x}_{t|t-1})
\end{align}
```

ただし、

```math
\begin{align}
D_{t|t-1}=H_tV_{t|t-1}\ ^TH_t+R_t
\end{align}
```

と置いた。パラメータの次元数を$M$と置くと、AICは

```math
\begin{align}
AIC=\log p({\bf y}_{1:T})-M
\end{align}
```

なので、これをもとにモデル選択を行う。

[次回：基礎から分かる時系列分析 3.カルマンフィルタ(ローカルレベルモデルの実装)](https://qiita.com/_R2_/items/d65957bcfa8c59248c33)