混合ガウス分布の最尤推定に一般のEMアルゴリズムを適用

はじめに

混合ガウス分布の最尤推定の記事で、混合ガウス分布の説明はしております。
使用する変数などはそちらをご覧ください。
また、本記事は混合ガウス分布の最尤推定に一般のEMアルゴリズムを適用する内容となっておりますので、
一般のEMアルゴリズムの記事をご覧になっていない方は、まずはそちらをご覧ください。

一応、本記事にも混合ガウス分布のグラフィカルモデルを貼り付けておきます。
f:id:olj611:20210915055429p:plain:w360

最尤推定

各データ点が混合ガウス分布から独立に生成されると仮定すると、尤度は以下のようになります。

$\begin{eqnarray} p({\bf X}|{\boldsymbol\pi},{\boldsymbol\mu},{\bf\Sigma})=\prod_{n=1}^N\sum_{k=1}^K\pi_k\mathcal{N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)\tag{1} \end{eqnarray}$

対数を取ります。

$\begin{eqnarray} \ln p({\bf X}|{\boldsymbol\pi},{\boldsymbol\mu},{\bf\Sigma})=\sum_{n=1}^N\ln\left(\sum_{k=1}^K\pi_k\mathcal{N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)\right)\tag{2} \end{eqnarray}$

式 $(2)$ はlog-sumになっており、解析的解くことができません。
そこでEMアルゴリズムを適用します。
一般のEMアルゴリズムと表記を合わせるために ${\boldsymbol\theta}=\{{\boldsymbol\pi},{\boldsymbol\mu},{\bf\Sigma}\}$ とおきます。

1.初期化

${\boldsymbol\pi^{old},{\boldsymbol\mu}^{old},{\bf\Sigma}^{old}}$ を初期化します。

$\begin{eqnarray} &&{\boldsymbol\theta^{old}}=\{{\boldsymbol\pi^{old},{\boldsymbol\mu}^{old},{\bf\Sigma}^{old}}\}\tag{3}\\ \end{eqnarray}$

2.Eステップ

Eステップでは、 ${\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})$ を計算します。
${\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})$ の計算には、 $p({\bf Z}|{\bf X},{\boldsymbol\theta}^{old})$ と $\ln p({\bf X},{\bf Z}|{\boldsymbol\theta})$ が必要です。

まず、 $\ln p({\bf X},{\bf Z}|{\boldsymbol\theta})$ の計算します。

以下の計算で使用する式を書いておきます。

$\begin{eqnarray} &&p({\bf z}_n|{\boldsymbol\theta})=p({\bf z}_n|{\boldsymbol\pi})=\prod_{k=1}^K\pi_k^{z_{nk}}\tag{4}\\ &&p({\bf x}_n|{\bf z}_n,{\boldsymbol\theta})=p({\bf x}_n|{\bf z}_n,{\boldsymbol\mu},{\bf\Sigma})=\prod_{k=1}^K{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)^{z_{nk}}\tag{5}\\ \end{eqnarray}$

式 $(4),(5)$ より以下が成り立ちます。

$\begin{eqnarray} p({\bf x}_n,{\bf z}_n|{\boldsymbol\theta})&=&p({\bf z}_n|{\boldsymbol\theta})p({\bf x}_n|{\bf z}_n,{\boldsymbol\theta})\\ &=&p({\bf z}_n|{\boldsymbol\pi})p({\bf x}_n|{\bf z}_n,{\boldsymbol\mu},{\bf\Sigma})\\ &=&\left(\prod_{k=1}^K\pi_k^{z_{nk}}\right)\left(\prod_{k=1}^K{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)^{z_{nk}}\right)\\ &=&\prod_{k=1}^K\pi_k^{z_{nk}}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)^{z_{nk}}\tag{6}\\ \end{eqnarray}$

完全データの尤度関数 $p({\bf X},{\bf Z}|{\boldsymbol\theta})$ を計算します。

$\begin{eqnarray} p({\bf X},{\bf Z}|{\boldsymbol\theta})&=&\prod_{n=1}^N p({\bf x}_n,{\bf z}_n|{\boldsymbol\theta})\\ &=&\prod_{n=1}^N \prod_{k=1}^K\pi_k^{z_{nk}}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)^{z_{nk}}\tag{7}\\ \end{eqnarray}$

完全データの対数尤度関数 $\ln p({\bf X},{\bf Z}|{\boldsymbol\theta})$ を計算します。

$\begin{eqnarray} \ln p({\bf X},{\bf Z}|{\boldsymbol\theta})&=&\ln\left(\prod_{n=1}^N \prod_{k=1}^K\pi_k^{z_{nk}}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)^{z_{nk}}\right)\\ &=&\sum_{n=1}^N \sum_{k=1}^K\ln\left(\pi_k^{z_{nk}}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)^{z_{nk}}\right)\\ &=&\sum_{n=1}^N \sum_{k=1}^Kz_{nk}\left(\ln\pi_k+\ln{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)\right)\tag{8}\\ \end{eqnarray}$

次に、 $p({\bf Z}|{\bf X},{\boldsymbol\theta})$ を計算します。
※本来は、 $p({\bf Z}|{\bf X},{\boldsymbol\theta}^{old})$ を計算する必要がありますが、 ${\boldsymbol\theta}$ か ${\boldsymbol\theta}^{old}$ の違いだけなので、 $p({\bf Z}|{\bf X},{\boldsymbol\theta})$ を計算すれば十分です。

$\begin{eqnarray} p({\bf Z}|{\bf X},{\boldsymbol\theta})&=&\frac{p({\bf X}, {\bf Z}|{\boldsymbol\theta})}{p({\bf X}|{\boldsymbol\theta})}\\ &=&\frac{p({\bf X}, {\bf Z}|{\boldsymbol\theta})}{\sum_{\bf Z}p({\bf X},{\bf Z}|{\boldsymbol\theta})}\\ &=&\frac{\prod_{n=1}^N\prod_{k=1}^K\pi_k^{z_{nk}}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k{\bf\Sigma}_k)^{z_{nk}}}{\sum_{\bf Z}\prod_{n=1}^N\prod_{k=1}^K\pi_k^{z_{nk}}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)^{z_{nk}}}\\ &\propto&\prod_{n=1}^N\prod_{k=1}^K\pi_k^{z_{nk}}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)^{z_{nk}}\tag{9}\\ \end{eqnarray}$

式 $(9)$ は $n$ についての積の形をしているので、事後分布の下では $\{{\bf z}_n\}$ は独立であることが分かります。

$\begin{eqnarray} p({\bf Z}|{\bf X},{\boldsymbol\theta})=\prod_{n=1}^Np({\bf z}_n|{\bf x}_n,{\boldsymbol\theta})\tag{10}\\ \end{eqnarray}$

式 $(10)$ の因数 $p({\bf z}_n|{\bf x}_n,{\boldsymbol\theta})$ を求めます。

$\begin{eqnarray} p({\bf z}_n|{\bf x}_n,{\boldsymbol\theta})&=&\frac{p({\bf x}_n,{\bf z}_n|{\boldsymbol\theta})}{p({\bf x}_n|{\bf z}_n,{\boldsymbol\theta})}\\ &=&\frac{p({\bf x}_n,{\bf z}_n|{\boldsymbol\theta})}{\sum_{{\bf z}_n}p({\bf x}_n,{\bf z}_n|{\boldsymbol\theta})}\\ &=&\frac{\prod_{k=1}^K\pi_k^{z_{nk}}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k{\bf\Sigma}_k)^{z_{nk}}}{\sum_{{\bf z}_n}\prod_{k=1}^K\pi_k^{z_{nk}}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)^{z_{nk}}}\tag{11}\\ \end{eqnarray}$

${\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})$ の計算時に $\langle z_{nk}\rangle_{p({\bf Z}|{\bf X},{\boldsymbol\theta}^{old})}$ が出てくるので、先に求めておきます。

$\begin{eqnarray} \langle z_{nk}\rangle_{p({\bf Z}|{\bf X},{\boldsymbol\theta}^{old})}&=&\langle z_{nk}\rangle_{p({\bf z}_n|{\bf x}_n,{\boldsymbol\theta}^{old})p({\bf Z}_{\backslash n}|{\bf X}_{\backslash n},{\boldsymbol\theta}^{old})}\\ &=&\left\langle\langle z_{nk}\rangle_{p({\bf Z}_{\backslash n}|{\bf X}_{\backslash n},{\boldsymbol\theta}^{old})}\right\rangle_{p({\bf z}_n|{\bf x}_n,{\boldsymbol\theta}^{old})}\\ &=&\langle z_{nk}\rangle_{p({\bf z}_n|{\bf x}_n,{\boldsymbol\theta}^{old})}\\ &=&\sum_{{\bf z}_n}z_{nk}\frac{\prod_{k'=1}^K(\pi_{k'}^{old})^{z_{nk'}}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_{k'}^{old},{\bf\Sigma}_{k'}^{old})^{z_{n{k'}}}} {\sum_{{\bf z}_n}\prod_{j=1}^K(\pi_j^{old})^{z_{nj}}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_j^{old},{\bf\Sigma}_j^{old})^{z_{nj}}}\\ &=&\frac{\pi_k^{old}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k^{old},{\bf\Sigma}_k^{old})} {\sum_{j=1}^K\pi_j^{old}{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_j^{old},{\bf\Sigma}_j^{old})}\\ &=&\gamma({z_{nk}})\tag{12}\\ \end{eqnarray}$

式 $(12)$ の4行目から5行目にかけての変形ですが
左側の $\sum_{{\bf z}_n}z_{nk}$ ですが、 $z_{nk}=1$ のときのみ1になります。
分子の $z_{nk'}=1$ になるのは、 $k'=k$ の時なので、分子は $\pi_k^{old},{\boldsymbol\mu}_k^{old},{\bf\Sigma}_k^{old}$ が残ります。
分母についてはこちらの記事の式 $(5)$ を参考にしてください。

${\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})$ を求めます。

$\begin{eqnarray} {\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})&=&\sum_{\bf Z}p({\bf Z}|{\bf X},{\boldsymbol\theta}^{old})\ln p({\bf X},{\bf Z}|{\boldsymbol\theta})\\ &=&\langle\ln p({\bf X},{\bf Z}|{\boldsymbol\theta})\rangle_{p({\bf Z}|{\bf X},{\boldsymbol\theta}^{old})}\\ &=&\left\langle\sum_{n=1}^N \sum_{k=1}^Kz_{nk}\left(\ln\pi_k+\ln{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)\right)\right\rangle_{p({\bf Z}|{\bf X},{\boldsymbol\theta}^{old})}\\ &=&\sum_{n=1}^N \sum_{k=1}^K\langle z_{nk}\rangle_{p({\bf Z}|{\bf X},{\boldsymbol\theta}^{old})}\left(\ln\pi_k+\ln{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)\right)\\ &=&\sum_{n=1}^N \sum_{k=1}^K\gamma(z_{nk})\left(\ln\pi_k+\ln{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)\right)\tag{13}\\ \end{eqnarray}$

Eステップでの目的は ${\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})$ を計算することであり、
${\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})$ の中で ${\boldsymbol\theta}^{old}$ に依存するのは、 $\gamma(z_{nk})=\langle z_{nk}\rangle_{p({\bf Z}|{\bf X},{\boldsymbol\theta}^{old})}$ なので、
Eステップで実質行うのは、 $\gamma(z_{nk})$ の計算です。

また、式 $(13)$ ですが、
「本当は完全データの対数尤度関数 $\ln p({\bf X},{\bf Z}|{\boldsymbol\theta})$ を使って計算したいが、 $z_{nk}$ は潜在変数なので観測できない。
だったら、 $z_{nk}$ の代わりに ${\bf z}$ の事後分布による期待値 $\gamma(z_{nk})=\langle z_{nk}\rangle_{p({\bf Z}|{\bf X},{\boldsymbol\theta}^{old})}$ で置き換えてしまおう。」
というノリです。

3. Mステップ

以下のMステップの全ての更新式で使う為、 $N_k$ を以下のようにおきます。

$\begin{eqnarray} N_k=\sum_{n=1}^N\gamma(z_{nk})\tag{14} \end{eqnarray}$

${\boldsymbol\mu}_k$ の更新式

準備として、 $\dfrac{\partial}{\partial{\boldsymbol\mu}_k}\ln\mathcal{N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)$ を求めます。

$\begin{eqnarray} \frac{\partial}{\partial{\boldsymbol\mu}_k}\ln\mathcal{N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)&=&\frac{\partial}{\partial{\boldsymbol\mu}_k}\left(-\frac{1}{2}({\bf x}_n-{\boldsymbol\mu}_k)^\top{\bf\Sigma}_k^{-1}({\bf x}_n-{\boldsymbol\mu}_k)\right)\\ &=&-\frac{1}{2}\frac{\partial}{\partial{\boldsymbol\mu}_k}\left({\bf x}_n^\top{\bf\Sigma}_k^{-1}{\bf x}_n-{\bf x}_n^\top{\bf\Sigma}_k^{-1}{\boldsymbol\mu}_k-{\boldsymbol\mu}_k^\top{\bf\Sigma}_k^{-1}{\bf x}_n+{\boldsymbol\mu}_k^\top{\bf\Sigma}_k^{-1}{\boldsymbol\mu}_k\right)\\ &=&-\frac{1}{2}\left(\frac{\partial}{\partial{\boldsymbol\mu}_k}{\bf x}_n^\top{\bf\Sigma}_k^{-1}{\bf x}_n-\frac{\partial}{\partial{\boldsymbol\mu}_k}{\bf x}_n^\top{\bf\Sigma}_k^{-1}{\boldsymbol\mu}_k-\frac{\partial}{\partial{\boldsymbol\mu}_k}{\boldsymbol\mu}_k^\top{\bf\Sigma}_k^{-1}{\bf x}_n+\frac{\partial}{\partial{\boldsymbol\mu}_k}{\boldsymbol\mu}_k^\top{\bf\Sigma}_k^{-1}{\boldsymbol\mu}_k\right)\\ &=&-\frac{1}{2}\left(-({\bf x}_n^\top{\bf\Sigma}_k^{-1})^\top-{\bf\Sigma}_k^{-1}{\bf x}_n+2{\bf\Sigma}_k^{-1}{\boldsymbol\mu}_k\right)\\ &=&-\frac{1}{2}\left(-{\bf\Sigma}_k^{-1}{\bf x}_n-{\bf\Sigma}_k^{-1}{\bf x}_n+2{\bf\Sigma}_k^{-1}{\boldsymbol\mu}_k\right)\\ &=&-\frac{1}{2}\left(2{\bf\Sigma}_k^{-1}({\boldsymbol\mu}_k-{\bf x}_n)\right)\\ &=&{\bf\Sigma}_k^{-1}({\bf x}_n-{\boldsymbol\mu}_k)\tag{15} \end{eqnarray}$

${\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})$ を ${\boldsymbol\mu}_k$ で微分して $={\bf 0}$ とおきます。

$\begin{eqnarray} &&\frac{\partial}{\partial{\boldsymbol\mu}_k}{\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})={\bf 0}\\ &&\Leftrightarrow\frac{\partial}{\partial{\boldsymbol\mu}_k}\sum_{n=1}^N\sum_{j=1}^K\gamma(z_{nj})(\ln\pi_j+\ln{\mathcal{N}({\bf x}_n|{\boldsymbol\mu}_j},{\bf\Sigma}_j)={\bf 0}\\ &&\Leftrightarrow\sum_{n=1}^N\sum_{j=1}^K\gamma(z_{nj})\frac{\partial}{\partial{\boldsymbol\mu}_k}(\ln\pi_j+\ln{\mathcal{N}({\bf x}_n|{\boldsymbol\mu}_j},{\bf\Sigma}_j)={\bf 0}\\ &&\Leftrightarrow\sum_{n=1}^N\gamma(z_{nj})\frac{\partial}{\partial{\boldsymbol\mu}_k}{\mathcal{N}({\bf x}_n|{\boldsymbol\mu}_k},{\bf\Sigma}_k)={\bf 0}\\ &&\Leftrightarrow {\bf\Sigma}_k^{-1}\sum_{n=1}^N\gamma(z_{nk})({\bf x}_n-{\boldsymbol\mu}_k)={\bf 0}\\ &&\Leftrightarrow \sum_{n=1}^N\gamma(z_{nk})({\bf x}_n-{\boldsymbol\mu}_k)={\bf 0}\\ &&\Leftrightarrow {\boldsymbol\mu}_k=\frac{\sum_{n=1}^N\gamma(z_{nk}){\bf x}_n}{\sum_{n=1}^N\gamma(z_{nk})}\\ &&\Leftrightarrow {\boldsymbol\mu}_k=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_{nk}){\bf x}_n\tag{16}\\ \end{eqnarray}$

${\bf\Sigma}_k$ の更新式

ここで使用する行列の微分やトレースの公式を書き出します。

$\begin{eqnarray} &&{\bf x}^\top{\bf A}{\bf x}={\rm Tr}({\bf A}{\bf x}{\bf x}^\top)\tag{17}\\ &&\frac{\partial}{\partial{\bf A}}\ln|{\bf A}|=({\bf A}^{-1})^\top\tag{18}\\ &&\frac{\partial}{\partial{\bf A}}{\rm Tr}({\bf A}^{-1}{\bf B})=-({\bf A}^{-1}{\bf B}{\bf A}^{-1})^\top\tag{19}\\ \end{eqnarray}$

準備として、 $\dfrac{\partial}{\partial{\bf\Sigma}_k}\ln{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)$ を求めます。

$\begin{eqnarray} \ln{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)&=&-\frac{D}{2}\ln2\pi-\frac{1}{2}\ln|{\bf\Sigma}_k|-\frac{1}{2}({\bf x}_n-{\boldsymbol\mu}_k)^\top{\bf\Sigma}_k^{-1}({\bf x}_n-{\boldsymbol\mu}_k)\\ &=&-\frac{D}{2}\ln2\pi-\frac{1}{2}\ln|{\bf\Sigma}_k|-\frac{1}{2}\mathrm{Tr}\left({\bf\Sigma}_k^{-1}({\bf x}_n-{\boldsymbol\mu}_k)({\bf x}_n-{\boldsymbol\mu}_k)^\top\right)\tag{20}\\ \end{eqnarray}$

式 $(20)$ の式変形に式 $(17)$ を用いました。
式 $(20)$ に対数を取ります。

$\begin{eqnarray} \frac{\partial}{\partial{\bf\Sigma}_k}\ln{\mathcal N}({\bf x}_n|{\boldsymbol\mu}_k,{\bf\Sigma}_k)&=&-\frac{1}{2}\frac{\partial}{\partial{\bf\Sigma}_k}\ln|{\bf\Sigma}_k|-\frac{1}{2}\frac{\partial}{\partial{\bf\Sigma}_k}\mathrm{Tr}\left({\bf\Sigma}_k^{-1}({\bf x}_n-{\boldsymbol\mu}_k)({\bf x}_n-{\boldsymbol\mu}_k)^\top\right)\\ &=&-\frac{1}{2}({\bf\Sigma}_k^{-1})^\top+\frac{1}{2}\left({\bf\Sigma}_k^{-1}({\bf x}_n-{\boldsymbol\mu}_k)({\bf x}_n-{\boldsymbol\mu}_k)^\top{\bf\Sigma}_k^{-1}\right)^\top\\ &=&-\frac{1}{2}{\bf\Sigma}_k^{-1}+\frac{1}{2}{\bf\Sigma}_k^{-1}({\bf x}_n-{\boldsymbol\mu}_k)({\bf x}_n-{\boldsymbol\mu}_k)^\top{\bf\Sigma}_k^{-1}\tag{21}\\ \end{eqnarray}$

式 $(21)$ の式変形には式 $(18),(19)$ を用いました。
${\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})$ を ${\bf\Sigma}_k$ で微分して $={\bf O}$ とおきます。

$\begin{eqnarray} &&\frac{\partial}{\partial{\bf\Sigma}_k}{\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})={\bf O}\\ &&\Leftrightarrow\frac{\partial}{\partial{\bf\Sigma}_k}\sum_{n=1}^N\sum_{j=1}^K\gamma(z_{nj})(\ln\pi_j+\ln{\mathcal{N}({\bf x}_n|{\boldsymbol\mu}_j},{\bf\Sigma}_j))={\bf O}\\ &&\Leftrightarrow\sum_{n=1}^N\sum_{j=1}^K\gamma(z_{nj})\frac{\partial}{\partial{\bf\Sigma}_k}(\ln\pi_j+\ln{\mathcal{N}({\bf x}_n|{\boldsymbol\mu}_j},{\bf\Sigma}_j))={\bf O}\\ &&\Leftrightarrow\sum_{n=1}^N\gamma(z_{nj})\frac{\partial}{\partial{{\bf\Sigma}_k}}{\mathcal{N}({\bf x}_n|{\boldsymbol\mu}_k},{\bf\Sigma}_k)={\bf O}\\ &&\Leftrightarrow \sum_{n=1}^N\gamma(z_{nk})\left(-\frac{1}{2}{\bf\Sigma}_k^{-1}+\frac{1}{2}{\bf\Sigma}_k^{-1}({\bf x}_n-{\boldsymbol\mu}_k)({\bf x}_n-{\boldsymbol\mu}_k)^\top{\bf\Sigma}_k^{-1}\right)={\bf O}\\ &&\Leftrightarrow {\bf\Sigma}_k^{-1}\sum_{n=1}^N\gamma(z_{nk})\left({\bf I}-({\bf x}_n-{\boldsymbol\mu}_k)({\bf x}_n-{\boldsymbol\mu}_k)^\top{\bf\Sigma}_k^{-1}\right)={\bf O}\\ &&\Leftrightarrow \sum_{n=1}^N\gamma(z_{nk})\left({\bf I}-({\bf x}_n-{\boldsymbol\mu}_k)({\bf x}_n-{\boldsymbol\mu}_k)^\top{\bf\Sigma}_k^{-1}\right)={\bf O}\\ &&\Leftrightarrow \sum_{n=1}^N\gamma(z_{nk})=\sum_{n=1}^N\gamma(z_{nk})({\bf x}_n-{\boldsymbol\mu}_k)({\bf x}_n-{\boldsymbol\mu}_k)^\top{\bf\Sigma}_k^{-1}\\ &&\Leftrightarrow {\bf\Sigma}_k=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_{nk})({\bf x}_n-{\boldsymbol\mu}_k)({\bf x}_n-{\boldsymbol\mu}_k)^\top\tag{22}\\ \end{eqnarray}$

$\pi_k$ の更新式

$\sum_{k=1}^K\pi_k=1\Leftrightarrow\sum_{k=1}^K\pi_k-1=0$ の制約の下で最大化します。
この時ラグランジュ関数 $G$ は以下のようになります。

$\begin{eqnarray} G&=&{\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})+\lambda\left(\sum_{j=1}^K\pi_j-1\right)\\ &=&\sum_{n=1}^N\sum_{j=1}^K\gamma(z_{nj})(\ln\pi_j+\ln{\mathcal{N}({\bf x}_n|{\boldsymbol\mu}_j},{\bf\Sigma}_j))+\lambda\left(\sum_{j=1}^K\pi_j-1\right)\tag{23}\\ \end{eqnarray}$

$G$ を $\pi_k$ で微分して $=0$ とおきます。

$\begin{eqnarray} &&\frac{\partial}{\partial\pi_k}G=0\\ &&\Leftrightarrow\sum_{n=1}^N\sum_{j=1}^K\gamma(z_{nj})\frac{\partial}{\partial\pi_k}(\ln\pi_j+\ln{\mathcal{N}({\bf x}_n|{\boldsymbol\mu}_j},{\bf\Sigma}_j))+\lambda\frac{\partial}{\partial\pi_k}\left(\sum_{k'=1}^K\pi_{k'}-1\right)=0\\ &&\Leftrightarrow\sum_{n=1}^N\gamma(z_{nj})\frac{1}{\pi_k}+\lambda={\bf 0}\\ &&\Leftrightarrow\frac{1}{\pi_k}\sum_{n=1}^N\gamma(z_{nk})+\lambda=0\\ &&\Leftrightarrow\frac{N_k}{\pi_k}+\lambda=0\\ &&\Leftrightarrow N_k=-\lambda\pi_k\tag{24}\\ \end{eqnarray}$

ところで、

$\begin{eqnarray} &&N=\sum_{k=1}^K N_k\tag{25}\\ \end{eqnarray}$

が成り立ちます。式 $(24)$ を式 $(25)$ に代入します。

$\begin{eqnarray} &&N=\sum_{k=1}^K(-\lambda\pi_k)\\ &&\Leftrightarrow N=-\lambda\sum_{k=1}^K\pi_k\\ &&\Leftrightarrow\lambda=-N\tag{26}\\ \end{eqnarray}$

式 $(26)$ を式 $(24)$ に代入して、 $\pi_k$ について解きます。

$\begin{eqnarray} \pi_k=\frac{N_k}{N}\tag{27}\\ \end{eqnarray}$

最後に

更新式は勿論、混合ガウス分布の最尤推定の記事と同じです。
混合ガウス分布の最尤推定の記事と比較すると、Mステップ以降の式変形が若干簡単になっています。

偉人の名言

f:id:olj611:20210305151651p:plain:w300
前進をしない人は、後退をしているのだ。
ゲーテ

参考文献

パターン認識と機械学習下巻
ベイズ推論による機械学習入門

動画

本記事は動画投稿後、大幅に書き直しました。

機械学習基礎理論独習

誤りがあればご指摘いただけると幸いです。数式が整うまで少し時間かかります。リンクフリーです。