표본 샘플링 방법(2) - 메트로폴리스-헤이스팅스 알고리즘(Metropolis-Hastings Algorithm)

저번 rejection sampling에 이어 이번 포스팅에서도 타겟 분포로부터 표본을 생성해낼 수 있는 샘플링 알고리즘 중 하나인 M-H Algorithm을 살펴보겠다.

 

2022.03.14 - [Statistics/Bayesian With Python] - 표본 샘플링 방법(1) - 기각 샘플링(Rejection Sampling)

표본 샘플링 방법(1) - 기각 샘플링(Rejection Sampling)

---

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

 

 메트로폴리스-헤이스팅스 알고리즘(이하 M-H)은 MCMC 방법 중 하나로, 이전에 살펴봤던 기각 샘플링과 마찬가지로 표본을 직접적으로 샘플링 할 수 없는 타겟 분포 $p(x)$에서 표본을 얻을 수 있도록 하는 샘플링 방법이다.

 

이 방법은 Metropolis et al.,1953에서 제안된 Metropolis algorithm을 1970년 Hastings가 응용하면서 제안된 방법이며,

베이지안 접근법을 활용한 연구에서 모수 추론을 위해 가장 널리 사용되는 방법이기도 하다.

 

M-H 방법은 앞서 살펴본 Rejection sampling 방법과 제안 분포 $q(x)$에 대한 조건과 Acceptance rate $\alpha$ 계산에서 약간의 차이가 존재하지만, 전체적인 알고리즘의 의사결정과정은 매우 유사하므로 이해하는 것에 어렵지 않을 것으로 생각된다.

 

---

$p(x)$는 타겟 분포(target density)로 표본을 직접적으로 생성해 낼 수 없는 복잡한 형태를 가지고 있는 확률 분포일 때,

우리는 쉽게 표본을 생성할 수 있는 제안 분포(proposal density) $q(x)$를 도입하여 $p(x)$의 표본을 생성한다.

 

M-H 방법은 MCMC 방법 중 하나로 제안 분포 $q(x)$가 markov chain의 성질을 가지고 있다. 즉, $t$-시점의 제안 분포는 $t-1$-시점에서 샘플링된 $x^{(t-1)}$에 의존한다 : $q(x|x^{(t-1)})$

M-H 방법은 제안분포에 대해 추가적인 조건을 두고있진 않지만(가령 $p(x)\leq Mq(x)$같은..), 제안 분포의 형태가 타겟 분포와 유사하도록 설정할 것을 권장한다. (→하지만, 이 부분은 타겟 분포의 차원 수가 커질 수록 만족시키기 어렵고, M-H 알고리즘 수렴에도 영향을 주게된다.)

 

M-H 방법은 기각 샘플링 방법과 마찬가지로 Acceptance rate $\alpha$를 계산하고, $U\sim Unif(0,1)$와 비교하여 candidate $x^{*}$를 표본으로 받아들일지를 결정하게 된다.

아래에서 의사결정과정을 차근차근 확인해보자.

---

Step 1. 초기 값 $x^{(0)}$을 설정한다.

 

   For iteration $t=1,\cdots, T$ : 

         Step 2. candidate $x^{*}$을 $q(x|x^{(t-1)})$로부터 샘플링한다.

         Step 3. Acceptance rate $\alpha= \cfrac{p(x^{*})/q(x^{*}|x^{(t-1)})}{p(x^{(t-1)})/q(x^{(t-1)}|x^{*})}$ 를 계산한다.

메트로폴리스 알고리즘은 제안 분포 $q(x|x^{(t-1)})$가 대칭 분포임을 가정하고 있으며($q(x^{*}|x^{(t-1)}) = q(x^{(t-1)}|x^{*})$), Acceptance rate $\alpha$는 $\cfrac{p(x^{*})}{p(x^{(t-1)})}$로 계산된다.

이 때, $\alpha$는 이미 타겟 분포의 표본으로 포함된 $x^{(t-1)}$에 비해 candidate $x^{*}$가 타겟 분포에 포함될 확률이 얼마나 되는지 그 비율을 나타낸 것으로 이해하면 좋다.
때문에, 여기서의 Acceptance rate은 1보다 큰 값을 가질 수 있다.

여기서 발전된 M-H 방법 같은 경우, 제안 분포가 대칭 분포가 아니어도 적용이 가능하도록 응용되었기 때문에, 그에 따른 보정이 acceptance rate $\alpha$에 포함되어 step 3에 나타난 $\alpha$의 형태가 되었다고 이해하면 쉽다.

         

         Step 4. $\alpha^{*}=\min\{\alpha,1\}$ 그리고 $U\sim Unif(0,1)$를 생성한다.

         Step 5. 다음과 같은 조건으로 $x^{(t)}$를 update한다.

                   - IF $U\leq \alpha^{*}$ then, $x^{(t)}=x^{*}$

                   - IF $U>\alpha^{*}$ then, $x^{(t)} = x^{(t-1)}$

---

M-H에서는 candidate $x^{*}$가 표본으로 reject될 경우, accept이 될 때까지 새로운 candidate을 샘플링 하는 것이 아니라, 이전에 표본으로 받아들여진 $x^{(t-1)}$를 다시 한번 표본으로 선택하며, Algorithm이 종료되면 하나의 chain을 얻게된다 : $\boldsymbol{x}=\{x^{(1)},x^{(2)},\cdots,x^{(T)}\}$

 

이러한 의사결정과정 때문에 M-H 방법에서는 얻어진 chain의 수렴이 잘 이루어지는 것이 매우 중요한데, 적절한 제안 분포를 설정하지 않았을 경우 많은 iteration에서 candidate이 reject 되며, $x^{(t)}$가 같은 값을 갖게되는 경우가 발생할 수 있다. 즉, M-H의 수렴에 문제가 발생할 수 있다.

 

이러한 수렴의 문제를 극복하기 위해 MCMC 기반 알고리즘들에서 사용되는 몇 가지 존재하는데, 여기서는 burn-in period만 간단히 언급하려고 한다. (→ 나중에 수렴에 대한 내용을 따로 다루려고한다.)

 

burn-in period란 말 그대로 표본을 태워서 버려버리는 구간을 나타낸다.

Iteration $t=1$부터 일정한 시점 $b$까지는 수렴이 되는 과정이라고 여겨 그 사이에 샘플링된 표본들은 버리고, $b$시점 이후부터 샘플링된 $x^{(b+1)},\cdots,x^{(T)}$만을 $p(x)$의 표본으로 채택하는 것이다.

일반적으로 전체 iteration $T$가 클수록 수렴은 더 잘되는 것으로 알려져있고, burn-in period를 나타내는 시점 $b$는 전체 시점의 절반 $b=T/2$로 설정된다.

---

이해를 돕기 위해 예를 살펴보자. 

다음과 같은 타겟 분포$p(x)$를 가정하자.

 

$$p(x) \propto 0.3e^{-0.2x^2}+0.7e^{-0.2(x-10)^2}$$

이 때, 제안 분포 $q(x^{*}|x^{(t-1)})$를 $N(x^{(t-1)},\sigma^2)$으로 설정하고, $\sigma$의 값을 바꿔가며 M-H algorithm을 적용해보자.

 

$\sigma=\{1,5,10\}$ 총 3가지 경우에 대한 시뮬레이션을 진행하였고, 모든 경우에 대하여 전체 iteration number $T=5000$으로 설정되었다. 

 

시뮬레이션 결과, 제안 분포의 형태에 따라 근사와 수렴 결과에 차이가 나타나는 것을 알 수 있다.

베이지안에서 수렴을 진단하는 방법으로는 여러가지 metrics가 고안되어져왔지만, 일반적으로 metrics들은 보조적인 도구로 사용되며, visual inspection이 가장 중요하다.

 

이 예제에서는 주어진 iteration안에서 $\sigma=10$이 가장 좋은 근사결과와 수렴 정도를 보이며, 더 많은 iteration을 수행하면 더 좋은 수렴결과를 보여줄 것으로 기대된다.

 

본 포스팅에서 사용된 파이썬 코드는 더보기에서 확인 할 수 있다.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as stats

def p(x):
    pr=0.3*np.exp(-0.2*(x**2))+0.7*np.exp(-0.2*((x-10)**2))
    return pr
    
######## density plot of p(x) ########

grid_x=np.arange(-10,20,0.01)
plt.plot(grid_x,p(grid_x))
ax = plt.gca()
ax.axes.yaxis.set_visible(False)
plt.title('Density plot of p(x)')

######## Simulation 1 - when sigma=1 ########

T=5000 # Total iteration number
sigma=1
MH_x_1=np.random.normal(size=1,loc=0,scale=1) # Set initial value

for i in range(1,T+1):
    # Draw candidate
    candi_x = np.random.normal(size=1,loc=MH_x_1[i-1],scale=sigma)
    
    # Compute acceptance rate
    nume = p(candi_x)/stats.norm(MH_x_1[i-1], sigma).pdf(candi_x)
    deno = p(MH_x_1[i-1])/stats.norm(candi_x,sigma).pdf(MH_x_1[i-1]) 
    alpha = min(nume/deno,1)
    
    #Compare acceptance rate with U
    U = np.random.uniform(size=1)
    if U<=alpha: 
        MH_x_1=np.append(MH_x_1,candi_x)
    else:
        MH_x_1=np.append(MH_x_1,MH_x_1[i-1])

MH_x_1=np.delete(MH_x_1,0) # Drop initial value

######## Histogram of samples ########

grid_x=np.arange(-10,20,0.01)
plt.plot(grid_x,p(grid_x)/sum(p(grid_x))*100)
plt.hist(MH_x_1,bins=40,density=True)
ax = plt.gca()
ax.axes.yaxis.set_visible(False)
plt.title(r'Histogram of samples when $\sigma=1$')

---

• Metropolis, Nicholas, et al. "Equation of state calculations by fast computing machines." The journal of chemical physics 21.6 (1953): 1087-1092.
• Hastings, W. Keith. "Monte Carlo sampling methods using Markov chains and their applications." (1970): 97-109.