TL;DR

同時分布がガウス分布に従うならば周辺分布もガウス分布に従う。
xに関する周辺分布$p(x)$と条件付き分布$p(y|x)$が与えられていれば
確率の法則等を使って、同時分布$p(x, y)$とyに関する周辺分布$p(y)$, xとyをひっくり返した
条件付き分布$p(x|y)$がどんなガウス分布に従うかも導出できる（まるでベイズの定理のように）

2.3.2 周辺ガウス分布

この節でいいたいことは

同時分布$ p(\mathbf{x} _a, \mathbf{x} _b) $がガウス分布であれば、下式周辺分布もガウス分布である

ということ。

$$
p( \mathbf{x}_a ) = \int p( \mathbf{x} _a, \mathbf{x} _b ) d \mathbf{x} _b \tag{2.83}
$$

そこで、条件付き分布と同様の指針で

($\mathbf{x} _b$は積分で消えるので)同時分布のexpの中身が$\mathbf{x} _a$の二次形式で書けることを示す ⇒ ガウス分布であることが示せた
$\mathbf{x} _a$の二次、一次の係数に注目して平均、分散を求める

方針で進める。

これから行う式変形の概要を先につかんでおこう。同時分布のexpの中身である式(2.70)から

積分消去したい$\mathbf{x} _b$に関する項を取り出して、平方完成する。すると、ガウス分布のexpの中身と同じ形の$\mathbf{x} _b$に関する項と$\mathbf{x} _a$に関する項に分けられる。
ガウス分布のexpの中身と同じ形の$\mathbf{x} _b$に関する項は、ガウス分布の積分が1になることを利用すれば$\mathbf{x} _a$にも$\mathbf{x} _b$にも依存しないことが分かるので、以後無視してよい。
結局、並行完成のときに出てきた$\mathbf{x} _a$に関する項と元々の(2.70)式から取り出した$\mathbf{x} _b$には依存しないが$\mathbf{x} _a$に関する項の2つについて考えればよい。
上記2つの項を考えると、$\mathbf{x} _a$の二次形式になっているので、ガウス分布になっていることが示せた！
さらに、$\mathbf{x} _a$の二次、一次の係数に注目して平均、分散を求まった！

Step1. 平方完成

式(2.70)から$\mathbf{x} _b$に関する項を取り出すと

$$
\frac{1}{2} x_b^T \Lambda_{bb} x_b + x_b^T m \tag{2.84}
$$
ただし、

$$
m = \Lambda_{bb} \mu_b - \Lambda_{ba} (x_a - \mu_a) \tag{2.85}
$$

と書ける。ここから式(2.84)を平方完成したい。ベクトルでなく、スカラーの以下の式を平方完成すると

\begin{align}

- \frac{1}{2} a x^2 + mx &= - \frac{a}{2} \left( x^2 - \frac{2m}{a} x \right) \\
&=  - \frac{1}{2} a \left( x - \frac{m}{a} \right)^2 - \frac{m^2}{a^2} \cdot - \frac{a}{2} \\
&= - \frac{1}{2} a \left( x - \frac{m}{a} \right)^2 + \frac{m^2}{2a} \\
&= - \frac{1}{2} a \left( x - a m^{-1} \right)^2 + \frac{1}{2} a^{-1} m^2 

\end{align}

となる。上式の$a$を$\Lambda_{bb}$に置き換えて、ベクトルである式(2.84)の平方完成を行うと以下のようになる（なりそうな感じが感じられればOKと思おう）。

$$

\frac{1}{2} ( x_b - \Lambda_{bb}^{-1} m )^T \Lambda_{bb} ( x_b - \Lambda_{bb}^{-1} m ) + \frac{1}{2} m^T \Lambda_{bb}^{-1} m \tag{2.84}
$$

(2.84)式をみると$ \mathbf{x}_b $に依存している第１項と$ \mathbf{x}_a $に依存している第２項に分けられた。

Step.2 平方完成で現れたx_bに依存している項

(2.84)式の第１項である$ - \frac{1}{2} ( x_b - \Lambda_{bb}^{-1} m )^T \Lambda_{bb} ( x_b - \Lambda_{bb}^{-1} m ) $は多次元ガウス分布のexpの中身である$ - \frac{1}{2} (x - \circ )^T 行列 (x - \circ )$ の形になっている！
なので、(2.84)式の第１項の積分には多次元ガウス分布の積分が1になるという定義式

$$
\int \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\} d \mathbf{x} = 1
$$

を変形して

$$
\int
\exp \left\{ - \frac{1}{2} (\mathbf{x} - \boldsymbol{\mu} )^T \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu} ) \right\} d \mathbf{x} =
\left( \frac{1}{ (2 \pi)^{D/2} } \frac{1}{| \Sigma |^{1/2}} \right)^{-1} \tag{A}
$$

となることが利用できる。式(A)に注目すると積分後の結果である右辺には（当然ながら）$x$が消えて共分散行列にのみ依存していることが分かる。そのため、(2.84)式の第１項の積分である

$$
\int \exp \left\{ - \frac{1}{2} ( x_b - \Lambda_{bb}^{-1} m )^T \Lambda_{bb} ( x_b - \Lambda_{bb}^{-1} m ) \right\} d x_b \tag{2.86}
$$

の結果は$x_a$, $x_b$にも依存せず、以後の式変形では無視してよいことが分かった。

Step.3 x_aの項

$ \mathbf{x} _a $に依存する項は$ \mathbf{x} _b $に関する積分では定数のように扱えるので以後は積分についての議論は無視して考える。

Step.2の平方完成で導出した$ \mathbf{x} _a $に依存する項である

$$\frac{1}{2} m^T \Lambda_{bb}^{-1} m = \frac{1}{2} [\Lambda_{bb} \mu_b - \Lambda_{ba} (x_a - \mu_a)]^T \Lambda_{bb}^{-1} [\Lambda_{bb} \mu_b - \Lambda_{ba} (x_a - \mu_a)] \tag{B}$$

と元々の(2.70)式から取り出した$\mathbf{x} _ b$には依存しないが$\mathbf{x} _a$に関する項について考えればよい。「PRML復活の呪文 part4」で使った線形代数のテクニックを使うと、$\mathbf{x} _ b$には依存しないが$\mathbf{x} _a$に関する項は

\begin{align}

& - \frac{1}{2} x_a^T \Lambda_{aa} x_a + 
\frac{1}{2} x_a^T \Lambda_{aa} \mu_a + \frac{1}{2} \mu_a^T \Lambda_{aa} x_a +
\frac{1}{2} x_a^T \Lambda_{ab} \mu_b + \frac{1}{2} \mu_b^T \Lambda_{ba} x_a \\

&=　- \frac{1}{2} x_a^T \Lambda_{aa} x_a +  x_a^T \Lambda_{aa} \mu_a + x_a^T \Lambda_{ab} \mu_b \\
&=　- \frac{1}{2} x_a^T \Lambda_{aa} x_a +  x_a^T \Lambda_{aa} \mu_a + x_a^T \Lambda_{ab} \mu_b \\
&=　- \frac{1}{2} x_a^T \Lambda_{aa} x_a +  x_a^T ( \Lambda_{aa} \mu_a + \Lambda_{ab} \mu_b )
\tag{C}

\end{align}

Step.4 x_aの項の式変形

式(B)と式(C)を合わせた項を変形すると以下の式が導出できる。

$$
\frac{1}{2} x_a^T ( \Lambda_{aa} - \Lambda_{ab}\Lambda_{bb}^{-1} \Lambda_{ba}) x_a + x_a^T ( \Lambda_{aa} - \Lambda_{ab}\Lambda_{bb}^{-1} \Lambda_{ba} ) \mu_a \tag{2.87}
$$

ここから式(2.83)によって周辺分布に変形してもexpの中身は$ \mathbf{x}_a $の2次形式となっているので、周辺分布もガウス分布であることが分かった！

Step.5 周辺分布の平均と分散

共分散

2.3.1節より、$ \mathbf{x}_a $の2次の項の係数は共分散行列の逆行列となっていたことから、
周辺分布$p( \mathbf{x} _ a )$の共分散は

$$
\Sigma_a = ( \Lambda_{aa} - \Lambda_{ab}\Lambda_{bb}^{-1} \Lambda_{ba})^{-1} \tag{2.88}
$$

ここでブロック行列表現の精度行列と共分散行列の関係である

\begin{align}

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

= 

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

\tag{2.90}

\end{align}

とシューア補行列を使うと

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

と書けるので結局は$ \Sigma_a = \Sigma_{aa} $であることが分かった。

平均

$ \mathbf{x}_a $の1次の項の係数は共分散行列の逆行列と平均の積となっていたことから、
周辺分布$p( \mathbf{x} _ a )$の平均$\hat{\mu _a}$ ¹は

$$
\Sigma_a^{-1} \hat{\mu_a} = ( \Lambda_{aa} - \Lambda_{ab}\Lambda_{bb}^{-1} \Lambda_{ba} ) \mu_a
$$

左から$\Sigma_a$、すなわち式(2.88)を掛けると

$$
\hat{\mu_a} = \mu_a
$$

と書けることが分かった。

2.3.3 ガウス変数に対するベイズの定理

この節では$p(x), p(y|x)$というガウス分布が与えられた際に、周辺ガウス分布$p(y)$と条件付きガウス分布$p(x|y)$がどうなるか求めよう。

2.3.1節, 2.3.2節で同時分布$p(x, y)$を式変形していくことで周辺分布と条件付き分布の性質を求めた。
そこで、本節でもまずは

同時分布を求めてから
$p(y), p(x|y)$を求める

指針となっている。
与えられる周辺分布$p(x)$は$\mu$を平均に関係するパラメータ、$\Lambda$を精度行列とおいて

$$
p(x) = N( x | \mu, \Lambda^{-1} ) \tag{2.99}
$$

とする。ここで2.3.1節で導出した条件付き分布の平均である

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

に着目すると、平均は観測値$\mathbf{x}_b $に関する一次式となっていることが分かる。
なので、与えられる条件付き分布$p(y|x)$は$A, b$を平均に関係するパラメータ、$L$を精度行列とおいて

$$
p(y | x) = N( y | Ax + b, L^{-1} ) \tag{2.100}
$$

とする。

同時分布を求める

\begin{align}

z = 

\begin{pmatrix}
x \\\
y
\end{pmatrix} \tag{2.101}

\end{align}

とおく。同時分布は$ p(x, y) = p(z) = p(x) p(y|x) $と書ける。また、ガウス分布のexpの外の部分は$x, y$に依存しないので無視してよい。対数をとって(2.99), (2.100)を使うと

\begin{align}

\ln p(z) &= \ln p(x) + \ln p(y|x) \\
&= -\frac{1}{2} (x - \mu)^T \Lambda (x - \mu) 
-\frac{1}{2} (y - Ax - b)^T L (y - Ax - b) \tag{2.102}

\end{align}

2.3.1節の議論より、$z$の2次の項の係数が同時分布の精度行列になっているので、式(2.102)を2次の項について整理すると、

\begin{align}

- \frac{1}{2}

\begin{pmatrix}
x \\\
y
\end{pmatrix}^T

\begin{pmatrix}
\Lambda + A^T L A  & - A^T L \\\
- LA               & L
\end{pmatrix}

\begin{pmatrix}
x \\\
y
\end{pmatrix}

= -\frac{1}{2} z R z \tag{2.103}

\end{align}

よって、同時分布の精度行列は

\begin{align}

R = 

\begin{pmatrix}
\Lambda + A^T L A  & - A^T L \\\
- LA               & L
\end{pmatrix} \tag{2.104}

\end{align}

となり、シューア補行列を使うことで精度行列の逆行列（共分散行列）も求めることができる。

\begin{align}

\text{cov} [z] = R^{-1} = 

\begin{pmatrix}
\Lambda^{-1}   & \Lambda^{-1} A^T \\\
A \Lambda^{-1} & L^{-1} + A \Lambda ^{-1} A^T
\end{pmatrix} \tag{2.105}

\end{align}

また、$z$の1次の項の係数から平均が求められるので、式(2.102)を1次の項について整理すると、

\begin{align}

\begin{pmatrix}
x \\\
y
\end{pmatrix}^T

\begin{pmatrix}
A \mu - A^T L b \\\
Lb
\end{pmatrix}
\tag{2.106}

\end{align}

よって、

\begin{align}

\mathrm{E}[z] = \mathrm{E} \left[
\begin{pmatrix}
x \\\
y
\end{pmatrix}

\right] = R^{-1}

\begin{pmatrix}
A \mu - A^T L b \\\
Lb
\end{pmatrix}

= 

\begin{pmatrix}
\mu \\\
A \mu + b
\end{pmatrix}
\tag{2.108}

\end{align}

周辺分布p(y)の平均と分散を求める

2.3.2節の議論を思い出そう。周辺分布の平均は$ \mathrm{E} [ \mathbf{x} _a ] = \mu _a \tag{2.92}$つまり（周辺化後の）変数$\mathbf{x} _a$に関する平均となる。
いま求めたいのは$\mathbf{x}$を周辺化した$\mathbf{y}$についての周辺分布の平均なので、式(2.108)から$\mathbf{y}$に関する平均を抜き出して

$$
\mathrm{E} [y] = A \mu + b \tag{2.109}
$$

となる。

次に周辺分布の分散は$ \text{cov} [ \mathbf{x} _a ] = \Sigma _ {aa} \tag{2.93} $
つまり（周辺化後の）変数$\mathbf{x} _a$のみについての共分散行列で書ける。
なので、式(2.105)から$\mathbf{y}$に関する共分散行列を抜き出して

$$
\text{cov} [ \mathbf{y} ] = L^{-1} + A \Lambda^{-1} A^T \tag{2.110}
$$

となることが分かった！

条件付き分布p(x|y)の平均と分散を求める

2.3.1節の議論より、条件付き分布の分散は

$$
\Sigma_{a|b} = \Lambda_{aa}^{-1} \tag{2.73}
$$

と精度行列を使えば簡単に書けるのであった。
いま求めたいのは$\Sigma_{x|y}$についても精度行列である式(2.104)から$\Lambda_{xx}$に相当する箇所を抜き出して

$$
\text{cov} [ \mathbf{x} | \mathbf{y} ] = ( \Lambda + A^T L A)^{-1} \tag{2.112}
$$

と書ける。条件付き分布の平均は式(2.75)より

$$
\mu _ {a|b} = \Sigma _ {a|b} { \Lambda _ {aa} \mu_a - \Lambda _ {ab} (x_b - \mu_b ) } \tag{2.75}
$$

と書ける。a, bをx, yに置き換えて、$\Sigma _ {x|y}$に相当する部分は式(2.112)で導出した。また、$\Lambda _ {xx}, \Lambda _ {xy}$は式(2.104)の精度行列から求まる。$\mu_x$は式(2.99)の定義から$\mu$になり、$x_b$は$y$, $\mu_b$は式(2.108)より$ A \mu + b$と求まった。これらを使って式変形を進める。

\begin{align}

\mathrm{E} [ \mathbf{x} | \mathbf{y} ] &= ( \Lambda + A^T L A)^{-1}
\left\{ ( \Lambda + A^T L A ) \mu + A^T L ( y - A \mu - b) \right\} \\

&= ( \Lambda + A^T L A)^{-1} ( \Lambda \mu + A^T L A \mu + A^T L y - A^T L A \mu - A^T L b ) \\

&= ( \Lambda + A^T L A)^{-1} \left\{ A^T L (y - b) + \Lambda \mu \right\} \tag{2.111}

\end{align}

とテキストの式を導出できた。

ベイズの定理？

ここまでの式変形で、$p(x), p(y|x)$が与えられたときに、条件である$x$と$y$を入れ替えた形の式を導出したりした。この操作は、あたかもベイズの定理を使って変数を入れ替えている例として見ることができる。

求めたい周辺分布の平均$\mu _ a$とこれまでの式変形で使っている$\mu_a$（つまり、式(2.65)で定義された$\mu_a$）が同じ記号で混乱するので、求めたい周辺分布の平均にハットをつけて表記した ↩