메트로폴리스 헤이스팅스 알고리즘(Metropolis-Hastings algorithm)

Metropolis-Hastings(이하 MH) 알고리즘에 대해 알아볼 것이다.

 

MH 알고리즘은 MCMC(Markov Chain Monte-Carlo)의 일반적인 형태로써 특정 분포로부터 정상분포로 갖는 체인을 발생시킬 수 있다. 이를 이용하여 특정 분포로부터 데이터를 생성할 수 있다.

 

다룰 내용으로는 다음과 같다.

 

1. MH 알고리즘

2. Random walk MH 알고리즘 구현

 

1. MH 알고리즘

 

MH 알고리즘은 다음과 같다.

 

 

▷ q는 제안 분포(Proposal distribution)를 의미하고, g는 우리의 목적 분포(Target distribution)에서 정규화 상수(Normalizing constant)를 제외한 부분이다. 즉, 목적 분포와 g(theta)는 비례 관계가 성립한다.

 

▷ 초기값 theta 0를 정하고, q(theta i-1 | theta star)로부터 데이터를 생성한 뒤, alpha를 구한다. alpha가 1 이상이면 theta star를 theta i로 정하고, 아닌 경우에는 기존의 theta i-1를 theta i로 정한다.

 

▷ 이 과정을 통해 생성한 theta 중 수렴하기 시작한 부분부터 목적 분포로부터 생성한 데이터로써 활용할 수 있다.

 

MH 알고리즘의 핵심인 alpha가 어떻게 만들어졌는지 알아보자.

 

 

▷ 위의 그림은 두 노드간의 체인을 나타낸 것이다. 목적 분포로부터 데이터를 생성하기 위해서는 마르코프 체인(Markov chain)이 정상 상태(Stationary state)가 되어야 한다. 이를 위해  미세 균형 조건(Detailed balance condition)을 만족하도록 제약을 준다. 제약으로써 위의 식과 같이 각 상태의 전이확률(Transition probability)에 r을 곱하여 같도록 만든다.

 

▷ alpha의 분자 q(theta i | theta star)×g(theta star)가 분모 q(theta star | theta i)×g(theta i)보다 크면 분자에 곱해지는 r은 작게, 분모에 곱해진 r은 크게 만들어 줘야 미세 균형 조건이 성립하게 된다. 따라서 r theta i -> theta star은 1로 r theta star -> theta i는 alpha의 역수로 둔다. 반대의 경우에는 r theta i -> theta star은 alpha로 r theta i -> theta star는 1로 둔다. 이를 정리하면 다음과 같다.

 

 

▷ 이전 상태로부터 다음 상태의 예측을 통해 데이터를 생성하므로,  r theta i -> theta star를 통해 다음 상태를 결정하여야 한다. 즉, alpha에 따른 r theta i -> theta star가 생성된 데이터를 받아들일지 말지(Accept or not)에 대해 결정한다. 이를 정리하면 다음과 같다.

 

 

2. Random walk MH 알고리즘 구현

 

문제)

 

가능도 함수가 정규분포를 따르고, 정규분포의 평균은 t 분포를 따른다한다. 데이터로 1.2, 1.4, -0.5, 0.3, 0.9, 2.3, 1.0, 0.1, 1.3, 1.9가 주어졌을 때, 평균의 사후분포로부터 데이터를 생성하여라.

 

풀이)

 

 

▷ 정규화 상수를 모른다고 가정하고, 사전분포와 가능도 함수의 곱을 이용하여 사후분포로부터 데이터를 생성할 것이다. Random walk MH 알고리즘을 사용하여 데이터를 생성해보자. 아래의 코드는 R을 이용한 Random walk MH 알고리즘을 구현한 것이다.

 

In:

 

log_g <- function(mu, n, x_bar) {
  mu2 <- mu^2
  
  return(n * (x_bar * mu - mu2 / 2.0) - log(1.0 + mu2))
}

mh_sampl <- function(n_data, x_bar, n_iter, mu_init, prop_sd) {
  mu_out <- numeric(n_iter)
  accpt_cnt <- 0
  
  mu_now <- mu_init
  log_g_now <- log_g(mu_now, n_data, x_bar)
  
  for (i in 1:n_iter) {
    mu_cand <- rnorm(1, mu_now, prop_sd)
    log_cand <- log_g(mu_cand, n_data, x_bar)
    log_alpha <- log_cand - log_g_now
    alpha <- exp(log_alpha)
    
    if (runif(1) < alpha) {
      mu_now <- mu_cand
      accpt_cnt <- accpt_cnt + 1
      log_g_now <- log_cand
    }
    
    mu_out[i] <- mu_now
  }
  
  return(list(mu = mu_out, accpt_rate = accpt_cnt/n_iter))
}

 

▷ 구현 과정에서 alpha를 구하기 위해 약간의 트릭을 사용하였다. alpha를 구할 때, 로그를 취하여 계산한 뒤, 다시 지수를 통해 원래의 alpha로 바꿔주었다. 이는 계산과정 중 너무 작은 숫자로 인해 발생하는 반올림 오류(Rounding error)를 막기 위함이다.

 

▷ 제안 분포로 정규분포를 사용하였는데, 이는 alpha를 계산을 편리하게 해준다. 왜냐하면 alpha의 분자 및 분모의 제안 분포에 따른 확률이 같아져 상쇄되기 때문이다.

 

In:

 

x <- c(1.2, 1.4, -0.5, 0.3, 0.9, 2.3, 1.0, 0.1, 1.3, 1.9)

x_bar <- mean(x)
n_data <- length(x)

library(coda)

post_sampl <- mh_sampl(n_data = n_data, x_bar = x_bar, n_iter = 1e4, mu_init = 0.0, prop_sd = 1)

traceplot(as.mcmc(post_sampl$mu))

 

Out:

 

 

▷ 위의 그림은 각 반복 수행의 진행에 따른 생성된 mu의 변화를 나타낸 것이다. 특정 구간에서 동일한 mu로 유지하는 것을 확인할 수 있다. 이는 alpha가 1보다 작은 경우, alpha의 확률로 수락 또는 거절하게(Accept or reject) 되는데, 거절이 반복되어 나타난 결과이다.

 

In:

 

str(post_sampl)

 

Out:

 

List of 2
 $ mu        : num [1:1000] 0.54 0.54 0.54 0.54 0.54 ...
 $ accpt_rate: num 0.347

 

▷ 허용 비율이 34.7%인 것으로 나타났다. 일반적으로 허용 비율은 23 ~ 50%가 적당하다고 한다.

 

MH 알고리즘으로부터 생성된 데이터와 실제 사후분포와의 차이를 확인해보자.

 

In:

 

plot(density(post_sampl$mu, adjust = 2.5), main = '', xlim = c(-1.0, 3.0), ylim = c(0, 1.2), xlab = expression(mu))

curve(dt(x = x, df = 1), lty = 2, add = T) # prior for mu

curve(0.017 * exp(log_g(mu = x, n = n_data, x_bar = x_bar)), from = -1.0, to = 3.0, add = TRUE, col = 'blue')

 

Out:

 

 

▷ 파랑선은 실제 사후분포를 나타내고, 검정선은 생성한 데이터로부터 밀도 추정(Density estimation)을 통해 생성한 선이다. 점선은 사전분포를 나타낸다.

 

▷ 비교적 파랑선과 검정선이 비슷한 것을 확인할 수 있다. 생성 데이터의 갯수를 늘린다면 실제 사후분포에 더 가까운 분포가 나타날 것이다.

 

---

Reference:

"Bayesian Statistics: From Concept to Data AnalysisTechniques and Models," Coursera, https://www.coursera.org/learn/bayesian-statistics/.