2개의 Gaussian mixture에 대한 EM algorithm

R · 케이플러스 한성탁 · 05 Jan 2017

EM 알고리즘은 관측되지 않는 잠재변수에 의존하는 확률 모델에서 최대가능도(maximum likelihood)나 최대사후확률(maximum a posteriori)을 갖는 매개변수를 찾는 반복적인 알고리즘입니다. (출처 : 위키백과)

EM 알고리즘에 대한 간단한 예제로 2개의 정규 분포를 일정 확률로 섞은 Gaussian mixture의 모수를 maximum likelihood로 추정해 보는 과정을 살펴보겠습니다. 아래 코드는 R. Tibshirani의 The Elements of Statistical Learning의 EM 알고리즘 설명을 R로 구현한 것입니다.

초기 설정

#최대 반복
set.seed(0)
maxiteration = 100
l0 = vector(length=maxiteration)

반복문의 최대 실행 횟수를 설정하고, log-likelihood 값을 저장할 벡터를 생성합니다.

# input data
data0=c(-0.39, 0.12, 0.94, 1.67, 1.76, 2.44, 3.72, 4.28, 4.92, 5.53,
        0.06, 0.48, 1.01, 1.68, 1.80, 3.25, 4.12, 4.60, 5.28, 6.22)
N=length(data0)

데이터를 입력하고, 데이터의 크기 \(N\)을 지정하였습니다.

1. 모수 초기화

첫 번째로 추정할 모수 \(\hat{\theta}^{(0)}\)를 초기화하는 단계입니다.

우선 \(\hat{\pi} = 0.5\)로 설정하고, \(\pi\)에 대한 베르누이 난수 \(Z_i\)를 발생시켜 원본 데이터를 두 그룹으로 나눕니다. 그리고 각각의 데이터 그룹에 대해 평균과 분산을 추정하고 이를 초기값으로 설정합니다.

2, 3. Expectation, Maximization

for(k in 1:maxiteration){
  
  # 2. Expectation Step : compute the responsivilities
  gamma.i = pi0*dnorm(data0, mu2, sigma2) / 
              ((1-pi0)*dnorm(data0, mu1, sigma1) + pi0*dnorm(data0, mu2, sigma2))
  
  # 3. Maximization Step : compute the weighted means and variances
  mu1 = sum((1-gamma.i)*data0)/sum(1-gamma.i)
  sigma1= sqrt( sum( (1-gamma.i)*(data0-mu1)^2 ) / sum(1-gamma.i) )
  mu2 = sum(gamma.i*data0)/sum(gamma.i)
  sigma2 = sqrt(sum(gamma.i*(data0-mu2)^2)/sum(gamma.i))
  
  pi0 = sum(gamma.i)/N
  l0[k] = sum(log( (1-pi0)*dnorm(data0, mu1, sigma1) + pi0*dnorm(data0, mu2, sigma2) ))
  
# 4. Iterate steps 2 and 3 until convergence
  #cat("(", k ,")", "pi :", pi0, "\n")
  
  if(k>1) {
    if(l0[k] - l0[k-1]<1e-8) {l0=l0[1:k]; break}
  }
}

이제 step 2에서 \(Q(\theta, \hat{\theta}^{(j)}) = E[l_0(\theta；Y)|X, \hat{\theta}^{(j)}] = \hat{\gamma}_j\)를 구합니다.

\[\hat{\gamma}_{i}=\frac{\hat{\pi}\phi_{\hat{\theta}_2}(y_i)}{(1-\hat{\pi})\phi_{\hat{\theta}_1}(y_i)+\hat{\pi}\phi_{\hat{\theta}_2}(y_i)}, i = 1, 2, ... , N\] 을 통해 계산할 수 있습니다.

또한 step 3에서는 앞에서 구한 \(\hat{\gamma}_i\)를 이용하여 \(\hat{\mu}_1, \hat{\mu}_2, \hat{\sigma}_1, \hat{\sigma}_2\)의 새로운 값을 구할 수 있습니다.

\[\begin{matrix} \hat{\mu}_1 = \frac{\sum_{i=1}^{N} (1-\hat{\gamma}_i)y_i}{\sum_{i=1}^{N} (1-\hat{\gamma}_i)}, & \hat{\sigma}_1 = \sqrt\frac{\sum_{i=1}^{N} (1-\hat{\gamma}_i)(y_i-\hat{\mu}_1)^2}{\sum_{i=1}^{N} (1-\hat{\gamma}_i)}, \\ \hat{\mu}_2 = \frac{\sum_{i=1}^{N} \hat{\gamma}_i y_i}{\sum_{i=1}^{N} \hat{\gamma}_i}, & \hat{\sigma}_2 = \sqrt\frac{\sum_{i=1}^{N} \hat{\gamma}_i(y_i-\hat{\mu}_2)^2}{\sum_{i=1}^{N} \hat{\gamma}_i}, \end{matrix}\]

한편 \(\hat{\pi} = \sum_{i=1}^{N}\hat{\gamma}_i/N\) 입니다.

4. 반복

위의 step 2, 3을 최대 반복수까지 반복합니다. 이 때 log-likelihood 값의 증가량이 \(10^{-8}\)보다 작으면 반복을 그만하도록 했습니다. log-likelihood 값은 대체로 100번 내로 수렴합니다.

plot(l0, type="b")

result = c(mu1, mu2, sigma1^2, sigma2^2, pi0)
names(result) = c("mu1", "mu2", "sigma1^2", "sigma2^2", "pi0")

result

##       mu1       mu2  sigma1^2  sigma2^2       pi0 
## 4.6558926 1.0831443 0.8188230 0.8113422 0.5545850

위의 결과는 예제 데이터에 대해 EM 알고리즘을 적용한 결과입니다. 책의 결과와 비슷한 것을 볼 수 있습니다.

원문에서 계산된 예측치는

\[\begin{matrix} \hat{\mu}_{1} = 4.62, && \hat{\sigma}_{1}^{2} = 0.87, \\ \hat{\mu}_2 = 1.06, && \hat{\sigma}_{1}^{2} = 0.77, \\ & \hat{\pi} = 0.546 \end{matrix}\]

이었습니다.