EMアルゴリズム（expectation maximization algorithm）とは、統計学において、確率モデルのパラメータを最尤推定する手法の一つであり、観測不可能な潜在変数に確率モデルが依存する場合に用いられる。
混合ガウス分布の何番目からデータが発生したのかを表す隠れ変数Zを事前に導入する。
混合ガウス分布
具体例で学ぶ数学 様 「混合ガウス分布の意味と役立つ例」の説明
https://mathwords.net/gmm
ガウス分布 \(N(x\mid\mu_k,\Sigma_k):k=1,2,\cdots,K\) に対し、
\(p(x)=\displaystyle\sum_{k=1}^K\pi_kN(x\mid\mu_k,\Sigma_k)\) \(\pi_k\) は混合率：\(\displaystyle\sum_{k=1}^K\pi_k=1\)
によって定まる分布を混合ガウス分布という。
［目的］ 与えられたデータに当てはまる混合ガウス分布を推定する。
データ \((x_1,x_2,\cdots,x_N)\) が与えられたとき、
最尤法により、 \(\displaystyle\prod_{i=1}^Np(x_i)\) を最大にする \((\pi_k,\mu_k,\Sigma_k:\ k=1,2,\cdots,K)\) を求める。
但し、 \(p(x_i)=\pi_1N(x_i\mid\mu_1,\Sigma_1)+\pi_2N(x_i\mid\mu_2,\Sigma_2)+\cdots+\pi_KN(x_i\mid\mu_K,\Sigma_K)\)
しかし、混合率 \(\pi_k\) とパラメーター \(\mu_k,\Sigma_k\) を同時に求めることは一般に困難である。
そこで、EMアルゴリズムという逐次近似法（混合率→パラメータ→混合率→パラメータ→‥）
により、収束する値を最適解とする。
［EMアルゴリズムの考え方］
数学的裏付けを離れてイメージで記す。（論理的正確さを犠牲にしています）
・ガウス分布 \(N(x\mid\mu_k,\Sigma_k)\) が決まっていても、混合率 \(\pi_k\) が変われば分布は変化する。
・混合率 \(\pi_k\) が決まっていても、パラメータ \(\mu_k,\Sigma_k\ :\ k=1,2,\cdots,K\) が変われば分布は変化する。
① 混合率 \(\pi_k\) の初期値を決め、 \((\mu_k,\Sigma_k:\ k=1,2,\cdots,K)\) を最尤法により推定する。
② 求まった \((\mu_k,\Sigma_k:\ k=1,2,\cdots,K)\) のもとで、混合率 \(\pi_k^{new}\) を最尤法で推定する。
③ 新しい混合率 \(\pi_k^{new}\) のもとで、パラメータ \((\mu_k,\Sigma_k:\ k=1,2,\cdots,K)^{new}\) を推定する。
これらを繰り返すことで、\((\pi_k,\mu_k,\Sigma_k:\ k=1,2,\cdots,K)\) が収束すれば、それを最適解とする。
［データ \(x_i\) の負担率］ \(\gamma_{ik}=\dfrac{\pi_kN(x_i\mid \mu_k,\Sigma_k)}{\displaystyle\sum_{k=1}^K \pi_kN(x_i\mid \mu_k,\Sigma_k)}\) \(p(x_i)\) のうち \(k\) 番目のガウス分布の割合
負担率は \(x_i\) に依存しており \(x_i\) によって値が変わる。
\(\pi_k^{new}=f_1(\gamma_{ik})=f_1(\gamma_{ik},\pi_k)\) 、\(\mu_k^{new}=f_2(\gamma_{ik},\mu_k)\) 、 \(\Sigma_k^{new}=f_3(\gamma_{ik},\Sigma_k)\)
\(\displaystyle \sum_{k=1}^K\gamma_{i,k}=1\) であるから、 \(\gamma_{i,k}\) は混合率と考えられる。
慶應義塾大学研究員 皆川 卓也 様（PRML読書会）の説明
https://www.slideshare.net/takmin/chapter9-2
［負担率］
データ x が与えられたときの z の条件付確率を負担率と呼ぶ
混合ガウス分布 \(\displaystyle p(x) = \sum_{k=1}^K\pi_k N(x|\mu_k , \Sigma_k)\)
尤度関数 \(p(X|\pi,\mu , \Sigma ) =\displaystyle \prod_{n=1}^N \sum_{k=1}^K\pi_k N(x_n|\mu_k , \Sigma_k)\)
［潜在変数の導入］
確率 \(\displaystyle p(x)=\sum_zp(x,z)=\sum_zp(z)p(x|z)\)
潜在変数 \(Z=\begin{pmatrix}
z_1\\
\vdots\\
z_K
\end{pmatrix}\) \(\left(z_k=0,1\ :\sum_{k=1}^Kz_k=1\right)\)
事前分布 \(p(z_k=1)=\pi_k\) 、 \(\displaystyle p(z)=\prod_{k=1}^K\pi_k^{z_k}\) （ \(\displaystyle 0<\pi_k<1:\sum_{k=1}^K\pi_k=1\) 混合率 ） (9.10)
事後分布 \(p(x|z_k=1) =N(x|\mu_k , \Sigma_k)\) 、\(p(x|z) =\displaystyle \prod_{k=1}^KN(x|\mu_k , \Sigma_k)^{z_k}\) (9.11)
このとき、 \(\displaystyle p(x) = \sum_{k=1}^K\pi_k N(x|\mu_k , \Sigma_k)\) (9.7)
\(\begin{split}\because\
p(x)&=\sum_zp(x,z)=\sum_zp(z)p(x|z)=\sum_z\left(\prod_{k=1}^K\pi_k^{z_k}\displaystyle \prod_{k=1}^KN(x|\mu_k , \Sigma_k)^{z_k}\right)\\
&=\sum_z\prod_{k=1}^K\left(\pi_kN(x|\mu_k , \Sigma_k)\right)^{z_k}
=\sum_z\pi_kN(x|\mu_k , \Sigma_k)
\end{split}\)
負担率：x が与えられたもとでの \(z_k\) の条件付き確率 （スライド11）
\(\gamma(z_k)=p(z_k=1|x)=\dfrac{p(z_k=1,x)}{p(x)}=\dfrac{p(z_k=1)p(x|z_k=1)}{\sum_{j=1}^Kp(z_j=1)p(x|z_j=1)}
=\dfrac{\pi_kN(x|\mu_k , \Sigma_k)}{\sum_{j=1}^K\pi_jN(x|\mu_j , \Sigma_j)}\)
（データ \(x\) が与えられたとき、データがグループ \(k\) 由来である確率）
・ \(\pi_k\) は事前分布 （ \(x\) が観測される前の）
・ \(\gamma(z_k)\) は事後分布 （ \(x\) が観測された後の）
x が与えられたとき、
\((\pi_k,\mu_k,\Sigma_k:\ k=1,2,\cdots,K)\rightarrow\gamma(z_k)\rightarrow
(\pi_k,\mu_k,\Sigma_k:\ k=1,2,\cdots,K)^{new}\)
リクルート 小山田 創哲 様 の説明
https://www.slideshare.net/sotetsukoyamada/em-36315279
EMアルゴリズムとは、隠れ変数を含むモデルの学習に使われるアルゴリズムである。
混合分布（Gaussian mixture）のモデリングには２通りの方法が考えられる。
［ 方法１］混合率 \(\pi_k\) を用いるモデル（隠れ変数を使わない観測変数 x だけのモデル）
\(p(x)=p(x|\pi,\mu,\Sigma)\) を確率関数として、最尤法により \((\pi_k,\mu_k,\Sigma_k:\ k=1,2,\cdots,K)\) を求める。
（ところが、この方法では最尤解が定まらないことがわかる）
［ 方法２］グループを決定する確率変数 z を用いるモデル（隠れ変数zを使うモデル）
\(\displaystyle p(x)=\sum_zp(x|z,\mu,\Sigma)p(z|\pi)\) \(p(z|\pi)\) は、混合率 のもとで が起こる確率
\(x_{1:N}=\{x_1,x_2,\cdots,x_N\}\) とおくと、
\(\displaystyle p(x_{1:N}|\pi,\mu,\Sigma)=\prod _{n=1}^N p(x_{n}|\pi,\mu,\Sigma)
=\displaystyle \prod_{n=1}^N \sum_{k=1}^K\pi_k N(x_n|\mu_k , \Sigma_k)\)
\(\displaystyle \log p(x_{1:N}|\pi,\mu,\Sigma)
=\sum_{n=1}^N \log\left( \sum_{k=1}^K\pi_k N(x_n|\mu_k , \Sigma_k)\right)\)
もし単体のガウス分布（k=1)であれば、 \(x_n\) が \(D\) 次元とすると
\(\begin{split}
\log p(x_{1:N}|\pi,\mu,\Sigma)
&=\sum_{n=1}^N \log N(x_n|\mu_k , \Sigma_k)\\
&=- \frac{ND}{2} \log (2 \pi ) - \frac{N}{2} \log \left[ \mathrm{det} \left( \Sigma \right) \right] - \sum_{n=1}^N \frac{1}{2} ( x_n - \mu )^T \Sigma^{-1} ( x_n - \mu )
\end{split}\)
\(\displaystyle \frac{ \partial }{ \partial \mu } \log p(x_{1:N}|\pi,\mu,\Sigma)
=-\sum_{n=1}^N \left(\Sigma^{-1} (x_n -\mu)\right)
=-\Sigma^{-1}\sum_{n=1}^N (x_n -\mu)
=
0\) より
\(\displaystyle \sum_{n=1}^N (x_n -\mu)=\sum_{n=1}^N x_n -N\mu=0\) \(\displaystyle \therefore\: \widehat\mu=\dfrac{1}{N}\sum_{n=1}^N x_n\) （平均の最尤推定量）
計算は難しいが、 \(\Sigma\) の最尤推定量は、 \(\displaystyle \widehat \Sigma = \frac{1}{N} \sum_{n=1}^N ( x_n - \hat{ \mu }) ( x_n - \hat{ \mu })^T \) となる。
混合ガウス分布の場合は
\(\displaystyle \dfrac{ \partial }{ \partial \mu_k } \log p(x_{1:N}|\pi,\mu,\Sigma)
=\sum_{n=1}^N \dfrac{\pi_k N(x_n|\mu_k , \Sigma_k)}{\sum_{k=1}^K\pi_k N(x_n|\mu_k , \Sigma_k)}\Sigma^{-1} (x_n -\mu_k)
=\sum_{n=1}^N \gamma_{n.k}\ \Sigma^{-1} (x_n -\mu_k)\)
合成関数の log 微分なので分数になり、負担率が現れるところがポイント。
\(\dfrac{ \partial }{ \partial \mu_k } \log p(x_{1:N}|\pi,\mu,\Sigma)=0\) とすると、
\(\displaystyle \Sigma^{-1} \left(\sum_{n=1}^N \gamma_{n,k}x_n - \sum_{n=1}^N\gamma_{n,k}\mu_k\right)=0\) \(\displaystyle \therefore \sum_{n=1}^N \gamma_{n,k}x_n - \sum_{n=1}^N\gamma_{n,k}\mu_k=0\)
\(\displaystyle \sum_{n=1}^N \gamma_{n,k}x_n - \mu_k\sum_{n=1}^N\gamma_{n,k}=0\) より、 \(\displaystyle \sum_{n=1}^N \gamma_{n,k}x_n - \mu_kN_k=0\) \(\displaystyle \therefore \mu_k=\dfrac{1}{N_k}\sum_{n=1}^N \gamma_{n,k}x_n\)
最尤推定値にも負担率が入り、 \(\displaystyle \widehat\mu_k=\dfrac{1}{N_k}\sum_{n=1}^N \gamma_{n,k}x_n\) となる。
負担率 \(\gamma(z_k)=p(z_k=1|x)
=\dfrac{\pi_kN(x|\mu_k , \Sigma_k)}{\sum_{j=1}^K\pi_jN(x|\mu_j , \Sigma_j)}\) は \(\pi_k,\mu_k,\Sigma_k,x\) によって決まる。
上の式では、 \(\gamma_{n,k}=p(z_k=1|x_n)
=\dfrac{\pi_kN(x_n|\mu_k , \Sigma_k)}{\sum_{j=1}^K\pi_jN(x_n|\mu_j , \Sigma_j)}\) である。
さて、
データ \((x_1,x_2,\cdots,x_N)\) が与えられたとする。（各 \(x_n\) は多次元データ）
\(n\) 番目のデータの \(k\) 群への負担率を \(\gamma(z_{n,k})\) とする。
最尤法により、
\(\displaystyle N_k^{new}=\sum_{n=1}^N\gamma(z_{n,k})\) k群のデータ数に対応し、 \(\pi_k^{new}=\dfrac{N_k}{N}\)
\(\displaystyle \mu_k^{new}=\dfrac{1}{N_k}\sum_{n=1}^N\gamma(z_{n,k})\cdot x_n\) \(\displaystyle \Sigma_k^{new}=\dfrac{1}{N_k}\sum_{n=1}^N\gamma(z_{n,k})(x_n-\mu_k)(x_n-\mu_k)^T\)
newパラメータから負担率 \(\gamma(z_{n,k})^{new}\) が決まり、負担率から次のパラメータが定まる。
これを繰り返すことによる収束先が最適解である。
最尤法であるから、与えられるデータ \((x_1,x_2,\cdots,x_N)\) によって負担率 \(\gamma(z_{n,k})\) は変わってよい。
| 1群 | 2群 | \(\cdots\) | K群 | 計 | |
| \(x_1\) | \(\gamma(z_{1,1})\) | \(\gamma(z_{1,2})\) | \(\cdots\) | \(\gamma(z_{1,K})\) | 1 |
| \(x_2\) | \(\gamma(z_{2,1})\) | \(\gamma(z_{2,2})\) | \(\cdots\) | \(\gamma(z_{2,K})\) | 1 |
| \(\vdots\) | \(\vdots\) | \(\vdots\) | \(\vdots\) | 1 | |
| \(x_N\) | \(\gamma(z_{N,1})\) | \(\gamma(z_{N,2})\) | \(\cdots\) | \(\gamma(z_{N,K})\) | 1 |
| 計 | \(N_1\) | \(N_2\) | \(\cdots\) | \(N_K\) | \(N\) |
\(\gamma(z_{n,k})\) は加重と同じ働きをすると考えれば
\(\displaystyle \mu_k=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_{n,k})x_n\) \(\displaystyle \Sigma_k=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_n,k)(x_n-\mu_k)(x_n-\mu_k)^T\) は理解できる。
［参考］
| 1群 | 2群 | \(\cdots\) | K群 | 計 | |
| \(x_1\) | \(\pi_1 N(x_1|\mu_1 , \Sigma_1)\) | \(\pi_2 N(x_1|\mu_2 , \Sigma_2)\) | \(\cdots\) | \(\pi_K N(x_1|\mu_K , \Sigma_K)\) | \(p(x_1)\) |
| \(x_2\) | \(\pi_1 N(x_2|\mu_1 , \Sigma_1)\) | \(\pi_2 N(x_2|\mu_2 , \Sigma_2)\) | \(\cdots\) | \(\pi_K N(x_2|\mu_K , \Sigma_K)\) | \(p(x_2)\) |
| \(\vdots\) | \(\vdots\) | \(\vdots\) | \(\vdots\) | \(\vdots\) | |
| \(x_N\) | \(\pi_1 N(x_N|\mu_1 , \Sigma_1)\) | \(\pi_2 N(x_N|\mu_2 , \Sigma_2)\) | \(\cdots\) | \(\pi_K N(x_N|\mu_K , \Sigma_K)\) | \(p(x_N)\) |
Re:ゼロから始めるML生活 様
https://www.nogawanogawa.com/entry/em_algorithm
［データの形成］
混合モデルだと仮定したとき、この試行において結果が出力されるまでの過程
・第一段階としてどちらの分布に従うかを選択し、
・第二段階として選択された分布の中で確率的にどこに位置するかが選択される
・潜在変数Zによって分布が選択され、（負担率）
・同時に観測変数によって分布が形成されることを表現できます。
\(\displaystyle p(x|\mu , \sigma, \pi) = \sum_j\pi_j N(x|\mu , \sigma)\)
観測できないπが含まれており、観測変数と潜在変数を同時に知ることができない不完全データ
潜在変数を確率変数として取り扱うことで尤度の期待値を最大化するというアプローチを考える
ことで、近似的にパラメータを推定する方法がEMアルゴリズムになります。
分布のパラメータと負担率は、一方が他方さえ決定されれば推定できるという状況になっている
・Eステップ(期待値計算)
分布のパラメータθを固定し、負担率を計算
・Mステップ(最大化)
負担率を固定して、対数尤度が最大化するパラメータθを求める。