[general] EM 알고리즘(Expectation and Maximization Algorithm)

이번 글에서는 최대 우도 추정법(MLE)을 사용할 수 없는 경우, 모수를 추정하기 위해 사용할 수 있는 EM 알고리즘에 대해 알아보겠습니다. EM 알고리즘은 일반적인 모델에서 쓰일 수 있지만, 이번 글에서는 가우시안 혼합 모델의 모수를 추정하기 위해 쓰는 경우에 대해 알아보겠습니다. 따라서 이전 글가우시안 혼합 모델의 확장(Extention ofGaussian Mixture)을 읽고 오시면 도움이 될 것입니다.

---

가우시안 혼합 모델의 주변 분포와 로그 우도는 다음과 같습니다.

$\begin {align} &p(x_n)=\displaystyle \sum_{k=1}^{K}{\pi_k N(x_n|\mu_k,~\Sigma_k)} \\ &ln~ p(\textbf {X}|\mathbf {\pi}, \mathbf {\mu}, \mathbf {\Sigma}) = \displaystyle \sum^N_{n=1}{ln~ \left \{ \displaystyle \sum^K_{k=1}{\pi_k N(x_n|\mu_k,~\Sigma_k)} \right \} }\end {align}$

최대 우도 추정법을 이용하여 가우시안 혼합 모델의 모수를 추정하는 것은 크게 두 가지 어려움이 있습니다. 

첫 번째는 미분 값이 닫힌 형식으로 표현되지 않을 수 있습니다.

두 번째는 만약 관측치가 하나밖에 없는 분포를 혼합한다면, 로그 우도의 식이 무한대로 발산하게 됩니다. 두 번째 경우를 식으로 나타내면 다음과 같습니다.

관측치가 하나밖에 없는 경우 평균은 그 관측치가 되고, 분산은 0이 되는데 그때의 가우시안 모델의 식은 다음과 같습니다.

$N(x_n|x_n,~\sigma^2_j\textbf {I})=\frac {1}{(2\pi)^{1/2}}\frac {1}{\sigma_j}$

따라서 가우시안 혼합 모델과 같이 최대 우도 추정법으로 모수를 추정할 수 없는 경우, 기울기 기반의 최적화 기법을 사용하거나, EM알고리즘을 사용하여 모수를 추정할 수 있습니다. 

---

EM알고리즘은 우리가 추정하고 싶은 모수의 초기값을 정한 후, 반복적인 추론을 통하여 모수를 갱신하는 방법입니다. 모수를 반복적으로 갱신하기 위해, 각각의 모수를 정리하여 나타내면 편할 것입니다.

먼저 위의 로그 우도 식을 우리가 찾고 싶은 모수$\mu_k$에 대해 미분하여 풀면 다음과 같습니다.

$\mu_k=\frac {1}{N_k}{\displaystyle \sum^N_{n=1}{\gamma(z_{nk})x_n}}$

여기서 $\gamma(z_{nk})$와 $N_k$는 다음과 같습니다.

$\begin {align} \gamma(z_k) \equiv p(z_k=1|\textbf {x}) &= \frac{{p(z_k=1) p(\textbf {x}|z_k=1)}}{\displaystyle \sum^{K}_{k=1}{p(z_j=1) p(\textbf {x}|z_j=1)}} \\ &= \frac {\pi_k N(\textbf {x}|\mu_k,~\Sigma_k)}{\displaystyle \sum^{K}_{j=1}{\pi_j N(\textbf {x}|\mu_j,~\Sigma_j)}} \end {align}$

$N_k=\displaystyle \sum^N_{n=1}{\gamma(z_{nk})}$

마찬가지로 위의 로그 우도 식을 우리가 찾고 싶은 모수 $\Sigma_k$에 대해 미분하여 풀면 다음과 같습니다.

$\Sigma_k=\frac {1}{N_k}{\displaystyle \sum^K_{n=1}{\gamma(z_{nk})(x_n-\mu_k)(x_n-\mu_k)^T}}$

그리고 mixing coefficients $\pi_k$의 합이 1이라는 것을 이용하여 라그랑주 승수법을 이용하여 구하면 다음과 같습니다.

$\pi_k=\frac {N_k}{N}$

위의 식들을 이용하여 모수를 반복적으로 갱신하게 되는데 EM알고리즘을 정리하면 다음과 같습니다.

---

1. $\mu_k,~\Sigma_k,~\pi_k$에 대한 초기값을 설정하고, 초기값들을 통해 로그 우도를 계산합니다.

2. E Step. 초기값을 이용하여 responsibility를 계산합니다.

$\begin {align} \gamma(z_k) \equiv p(z_k=1|\textbf {x}) &= \frac{{p(z_k=1) p(\textbf {x}|z_k=1)}}{\displaystyle \sum^{K}_{k=1}{p(z_j=1) p(\textbf {x}|z_j=1)}} \\ &= \frac {\pi_k N(\textbf {x}|\mu_k,~\Sigma_k)}{\displaystyle \sum^{K}_{j=1}{\pi_j N(\textbf {x}|\mu_j,~\Sigma_j)}} \end {align}$

3. M Step. responsibility 값을 통해 새로운 모수를 계산합니다.

$ \begin {align} &\mu_k^{\textrm {new}}=\frac {1}{N_k}{\displaystyle \sum^N_{n=1}{\gamma(z_{nk})x_n}} \\ &\Sigma_k=\frac {1}{N_k}{\displaystyle \sum^K_{n=1}{\gamma(z_{nk})(x_n-\mu_k^{\textrm {new}})(x_n-\mu_k^{\textrm {new}})^T}} \\ &\pi_k^{\textrm {new}}=\frac {N_k}{N} \end {align}$

4. 로그 우도 값을 계산한 후, 수렴하면 알고리즘을 멈추고, 아닌 경우 다시 2. 번으로 돌아갑니다.

---

가우시안 혼합 모델에서 EM알고리즘으로 갱신하는 과정을 그림으로 나타내면 다음과 같습니다.

1. 과정에서 설정된 초기값을 통해, responsibility(2. 과정)를 계산하면 (b)와 같습니다. 그리고 responsibility를 통해 다시 모수를 갱신하여 나타내면 그림 (c)와 같습니다. 이 과정을 각각 2번, 5번, 20번 반복한 그림은 (d), (e), (f)와 같습니다.

과정이 진행될수록 가우시안 혼합 모델의 모수들이 잘 적합되는 모습을 볼 수 있습니다.