메트로폴리스-해스팅스(Metropolis-Hastings) 난수추출 알고리즘에 대하여

이번 포스팅에서는 베이지안 통계에서 많이 사용되는 알고리즘인 Metropolis-Hastings(이하 MH) 알고리즘에 대하여 간단히 알아보고자 한다.

 

MH 알고리즘은 마르코프 체인 몬테칼로(Markov Chain Monte Carlo) 방법 중 하나로, MH 알고리즘을 사용하여 타겟 분포를 정상분포(stationary distribution)로 갖는 마코브 체인을 발생시킬수 있기 때문이다. 이렇게 발생시킨 체인은 타켓 분포에서 발생시킨 상관성이 존재하는 표본으로 생각할 수 있는데, 이 표본들을 가지고 몬테칼로 방법을 적용하게 된다.

 

이러한 전반적인 과정을 머리속에 생각하면서, 오늘은 우리가 생각하는 특정 분포을 정상분포로 갖는 체인을 발생시키는 방법에 대하여 알아보자.

 

 

 

Metropolis-Hastings 알고리즘

MH 알고리즘은 분포 \(f\)를 정상분포로 갖는 마르코프 체인 \(\{X_0, X_1,...,X_t,...\}\)을 다음의 과정을 통하여 발생시킨다.

1. 우리가 잘 알고있으며, 표본을 쉽게 추출할 수 있는 분포 \(g(.|X_t)\)를 하나 정한다.
2. 마르코프 체인의 시작점 \(X_0\)를 정한다.
3. 체인이 수렴할 때까지 다음의 과정을 반복한다. 현재의 체인이 \(X_t=x_t\)일 경우,

1. \(X_{t+1}\)이 될 수 있는 후보값 Y를 분포 \(g(.|X_t=x_t)\)로부터 하나 추출한다.
2. 1)번에서 뽑힌 Y의 값, y,를 사용하여 다음의 값을 계산한다. \[ \frac{f\left(y\right)g\left(X_{t}=x_{t}|y\right)}{f\left(x_{t}\right)g\left(y|X_{t}=x_{t}\right)} \]
3. 균등분포에서 표본 하나, \(U\),를 뽑아서 다음과 같은 과정을 통하여 \(X_{t+1}\)을 설정한다. \[ X_{t+1}=\begin{cases} y & U\leq min\left(1,\frac{f\left(y\right)g\left(X_{t}=x_{t}|y\right)}{f\left(x_{t}\right)g\left(y|X_{t}=x_{t}\right)}\right)\\ x_{t} & Otherwise \end{cases} \]

3)의 과정은 마르코프 체인 \(X_t\)의 값이 \(x_t\)인 상황에서 \(X_{t+1}\)가 y가 될 확률이

\[ \alpha\left(x_{t},y\right):=min\left(1,\frac{f\left(y\right)g\left(X_{t}=x_{t}|y\right)}{f\left(x_{t}\right)g\left(y|X_{t}=x_{t}\right)}\right) \]

 

가 됨을 의미한다. 이것을 acceptance probability라고 부른다. 만약 \(\alpha\) 값이 1로 고정되어 있다면, 발생되는 체인은 우리가 잘 알고 있는 \(g(.|X_t)\)에 의하여 발생되는 체인이 아무런 제약을 받지않고 움직이는 상황이 될 것이다. 하지만 \(\alpha\) 값이 1보다 작다면, \(g(.|X_t)\)에 의하여 발생되는 체인에 일정한 제약이 걸린채로 움직이는 체인이 된다.

 

실제 예

Target density

MH 알고리즘을 적용하여, 실제 분포에서 표본을 추출하여 보도록 하자. 먼저 우리가 표본을 추출할 분포는 다음과 같다.

\[ f(x;\sigma )={\frac {x}{\sigma ^{2}}}e^{-x^{2}/(2\sigma ^{2})},\quad x\geq 0, \]

 

위의 분포는 Rayleigh 분포의 pdf로 \(\sigma>0\)는 스케일 모수의 역할을 한다.

 

Rayleigh 분포의 pdf - 위키피디아

후보값 제안분포 (Proposal distribution) 설정

알고리즘의 적용을 위하여 후보값을 제안하는 분포를 정해야 하는데, 보통 목표 분포와 같은 정의역을 갖는 분포를 사용한다. 따라서 우리는 0보다 큰 값에서 밀도함수값이 정의되는 감마분포를 사용하도록 하겠다.

\[ \begin{align*} Y & \sim Gamma\left(X_{t},1\right)\\ g\left(x_{t}|y\right) & =p\left(y|X_{t}=x_{t},\beta=1\right)\\ & =\frac{1}{\Gamma(x_{t})}y^{x_{t}-1}e^{-y} \end{align*} \]

채택 확률(Acceptance prob.) 구하기

앞에서 설명한 채택 확률을 구하기 위해서 위의 주어진 감마분포 밀도함수와 Rayleigh 분포의 밀도함수를 사용하여 다음의 값을 조금 정리하면 다음과 같다.

\[ \begin{align*} & \frac{f\left(y\right)g\left(X_{t}=x_{t}|y\right)}{f\left(x_{t}\right)g\left(y|X_{t}=x_{t}\right)}\\ = & \frac{\frac{y}{\sigma^{2}}e^{-y^{2}/(2\sigma^{2})}}{\frac{x_{t}}{\sigma^{2}}e^{-x_{t}^{2}/(2\sigma^{2})}}\times\frac{\frac{1}{\Gamma(y)}x_{t}^{y-1}e^{-x_{t}}}{\frac{1}{\Gamma(x_{t})}y^{x_{t}-1}e^{-y}}\\ = & \frac{\Gamma(x_{t})}{\Gamma(y)}\frac{y\left(x_{t}^{y-1}\right)exp\left(-\frac{y^{2}}{2\sigma^{2}}-x_{t}\right)}{x_{t}\left(y^{x_{t}-1}\right)exp\left(-\frac{x_{t}^{2}}{2\sigma^{2}}-y\right)} \end{align*} \]

 

이렇게 설정한 이유는 분모와 분자의 꼴이 정확히 대칭되는 꼴이므로 다음과 같이 R에서 함수를 쉽게 정의할 수 있기 때문이다.

deno <- function(fir, sec, sig) {
	gamma(fir) * sec * fir^(sec-1) * exp(-(sec^2)/(2*sig^2)-fir)
}
accept_p <- function(x_t, y, sigma) {
	min(1, (deno(x_t, y, sigma)/deno(y, x_t, sigma)))
}

R에서 알고리즘 구현하기

다음은 Rayleigh(\(\sigma=2\)) 분포를 정상분포로 갖는 마르코프 체인을 MH 알고리즘을 통해 발생시키는 과정을 R에서 구현한 것이다.

1. 초기값을 설정한다.

set.seed(2017) # Initialization
m <- 10000
mchain <- numeric(m)
my_sigma <- 2
rej_num <- 0

1. 먼저 시작점 \(X_0\)값을 감마분포에서 하나 뽑도록 하자.

mchain[1] <- rgamma(1, shape = 2)

1. 체인이 수렴할 때까지 다음의 과정을 반복한다.

for (i in 2:m) {
	xt <- mchain[i-1]
    y  <- rgamma(1, shape = xt, rate = 1)
    if (runif(1) <= accept_p(xt, y, my_sigma)){
    	mchain[i] <- y
        } else {
        mchain[i] <- xt
        rej_num <- rej_num + 1
    }
}
rej_num

## [1] 3989

head(mchain)

## [1] 3.770764 3.132473 3.132473 3.416740 2.631868 3.045583

rej_num에 비추어 보았을 때, 기각률이 40%정도 되는 것은 체인이 그닥 효율적이지 못하다는 것을 의미한다. 이상적인 체인은 25%정도의 기각률이라고 알려져있다(Roberts, Gelman and Gilks, 1997). 일단, 우리가 발생시킨 마르코프 체인을 그래프로 나타내보도록 하자.

library(ggplot2)
library(RColorBrewer)
mchain_part <- mchain[5000:6000]

p <- ggplot(data = data.frame(
					mc = mchain_part,
					time = c(5000:6000))) +
     geom_line(mapping = aes(y = mc, x = time)) + 
     theme_light()

p + labs(title = "Sample path of MC generated by MH algorithm",
		 y = "Markov chain")

위의 그래프에서 볼 수 있다시피 중간중간 체인이 걸려있는 모습을 보인다. 이러한 이유는 후보값이 계속 기각을 당해서 체인이 그 자리에 계속 있는 것이다.

표본 추출 결과확인

우리의 목표는 우리가 생성한 마르코프 체인의 값들이 목표 분포인 Rayleigh(\(\sigma=2\)) 분포를 따르게 만드는 것이었다. 이것을 눈으로 확인하기 위하여 QQ plot과 Rayleigh 분포의 밀도함수 값과 표본을 바탕으로 추정된 밀도함수의 겹쳐서 그려보자.

rayleigh_pdf <- function(x, sigma){
	(x / sigma^2) * exp(-x^2 / (2*sigma^2)) * (x > 0)
}

my_color <- brewer.pal(5, "Accent")
my_sample <- mchain[3000:10000] 

p <- ggplot(data = as.data.frame(my_sample)) + 
	geom_hline(yintercept=0) + 
	geom_vline(xintercept=0) +  
	theme_light() 
    
p <- p + 
	stat_density(mapping = aes(x = my_sample),      
							fill = my_color[1], 
                            alpha = 0.3) +   
	geom_rug(aes(x = my_sample, y = 0),
			 position = position_jitter(height = 0),
			 color = my_color[3]) +    
	stat_function(
    		fun = rayleigh_pdf,
    		args = list(sigma = 2),
			color = my_color[5],
            size = 0.7)

p  + labs(title = "Theoretical density vs. estimated density from MC",   
			  x = "Sample from Markov Chain",
		caption = "Fig. 1")

MCMC 결과 시각화

QQ plot은 이론적인 Quantile값과 우리가 추출한 표본의 Quantile값을 비교해보는 그래프로 목표 분포에 가까울수록 직선에 가까운 모양이 나타난다. Rayleigh 분포의 누적분포 함수는 다음과 같다.

\[ F(x;\sigma)={\displaystyle 1-e^{-x^{2}/\left(2\sigma ^{2}\right)}} \]

 

따라서, 0에서 1사이의 값, p,에 대응하는 Rayleigh 분포의 quantile 함수는 다음과 같다.

 

\[ {\displaystyle Q(p;\sigma )=\sigma {\sqrt {-2\ln(1-p)}}} \]

# Quantile function define 
q_rayleigh <- function(p, sigma){
	sigma * sqrt(-2 * log(1-p)) 
}

qq <- ggplot(data = as.data.frame(my_sample),
			 aes(sample = my_sample)) +
	geom_hline(yintercept=0) +
	geom_vline(xintercept=0) + 
	xlim(0, 10) +
    ylim(0, 10) +  
	coord_fixed(ratio = 1) +   
	theme_light() 

qq <- qq + 
	geom_abline(
		slope = 1,
    	intercept = 0, 
		color = my_color[1]) +    
	stat_qq(distribution = q_rayleigh,
			dparams = list(sigma = 2),
			size = 0.6)

qq + labs(title = "Q-Q plot of Rayleigh distribution",  
			x = "Theoretical Quantile", 
			y = "Sample Quantiles",   
			caption = "Fig. 2")

위의 QQ plot에서 볼 수 있듯이 MH를 통하여 발생시킨 체인이 목표 분포인 Rayleigh 분포를 잘 따르고 있는 것을 확인 할 수 있다. 하지만 마지막 꼬리부분에서는 샘플이 잘 뽑히지 않아서 약간 틀어지는 모습을 보인다.

마치며

이론적인 부분

이번 포스팅에서는 MH 알고리즘의 구현에만 초점을 맞췄다. MH를 통하여 발생시킨 체인이 이론적으로 목표 분포를 정상분포로 하는가에 대한 내용에 대하여는 추후 포스팅을 해보도록 하겠다. 또한 MH 알고리즘의 사용의 가장 중요한 포인트는 바로 제안분포(proposal dist.)를 어떻게 설정해야하는가이다. MH 알고리즘을 통하여 발생된 체인에 몬테카를로 방법을 적용하려면, irreducibility, positive recurrence, aperiodicity의 성질을 만족해야만 한다. 필자의 본문에 언급한 것처럼 목표 분포와 동일한 정의역을 가진 제안 분포를 사용하면, 대부분의 경우 발생시킨 체인이 이러한 조건들을 만족한다(Rizzo 2007, p248).

Rayleigh 분포 표본 추출

이번 포스팅에 사용된 Rayleigh에서 표본을 추출은 한가지 예일 뿐이다. 포스팅에서 구했든 quantitle 함수를 그냥 구할 수 있는 분포의 경우는 inverse transform을 이용하여 표본을 추출하는 것이 훨씬 효율적이다.

Reference

[1] Roberts, Gareth O., Andrew Gelman, and Walter R. Gilks. “Weak convergence and optimal scaling of random walk Metropolis algorithms.” The annals of applied probability 7.1 (1997): 110-120.

[2] Rizzo, Maria L. Statistical computing with R. CRC Press, 2007.