Summary

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

출처: Lesson 5.6 MCMC example 2

R 의 내장된 데이터인 faithful 을 이용해서, 간헐천의 분화시간과 분화시간 사이의 대기시간을 이용해서, 군집화를 시도한다. 혼합모델의 파라미터 추정은 MCMC를 사용한다.

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)

MCMC 알고리즘 수행

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

library(mvtnorm)
library(ellipse)
library(MCMCpack)


rm(list=ls()) 
set.seed(2022)     # For reproducibility

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

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

## Initialize the parameters
# -----------------------------------------------------
# 각 성분에 대해서 동일한 가중치 1/3 씩을 할당한다.
# 분산이 너무 퍼져서, 느리게 수렴하는것을 막는다.
w <- rep(1,KK)/KK  #Assign equal weight to each component to start with

mu <- rmvnorm(KK, apply(x,2,mean), var(x))   #RandomCluster centers randomly spread over the support of the data

# 각 성분의 sigma 값은 데이터의 분산을 이용한다.
# 성분의 개수 KK 로 나누고 있는데, 너무 퍼져 있는 성분끼리에서 발생하는 문제를 막기 위해서이다.
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
# Sigma[3,,] = var(x)/KK

cc <- sample(1:KK, n, replace=TRUE, prob=w)

par(mfrow=c(1,1))
plot(x[,1], x[,2], 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)
}
title(main="Initial estimate + Observations")

MCMC 알고리즘을 수행해보자. 수렴여부를 확인하기 위해서 따로 번인을 수행하지는 않는다.

# Priors
aa <- rep(1, KK) # 사전분포로 가중치에 대해서는 동등하게 설정한다.

# 이변량 데이터의 평균을 사전분포의 평균으로 사용하는거야.
dd <- apply(x,2,mean) # 이변량 정규분포를 mu 들에 대한 사전분포로 필요하다.

# 비정보적 사전분포(flat 사전분포)를 사용하지 않는다는 의미이다.
DD <- 10*var(x) # 데이터의 분산으로 10배 스케일 되었다. 

nu <- p # 역-위샤트 사전분포는 j 파라미터를 가진다. 여기서는 2를 값으로 가진다.
SS <- var(x)/KK # scale matrix 를 가진다. 데이터의 분산을 3으로 나누고 있다. 경험적인 설정이댜. 

# 번인을 하지 않는데, 알고리즘이 어떻게 수렴하는지를 설명하기 위해서이다.
# Number of iteration of the sampler
rrr <- 40

# 샘플들의 값들을 저장한다.
cc.out <- array(0, dim=c(rrr, n))
w.out <- array(0, dim=c(rrr, KK))
mu.out <- array(0, dim=c(rrr, KK, p))
Sigma.out <- array(0, dim=c(rrr, KK, p, p))
logpost <- rep(0, rrr)


for(s in 1:rrr){
  # Sample the indicators
  for(i in 1:n){
    v = rep(0,KK)
    for(k in 1:KK){
      v[k] <- log(w[k]) + mvtnorm::dmvnorm(x[i,], mu[k,], Sigma[k,,], log=TRUE)  #Compute the log of the weights
    }
    v <- exp(v - max(v))/sum(exp(v - max(v)))
    
    # 각 관측치가 어느 성분에서 나왔는지를 말해준다.
    # 즉, 성분지시자이다.
    cc[i] <- sample(1:KK, 1, replace=TRUE, prob=v) 
  }
  
  # Sample the weights
  # 가중치 갱신, 이 강좌 전체적으로 혼합모델에서는 동일한 아래 코드가 사용된다.
  w <- as.vector(rdirichlet(1, aa + tabulate(cc)))
  
  # Sample the means
  # 평균을 샘플링한다.
  DD.st = matrix(0, nrow=p, ncol=p)
  for(k in 1:KK){
    mk <- sum(cc==k)
    xsumk <- apply(x[cc==k,], 2, sum)
    DD.st <- solve(mk*solve(Sigma[k,,]) + solve(DD)) # 역행렬 게산
    dd.st <- DD.st%*%(solve(Sigma[k,,])%*%xsumk + solve(DD)%*%dd) # 행렬곱셈, 역행렬
    mu[k,] <- as.vector(rmvnorm(1,dd.st,DD.st)) # 다변량 정규분포에서의 샘플링 등
  }
  
  # Sample the variances
  xcensumk <- array(0, dim=c(KK,p,p))
  for(i in 1:n){
    xcensumk[cc[i],,] <- xcensumk[cc[i],,] + (x[i,] - mu[cc[i],])%*%t(x[i,] - mu[cc[i],])
  }
  for(k in 1:KK){
    # 각 성분에 대해서 서로 다른 분산을 가지고
    # 역-위샤트 분포에서 샘플링 된다.
    Sigma[k,,] <- riwish(nu + sum(cc==k), SS + xcensumk[k,,])
  }
  
  # Store samples
  cc.out[s,] <- cc
  w.out[s,] <- w
  mu.out[s,,] <- mu
  Sigma.out[s,,,] <- Sigma
  for(i in 1:n){
    logpost[s] <- logpost[s] + log(w[cc[i]]) + mvtnorm::dmvnorm(x[i,], mu[cc[i],], Sigma[cc[i],,], log=TRUE)
  }
  logpost[s] = logpost[s] + ddirichlet(w, aa)
  for(k in 1:KK){
    logpost[s] <- logpost[s] + mvtnorm::dmvnorm(mu[k,], dd, DD, log=TRUE)
    logpost[s] <- logpost[s] + log(diwish(Sigma[k,,], nu, SS))
  }
  
  if(s/250==floor(s/250)){
    print(paste("s = ", s))
  }
  
}

## Plot the logposterior distribution for various samples
par(mfrow=c(1,1))
par(mar=c(4,4,1,1)+0.1)
plot(logpost, type="l", xlab="Iterations", ylab="Log posterior")

## Plot the density estimate for the last iteration of the MCMC
par(mfrow=c(1,1))
par(mar=c(4,4,2,1)+0.1)
plot(x[,1], x[,2], main=paste("s =",s,"   logpost =", round(logpost[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)
}