はじめに
PRML解答例まとめを参照
演習 2.21
大きさがD\times Dの実対称行列の独立なパラメータは、D(D+1)/2個であることを示せ．
大きさがD\times Dの実対称行列の全成分の個数は当然D^2個である。このうち、対角成分のD個を除いて残りのパラメータ（非対角成分）は対角成分に対して対称な値になっていなければならないので、そのパラメータの自由度は\displaystyle \frac{D^2-D}{2}個である。これにD個を足して
\frac{D^2-D}{2}+D = \frac{D(D+1)}{2}
つまり独立なパラメータは\displaystyle \frac{D(D+1)}{2}個である。
演習 2.22
対称行列の逆行列も対称であることを示せ．
ある任意の対称行列\mathbf{A}があり、逆行列が存在する場合それを\mathbf{A}^{-1}とすると、
\mathbf{A}\mathbf{A}^{-1} = \mathbf{I}
となる（\mathbf{I}は単位行列）。両辺の転置を取り、\mathbf{A}は対称行列なので\mathbf{A} = \mathbf{A}^{\mathrm{T}}であることに注意すると
(\mathbf{A}\mathbf{A}^{-1})^{\mathrm{T}} = (\mathbf{A}^{-1})^{\mathrm{T}}\mathbf{A}^{\mathrm{T}} = (\mathbf{A}^{-1})^
{\mathrm{T}}\mathbf{A} = \mathbf{I}
ここで第3項について、逆行列の定義から
(\mathbf{A}^{-1})^{\mathrm{T}} = \mathbf{A}^{-1}
とならなければならないことがわかる。これは対称行列の逆行列\mathbf{A}^{-1}も対称行列となっていることを表している。
演習 2.23
\mathbf{\Sigma}=\sum_{i=1}^{D} \lambda_{i} \mathbf{u}_{i} \mathbf{u}_{i}^{\mathrm{T}} \tag{2.48}
の固有ベクトル展開を用いて座標系を対角化することで，マハラノビス距離\Deltaが定数になる超楕円体の内部の体積が，
V_{D}|\mathbf{\Sigma}|^{1 / 2} \Delta^{D} \tag{2.286}
になることを示せ．ただし，V_{D}はD次元単位球の体積で，マハラノビス距離は (2.44)で定義される．
※ PRML第2章 2.3 ガウス分布の議論にあるガウス分布の幾何的な形状についての理解を深めるための問題です。
※ そもそも超楕円体って何？って調べてみても意外とGoogleでヒットしないのですが、以下の定義を使います。
楕円体は，2次曲面の一種です．2次元において，次の方程式：
\frac{x^2}{a^2}+\frac{y^2}{b^2}=1
で表現される図形を楕円と呼びますが，これのn次元へ拡張したものと捉えて問題ありません．より厳密な呼び分けとしては，n=3のときのみ楕円体と呼び，n\ge4のとき超楕円体と呼ぶ場合もあるようです．
http://ssr-yuki.hatenablog.com/entry/2020/04/26/230647
マハラノビス距離の2乗\displaystyle \Delta^2 = (\mathbf{x}-\mathbf{\mu})^{\mathrm{T}}\mathbf{\Sigma}^{-1}(\mathbf{x}-\mathbf{\mu})は、P.78の手続きから固有ベクトル展開を用いて座標系を対角化することで(2.50)式\displaystyle \Delta^2 = \sum_{i=1}^{D}\frac{y_i^2}{\lambda_i}と書くことができる。例としてD=2であれば
\frac{y_1^2}{\lambda_1}+\frac{y_2^2}{\lambda_2}=\Delta^2
と書ける。これは平面図形の楕円である（P.79の図2.7のイメージ）。ちなみにD=3では楕円体（Wikipediaの楕円体を参照）を表す式になり、D \ge 4では超楕円体を表す。
D次元の超楕円体の体積V_eは以下の式で定義される。
V_e = \int\int\cdots\int dy_1dy_2\cdots dy_D
これは3次元の場合の式V_3 = \int\int\int dxdydzの拡張です。この辺についての説明は 楕円の面積と楕円体の体積の求め方のページも参考にしてみてください。
今、マハラノビス距離\Deltaは定数ということになっているので、超楕円体はa_i^2=y_i^2/\lambda_iの変数変換を行うことで、半径\Deltaの超球へと変換させることができる。つまりヤコビアン\mathbf{J}を使って表現すると
\begin{aligned}
V_e &= \int\int\cdots\int dy_1dy_2\cdots dy_D \\
&= \int\int\cdots\int |\mathbf{J}|da_1da_2\cdots da_D
\end{aligned}
となる。ここでヤコビアンは(2.53)-(2.55)での議論から
\mathbf{J}=\left(\begin{array}{ccc}\frac{\partial y_{1}}{\partial a_{1}} & \cdots & \frac{\partial y_{1}}{\partial a_{D}} \\ \vdots & \ddots & \vdots \\ \frac{\partial y_{D}}{\partial a_{1}} & \cdots & \frac{\partial y_{0}}{\partial a_{2}}\end{array}\right)=\left(\begin{array}{ccc}
\sqrt{\lambda_{1}} & & & 0 \\
& \sqrt{\lambda_{2}} \\
& & \ddots & \\
0 & & & \sqrt{\lambda_{D}}
\end{array}\right)
なので、|\mathbf{J}| = \prod_{i=1}^{D}\lambda_i^{1/2} =|\mathbf{\Sigma}|^{1/2}となる。
一方、\displaystyle \int\int\cdots\int da_1da_2\cdots da_D = \int \prod_{i=1}^{D}da_i部分は、半径\DeltaのD次元超球の体積を表しているので（超球は演習問題1.18でも登場。各変数a_iの定義域は-\Delta \le a_i \le \Deltaである）、問題文の通りにV_DをD次元単位球の体積とすると、
\int\int\cdots\int da_1da_2\cdots da_D = V_D\Delta^D
となる。
以上から求める超楕円体の内部の体積V_eは(2.286)式の通りに
V_e = V_D|\mathbf{\Sigma}|^{1/2}\Delta^D
となることが示された。
マハラノビス距離の直感的な理解としては、例えば https://mathwords.net/mahalanobis などのサイトの説明を読んでください。多次元からなるデータ群の中で例えば外れ値を検出したい場合、データの各次元への分散まで考慮したデータ群からの距離を考える必要があります。これを実現するのがマハラノビス距離です。マハラノビス距離が大きい → その点での確率密度が小さい → 異常度が高いと考えることができます。
演習 2.24
\left(\begin{array}{cc}\mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D}\end{array}\right)^{-1}=\left(\begin{array}{cc}\mathbf{M} & -\mathbf{M B D}^{-1} \\ -\mathbf{D}^{-1} \mathbf{C M} & \mathbf{D}^{-1}+\mathbf{D}^{-1} \mathbf{C M B D}^{-1}\end{array}\right) \tag{2.76}
の両辺に次の行列\displaystyle \left(\begin{array}{ll}\mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D}\end{array}\right) \ (2.287)を掛け，また，
\mathbf{M}=\left(\mathbf{A}-\mathbf{B} \mathbf{D}^{-1} \mathbf{C}\right)^{-1} \tag{2.77}
の定義を用いることで，(2.76)の恒等式を証明せよ．
指示通り(2.76)の右辺に左から\displaystyle \left(\begin{array}{cc}\mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \end{array}\right)をかけたものを\mathbf{X}とおく。これが左辺に左から\displaystyle \left(\begin{array}{cc}\mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \end{array}\right)をかけたもの、すなわち単位行列\mathbf{I}になっていることを示せば良い。
\begin{aligned} \mathbf{X}
&=\left(\begin{array}{cc}\mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \end{array}\right)
\left(\begin{array}{cc}\mathbf{M} & -\mathbf{MBD}^{-1} \\ -\mathbf{D}^{-1} \mathbf{CM} & \mathbf{D}^{-1}+\mathbf{D}^{-1} \mathbf{CMBD}^{-1}\end{array}\right) \\
&=\left(\begin{array}{cc}\mathbf{AM}-\mathbf{BD}^{-1}\mathbf{CM} & -\mathbf{AMBD}^{-1}+\mathbf{B}\left(\mathbf{D}^{-1}+\mathbf{D}^{-1} \mathbf{CMBD}^{-1}\right) \\ \mathbf{CM}-\mathbf{DD}^{-1} \mathbf{CM} & -\mathbf{CMBD}^{-1}+\mathbf{D}\left(\mathbf{D}^{-1}+\mathbf{D}^{-1} \mathbf{CMBD}^{-1}\right)\end{array}\right)
\end{aligned}
このそれぞれの部分行列成分について計算していくと
\begin{aligned}
\mathbf{X}_{11}
&=\mathbf{AM}-\mathbf{BD}^{-1} \mathbf{CM} \\
&=\left(\mathbf{A}-\mathbf{BD}^{-1} \mathbf{C}\right) \mathbf{M} \\ &=\left(\mathbf{A}-\mathbf{B D}^{-1} \mathbf{C}\right) \left(\mathbf{A}-\mathbf{BD}^{-1} \mathbf{C}\right)^{-1} \\ &=\mathbf{I}
\end{aligned}
\begin{aligned}
\mathbf{X}_{12}
&=-\mathbf{A M B D}^{-1}+\mathbf{B}\left(\mathbf{D}^{-1}+\mathbf{D}^{-1} \mathbf{C M B D}^{-1}\right) \\
&=-\mathbf{A M B D}^{-1}+\mathbf{B D}^{-1}+\mathbf{B D}^{-1}\mathbf{C M B D}^{-1} \\
&=-\left(\mathbf{A}-\mathbf{B D}^{-1} \mathbf{C}\right) \mathbf{M B D}^{-1}+\mathbf{B D}^{-1} \\
&=-\left(\mathbf{A}-\mathbf{B D}^{-1} \mathbf{C}\right)\left(\mathbf{A}-\mathbf{B D}^{-1} \mathbf{C}\right)^{-1} \mathbf{B D}^{-1}+\mathbf{B D}^{-1} \\
&=-\mathbf{B D}^{-1}+\mathbf{B D}^{-1} \\
&=\mathbf{O}
\end{aligned}
\begin{aligned}
\mathbf{X}_{21}
&=\mathbf{CM}-\mathbf{DD}^{-1}\mathbf{CM} \\
&=\mathbf{O}
\end{aligned}
\begin{aligned}
X_{22}
&=-\mathbf{C M B D}^{-1}+\mathbf{D}\left(\mathbf{D}^{-1}+\mathbf{D}^{-1} \mathbf{C M B D}^{-1}\right) \\
&=-\mathbf{C M B D}^{-1}+\mathbf{I}+\mathbf{C M B D}^{-1} \\ &=\mathbf{I}
\end{aligned}
よって全体として\mathbf{X} = \mathbf{I}となっていることが示せたので、(2.76)の恒等式は示された。
演習 2.25
2.3.1節や2.3.2節では，多変量ガウス分布の条件付き分布や周辺分布について考察した．より一般的に，\mathbf{x}の要素を\mathbf{x}_a, \mathbf{x}_b,および\mathbf{x}_cの3つに分けることを考える．この分割により，対応する平均ベクトル\boldsymbol{\mu}と共分散行列\mathbf{\Sigma}は
\boldsymbol{\mu}=\left(\begin{array}{c}\boldsymbol{\mu}_{a} \\ \boldsymbol{\mu}_{b} \\ \boldsymbol{\mu}_{c}\end{array}\right), \quad \mathbf{\Sigma}=\left(\begin{array}{ccc}\mathbf{\Sigma}_{a a} & \mathbf{\Sigma}_{a b} & \mathbf{\Sigma}_{a c} \\ \mathbf{\Sigma}_{b a} & \mathbf{\Sigma}_{b b} & \mathbf{\Sigma}_{b c} \\ \mathbf{\Sigma}_{c a} & \mathbf{\Sigma}_{c b} & \mathbf{\Sigma}_{c c}\end{array}\right)
のように分割される. 2.3節の結果を用いて\mathbf{x}_cを周辺化で消去した条件付き分布p(\mathbf{x}_a|\mathbf{x}_b)の式を求めよ．
\mathbf{x}_cを消去したときの同時分布p(\mathbf{x}_a,\mathbf{x}_b)は、平均ベクトルと共分散行列が
\boldsymbol{\mu}=\left(\begin{array}{c}\boldsymbol{\mu}_{a} \\ \boldsymbol{\mu}_{b} \end{array}\right), \quad \mathbf{\Sigma}=\left(\begin{array}{ccc}\mathbf{\Sigma}_{a a} & \mathbf{\Sigma}_{a b} \\ \mathbf{\Sigma}_{b a} & \mathbf{\Sigma}_{b b} \end{array}\right)
のガウス分布となる。よって条件付き分布p(\mathbf{x}_a|\mathbf{x}_b)もガウス分布であり、その平均は(2.81)(2.82)で与えられ、式は(2.96)となる。
演習 2.26
非常に有用な線形代数の結果であるWoodbury行列反転公式(Woodbury matrix inversion formula)は
(\mathbf{A}+\mathbf{B C D})^{-1}=\mathbf{A}^{-1}-\mathbf{A}^{-1} \mathbf{B}\left(\mathbf{C}^{-1}+\mathbf{D} \mathbf{A}^{-1} \mathbf{B}\right)^{-1} \mathbf{D} \mathbf{A}^{-1} \tag{2.289}
である．この両辺に(\mathbf{A}+\mathbf{B C D})を掛けて，この公式を証明せよ.
右辺に(\mathbf{A}+\mathbf{B C D})を掛けると
\begin{aligned}
& (\mathbf{A}+\mathbf{B C D}) (\mathbf{A}^{-1}-\mathbf{A}^{-1} \mathbf{B}\left(\mathbf{C}^{-1}+\mathbf{D} \mathbf{A}^{-1} \mathbf{B}\right)^{-1} \mathbf{D} \mathbf{A}^{-1}) \\
&= \mathbf{I} - \mathbf{B}\left(\mathbf{C}^{-1}+\mathbf{D} \mathbf{A}^{-1} \mathbf{B}\right)^{-1} \mathbf{D} \mathbf{A}^{-1} + \mathbf{B C D}\mathbf{A}^{-1} - \mathbf{B C D} \mathbf{A}^{-1} \mathbf{B}\left(\mathbf{C}^{-1}+\mathbf{D} \mathbf{A}^{-1} \mathbf{B}\right)^{-1} \mathbf{D} \mathbf{A}^{-1}\\
&= \mathbf{I} + \mathbf{B C D}\mathbf{A}^{-1} - \mathbf{B}\left(\mathbf{I}+\mathbf{C D} \mathbf{A}^{-1} \mathbf{B}\right) \left(\mathbf{C}^{-1}+\mathbf{D} \mathbf{A}^{-1} \mathbf{B}\right)^{-1} \mathbf{D} \mathbf{A}^{-1} \\
&= \mathbf{I} + \mathbf{B C D}\mathbf{A}^{-1} - \mathbf{B C}\left(\mathbf{C}^{-1}+\mathbf{D} \mathbf{A}^{-1} \mathbf{B}\right) \left(\mathbf{C}^{-1}+\mathbf{D} \mathbf{A}^{-1} \mathbf{B}\right)^{-1} \mathbf{D} \mathbf{A}^{-1} \\
&= \mathbf{I} + \mathbf{B C D}\mathbf{A}^{-1} - \mathbf{B C D}\mathbf{A}^{-1} \\
&= \mathbf{I}
\end{aligned}
となり示された。
演習 2.27
\mathbf{x}と\mathbf{z}を２つの独立な確率ベクトル，すなわち，p(\mathbf{x, z}) = p(\mathbf{x})p(\mathbf{z})であるとする．これらの和\mathbf{y}=\mathbf{x}+\mathbf{z}の平均が，それぞれの変数について個別に求めた平均の和となることを示せ．同様に，\mathbf{y}の共分散行列が，\mathbf{x}と\mathbf{z}それぞれの共分散行列の和であることを示せ．これが，演習問題1.10の結果と一致することを確認せよ．
\mathbb{E}[\mathbf{y}] = \mathbb{E}[\mathbf{x}+\mathbf{z}]=\int\int(\mathbf{x}+\mathbf{z})p(\mathbf{x},\mathbf{z})d\mathbf{x}d\mathbf{z}
\mathbf{x}と\mathbf{z}は独立であるから
\begin{aligned}
\int\int(\mathbf{x}+\mathbf{z})p(\mathbf{x},\mathbf{z})d\mathbf{x}d\mathbf{z} &= \int\int(\mathbf{x}+\mathbf{z})p(\mathbf{x})p(\mathbf{z})d\mathbf{x}d\mathbf{z}\\
&= \int\int \mathbf{x}p(\mathbf{x})p(\mathbf{z})d\mathbf{x}d\mathbf{z} + \int\int \mathbf{z}p(\mathbf{z})p(\mathbf{x})d\mathbf{x}d\mathbf{z} \\
&=\int \mathbf{x}p(\mathbf{x})d\mathbf{x} + \int \mathbf{z}p(\mathbf{z})d\mathbf{z} \\
&=\mathbb{E}[\mathbf{x}]+\mathbb{E}[\mathbf{z}]
\end{aligned}
以上より\mathbb{E}[\mathbf{y}]=\mathbb{E}[\mathbf{x}]+\mathbb{E}[\mathbf{z}]
共分散行列の定義より
\begin{aligned}
cov[y]
&= \mathbb{E}[(y-\mathbb{E}[y])(y-\mathbb{E}[y])^T]\\
&= \mathbb{E}[(x-\mathbb{E}[\mathbf{x}]+z-\mathbb{E}[z])(x-\mathbb{E}[\mathbf{x}]+z-\mathbb{E}[z])^T]\\
\end{aligned}
\mathbf{x}-\mathbb{E}[\mathbf{x}]=A、\mathbf{z}-\mathbb{E}[\mathbf{z}]=Bと置くと
\begin{aligned}
cov[\mathbf{y}]
&= \mathbb{E}[(A+B)(A+B)^T]\\
&= \mathbb{E}[AA^T]+\mathbb{E}[AB^T]+\mathbb{E}[BA^T]+\mathbb{E}[BB^T]
\end{aligned}
独立な変数同士の共分散は0になることから
\operatorname{cov}[\mathbf{y}] = \mathbb{E}[AA^T]+\mathbb{E}[BB^T] = \operatorname{cov}[\mathbf{x}]+\operatorname{cov}[\mathbf{z}]
演習 2.28
平均と共分散がそれぞれ
\mathbb{E}[\mathbf{z}]=
\left(\begin{array}{c}\boldsymbol{\mu} \\ \mathbf{A} \boldsymbol{\mu}+\mathbf{b}\end{array}\right) \tag{2.108}
と
\operatorname{cov}[\mathbf{z}]=\mathbf{R}^{-1}=\left(\begin{array}{cc}\mathbf{\Lambda}^{-1} & \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathbf{T}} \\ \mathbf{A} \mathbf{\Lambda}^{-1} & \mathbf{L}^{-1}+\mathbf{A} \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathbf{T}}\end{array}\right) \tag{2.105}
であるような，次の変数上の同時分布を考える．
\mathbf{z}=\left(\begin{array}{l}\mathbf{x} \\ \mathbf{y}\end{array}\right) \tag{2.190}
(2.92)\ \mathbb{E}[\mathbf{x}_a] = \boldsymbol{\mu}_aと(2.93)\ \operatorname{cov}[\mathbf{x}_a] = \mathbf{\Sigma}_{aa}の結果を用いて，周辺分布p(\mathbf{x})が\mathcal{N}(\mathbf{x}|\boldsymbol{\mu}, \mathbf{\Lambda}^{-1})\ (2.99)となることを示せ．
同様に，
\boldsymbol{\mu}_{a \mid b}=\boldsymbol{\mu}_{a}+\mathbf{\Sigma}_{ab} \mathbf{\Sigma}_{bb}^{-1}\left(\mathbf{x}_{b}-\boldsymbol{\mu}_{b}\right) \tag{2.81}
と
\mathbf{\Sigma}_{a \mid b}=\mathbf{\Sigma}_{a a}-\mathbf{\Sigma}_{a b} \mathbf{\Sigma}_{b b}^{-1} \mathbf{\Sigma}_{b a} \tag{2.82}
の結果を用いて，条件付き分布p(\mathbf{y}|\mathbf{x})が
p(\mathbf{y|x}) = \mathcal{N}(\mathbf{y}|\mathbf{Ax+b}, \mathbf{L}^{-1}) \tag{2.100}
となることを示せ．
(2.92)と(2.93)を(2.98)に代入すると，p(\mathbf{x})=\mathcal{N}(\mathbf{x}|\boldsymbol{\mu}, \mathbf{\Lambda}^{-1})であることが簡単に示された．
また
\begin{aligned}
\mathbb{E}[\mathbf{z}] &= \left(\begin{array}{c}\mu\\\mathbf{A}\mu+\mathbf{b}\\\end{array}\right)\\ cov[\mathbf{z}] &= \left(\begin{array}{cc}\Lambda^{-1} & \Lambda^{-1}\mathbf{A}^\top \\\mathbf{A}\Lambda^{-1} & \mathbf{L}^{-1} + \mathbf{A}\Lambda^{-1}\mathbf{A}^\top\\\end{array}\right)\end{aligned}
を(2.81)と(2.82)に代入すると
\begin{aligned}
\mu_{\mathbf{y}|\mathbf{x}} &= \mu_{\mathbf{y}} + \Lambda_{\mathbf{yx}}\Lambda_{\mathbf{xx}}^{-1}(\mathbf{x} - \mu_{\mathbf{x}})\\
&= \mathbf{A}\mu + \mathbf{b} + \mathbf{A}\Lambda^{-1}\Lambda(\mathbf{x}-\mu)\\
&= \mathbf{A}\mu + \mathbf{b} + \mathbf{Ax} - \mathbf{A\mu}\\
&= \mathbf{Ax+b}\\
\Sigma_{\mathbf{y}|\mathbf{x}} &= \Sigma_{\mathbf{yy}} - \Sigma_{\mathbf{yx}}\Sigma_{xx}^{-1}\Sigma_{\mathbf{xy}}\\
&= \mathbf{L}^{-1} + \mathbf{A}\Lambda^{-1}\mathbf{A}^\top - \mathbf{A}\Lambda^{-1}\Lambda\Lambda^{-1}\mathbf{A}^\top \\
&= \mathbf{L}^{-1}\\
\end{aligned}
となって，p(\mathbf{y|x}) = \mathcal{N}(\mathbf{y}|\mathbf{Ax+b}, \mathbf{L}^{-1})も示された．
演習 2.29
分割行列の逆行列の公式
\left(\begin{array}{cc}\mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D}\end{array}\right)^{-1}=\left(\begin{array}{cc}\mathbf{M} & -\mathbf{M B D}^{-1} \\ -\mathbf{D}^{-1} \mathbf{C M} & \mathbf{D}^{-1}+\mathbf{D}^{-1} \mathbf{C M B D}^{-1}\end{array}\right) \tag{2.76}
を用いて，精度行列
\mathbf{R}=\left(\begin{array}{cc}\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{L} \mathbf{A} & -\mathbf{A}^{\mathrm{T}} \mathbf{L} \\ -\mathbf{LA} & \mathbf{L}\end{array}\right) \tag{2.104}
の逆行列が，共分散行列
\operatorname{cov}
[\mathbf{z}]=\mathbf{R}^{-1}=\left(\begin{array}{cc}\mathbf{\Lambda}^{-1} & \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}} \\ \mathbf{A} \mathbf{\Lambda}^{-1} & \mathbf{L}^{-1}+\mathbf{A} \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}}\end{array}\right) \tag{2.105}
となることを示せ．
\mathbf{M}=\left(\mathbf{A}-\mathbf{B} \mathbf{D}^{-1} \mathbf{C}\right)^{-1}
(2.76)より、\mathbf{R}^{-1}について
\begin{aligned}( \text{左上} ) &= ( \mathbf{\Lambda} + \mathbf{A}^T \mathbf{L} \mathbf{A} - (- \mathbf{A}^T \mathbf{L}) \mathbf{L}^{-1} \mathbf{L} \mathbf{A})^{-1} \\&= ( \mathbf{\Lambda} + \mathbf{A}^T \mathbf{L} \mathbf{A} - \mathbf{A}^T \mathbf{L} \mathbf{A} )^{-1} \\&= \mathbf{\Lambda}^{-1} \end{aligned}
\begin{aligned}( \text{右上} ) &= - \mathbf{\Lambda}^{-1} (- \mathbf{A}^T \mathbf{L}) \mathbf{L}^{-1} \\&= \mathbf{\Lambda}^{-1} \mathbf{A}^T \end{aligned}
\begin{aligned}( \text{左下} ) &= - \mathbf{L}^{-1} (- \mathbf{L} \mathbf{A}) \mathbf{\Lambda}^{-1} \\&= \mathbf{A} \mathbf{\Lambda}^{-1} \end{aligned}
\begin{aligned}( \text{右下} ) &= \mathbf{L}^{-1} + \mathbf{L}^{-1} (- \mathbf{L} \mathbf{A}) \mathbf{\Lambda}^{-1} (- \mathbf{A}^T \mathbf{L}) \mathbf{L}^{-1} \\&= \mathbf{L}^{-1} + \mathbf{A} \mathbf{\Lambda}^{-1} \mathbf{A}^T \end{aligned}
これは共分散行列(2.105)と一致する。
演習 2.30
\mathbb{E}[\mathbf{z}]=\mathbf{R}^{-1}\left(\begin{array}{c}\mathbf{\Lambda} \boldsymbol{\mu}-\mathbf{A}^{\mathrm{T}} \mathbf{L} \mathbf{b} \\ \mathbf{Lb}\end{array}\right) \tag{2.107}
に，
\operatorname{cov}
[\mathbf{z}]=\mathbf{R}^{-1}=\left(\begin{array}{cc}\mathbf{\Lambda}^{-1} & \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}} \\ \mathbf{A} \mathbf{\Lambda}^{-1} & \mathbf{L}^{-1}+\mathbf{A} \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}}\end{array}\right) \tag{2.105}
の結果を用いて，
\mathbb{E}[\mathbf{z}]=\left(\begin{array}{c}\boldsymbol{\mu} \\
\mathbf{A} \boldsymbol{\mu}+\mathbf{b}\end{array}\right) \tag{2.108}
を確かめよ．
単純に計算するだけ
\begin{aligned}
\mathbb{E}[\mathbf{z}] &= \mathbf{R}^{-1}\left(\begin{array}{c}\mathbf{\Lambda} \boldsymbol{\mu}-\mathbf{A}^{\mathrm{T}} \mathbf{L} \mathbf{b} \\ \mathbf{Lb}\end{array}\right) \\
&= \left(\begin{array}{cc}\mathbf{\Lambda}^{-1} & \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}} \\ \mathbf{A} \mathbf{\Lambda}^{-1} & \mathbf{L}^{-1}+\mathbf{A} \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}}\end{array}\right)\left(\begin{array}{c}\mathbf{\Lambda} \boldsymbol{\mu}-\mathbf{A}^{\mathrm{T}} \mathbf{L b} \\ \mathbf{L b}\end{array}\right) \\
&=\left(\begin{array}{c}\mathbf{\Lambda}^{-1}\left(\mathbf{\Lambda} \boldsymbol{\mu}-\mathbf{A}^{\mathrm{T}} \mathbf{L b}\right)+\mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{L b} \\
\mathbf{A} \mathbf{\Lambda}^{-1}\left(\mathbf{\Lambda} \boldsymbol{\mu}-\mathbf{A}^{\mathrm{T}} \mathbf{L} \mathbf{b}\right)+\left(\mathbf{L}^{-1}+\mathbf{A} \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}}\right) \mathbf{L b}\end{array}\right) \\
&=\left(\begin{array}{c}\boldsymbol{\mu}-\mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{L} \mathbf{b}+\mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{L} \mathbf{b}
\\ \mathbf{A} \boldsymbol{\mu}-\mathbf{A} \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{L} \mathbf{b}+\mathbf{b}+\mathbf{A} \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{L} \mathbf{b}\end{array}\right) \\
&=\left(\begin{array}{c}\boldsymbol{\mu} \\
\mathbf{A} \boldsymbol{\mu}+\mathbf{b}\end{array}\right)
\end{aligned}
演習 2.31
2つの多次元確率ベクトル\mathbf{x}と\mathbf{z}を考える．これらは，それぞれガウス分布p(\mathbf{x}) = \mathcal{N}(\mathbf{x}| \boldsymbol{\mu}_{\mathbf{x}}, \mathbf{\Sigma_x})とp(\mathbf{z}) = \mathcal{N}(\mathbf{z}| \boldsymbol{\mu}_{\mathbf{z}}, \mathbf{\Sigma_z})に従い，これらの和は\mathbf{y} = \mathbf{x} + \mathbf{z}であるとする．周辺分布p(\mathbf{x})と条件付き分布p(\mathbf{y|x})の積からなる線形ガウスモデルを用いて，
\mathbb{E}[\mathbf{y}] =\mathbf{A} \boldsymbol{\mu}+\mathbf{b} \tag{2.109}
\operatorname{cov}[\mathbf{y}] = \mathbf{L}^{-1}+\mathbf{A}\mathbf{\Lambda}^{-1}\mathbf{A}^{\mathrm{T}} \tag{2.110}
から周辺分布p(\mathbf{y})についての式を求めよ．
条件付き確率分布p(\mathbf{y|x})は\mathbf{x}が与えられた（つまり定数とした）中での\mathbf{y} = \mathbf{x}+\mathbf{z}の確率分布を表しているので、\mathbf{y}の平均は\boldsymbol{\mu}_{\mathbf{z}}に\mathbf{x}を足したもの、\mathbf{y}の分散は\mathbf{\Sigma}_{\mathbf{z}}となるので、
p(\mathbf{y|x}) = \mathcal{N}(\mathbf{y}| \boldsymbol{\mu}_{\mathbf{z}}+\mathbf{x}, \mathbf{\Sigma_z})
と表すことができる。
以降、pp.88〜90の線形ガウスモデルの議論とまとめ(2.113)〜(2.117)と、この問題設定の
\begin{aligned}
p(\mathbf{x}) &= \mathcal{N}(\mathbf{x}|\boldsymbol{\mu}_{\mathbf{x}},\mathbf{\Sigma}_{\mathbf{x}}) \\
p(\mathbf{y|x}) &= \mathcal{N}(\mathbf{y}|\mathbf{x}+\boldsymbol{\mu}_{\mathbf{z}},\mathbf{\Sigma}_{\mathbf{z}})
\end{aligned}
を比較すると、
\boldsymbol{\mu} \to \boldsymbol{\mu}_{\mathbf{x}}, \mathbf{\Lambda}^{-1} \to \mathbf{\Sigma}_{\mathbf{x}}, \mathbf{A} \to \mathbf{I}, \mathbf{b} \to \boldsymbol{\mu}_{\mathbf{z}}, \mathbf{L}^{-1} \to \mathbf{\Sigma}_{\mathbf{z}}
と置換することで、(2.115)式から
p(\mathbf{y}) = \mathcal{N}(\mathbf{y}|\boldsymbol{\mu}_{\mathbf{x}}+\boldsymbol{\mu}_{\mathbf{z}}, \mathbf{\Sigma}_{\mathbf{x}}+\mathbf{\Sigma}_{\mathbf{z}})
となる。これは演習2.27の結果と同じである。
演習 2.32
これと次の演習問題で，線形ガウスモデル中の二次形式の操作を練習し，また，本文中で導いた結果も検証する．周辺分布
p(\mathbf{x})=\mathcal{N}\left(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}^{-1}\right) \tag{2.99}
と条件付き分布
p(\mathbf{y} \mid \mathbf{x})=\mathcal{N}\left(\mathbf{y} \mid \mathbf{Ax}+\mathbf{b}, \mathbf{L}^{-1}\right) \tag{2.100}
で定義される同時確率p(\mathbf{x},\mathbf{y})を考える．同時確率の指数部分の二次形式に，2.3節で述べた平方完成の技法を適用して，変数\mathbf{x}を積分消去した周辺分布p(\mathbf{y})の平均と共分散の式を求めよ．これには，Woodbury行列反転公式
(\mathbf{A}+\mathbf{BCD})^{-1}=\mathbf{A}^{-1}-\mathbf{A}^{-1} \mathbf{B}\left(\mathbf{C}^{-1}+\mathbf{DA}^{-1} \mathbf{B}\right)^{-1} \mathbf{DA}^{-1} \tag{2.289}
を用いる．これらが，2章の結果を用いて得た結果
\mathbb{E}[\mathbf{y}] =\mathbf{A} \boldsymbol{\mu}+\mathbf{b} \tag{2.109}
\operatorname{cov}[\mathbf{y}] = \mathbf{L}^{-1}+\mathbf{A}\mathbf{\Lambda}^{-1}\mathbf{A}^{\mathrm{T}} \tag{2.110}
と一致することを確かめよ．
※ 2.3.3 ガウス変数に対するベイズの定理における導出の流れに則りつつも、問題文の指示によれば途中から二次形式を使ったやり方で導出するよう求めているため、2.3.2 周辺ガウス分布で行ったような二次形式を使ったやり方で求めていきます。結果は当然変わらないですが、こちらの方が2.3.3の行列形式を使ったやり方よりも計算量がものすごく増えます。
演習2.31の内容と同様に、p(\mathbf{x,y}) = p(\mathbf{y|x})p(\mathbf{x})から
\begin{aligned}
p(\mathbf{y}) &= \int p(\mathbf{x}, \mathbf{y}) d\mathbf{x} \\
&= \int \mathcal{N}\left(\mathbf{y} \mid \mathbf{Ax}+\mathbf{b}, \mathbf{L}^{-1}\right) \mathcal{N}\left(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}^{-1}\right) d\mathbf{x}
\end{aligned}
となる。正規化係数を無視して指数部分だけを考えると
\exp \left\{-\frac{1}{2}\left(\mathbf{y}-\mathbf{Ax}-\mathbf{b}\right)^{\mathrm{T}}\mathbf{L}(\mathbf{y}-\mathbf{Ax}-\mathbf{b})-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\mathrm{T}}\mathbf{\Lambda}(\mathbf{x}-\boldsymbol{\mu})\right\}
となるので、二次形式を展開すると
\begin{aligned}
&-\frac{1}{2}\left\{\mathbf{y}^{\mathrm{T}}\mathbf{L}\mathbf{y}-2\mathbf{y}^{\mathrm{T}}\mathbf{L}(\mathbf{Ax}+\mathbf{b})+(\mathbf{Ax}+\mathbf{b})^{\mathrm{T}}\mathbf{L}(\mathbf{Ax}+\mathbf{b})+\mathbf{x}^{\mathrm{T}}\mathbf{\Lambda}\mathbf{x}-2\boldsymbol{\mu}^{\mathrm{T}} \mathbf{\Lambda} \mathbf{x}+\boldsymbol{\mu}^{\mathrm{T}} \mathbf{\Lambda} \boldsymbol{\mu} \right\} \\
=&-\frac{1}{2}\left\{\mathbf{y}^{\mathrm{T}} \mathbf{L} \mathbf{y}-2 \mathbf{y}^{\mathrm{T}} \mathbf{L} \mathbf{Ax}-2 \mathbf{y}^{\mathrm{T}} \mathbf{L} \mathbf{b}+\mathbf{x}^{\mathrm{T}} \mathbf{A}^{\mathrm{T}} \mathbf{LAx}+2 \mathbf{b}^{\mathrm{T}} \mathbf{LAx}+\mathbf{b}^{\mathrm{T}} \mathbf{Lb}\right. \\
&+\left.\mathbf{x}^{\mathrm{T}} \mathbf{\Lambda} \mathbf{x}-2 \boldsymbol{\mu}^{\mathrm{T}} \mathbf{\Lambda} \mathbf{x}+\boldsymbol{\mu}^{\mathrm{T}} \mathbf{A} \boldsymbol{\mu} \right\}
\end{aligned}
ここで2.3.2の技法に則り、\mathbf{x}を積分消去することを目標として\mathbf{x}の関わる項をまとめてから積分を容易にするために平方完成の技法を使う。
\mathbf{x}の関わる項を上式から抽出すると
\begin{aligned}
&-\frac{1}{2}\left\{\mathbf{x}^{\mathrm{T}}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right) \mathbf{x}-2\left(\left(\mathbf{y}^{\mathrm{T}}-\mathbf{b}^{\mathrm{T}}\right) \mathbf{LA}+\boldsymbol{\mu}^{\mathrm{T}} \mathbf{\Lambda}\right) \mathbf{x}\right\} \\
=&-\frac{1}{2}(\mathbf{x}-\mathbf{m})^{\mathrm{T}}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)(\mathbf{x}-\mathbf{m})+\frac{1}{2} \mathbf{m}^{\mathrm{T}}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right) \mathbf{m}
\end{aligned}
ここで\mathbf{m}=\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1}(\mathbf{\Lambda}^{\mathrm{T}} \boldsymbol{\mu}+\mathbf{A}^{\mathrm{T}}\mathbf{L}(\mathbf{y}-\mathbf{b}))とした。
指数内に戻して考えると
\int\exp\left\{
-\frac{1}{2}(\mathbf{x}-\mathbf{m})^{\mathrm{T}}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)(\mathbf{x}-\mathbf{m})
\right\}d\mathbf{x}
この積分は正規化されていないガウス分布の標準的な二次形式になっているので、正規化係数の逆数になる。その逆数は|(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}}\mathbf{L A})^{-1}|、つまり共分散行列のみに依存することがガウス分布の性質から分かるので\mathbf{x}には依存しなくなっている。つまりこれは\mathbf{x}が積分によって消去されていることを意味する。
まだ\frac{1}{2} \mathbf{m}^{\mathrm{T}}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right) \mathbf{m}部分が残っているので、これを\mathbf{y}に依存する項とあわせて再び計算すると（\mathbf{y}に依存しない項は\text{const.}とした）
\begin{aligned}
&=-\frac{1}{2}\left\{\mathbf{y}^{\mathrm{T}} \mathbf{Ly}-2 \mathbf{y}^{\mathrm{T}} \mathbf{Lb}+\mathbf{m}^{\mathrm{T}}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right) \mathbf{m}\right\}+\text { const. } \\
&=-\frac{1}{2}\left\{\mathbf{y}^{\mathrm{T}} \mathbf{Ly}-2 \mathbf{y}^{\mathrm{T}} \mathbf{Lb}+\mathbf{y}^{\mathrm{T}} \mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{Ly}-2 \mathbf{y}^{\mathrm{T}}\left[\mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1}\left(\mathbf{\Lambda} \mu-\mathbf{A}^{\mathrm{T}} \mathbf{Lb}\right)\right] + \text{const.} \right\} \\
&=-\frac{1}{2} \mathbf{y}^{\mathrm{T}}\left[\mathbf{L}-\mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{L}\right] \mathbf{y}+\mathbf{y}^{\mathrm{T}}\left[\left(\mathbf{L}-\mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{L}\right)\mathbf{b}+ \mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{T}\mathbf{LA}\right)^{-1}\mathbf{\Lambda} \boldsymbol{\mu}\right]
\end{aligned}
これが周辺分布p(\mathbf{y})の指数部分となっているので、二次形式部分の\mathbf{L}-\mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{L}がp(\mathbf{y})の精度行列に相当することがわかる。
Woodburyの反転公式を使うと、共分散行列\operatorname{cov}[\mathbf{y}]は
\operatorname{cov}[\mathbf{y}] = (\mathbf{L}-\mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{L})^{-1} = \mathbf{L}^{-1}+\mathbf{A}\mathbf{\Lambda}^{-1}\mathbf{A}^{\mathrm{T}}
となる。
最後に平均\mathbb{E}[\mathbf{y}]は、(2.71), (2.75), (2.89)で議論したように、\mathbf{y}^{\mathrm{T}}の係数\left[\left(\mathbf{L}-\mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{L}\right)\mathbf{b}+\mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{T} \mathbf{LA}\right)^{-1}\mathbf{\Lambda} \boldsymbol{\mu}\right]が\{\operatorname{cov}[\mathbf{y}]\}^{-1}\mathbb{E}[\mathbf{y}]に一致することから求められる。
\begin{aligned}
\mathbb{E}[\mathbf{y}] &= \operatorname{cov}[\mathbf{y}]\left[\left(\mathbf{L}-\mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{L}\right)\mathbf{b}+\mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{T} \mathbf{LA}\right)^{-1}\mathbf{\Lambda} \boldsymbol{\mu}\right] \\
&=\left(\mathbf{L}^{-1}+\mathbf{A}\mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}}\right)\left(\left(\mathbf{L}-\mathbf{LA}(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA})^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{L}\right) \mathbf{b}+\mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \mathbf{\Lambda} \boldsymbol{\mu}\right) \\
&=\left(\mathbf{L}^{-1}+\mathbf{A}\mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}}\right)\left(\left(\mathbf{L}^{-1}+\mathbf{A}\mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}}\right)^{-1} \mathbf{b}+\mathbf{LA}\left(\Lambda+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \mathbf{\Lambda} \boldsymbol{\mu}\right) \\
&=\mathbf{b}+\left(\mathbf{L}^{-1}+\mathbf{A} \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}}\right) \mathbf{LA}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \mathbf{\Lambda} \boldsymbol{\mu} \\
&=\mathbf{b}+\left(\mathbf{A}+\mathbf{A} \mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)\left(\mathbf{\Lambda}\left(\mathbf{I}+\mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)\right)^{-1} \mathbf{\Lambda} \boldsymbol{\mu} \\
&=\mathbf{b}+\left\{\mathbf{A}\left(\mathbf{I}+\mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)\left(\mathbf{I}+\mathbf{\Lambda}^{-1} \mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \mathbf{\Lambda}^{-1} \mathbf{\Lambda} \boldsymbol{\mu}\right\} \\
&=\mathbf{A}\boldsymbol{\mu}+\mathbf{b}
\end{aligned}
よって(2.109)の平均が求められた。
演習 2.33
演習問題2.32と同じ同時分布について考えるが，今度は，平方完成の技法によって，条件付き分布p(\mathbf{x|y})の平均と共分散の式を求める．この結果も，
\mathbb{E}[\mathbf{x} \mid \mathbf{y}] =\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1}\left\{\mathbf{A}^{\mathrm{T}} \mathbf{L}(\mathbf{y}-\mathbf{b})+\mathbf{\Lambda} \boldsymbol{\mu}\right\} \tag{2.111}
\operatorname{cov}[\mathbf{x} \mid \mathbf{y}] =\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \tag{2.112}
の式と一致することを確かめよ．
※条件付き分布p(\mathbf{x\mid y})の場合、\mathbf{y}は固定である。演習2.32と同様に考えると
\begin{aligned}
p(\mathbf{x}\mid \mathbf{y}) &= \frac{p(\mathbf{x},\mathbf{y})}{p(\mathbf{y})}\\
&= \frac{p(\mathbf{x}, \mathbf{y})}{\int p(\mathbf{x}, \mathbf{y}) d\mathbf{x}}
\end{aligned}
であり、分母は\mathbf{y}が固定で\mathbf{x}について積分するので定数となる。よって、\mathbf{y}が固定された状態（定数とみなせる状態）でp(\mathbf{x},\mathbf{y})について平方完成し、指数部分を調べれば共分散\operatorname{cov}[\mathbf{x} \mid \mathbf{y}]が求まり、そこから平均\mathbb{E}[\mathbf{x} \mid \mathbf{y}]も求まる。
ここについては演習2.32と同じ流れで
\begin{aligned}
p(\mathbf{x}\mid \mathbf{y}) &= \frac{p(\mathbf{x}, \mathbf{y})}{\int p(\mathbf{x}, \mathbf{y}) d\mathbf{x}} \\
&= C_1 \mathcal{N}\left(\mathbf{y} \mid \mathbf{Ax}+\mathbf{b}, \mathbf{L}^{-1}\right) \mathcal{N}\left(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Lambda}^{-1}\right) \\
&= C_2 \exp\left\{-\frac{1}{2}\left\{\mathbf{x}^{\mathrm{T}}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right) \mathbf{x}-2\left(\left(\mathbf{y}^{\mathrm{T}}-\mathbf{b}^{\mathrm{T}}\right) \mathbf{LA}+\boldsymbol{\mu}^{\mathrm{T}} \mathbf{\Lambda}\right) \mathbf{x} + C_3 \right\} \right\}\\
&= C_4 \exp \left\{ -\frac{1}{2}(\mathbf{x}-\mathbf{m})^{\mathrm{T}}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)(\mathbf{x}-\mathbf{m})+\frac{1}{2} \mathbf{m}^{\mathrm{T}}\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right) \mathbf{m} \right\}
\end{aligned}
ここでC_1, C_2, C_3, C_4はそれぞれ定数で、\mathbf{m}=\left(\mathbf{\Lambda}+\mathbf{A}^{\top} \mathbf{LA}\right)^{-1}(\mathbf{\Lambda} \boldsymbol{\mu}+\mathbf{A}^{\mathrm{T}}\mathbf{L}(\mathbf{y}-\mathbf{b}))であるから、結局p(\mathbf{x}\mid \mathbf{y})の共分散は
\operatorname{cov}[\mathbf{x} \mid \mathbf{y}] =\left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1} \tag{2.112}
となり、
\mathbb{E}[\mathbf{x} \mid \mathbf{y}] = \mathbf{m} = \left(\mathbf{\Lambda}+\mathbf{A}^{\mathrm{T}} \mathbf{LA}\right)^{-1}\left\{\mathbf{A}^{\mathrm{T}} \mathbf{L}(\mathbf{y}-\mathbf{b})+\mathbf{\Lambda} \boldsymbol{\mu}\right\}
を得る。
演習 2.34
多変量ガウス分布の共分散行列の最尤推定解を求めるには，共分散行列が対称で正定値である制約の下で\mathbf{\Sigma}について対数尤度関数
\ln p(\mathbf{x} \mid \boldsymbol{\mu}, \mathbf{\Sigma})=-\frac{N D}{2} \ln (2 \pi)-\frac{N}{2} \ln |\mathbf{\Sigma}|-\frac{1}{2} \sum_{n=1}^{N}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)^{\mathrm{T}} \mathbf{\Sigma}^{-1}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right) \tag{2.118}
を最大化しなくてはならない．ここでは，こうした制約を無視して，ただ最大化することにする．付録Cの(C.21)，(C.26)，および(C.28)の結果を用いて，対数尤度関数(2.118)を最大化する共分散行列\mathbf{\Sigma}が，サンプル共分散
\mathbf{\Sigma}_{\mathbf{ML}}=\frac{1}{N} \sum_{n=1}^{N}\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{\mathbf{ML}}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{\mathbf{ML}}\right)^{\mathrm{T}} \tag{2.122}
となることを示せ．なお，（サンプル共分散が非特異なら）最終結果は対称かつ，正定値である必要がある．
※ ※ 公式解答集では共分散行列\mathbf{\Sigma}が対称行列であること（\mathbf{\Sigma}^{-1}=({\mathbf{\Sigma}^{-1}})^{\mathrm{T}}）の制約を無視して〜と書いてあるのに、その対称性を利用した解答例になっているが、それを利用しなくても答えを出すことはできる。
\begin{aligned}
\frac{\partial}{\partial \mathbf{\Sigma}} \ln p(\mathbf{X} | \boldsymbol{\mu}, \mathbf{\Sigma}) &=-\frac{N}{2} \frac{\partial}{\partial \mathbf{\Sigma}} \ln |\mathbf{\Sigma}|-\frac{1}{2} \frac{\partial}{\partial \mathbf{\Sigma}}\left\{\sum_{n=1}^{N}\left(\mathbf{x}_n-\boldsymbol{\mu}\right)^{\mathrm{T}} \Sigma^{-1} \left(\mathbf{x}_n-\boldsymbol{\mu} \right)\right\} \\
&=-\frac{N}{2}\left(\mathbf{\Sigma}^{-1}\right)^{\mathrm{T}}-\frac{1}{2} \frac{\partial}{\partial \mathbf{\Sigma}}\left\{\sum_{n=1}^{N} \operatorname{Tr}\left(\mathbf{\Sigma}^{-1}\left(\mathbf{x}_n-\boldsymbol{\mu}\right)\left(\mathbf{x}_n-\boldsymbol{\mu} \right)^{\mathrm{T}}\right)\right\} \\
&=-\frac{N}{2}\left(\mathbf{\Sigma}^{-1}\right)^{\mathrm{T}} - \frac{1}{2} \frac{\partial}{\partial \mathbf{\Sigma}} \operatorname{Tr}\left( \mathbf{\Sigma}^{-1} \sum_{n=1}^{N}\left(\mathbf{x}_n-\boldsymbol{\mu}\right)\left(\mathbf{x}_n-\boldsymbol{\mu} \right)^{\mathrm{T}}\right)
\end{aligned}
ここで第2項について、\mathbf{A}=\mathbf{\Sigma}, \mathbf{B}=\sum_{n=1}^{N}\left(\mathbf{x}_n-\boldsymbol{\mu}\right)\left(\mathbf{x}_n-\boldsymbol{\mu} \right)^{\mathrm{T}}とすると、求めるべきは\displaystyle \frac{\partial}{\partial \mathbf{A}}{\mathrm{Tr}}(\mathbf{A}^{-1}\mathbf{B})である。
\begin{aligned}
\frac{\partial}{\partial A_{ij}}{\mathrm{Tr}}(\mathbf{A}^{-1}\mathbf{B}) =& \operatorname{Tr}\left(\left(\frac{\partial}{\partial A_{i j}} \mathbf{A}^{-1} \right) \mathbf{B}\right) \\
=&- \operatorname{Tr}\left(\mathbf{A}^{-1}\left(\frac{\partial}{\partial A_{ij}} \mathbf{A}\right) \mathbf{A}^{-1} \mathbf{B}\right)\hspace{1em}(\because 付録(C.21))\\
=&-\operatorname{Tr}\left(\left(\frac{\partial}{\partial A_{ij}} \mathbf{A}\right) \mathbf{A}^{-1} \mathbf{B} \mathbf{A}^{-1}\right)\hspace{1em}(\because \operatorname{Tr}(\mathbf{ABCD}) = \operatorname{Tr}(\mathbf{BCDA}))
\end{aligned}
なので、\mathbf{C}=\mathbf{A^{-1}BA^{-1}}とすると
\begin{aligned} \operatorname{Tr}\left(\left(\frac{\partial}{\partial A_{ij}} \mathbf{A}\right) \mathbf{C}\right) &=\sum_{s}\left(\left(\frac{\partial}{\partial A_{ij}} \mathbf{A}\right) \mathbf{C}\right)_{ss} \\
&= \sum_{s}\left(\sum_{t}\left(\frac{\partial}{\partial A_{ij}} \mathbf{A}\right)_{s t} c_{ts}\right) \\
&=\sum_{s, t} \delta_{is} \delta_{jt} c_{ts}=c_{ji}
\end{aligned}
最後はs=iかつt=jのときのみクロネッカーのデルタ\delta_{is}\delta_{jt}が1になり、その他は0となることを表している。
ここについては、付録(C.23)に書かれてある定理に則って
\begin{aligned} \operatorname{Tr}\left(\left(\frac{\partial}{\partial A_{ij}} \mathbf{A}\right) \mathbf{C}\right)
&=\frac{\partial}{\partial A_{ij}} \operatorname{Tr}\left(\mathbf{AC} \right) \\
&=c_{ji}
\end{aligned}
としても良い。
（※ 一般に正方行列\mathbf{F}の行列のij成分f_{ij}について\operatorname{Tr}(\mathbf{F})=\sum_{i=1}f_{ii}より、
\frac{\partial}{\partial f_{ij}}\operatorname{Tr}(\mathbf{F}) = \frac{\partial f_{11}}{\partial f_{ij}}+\frac{\partial f_{22}}{\partial f_{ij}}+\cdots= \operatorname{Tr}\left(\frac{\partial}{\partial f_{ij}} \mathbf{F}\right)
が成立する）
これより、\displaystyle \frac{\partial}{\partial A_{ij}}{\mathrm{Tr}}(\mathbf{A}^{-1}\mathbf{B})=-c_{ji}となるので、付録(C.24)のように
\begin{aligned}
\frac{\partial}{\partial \mathbf{A}}{\mathrm{Tr}}(\mathbf{A}^{-1}\mathbf{B}) &= -\mathbf{C}^{\mathrm T} \\
&=-(\mathbf{A^{-1}BA^{-1}})^{\mathrm T} \\
&=-(\mathbf{\Sigma^{-1}B\Sigma^{-1}})^{\mathrm T}
\end{aligned}
よって、最大値を求めるために\frac{\partial}{\partial \mathbf{\Sigma}} \ln p(\mathbf{X} | \boldsymbol{\mu}, \mathbf{\Sigma}) = 0とすると
\begin{aligned}
-\frac{N}{2}\left(\mathbf{\Sigma}^{-1}\right)^{\mathrm{T}} - \frac{1}{2} \frac{\partial}{\partial \mathbf{\Sigma}} \operatorname{Tr}\left( \mathbf{\Sigma}^{-1} \sum_{n=1}^{N}\left(\mathbf{x}_n-\boldsymbol{\mu}\right)\left(\mathbf{x}_n-\boldsymbol{\mu} \right)^{\mathrm{T}}\right)
=0 \\
\frac{N}{2}\left(\mathbf{\Sigma}^{-1}\right)^{\mathrm{T}} = \frac{1}{2} \left( \mathbf{\Sigma^{-1}} \sum_{n=1}^{N}\left(\mathbf{x}_n-\boldsymbol{\mu}\right)\left(\mathbf{x}_n-\boldsymbol{\mu} \right)^{\mathrm{T}} \mathbf{\Sigma^{-1}} \right) ^{\mathrm T}
\end{aligned}
互いの転置をとって移項すると、
\mathbf{\Sigma} = \frac{1}{N}\sum_{n=1}^{N}\left(\mathbf{x}_n-\boldsymbol{\mu}\right)\left(\mathbf{x}_n-\boldsymbol{\mu} \right)^{\mathrm{T}}
となり、これは式(2.122)となる。また、これによって\mathbf{\Sigma}についての対称性を仮定せずに\mathbf{\Sigma}が対称行列になることがわかった。
演習 2.35
\mathbb{E}[\mathbf{x}] = \boldsymbol{\mu} \tag{2.59}
の結果を用いて，
\mathbb{E}[\mathbf{xx}^{\mathrm{T}}] = \boldsymbol{\mu\mu}^{\mathrm{T}}+\mathbf{\Sigma} \tag{2.62}
を証明せよ．そして，(2.59)と(2.62)の結果から，
\mathbb{E}\left[\mathbf{x}_{n} \mathbf{x}_{m}^{\mathrm{T}}\right]=\boldsymbol{\mu} \boldsymbol{\mu}^{\mathrm{T}}+I_{n m} \mathbf{\Sigma} \tag{2.291}
を示せ．ただし，\mathbf{x}_{n}は，平均が\boldsymbol{\mu}で共分散が\mathbf{\Sigma}のガウス分布からサンプリングされたデータ点，I_{nm}は単位行列の(n,m)要素を表す．これから，
\mathbb{E}[\mathbf{\Sigma}_{\mathrm{ML}}] = \frac{N-1}{N}\mathbf{\Sigma} \tag{2.124}
の結果を証明せよ.
(2.59)と(2.62)から、m=nであれば\mathbb{E}\left[\mathbf{x}_{n} \mathbf{x}_{m}\right]=\boldsymbol{\mu} \boldsymbol{\mu}^{\mathrm{T}}+\mathbf{\Sigma}となる。
一方で、m \ne nだったときには2つのデータセット\mathbf{x}_nと\mathbf{x}_mは独立なので\mathbb{E}\left[\mathbf{x}_{n} \mathbf{x}_{m}\right]=\boldsymbol{\mu} \boldsymbol{\mu}^{\mathrm{T}}となる。
よって、I_{nm}を使うと\mathbb{E}[\mathbf{x}_n\mathbf{x}_m^{\mathrm T}] = \boldsymbol{\mu}\boldsymbol{\mu}^{\mathrm{T}} + I_{nm}\mathbf{\Sigma}となる。
(2.124)式 \displaystyle \mathbb{E}\left[\mathbf{\Sigma}_{\mathrm{ML}}\right]=\frac{N-1}{N} \mathbf{\Sigma}を示す。
(2.122)式より
\begin{align}
\mathbb{E}\left[\mathbf{\Sigma}_{\mathrm{ML}}\right] &=\frac{1}{N} \sum_{n=1}^{N} \mathbb{E}\left[\left(\mathbf{x}_{n}-\frac{1}{N} \sum_{m=1}^{N} \mathbf{x}_{m}\right)\left(\mathbf{x}_{n}^{\mathrm{T}}-\frac{1}{N} \sum_{l=1}^{N} \mathbf{x}_{l}^{\mathrm{T}}\right)\right] \\
&=\frac{1}{N} \sum_{n=1}^{N} \mathbb{E}\left[\mathbf{x}_{n} \mathbf{x}_{n}^{\mathrm{T}}-\frac{2}{N} \mathbf{x}_{n} \sum_{m=1}^{N} \mathbf{x}_{m}^{\mathrm{T}}+\frac{1}{N^{2}} \sum_{m=1}^{N} \sum_{l=1}^{N} \mathbf{x}_{m} \mathbf{x}_{l}^{\mathrm{T}}\right] \\
&=\frac{1}{N}\sum_{n=1}^{N}\left\{ \mathbb{E}[\mathbf{x}_n\mathbf{x}_n^{\mathrm{T}}]-\frac{2}{N}\mathbb{E}\left[ \mathbf{x}_{n} \sum_{m=1}^{N} \mathbf{x}_{m}^{\mathrm{T}} \right] +\frac{1}{N^2}\mathbb{E}\left[ \sum_{m=1}^{N} \sum_{l=1}^{N} \mathbf{x}_{m} \mathbf{x}_{l}^{\mathrm{T}} \right] \right\} \\
&=\frac{1}{N}\sum_{n=1}^{N}\left\{\boldsymbol{\mu} \boldsymbol{\mu}^{\mathrm{T}}+\mathbf{\Sigma}-2\left(\boldsymbol{\mu} \boldsymbol{\mu}^{\mathrm{T}}+\frac{1}{N} \mathbf{\Sigma}\right)+\boldsymbol{\mu} \boldsymbol{\mu}^{\mathrm{T}}+\frac{1}{N} \mathbf{\Sigma}\right\} \\
&=\frac{N-1}{N} \mathbf{\Sigma}
\end{align}
演習 2.36
\begin{align}
\boldsymbol{\mu}_{\mathrm{ML}}^{(N)} &=\frac{1}{N} \sum_{n=1}^{N} \mathbf{x}_{n} \\
&=\frac{1}{N} \mathbf{x}_{N}+\frac{1}{N} \sum_{n=1}^{N-1} \mathbf{x}_{n} \\
&=\frac{1}{N} \mathbf{x}_{N}+\frac{N-1}{N} \boldsymbol{\mu}_{\mathrm{ML}}^{(N-1)} \\
&=\boldsymbol{\mu}_{\mathrm{ML}}^{(N-1)}+\frac{1}{N}\left(\mathbf{x}_{N}-\boldsymbol{\mu}_{\mathrm{ML}}^{(N-1)}\right) \tag{2.126}
\end{align}
を得たのと同様の手続きで，次の最尤推定の式から，1変数ガウス分布の分散の逐次推定の式を導出せよ．
\sigma_{\mathrm{ML}}^{2}=\frac{1}{N} \sum_{n=1}^{N}\left(x_{n}-\mu\right)^{2} \tag{2.292}
Robbins-Monroの逐次推定の式
\theta^{(N)}=\theta^{(N-1)}-a_{N-1} \frac{\partial}{\partial \theta^{(N-1)}}\left[-\ln p\left(x_{N} \mid \theta^{(N-1)}\right)\right] \tag{2.135}
にガウス分布を代入すると，これと同じ式になることを確かめ，ここから対応する係数a_Nの式を求めよ．
※ この問題の主題は(2.292)式から(2.126)式のような手続きによって1変数ガウス分布の分散を求めるというのが前半で、後半はこの式とRobbins-Monroアルゴリズムによる(2.135)式が係数a_Nを適切に定めることで同型になることを示すことである。
まずは(2.292)式から分散の逐次更新式を求める。
\begin{aligned}
\sigma^2_{(N)} &=\frac{1}{N} \sum_{n=1}^{N}\left(x_{n}-\mu\right)^{2} \\
&=\frac{1}{N}\left\{\sum_{n=1}^{N-1}\left(x_{n}-\mu\right)^{2}+\left(x_{N}-\mu\right)^{2}\right\} \\
&=\frac{N-1}{N} \cdot \frac{1}{N-1} \sum_{n=1}^{N-1}\left(x_{n}-\mu \right)^{2}+\frac{1}{N}\left(x_{N}-\mu\right)^{2} \\
&=\frac{N-1}{N} \sigma^2_{(N-1)}+\frac{1}{N}\left(x_{N}-\mu \right)^{2} \\
&=\sigma^2_{(N-1)}+\frac{1}{N}\left( (x_N-\mu)^2-\sigma^2_{(N-1)}\right) \hspace{2em} \cdots (1)
\end{aligned}
一方、Robbins-Monroの式ガウス分布の式を代入すると
\begin{aligned}
\sigma_{(N)}^{2}
&=\sigma_{(N-1)}^{2}-a_{N-1} \frac{\partial}{\partial \sigma_{(N-1)}^{2}}\left[-\ln \left\{\frac{1}{\sqrt{2 \pi \sigma_{(N-1)}^{2}}} \exp \left(-\frac{\left(x_{N}-\mu\right)^{2}}{2 \sigma_{(N-1)}^{2}}\right)\right\}\right]\\
&=\sigma_{(N-1)}^{2}-a_{N-1} \frac{\partial}{\partial \sigma_{(N-1)}^{2}}\left[\frac{1}{2} \ln \left(2 \pi \sigma_{(N-1)}^{2}\right)+\frac{\left(x_{N}-\mu\right)^{2}}{2 \sigma_{(N-1)}^{2}}\right] \\
&=\sigma_{(N-1)}^{2}-a_{N-1}\left\{\frac{1}{2} \cdot \frac{1}{\sigma_{(N-1)}^{2}}-\frac{\left(x_{N}-\mu\right)^{2}}{2 \sigma_{(N-1)}^{4}}\right\} \\
&=\sigma_{(N-1)}^{2}+\frac{a_{N-1}}{2 \sigma_{(N-1)}^{4}}\left\{-\sigma^2_{(N-1)}+\left(x_{N}-\mu\right)^{2}\right\} \hspace{2em} \cdots (2)
\end{aligned}
(1), (2)式を比較すると\displaystyle a_{N-1} = \frac{2\sigma^4_{(N-1)}}{N} \hspace{1em} \left( a_{N} = \frac{2\sigma^4_{(N)}}{N+1} \right)が得られる。
演習 2.37
(2.126)を得たのと同様の手続きで、最尤推定の式
\mathbf{\Sigma}_{\mathbf{ML}}=\frac{1}{N} \sum_{n=1}^{N}\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{\mathbf{ML}}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{\mathbf{ML}}\right)^{\mathrm{T}} \tag{2.122}
から，多変量ガウス分布の共分散の逐次推定の式を導出せよ．Robbins-Monroの逐次推定の式(2.135)にガウス分布を代入すると，これと同じ式になることを確かめ，ここから係数a_Nの式を求めよ．
※問題設定が悪く、解けない可能性が高い。
(2.122)式
\Sigma_{\mathrm{ML}} = \frac{1}{N}\sum_{n=1}^{N}\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{\mathrm{ML}}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{\mathrm{ML}}\right)^{\mathrm{T}}
から多次元ガウス分布の分散の逐次更新式を求める。
\begin{aligned}
\Sigma_{(N)} &= \frac{1}{N}\sum_{n=1}^{N}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)^{\mathrm{T}} \\
&= \frac{N-1}{N}\cdot\frac{1}{N-1}\left\{ \sum_{n=1}^{N-1}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)^{\mathrm{T}}+\left(\mathbf{x}_{N}-\boldsymbol{\mu}\right)\left(\mathbf{x}_{N}-\boldsymbol{\mu}\right)^{\mathrm{T}} \right\} \\
&= \frac{N-1}{N}\Sigma_{(N-1)}+\frac{1}{N}(\mathbf{x}_{N}-\boldsymbol{\mu})(\mathbf{x}_{N}-\boldsymbol{\mu})^{\mathrm{T}} \\
&= \Sigma_{(N-1)}+\frac{1}{N}\left\{ (\mathbf{x}_{N}-\boldsymbol{\mu})(\mathbf{x}_{N}-\boldsymbol{\mu})^{\mathrm{T}} - \Sigma_{(N-1)}\right\} \hspace{2em} \cdots (1)
\end{aligned}
一方、Robbins-Monroの式に多変量ガウス分布を当てはめると、
\Sigma_{(N)} = \Sigma_{(N-1)}-a_{N-1}\frac{\partial}{\partial \Sigma_{(N-1)}}\left[ -\ln \left\{ \frac{1}{(2\pi)^{D/2}}\frac{1}{|\Sigma_{(N-1)}|^{1/2}}\exp\left\{ -\frac{1}{2}(\mathbf{x}_N-\boldsymbol{\mu})^{\mathrm{T}}\Sigma_{(N-1)}(\mathbf{x}_N-\boldsymbol{\mu}) \right\} \right\} \right]\hspace{2em} \cdots (2)
この微分項について考える。
\begin{aligned}
&\frac{\partial}{\partial \Sigma_{(N-1)}}\left[ -\ln \left\{ \frac{1}{(2\pi)^{D/2}}\frac{1}{|\Sigma_{(N-1)}|^{1/2}}\exp\left\{ -\frac{1}{2}(\mathbf{x}_N-\boldsymbol{\mu})^{\mathrm{T}}\Sigma_{(N-1)}(\mathbf{x}_N-\boldsymbol{\mu}) \right\} \right\} \right] \\
=& \frac{1}{2}\Sigma_{(N-1)}^{-1}-\frac{1}{2}\left\{ \Sigma_{(N-1)}^{-1}(\mathbf{x}_{N}-\boldsymbol{\mu})(\mathbf{x}_{N}-\boldsymbol{\mu})^{\mathrm{T}}\Sigma_{(N-1)}^{-1} \right\} \\
=& -\frac{1}{2}\Sigma_{(N-1)}^{-1}\left\{ (\mathbf{x}_{N}-\boldsymbol{\mu})(\mathbf{x}_{N}-\boldsymbol{\mu})^{\mathrm{T}} -\Sigma_{(N-1)} \right\}\Sigma_{(N-1)}^{-1}
\end{aligned}
(2)式に当てはめれば
\Sigma_{(N)} = \Sigma_{(N-1)}+\frac{a_{N-1}\Sigma_{(N-1)}^{-1}}{2}\left\{ (\mathbf{x}_{N}-\boldsymbol{\mu})(\mathbf{x}_{N}-\boldsymbol{\mu})^{\mathrm{T}} -\Sigma_{(N-1)} \right\}\Sigma_{(N-1)}^{-1}
ここまではわかったけれど係数a_{N-1}がきれいな式にならない……。解答例では、追加で共分散行列がdiagonal（対角行列？）ならば\displaystyle a_{N-1} = \frac{2}{N}\left( \Sigma_{(N-1)}\right)^2とすれば良いと書いてあるが、仮にdiagonalでも合わない気がする。Googleで調べてみたけれどよくわからない。
参考：https://math.stackexchange.com/questions/1558016/deriving-mle-for-covariance-matrix-using-robbins-monro
演習 2.38
二次形式を平方完成することで，
\mu_{N} =\frac{\sigma^{2}}{N \sigma_{0}^{2}+\sigma^{2}} \mu_{0}+\frac{N \sigma_{0}^{2}}{N \sigma_{0}^{2}+\sigma^{2}} \mu_{\mathrm{ML}} \tag{2.141}
\frac{1}{\sigma_{N}^{2}} =\frac{1}{\sigma_{0}^{2}}+\frac{N}{\sigma^{2}} \tag{2.142}
の結果を導出せよ．
p(\mu \mid \mathbf{x}) \propto p(\mathbf{x} \mid \mu) p(\mu) \tag{2.139}
について、
p(\mathbf{x} \mid \mu)=\prod_{n=1}^{N} p\left(x_{n} \mid \mu\right)=\frac{1}{\left(2 \pi \sigma^{2}\right)^{N / 2}} \exp \left\{-\frac{1}{2 \sigma^{2}} \sum_{n=1}^{N}\left(x_{n}-\mu\right)^{2}\right\} \tag{2.137}
とp(\mu) = \mathcal{N}(\mu | \mu_0, \sigma_0^2) \hspace{1em}(2.138)式で変形する。
\muについての式としてのp(\mathbf{X}|\mu)p(\mu)の指数部分について考えると
\begin{aligned}
& \exp \left\{ -\frac{1}{2\sigma^2}\sum_{n=1}^{N}(x_n-\mu)^2 \right\}\cdot \exp \left\{ -\frac{1}{2\sigma_0^2}(\mu-\mu_0)^2 \right\} \\
=&\exp \left \{ -\frac{1}{2\sigma^2} \sum_{n=1}^{N}(x_n^2 -2x_n \mu + \mu^2) - \frac{1}{2\sigma_0^2} (\mu^2 - 2\mu\mu_0 + \mu_0^2)\right\} \\
=&\exp \left\{ -\frac{1}{2}\left( \frac{N}{\sigma^2}+\frac{1}{\sigma_0^2}\right)\mu^2 + \left( \frac{1}{\sigma^2}\sum_{n=1}^{N}x_n+\frac{\mu_0}{\sigma_0^2} \right)\mu + \cdots \right\}
\end{aligned}
となる。ここで、\cdotsの部分は\muに依存しない項なので考えなくても良い。一方で、(2.140)の右辺の指数部分について考えると、
\begin{aligned}
& \exp \left\{ -\frac{1}{2\sigma_N^2}(\mu - \mu_N)^2 \right\} \\
=& \exp \left\{ -\frac{1}{2}\left( \frac{1}{\sigma_N^2} \right)\mu^2+\frac{\mu_N}{\sigma_N^2}\mu+\cdots \right\}
\end{aligned}
\displaystyle \sum_{n=1}^{N}x_n = N\mu_{\mathrm{ML}}に注意して、これら2つの式を比較すると
\begin{aligned}
\frac{1}{\sigma_N^2}&=\frac{1}{\sigma_0^2}+\frac{N}{\sigma^2},\hspace{2em}\left( \sigma_N^2 = \frac{\sigma^2\sigma_0^2}{\sigma^2+N\sigma_0^2} \right) \\
\frac{\mu_N}{\sigma_N^2}&=\frac{\mu_0}{\sigma_0^2}+\frac{N}{\sigma^2}\mu_{\mathrm{ML}}
\end{aligned}
2式めを変形して
\begin{aligned}
\mu_{N}&=\left( \frac{\mu_0}{\sigma_0^2}+\frac{N\mu_{\mathrm{ML}}}{\sigma^2} \right)\sigma_{N}^2 \\
&=\left( \frac{\mu_0}{\sigma_0^2}+\frac{N\mu_{\mathrm{ML}}}{\sigma^2} \right) \left( \frac{\sigma^2\sigma_0^2}{\sigma^2+N\sigma_0^2} \right) \\
&=\frac{\sigma^{2}}{N \sigma_{0}^{2}+\sigma^{2}} \mu_{0}+\frac{N \sigma_{0}^{2}}{N \sigma_{0}^{2}+\sigma^{2}} \mu_{\mathrm{ML}}
\end{aligned}
ベイズ推定では、事前知識を決め打ちではなく分布の形で表す。
演習 2.39
ガウス確率変数の平均の事後分布についての結果
\mu_{N} =\frac{\sigma^{2}}{N \sigma_{0}^{2}+\sigma^{2}} \mu_{0}+\frac{N \sigma_{0}^{2}}{N \sigma_{0}^{2}+\sigma^{2}} \mu_{\mathrm{ML}} \tag{2.141}
\frac{1}{\sigma_{N}^{2}} =\frac{1}{\sigma_{0}^{2}}+\frac{N}{\sigma^{2}} \tag{2.142}
を元に，最初のN-1個のデータ点の影響を分離し\mu_{N}と\sigma^2_{N}の逐次更新の式を求めよ．そして，事後分布p\left(\mu | x_{1}, \ldots, x_{N-1}\right)=\mathcal{N}\left(\mu | \mu_{N-1}, \sigma_{N-1}^{2}\right)に尤度関数p(x_N|\mu) = \mathcal{N}(x_N|\mu, \sigma^2)を掛けた後,平方完成と正規化をすることで，N個の観測値を得た後の事後分布と同じ結果を導出せよ．
問題文の通り、N-1個のデータ点との影響を分離する。ただし、分散\sigma^2は既知であることに注意する。
(2.141)と(2.142)式と、演習問題2.38の結果から、\mu_0 \to \mu_{N-1}, \sigma_0 \to \sigma_{N-1}, \mu_{\mathrm{ML}}=x_N, N=1の場合に当てはまるので、\displaystyle \mu_N = \left( \frac{\mu_{N-1}}{\sigma_{N-1}^2}+\frac{x_N}{\sigma^2} \right)\sigma_N^2の形になるはずである。これを示す。
\begin{aligned}
\frac{1}{\sigma_N^2}&=\frac{1}{\sigma_0^2}+\frac{N}{\sigma^2}=\frac{1}{\sigma_0^2}+\frac{N-1}{\sigma^2}+\frac{1}{\sigma^2} = \frac{1}{\sigma_{N-1}^2}+\frac{1}{\sigma^2} \\
\sigma_N^2 &= \frac{\sigma^2+\sigma_{N-1}^2}{\sigma^2\sigma_{N-1}^2}
\end{aligned}
平均\mu_Nについて
\begin{aligned}
\mu_N &= \frac{\sigma^{2}}{N \sigma_{0}^{2}+\sigma^{2}} \mu_{0}+\frac{N \sigma_{0}^{2}}{N \sigma_{0}^{2}+\sigma^{2}} \mu_{\mathrm{ML}} \\
&= \frac{\mu_0}{\sigma_0^2}\sigma_N^2+\frac{\sigma_N^2}{\sigma^2}\left( \sum_{n=1}^{N-1}x_n + x_N \right) \\
&= \left( \frac{\mu_0}{\sigma_0^2}+\frac{1}{\sigma^2}\sum_{n=1}^{N-1}x_n \right)\sigma_N^2 + \frac{x_N}{\sigma^2}\sigma_N^2
\end{aligned}
ここで\displaystyle \frac{\mu_0}{\sigma_0^2}+\frac{1}{\sigma^2}\sum_{n=1}^{N-1}x_nについて、演習2.38の途中式で出てきたように
\frac{\mu_0}{\sigma_0^2}+\frac{1}{\sigma^2}\sum_{n=1}^{N-1}x_n = \frac{\mu_{N-1}}{\sigma_{N-1}^2}
と書けるので、
\mu_N = \left( \frac{\mu_{N-1}}{\sigma_{N-1}^2} + \frac{x_N}{\sigma^2} \right)\sigma_N^2\, \hspace{2em}\left( \frac{\mu_N}{\sigma_N^2} = \frac{\mu_{N-1}}{\sigma_{N-1}^2}+\frac{x_N}{\sigma^2} \right)
となる。これが逐次更新の式となる。
事後分布に尤度関数をかけて、平方完成と正規化〜の部分は演習問題2.38の解き方と同じであり、得られる結果は上記となる。
演習 2.40
D次元ガウス確率変数\mathbf{x}を考える．この分布\mathcal{N}(\mathbf{x}|\boldsymbol{\mu}, \mathbf{\Sigma})の共分散\mathbf{\Sigma}は既知としたとき，観測値集合\mathbf{X} = \{ \mathbf{x}_1, \ldots, \mathbf{x}_N \}から平均\boldsymbol{\mu}を推定したいとする．事前分布p(\boldsymbol{\mu})=\mathcal{N}\left(\boldsymbol{\mu} | \boldsymbol{\mu}_{0}, \mathbf{\Sigma_{0}} \right)について、これに対応する事後分布p(\boldsymbol{\mu} | \mathbf{X})を求めよ．
事後分布p(\boldsymbol{\mu} | \mathbf{X})は事前分布\prod_{n=1}^{N} p\left(\mathbf{x}_{n} | \boldsymbol{\mu}, \mathbf{\Sigma}\right)に尤度関数p(\boldsymbol{\mu})をかけたものに比例するので、
p(\boldsymbol{\mu} | \mathbf{X}) \propto \prod_{n=1}^{N} p\left(\mathbf{x}_{n} | \boldsymbol{\mu}, \mathbf{\Sigma}\right) p(\boldsymbol{\mu})
先の演習問題と同様に指数部分を考えると
\begin{aligned}
&-\frac{1}{2}\left(\boldsymbol{\mu}-\boldsymbol{\mu}_{0}\right)^{\mathrm{T}} \mathbf{\Sigma}_{0}^{-1}\left(\boldsymbol{\mu}-\boldsymbol{\mu}_{0}\right)-\frac{1}{2} \sum_{n=1}^{N}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right)^{\mathrm{T}} \mathbf{\Sigma}^{-1}\left(\mathbf{x}_{n}-\boldsymbol{\mu}\right) \\
=&-\frac{1}{2} \boldsymbol{\mu}^{\mathrm{T}}\left(\mathbf{\Sigma}_{0}^{-1}+N \mathbf{\Sigma}^{-1}\right) \boldsymbol{\mu}+\boldsymbol{\mu}^{\mathbf{T}}\left(\mathbf{\Sigma}_{0}^{-1} \boldsymbol{\mu}_{0}+\mathbf{\Sigma}^{-1} \sum_{n=1}^{N} \mathbf{x}_{n}\right)+\mathrm{const.}
\end{aligned}
ここで、\mathrm{const.}は\boldsymbol{\mu}と独立な項であり、正規化時に吸収される。
求める事後分布の形を\mathcal{N}(\boldsymbol{\mu}|\boldsymbol{\mu}_N,\mathbf{\Sigma}_{N})とすると、この指数部分は
-\frac{1}{2}(\boldsymbol{\mu}-\boldsymbol{\mu}_N)^{\mathrm{T}}\mathbf{\Sigma}_N^{-1}(\boldsymbol{\mu}-\boldsymbol{\mu}_N)=-\frac{1}{2}\boldsymbol{\mu}^{\mathrm{T}}\mathbf{\Sigma}_N^{-1}\boldsymbol{\mu}+\boldsymbol{\mu}^{\mathrm{T}}\mathbf{\Sigma}_N^{-1}\boldsymbol{\mu}_N+\textrm{const.}
となるので、まず\boldsymbol{\mu}の二次形式部分について
\mathbf{\Sigma}_{N}^{-1} =\mathbf{\Sigma}_{0}^{-1}+N \mathbf{\Sigma}^{-1}
となり、\boldsymbol{\mu}^{\mathrm{T}}の係数について
\begin{aligned}
\boldsymbol{\mu}_{N} &=\mathbf{\Sigma}_N\left(\mathbf{\Sigma}_{0}^{-1} \boldsymbol{\mu}_{0}+\mathbf{\Sigma}^{-1} N \boldsymbol{\mu}_{\mathrm{ML}}\right) \\
&=\left(\mathbf{\Sigma}_{0}^{-1}+N \boldsymbol{\Sigma}^{-1}\right)^{-1}\left(\mathbf{\Sigma}_{0}^{-1} \boldsymbol{\mu}_{0}+\mathbf{\Sigma}^{-1} N \boldsymbol{\mu}_{\mathrm{ML}}\right)
\end{aligned}
となる。ただし、ここで\displaystyle \sum_{n=1}^{N}\mathbf{x}_n=N\boldsymbol{\mu}_{\mathrm{ML}}とした。
すなわち、事後分布p(\boldsymbol{\mu} | \mathbf{X})は\mathcal{N}(\boldsymbol{\mu}|\boldsymbol{\mu}_N,\mathbf{\Sigma}_{N}) = \mathcal{N}(\boldsymbol{\mu}\mid \left(\mathbf{\Sigma}_{0}^{-1}+N \boldsymbol{\Sigma}^{-1}\right)^{-1}\left(\mathbf{\Sigma}_{0}^{-1} \boldsymbol{\mu}_{0}+\mathbf{\Sigma}^{-1} N \boldsymbol{\mu}_{\mathrm{ML}}\right), (\mathbf{\Sigma}_{0}^{-1}+N \mathbf{\Sigma}^{-1})^{-1})の形で書ける。
Discussion