EM 알고리즘에 대하여 Optimization 강의 1

이번 포스팅에서는 머신러닝(machine learning) 분야에서 많이 등장하는 EM 알고리즘(Expectation–Maximization algorithm)에 대하여 알아보도록 하겠다. 요즘 머신러닝이 유행하면서 EM 알고리즘이 뭔가 새롭게 등장한 것 같지만, EM 알고리즘은 어찌보면 전통 통계기법이라 할 만큼 그 역사가 길다.

Arthur P. Dempster (from Wiki)

EM 알고리즘의 이름이 처음 등장한 것은 1977년 Arthur Dempster, Nan Laird, and Donald Rubin이 쓴 논문에 처음 등장하지만, 이 논문의 요지는 ’이제까지 특정 분야들에서 여러 사람들에 의하여 제안되었던 방법들이 하나의 공통적인 과정을 따르고 있다.’라는 것이었다. 즉, 1977년 이전부터 쭉 사람들이 EM 알고리즘을 사용하여 문제를 해결했지만, 각자의 상황들이 달라서 같은 알고리즘이라도 인식하지 못해왔다는 것이다. 그렇다면 1977년의 논문은 실제로 독창적인 내용이 담겨져있는 것이 아니라고 생각할 수도 있겠지만, 예를들어 논문 A와 논문 B가 같은 알고리즘을 사용했다 할 지라도, 실제로 논문마다 주어진 상황이 다르고, 각자가 알고리즘을 설명하는 방식이 다르다면, A와 B의 알고리즘이 결국 큰 틀안에서 보면 동일하다는 것을 인지하기는 결코 쉬운일이 아니다. 이건 마치 수학에서의 일반화가 어려운 것 처럼 컴퓨터 사이언스 분야에서의 일반화 역시 동일하다고 생각하면 좋을 것 같다. 1977년의 논문은 실제로 엄청난 학계의 관심을 불러일으켰으며, 이 방식은 오늘날까지도 여전히 여러 분야에서 사용되고 있다.

 

EM 알고리즘은 언제 사용할 수 있을까

EM 알고리즘은 결국 우리가 앞에서 다루었던 우도함수를 최대화 시키는 모수를 찾는 방법, 즉, MLE를 찾는 방법 중 하나이다. 예를들어, 우리가 확률변수 \(Y\sim g_{\theta}\left(y\right)\)를 관찰하고 있다고 생각하자. 여기서 \(g\)는 \(Y\)의 확률밀도함수이다. 따라서 우리가 \(y_1,...,y_n\)의 데이터를 가지고 있다면, \(Y\)의 모수 \(\theta\)에 대한 MLE는 이전 포스팅에서 살펴보았듯, 로그우도 함수를 최대로 만드는 \(\theta\)의 값이 된다.

\[ \hat{\theta}=\underset{\theta}{argmax}\,\sum_{i=1}^{n}log\,g_{\theta}\left(y_{i}\right) \]

하지만 위의 로그우도 함수를 직접 최대화하기 힘든 경우가 많다. 이럴 경우 하나의 대안책으로 EM 알고리즘을 사용할 수 있다. 알고리즘을 적용할 수 있는 조건은 다음과 같다.

1. 우리가 관찰하는 확률변수 Y가 우리가 잘 알고있으며 같은 모수를 사용하는 좀 더 단순한 확률변수 \(X\sim p_{\theta}\left(x\right)\)로부터 나왔다고 볼 수 있는 경우, 즉, \(Y=T(X)\)같은 어떤 함수 T를 만들 수 있는 경우.
2. 확률변수 X를 사용한 우도 함수의 최대화를 Y를 사용하는 위의 경우보다 비교적 쉽게 할 수 있는 경우.

우리가 풀려고 하는 문제가 위의 두 가지 조건을 만족할 때 EM 알고리즘은 다음과 같은 방법을 사용하여 모수 \(\theta\)의 추정 할 수 있다.

EM 알고리즘

주어진 데이터를 \(\underline{y}\), \(\theta\)에 대한 초기값을 \(\theta^{\left(0\right)}\)라고 하자. 그렇다면 \(k=1,2,...\)에 대하여 다음의 두 과정을 반복한다.

-E-step: 주어진 데이터 \(\underline{y}\)와 현재 \(\theta\)의 추정치 \(\theta^{\left(k\right)}\)를 사용하여 다음과 같은 함수를 계산한다.

\[ Q\left(\theta;\theta^{\left(k\right)}\right)=\mathbb{E}\left[log\,p_{\theta}\left(X\right)|Y=\underline{y},\theta^{\left(k\right)}\right] \]

-M-step: 위에서 계산한 \(Q\)함수를 \(\theta\)에 대하여 최대화한다.

\[ \theta^{\left(k+1\right)}=\underset{\theta}{argmax}\,Q\left(\theta;\theta^{\left(k\right)}\right) \]

예제를 통한 이해

EM 알고리즘을 위와 같이 식으로만 써져있는 것을 보고 이해가 안된다고 좌절하면 안된다. 위의 식을 보고 바로 알고리즘을 이해하기란 결코 쉬운 일이 아니다. 다음의 예제를 따라가면서 위의 식을 따라가보자. 문제에서 주어진 상황을 완벽히 이해하고, 예제의 계산과정이 위의 알고리즘 설명에서 어떠한 부분을 의미하는지 계속해서 번갈아가면서 확인해야 완벽히 이해가 가능할 것이라 생각한다.

---

예제(Dempster et al. 1977 변형)

길동이가 운영하는 아이스크림 가게는 초콜렛바닐라, 딸기, 민트, 체리 이렇게 4가지의 아이스크림을 팔고있다. 길동이는 가게에 오는 손님들이 4가지 아이스크림을 다음과 같은 확률로 고른다고 생각하고 있다.

\[ \left(\frac{1}{2}+\frac{1}{4}\theta,\frac{1}{4}\left(1-\theta\right),\frac{1}{4}\left(1-\theta\right),\frac{1}{4}\theta\right) \]

길동이는 하루동안의 매출 기록을 이용하여 \(\theta\)의 값을 알아보고자 한다. 만약 하루동안 n명의 손님이 가게를 다녀갔을때 각 종류별로 팔린 아이스크림의 갯수, 확률변수 \(Y\), 는 다음과 같은 다항분포를 따를 것이다.

\[ \begin{align*} Y & \sim multi\left(n,\underline{p}\right)\\ \underline{p} & =\left(\frac{1}{2}+\frac{1}{4}\theta,\frac{1}{4}\left(1-\theta\right),\frac{1}{4}\left(1-\theta\right),\frac{1}{4}\theta\right) \end{align*} \]

따라서 Y에 대한 실제 모수는 \(\theta\)가 되며, p 벡터는 확률벡터이므로 적절한 \(\theta\)의 값은 \(0 \leq \theta \leq 1\)이 된다. 하루동안 길동이의 가게에서 팔린 아이스크림 별 갯수 \(y\)가 다음과 같을 때 \(\theta\)의 MLE값을 추정해보자.

\[ \begin{align*} \underline{y} & =\left(y_{1},y_{2},y_{3},y_{4}\right)\\ & =\left(125,18,20,34\right) \end{align*} \]

---

위의 예제에서 \(\theta\)의 MLE값을 구하기 위해서 우도 함수를 정의하면 다음과 같다.

\[ g_{\theta}\left(\underline{y}\right)=\frac{n!}{y_{1}!y_{2}!y_{3}!y_{4}!}\left(\frac{1}{2}+\frac{1}{4}\theta\right)^{y_{1}}\left(\frac{1-\theta}{4}\right)^{y_{2}}\left(\frac{1-\theta}{4}\right)^{y_{3}}\left(\frac{1}{4}\theta\right)^{y_{4}}, \]

여기서 \(n\)은 \(\sum_{i=1}^{4}y_{i}\)과 같다. 따라서 로그우도함수는 다음과 같다.

\[ \tag{1} l\left(\theta;\underline{y}\right)=y_{1}log\left(\frac{1}{2}+\frac{1}{4}\theta\right)+y_{2}log\left(\frac{1-\theta}{4}\right)+y_{3}log\left(\frac{1-\theta}{4}\right)+y_{4}log\left(\frac{1}{4}\theta\right) \]

즉, 우리의 목표는 이 함수를 최대로 만드는 \(\theta\)를 찾으면 되는 것이다. 하지만, 처음 로그 안의 \(1/2+1/4\theta\)부분 때문에 계산이 복잡해지는 것을 알 수 있다. 만약 \(1/2\) 부분이 없다면 아주 깔끔하게 계산이 떨어질텐데 말이다. 이 부분이 바로 EM 알고리즘을 적용하기 좋은 부분이 된다.

 

완전한(Complete) 데이터를 사용한 MLE

만약 확률변수 \(X=\left(X_{1},X_{2},X_{3},X_{4},X_{5}\right)\)가 다섯개의 카테고리를 갖는 다항분포를 따른다고 가정하자.

\[ \begin{align*} X & \sim multi\left(n,\underline{p}\right)\\ \underline{p} & =\left(\frac{1}{2},\frac{1}{4}\theta,\frac{1}{4}\left(1-\theta\right),\frac{1}{4}\left(1-\theta\right),\frac{1}{4}\theta\right) \end{align*} \]

그렇다면 확률변수 Y는 X와 다음과 같은 관계를 갖는다고 생각할 수 있다.

\[ \begin{align*} Y & =\left(Y_{1},Y_{2},Y_{3},Y_{4}\right)\\ & =T\left(X\right)\\ & =\left(X_{1}+X_{2},X_{3},X_{4},X_{5}\right) \end{align*} \]

이것을 길동의 예제로 옮겨놓으면 다음과 같다. 만약 고객의 원래 선호도는 초콜렛, 바닐라, 딸기, 민트, 체리와 같은 다섯개의 카테고리로 이루어져있는데, 초콜렛과 바닐라를 좋아하는 고객이 초콜렛바닐라 맛 아이스크림을 사가는 것이라고 생각할 수 있는 것이다. 즉, 길동의 데이터는 다섯개의 카테고리를 갖는 확률변수 X의 데이터를 완전한 데이터라고 한다면 Y의 데이터는 완전 데이터의 일부분이 생략된 불완전 데이터를 관찰하고 있는 상태라고 생각할 수 있다.

만약 길동이 확률변수 X에 대한 데이터, 즉 다섯가지 맛에 대한 데이터를 관찰하고 있다면, 로그우도 함수를 최대로 만드는 \(\theta\)를 구하는 것은 아주 쉬워진다.

확률변수 \(X\)를 따르는 데이터 \(\underline{x}\)가 주어졌을 때, 우도 함수를 정의하면 다음과 같다.

\[ p_{\theta}\left(\underline{x}\right)=\frac{n!}{x_{1}!x_{2}!x_{3}!x_{4}!x_{5}!}\left(\frac{1}{2}\right)^{x_{1}}\left(\frac{1}{4}\theta\right)^{x_{2}}\left(\frac{1-\theta}{4}\right)^{x_{3}}\left(\frac{1-\theta}{4}\right)^{x_{4}}\left(\frac{1}{4}\theta\right)^{x_{5}} \]

따라서 로그우도 함수는 다음과 같다.

\[ \begin{eqnarray*} l\left(\theta;\underline{x}\right) & = & logp_{\theta}\left(\underline{x}\right)\\ & = & h\left(\underline{x},n\right)+x_{2}log\left(\frac{1}{4}\theta\right)+x_{3}log\left(\frac{1-\theta}{4}\right)+x_{4}log\left(\frac{1-\theta}{4}\right)+x_{5}log\left(\frac{1}{4}\theta\right)\\ \tag{2} & = & h\left(\underline{x},n\right)+\left(x_{2}+x_{5}\right)log\left(\frac{1}{4}\theta\right)+\left(x_{3}+x_{4}\right)log\left(\frac{1-\theta}{4}\right) \end{eqnarray*} \]

그러므로, \(\theta\)의 MLE는 다음과 같이 쉽게 구해진다.

\[ \begin{align*} \frac{\partial l}{\partial\theta} & =\left(x_{2}+x_{5}\right)\frac{1}{\theta}-\left(x_{3}+x_{4}\right)\frac{1}{1-\theta}\overset{set}{=}0\\ \tag{3} \hat{\theta} & =\frac{x_{2}+x_{5}}{x_{2}+x_{3}+x_{4}+x_{5}} \end{align*} \]

E-step: Q함수 계산

자 이제 EM 알고리즘을 적용하기 위하여 Q함수를 구해보자.

\[ Q\left(\theta,\theta^{\left(k\right)}\right)=\mathbb{E}\left[l\left(\theta;\underline{X}\right)\vert\underline{Y}=\underline{y},\theta^{\left(k\right)}\right] \]

Q 함수를 구하기 위해서 먼저 주목해야하는 부분은 \(X_3, X_4, X_5\)의 값이 \(Y_2, Y_3, Y_4\)와 같다는 사실이다. 즉, Y의 관측값이 결정되면 우리는 \(X_3, X_4, X_5\) 값을 알고 있는 것과 마찬가지이다. 길동의 예제에 적용해서 다시 말하면, 길동의 데이터를 보고 우리는 \(X_3, X_4, X_5\)의 관측값이 \(18, 20, 34\)라는 사실을 알 수 있다. 따라서 다음의 식이 성립한다.

\[ \mathbb{E}\left[X_{i}\vert\underline{Y}=\underline{y},\theta^{\left(k\right)}\right]=y_{i-1},\mbox{ where i=3,4,5} \]

반면, 우리가 \(Y_1\)의 관측값을 관찰했다고 하더라도 \(X_1\)과 \(X_2\)의 경우는 알 수 없다. 하지만 \(Y_1\)값을 알고있다면, \(X_1\)가 이항분포를 따르게 된다는 것을 짐작할 수 있다. 데이터 \(\underline{y}\)와 \(\theta^{(k)}\)가 주어졌을때 \(X_1\) 조건부 확률질량 함수를 구해보면 다음과 같다.

\[ \begin{eqnarray*} p\left(x_{1}|\underline{y},\theta^{\left(k\right)}\right) & = & \frac{p\left(x_{1},y|\theta^{\left(k\right)}\right)}{p\left(y|\theta^{\left(k\right)}\right)}\\ & = & \frac{\frac{n!}{x_{1}!\left(y_{1}-x_{1}\right)!y_{2}!y_{3}!y_{4}!}\left(\frac{1}{2}\right)^{x_{1}}\left(\frac{1}{4}\theta^{\left(k\right)}\right)^{y_{1}-x_{1}}}{\frac{n!}{y_{1}!y_{2}!y_{3}!y_{4}!}\left(\frac{1}{2}+\frac{1}{4}\theta^{\left(k\right)}\right)^{y_{1}}}\\ & = & \frac{y_{1}!}{x_{1}!\left(y_{1}-x_{1}\right)!}\left(\frac{\frac{1}{2}}{\frac{1}{2}+\frac{1}{4}\theta^{\left(k\right)}}\right)^{x_{1}}\left(\frac{\frac{1}{4}\theta^{\left(k\right)}}{\frac{1}{2}+\frac{1}{4}\theta^{\left(k\right)}}\right)^{y_{1}-x_{1}}. \end{eqnarray*} \]

또한, X_2의 경우 \(X_2 = y_1 - X_1\)의 관계를 만족시키므로,

\[ \begin{align*} X_{1}|\underline{Y} & =\underline{y},\theta^{\left(k\right)}\sim B\left(y_{1},\frac{\frac{1}{2}}{\frac{1}{2}+\frac{1}{4}\theta^{\left(k\right)}}\right)\\ X_{2}|\underline{Y} & =\underline{y},\theta^{\left(k\right)}\sim B\left(y_{1},\frac{\frac{1}{4}\theta^{\left(k\right)}}{\frac{1}{2}+\frac{1}{4}\theta^{\left(k\right)}}\right) \end{align*} \]

따라서 이항분포의 기댓값 공식에 따라서 다음의 조건부 기대값을 얻을 수 있다.

\[ \begin{eqnarray*} \mathbb{E}\left[X_{1}\vert\underline{Y}=\underline{y},\theta^{\left(k\right)}\right] & = & y_{1}\frac{\frac{1}{2}}{\frac{1}{2}+\frac{1}{4}\theta^{\left(k\right)}},\\ \mathbb{E}\left[X_{2}\vert\underline{Y}=\underline{y},\theta^{\left(k\right)}\right] & = & y_{1}\frac{\frac{1}{4}\theta^{\left(k\right)}}{\frac{1}{2}+\frac{1}{4}\theta^{\left(k\right)}}. \end{eqnarray*} \]

따라서 \(\theta^{(k)}\)에 대한 Q함수는 다음과 같이 구할 수 있다.

\[ \begin{eqnarray*} Q\left(\theta,\theta^{\left(k\right)}\right) & = & \mathbb{E}\left[l\left(\theta;\underline{X}\right)\vert\underline{Y}=\underline{y},\theta^{\left(k\right)}\right]\\ & = & \left(y_{1}\frac{\frac{1}{4}\theta^{\left(k\right)}}{\frac{1}{2}+\frac{1}{4}\theta^{\left(k\right)}}+y_{4}\right)log\left(\frac{1}{4}\theta\right)+\\ & & \left(y_{2}+y_{3}\right)log\left(\frac{1-\theta}{4}\right)+h\left(\underline{y},\theta^{\left(k\right)}\right) \end{eqnarray*} \]

M-step: Q함수 최대화

위에서 구한 Q함수의 모양을 잘 들여다보면 위에서 구한 완전 데이터를 사용한 로그우도 함수 식 (2)과 형태가 똑같다는 사실을 알 수 있다. 따라서 위의 Q함수를 최대로 만드는 \(\theta\)의 값은 식 (3)를 이용하여 쉽게 얻을 수 있다.

\[ \theta^{\left(k+1\right)}=\frac{\left(y_{1}\frac{\frac{1}{4}\theta^{\left(k\right)}}{\frac{1}{2}+\frac{1}{4}\theta^{\left(k\right)}}+y_{4}\right)}{y_{1}\frac{\frac{1}{4}\theta^{\left(k\right)}}{\frac{1}{2}+\frac{1}{4}\theta^{\left(k\right)}}+y_{2}+y_{3}+y_{4}} \]

R 코드를 사용한 구현

이제 EM 알고리즘을 R에서 구현할 준비가 끝났다. 다음의 EM 함수는 E스텝과 M스텝을 모수 \(\theta\)의 값의 변화가 거의 없을때까지 반복하도록 되어 있다.

EM <- function(theta, y){
	# Initial values
    iter <- 0
    value <- NULL
    
    # repeat E and M steps
    repeat{
    	## E step
        p <- 0.5 / ( 0.5 + 0.25 * theta )
        x <- c( y[1]*p , y[1]*(1-p), y[2:4] )
        ## M step
        theta_new <- ( x[2] + x[5] ) / ( x[2] + x[3] + x[4] + x[5] )
        ## update old to new
        iter <- iter + 1
        value[iter] <- theta_new
        epsilon  <- abs( theta_new - theta )
        theta    <- theta_new
        ## decision based on
        if (epsilon < .Machine$double.eps){ break }
	}
    # pack up result by converting theta to r
    list(theta = theta_new,
    	 step  = iter,
         whole = value )
}

위의 길동의 데이터를 사용하여 위의 함수를 실행시켜보자. \(\theta\)의 초기값은 임의로 0.5를 부여하였다.

y <- c( 125, 18, 20, 34 )
resultEM <- EM( 0.5, y )
resultEM$theta

## [1] 0.6268215

resultEM$step

## [1] 19

따라서 EM 알고리즘을 통하여 구한 \(\theta\)의 MLE 값은 0.6268215이고 19의 반복을 통하여 찾았다는 것을 알 수 있다.

Newton 방법을 사용한 MLE 구하기

눈치 챈 독자들도 있겠지만, 위의 길동의 예제는 \(\theta\)의 MLE를 찾기 위해서 꼭 EM 알고리즘 만을 사용할 수 있는 것은 아니다. 위의 식 (1)의 로그우도 함수를 직접 뉴턴 방법을 사용해서 최대화를 할 수도 있기 때문이다.

뉴턴 방법을 적용하기 위해서는 식 (1)의 1차 도함수와 2차 도함수가 필요하므로 식 (1)을 \(\theta\)에 대하여 미분하도록 하자. 뉴턴 방법을 이용한 최대화에 대한 이해가 필요한 독자의 필자의 이전 포스팅을 참고하도록 하자.

\[ \frac{\partial l}{\partial\theta}=\frac{y_{1}}{2+\theta}-\frac{\left(y_{2}+y_{3}\right)}{1-\theta}+\frac{y_{4}}{\theta} \]

\[ \frac{\partial^{2}l}{\partial\theta^{2}}=-\frac{y_{1}}{\left(2+\theta\right)^{2}}-\frac{\left(y_{2}+y_{3}\right)}{\left(1-\theta\right)^{2}}-\frac{y_{4}}{\theta^{2}} \]

dloglik <- function(theta, y){
	y[1]/(2 + theta) - (y[2] + y[3])/(1 - theta) + y[4]/theta 
}
d2loglik <- function(theta, y){
	-y[1]/(2 + theta)^2 - (y[2] + y[3]) / (1 - theta)^2 - y[4] / theta^2 
}

길동의 데이터를 사용하여 1차 도함수와 2차 도함수를 정의하고, 이전에 필자가 정의해 놓은 함수를 사용해서 MLE를 찾아보면 다음과 같다. 결과값을 보면 EM 알고리즘을 통하여 찾은 \(\theta\)의 값과 동일하다는 것을 확인할 수 있다.

y <- c(125, 18, 20, 34)
dl <- function(theta){
	dloglik(theta, y)
}
d2l <- function(theta){
	d2loglik(theta, y)
}

library(r4issactoast)
newton_method(dl, d2l, 0.5, left_b = 0.25, right_b = 1)

## Final value of function: 7.105427e-15

## [1] 0.6268215

이 포스팅은 필자가 수강한 The University of Iowa의 Jian Huang 교수의 Advanced inference II 수업을 들으면서 필자가 했던 숙제를 정리한 것임을 밝혀둔다.

 

Reference

[1] Dempster, Arthur P., Nan M. Laird, and Donald B. Rubin. “Maximum likelihood from incomplete data via the EM algorithm.” Journal of the royal statistical society. Series B (methodological) (1977): 1-38.