混合ベルヌーイ分布の最尤推定に一般のEMアルゴリズムを適用

混合ベルヌーイ分布とは

${\bf x}$ を $D$ 次元ベクトル ${\bf x}=(x_1,\ldots,x_D)^\top,x_i\in\{0,1\}$ ,
${\boldsymbol\mu}$ を $D$ 次元ベクトル ${\boldsymbol\mu}=(\mu_1,\ldots,\mu_D)^\top,\mu_i\in(0,1)$ とします。
${\bf x}$ は以下の分布に従うと仮定します。

$\begin{eqnarray} p({\bf x}|{\boldsymbol\mu})&=&\prod_{i=1}^D{\rm Bern}(x_i|\mu_i)\\ &=&\prod_{i=1}^D\mu_i^{x_i}(1-\mu_i)^{(1-x_i)}\tag{1} \end{eqnarray}$

式 $(1)$ は積の形をしているので、 ${\boldsymbol\mu}$ が与えられているとき各変数 $x_i$ は独立であることが分かります。

確率分布 $(1)$ の平均を求めます。

$\begin{eqnarray} \langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu})}&=&(\langle x_1\rangle_{p({\bf x}|{\boldsymbol\mu})},\ldots,\langle x_D\rangle_{p({\bf x}|{\boldsymbol\mu})})^\top\\ &=&(\langle x_1\rangle_{p(x_1|\mu_1)},\ldots,\langle x_D\rangle_{p(x_D|\mu_D)})^\top\\ &=&(\mu_1,\ldots,\mu_D)^\top\\ &=&{\boldsymbol\mu}\tag{2} \end{eqnarray}$

確率分布 $(1)$ の共分散を求めるために、先に $x_ix_j$ の期待値を $i\not=j,i=j$ と場合分けして求めます。
まず、 $i\not=j$ の場合です。

$\begin{eqnarray} \langle x_ix_j\rangle_{p({\bf x}|{\boldsymbol\mu})}&=&\langle x_ix_j\rangle_{p(x_i|\mu_i)p(x_j|\mu_j)}\\ &=&\langle x_i\rangle_{p(x_i|\mu_i)}\langle x_j\rangle_{p(x_i|\mu_j)}\\ &=&\mu_i\mu_j\tag{2} \end{eqnarray}$

次に、 $i=j$ の場合です。

$\begin{eqnarray} \langle x_ix_j\rangle_{p({\bf x}|{\boldsymbol\mu})}&=&\langle x_i^2\rangle_{p({\bf x}|{\boldsymbol\mu})}\\ &=&\langle x_i^2\rangle_{p(x_i|\mu_i)}\\ &=&\mu_i\tag{3} \end{eqnarray}$

式 $(2),(3)$ をまとめて書くと、以下のようになります。

$\begin{eqnarray} \langle x_ix_j\rangle_{p({\bf x}|{\boldsymbol\mu})}= \left\{ \begin{array}{l} \mu_i\mu_j\ (i\not=j)\\ \mu_i\ \ \ \ \ (i=j) \end{array}\tag{4} \right. \end{eqnarray}$

確率分布 $(1)$ の共分散を求めます。

$\begin{eqnarray} \langle({\bf x}-{\boldsymbol\mu})({\bf x}-{\boldsymbol\mu})^\top\rangle_{p({\bf x}|{\boldsymbol\mu})}&=&\langle({\bf x}-{\boldsymbol\mu})({\bf x}^\top-{\boldsymbol\mu}^\top)\rangle_{p({\bf x}|{\boldsymbol\mu})}\\ &=&\langle{\bf x}{\bf x}^\top-{\bf x}{\boldsymbol\mu}^\top-{\boldsymbol\mu}{\bf x}^\top+{\boldsymbol\mu}{\boldsymbol\mu}^\top\rangle_{p({\bf x}|{\boldsymbol\mu})}\\ &=&\langle{\bf x}{\bf x}^\top\rangle_{p({\bf x}|{\boldsymbol\mu})}-\langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu})}{\boldsymbol\mu}^\top-{\boldsymbol\mu}\langle{\bf x}^\top\rangle_{p({\bf x}|{\boldsymbol\mu})}+{\boldsymbol\mu}{\boldsymbol\mu}^\top\\ &=&\langle{\bf x}{\bf x}^\top\rangle_{p({\bf x}|{\boldsymbol\mu})}-{\boldsymbol\mu}{\boldsymbol\mu}^\top-{\boldsymbol\mu}{\boldsymbol\mu}^\top+{\boldsymbol\mu}{\boldsymbol\mu}^\top\\ &=&\langle{\bf x}{\bf x}^\top\rangle_{p({\bf x}|{\boldsymbol\mu})}-{\boldsymbol\mu}{\boldsymbol\mu}^\top\\ &=&\begin{pmatrix} \mu_1-\mu_1^2 & \cdots &\mu_1\mu_i-\mu_1\mu_i & \cdots & \mu_1\mu_D-\mu_1\mu_D\\ \vdots & & & & \vdots \\ \mu_i\mu_1-\mu_i\mu_1 & \cdots &\mu_i-\mu_i^2 &\cdots & \mu_i\mu_D-\mu_i\mu_D\\ \vdots & & & & \vdots \\ \mu_D\mu_1-\mu_D\mu_1 & \cdots &\mu_D\mu_i-\mu_D\mu_i &\cdots & \mu_D-\mu_D^2\\ \end{pmatrix}\\ &=&\begin{pmatrix} \mu_1(1-\mu_1) & \cdots & 0 & \cdots & 0\\ \vdots & & & & \vdots \\ 0 & \cdots &\mu_i(1-\mu_i) &\cdots & 0\\ \vdots & & & & \vdots \\ 0 & \cdots & 0 &\cdots & \mu_D(1 -\mu_D)\\ \end{pmatrix}\\ &=&{\rm diag}(\mu_i(1-\mu_i))\tag{5} \end{eqnarray}$

式 $(5)$ の下から3行目の変形は、式 $(4)$ を利用しました。

ここで、確率分布 $(1)$ の混合分布を考えます。

$\begin{eqnarray} p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})=\sum_{k=1}^K\pi_kp({\bf x}|{\boldsymbol\mu}_k)\tag{6} \end{eqnarray}$

ただし、

$\begin{eqnarray} &&{\boldsymbol\mu}=\{{\boldsymbol\mu}_1,\ldots,{\boldsymbol\mu}_K\},\\ &&{\boldsymbol\pi}=\{{\pi_1,\ldots,\pi_K}\},\ 0\leqslant \pi_i\leqslant 1,\ \sum_{k=1}^K\pi_k=1,\\ &&p({\bf x}|{\boldsymbol\mu}_k)=\prod_{i=1}^D\mu_{ki}^{x_i}(1-\mu_{ki})^{(1-x_i)}\tag{7} \end{eqnarray}$

とします。

式 $(6)$ が混合ベルヌーイ分布です。

混合ベルヌーイ分布の平均と共分散

混合ベルヌーイ分布 $(6)$ の平均を求めます。

$\begin{eqnarray} \langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}&=&\sum_{\bf x}{\bf x}p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})\\ &=&\sum_{\bf x}{\bf x}\sum_{k=1}^K\pi_kp({\bf x}|{\boldsymbol\mu}_k)\\ &=&\sum_{k=1}^K\pi_k\sum_{\bf x}{\bf x}p({\bf x}|{\boldsymbol\mu}_k)\\ &=&\sum_{k=1}^K\pi_k\langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu}_k)}\\ &=&\sum_{k=1}^K\pi_k{\boldsymbol\mu}_k\tag{8} \end{eqnarray}$

混合ベルヌーイ分布 $(6)$ の共分散を求めます。

$\begin{eqnarray} &&\langle({\bf x}-\langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})})({\bf x}-\langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})})^\top\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}\\ &=&\langle({\bf x}-\langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})})({\bf x}^\top-\langle{\bf x}\rangle^\top_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})})\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}\\ &=&\langle{\bf x}{\bf x}^\top-{\bf x}\langle{\bf x}\rangle^\top_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}-\langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}{\bf x}^\top+\langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}\langle{\bf x}\rangle^\top_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}\\ &=&\langle{\bf x}{\bf x}^\top\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}-\langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}\langle{\bf x}\rangle^\top_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}-\langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}\langle{\bf x}\rangle^\top_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}+\langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}\langle{\bf x}\rangle^\top_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}\\ &=&\langle{\bf x}{\bf x}^\top\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}-\langle{\bf x}\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}\langle{\bf x}\rangle^\top_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}\tag{9} \end{eqnarray}$

ここで、式 $(9)$ を変形するため、 $\langle{\bf x}{\bf x}^\top\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}$ を計算します。

$\begin{eqnarray} \langle{\bf x}{\bf x}^\top\rangle_{p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})}&=&\sum_{\bf x}{\bf x}{\bf x}^\top p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})\\ &=&\sum_{\bf x}{\bf x}{\bf x}^\top\sum_{k=1}^K\pi_kp({\bf x}|{\boldsymbol\mu}_k)\\ &=&\sum_{k=1}^K\pi_k\sum_{\bf x}{\bf x}{\bf x}^\top p({\bf x}|{\boldsymbol\mu}_k)\\ &=&\sum_{k=1}^K\pi_k\langle{\bf x}{\bf x}^\top\rangle_{p({\bf x}|{\boldsymbol\mu}_k)}\tag{10} \end{eqnarray}$

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

$\begin{eqnarray} \langle({\bf x}-{\boldsymbol\mu}_k)({\bf x}-{\boldsymbol\mu}_k)^\top\rangle_{p({\bf x}|{\boldsymbol\mu}_k)}&=&\langle{\bf x}{\bf x}^\top\rangle_{p({\bf x}|{\boldsymbol\mu}_k)}-{\boldsymbol\mu}_k{\boldsymbol\mu}_k^\top\tag{11} \end{eqnarray}$

ここで、 ${\bf\Sigma}_k=\langle({\bf x}-{\boldsymbol\mu}_k)({\bf x}-{\boldsymbol\mu}_k)^\top\rangle_{p({\bf x}|{\boldsymbol\mu}_k)}$ とおくと、式 $(11)$ より以下のように書けます。

$\begin{eqnarray} \langle{\bf x}{\bf x}^\top\rangle_{p({\bf x}|{\boldsymbol\mu}_k)}={\bf\Sigma}_k+{\boldsymbol\mu}_k{\boldsymbol\mu}_k^\top\tag{12} \end{eqnarray}$

よって、式 $(10),(12)$ より、式 $(9)$ は次のように書けます。

潜在変数を潜り込ませる

データ集合 ${\bf X}=\{{\bf x}_1,\ldots,{\bf x}_N\}$ が与えられたとき、尤度関数は次式で表されます。

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

式 $(14)$ より、対数尤度関数は次のようになります。

$\begin{eqnarray} \ln p({\bf X}|{\boldsymbol\mu},{\boldsymbol\pi})&=&\sum_{n=1}^N\ln\sum_{k=1}^K\pi_kp({\bf x}|{\boldsymbol\mu}_k)\tag{15} \end{eqnarray}$

式 $(15)$ はlog-sumの形をしているので、解析的に解くのが難しいです。
そこで、混合ガウス分布同様、潜在変数を潜り込ませて、EMアルゴリズムを用いることにします。

まず、潜在変数を潜り込ませます。
${\bf x}$ に対応する $K$ 次元の2値確率変数 ${\bf z}=(z_1,\ldots,z_K)^\top$ を導入します。
これは1-ofK符号化を取るとします。 $z_k\in\{0,1\}$ かつ $\sum_{k=1}^Kz_k=1$ を満たします。
${\bf z}$ の事前分布を

$\begin{eqnarray} p({\bf z}|{\boldsymbol\pi})=\prod_{k=1}^K\pi_k^{z_k}\tag{16} \end{eqnarray}$

とします。
$z_k=1$ のとき、 ${\bf x}$ のクラスタが $k$ であることを表すので、以下の式が成り立ちます。

$\begin{eqnarray} p({\bf x}|z_k=1,{\boldsymbol\mu})=p({\bf x}|{\boldsymbol\mu}_k)\tag{17} \end{eqnarray}$

式 $(16),(17)$ より、 ${\bf x}$ の条件付き分布は

$\begin{eqnarray} p({\bf x}|{\bf z},{\boldsymbol\mu})=\prod_{k=1}^Kp({\bf x}|{\boldsymbol\mu}_k)^{z_k}\tag{18} \end{eqnarray}$

となります。

また、式 $(16),(18)$ から次式が成り立ちます。

$\begin{eqnarray} p({\bf x},{\bf z}|{\boldsymbol\mu},{\boldsymbol\pi})&=&p({\bf x}|{\bf z},{\boldsymbol\mu})p({\bf z}|{\boldsymbol\pi})\\ &=&\prod_{k=1}^Kp({\bf x}|{\boldsymbol\mu}_k)^{z_k}\prod_{k=1}^K\pi_k^{z_k}\\ &=&\prod_{k=1}^K\left(\pi_kp({\bf x}|{\boldsymbol\mu}_k)\right)^{z_k}\tag{19} \end{eqnarray}$

式 $(19)$ を ${\bf z}$ について周辺化します。

$\begin{eqnarray} p({\bf x}|{\boldsymbol\mu},{\boldsymbol\pi})&=&\sum_{\bf z}\prod_{k=1}^K\left(\pi_kp({\bf x}|{\boldsymbol\mu}_k)\right)^{z_k}\\ &=&\sum_{\bf z}(\pi_kp({\bf x}|{\boldsymbol\mu}_k))^{z_1}(\pi_kp({\bf x}|{\boldsymbol\mu}_k))^{z_2}\cdots(\pi_kp({\bf x}|{\boldsymbol\mu}_k))^{z_K}\\ &=&(\pi_kp({\bf x}|{\boldsymbol\mu}_1))^1(\pi_kp({\bf x}|{\boldsymbol\mu}_2))^0\cdots(\pi_kp({\bf x}|{\boldsymbol\mu}_K))^0\\ &+&(\pi_kp({\bf x}|{\boldsymbol\mu}_1))^0(\pi_kp({\bf x}|{\boldsymbol\mu}_2))^1\cdots(\pi_kp({\bf x}|{\boldsymbol\mu}_K))^0\\ &+&\cdots\\ &+&(\pi_kp({\bf x}|{\boldsymbol\mu}_1))^0(\pi_kp({\bf x}|{\boldsymbol\mu}_2))^0\cdots(\pi_kp({\bf x}|{\boldsymbol\mu}_K))^1\\ &=&\pi_kp({\bf x}|{\boldsymbol\mu}_1)+\pi_kp({\bf x}|{\boldsymbol\mu}_2)+\cdots+\pi_kp({\bf x}|{\boldsymbol\mu}_K)\\ &=&\sum_{k=1}^K\pi_kp({\bf x}|{\boldsymbol\mu}_k)\tag{20} \end{eqnarray}$

式 $(20)$ は混合ベルヌーイ分布の式 $(6)$ と一致するので潜在変数 ${\bf z}$ を潜り込ませてもよいことが分かります。

データ集合 ${\bf X}$ に対応する潜在変数の集合を ${\bf Z}=\{{\bf z}_1,\ldots,{\bf z}_N\}$ とします。

f:id:olj611:20210915062346p:plain:w280

さて、ここからは混合ガウス分布の最尤推定に一般のEMアルゴリズムを適用の記事同様、EMアルゴリズムを適用していきます。
一般のEMアルゴリズムと表記を合わせるために ${\boldsymbol\theta}=\{{\boldsymbol\mu},{\boldsymbol\pi}\}$ とおきます。

1. 初期化

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

$\begin{eqnarray} &&{\boldsymbol\theta^{old}}=\{{\boldsymbol\mu}^{old},{\boldsymbol\pi^{old}}\}\tag{21}\\ \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})$ の計算します。

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

$\begin{eqnarray} p({\bf X},{\bf Z}|{\boldsymbol\theta})&=&\prod_{n=1}^Np({\bf x}_n,{\bf z}_n|{\boldsymbol\theta})\\ &=&\prod_{n=1}^Np({\bf x}_n,{\bf z}_n|{\boldsymbol\mu},{\boldsymbol\pi})\\ &=&\prod_{n=1}^N\prod_{k=1}^K(\pi_kp({\bf x}_n|{\boldsymbol\mu}_k))^{z_{nk}}\\ &=&\prod_{n=1}^N\prod_{k=1}^K\left(\pi_k\prod_{i=1}^D\mu_{ki}^{x_{ni}}(1-\mu_{ki})^{(1-x_{ni})}\right)^{z_{nk}}\tag{22} \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\left(\pi_k\prod_{i=1}^D\mu_{ki}^{x_{ni}}(1-\mu_{ki})^{(1-x_{ni})}\right)^{z_{nk}}\right)\\ &=&\sum_{n=1}^N\sum_{k=1}^Kz_{nk}\left(\ln\pi_k+\sum_{i=1}^D\left(x_{ni}\ln\mu_{ki}+(1-x_{ni})\ln(1-\mu_{ki})\right)\right)\tag{23} \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\left(\pi_k\prod_{i=1}^D\mu_{ki}^{x_{ni}}(1-\mu_{ki})^{(1-x_{ni})}\right)^{z_{nk}}}{\sum_{\bf Z}\prod_{n=1}^N\prod_{k=1}^K\left(\pi_k\prod_{i=1}^D\mu_{ki}^{x_{ni}}(1-\mu_{ki})^{(1-x_{ni})}\right)^{z_{nk}}}\\ &\propto&\prod_{n=1}^N\prod_{k=1}^K\left(\pi_k\prod_{i=1}^D\mu_{ki}^{x_{ni}}(1-\mu_{ki})^{(1-x_{ni})}\right)^{z_{nk}}\tag{24}\\ \end{eqnarray}$

式 $(24)$ は $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{25}\\ \end{eqnarray}$

式 $(25)$ の因数 $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\left(\pi_k\prod_{i=1}^D\mu_{ki}^{x_{ni}}(1-\mu_{ki})^{(1-x_{ni})}\right)^{z_{nk}}}{\sum_{{\bf z}_n}\prod_{k=1}^K\left(\pi_k\prod_{i=1}^D\mu_{ki}^{x_{ni}}(1-\mu_{ki})^{(1-x_{ni})}\right)^{z_{nk}}}\tag{26} \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\left(\pi_{k'}^{old}\prod_{i=1}^D(\mu_{{k'}i}^{old})^{x_{ni}}(1-\mu_{{k'}i}^{old})^{(1-x_{ni})}\right)^{z_{n{k'}}}} {\sum_{{\bf z}_n}\prod_{j=1}^K\left(\pi_j^{old}\prod_{i=1}^D(\mu_{ji}^{old})^{x_{ni}}(1-\mu_{ji}^{old})^{(1-x_{ni})}\right)^{z_{nj}}}\\ &=&\frac{\pi_k^{old}\prod_{i=1}^D(\mu_{ki}^{old})^{x_{ni}}(1-\mu_{ki}^{old})^{(1-x_{ni})}} {\sum_{j=1}^K\pi_j^{old}\prod_{i=1}^D(\mu_{ji}^{old})^{x_{ni}}(1-\mu_{ji}^{old})^{(1-x_{ni})}}\\ &=&\gamma({z_{nk}})\tag{27} \end{eqnarray}$

式 $(27)$ の4行目から5行目にかけての変形ですがこちらの記事の式 $(12)$ の説明を参照してください。

${\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+\sum_{i=1}^D\left(x_{ni}\ln\mu_{ki}+(1-x_{ni})\ln(1-\mu_{ki})\right)\right)\right\rangle_{p({\bf Z}|{\bf X},{\boldsymbol\theta}^{old})}\\ &=&\sum_{n=1}^N\sum_{k=1}^K\left\langle z_{nk}\right\rangle_{p({\bf Z}|{\bf X},{\boldsymbol\theta}^{old})}\left(\ln\pi_k+\sum_{i=1}^D\left(x_{ni}\ln\mu_{ki}+(1-x_{ni})\ln(1-\mu_{ki})\right)\right)\\ &=&\sum_{n=1}^N\sum_{k=1}^K\gamma(z_{nk})\left(\ln\pi_k+\sum_{i=1}^D\left(x_{ni}\ln\mu_{ki}+(1-x_{ni})\ln(1-\mu_{ki})\right)\right)\tag{28} \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})$ の計算です。

また、式 $(28)$ ですが、
「本当は完全データの対数尤度関数 $\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{29} \end{eqnarray}$

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

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

$\begin{eqnarray} &&\frac{\partial}{\partial\mu_{ki}}{\mathcal{Q}}({\boldsymbol\theta},{\boldsymbol\theta}^{old})=0\\ &&\Leftrightarrow \frac{\partial}{\partial\mu_{ki}}\sum_{n=1}^N\sum_{k'=1}^K\gamma(z_{nk'})\left(\ln\pi_{k'}+\sum_{i'=1}^D\left(x_{ni'}\ln\mu_{k'i'}+(1-x_{ni'})\ln(1-\mu_{k'i'})\right)\right)=0\\ &&\Leftrightarrow \sum_{n=1}^N\gamma(z_{nk})\left(x_{ni}\frac{\partial}{\partial\mu_{ki}}\ln\mu_{ki}+(1-x_{ni})\frac{\partial}{\partial\mu_{ki}}\ln(1-\mu_{ki})\right)=0\\ &&\Leftrightarrow \sum_{n=1}^N\gamma(z_{nk})\left(\frac{x_{ni}}{\mu_{ki}}-\frac{1-x_{ni}}{1-\mu_{ki}}\right)=0\\ &&\Leftrightarrow \sum_{n=1}^N\gamma(z_{nk})\frac{x_{ni}}{\mu_{ki}}-\sum_{n=1}^N\gamma(z_{nk})\frac{1-x_{ni}}{1-\mu_{ki}}=0\\ &&\Leftrightarrow \sum_{n=1}^N\gamma(z_{nk})x_{ni}(1-\mu_{ki})-\sum_{n=1}^N\gamma(z_{nk})(1-x_{ni})\mu_{ki}=0\\ &&\Leftrightarrow \mu_{ki}\sum_{n=1}^N\gamma(z_{nk})=\sum_{n=1}^N\gamma(z_{nk})x_{ni}\\ &&\Leftrightarrow \mu_{ki}=\frac{\sum_{n=1}^N\gamma(z_{nk})x_{ni}}{\sum_{n=1}^N\gamma(z_{nk})}\\ &&\Leftrightarrow \mu_{ki}=\frac{\bar{x}_{ni}}{N_k}\tag{30} \end{eqnarray}$

式 $(30)$ で

$\begin{eqnarray} \bar{x}_{ki}=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_{nk})x_{ni}\tag{31} \end{eqnarray}$

とおきました。
また、 $\bar{\bf x}_k=(\bar{x}_{k1},\ldots,\bar{x}_{kD})^\top$ とおくことにより、次の更新式が得られます。

$\begin{eqnarray} {\boldsymbol\mu}_k=\bar{\bf x}_k\tag{32} \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_{k=1}^K\gamma(z_{nk})\left(\ln\pi_{k}+\sum_{i=1}^D\left(x_{ni}\ln\mu_{ki}+(1-x_{ni})\ln(1-\mu_{ki})\right)\right)+\lambda\left(\sum_{j=1}^K\pi_j-1\right)\tag{33}\\ \end{eqnarray}$

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

$\begin{eqnarray} &&\frac{\partial}{\partial\pi_k}G=0\\ &&\Leftrightarrow\sum_{n=1}^N\sum_{k'=1}^K\gamma(z_{nk'})\frac{\partial}{\partial\pi_k}\left(\ln\pi_{k'}+\sum_{i=1}^D\left(x_{ni}\ln\mu_{k'i}+(1-x_{ni})\ln(1-\mu_{k'i})\right)\right)+\lambda\frac{\partial}{\partial\pi_k}\left(\sum_{j'=1}^K\pi_{j'}-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{34} \end{eqnarray}$

ところで、

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

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

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

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

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

アルゴリズム

1.初期化
${\boldsymbol\mu}_k,\pi_k$ を初期化します。 $(k=1,\ldots,K)$

2.Eステップ
$\gamma(z_{nk})$ を更新します。 $(n=1,\ldots,N,\ k=1,\ldots,K)$

$\begin{eqnarray} &&\gamma(z_{nk})\leftarrow\frac{\pi_k\prod_{i=1}^D\mu_{ki}^{x_{ni}}(1-\mu_{ki})^{(1-x_{ni})}} {\sum_{j=1}^K\pi_j\prod_{i=1}^D\mu_{ji}^{x_{ni}}(1-\mu_{ji})^{(1-x_{ni})}}\tag{38}\\ \end{eqnarray}$

3.Mステップ
${\boldsymbol\mu}_k,\pi_k$ を更新します。 $(k=1,\ldots,K)$

$\begin{eqnarray} &&{\boldsymbol\mu}_k\leftarrow\bar{\bf x}_k\tag{39}\\ &&\pi_k\leftarrow\frac{N_k}{N}\tag{40} \end{eqnarray}$

ただし、

$\begin{eqnarray} &&N_k=\sum_{n=1}^N\gamma(z_{nk})\tag{41}\\ &&\bar{x}_{ki}=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_{nk})x_{ni}\tag{42}\\ &&\bar{\bf x}_k=(\bar{x}_{k1},\ldots,\bar{x}_{kD})^\top\tag{43} \end{eqnarray}$

とします。

4.収束確認
対数尤度を再計算し、前回との差分があらかじめ設定していた収束条件を満たしていなければ2に戻り、満たしていれば終了します。

偉人の名言

f:id:olj611:20210915090612p:plain:w300
大事は寄せ集められた小事によってなされる。
フィンセント・ファン・ゴッホ

参考文献

パターン認識と機械学習下巻 p160-p163
パターン認識と機械学習の学習 p60-p62

参考リンク

9.3.3：混合ベルヌーイ分布のEMアルゴリズム【PRMLのノート】

動画

なし

機械学習基礎理論独習

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