The Metropolis-Hastings algorithm

이전 포스팅인 Rejection sampling 이나 Sampling-Importance-Resampling algorithm과 같이 특정 확률분포 함수로부터 샘플링을 추출해 확률값을 근사시킬 수 있는 알고리즘이다. 다만 마르코프 체인의 개념을 이용해 이전확률값을 기준으로 현재 확률값을 평가해 이 기준에 만족하는 경우 샘플을 수용하고 아닐 경우 과거 샘플을 다시 넣는 방식을 사용한다. 흡사 Rejection Sampling과 개념은 비슷하나 Rejection 샘플링의 경우 샘플을 버리는 방식을 사용한다는 차이점이 있다. 이 알고리즘의 경우 다차원 분포의 샘플링을 위해 주로 사용되며 많은 MCMC 패키지들이 이 알고리즘을 포함하고 있고 실제 많이 사용된다.

마르코프 체인은 현재 상태가 과거 상태의 영향을 받는다는 것을 의미하며, 아래와 같은 코드와 그래프가 이전 상태의 영향을 받는 마르코프 체인의 한 예를 보여준다.

library(ggplot2)
set.seed(1234)
N <- 10000
x <- vector(N, mode = "numeric")
x[1] <- rnorm(1, mean = 0, sd = 2)
for (i in 1:(N - 1)) {
    x[i + 1] <- x[i] + rnorm(1, mean = 0, sd = 2)
}

ggplot(data.frame(x = 1:N, y = x), aes(x, y)) + geom_line()

일단 random wark를 만족하는 기본적인 알고리즘을 구현하기 위해서는 대칭 분포를 가지는(q(x'|x_t) == q(x_t|x')) 임의의 확률 분포 함수를 제공해야 되는데 주로 가우시언분포 함수를 사용한다.

아래 코드는 mu가 0이고 sigma가 5인 임의 확률 분포 함수를 기반으로 샘플링을 하고 이 값이 이전 p()값을 기준으로 alpha값이 유의미하게 크면 샘플링을 수용하고 아니면 과거 샘플링값을 재사용 하는 프로세스를 수행한다.

아래 코드에서 mu와 sigma를 조절해 가면서 시뮬레이션 해보면 재미난 결과가 나온다.

sigma가 너무 작으면 random walk가 제대로 동작하지 않아 샘플링이 제대로 수행되지 않으니 테스트 해보길 바란다.

역시 책을 참고했다.

p <- function(x) {
    mu1 <- 3
    mu2 <- 10
    v1 <- 10
    v2 <- 3
    0.3 * exp(-(x - mu1)^2/v1) + 0.7 * exp(-(x - mu2)^2/v2)
}


mh <- function(N) {
    u <- runif(N)
    mu <- 0
    sigma <- 5
    y <- vector(N, mode = "numeric")
    y[1] <- rnorm(1, mean = mu, sd = sigma)
    for (i in 1:N) {
        # 마르코프 체인
        ynew <- y[i] + rnorm(1, mean = mu, sd = sigma)
        alpha <- min(1, p(ynew)/p(y[i]))
        if (u[i] < alpha) {
            y[i + 1] <- ynew
        } else {
            y[i + 1] <- y[i]
        }
    }
    return(y)
}

그래프

N <- 5000
real <- data.frame(real_y = p(seq(-10, 20, by = 0.8)), real_x = seq(-10, 20, 
    by = 0.8))
real$count <- with(real, real_y * N/sum(real_y))
mh_sampling <- data.frame(mh = mh(N))

ggplot(mh_sampling, aes(mh)) + geom_histogram() + geom_line(aes(x = real_x, 
    y = count), data = real, colour = "red2") + scale_fill_gradient2(guide = guide_legend(reverse = TRUE))

The Metropolis-Hastings algorithm by from __future__ import dream is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.