Summary

본 내용은 코세라 BAYES3 과정의 Lesson 4.5 EM example 2 의 코드를 사용했습니다. 자세한 내용은 해당 강의를 참조하세요.

출처: Lesson 4.5 EM example 2

R 의 내장된 데이터인 faithful 을 이용해서, 간헐천의 분화시간과 분화시간 사이의 대기시간을 이용해서, 군집화를 시도한다. 군집화 알고리즘으로 EM(Expectation Maximization) 을 사용한다.

faithful 데이터는, 미국 와이오밍주 옐로스톤 국립공원에 있는 Old faithful 간헐천(geyser) 에서의 분화지속시간(eruptions, min) 과 분화지속시간 사이의 대기사간(waiting, min)으로 측정한 자료이다.

크게 2개의 군집이 있을것으로 예상한다. 분화지속시간이 길수록 대기시간이 길거나, 분화지속시간이 짧을수록 대기시간이 짧거나 하는 부류이다.

기본 데이터 확인

# 데이터 구조 확인
str(faithful)
## 'data.frame':    272 obs. of  2 variables:
##  $ eruptions: num  3.6 1.8 3.33 2.28 4.53 ...
##  $ waiting  : num  79 54 74 62 85 55 88 85 51 85 ...
# 데이터 앞 부분 확인
head(faithful)
##   eruptions waiting
## 1     3.600      79
## 2     1.800      54
## 3     3.333      74
## 4     2.283      62
## 5     4.533      85
## 6     2.883      55
# 분화시간 대비 대기시간에 대한 플롯
plot(faithful)

EM 알고리즘 수행

이변량 가우시안 성분을 가진 혼합모델(Mixture model) 을 피팅하기 위해 EM 알고리즘을 사용한다. 초기 그래프는 아래와 같다.

library(mvtnorm)    # Multivariate normals are not default in R
library(ellipse)    # Required for plotting
## 
## 다음의 패키지를 부착합니다: 'ellipse'
## The following object is masked from 'package:graphics':
## 
##     pairs
rm(list=ls()) 
set.seed(2022)     # For reproducibility

x <- data.matrix(faithful) # XXX: 데이터를 행렬형태로 변환한다.
n <- nrow(faithful)

KK <- 2 # 혼합모델에서 2개의 성분으로 작업한다.
p <- 2 # 이변량 정규분포, 2는 bivariate를 의미함 

## EM 알고리즘의 파라미터 초기값 지정
# 각 성분에 대한 가중치를 1/2로 설정한다.
w <- rep(1,KK)/KK  

# 평균은 무작위로 2개를 생성한다.
# apply(x, 2, mean) 은 각 세로축(컬럼)에 대해서 평균을 구한다.
# var(x) = cov(x)
mu <- rmvnorm(KK, apply(x, 2, mean), var(x))

# 분산-공분산 행렬에 대해서 초기값 설정
Sigma <- array(0, dim=c(KK,p,p))  #Initial variances are assumed to be the same
for (k in 1:KK) {
  Sigma[k,,] <- var(x)/KK  
}
# Sigma[1,,] = var(x)/KK  
# Sigma[2,,] = var(x)/KK

# 초기 설정값에 따른 그래프를 플롯함
par(mfrow=c(1,1))
plot(x[,1], x[,2], xlab=expression(x[1]), ylab=expression(x[2]), 
     xlim = c(0, 7), ylim = c(30, 100))
for(k in 1:KK){
  lines(ellipse(x=Sigma[k,,], centre=mu[k,], level=0.50), col=k, lty=2, lwd=2)
  lines(ellipse(x=Sigma[k,,], centre=mu[k,], level=0.82), col=k, lty=2, lwd=2)
  lines(ellipse(x=Sigma[k,,], centre=mu[k,], level=0.95), col=k, lty=2, lwd=2)
}
title(main="Initial estimate + Observations")

EM 알고리즘을 수행한 후, 수렴된 Mixture 모델의 그래프를 확인해보자.

s <- 0
sw <- FALSE # 알고리즘 중지 변수, 솔루션에 충분히 가까워진다면 알고리즘을 중지한다.
QQ <- -Inf # XXX: 강의 영상에서는 KL 로 소개하고 있음
QQ.out <- NULL # XXX: 강의 영상에서는 KL.out 으로 소개하고 있음
epsilon <- 10^(-5)

while(!sw){
  ## E step
  v <- array(0, dim=c(n,KK))
  for(k in 1:KK){
    # 다변량 정규분포 dmvnorm 으로 사용한다.
    v[,k] <- log(w[k]) + dmvnorm(x, mu[k,], Sigma[k,,],log=TRUE)  #Compute the log of the weights
  }
  for(i in 1:n){
    v[i,] <- exp(v[i,] - max(v[i,]))/sum(exp(v[i,] - max(v[i,])))  #Go from logs to actual weights in a numerically stable manner
  }
  
  ## M step
  w <- apply(v,2,mean) # omega 추정치 계산
  
  # mu 추정치 계산
  # XXX: 참고로 강의에서는 동일 루프 안에 mu, sigma 를 같이 계산하고 있다.
  mu <- array(0, dim=c(KK, p)) 
  # 각각의 성분에 대해서
  for(k in 1:KK){
    # 각각의 관측치에 대해서
    for(i in 1:n){
      # mu 가 벡터로 처리되고 있다.
      mu[k,] <- mu[k,] + v[i,k]*x[i,]
    }
    mu[k,] <- mu[k,]/sum(v[,k])
  }
  
  # 분산 추정치 계산
  Sigma <- array(0, dim=c(KK, p, p))
  for(k in 1:KK){
    for(i in 1:n){
      Sigma[k,,] <- Sigma[k,,] + v[i,k]*(x[i,] - mu[k,])%*%t(x[i,] - mu[k,])
    }
    Sigma[k,,] <- Sigma[k,,]/sum(v[,k])
  }
  
  ##Check convergence
  # 수렴되었는지를 체크하기 위한, Q함수값 계산
  QQn = 0
  for(i in 1:n){
    for(k in 1:KK){
      QQn <- QQn + v[i,k]*(log(w[k]) + dmvnorm(x[i,],mu[k,],Sigma[k,,],log=TRUE))
    }
  }
  if(abs(QQn-QQ)/abs(QQn)<epsilon){
    sw <- TRUE
  }
  QQ <- QQn 
  QQ.out <- c(QQ.out, QQ)
  s <- s + 1
  
  # print(paste(s, QQn))
  
  #Plot current components over data
  # 현재 추정된 파라미터 기초로 그래프를 플롯한다.
  # layout(matrix(c(1,2),2,1), widths=c(1,1), heights=c(1.3,3))
  # par(mar=c(3.1,4.1,0.5,0.5))
  # plot(QQ.out[1:s],type="l", xlim=c(1,max(10,s)), las=1, ylab="Q")
  # 
  # par(mar=c(5,4,1,0.5))
  # plot(x[,1], x[,2], main=paste("s =",s,"   Q =", round(QQ.out[s],4)),
  #      xlab=expression(x[1]), ylab=expression(x[2]), lwd=2)
  # for(k in 1:KK){
  #   lines(ellipse(x=Sigma[k,,], centre=mu[k,], level=0.50), col=k, lty=2, lwd=2)
  #   lines(ellipse(x=Sigma[k,,], centre=mu[k,], level=0.82), col=k, lty=2, lwd=2)
  #   lines(ellipse(x=Sigma[k,,], centre=mu[k,], level=0.95), col=k, lty=2, lwd=2)
  # }
  
}

#Plot current components over data
layout(matrix(c(1,2),2,1), widths=c(1,1), heights=c(1.3,3))
par(mar=c(3.1,4.1,0.5,0.5))
plot(QQ.out[1:s],type="l", xlim=c(1,max(10,s)), las=1, ylab="Q", lwd=2)

par(mar=c(5,4,1,0.5))
plot(x[,1], x[,2], main=paste("s =",s,"   Q =", round(QQ.out[s],4)), xlab=expression(x[1]), ylab=expression(x[2]))
for(k in 1:KK){
  lines(ellipse(x=Sigma[k,,], centre=mu[k,], level=0.50), col=k, lty=2, lwd=2)
  lines(ellipse(x=Sigma[k,,], centre=mu[k,], level=0.82), col=k, lty=2, lwd=2)
  lines(ellipse(x=Sigma[k,,], centre=mu[k,], level=0.95), col=k, lty=2, lwd=2)
}