$\gdef \vec#1{\boldsymbol{#1}} \\ \gdef \rank {\mathrm{rank}} \\ \gdef \det {\mathrm{det}} \\ \gdef \Bern {\mathrm{Bern}} \\ \gdef \Bin {\mathrm{Bin}} \\ \gdef \Mn {\mathrm{Mn}} \\ \gdef \Cov {\mathrm{Cov}} \\ \gdef \Po {\mathrm{Po}} \\ \gdef \HG {\mathrm{HG}} \\ \gdef \Geo {\mathrm{Geo}}\\ \gdef \N {\mathrm{N}} \\ \gdef \LN {\mathrm{LN}} \\ \gdef \U {\mathrm{U}} \\ \gdef \t {\mathrm{t}} \\ \gdef \F {\mathrm{F}} \\ \gdef \Exp {\mathrm{Exp}} \\ \gdef \Ga {\mathrm{Ga}} \\ \gdef \Be {\mathrm{Be}} \\ \gdef \NB {\mathrm{NB}} \\ \gdef \indep {\mathop{\perp\!\!\!\!\perp}} \\ \gdef \tr {\mathrm{tr}}$
EMアルゴリズム扱う場面が『欠損データが存在する時』である場合を本記事（例題1）では扱いました。
今回はそれが『潜在変数が想定される時』である場合を扱います。
（扱う問題）
まず、確率変数$\vec X, \vec Z$を
・$\vec X$：観測データ
・$\vec Z$：潜在変数データ
とします。
具体例として、東京都で1000人をランダムサンプルして身長データ（$X_i$）を測定したものとすると、
$\vec X = (X_1,\cdots,X_{1000}), \vec Z = (Z_1, \cdots, Z_{1000})$
となります。
この時、
$$\begin{aligned} X_i &\overset{i.i.d} \sim \pi \cdot \N(\mu_1,\sigma_1^2) + (1-\pi) \cdot \N(\mu_2,\sigma_2^2) \\ &{\scriptsize(i=1,\cdots,1000; 0 \leqq \pi \leqq 1)} \end{aligned}$$と混合正規分布でモデリングしたものとしましょう。
イメージとしては$\pi$が男性の割合、$\N(\mu_1,\sigma_1^2)$が男性の身長の分布、$(1-\pi)$が女性の割合、$\N(\mu_2,\sigma_2^2)$が女性の身長の分布、といった具合ですね。
グルーピングは性別に限りませんが、性別によってグルーピングするという状況を想定してその様にイメージしていただいても構いません。
（上記問題を最尤法で直接的に解こうとすると）
$$\begin{aligned} p(\vec X | \theta_1, \theta_2, \pi) &= \prod_{i} \{ \pi f(X_i | \theta_1) + (1-\pi) f(X_i | \theta_2) \} \\ &{\scriptsize(但し、\theta_i=(\mu_i,\sigma_i)とまとめた。fは正規分布の確率密度関数を表す。)} \\[10px] \Rightarrow \log p(\vec X | \theta_1, \theta_2, \pi) &= \sum_i \log \{ \pi f(X_i | \theta_1) + (1-\pi) f(X_i | \theta_2) \} \end{aligned}$$となります。
この時、$\log p(\vec X | \theta_1, \theta_2, \pi)$を$\theta_1,\theta_2$でそれぞれ偏微分して$=0$とした連立方程式を解くことになります。
しかしながら、$\log p(\vec X | \theta_1, \theta_2, \pi)$を$\pi$に注目してみると$\log(\square+\triangle)$の形（いわゆるlog-sumの形）になっており、その計算は困難です。
（上記問題をEMアルゴリズムで解いてみる）
上記モデルは、$Z_i$を潜在変数として、
・$X_i | Z_i = 1 \overset{i.i.d} \sim \N(\mu_1, \sigma_1^2)$
・$X_i | Z_i = 2 \overset{i.i.d} \sim \N(\mu_2, \sigma_2^2)$
・$p(Z_i=1) = \pi$
・$p(Z_i=2) = 1-\pi$
というモデルと同等です。
この時、
$$\begin{aligned} p(Z_i=1 | X_i,\theta_1,\theta_2,\pi) &= \frac{ p(Z_i=1,X_i | \theta_1,\theta_2,\pi) }{ p(X_i | \theta_1, \theta_2, \pi) } \\[10px] &= \frac{ \pi f(X_i | \theta_1) }{ \pi f(X_i | \theta_1) + (1-\pi) f(X_i | \theta_2) } \\[20px] p(Z_i=2 | X_i,\theta_1,\theta_2,\pi) &= \frac{ (1-\pi) f(X_i | \theta_2) }{ \pi f(X_i | \theta_1) + (1-\pi) f(X_i | \theta_2) } \\ {\scriptsize(こちらは上記と同様なので最初の行は省略)} \end{aligned}$$となります。
よって、Q関数$Q((\theta_1,\theta_2,\pi) | (\theta_1^*,\theta_2^*,\pi^*))$は、
$$\begin{aligned} Q((\theta_1,\theta_2,\pi) | (\theta_1^*,\theta_2^*,\pi^*)) &= E_{\vec Z | \vec X,\theta_1^*,\theta_2^*,\pi^*} [\log p(\vec X, \vec Z | \theta_1, \theta_2, \pi)] \\[10px] &= E_{\vec Z | \vec X,\theta_1^*,\theta_2^*,\pi^*} [\sum_i \log p(X_i, Z_i | \theta_1, \theta_2, \pi)] \\ &{\scriptsize(独立性より)} \\[10px] &= \underbrace{\sum_i \log p(X_i,Z_i | \theta_1, \theta_2, \pi)}_{(*)} \times p(Z_i | X_i, \theta_1^*, \theta_2^*, \pi^*) \\ &{\scriptsize((*)はZ_j(j \neq i)に依存せず、また、p(Z_i | \vec X, \theta_1^*, \theta_2^*, \pi^*) = p(Z_i | X_i, \theta_1^*, \theta_2^*, \pi^*)であることより)} \\[10px] &= \sum_i \{ \underbrace{p(Z_i=1 | X_i, \theta_1^*, \theta_2^*, \pi^*)}_{:=C_i^{(1)}} \cdot \log p(X_i,Z_i=1 | \theta_1,\theta_2,\pi) \\ &+ \underbrace{p(Z_i=2 | X_i, \theta_1^*, \theta_2^*, \pi^*)}_{:=C_i^{(2)}} \cdot \log p(X_i,Z_i=2 | \theta_1,\theta_2,\pi) \} \\[10px] &= \sum_i \{ C_i^{(1)} \cdot \log p(X_i,Z_i=1 | \theta_1,\theta_2,\pi) + C_i^{(2)} \cdot \log p(X_i,Z_i=2 | \theta_1,\theta_2,\pi) \} \\[10px] &= \sum_i \{ C_i^{(1)} \cdot \log \pi f(X_i | \theta_1) + C_i^{(2)} \cdot \log (1-\pi) f(X_i | \theta_2) \} \end{aligned}$$となります。
これより、
$$\begin{aligned} \frac{\partial Q}{\partial \pi} &= \sum_i \{\frac{C_i^{(1)}}{\pi}-\frac{C_i^{(2)}}{1-\pi}\} = 0 (として) \\[10px] \Rightarrow \pi &= \frac{ \sum_i C_i^{(1)} }{ \sum_i (C_i^{(1)} + C_i^{(2)}) } \\[10px] &= \frac{C_i^{(1)}}{n} \\ &{\scriptsize(n=1000である)} \end{aligned}$$となります。
これは$\theta_1^*,\theta_2^*,\pi^*$が与えられた下での$Z_i=1$となる事後確率に相当します。
また、
$$\begin{aligned} \frac{\partial Q}{\partial \mu_1} &= \sum_i (\frac{X_i-\mu_1}{\sigma_1}) C_i^{(1)} = 0 (として) \\[10px] \Rightarrow \mu_1 &= \frac{\sum_i C_i^{(1)} X_i}{\sum_i C_i^{(1)}} \\[40px] \frac{\partial Q}{\partial \sigma_1} &= \sum_i (-\frac{1}{\sigma_1} + \frac{(X_i-\mu_1)^2}{\sigma_1^3}) C_i^{(1)} = 0 (として) \\[10px] \sigma_1^2 &= \frac{ \sum_i C_i^{(1)} (X_i-\mu_1)^2 }{ \sum_i C_i^{(1)} } \end{aligned}$$となります。
（$\mu_2,\sigma_2$は同様なので省略）
求めた$\mu_1,\mu_2$はそれぞれの事後分布における重み付け平均になります。
また、求めた$\sigma_1^2,\sigma_2^2$はそれぞれの事後分布における重み付き分散になります。
ここで計算してきた様に、潜在変数を想定できる場合にEMアルゴリズムを用いることでその計算過程が非常にシンプルな結果で表されます。
ただしもちろん問題設定によっては複雑になることもあります。
