1 Introduction

EM 알고리즘은 latent variable이 존재하는 statistical model의 maximum likelihood 혹은 maximum a posterior(MAP) 추정값를 구하기 위한 iterative method이다.

EM 알고리즘은 E step과 M step을 번갈아 적용한다.

Expectation (E) step에서는 parameter(\(\theta\))에 대한 현재 추정값(\(\theta^{(t)}\))를 이용해서 log likelihood의 기댓값을 계산하기 위한 함수(\(Q(\theta, \theta^{(t)}), Y\))를 만든다.

Maximization (M) step에서는 E step에서 계산된 expected log likelihood를 최대화하는 parameter 추정값들을 계산한다. 이 추정값들은 다음 E step에서 latent variable들의 분포를 결정하기 위해 사용된다.[1]

2 EM Algorithm

데이터 구조 아래와 같이 정의된다.

\(Y = (Y_1, \cdots, Y_n) \quad \leftarrow \quad observed \quad \overset{iid}{\sim} \quad g(y; \theta)\)

\(Z=(Z_1, \cdots, Z_n) \quad \leftarrow \quad unobserved\)

\(X=(Y, Z) \overset{iid}{\sim} f(x; \theta) \quad \leftarrow \quad complete\)

목표는 \(\hat{\theta}_{MLE}= \underset{\theta}{arg \;max} \; g(y;\theta)\)를 구하는 것이다.

2.1 E-step

log likelihood의 기대값 계산

\[Q(\theta, \theta^{(t)}, Y)=E_{\theta^{(t)}}[logf(X;\theta)|Y]\]

  • \(\theta^{(t)}\): current value

  • \(logf(X;\theta)\): complete data의 log likelihood

  • Z를 모르기 때문에 Joint likelihood(\(logf(X;\theta)\))를 관측된 데이터(\(Y\))로 conditional expectation 처리

2.2 M-step

이제 아래 식이 converge할 때 까지 돌리면 된다.

\[\theta^{(t+1)}= \underset{\theta}{arg\;max} Q(\theta, \theta^{(t)}, Y)\]

  • \(Q(\theta, \theta^{(t)}, Y)\)가 \(log(g(y; \theta))\)를 minorize함을 보이면 MM 알고리즘으로 설명가능

2.2.1 Jensen’s Inequality

http://people.duke.edu/~ccc14/sta-663-2016/14_ExpectationMaximization.html

http://people.duke.edu/~ccc14/sta-663-2016/14_ExpectationMaximization.html

For a convex function h(.) \[E(h(x)) \ge h(E(x))\]

  • Jensen’s Inequality를 이용하여 다음의 information inequality를 증명할 수 있다.

2.2.2 Information inequality

Let g, f are density, \[E_f (log(f)) \ge E_f (log(g))\]

[proof]

KL-divergence = \(E_f(-log \frac{g}{f})\)

\[E_f (-log \frac{g}{f}) \ge -logE_f (g/f) = -log \int \frac{g}{f} f dx= -log \int g dx=0\]

\[\underset{by \; Jensen's\; inequality}{\uparrow} \quad \quad \Leftrightarrow \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad\]

\[E_f(log(f)) \ge E_f(log(g))\]

2.2.3 Q가 g를 Minorize 함을 보임

\(log(g(y;\theta))\ge Q(\theta, \theta^{(t)}, Y)\)임을 보이면 Q가 g를 minorize함을 보이는 것이다.

\[< \;proof\; >\]

\[log(g(y;\theta))\ge Q(\theta, \theta^{(t)}, Y)\]

임을 보여야 하는데

Minorization의 정의에 의해

\[log(g(y;\theta^{(t)}))=Q(\theta^{(t)}, \theta^{(t)}, Y)\]임을 알고 있다.

이를 이용해 부등식을 쓰면

\[log(g(y;\theta))\ge Q(\theta, \theta^{(t)}, Y) + log(g(y;\theta^{(t)})) - Q(\theta^{(t)}, \theta^{(t)}, Y)\]

로 표현할 수 있고, 부등식을 정리하면 아래와 같다.

\[Q(\theta, \theta^{(t)}, Y)-log(g(y;\theta)) \le Q(\theta^{(t)}, \theta^{(t)}, Y)-log(g(y;\theta^{(t)}))\]

이제 위의 부등식이 성립함을 보이면 되는데 이는 information inequality로 간단하게 증명된다.

\[Q(\theta, \theta^{(t)}, Y)-log(g(y;\theta))=E_{\theta^{(t)}}(log\frac{f(y,z;\theta)}{g(y;\theta)}|Y)\]

가 되는데, Information inequality에 의해

\[Q(\theta, \theta^{(t)}, Y)-log(g(y;\theta)) =E_{\theta^{(t)}}(log\frac{f(y,z;\theta)}{g(y;\theta)}|Y) \le E_{\theta^{(t)}}(log\frac{f(y,z;\theta^{(t)})}{g(y;\theta^{(t)})}|Y)\\ =Q(\theta^{(t)}, \theta^{(t)}, Y)-log(g(y;\theta^{(t)}))=0\]

가 성립한다.

3 Example 1

3.1 Gaussian Mixture

\[y_1, \cdots, y_n \sim pf_1+(1-p)f_2, \quad \quad f_1 \sim N(\mu_1, \sigma_1^2), f_2 \sim N(\mu_2, \sigma_2^2)\]

\[\theta = (\mu_1, \mu_2, \sigma_1^2, \sigma_2^2, p)\]

\(\hat{\theta}_{MLE}\)를 EM 알고리즘을 통해 구해보자.

EM 알고리즘을 적용하기 위해 missing(unobserved)을 만들어내면

\[y_i = z_i x_{i1}+(1-z_i)x_{i2}, \quad where \;\; z_i \sim Bern(p), \;\; x_{i1} \sim N(\mu_1, \sigma_1^2), \;\; x_{i2} \sim N(\mu_2, \sigma_2^2)\]

1) E-step : complete likelihood 계산

\[f(x;\theta)=f_1(y_i ; \mu_1, \sigma_1^2)\times p, \quad \quad if \; z_i=1 \\ \quad \quad \quad \quad \quad= f_2(y_i ; \mu_2, \sigma_2^2)\times (1-p), \quad \quad if \; z_i=0\]

이를 다시 쓰면

\[f(x;\theta)=[f_1(y_i ; \mu_1, \sigma_1^2)\times p]^{z_i} [f_2(y_i ; \mu_2, \sigma_2^2)\times (1-p)]^{1-z_i}\]

양변에 로그를 취하면

\[logf(x;\theta)=z_ilogf_1(y_i ; \mu_1, \sigma_1^2) + z_ilogp + (1-z_i)logf_2(y_i ; \mu_2, \sigma_2^2) + (1-z_i)log(1-p) \]

이제 conditional expectation을 취해야 하는데, expectation은 \(z_i\)에만 걸린다.

\[E_{\theta^{(t)}}[logf(x;\theta)|Y=y_i]=E_{\theta^{(t)}}[z_i|Y=y_i]logf_1(y_i ; \mu_1, \sigma_1^2) + E_{\theta^{(t)}}[z_i|Y=y_i]logp \\ \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad -E_{\theta^{(t)}}[z_i|Y=y_i]logf_2(y_i ; \mu_2, \sigma_2^2) -E_{\theta^{(t)}}[z_i|Y=y_i]log(1-p) \\ \quad + logf_2(y_i ; \mu_2, \sigma_2^2) + log(1-p) \\ \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \cdots \quad (1)\]

이제 \(E_{\theta^{(t)}}[z_i|Y=y_i]\)을 구해서 대입해주면 된다. \(E_{\theta^{(t)}}[z_i|Y=y_i]\)는 Bayes rule을 이용해서 구할 수 있다.

\[E_{\theta^{(t)}}[z_i|Y=y_i] = P(z_i=1|Y=y_i) = w_i^{(t)} \\ = \frac{P(Y=y_i|z_i=1)P(z_i=1)}{P(Y=y_i|z_i=1)P(z_i=1) + P(Y=y_i|z_i=0)P(z_i=0)} \\ =\frac{f_1(y_i; \mu_1^{(t)}, \sigma_1^{2(t)})p^{(t)}}{f_1(y_i; \mu_1^{(t)}, \sigma_1^{2(t)})p^{(t)} + f_2(y_i; \mu_2^{(t)}, \sigma_2^{2(t)})(1-p^{(t)})}\]

\(w_i^{(t)}\)를 (1)에 대입하면

\[Q(\theta, \theta^{(t)}, Y)=E_{\theta^{(t)}}[logf(x;\theta)|Y=y_i] \\ = \sum_{i=1}^n[{w_i^{(t)} logf_1(y_i;\mu_1,\sigma_1^2)} + (1-w_i^{(t)})logf_2(y_i;\mu_2,\sigma_2^2)+w_i^{(t)}p+(1-w_i^{(t)})(1-p)]\]

2) M-step

\(Q(\theta, \theta^{(t)}, Y)\)를 각 모수로 편미분하여 0으로 놓고 풀면 다음의 5개 식을 얻을 수 있다.

\[\mu_1^{(t+1)}=\frac{\sum_{i=1}^nw_i^{(t)}y_i}{\sum_{i=1}^nw_i^{(t)}}\]

\[\mu_2^{(t+1)}=\frac{\sum_{i=1}^n(1-w_i^{(t)})y_i}{\sum_{i=1}^n(1-w_i^{(t)})}\]

\[\sigma_1^{2(t+1)} = \frac{\sum_{i=1}^nw_i^{(t)}(y_i-\mu_1^{(t)})^2}{\sum_{i=1}^n \mu_1^{(t)}}\]

\[\sigma_2^{2(t+1)} = \frac{\sum_{i=1}^n(1-w_i^{(t)})(y_i-\mu_2^{(t)})^2}{\sum_{i=1}^n \mu_2^{(t)}}\]

\[p^{(t+1)}=\frac{1}{n}\sum_{i=1}^nw_i^{(t)}\]

3.2 R code

set.seed(3)
n <- 500 # sample size

p <- 0.5 # true, 0.3
m1 <- -3 # mean
s1 <- 1 # sd

m2 <- 3 # mean
s2 <- 2 # sd
true.theta <- c(m1, m2, s1, s2, p) # 5차원
# data generation
z <- rbinom(n, 1, p)
x1 <- rnorm(n, m1, s1)
x2 <- rnorm(n, m2, s2)

y <- z * x1 + (1-z) * x2
plot(density(y))

max.iter <- 100
eps <- 1.0e-5

# initialization
m1 <- 1
m2 <- 2
s1 <- 1
s2 <- 2
p <- 0.3

theta <- c(m1, m2, s1, s2, p)
print(theta)
[1] 1.0 2.0 1.0 2.0 0.3
for (iter in 1:max.iter)
{
  p1 <- p
  p2 <- 1-p1

  # compute weight
  f1 <- exp(log(p1) + dnorm(y, m1, s1, log = T)) 
  # log transformation => underflow 막기 위해서, => 작은 값들 보존
  f2 <- exp(log(p2) + dnorm(y, m2, s2, log = T))

  w1 <- f1/(f1 + f2)
  w2 <- f2/(f1 + f2) # 1- w1

  # update mixture proportion
  p <- mean(w1)

  # update mean
  m1 <- sum(w1 * y)/sum(w1)
  m2 <- sum(w2 * y)/sum(w2)

  # update standard deviation
  s1 <- sqrt(sum(w1 * (y - m1)^2)/sum(w1))
  s2 <- sqrt(sum(w2 * (y - m2)^2)/sum(w2))

  new.theta <- c(m1, m2, s1, s2, p)

  if (max(abs(new.theta - theta)) < eps) break
  theta <- new.theta
}

print(cbind(true.theta, new.theta)) # 어떤게 mu1, mu2인지 지정되지 않아서 순서가 섞여있음
     true.theta  new.theta
[1,]       -3.0  3.0379737
[2,]        3.0 -3.0498538
[3,]        1.0  1.9862645
[4,]        2.0  0.9882122
[5,]        0.5  0.4872378