PRML復活の呪文 part4 (2.3 - 2.3.1)

### おことわり

PRMLの英語版がフリーで公開されていますので、必要に応じて英語版の図を引用させていただいてます。

# TL;DR

2つの変数集合$\mathbf{x}, \mathbf{y}$の同時分布$ p(\mathbf{x}, \mathbf{y}) $がガウス分布に従うならば、条件付き分布$ p(\mathbf{x} | \mathbf{y}) $や周辺分布$ p(\mathbf{x} ) $もガウス分布になる。

# 2.3 ガウス分布

連続変数の分布のモデルとしてよく用いられる。
1次元の場合は

$$
N(x | \mu, \sigma^2 ) = \frac{1}{ (2 \pi \sigma^2 )^{1/2} } 
\exp \left\\{ - \frac{1}{2 \sigma^2} ( x - \mu ) ^2 \right\\} \tag{2.42}
$$

D次元の場合は
$$
N( \mathbf{x} | \boldsymbol{\mu}, \Sigma ) = \frac{1}{ (2 \pi)^{D/2} } \frac{1}{| \Sigma |^{1/2}} 
\exp \left\\{ - \frac{1}{2} (\mathbf{x} - \boldsymbol{\mu} )^T \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu} )  \right\\} \tag{2.43}
$$

expの前にかかっているのは確率分布の足し合わせが1になるようにかけている正規化係数。
行列式$|\Sigma|$の解釈についてはpart2参照。ベクトルに行列$ \Sigma^{-1} $をかけると
ベクトルが回転し、長さが変わる。長さが変わると大きさが変わってしまうので、
その逆数の$ \frac{1}{| \Sigma |^{1/2}} $ をかけて大きさを1に戻しているイメージ。
参照：[part2, MAP推定のところ](https://qiita.com/yosshun/items/c31cf62fb4668b89c0e9#map%E6%8E%A8%E5%AE%9A)

### マハラノビス距離

多次元ガウス分布 - 式(2.43) - のexpの中身が何を意味しているのかを調べる。
結果から言うと以下の量はマハラノビス距離と呼ばれる。
$$
\Delta^2 = (\mathbf{x} - \boldsymbol{\mu} )^T \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu} )
\tag{2.44}
$$

(2.44)式はスカラーとなることに留意。
日常使う距離（ユークリッド距離）は等高線を青線で描くと、下図のように原点とした点からの等距離となる位置に線が引かれるので**真円**となる。
一方、マハラノビス距離の等高線は**楕円**となる。右上がり（左下から右上）の方向は距離のスケールが長いが、
右下がり（左上から右下）の方向の距離のスケールは短い。方向によって距離のスケールが変わる。
![キャプチャ.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/140135/5cb8e759-5bc2-4f49-f560-3dc7c80760ea.png)

(補足) 
　マハラノビス距離は異常検出などにも応用できる。例として、人の身長、体重の2次元データが大量にあるとして、そこから珍しい体系（異常値）のデータを見つけたいとしよう。2次元のデータをプロットすると、身長と体重は相関関係にあるので、下図①のようになるだろう。
　さて、②の点A, 点Bのどちらが珍しい体系と思うだろうか？多くの人は点Aと思うだろう。データの中心である黒点からユークリッド距離が同じである等高線を青線で引いた。点A, 点Bとも黒点からのユークリッド距離は等しいので、ユークリッド距離ではどちらが珍しいか判断できない。
　一方、マハラノビス距離で考えると③の等高線のようになる。データは右上がりの方向に分布しているので右上がり方向にの距離のスケールは長い（離れていても値は小さい）と考える。逆に、右下がりの方向は距離のスケールは短い（少しでも離れていれば値が大きくなる）と考える。この結果、黒点からのマハラノビス距離は点Aの方が点Bより遠くなり、珍しいデータであると判定できる。

![キャプチャ2.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/140135/d2ac72c2-23d0-15ba-0d7f-5905574afe59.png)

(補足ここまで)

　ここから、より詳細に式(2.44)の幾何的性質を調べよう。
共分散行列$\Sigma$はその定義から（実数の）対称行列であり[^1]、実数の対数行列の固有値は実数になる（線形代数の性質でそういうものだと考える。納得できない場合はぐぐろう）。
$\Sigma$の固有ベクトルを$\mathbf{u}_i$、固有値を$ \lambda _{i} $とおく。$\Sigma, \Sigma^{-1}$は固有値と固有ベクトルの積の足し合わせで書き直すことができる（線形代数のテクニックでそうなるのだと思おう）。

$$
\Sigma = \sum_{i=1}^D \lambda_i \mathbf{u}_{i} \mathbf{u} _{i}^T \tag{2.48}
$$

$$
\Sigma^{-1} = \sum_{i=1}^D \frac{1}{\lambda_i} \mathbf{u}_{i} \mathbf{u} _{i}^T \tag{2.49}
$$

(2.49)を(2.44)に代入して

```math
\begin{align}

\Delta^2 &= (\mathbf{x} - \boldsymbol{\mu} )^T \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu} ) \\
&= (\mathbf{x} - \boldsymbol{\mu} )^T \left( \sum_{i=1}^D \frac{1}{\lambda_i}  \mathbf{u}_{i} \mathbf{u} _{i}^T \right) (\mathbf{x} - \boldsymbol{\mu} ) \\
スカラー（\lambda）の部分を頭にもってきて \\
&= \sum_{i=1}^D \frac{1}{\lambda_i} (\mathbf{x} - \boldsymbol{\mu} )^T \mathbf{u}_{i} \mathbf{u} _{i}^T (\mathbf{x} - \boldsymbol{\mu} ) \tag{A}

\end{align}
```

ここで$ y_i =\mathbf{u}_i^T (\mathbf{x} - \boldsymbol{\mu} ) $と定義する。これはベクトル$\mathbf{u} _i$とベクトル$(\mathbf{x} - \boldsymbol{\mu} )$の内積であり、$ y$はスカラーであることに留意。内積はかける順番を変えても結果は同じなので、

$$
y_i =\mathbf{u}_i^T (\mathbf{x} - \boldsymbol{\mu} )
 = (\mathbf{x} - \boldsymbol{\mu} )^T \mathbf{u} _i \tag{2.51}
$$

よって式(A)は
$$
\Delta^2 = \sum_{i=1}^D \frac{y_i^2}{\lambda _i} \tag{2.50}
$$
と書ける。

内積(2.51)が下図のようにベクトル$ \mathbf{x} - \boldsymbol{\mu} $からベクトル$ \mathbf{u} $への写像であることを意識する。

![キャプチャ3.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/140135/6887cc1f-4ade-2728-c143-6c78f5755f4b.png)

すると式(2.50), (2.51)は$\mathbf{x}$を$\mathbf{ \mu }$だけ平行移動したベクトルを$\mathbf{u}$に写像したので、座標系は$\mathbf{u}$に変わっていると解釈できる。
さらに、$\mathbf{u} _1, \mathbf{u} _2$の距離のスケールは式(2.50)の$\lambda$でコントロールされているのでテキストの図2.7のようになり、距離の等高線を赤線で書くと  

- 座標系が$\mathbf{u}$になったのでx軸に対して傾いており
- 距離のスケールが$u_1, u_2$で異なるので楕円

になっている。

![キャプチャ4.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/140135/a64c4a67-da8e-2ff5-9c9b-0e7bdd9d652c.png)

[^1]: 共分散行列の(i,j)成分はi番目の要素とj番目の要素の共分散。一方、(j,i)成分はj番目の要素とi番目の要素の共分散。つまり、(i,j)成分と(j,i)成分は同じ値。


## 2.3.1 条件付きガウス分布

　ここから行列の計算がゴリゴリ続く。2.3.1節で言いたいことは
- 2つの多次元確率変数$ \mathbf{x} _a, \mathbf{x} _b$の同時分布$ p(\mathbf{x} _a, \mathbf{x} _b) $が多次元ガウス分布に従うなら、条件付き分布$ p(\mathbf{x} _a | \mathbf{x} _b) $あるいは$ p(\mathbf{x} _b | \mathbf{x} _a) $もガウス分布に従う、ということを計算で示す。
- 行列の計算の際、線形代数の性質により**式2.67のように1つの行列をいくつかのブロックに分割しても計算ができる**。
条件付き分布（ガウス分布）のパラメータである平均、分散が共分散行列$\Sigma$およびその逆行列$\Sigma^{-1}$のブロック行列（式2.67, 2.68）でどのように書けるかを計算する

### ブロック行列で計算？

具体例でまず理解しよう。4×4のサイズの行列A, Bの積を求めたいとしよう。このとき、行列A, Bを下図のように2×2のブロックに分割し、$A_{11}, \cdots, A_{22}$と置き直そう。分割された行列を**ブロック行列**とよぶ。

![キャプチャ7.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/140135/9c6446c1-4e49-35bc-d4fe-4d9455a17561.png)

すると行列ABの積は、あたかもブロック行列$ A_{11}, \cdots, B_{22} $がスカラーであるように計算でき、

```math
\begin{align}

AB = 

\begin{pmatrix}
A_{11} & A_{12} \\\
A_{21} & A_{22}
\end{pmatrix}

\begin{pmatrix}
B_{11} & B_{12} \\\
B_{21} & B_{22}
\end{pmatrix}

=

\begin{pmatrix}
A_{11} B_{11} + A_{12} B_{21} & A_{11} B_{12} + A_{12} B_{22} \\\
A_{21} B_{11} + A_{22} B_{21} & A_{21} B_{12} + A_{22} B_{22}
\end{pmatrix}

\end{align}
```
ここで左辺の行列の各要素である$A_* B_* + A_* B_*$は2×2の行列積の足し算なので、サイズは2×2。上下左右にサイズ2×2のブロック行列があるので、結合させるとサイズは8×8となっており、普通にABの積を求めたときとサイズが同じになっている。また、証明はしないが値も同じになる。

　このように行列を適当な位置でブロック行列に分割して計算するテクニックがあることを知っておこう。

### 条件付き分布が本当にガウス分布に従うのか？

　多次元ガウス分布（式2.43）の$\mathbf{x}$に依存している項はexpの中身の$
(\mathbf{x} - \boldsymbol{\mu} )^T \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu} )
$だけであり、ベクトルの2乗である$ \mathbf{x}^T \mathbf{x}$に$\boldsymbol{\mu}$やら$\Sigma^{-1}$やら係数がかかっているものとみなせる（二次形式という）。
これは（ベクトルでなく）一次元の変数でかくと$a x^2 + b x + c$と同じことを表している。
　なので、条件付き分布$ p(\mathbf{x} _a | \mathbf{x} _b) $がガウス分布に従うということを示すにはexpの中身が$\mathbf{x} _a$の二次形式で書けることを示せばよい！
(観測値$\mathbf{x} _b $を固定値＝定数とみなせる)
　
　計算にあたり、数学テクニックとしてブロック行列を使っていく。$\mathbf{x} _a$と$\mathbf{x} _b$をまとめたベクトル$\mathbf{x}$を考えて、このベクトルをブロック行列の計算により、$\mathbf{x} _a$の形の式に書き換える。ここで、$\Lambda$は共分散行列$\Sigma$の逆行列であり、精度行列 (precision matrix) と呼ばれる。こいつを定義しておくと、後々で式が簡単に書けるようになるので導入されている。
また、共分散行列の中のブロック $ \Sigma _ {aa} $ はM次元ベクトル$\mathbf{x} _a$どおしの共分散行列となっている。なので、行列の大きさは$ M \times M$。他の共分散行列ブロックも同じ考え。

![キャプチャ5.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/140135/beb56a63-5a7e-275c-8066-5e5229963709.png)

　さて、product ruleにより

```math
\begin{align}

p(x_a, x_b) &= p(x_a | x_b) p(x_b) \\
p(x_a | x_b) &= \frac{ p(x_a, x_b) }{ p(x_b) }

\end{align}
```
と書ける。条件付き分布$ p(x_a | x_b) $は観測値$x_b$を固定した$x_a$の関数とみなせる。なので、上式より条件付き分布は同時分布$ p(x_a , x_b) $に比例しているといえる。そこで、同時分布の式変形を行って、条件付き分布がガウス分布に従うことを示そう。
　同時分布$ p(x_a , x_b) $の多次元ガウス分布のexpの中身をブロック行列の考えで変形すると以下の式2.70が導出できる（納得できない場合は、以下の「補足」を見てみると腹落ちできるかも）

![キャプチャ6.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/140135/47151c0d-b337-07d9-9531-50908826eddc.png)

　expの中身が$\mathbf{x} _a$の二次形式で書けることが示せたので、条件付き分布もガウス分布に従うことが分かった！


### 条件付き分布のガウス分布の平均、分散パラメータはどうなるのか？

平均、分散を求めるための指針を確認しよう。
（条件付き分布でない）普通のガウス分布のexpの中身を$\mathbf{x}$について整理する。つまり、
$\mathbf{x}$についての2次式$\mathbf{x}^T (係数) \mathbf{x}$と
$\mathbf{x}$についての1次式$\mathbf{x}^T (係数)$と
$\mathbf{x}$についての0次式に書き直す。

$$
- \frac{1}{2} ( \mathbf{x} - \mathbf{\mu} )^T \Sigma^{-1} ( \mathbf{x} - \mathbf{\mu} ) = 
- \frac{1}{2} \mathbf{x}^T \Sigma^{-1} \mathbf{x} + \mathbf{x}^T \Sigma^{-1} \mu + 
( \mathbf{x} によらない項 ) \tag{2.71}
$$

右辺に注目すると
- $\mathbf{x}$の2次の項には共分散行列の逆行列と$ - \frac{1}{2} $ が係数としてかかっている
- $\mathbf{x}$の1次の項には共分散行列の逆行列と平均ベクトルがかかっている
ことが分かる。

なので観測値$\mathbf{x} _b $を固定値＝定数とみなして、(2.70)式を$\mathbf{x} _a $について整理して$\mathbf{x} _a $の2次の項、1次の項の係数に注目して平均、分散を求めよう！

#### 条件付き分布の分散

上記の(2.70)式の図より$\mathbf{x} _a $の2次の項を取り出すと

$$
- \frac{1}{2} \mathbf{x} _ a^T  \Lambda_{aa} \mathbf{x} _ a
$$
係数$\Lambda_{aa}$が共分散行列の逆行列になっているはずなので、条件付き分布の共分散行列は
$$
\Sigma_{a|b} = \Lambda_{aa}^{-1} \tag{2.72}
$$
と書ける。精度行列を使うことで非常に簡単に書くことができた！（これが精度行列を導入したメリット）

精度行列でなく、共分散行列で書くにはシューア補行列と呼ばれるブロック行列の逆行列を計算する公式を使う（証明は下記の演習2.24で）。

```math
\begin{align}

\begin{pmatrix}
A & B \\\
C & D
\end{pmatrix}^{-1}

&=

\begin{pmatrix}
M           & -MBD^{-1} \\\
-D^{-1} C M & D^{-1}+D^{-1} C M B D^{-1}
\end{pmatrix} \\ 

\tag{2.76}

\end{align}
```

$$
M = (A - B D^{-1} C )^{-1} \tag{2.77}
$$

これを使うと、

```math
\begin{align}

\Lambda_{aa} &= ( \Sigma_{aa} - \Sigma_{ab} \Sigma_{bb}^{-1} \Sigma_{ba} ) ^ {-1}

\end{align}
```
すなわち、

$$
\Lambda_{aa}^{-1} = ( \Sigma_{aa} - \Sigma_{ab} \Sigma_{bb}^{-1} \Sigma_{ba} )
$$

と書けるので、

```math
\begin{align}

\Sigma_{a|b} &=  \Sigma_{aa} - \Sigma_{ab} \Sigma_{bb}^{-1} \Sigma_{ba} \tag{2.82}

\end{align}
```
となる。

#### 条件付き分布の平均

上記の(2.70)式の図より$\mathbf{x} _a $の1次の項を取り出すと

```math
\begin{align}

\frac{1}{2} \mathbf{x}_a^T \Lambda_{aa} \mu_a 
+ \frac{1}{2} \mu_a^T \Lambda_{aa} \mathbf{x}_a
- \frac{1}{2} \mathbf{x}_a^T \Lambda_{ab} \mathbf{x}_b
+ \frac{1}{2} \mathbf{x}_a^T \Lambda_{ab} \mu_b
- \frac{1}{2} \mathbf{x}_b^T \Lambda_{ba} \mathbf{x}_a
+ \frac{1}{2} \mu_b^T \Lambda_{ba} \mathbf{x}_a \tag{A}

\end{align}
```

ここで、式変形に必要となる線形代数の公式を確認しよう。$a, b$がベクトル。$A$が対称行列とする。

- 対称行列の逆行列も対称行列（そういうものだと思おう）。$ (A^{-1})^T = A^{-1}$
- $ a^T A b = b^T A^T a$（なんでこうなるか気になったら、補足参照）
- さらに$A$が対称行列ならば$A^T = A$なので$　a^T A b = b^T A a $

$\Lambda_{ab} = \Lambda_{ba}^T$であることに注意すると、この公式より上式(A)が

```math
\begin{align}

\frac{1}{2} \mathbf{x}_a^T \Lambda_{aa} \mu_a 
+ \frac{1}{2} \mu_a^T \Lambda_{aa} \mathbf{x}_a
= \mathbf{x}_a^T \Lambda_{aa} \mu_a \\

- \frac{1}{2} \mathbf{x}_a^T \Lambda_{ab} \mathbf{x}_b
- \frac{1}{2} \mathbf{x}_b^T \Lambda_{ba} \mathbf{x}_a
= \mathbf{x}_a^T \Lambda_{ab} \mathbf{x}_b \\

\frac{1}{2} \mathbf{x}_a^T \Lambda_{ab} \mu_b
+ \frac{1}{2} \mu_b^T \Lambda_{ba} \mathbf{x}_a
= \mathbf{x}_a^T \Lambda_{ab} \mu_b \\

\end{align}
```

すなわち、

```math
\begin{align}

(A) &= \mathbf{x}_a^T \Lambda_{aa} \mu_a + \mathbf{x}_a^T \Lambda_{ab} \mathbf{x}_b 
+ \mathbf{x}_a^T \Lambda_{ab} \mu_b \\
&= \mathbf{x}_a^T ( \Lambda_{aa} \mu_a - \Lambda_{ab} ( \mathbf{x}_b - \mu_b ) )

\end{align}
```

となる。
$\mathbf{x} _ a$の係数である$( \Lambda_{aa} \mu_a - \Lambda_{ab} ( \mathbf{x} _ b - \mu_b ) )$が$ \Sigma_{a|b}^{-1} \mu_{a|b}$と等しくなるので

```math
\begin{align}

\Sigma_{a|b}^{-1} \mu_{a|b} = \Lambda_{aa}^{-1} \mu_{a|b} 
= \Lambda_{aa} \mu_a - \Lambda_{ab} ( \mathbf{x}_b - \mu_b ) \\

\mu_{a|b} = \mu_a - \Lambda_{aa}^{-1} \Lambda_{ab} ( \mathbf{x}_b - \mu_b ) \tag{2.75}
\end{align}
```

共分散を求めた時と同様、シューア補行列を使うと$\mu_{a|b}$をブロック共分散行列で書くこともできる

$$
\mu_{a|b} = \mu_a + \Sigma_{ab} \Sigma_{bb}^{-1} ( \mathbf{x}_b  - \mu_b ) \tag{2.81}
$$


# 補足（分かったつもりになるために）

### 式2.70のブロック行列の変形は本当に正しいの？

$(\mathbf{x} - \boldsymbol{\mu} )^T \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu} )$において、$ (\mathbf{x} - \boldsymbol{\mu} ) $はベクトル、$ \Sigma^{-1}  $は行列なので簡単化して2次元のベクトルと2×2の行列に置き換えて計算してみよう。


$$
x - \mu = 
\begin{pmatrix}
a \\\ 
b
\end{pmatrix}
, 
\Sigma^{-1} = 
\begin{pmatrix}
\Lambda_{aa} & \Lambda_{ab} \\\ 
\Lambda_{ba} & \Lambda_{bb}
\end{pmatrix}
$$
と置き直すと

```math
\begin{align}

\begin{pmatrix}
a & b
\end{pmatrix}

\begin{pmatrix}
\Lambda_{aa} & \Lambda_{ab} \\\ 
\Lambda_{ba} & \Lambda_{bb}
\end{pmatrix}

\begin{pmatrix}
a \\\
b
\end{pmatrix}

=

\begin{pmatrix}
a & b
\end{pmatrix}

\begin{pmatrix}
\Lambda_{aa} a + \Lambda_{ab} b \\\ 
\Lambda_{ba} a + \Lambda_{bb} b
\end{pmatrix}

= \Lambda_{aa} a^2 + \Lambda_{ab} ab + \Lambda_{ba} ab + \Lambda_{bb} b^2

\end{align}
```

式(2.70)は結局スカラーの計算ということを意識すれば、上式と(2.70)はよく似た形になっているので、なんか(2.70)の変形は成り立ちそうな気がしないだろうか...


### 以下の線形代数公式は本当に成り立つのか？[^2]

[^2]: 見出し内に数式が書けない

$ a^T A b = b^T A^T a $

2次元ベクトル、2×2サイズの行列で要素の積に実際に書き下して、$ a, b $について整理した式を比較してみると成り立ちそうな感じを納得できる。

```math
\begin{align}

a =
\begin{pmatrix}
a_1 \\\
a_2
\end{pmatrix}

, A = 
\begin{pmatrix}
A_{11} & A_{12} \\\ 
A_{21} &  A_{22}
\end{pmatrix}

, b =
\begin{pmatrix}
b_1 \\\
b_2
\end{pmatrix}

\end{align}
```
とおくと

```math
\begin{align}

a^T A b &= a_1 ( A_{11} b_1 + A_{12} b_2 ) + a_2 ( A_{21} b_1 + A_{22} b_2 ) \\\
&= b_1 ( A_{11} a_1 + A_{21} a_2 ) + b_2 ( A_{12} a_1 + A_{22} a_2 )

\end{align}
```

となる。$a$と$b$を入れ替えて式が書けること、対応する行列の要素が転置になっている（$ A_{12} $が $A_{21} $になっていたり）ことを見ると成り立ちそうな気がしないだろうか...

## 演習2.24

式(2.76)の左辺に左から

```math
\begin{align}

\begin{pmatrix}
A & B \\\
C & D
\end{pmatrix}

\end{align}
```
をかければ左辺は(ABCDっていう行列)×(その行列の逆行列)なので定義通り単位行列となる。
なので右辺の左からもこの行列をかけて、その積が単位行列となることを示せばよい。

```math
\begin{align}

\begin{pmatrix}
A & B \\\
C & D
\end{pmatrix}

\begin{pmatrix}
M           & -MBD^{-1} \\\
-D^{-1} C M & D^{-1}+D^{-1} C M B D^{-1}
\end{pmatrix} 

\\\

= 
\begin{pmatrix}
AM - BD^{-1}CM         & -AMBD^{-1} + BD^{-1} + BD^{-1}CMBD^{-1} \\\
CM - CM                & -CMBD^{-1} + CMBD^{-1} + I
\end{pmatrix} 

\\\

= 
\begin{pmatrix}
(A - BD^{-1}C) M & (AM - BD^{-1}CM ) (- BD^{-1} ) + BD^{-1} \\\
0                & I
\end{pmatrix} \tag{B}

\end{align}
```

ここで行列Mの定義である式(2.77)より
$$
(A - BD^{-1}C) M = (A - BD^{-1}C) (A - BD^{-1}C)^{-1} = I
$$

$$
(AM - BD^{-1}CM ) (- BD^{-1} ) + BD^{-1} = I (- BD^{-1} ) + BD^{-1} = 0
$$

よって式(B)は

```math
\begin{align}

\begin{pmatrix}
I & 0 \\\
0 & I
\end{pmatrix}

\end{align}
```

となり、単位行列となることが示せたので式(2.76)は成り立つ。