はじめに
PRML解答例まとめを参照
演習 9.1
9.1節で議論したK-meansアルゴリズムを考える．離散的な支持変数r_{nk}の状態が有限通りしかないこと，そしてそのような状態の各々について，\{\boldsymbol{\mu}_k\}のただ1つの最適値があることの帰結としてK-meansアルゴリズムは有限回の繰り返しで収束することを示せ．
※教科書pp.140–142参照。
教科書にあるように、K-meansアルゴリズムはEステップとMステップに分かれている。Eステップではデータ点\mathbf{x}_nと\boldsymbol{\mu}_kは固定されている。各nについて、r_{nk}=1としたときに\|\mathbf{x}_n - \boldsymbol{\mu}_k\|^2が最小となるようなk\ (k=1\ldots K)を選んでr_{nk}=1とし、残りを0とする。式で書くと
r_{n k}=\left\{\begin{array}{ll}1 & k=\arg \min _{j}\left\|\mathrm{x}_{n}-\mu_{j}\right\|^{2} \text { のとき } \\ 0 & \text { それ以外. }\end{array}\right. \tag{9.2}
言い換えると、各nについてr_{nk}=1の選び方はK通り存在するので、N個のデータ点について言えば全部でK^N通りである。一方で、Eステップの操作によって必ず(9.1)式の歪み尺度(distortion measure)Jは小さくなる(図9.2を参照)。
Mステップについてはr_{nk}を固定した状況下で歪み尺度Jを最小化する。このとき最適な\boldsymbol{\mu}_kは(9.3)式にしたがって解析的に求まる。言い換えれば、前回のMステップからr_{nk}の割り当てが変わっている場合、Mステップでの\boldsymbol{\mu}_kの更新操作で必ずJは低下するが、Eステップでのr_{nk}の割り当てが変わらない場合はMステップでの\boldsymbol{\mu}_kも変化しない。このときK-meansアルゴリズムは収束解を得る。
したがって、すべてのデータ点\mathbf{x}_nについて最適なクラス割り当てが見つかるまでEステップとMステップが繰り返されるわけだが、これは高々K^N回の繰り返しの中で必ず収束することとなる。（※ただし、この収束解は大域的な最適解になるとは限らないことに注意。）
演習 9.2
2.3.5節で述べた Robbins-Monroの逐次的推定法を
J=\sum_{n=1}^{N} \sum_{k=1}^{K} r_{n k}\left\|\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right\|^{2} \tag{9.1}
で与えたJの\boldsymbol{\mu}_kによる微分で得られる回帰関数の根を求める問題に適用せよ．また，これから各データ点\mathbf{x}_nについて，最も近いプロトタイプ\boldsymbol{\mu}_kが
\boldsymbol{\mu}_{k}^{\text {new }}=\boldsymbol{\mu}_{k}^{\text {old }}+\eta_{n}\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}^{\text {old }}\right) \tag{9.5}
によって更新されるような確率的K-meansアルゴリズムが導かれることを示せ．
\frac{\partial J}{\partial \boldsymbol{\mu} _k}=0となるように、Jに含まれるパラメータr_{nk}, \boldsymbol{\mu}_kを最適化する問題を、データ系列の数Nについて逐次的に解く。
データ点\{ \mathbf{x}_1, \mathbf{x}_2,\cdots \mathbf{x}_{N-1} \}に基づいてパラメータ\{r_{nk}, \boldsymbol{\mu}_k\}_{k=1,2,\cdots,K, n=1,2,\cdots, N-1}が最適化されている（つまり、\left(\frac{\partial J}{\partial \boldsymbol{\mu} _k}\right)^{(N-1)}=0となっている）とき、データ点\mathbf{x}_Nを追加で入手した時にパラメータ\{r_{nk}, \boldsymbol{\mu}_k\}_{k=1,2,\cdots,K, n=1,2,\cdots, N}を最適化する（つまり、\left(\frac{\partial J}{\partial \boldsymbol{\mu} _k}\right)^{(N)}=0とする）ことを考える。
\begin{aligned}
\left( \frac{\partial J}{\partial \boldsymbol{\mu} _k}
\right)^{(N)}
&=\frac{\partial}{\partial \boldsymbol{\mu} _k}
\left( \sum_{n=1}^{N} \sum_{k'=1}^{K} r_{n k'}\left\|\mathbf{x}_{n}-\boldsymbol{\mu}_{k'}\right\|^{2} \right) \\
&=\sum_{n=1}^{N} \sum_{k'=1}^{K} r_{n k'} (-2) (\mathbf{x}_{n}-\boldsymbol{\mu}_{k'}) \delta_{kk'} \\
&=-2 \sum_{n=1}^{N} r_{n k} (\mathbf{x}_{n}-\boldsymbol{\mu}_{k})\\
&=
\left( \frac{\partial J}{\partial \boldsymbol{\mu} _k}
\right)^{(N-1)}
-2 r_{N k} (\mathbf{x}_{N}-\boldsymbol{\mu}_{k})
\end{aligned}
\mathbf{x}_{N}に最も近い\{\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2}\cdots\boldsymbol{\mu}_{K} \}が\boldsymbol{\mu}_{k}と異なる場合、Eステップにおいてr_{Nk}=0とする。この場合、Mステップにおいて\{\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2}\cdots\boldsymbol{\mu}_{K} \}の値は更新されない。
\mathbf{x}_{N}に最も近い\{\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2}\cdots\boldsymbol{\mu}_{K} \}が\boldsymbol{\mu}_{k}に一致している場合、Eステップにおいてr_{Nk}=1とする。この場合、Mステップにおける更新式（上式）は、
\begin{aligned}
\left( \frac{\partial J}{\partial \boldsymbol{\mu} _k}
\right)^{(N)}
&=
\left( \frac{\partial J}{\partial \boldsymbol{\mu} _k}
\right)^{(N-1)}
-2 (\mathbf{x}_{N}-\boldsymbol{\mu}_{k})
\end{aligned}
となる。
Robbins-Monroの逐次推定の(2.129)式
\theta^{(N)}=\theta^{(N-1)}-a_{N-1} z\left(\theta^{(N-1)}\right)
に代入して、
\begin{aligned}
\boldsymbol{\mu}_k ^{(N)}
= \boldsymbol{\mu}_k ^{(N-1)} + 2 a_{N-1} (\mathbf{x}_{N}-\boldsymbol{\mu}_{k}^{(N-1)})
\end{aligned}
となる。
2 a_{N-1} = \eta_Nと変数を変換すれば、(9.5)式を得る。
演習 9.3
混合ガウスモデルを考え.潜在変数の周辺分布p(\mathbf{z})が
p(\mathbf{z})=\prod_{k=1}^{K} \pi_{k}^{z_{k}} \tag{9.10}
で，観測変数の条件付き分布
p(\mathbf{x}\mid\mathbf{z})が
p(\mathbf{x} \mid \mathbf{z})=\prod_{k=1}^{K} \mathcal{N}\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)^{z_{k}} \tag{9.11}
でそれぞれ与えられると仮定する．
p(\mathbf{z})p(\mathbf{x}\mid\mathbf{z})を
\mathbf{z}の可能な値について足して得られる周辺分布
p(\mathbf{x})が
p(\mathbf{x})=\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x} \mid \mu_{k}, \mathbf{\Sigma}_{k}\right) \tag{9.7}
の形の混合ガウス分布になることを示せ．
導入に従ってとけばよい。
\begin{aligned}
p(\mathbf{x}) &= \sum_{\mathbf{z}} p(\mathbf{z})p(\mathbf{x}|\mathbf{z}) \\
&=\sum_{\mathbf{z}}\prod_{k} \pi_k^{z_k} \prod_k N(\mathbf{x}|\mathbf{\mu}_k, \mathbf{\Sigma_k})^{z_k} &(\because (9.10), (9.11)))\\
&=\sum_{\mathbf{z}} \prod_{k} \pi_k^{z_k} N(\mathbf{x}|\mathbf{\mu}_k, \mathbf{\Sigma_k})^{z_k}\\
&= \sum_{\mathbf{z}}\prod_{k}( \pi_k N(\mathbf{x}|\mathbf{\mu}_k, \mathbf{\Sigma_k}))^{z_k}\\
&= \sum_{k} \pi_k N(\mathbf{x}|\mathbf{\mu}_k, \mathbf{\Sigma_k}) &(\because \text{1-of-K})
\end{aligned}
演習 9.4
潜在変数を持つモデルに関して，データ集合\mathbf{X}を観測した下での事後分布p(\boldsymbol{\theta}\mid \mathbf{X})を，EMアルゴリスムを用いて\boldsymbol{\theta}について最大化する問題を考える．このとき，Eステップは最尤推定問題の場合と同じであるのに対し，Mステップでは，最大化すべき量が\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }})+\ln p(\boldsymbol{\theta})で与えられることを示せ．ただし，\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }})は
\mathcal{Q}\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}\right)=\sum_{\mathbf{Z}} p\left(\mathbf{Z} \mid \mathbf{X}, \theta^{\text {old }}\right) \ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta}) \tag{9.30}
で定義されている．
ベイズの定理より
p(\mathbf{\theta}|\mathbf{X})=\frac{p(\mathbf{X}|\mathbf{\theta})p(\mathbf{\theta})}{p(\mathbf{X})}
である．対数尤度を取って，
\begin{aligned}
\ln p(\mathbf{\theta}|\mathbf{X})&= \ln p(\mathbf{X}|\mathbf{\theta}) + \ln p(\mathbf{\theta}) + \mathrm{const.}\\
&=\ln\sum_{\mathbf{Z}} p(\mathbf{X},\mathbf{Z}|\theta) + \ln p(\mathbf{\theta}) + \mathrm{const.}\\
\end{aligned}
ここで第一項については教科書の(9.29), (9.30)の議論と同様に対数の中の総和への処理，および潜在変数\mathbf{Z}についての事後確率の期待値を考えるといった工夫をすることで\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }})を用いて表すことができて，
\begin{aligned}
\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}) + \ln p(\mathbf{\theta})
\end{aligned}
を最大化すればいいことがわかる。
演習 9.5
図9.6で示すような，混合ガウスモデルの有向グラフを考える．8.2節で議論した有向分離の規準を利用して，潜在変数の事後分布が次式のように各データ点ごとの事後分布の積になることを示せ．
p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\mu}, \mathbf{\Sigma}, \boldsymbol{\pi})=\prod_{n=1}^{N} p\left(\mathbf{z}_{n} \mid \mathbf{x}_{n}, \boldsymbol{\mu}, \mathbf{\Sigma}, \boldsymbol{\pi}\right) \tag{9.80}
まず、図9.9より、\mathbf{\pi}と\mathbf{\mu}、\mathbf{\mu} と\mathbf{\Sigma}、\mathbf{\Sigma}と\mathbf{\pi}は、tail-to-tailなので、\mathbf{Z}と\mathbf{X}が観察された時、遮断されている。すなわち条件付き独立が成り立つ。よって、以下では、\mathbf{X}と\mathbf{Z}について考えていく。
今、図9.9より、n \in \{1, 2, \cdots N\}について、それぞれ\mathbf{z}_n \rightarrow \mathbf{x}_nが成り立つ。また、i.i.dデータ集合であることが仮定されているので、p(\mathbf{X}) = \prod_n p(\mathbf{x}_n)である。
以上から、
\begin{aligned}
p(\mathbf{Z}| \mathbf{X}, \boldsymbol{\mu}, \mathbf{\Sigma}, \boldsymbol{\pi}) &=
p(\mathbf{Z}|\mathbf{X}) \\
&= \prod_n p(\mathbf{z}_n|\mathbf{x}_n)
\end{aligned}
が成り立つ。最後に、条件付き独立であるパラメータを加えて、
p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\mu}, \mathbf{\Sigma}, \boldsymbol{\pi})=\prod_{n=1}^{N} p\left(\mathbf{z}_{n} \mid \mathbf{x}_{n}, \boldsymbol{\mu}, \mathbf{\Sigma}, \boldsymbol{\pi}\right) \tag{9.80}
が導かれる。
演習 9.6
混合ガウスモデルについて，各混合要素の共分散行列\mathbf{\Sigma}_kすべてが共通の値\mathbf{\Sigma}に制限された特別な場合を考える．そのようなモデルにおいて，尤度関数を最大化するEM方程式を導け．
(9.17),(9.18),(9.19),(9.22)に\mathbf{\Sigma}_{k}=\mathbf{\Sigma}_{j}=\mathbf{\Sigma}を代入する。
\begin{aligned}
&N_{k}=\sum_{n=1}^{N} \gamma\left(z_{n k}\right) \\
&\left(\gamma\left(z_{n k}\right)=\frac{\pi_{k} \mathcal{N}\left(X_{n} \mid \boldsymbol{\mu}_{k}, \mathbf{\Sigma}\right)}{\sum_{j} \pi_j \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{j}, \mathbf{\Sigma}\right)}\right) \\
&\boldsymbol{\mu}_{k}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right) \mathbf{x}_{n} \\
&\mathbf{\Sigma}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)^{\mathrm{T}} \\
&\pi_{k}=\frac{N_{k}}{N}
\end{aligned}
演習 9.7
混合ガウスモデルに関する完全データ対数尤度関数
\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})=\sum_{n=1}^{N} \sum_{k=1}^{K} z_{n k}\left\{\ln \pi_{k}+\ln \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)\right\} \tag{9.36}
を最大化することは，各混合要素について独立に平均と共分散をその混合要素に対応するデータ点集合にあてはめ，かつ混合係数を各グループに属するデータの数の割合に一致させるという結果になることを示せ．
(9.36)を\mathbf{\mu}_{k}で微分すると、
\begin{aligned}
\frac{\partial}{\partial\mathbf{\mu}_k} \ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \mathbf{\Sigma}, \boldsymbol{\pi}) &= \sum_{n=1}^{N} \sum_{k=1}^{K} \frac{\partial}{\partial\mathbf{\mu}_k}{ z_{n k}\left\{\ln \pi_{k}+\ln \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)\right\}}\\
&= \sum_{x_n \in C_k} \frac{\partial}{\partial\mathbf{\mu}_k} \{\ln \pi_{k} + \ln \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)\}\\
&= \sum_{x_n \in C_k} \frac{\partial}{\partial\mathbf{\mu}_k} \{\ln \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)\}
\end{aligned}
これは、「ガウス分布でクラスタC_kに属するデータ点を利用した場合の
パラメータ最適化」と同じである。
\pi_iについては、\sum_{k=1}^{K}\pi_k=1...(A)という制限があるため、ラグランジュの未定乗数法を用いて
L=\ln p(\mathbf{X},\mathbf{Z}|\mathbf{\mu},\mathbf{\Sigma},\mathbf{\pi}) + \lambda(\mathbf{\Sigma}_{k}\pi_{k}-1) ...\mathrm{(B)}
\frac{\partial L}{\partial \pi_{k}} = \sum_{n=1}^{N}\frac{z_{nk}}{\pi_k} + \lambda=0
両辺に\pi_kをかけて\mathbf{\Sigma}_{k}をとると、
\sum_{k=1}^K \sum_{n=1}^N z_{nk} + \sum_{k=1}^K \lambda\pi_k=0
制約条件(A),
\mathbf{\Sigma}_{k}\mathbf{\Sigma}_{n}z_{nk}=Nより、
\sum_{k=1}^K \sum_{n=1}^N z_{nk} + \lambda=0
N + \lambda = 0, \lambda = -N(B)
∴\pi_{k} = \frac{1}{N}\sum_{n=1}^N z_{nk}
演習 9.8
もし負担率\gamma(z_{nk})を固定した下で，
\mathbb{E}_{\mathbf{Z}}[\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})]=\sum_{n=1}^{N} \sum_{k=1}^{K} \gamma\left(z_{n k}\right)\left\{\ln \pi_{k}+\ln \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right\} \tag{9.40}
を
\boldsymbol{\mu}_kについての最大化をしようとすると，
(9.17)で与えられる陽な解が得られることを示せ．
\begin{aligned}
\mathbf{\mu}_k = \frac{1}{N_k}\sum_n \gamma(z_{nk})\mathbf{x}_n \tag{9.17}
\end{aligned}
ただし、
\begin{aligned}
N_k = \sum_n \gamma(z_{nk}) \tag{9/18}
\end{aligned}
(9.40)式を\mathbf{\mu}_kについて微分し、0となる\mathbf{\mu}_kを求めれば良い。
\begin{aligned}
&\frac{d}{d\mathbf{\mu}_k} \mathbb{E}_{\mathbf{Z}}[\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})] =
\sum_n \gamma (z_{nk})\Sigma^{-1}(\mathbf{x}_n-\mathbf{\mu}_k) = 0 \\
&\Rightarrow \sum_n \gamma (z_{nk})(\mathbf{x}_n-\mathbf{\mu}_k) = 0 \\
&\Rightarrow \mathbf{\mu}_k = \frac{1}{N_k}\sum_n \gamma(z_{nk})\mathbf{x}_n &(\because (9.18))
\end{aligned}
演習 9.9
もし負担率\gamma(z_{nk})を固定した下で，
\mathbb{E}_{\mathbf{Z}}[\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})]=\sum_{n=1}^{N} \sum_{k=1}^{K} \gamma\left(z_{n k}\right)\left\{\ln \pi_{k}+\ln \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right\} \tag{9.40}
を
\boldsymbol{\Sigma}_{k}と
\pi_kについての最大化をしようとすると，
\mathbf{\Sigma}_{k}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right)\left(\mathbf{x}_{n}-\mu_{k}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)^{\mathrm{T}} \tag{9.19}
と
\pi_k = \frac{N_k}{N} \tag{9.22}
で与えられる陽な解が得られることを示せ．
演習(9.8)と同様に最適化したい変数の微分を0として解くことで目的の解を得る．
\Sigma_k^{-1}について(9.40)を偏微分する
\begin{aligned}
\frac{\partial}{\partial\Sigma_k^{-1}}\mathbb{E}_{\mathbf{Z}}[\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})]&=\frac{\partial}{\partial\Sigma_k^{-1}}\sum_{n=1}^{N} \sum_{k=1}^{K} \gamma\left(z_{n k}\right)\left\{\ln \pi_{k}+\ln \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right\}\\
&=\frac{\partial}{\partial\Sigma_k^{-1}}\sum_{n=1}^{N} \gamma\left(z_{n k}\right)\ln \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\\
&=\frac{\partial}{\partial\Sigma_k^{-1}}\sum_{n=1}^{N} \gamma\left(z_{n k}\right)\left\{-\frac{|D|}{2}\ln 2\pi+\frac{1}{2}\ln |\Sigma_k|^{-1}-\frac{1}{2}(\mathbf{x}_n-\mathbf{\mu}_k)^T\Sigma_k^{-1}(\mathbf{x}_n-\mathbf{\mu}_k)\right\}\\
&=\frac{\partial}{\partial\Sigma_k^{-1}}\sum_{n=1}^{N} \frac{1}{2}\gamma\left(z_{n k}\right)\left\{\ln |\Sigma_k^{-1}|-(\mathbf{x}_n-\mathbf{\mu}_k)^T\Sigma_k^{-1}(\mathbf{x}_n-\mathbf{\mu}_k)\right\}\\
&=\sum_{n=1}^{N} \frac{1}{2}\gamma\left(z_{n k}\right)\left\{\Sigma_k^T-(\mathbf{x}_n-\mathbf{\mu}_k)(\mathbf{x}_n-\mathbf{\mu}_k)^T\right\}\\
&=\sum_{n=1}^{N} \frac{1}{2}\gamma\left(z_{n k}\right)\left\{\Sigma_k-(\mathbf{x}_n-\mathbf{\mu}_k)(\mathbf{x}_n-\mathbf{\mu}_k)^T\right\}\\
\end{aligned}
|A^{-1}|=|A|^{-1},\frac{\partial}{\partial{A}}\ln|A|=(A^{-1})^Tを用いた．
これを0として
\begin{aligned}
\sum_{n=1}^{N} \frac{1}{2}\gamma\left(z_{n k}\right)\left\{\mathbf{\Sigma}_k-(\mathbf{x}_n-\mathbf{\mu}_k)^T(\mathbf{x}_n-\mathbf{\mu}_k)\right\}=0
\end{aligned}
\begin{aligned}
\mathbf{\Sigma}_k\sum_{n=1}^{N} \gamma\left(z_{n k}\right)=\sum_{n=1}^{N} \gamma\left(z_{n k}\right)(\mathbf{x}_n-\mathbf{\mu}_k)^T(\mathbf{x}_n-\mathbf{\mu}_k)
\end{aligned}
\begin{aligned}
\mathbf{\Sigma}_k=\frac{1}{N}\sum_{n=1}^{N} \gamma\left(z_{n k}\right)(\mathbf{x}_n-\mathbf{\mu}_k)^T(\mathbf{x}_n-\mathbf{\mu}_k)
\end{aligned}
が得られる．
次に\pi_kについて考える．これは演習(9.7)を応用して考えることができて，\mathbf{z}_{nk}を\gamma(\mathbf{z}_{nk})で置き換え，同様にラグランジュ未定乗数法を用いて計算することで
\begin{aligned}
\boldsymbol{\pi}_k&=\frac{1}{N}\sum_{n=1}^N\gamma(\mathbf{z}_{nk})\\
&=\frac{N_k}{N}
\end{aligned}
以上から目的の解が得られた
演習 9.10
p(\mathbf{x})=\sum_{k=1}^{K} \pi_{k} p(\mathbf{x} \mid k) \tag{9.81}
の形の混合分布で与えられる密度のモデルを考え，ベクトル\mathbf{x}を\mathbf{x} = (\mathbf{x}_a, \mathbf{x}_b)のように2つの部分に分解すると仮定する．このとき，条件付き密度p(\mathbf{x}_b\mid \mathbf{x}_a)自体が混合分布であることを示し，混合係数と各混合要素の表式を求めよ．
(9.81)より
p\left(\mathbf{x}_{a}, \mathbf{x}_{b}\right)=\sum_{k=1}^{K} \pi_{k} p\left(\mathbf{x}_{a}, \mathbf{x}_{b} \mid k\right)
であることを用いて、
\begin{aligned}
p\left(\mathbf{x}_{b} \mid \mathbf{x}_{a}\right)&=\frac{p\left(\mathbf{x}_{a}, \mathbf{x}_{b}\right)}{p\left(\mathbf{x}_{a}\right)}
\\&=\frac{\sum_{k=1}^{K} \pi_{k} p\left(\mathbf{x}_{a}, \mathbf{x}_{b} \mid k\right)}{p\left(\mathbf{x}_{a}\right)}\\&=\frac{\sum_{k=1}^{K} \pi_{k} p\left(\mathbf{x}_{b} \mid \mathbf{x}_{a}, k\right) p\left(\mathbf{x}_{a} \mid k\right)}{p\left(\mathbf{x}_{a}\right)}
\\&=\sum_{k=1}^{K} \lambda_{k} p\left(\mathbf{x}_{b} \mid \mathbf{x}_{a}, k\right)
\end{aligned}
ただし混合係数は
\begin{aligned}
\lambda_{k} &\equiv \frac{\pi_{k} p\left(\mathbf{x}_{a} \mid k\right)}{p\left(\mathbf{x}_{a}\right)}\\&=\frac{\pi_{k} p\left(\mathbf{x}_{a} \mid k\right)}{\sum_{k=1}^{K} \pi_{k} p\left(\mathbf{x}_{a} \mid k\right)}
\end{aligned}
である。
演習 9.11
9.3.2節においてすべての混合要素の共分散が\epsilon\mathbf{I}である混合ガウス分布を考えることで，K-means法とEMアルゴリズムの関係を導いた．ここに，\epsilon \to 0の極限において，このモデルの期待完全データ対数尤度
\mathbb{E}_{\mathbf{Z}}[\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\pi})]=\sum_{n=1}^{N} \sum_{k=1}^{K} \gamma\left(z_{n k}\right)\left\{\ln \pi_{k}+\ln \mathcal{N}\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right\} \tag{9.40}
の最大化が，K-means法の歪み尺度
J=\sum_{n=1}^{N} \sum_{k=1}^{K} r_{n k}\left\|\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right\|^{2} \tag{9.1}
の最小化と等価になることを示せ．
式 (9.37) より\pi_kは常に非負。従って、\epsilon \to 0となるとき、\gamma(z_{nk}) \to r_{nk}である (p.160) ことから、式 (9.40) の最大化は
\begin{aligned}
\sum_{n = 1}^N \sum_{k = 1}^K r_{nk} \left(-\frac{1}{2\epsilon}\left\|\mathbf{x}_n - \mathbf{\mu}_k\right\|^2\right)
\end{aligned}
の最大化に等しい。上式に最適化の対象である\mathbf{\mu}_kを含まないscaling factor \left(-2\epsilon\right) をかけると式 (9.1) に等しくなることから、題意が示される。
演習 9.12
p(\mathbf{x})=\sum_{k=1}^{K} \pi_{k} p(\mathbf{x} \mid k) \tag{9.82}
の形の混合分布を考える．ここに，\mathbf{x}は離散，連続，その組み合わせのいずれでもよい．p(\mathbf{x}\mid k)の平均と共分散をそれぞれ\boldsymbol{\mu}_kと\mathbf{\Sigma}_kで表す．混合分布の平均と共分散がそれぞれ
\mathbb{E}[\mathbf{x}] =\sum_{k=1}^{K} \pi_{k} \boldsymbol{\mu}_{k} \tag{9.49}
\operatorname{cov}[\mathbf{x}] =\sum_{k=1}^{K} \pi_{k}\left\{\mathbf{\Sigma}_{k}+\boldsymbol{\mu}_{k} \boldsymbol{\mu}_{k}^{\mathrm{T}}\right\}-\mathbb{E}[\mathbf{x}] \mathbb{E}[\mathbf{x}]^{\mathrm{T}} \tag{9.50}
で与えられることを示せ．
問題設定より、p(\mathbf{x}\mid k)（の分布の下での\mathbf{x}）の平均と共分散がそれぞれ\boldsymbol{\mu}_k, \mathbf{\Sigma}_kとなるので、上巻pp.19-20, (1.36)の記法にしたがって
\mathbb{E}_{k}[\mathbf{x}]=\int \mathbf{x} p(\mathbf{x} \mid k) d \mathbf{x}=\boldsymbol{\mu}_{k},\ \operatorname{cov}_{k}[\mathbf{x}]=\mathbf{\Sigma}_{k}
と書ける。
\mathbb{E}[\mathbf{x}]について、期待値の定義から
\begin{aligned}
\mathbb{E}[\mathbf{x}] &=\int \mathbf{x} p(\mathbf{x}) d \mathbf{x} \\
&=\int \mathbf{x} \sum_{k=1}^{K} \pi_{k} p(\mathbf{x} \mid k) d \mathbf{x} \\
&=\sum_{k=1}^{K} \pi_{k} \int \mathbf{x} p(\mathbf{x} \mid k) d \mathbf{x} \\
&=\sum_{k=1}^{K} \pi_{k} \boldsymbol{\mu}_{k} \end{aligned}
となる。次に共分散\operatorname{cov}[\mathbf{x}]について、
\begin{aligned} \operatorname{cov}[\mathbf{x}]
&=\mathbb{E}\left[(\mathbf{x}-\mathbb{E}[\mathbf{x}])(\mathbf{x}-\mathbb{E}[\mathbf{x}])^{\mathrm{T}}\right] \quad(\because(1.42)) \\
&=\mathbb{E}\left[\mathbf{xx}^{\mathrm{T}}-2 \mathbb{E}[\mathbf{x}]\mathbb{E}[\mathbf{x}]^{\mathrm{T}}+\mathbb{E}[\mathbf{x}] \mathbb{E}[\mathbf{x}]^{\mathrm{T}}\right] \\
&=\mathbb{E}\left[\mathbf{xx}^{\mathrm{T}}\right]-\mathbb{E}[\mathbf{x}] \mathbb{E}[\mathbf{x}]^{\mathrm{T}} \\
&=\int \mathbf{xx}^{\mathrm{T}} \sum_{k=1}^{K} \pi_{k} p(\mathbf{x} \mid k) d x-\mathbb{E}[\mathbf{x}] \mathbb{E}[\mathbf{x}]^{\mathrm{T}} \\
&=\sum_{k=1}^{K} \pi_{k} \int \mathbf{xx}^{\mathrm{T}} p(\mathbf{x} \mid k) d x-\mathbb{E}[\mathbf{x}] \mathbb{E}[\mathbf{x}]^{\mathrm{T}} \\
&=\sum_{k=1}^{K} \pi_{k} \mathbb{E}_{k}\left[\mathbf{xx}^{\mathrm{T}}\right]-\mathbb{E}[\mathbf{x}] \mathbb{E}[\mathbf{x}]^{\mathrm{T}}
\end{aligned}
ここで\mathbb{E}_{k}[\mathbf{xx}^{\mathrm{T}}]について、(1.42)の共分散行列の定義から
\begin{aligned} \mathbf{\Sigma}_{k}=\operatorname{cov}_{k}[\mathbf{x}]
&=\mathbb{E}_{k}\left[\left(\mathbf{x}-\mathbb{E}_{k}[\mathbf{x}]\right)\left(\mathbf{x}-\mathbb{E}_{k}[\mathbf{x}]\right)^{\mathrm{T}}\right] \\
&=\mathbb{E}_{k}\left[\left(\mathbf{x}-\boldsymbol{\mu}_{k}\right)\left(\mathbf{x}-\boldsymbol{\mu}_{k}\right)^{\mathrm{T}}\right] \\
&=\mathbb{E}_{k}\left[\mathbf{xx}^{\mathrm{T}}\right]-\boldsymbol{\mu}_{k} \boldsymbol{\mu}_{k}^{\mathrm{T}} \end{aligned}
であるから、\mathbb{E}_{k}\left[\mathbf{xx}^{\mathrm{T}}\right] = \mathbf{\Sigma}_{k}+\boldsymbol{\mu}_{k} \boldsymbol{\mu}_{k}^{\mathrm{T}}が得られる。以上を代入して、
\operatorname{cov}[\mathbf{x}] =\sum_{k=1}^{K} \pi_{k}\left\{\mathbf{\Sigma}_{k}+\boldsymbol{\mu}_{k} \boldsymbol{\mu}_{k}^{\mathrm{T}}\right\}-\mathbb{E}[\mathbf{x}] \mathbb{E}[\mathbf{x}]^{\mathrm{T}}
となる。
演習 9.13
ベルヌーイ分布の混合は，パラメータを最尤解のときの値に一致させた場合
\mathbb{E}[\mathrm{x}]=\frac{1}{N} \sum_{n=1}^{N} \mathrm{x}_{n} \equiv \overline{\mathrm{x}} \tag{9.83}
を満たすことをEMアルゴリズムの更新式を用いて示せ．さらにこのモデルにおいて,すべての混合要素の平均値を\boldsymbol{\mu}_{k}=\widehat{\boldsymbol{\mu}}\ (k=1, \ldots, K)と同ーの値になるようにパラメータを初期化すると，混合係数の初期値が何であっても， EMアルゴリズムは一度の繰り返しで収束し，\boldsymbol{\mu}_k = \overline{\mathbf{x}}となることを示せ．これは，すべての混合要素が同ーになる縮退したケースを表しており，実用上，我々は適切な初期化によってこのような解を避ける努力をする．
\begin{aligned}
p(\mathbf{x}|\boldsymbol{\mu}, \boldsymbol{\pi}) = \sum_{k=1}^K \pi_k p(\mathbf{x}|\boldsymbol{\mu}_k ) \tag{9.47}
\end{aligned}
なので、
\begin{aligned}
\mathbb{E}[\mathbf{x}] &= \sum_{k=1}^K \pi_k \boldsymbol{\mu}_k \\
&= \sum_{k=1}^K \pi_k \overline{\mathbf{x}_k} \\
&= \sum_{k=1}^K \pi_k \frac{1}{N_k} \sum_{n=1}^N \gamma (z_{nk}) \mathbf{x}_n \\
&= \frac{1}{N} \sum_{n=1}^N \sum_{k=1}^K \gamma (z_{nk}) \mathbf{x}_n \\
&= \frac{1}{N} \sum_{n=1}^N \mathbf{x}_n \\
\end{aligned}
となる。ここで、\gamma (z_{nk})は負担率であることから、kについて和をとると1になることを用いた。
次に、後半部分を証明する。今、\boldsymbol{\mu}_{k} = \widehat{\boldsymbol{\mu}}\ (k=1, \ldots, K)と設定された時、Eステップの更新式は(9.56)式で表され、
\begin{aligned}
\gamma (z_{nk})
&= \frac{\pi _k p(\mathbf{x}_n | \boldsymbol{\mu}_k = \widehat{\boldsymbol{\mu}})}{\sum_{j=1}^K \pi _j p(\mathbf{x}_n | \boldsymbol{\mu}_j = \widehat{\boldsymbol{\mu}})} \\
&= \frac{\pi _k }{\sum_{j=1}^K \pi _j} \\
&= \pi_k \left(\because \sum_{j=1}^K \pi _j = 1 \right)
\end{aligned}
となり、nによらない。
Mステップの更新式は(9.59)式で表され、\gamma (z_{nk}) = \pi_{k}を用いると
\begin{aligned}
\boldsymbol{\mu}_k
&= \overline{x_k} \\
&= \frac{1}{N_k} \sum_{n=1}^N \gamma (z_{nk}) \mathbf{x}_n \\
&= \frac{1}{N_k} \pi_k \sum_{n=1}^N \mathbf{x}_n \\
&= \frac{1}{N} \sum_{n=1}^N \mathbf{x}_n \\
\end{aligned}
となり、kによらない。
従って、逐次的に更新式を適用しても2回目以降は値が変化しない。結果として、純粋な（混合分布ではない）ベルヌーイ分布と同じ確率分布しか表現できない。
演習 9.14
p(\mathbf{x} \mid \mathbf{z}, \boldsymbol{\mu})=\prod_{k=1}^{K} p\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}\right)^{z_{k}} \tag{9.52}
と
p(\mathbf{z} \mid \boldsymbol{\pi})=\prod_{k=1}^{K} \pi_{k}^{z_{k}} \tag{9.53}
の積でベルヌーイ分布の潜在変数と観測変数の同時分布を構成しよう．この同時分布を\mathbf{z}について周辺化すると
p(\mathbf{x} \mid \boldsymbol{\mu}, \boldsymbol{\pi})=\sum_{k=1}^{K} \pi_{k} p\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}\right) \tag{9.47}
が得られることを示せ．
乗法定理から
\begin{aligned}
p(\mathbf{x},\mathbf{z},\mid \boldsymbol{\mu},\boldsymbol{\pi})&=p(\mathbf{x} \mid \mathbf{z}, \boldsymbol{\mu})p(\mathbf{z} \mid \boldsymbol{\pi})
\end{aligned}
である．これを\mathbf{z}について周辺化し，(9.52), (9.53)を代入すると
\begin{aligned}
p(\mathbf{x}\mid \boldsymbol{\mu},\boldsymbol{\pi})&=\sum_\mathbf{z} p(\mathbf{x} \mid \mathbf{z}, \boldsymbol{\mu})p(\mathbf{z} \mid \boldsymbol{\pi})\\
&=\sum_{\mathbf{z}}\prod_{k=1}^{K} \pi_{k}^{z_{k}}\prod_{k=1}^{K} p\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}\right)^{z_{k}}\\
&=\sum_{\mathbf{z}}\prod_{k=1}^{K} \pi_{k}^{z_{k}}p\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}\right)^{z_{k}}\\
&=\sum_{k=1}^K \pi_{k}p\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}\right)
\end{aligned}
が得られる．
ここでz_kはone-hot-vectorであるため0か1の値を取る．したがってz_k=0のとき\pi_{k}^{0}p\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}\right)^{0}=1であり，積を考える場合にはz_k=1となる場合に注目すれば良いことがわかる．以上により最後の行の式に変形できる．
演習 9.15
混合ベルヌーイ分布についての期待完全データ対数尤度関数
\mathbb{E}_{\mathbf{Z}}[\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \boldsymbol{\pi})]=\sum_{n=1}^{N} \sum_{k=1}^{K} \gamma\left(z_{n k}\right)\Biggl\{
\ln \pi_{k} \\ \left.+\sum_{i=1}^{D}\left[x_{n i} \ln \mu_{k i}+\left(1-x_{n i}\right) \ln \left(1-\mu_{k i}\right)\right]\right\} \tag{9.55}
を\boldsymbol{\mu}_kについて最大化すると，Mステップ方程式
\boldsymbol{\mu}_{k}=\overline{\mathbf{x}}_{k} \tag{9.59}
が得られることを示せ．
問題文の通り、\mathbb{E}_{\mathbf{Z}}[\ln p]を\boldsymbol{\mu}_kについて最大化する。このために\mu_{ki}での微分を考える。
\begin{aligned} \frac{\partial}{\partial \mu_{k i}} \mathbb{E}_{\mathbf{Z}}[\ln p]
&=\sum_{n=1}^{N} \gamma\left(z_{n k}\right)\left(\frac{x_{n i}}{\mu_{k i}}-\frac{1-x_{n i}}{1-\mu_{k i}}\right) \\
&=\frac{\sum_{n=1}^{N} \gamma\left(z_{n k}\right) x_{n i}-\sum_{n=1}^{N} \gamma\left(z_{n k}\right) \mu_{k i}}{\mu_{k i}\left(1-\mu_{k i}\right)} \end{aligned}
この式を0とおいて、
\begin{aligned}& \sum_{n=1}^{N} \gamma\left(z_{n k}\right) \mu_{k i}=\sum_{n=1}^{N} \gamma\left(z_{nk}\right) x_{n i} \\
& \mu_{k i}=\frac{\sum_{n=1}^{N} \gamma\left(z_{n k}\right) x_{n i}}{\sum_{n=1}^{N} \gamma\left(z_{n k}\right)}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right) x_{n i}\hspace{1em}(\because (9.57))\end{aligned}
を得る。これを踏まえると、
\begin{aligned} \frac{\partial \mathbb{E}_{\mathbf{Z}}[\ln p]}{\partial \boldsymbol{\mu}_{k}} &=\left(\frac{\partial \mathbb{E}_{\mathbf{Z}}[\ln p]}{\partial \mu_{k 1}}, \frac{\partial \mathbb{E}_{\mathbf{Z}}[\ln p]}{\partial \mu_{k 2}}, \cdots\right)^{\mathrm{T}} \\
&=\frac{1}{N_{k}}\left(\sum_{n=1}^{N} \gamma \left(z_{n k}\right) x_{n 1}, \sum_{n=1}^{N} \gamma\left(z_{n k}\right) x_{n 2}, \cdots\right)^{\mathrm{T}} \\
&=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right) \mathbf{x}_{n} \\
&=\overline{\mathbf{x}}_{k}
\end{aligned}
したがって、Mステップでは\boldsymbol{\mu}_{k}=\overline{\mathbf{x}}_{k}となる。
演習 9.16
変数混合ベルヌーイ分布についての期待完全データ対数尤度関数
\mathbb{E}_{\mathbf{Z}}[\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \boldsymbol{\pi})]=\sum_{n=1}^{N} \sum_{k=1}^{K} \gamma\left(z_{n k}\right)\left\{\ln \pi_{k}\right. \\ \left.+\sum_{i=1}^{D}\left[x_{n i} \ln \mu_{k i}+\left(1-x_{n i}\right) \ln \left(1-\mu_{k i}\right)\right]\right\} \tag{9.55}
を，ラグランジュ未定乗数法を用いて混合係数\pi_kの総和を一定に保ちつつ\pi_kについて最大化すると，Mステップ方程式
\pi_{k} = \frac{N_k}{N} \tag{9.60}
が得られることを示せ．
ラグランジュの未定乗数法を用いると
L=\mathbb{E}_{\mathbf{Z}}[\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\mu}, \boldsymbol{\pi})]+\lambda\left(\sum_{k=1}^{K} \pi_{k}-1\right)
これを\pi_{k}で微分してゼロとおくと
\frac{\partial L}{\partial \pi_{k}}=\sum_{n=1}^{N} \frac{\gamma\left(z_{n k}\right)}{\pi_{k}}+\lambda=0-①
この式の両辺に\pi_{k}をかけると
\sum_{n=1}^{N} \gamma\left(z_{n k}\right)+\lambda \pi_{k}=0
kについて総和を取れば
\sum_{n=1}^{N} \sum_{k=1}^{K} \gamma\left(z_{n k}\right)+\lambda=0 \Leftrightarrow \lambda=-\sum_{n=1}^{N} \sum_{k=1}^{K} \gamma\left(z_{n k}\right)
これを①に代入すれば
\pi_{k}=\frac{\sum_{n} \gamma\left(z_{n k}\right)}{\sum_{n, k} \gamma\left(z_{n k}\right)}=\frac{N_{k}}{N}
※ (9.57)を使用
演習 9.17
混合ベルヌーイ分布については離散変数\mathbf{x}_nについて0 \leqslant p\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}\right) \leqslant 1という制約があるために，不完全データ対数尤度関数は上に有界であり，よって尤度関数が発散する特異点が存在しないことを示せ．
不完全データ対数尤度関数\ln p\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}, \boldsymbol{\pi}\right)は、(9.51)で与えられる。
0\leqslant p\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}\right) \leqslant 1より、
\begin{aligned} \ln p\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}, \boldsymbol{\pi}\right) &=\sum_{n=1}^{N} \ln \left(\sum_{k=1}^{K} \pi_{k} p\left(\mathbf{x}_{n} \mid \boldsymbol{\mu}_{k}\right)\right) ...(9.51)\\ & \leqslant \sum_{n=1}^{N} \ln \left(\sum_{k=1}^{K} \pi_{k}\right) \\ &=\sum_{n=1}^{N} \ln 1=0 \end{aligned}
したがって、対数尤度関数は上に有界であり、発散することはない。つまり、そのような特異点は存在しない。
※ 連続変数である混合ガウス分布のときは発散していた（pp.149-150）との対比
演習 9.18
9.3.3節で論じたように，混合ベルヌーイモデルに事前分布を仮定する．ただし各パラメータベクトル\boldsymbol{\mu}_kの分布ρ(\boldsymbol{\mu}_k\mid a_k, b_k)において，\boldsymbol{\mu}_kの各成分はそれぞれ独立にパラメータa_kとb_kを持つベータ分布
\operatorname{Beta}(\mu \mid a, b)=\frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} \mu^{a-1}(1-\mu)^{b-1} \tag{2.13}
に従うものとし，また，p(\boldsymbol{\pi}\mid\boldsymbol{\alpha})はディリクレ分布
\operatorname{Dir}(\mu \mid \alpha)=\frac{\Gamma\left(\alpha_{0}\right)}{\Gamma\left(\alpha_{1}\right) \cdots \Gamma\left(\alpha_{K}\right)} \prod_{k=1}^{K} \mu_{k}^{\alpha_{k}-1} \tag{2.38}
であるとする．このとき，事後分布p(\boldsymbol{\mu},\boldsymbol{\pi}\mid \mathbf{X})を最大化するEMアルゴリズムを導け．
演習9.4より、対数尤度関数ではなく対数事後分布について最大化する場合は、
\begin{aligned}
\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}) = \mathcal{Q}_l(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}) + \ln p(\mathbf{\theta})
\end{aligned}
を最大化すれば良い。ただし、ここでは、
\begin{aligned}
p(\theta) = p(\mathbf{\mu}|\mathbf{a, b})p(\mathbf{\pi}|\mathbf{\alpha})
\end{aligned}
であり、\mathcal{Q}_l(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }})は(9.55)式で表される。
今、(2.13)より、
\begin{aligned}
\ln p(\mathbf{\mu}|\mathbf{a, b}) &= \sum_{i=1}^{D} \ln p(\mathbf{\mu}_k|a_i, b_i) \\
&= \sum_{i=1}^{D} \ln \prod_k \mu_{ki}^{a_k-1}(1-\mu_{ki})^{b_k-1} + \mathrm{const.} \\
&= \sum_{i=1}^{D} \sum_{k=1}^{K} \{(a_k -1) \ln \mu_{ki}+(b_k-1)\ln (1-\mu_{ki})\} + \mathrm{const.}
\end{aligned}
そして、同様に、(2.14)より、
\begin{aligned}
\ln p(\mathbf{\pi}|\mathbf{\alpha})
&= \sum_{k=1}^{K} (\alpha_k -1)\ln \pi_k + \mathrm{const.}
\end{aligned}
したがって、
\begin{aligned}
\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}) = E_{\mathbf{z}}[\ln p] + \sum_{i=1}^{D} \sum_{k=1}^{K} \{(a_k -1) \ln \mu_{ki}+(b_k-1)\ln (1-\mu_{ki})\} + \sum_{k=1}^{K} (\alpha_k -1)\ln \pi_k \tag{9.Ex18.1}
\end{aligned}
を最大化する問題を考えれば良い。よって、\mathcal{Q}を\mu_{ki}について微分を、\pi_kについて未定乗数法を用いて更新式を導出する。
\begin{aligned}
\frac{\partial \mathcal{Q}}{\partial \mu_{ki}} &= \frac{\partial E_{z}[\ln p] }{\partial \mu_{ki}} + \frac{a_i -1}{\mu_{ki}} - \frac{b_i -1 }{1 - \mu_{ki}} \\
&= \sum_n \gamma(z_{nk})\left(\frac{x_{ni}}{\mu_{ki}} -\frac{1-x_{ni}}{1-\mu_{ki}}\right) +\frac{a_i -1}{\mu_{ki}} - \frac{b_i -1 }{1 - \mu_{ki}} &(\because Ex. 9.15)\\
&= \frac{N_k \bar{x}_{ki}+a_i -1 }{\mu_{ki}} - \frac{N_k - N_k \bar{x}_{ki}+b_i -1}{1 - \mu_{ki}} &(\because (9.57), (9.58))
\end{aligned}
よって、整理して
\begin{aligned}
\mu_{ki} = \frac{N_k \bar{x}_{ki}+a_i -1 }{N_k +a_i -1 +b_i -1}
\end{aligned}
次に、\pi_kについて、ラグランジュ未定乗数法を用いて\pi_kを求める。以下の式を考える。
\begin{aligned}
L \propto E_z + \sum_{k=1}^{K} (\alpha_k -1)\ln \pi_k +\lambda \left(\sum_{k=1}^{K} \pi_k -1\right)
\end{aligned}
これを\pi_kについて微分して、
\begin{aligned}
\frac{\partial L}{\partial \pi_k} = \sum_{n=1}^{N} \frac{\gamma(z_{nk})}{\pi_k} + \frac{a_k-1}{\pi_k}+ \lambda = 0 &&(\because Ex. 9.16)
\end{aligned}
そして、両辺に\pi_kを乗じ、kについて和をとると、
\begin{aligned}
\sum_{k=1}^{K} \sum_{n=1}^{N} \gamma(z_{nk}) + \sum_{k=1}^{K}(\alpha_k -1)+ \lambda = 0 &&(\because \sum_k \pi_k = 1)
\end{aligned}
これを整理すると、
\begin{aligned}
\lambda = -N -\alpha_0 + K &&(\because \alpha_0 = \sum_k \alpha_k)
\end{aligned}
これを微分式に代入して整理すると
\begin{aligned}
\pi_k = \frac{\sum_n \gamma(z_{nk})+\alpha_k -1}{-\lambda} = \frac{N_k + \alpha_k -1}{N + \alpha_0 - K} &&(\because (9.57))
\end{aligned}
演習 9.19
\mathbf{x}を，各成分がx_{i j}(i=1, \ldots, D, j=1, \ldots, M)であるD\times M次元のベクトル値の確率変数とする．ただし，x_{ij}の各々は\{0,1\}に値を取り，すべてのIに対し制約条件\sum_{j}x_{ij}=1を満たすものとする．\mathbf{x}の分布は，以下のように，2.2節で議論した(N=1の)多項分布の混合であるとする．
p(\mathbf{x})=\sum_{k=1}^{K} \pi_{k} p\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}\right) \tag{9.84}
ただし
p\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}\right)=\prod_{i=1}^{D} \prod_{j=1}^{M} \mu_{k i j}^{x_{i j}} \tag{9.85}
である．パラメータ\mu_{kij}は，確率p(x_{ij} =1 \mid \boldsymbol{\mu}_k)を表し，0 \leqslant \mu_{kij} \leqslant 1と，各k, iの値について制約条件\sum_{j}\mu_{kij} = 1を満たさなくてはならない．観測データ集合\left\{\mathbf{x}_{n}\right\}(n=1, \ldots, N)が与えられた下で，この分布の混合係数\pi_kと要素のパラメータ\mu_{kij}を最尤推定するEMアルゴリズムの，EステップとMステップの更新式を求めよ．
観測変数\mathbf{x}に対応する潜在変数\mathbf{z}を導入する．このとき観測データセットの条件付き分布は
\begin{aligned}
p(\mathbf{X} \mid\mathbf{Z},\boldsymbol{\mu})&=\prod_{n=1}^N\prod_{k=1}^Kp(\mathbf{x}_n\mid\boldsymbol{\mu}_k)^{z_{nk}}\\
&=\prod_{n=1}^N\prod_{k=1}^K\prod_{i=1}^{D} \prod_{j=1}^{M} \mu_{k i j}^{x_{i j}z_{nk}}
\end{aligned}
と表せ，同様に潜在変数の分布は
p(\mathbf{Z}\mid\pi)=\prod_{n=1}^N\prod_{k=1}^K\pi^{z_{nk}}_k
以上から完全データ尤度関数は
\begin{aligned}
p(\mathbf{X},\mathbf{Z}\mid\boldsymbol{\mu},\pi)&=p(\mathbf{X} \mid\mathbf{Z},\boldsymbol{\mu})p(\mathbf{Z}\mid\pi)\\
&=\prod_{n=1}^N\prod_{k=1}^K\prod_{i=1}^{D} \prod_{j=1}^{M} \mu_{k i j}^{x_{i j}z_{nk}}\prod_{n=1}^N\prod_{k=1}^K\pi^{z_{nk}}_k\\
&=\prod_{n=1}^N\prod_{k=1}^K\prod_{i=1}^{D} \prod_{j=1}^{M} (\pi_k\mu_{k i j}^{x_{i j}})^{z_{nk}}
\end{aligned}
対数をとって添字に注意しながら整理すると
\ln p(\mathbf{X},\mathbf{Z}\mid\boldsymbol{\mu},\pi)=\sum_{n=1}^{N} \sum_{k=1}^{K} z_{n k}\left\{\ln \pi_{k}+\sum_{i=1}^{D} \sum_{j=1}^{M}\mathbf{x}_{nij}\ln\mu_{kij}\right\}
以上のような完全データ対数尤度関数が得られる．
E-stepでは以下の式でp(\mathbf{Z}\mid\pi)によるz_{nk}の期待値計算を行う．
\begin{aligned}
\mathbb{E}[z_{nk}]=\frac{\pi_kp\left(\mathbf{x} \mid \boldsymbol{\mu}_{k}\right)}{\sum_j\pi_jp\left(\mathbf{x} \mid \boldsymbol{\mu}_{j}\right)}=\gamma(z_{nk})
\end{aligned}
M-stepではラグランジュ未定乗数法を用いて各種パラメータの最大化を行う.ラグランジアンが制約条件と係数\lambda_1,\lambda_2を用いて以下のようにかけて
\mathcal{L}=\sum_{n=1}^{N} \sum_{k=1}^{K} \gamma(z_{n k})\left\{\ln \pi_{k}+\sum_{i=1}^{D} \sum_{j=1}^{M}\mathbf{x}_{nij}\ln\mu_{kij}\right\}+\lambda_1(\sum_{k=1}^K\pi_k-1)+\sum_{k=1}^K\sum_{i=1}^D\lambda_{2ki}(\sum_{j=1}^M\mu_{kij}-1)
これを\pi_k,\mu_{kij}について微分しものを0とおいてそれぞれ最大化を行う．
まず\pi_kについて
\begin{aligned}
\frac{\partial\mathcal{L}}{\partial\pi_k}&=\sum_{n=1}^N\frac{\gamma(z_{nk})}{\pi_k}+\lambda_1=0
\end{aligned}
\sum_n\gamma(z_{nk})=N_k,\sum_k\pi_k=1を用いて\lambda_1についてこれを解くと
\frac{N_k}{\pi_k}+\lambda_1=0
よって
が得られる．次に\mu_{kij}について偏微分を0として
\begin{aligned}
\frac{\partial\mathcal{L}}{\partial\mu_{kij}}&=\sum_{n=1}^N\frac{\gamma(z_{nk})\mathbf{x}_{nij}}{\mu_{kij}}+\lambda_{2ki}=0
\end{aligned}
\sum_n\gamma(z_{nk})=N_k,\sum_{j}\mu_{kij}=1を用いて\lambda_{2ki}について解くと
\sum_{n=1}^N\gamma(z_{nk})\mathbf{x}_{nij}+\lambda_{2ki}\mu_{kij}=0
\sum_{n=1}^N\sum_{j=1}^M\gamma(z_{nk})\mathbf{x}_{nij}+\lambda_{2ki}\sum_{j=1}^M\mu_{kij}=0
\begin{aligned}
\lambda_{2ki}=-\sum_{n=1}^N\gamma(z_{nk})\sum_{j=1}^M\mathbf{x}_{nij}=-N_k
\end{aligned}
これを\lambda_{2ki}に代入して式を整理すると
\mu_{kij}=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_{nk})\mathbf{x}_{nij}
となる．
演習 9.20
ベイズ線形回帰モデルについて，期待完全データ対数尤度関数
\begin{aligned} \mathbb{E}[\ln p(\mathbf{t}, \mathbf{w} \mid \alpha, \beta)]=& \frac{M}{2} \ln \left(\frac{\alpha}{2 \pi}\right)-\frac{\alpha}{2} \mathbb{E}\left[\mathbf{w}^{\mathrm{T}} \mathbf{w}\right]+\frac{N}{2} \ln \left(\frac{\beta}{2 \pi}\right) \\ &-\frac{\beta}{2} \sum_{n=1}^{N} \mathbb{E}\left[\left(t_{n}-\mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right)^{2}\right] \end{aligned} \tag{9.62}
の最大化は\alphaに関するMステップの更新式
\alpha=\frac{M}{\mathbb{E}\left[\mathbf{w}^{\mathrm{T}} \mathbf{w}\right]}=\frac{M}{\mathbf{m}_{N}^{\mathrm{T}} \mathbf{m}_{N}+\operatorname{Tr}\left(\mathbf{S}_{N}\right)} \tag{9.63}
を導くことを示せ．
\mathbf{w}の事後分布は(3.49)より以下になる。
p(\mathbf{w} \mid \mathbf{t})=\mathcal{N}\left(\mathbf{w} \mid \mathbf{m}_{N}, \mathbf{S}_{N}\right)(9.62)
\frac{\partial \mathbb{E}}{\partial \alpha}=\frac{M}{2} \frac{1}{\alpha}-\frac{1}{2} \mathbb{E}\left[\mathbf{w}^{\mathrm{T}} \mathbf{w}\right]=0 \\ \therefore \alpha=\frac{M}{\mathbb{E}\left[\mathbf{w}^{\mathrm{T}} \mathbf{w}\right]}
\begin{aligned} \mathbb{E}\left[\mathbf{w}^{\mathrm{T}} \mathbf{w}\right] &=\mathbb{E}\left[\mathbf{w}_{1}^{2}+\mathbf{w}_{2}^{2}+\cdots \cdot\right] \\ &=\mathbb{E}\left[\mathbf{w}_{1}^{2}\right]+\mathbb{E}\left[\mathbf{w}_{2}^{2}\right]+\cdots \cdot \\ &=\operatorname{Tr}\left(\mathbb{E}\left[\mathbf{w} \mathbf{w}^{\mathrm{T}}\right]\right) \quad \because \operatorname{Tr}\left(\mathbb{E}\left[\mathbf{w} \mathbf{w}^{\mathrm{T}}\right]\right)=E\left[\mathbf{w}_{1}^{2}\right]+\cdots \cdot\\ &=\operatorname{Tr}\left(\mathbf{m}_{N} \mathbf{m}_{N}^{\mathrm{T}}+\mathbf{S}_{N}\right)\hspace{1em} \because(2.62) E\left[\mathbf{w} \mathbf{w}^{\mathrm{T}}\right]=\mathbf{m}_{N} \mathbf{m}_{N}^{\mathrm{T}}+\mathbf{S}_{N} \\ &=\mathbf{m}_{N}^{\mathrm{T}} \mathbf{m}_{N}+\operatorname{Tr}\left(\mathbf{S}_{N}\right) \end{aligned}
ゆえに(9.63)が導かれる。
\alpha=\frac{M}{\mathbb{E}\left[\mathbf{\mathbf{w}}^{\mathrm{T}} \mathbf{\mathbf{w}}\right]}=\frac{M}{\mathbf{m}_{N}^{\mathrm{T}} \mathbf{m}_{N}+\operatorname{Tr}\left(\mathbf{S}_{N}\right)}
演習 9.21
ベイズ線形回帰モデルについて，3.5節におけるエビデンスの枠組みを用いて，パラメータ\alphaに関する
\alpha=\frac{M}{\mathbb{E}\left[\mathbf{w}^{\mathrm{T}} \mathbf{w}\right]}=\frac{M}{\mathbf{m}_{N}^{\mathrm{T}} \mathbf{m}_{N}+\operatorname{Tr}\left(\mathbf{S}_{N}\right)} \tag{9.63}
と同様の，パラメータ\betaに関するMステップ更新式を導け．
※教科書P.164より、この目的は第3章のエビデンス近似のEMアルゴリズムを用いた\alpha, \betaの計算法の再定義。
Eステップではパラメータ\alpha, \betaの現在の値に基づく\mathbf{w}の事後分布の計算を行う。(3.49)ですでに\mathbf{w}の事後分布を求めてあるので
p(\mathbf{w} \mid t)=\mathcal{N}\left(\mathbf{w} \mid \mathbf{m}_{N}, \mathbf{S}_{N}\right) \tag{3.49}
Mステップではこの量を\alpha, \betaについて最大化する。完全データ対数尤度
\ln p(\mathbf{t}, \mathbf{w} \mid \alpha, \beta)=\ln p(\mathbf{t} \mid \mathbf{w}, \beta)+\ln p(\mathbf{w} \mid \alpha) \tag{5.61}
\mathbf{w}について期待値を取ると、
\begin{aligned} \mathbb{E}[\ln p(\mathbf{t}, \mathbf{w} \mid \alpha, \beta)]=& \frac{M}{2} \ln \left(\frac{\alpha}{2 \pi}\right)-\frac{\alpha}{2} \mathbb{E}\left[\mathbf{w}^{\mathrm{T}} \mathbf{w}\right]+\frac{N}{2} \ln \left(\frac{\beta}{2 \pi}\right) \\ &-\frac{\beta}{2} \sum_{n=1}^{N} \mathbb{E}\left[\left(t_{n}-\mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right)^{2}\right] \end{aligned} \tag{9.62}
これを\betaで微分して0とおくと
\frac{\partial \mathbb{E}}{\partial \beta}=\frac{N}{2} \frac{1}{\beta}-\frac{1}{2} \sum_{n} \mathbb{E}\left[\left(t_{n}-\mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right)^{2}\right]=0
\begin{aligned}
\mathbb{E}\left[\left(t_{n}-\mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right)^{2}\right]
&=\mathbb{E}\left[t_{n}^{2}-2 t_{n} \mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}+\boldsymbol{\phi}_{n}^{\mathrm{T}} \mathbf{w} \mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right] \\
&=\mathbb{E}[t_{n}]^{2}-2 t_{n} \mathbb{E}[\mathbf{w}]^{\mathrm{T}} \boldsymbol{\phi}_{n}+\boldsymbol{\phi}_{n}^{\mathrm{T}} \mathbb{E}\left[\mathbf{w}\mathbf{w}^{\mathrm{T}}\right] \boldsymbol{\phi}_{n} \\
&=t_{n}^{2}-2 t_{n} \mathbf{m}_{N}^{\mathrm T} \boldsymbol{\phi}_{n}+\boldsymbol{\phi}_{n}^{\mathrm T}\left(\mathbf{m}_{N} \mathbf{m}_{N}^{\mathrm{T}}+\mathbf{S}_{N}\right) \boldsymbol{\phi}_{n}\qquad \because (2.59),(2.62)\\
&=\left(t_{n}-\mathbf{m}_{N}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right)^{2}+\boldsymbol{\phi}_{n}^{\mathrm{T}} \mathbf{S}_{N} \boldsymbol{\phi}_{n}
\end{aligned}
なので以下となる。
\frac{N}{\beta}-\sum_{n=1}^N\left\{\left(t_{n}-\mathbf{m}_{N}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right)^{2}+\boldsymbol{\phi}_{n}^{\mathrm{T}} \mathbf{S}_{N} \boldsymbol{\phi}_{n}\right\}=0
\begin{aligned}
\frac{1}{\beta} &=\frac{1}{N} \sum_{n=1}^N\left\{\left(t_{n}-\mathbf{m}_{N}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right)^{2}+\boldsymbol{\phi}_{n}^{\mathrm{T}} \mathbf{S}_{N} \boldsymbol{\phi}_{n}\right\} \\
&=\frac{1}{N}\left\{\sum_{n=1}^N\left(t_{n}-\mathbf{m}_{N}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right)^{2}+\operatorname{Tr}\left(\mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi} \mathbf{S}_{N}\right)\right\}
\end{aligned}
※最後の変形は以下のようになる。N \times Mの計画行列\mathbf{\Phi}\ (3.16)を利用して
\mathbf{\Phi} = \begin{pmatrix} \boldsymbol{\phi}_1^{\mathrm T} \\ \vdots \\ \boldsymbol{\phi}_N^{\mathrm T} \end{pmatrix} \\
\mathbf{\Phi S}_N\mathbf{\Phi}^{\mathrm T} = \begin{pmatrix} \boldsymbol{\phi}_1^{\mathrm T} \\ \vdots \\ \boldsymbol{\phi}_N^{\mathrm T} \end{pmatrix} \mathbf{S}_N \begin{pmatrix} \boldsymbol{\phi}_1 & \ldots \ \boldsymbol{\phi}_N \end{pmatrix} =
\begin{pmatrix} \boldsymbol{\phi}_1^{\mathrm T}\mathbf{S}_N\boldsymbol{\phi}_1 & \cdots & \boldsymbol{\phi}_1^{\mathrm T}\mathbf{S}_N\boldsymbol{\phi}_N \\ \vdots & \ddots & \vdots \\ \boldsymbol{\phi}_N^{\mathrm T}\mathbf{S}_N\boldsymbol{\phi}_1 & \cdots & \boldsymbol{\phi}_N^{\mathrm T}\mathbf{S}_N\boldsymbol{\phi}_N \end{pmatrix}
この結果から、\sum_{n=1}^{N}\boldsymbol{\phi}_n^{\mathrm T}\mathbf{S}_N\boldsymbol{\phi}_n = \operatorname{Tr}\left(\mathbf{\Phi} \mathbf{S}_{N} \mathbf{\Phi}^{\mathrm{T}}\right) = \operatorname{Tr}\left(\mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi} \mathbf{S}_{N}\right)となることがわかる。
演習 9.22
期待完全データ対数尤度
\mathbb{E}_{\mathbf{w}}[\{\ln p(\mathbf{t} \mid \mathbf{X}, \mathbf{w}, \beta) p(\mathbf{w} \mid \boldsymbol{\alpha})\}] \tag{9.66}
を最大化することで.回帰問題のためのRVM(関連ベクトルマシン)の超パラメータについてのMステップ更新式
\alpha_{i}^{\text {new }} =\frac{1}{m_{i}^{2}+\Sigma_{i i}} \tag{9.67}
\left(\beta^{\text {new }}\right)^{-1} =\frac{\|\mathbf{t}-\Phi \mathbf{m}\|^{2}+\beta^{-1} \sum_{i} \gamma_{i}}{N} \tag{9.68}
を導け．
p(t \mid \mathbf{x}, \mathbf{w}, \beta)=\mathcal{N}\left(t \mid y(\mathbf{x}), \beta^{-1}\right) \tag{7.76}
p(\mathbf{t} \mid \mathbf{X}, \mathbf{w}, \beta)=\prod_{n=1}^{N} p\left(t_{n} \mid \mathbf{x}_{n}, \mathbf{w}, \beta\right) \tag{7.79}
p(\mathbf{w} \mid \boldsymbol{\alpha})=\prod_{i=1}^{M} \mathcal{N}\left(w_{i} \mid 0, \alpha_{i}^{-1}\right) \tag{7.80}
よりまずは対数尤度関数を計算する。
\begin{aligned}
\ln p(\mathbf{t} \mid \mathbf{X}, \mathbf{w}, \beta) p(\mathbf{w} \mid \boldsymbol{\alpha}) &=\ln p(\mathbf{t} \mid \mathbf{X}, \mathbf{w}, \beta)+\ln p(\mathbf{w} \mid \boldsymbol{\alpha}) \\
&=\sum_{n=1}^{N} \ln p\left(t_{n} \mid x_{n}, \mathbf{w}, \beta^{-1}\right)+\sum_{i=1}^{M} \ln \mathcal{N}\left(w_{i} \mid 0, \alpha_{i}^{-1}\right) \\
&=\sum_{n=1}^{N} \ln \mathcal{N}\left(t_{n} \mid \mathbf{w}^{\mathrm T} \boldsymbol{\phi}_{n}, \beta^{-1}\right)+\sum_{i=1}^{M} \ln \mathcal{N}\left(w_{i} \mid 0, \alpha_{i}^{-1}\right) \\
&=\frac{N}{2} \ln \frac{\beta}{2 \pi}-\frac{\beta}{2} \sum_{n=1}^{N}\left(t_{n}-\mathbf{w}^{\mathrm T} \boldsymbol{\phi}_{n}\right)^{2}+\frac{1}{2} \sum_{i=1}^{M} \ln \frac{\alpha_{i}}{2 \pi}-\sum_{i=1}^{M} \frac{\alpha_{i}}{2} w_{i}^{2}
\end{aligned}
Mステップでは、この対数尤度関数の潜在変数\mathbf{w}についての期待値を取る。
期待値を取る\mathbf{w}の事後分布は、
p(\mathbf{w} \mid \mathbf{t}, \mathbf{X}, \boldsymbol{\alpha}, \beta)=\mathcal{N}(\mathbf{w} \mid \mathbf{m}, \mathbf{\Sigma})\tag{7.81}
を利用する。
\mathbb{E}_{\mathbf{w}}[\ln p]=\frac{N}{2} \ln \frac{\beta}{2 \pi}-\frac{\beta}{2} \sum_{n=1}^{N} \mathbb{E}_{\mathbf{w}}\left[\left(t_{n}-\mathbf{w}^{\mathrm T} \boldsymbol{\phi}_{n}\right)^{2}\right]+\frac{1}{2} \sum_{i=1}^{M} \ln \frac{\alpha_{i}}{2 \pi}-\sum_{i=1}^{M} \frac{\alpha_{i}}{2} \mathbb{E}_{\mathbf{w}}\left[w_{i}^{2}\right]
\alpha_{i}に関して微分したものを0とおくと
\frac{1}{2} \cdot \frac{2 \pi}{\alpha_{i}} \cdot \frac{1}{2 \pi}-\frac{1}{2} \mathbb{E}_{\mathbf{w}}\left[w_{i}^{2}\right] =0
\begin{aligned}
\alpha_{i}^{\textrm{new}} &=\frac{1}{\mathbb{E}_{\mathbf{w}}\left[w_{i}^{2}\right]} \\
&=\frac{1}{\mathbb{E}_{\mathbf{w}}\left[\mathbf{ww}^{\mathrm T}\right]_{ii} }\\
\end{aligned}
ここで、\mathbb{E}_{\mathbf{w}}\left[\mathbf{ww}^{\mathrm T}\right]_{ii}は\mathbb{E}_{\mathbf{w}}\left[\mathbf{ww}^{\mathrm T}\right]の行列ii成分を表している。今、(7.81)式から\mathbf{w}は\mathcal{N}(\mathbf{w}\mid \mathbf{m}, \mathbf{\Sigma})から生成されているので（つまり\mathbf{w} \sim \mathcal{N}(\mathbf{w}\mid \mathbf{m}, \mathbf{\Sigma})）、(2.62)の結果から
\mathbb{E}\left[\mathbf{w} \mathbf{w}^{\mathrm T}\right]=\mathbf{m} \mathbf{m}^{\mathrm T}+\mathbf{\Sigma}
これを用いて
\alpha_{i}=\frac{1}{\mathbb{E}_{\mathbf{w}}\left[\mathbf{w} \mathbf{w}^{\mathrm T}\right]_{ii}}=\frac{1}{\left(\mathbf{m m}^{\mathrm T}+\mathbf{\Sigma}\right)_{ii}}=\frac{1}{m_{i}^{2}+\Sigma_{i i}}
同様に\betaを微分したものを0とおくと
\frac{N}{2} \cdot \frac{2 \pi}{\beta} \cdot \frac{1}{2 \pi}-\frac{1}{2} \sum_{n=1}^{N} \mathbb{E}_{\mathbf{w}}\left[\left(t_{n}-\mathbf{w}^{\mathrm T} \boldsymbol{\phi}_{n}\right)^{2}\right]=0
\begin{aligned}\beta^{\textrm{new}} &= \frac{N}{\sum_{n=1}^{N} \mathbb{E}_{\mathbf{w}}\left[\left(t_{n}-\mathbf{w}^{\mathrm T} \boldsymbol{\phi}_{n}\right)^{2}\right]} \\
\left(\beta^{\textrm{new}}\right)^{-1} &=\frac{1}{N}\sum_{n=1}^{N} \mathbb{E}_{\mathbf{w}}\left[\left(t_{n}-\mathbf{w}^{\mathrm T} \boldsymbol{\phi}_{n}\right)^{2}\right] \\
&=\frac{1}{N} \sum_{n=1}^{N}\left\{\left(t_{n}-\mathbf{m}^{\mathrm T} \boldsymbol{\phi}_{N}\right)^{2}+\operatorname{Tr}\left[\boldsymbol{\phi}_{n}^{\mathrm{T}} \mathbf{\Sigma} \boldsymbol{\phi}_{n}\right]\right\} \\
&=\frac{1}{N}\left\{\|\mathbf{t}-\mathbf{\Phi} \mathbf{m}\|^{2}+\operatorname{Tr}\left[\mathbf{\Phi}^{\mathrm T} \mathbf{\Phi} \mathbf{\Sigma}\right]\right\} \end{aligned}
ここで\sum_{n=1}^{N} \mathbb{E}_{\mathbf{w}}\left[\left(t_{n}-\mathbf{w}^{\mathrm T} \boldsymbol{\phi}_{n}\right)^{2}\right]部分については演習9.21と同様に変形した。
以下では\operatorname{Tr}\left[\mathbf{\Phi}^{\mathrm T} \mathbf{\Phi} \mathbf{\Sigma}\right] = \beta^{-1} \sum_{i} \gamma_{i}を示す。RVMでは\mathbf{\Phi} , \mathbf{\Sigma}に関する情報は
\mathbf{\Sigma}=\left(\mathbf{A}+\beta \mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi}\right)^{-1} \tag{7.83}
で与えられる。この両辺にまず\left(\mathbf{A}+\beta \mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi}\right)をかけて
\begin{aligned} \left(\mathbf{A}+\beta \mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi}\right)\mathbf{\Sigma} &= \mathbf{I} \\
\mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi}\mathbf{\Sigma} &= \beta^{-1}\left(\mathbf{I}-\mathbf{A\Sigma}\right)
\end{aligned}
これを\operatorname{Tr}\left(\mathbf{\Phi} \mathbf{\Sigma} \mathbf{\Phi}^{\mathrm T}\right)=\operatorname{Tr}\left(\mathbf{\Phi}^{\mathrm T} \mathbf{\Phi} \mathbf{\Sigma} \right)に代入すれば
\begin{aligned} \operatorname{Tr}\left(\mathbf{\Phi}^{\mathrm T} \mathbf{\Phi} \mathbf{\Sigma} \right) &=\beta^{-1} \operatorname{Tr}(\mathbf{I}-\mathbf{A\Sigma}) \\
&=\beta^{-1} \sum_{i}\left(1-\alpha_{i} \Sigma_{ii}\right) \\ &=\beta^{-1} \sum_{i} \gamma_{i} \end{aligned}
以上で(9.68)式が導かれた。
演習 9.23
7.2.1節において回帰問題のためのRVMの超パラメータ\boldsymbol{\alpha}と\betaの更新式(7.87)と(7.88)を導くのに周辺尤度を直接最大化した．同様に9.3.4節において，同じ周辺尤度を最大化するのにEMアルゴリズムを用いて更新式
\alpha_{i}^{\textrm{new}} =\frac{1}{m_{i}^{2}+\Sigma_{i i}} \tag{9.67}
\left(\beta^{\textrm{new}}\right)^{-1} =\frac{\|\mathbf{t}-\mathbf{\Phi} \mathbf{m}\|^{2}+\beta^{-1} \sum_{i} \gamma_{i}}{N} \tag{9.68}
を求めた．任意の停留点においてこれら2組の更新式が厳密に等価であることを示せ．
\begin{aligned}
\alpha_{i}^{\textrm{new}} &=\frac{\gamma_{i}}{m_{i}^{2}}\qquad (7.87) \\
\left(\beta^{\mathrm{new}}\right)^{-1} &=\frac{\|\mathbf{t}-\mathbf{\Phi} \mathbf{m}\|^{2}}{N-\sum_{i} \gamma_{i}}\quad (7.88)\\
\end{aligned}
\gamma_{i}=1-\alpha_{i} \Sigma_{i i}\quad (7.89)
この問題では「(7.87)かつ(7.88)」⇔「(9.67)かつ(9.68)」を示せばよい。(7.87)に(7.89)を代入すれば
\begin{aligned}
\alpha_{i}^{\textrm{new}} &=\frac{1-\alpha_{i}^{\textrm{new}} \Sigma_{i i}}{m_{i}^{2}} \\
\Leftrightarrow \alpha_{i}^{\textrm{new}}&=\frac{1}{m_{i}^{2}+\sum_{i i}}\quad (\Leftrightarrow(9.67))
\end{aligned}
βに関しては、逆に導出する。(9.68)において\betaを\beta^{\textrm{new}}に置き換えると
\begin{aligned}
\left(\beta^{\textrm{new}}\right)^{-1}&=\frac{\|\mathbf{t}-\mathbf{\Phi} \mathbf{m}\|^{2}+\left(\beta^{\textrm{new}}\right)^{-1} \sum_{i} \gamma_{i}}{N} \\
\Leftrightarrow\left(\beta^{\textrm{new}}\right)^{-1}&=\frac{\|\mathbf{t}-\mathbf{\Phi} \mathbf{m}\|^{2}}{N-\sum_{i} \gamma_{i}}
\end{aligned}
演習 9.24
\ln p(\mathbf{X} \mid \boldsymbol{\theta})=\mathcal{L}(q, \theta)+\mathrm{KL}(q \| p) \tag{9.70}
を示せ．ただし\mathcal{L}(q, \boldsymbol{\theta})と\mathrm{KL}(q \| p)はそれぞれ
\mathcal{L}(q, \boldsymbol{\theta})=\sum_{\mathbf{Z}} q(\mathbf{Z}) \ln \left\{\frac{p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})}{q(\mathbf{Z})}\right\} \tag{9.71}
\mathrm{KL}(q \| p)=-\sum_{\mathbf{Z}} q(\mathbf{Z}) \ln \left\{\frac{p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})}{q(\mathbf{Z})}\right\}\tag{9.72}
で定義されている．
乗法定理より、\ln p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})=\ln p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})+\ln p(\mathbf{X} \mid \boldsymbol{\theta})であるから
\begin{aligned} \mathcal{L}(q, \boldsymbol{\theta}) &=\sum_{\mathbf{Z}} q(\mathbf{Z})\{\ln p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})+\ln p(\mathbf{X} \mid \boldsymbol{\theta})-\ln q(\mathbf{Z})\} \\ &=\ln p(\mathbf{X} \mid \boldsymbol{\theta}) \sum_{\mathbf{Z}} q(\mathbf{Z})+\sum_{\mathbf{Z}} q(\mathbf{Z}) \ln \frac{p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})}{q(\mathbf{Z})} \\ &=\ln p(\mathbf{X} \mid \boldsymbol{\theta})-\operatorname{KL}(q \| p) \end{aligned}
より、(9.70)が示された。
演習 9.25
対数尤度関数\ln p(\mathbf{X}\mid \boldsymbol{\theta})とその下界
\mathcal{L}(q, \boldsymbol{\theta})=\sum_{\mathbf{Z}} q(\mathbf{Z}) \ln \left\{\frac{p(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})}{q(\mathbf{Z})}\right\} \tag{9.71}
ただしq(\mathbf{Z}) = p(\mathbf{Z}\mid \mathbf{X}, \boldsymbol{\theta}^{\textrm{old}}) は，点\boldsymbol{\theta} = \boldsymbol{\theta}^{\textrm{old}}において同じ勾配を持つことを示せ．
前提条件である、q(\mathbf{Z})=p\left(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{\text {old }}\right) の場合、
\mathrm{KL}(q \| p)=-\sum_{\mathbf{Z}} q(\mathbf{Z}) \ln \left\{\frac{p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta})}{q(\mathbf{Z})}\right\} = 0
となっている。また常に\mathrm{KL}(q \| p) \ge 0なので、点\boldsymbol{\theta} = \boldsymbol{\theta}^{\textrm{old}}において
\left.\frac{\partial}{\partial \boldsymbol{\theta}} \mathrm{KL}(q \| p)\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}_{\textrm{old}}} = \mathbf{0}
となる。
\ln p(\mathbf{X} \mid \boldsymbol{\theta})=\mathcal{L}(q, \theta)+\mathrm{KL}(q \| p) \tag{9.70}
式の\boldsymbol{\theta}についての両辺の微分をとり、点\boldsymbol{\theta} = \boldsymbol{\theta}^{\textrm{old}}において右辺の第2項は\mathbf{0}になるので、左辺の微分（勾配）と右辺の第1項の微分（勾配）は同じになる。すなわち、
\left.\frac{\partial}{\partial \boldsymbol{\theta}} \operatorname{ln}p(\mathbf{X} \mid \boldsymbol{\theta})\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}_{\textrm{old}}} = \left.\frac{\partial}{\partial \boldsymbol{\theta}} \mathcal{L}(q, \boldsymbol{\theta})\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}_{\textrm{old}}}
を得る。
演習 9.26
混合ガウス分布について，負担率を1つのデータ点\mathbf{x}_mのみについてしか更新しない逐次型EMアルゴリズムを考える．Mステップの式(9.17)と(9.18)から始めて，混合要素の平均を更新する式
\boldsymbol{\mu}_{k}^{\text {new }}=\boldsymbol{\mu}_{k}^{\text {old }}+\left(\frac{\gamma^{\text {new }}\left(z_{m k}\right)-\gamma^{\text {old }}\left(z_{m k}\right)}{N_{k}^{\text {new }}}\right)\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\text {old }}\right) \tag{9.78}
N_{k}^{\text {new }}=N_{k}^{\text {old }}+\gamma^{\text {new }}\left(z_{m k}\right)-\gamma^{\text {old }}\left(z_{m k}\right)\tag{9.79}
を導け．
※P.171の流れの通り、Eステップでn=mとなる1つのデータ点mのみについて\gamma(z_{nk})を更新する。混合分布が指数型分布族の場合は、負担率が十分統計量となるので、n=mについての更新差分だけ考慮すれば良い。
N_{k}=\sum_{n=1}^{N} \gamma\left(z_{n k}\right) \tag{9.18}
より
\begin{aligned}
N_{k}^{\textrm{new}} &=\gamma^{\textrm{old}}\left(z_{1 k}\right)+\gamma^{\textrm{old}}\left(z_{2 k}\right)+\cdots+\gamma^{\textrm{new}}\left(z_{m k}\right)+\cdots+\gamma^{\textrm{old}}\left(z_{N k}\right) \\
&=\sum_{n=1}^{N} \gamma^{\textrm{old}}\left(z_{n k}\right)-\gamma^{\textrm{old}}\left(z_{m k}\right)+\gamma^{\textrm{new}}\left(z_{m k}\right) \\
&=N_{k}^{\textrm{old}}+\gamma^{\textrm{new}}\left(z_{m k}\right)-\gamma^{\textrm{old}}\left(z_{m k}\right)
\end{aligned}
以上から(9.79)が示された。
\mu_{k}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right) \mathbf{x}_{n} \tag{9.17}
より
\begin{aligned}
\boldsymbol{\mu}_{k}^{\textrm{new}}&=\frac{1}{N_{k}^{\textrm{new}}}\left\{\gamma^{\textrm{old}}\left(z_{1 k}\right) \mathbf{x}_{1}+\gamma^{\textrm{old}}\left(z_{2 k}\right) \mathbf{x}_{2}+\cdots+\gamma^{\textrm{new}}\left(z_{m k}\right) \mathbf{x}_{m}+\cdots+\gamma^{\textrm{old}}\left(z_{N k}\right) \mathbf{x}_{N}\right\}\\
&=\frac{1}{N_{k}^{\textrm{new}}}\left\{\sum_{n=1}^{N} \gamma^{\textrm{old}}\left(z_{n k}\right) \mathbf{x}_{n}-\gamma^{\textrm{old}}\left(z_{m k}\right) \mathbf{x}_{m}+\gamma^{\textrm{new}}\left(z_{m k}\right) \mathbf{x}_{m}\right\}\\
&=\frac{1}{N_{k}^{\textrm{new}}}\left\{N_{k}^{\textrm{old}} \boldsymbol{\mu}_{k}^{\textrm{old}}-\gamma^{\textrm{old}}\left(z_{m k}\right) \mathbf{x}_{m}+\gamma^{\textrm{new}}\left(z_{m k}\right) \mathbf{x}_{m}\right\}\\
&=\frac{1}{N_{k}^{\textrm{new}}}\left[\underbrace{\left\{N_{k}^{\textrm{new}}+\gamma^{\textrm{old}}\left(z_{m k}\right)-\gamma^{\textrm{new}}\left(z_{m k}\right)\right.}_{(\because(9.79))}\} \boldsymbol{\mu}_{k}^{\textrm{old}}-\gamma^{\textrm{old}}\left(z_{m k}\right) \mathbf{x}_{m}+\gamma^{\textrm{new}}\left(z_{m k}\right) \mathbf{x}_{m}\right] \\
&=\frac{1}{N_{k}^{\textrm{new}}}\left[N_{k}^{\textrm{new}} \boldsymbol{\mu}_{k}^{\textrm{old}}-\left\{\gamma^{\textrm{new}}\left(z_{m k}\right)-\gamma^{\textrm{old}}\left(z_{m k}\right)\right\} \boldsymbol{\mu}_{k}^{\textrm{old}}+\left\{\gamma^{\textrm{new}}\left(z_{m k}\right)-\gamma^{\textrm{old}}\left(z_{m k}\right)\right\} \mathbf{x}_{m}\right]\\
&=\boldsymbol{\mu}_{k}^{\textrm{old}}+\frac{1}{N_{k}^{\textrm{new}}}\left\{\gamma^{\textrm{new}}\left(z_{m k}\right)-\gamma^{\textrm{old}}\left(z_{m k}\right)\right\}\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\textrm{old}}\right)
\end{aligned}
以上から(9.78)が示された。
演習 9.27
平均の更新式
\boldsymbol{\mu}_{k}^{\text {new }}=\boldsymbol{\mu}_{k}^{\text {old }}+\left(\frac{\gamma^{\text {new }}\left(z_{m k}\right)-\gamma^{\text {old }}\left(z_{m k}\right)}{N_{k}^{\text {new }}}\right)\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\text {old }}\right) \tag{9.78}
と同様に，負担率を逐次型で更新する場合の混合ガウス分布の共分散行列と混合係数を更新するMステップの更新式を求めよ．
教科書P.154の混合ガウス分布のためのEMアルゴリズムから、(9.24)-(9.26)を利用して、Mステップでの更新後の混合係数は
\begin{aligned}
\pi^{\textrm{new}} &=\frac{N^{\textrm{new}}}{N}(\because(9.26)) \\ &=\frac{N_{k}^{\text {old}}+\gamma^{\text {new}}\left(z_{m k}\right)-\gamma^{\text {old}}\left(z_{m k}\right)}{N} \\
&=\pi^{\textrm{old}} + \frac{\gamma^{\text {new}}\left(z_{m k}\right)-\gamma^{\text {old}}\left(z_{m k}\right)}{N}
\end{aligned}
となる。
次に更新後の共分散行列は、古い値は(9.19)を利用して
\mathbf{\Sigma}_{k}^{\textrm{old}}=\frac{1}{N_{k}^{\textrm{old}}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}^{\textrm{old}}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}^{\textrm{old}}\right)^{\mathrm{T}}
である。これをn=mについて更新すると
\begin{aligned}\mathbf{\Sigma}_{k}^{\textrm{new}}&=\frac{1}{N_{k}^{\textrm{new}}}\left(\sum_{n \neq m} \gamma^{\textrm{old}}\left(z_{n k}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}^{\textrm{old}}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}^{\textrm{old}}\right)^{\mathrm{T}} +\gamma^{\textrm{new}}\left(z_{m k}\right)\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\textrm{new}}\right)\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\textrm{new}}\right)^{\mathrm{T}}\right) \\
&=\frac{1}{N_{k}^{\textrm{new}}}\left(N_{k}^{\textrm{old}} \mathbf{\Sigma}_{k}^{\textrm{old}}-\gamma^{\textrm{old}}\left(z_{m k}\right)\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\textrm{old}}\right)\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\textrm{old}}\right)^{\mathrm{T}}+\gamma^{\textrm{new}}\left(z_{m k}\right)\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\textrm{new}}\right)\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\textrm{new}}\right)^{\mathrm{T}}\right) \\
&=\mathbf{\Sigma}_{k}^{\mathrm{old}} -\frac{\gamma^{\textrm{old}}\left(z_{m k}\right)}{N_{k}^{\mathrm{new}}}\left(\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\textrm{old}}\right)\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\textrm{old}}\right)^{\mathrm{T}}-\mathbf{\Sigma}_{k}^{\textrm{old}}\right) +\frac{\gamma^{\textrm{new}}\left(z_{m k}\right)}{N_{k}^{\textrm{new}}}\left(\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\textrm{new}}\right)\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{k}^{\textrm{new}}\right)^{\mathrm{T}}-\mathbf{\Sigma}_{k}^{\textrm{old}}\right)
\end{aligned}
となる。途中で(9.79)式の変形N_{k}^{\textrm{old}} = N_{k}^{\textrm{new}} - \gamma^{\text {new }}\left(z_{m k}\right)+ \gamma^{\text {old }}\left(z_{m k}\right)を利用して\mathbf{\Sigma}_{k}^{\mathrm{old}}とその更新差分に分解した。
Discussion
演習9.4の解答部分の5行目についてですが、Σの下の記号はkではなくzではないかと思いますが、いかがでしょうか？
長い間返信できておらずすみません。ご指摘ありがとうございました。修正いたしました。
演習9.13の解答部分の下から11行目についてですが、μハットのところは、P(Xn|μハット）ではないかと思いますがいかがでしょうか？
私の思い違いであれば、スルーしてくださいませ。
ありがとうございます。間違えておりましたので、修正いたしました。
演習9.18の解答部分の下から11行目についてですが、分母のNk・xkiはNkではないかと思います。
また、７行目以降のΣzは、iについての和だと思いますので、zではない文字（D等）の方が良いように思います。
修正しました。