저번 rejection sampling에 이어 이번 포스팅에서도 타겟 분포로부터 표본을 생성해낼 수 있는 샘플링 알고리즘 중 하나인 M-H Algorithm을 살펴보겠다.
2022.03.14 - [Statistics/Bayesian With Python] - 표본 샘플링 방법(1) - 기각 샘플링(Rejection Sampling)
표본 샘플링 방법(1) - 기각 샘플링(Rejection Sampling)
Bayesian과 Frequentist들의 가장 큰 차이점은 모수, parameter에 대한 관점의 차이이다. Frequentist들은 모수 θ를 Unknown constant로 가정하며, Bayesian은 모수 θ를 Random variable로 가정한다..
harang3418.tistory.com
● 메트로폴리스-헤이스팅스 알고리즘(Metropolis-Hastings Algorithm)
메트로폴리스-헤이스팅스 알고리즘(이하 M-H)은 MCMC 방법 중 하나로, 이전에 살펴봤던 기각 샘플링과 마찬가지로 표본을 직접적으로 샘플링 할 수 없는 타겟 분포 p(x)에서 표본을 얻을 수 있도록 하는 샘플링 방법이다.
이 방법은 Metropolis et al.,1953에서 제안된 Metropolis algorithm을 1970년 Hastings가 응용하면서 제안된 방법이며,
베이지안 접근법을 활용한 연구에서 모수 추론을 위해 가장 널리 사용되는 방법이기도 하다.
M-H 방법은 앞서 살펴본 Rejection sampling 방법과 제안 분포 q(x)에 대한 조건과 Acceptance rate α 계산에서 약간의 차이가 존재하지만, 전체적인 알고리즘의 의사결정과정은 매우 유사하므로 이해하는 것에 어렵지 않을 것으로 생각된다.
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)≤Mq(x)같은..), 제안 분포의 형태가 타겟 분포와 유사하도록 설정할 것을 권장한다. (→하지만, 이 부분은 타겟 분포의 차원 수가 커질 수록 만족시키기 어렵고, M-H 알고리즘 수렴에도 영향을 주게된다.)
M-H 방법은 기각 샘플링 방법과 마찬가지로 Acceptance rate α를 계산하고, U∼Unif(0,1)와 비교하여 candidate x∗를 표본으로 받아들일지를 결정하게 된다.
아래에서 의사결정과정을 차근차근 확인해보자.
Step 1. 초기 값 x(0)을 설정한다.
For iteration t=1,⋯,T :
Step 2. candidate x∗을 q(x|x(t−1))로부터 샘플링한다.
Step 3. Acceptance rate α=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 α는 p(x∗)p(x(t−1))로 계산된다.
이 때, α는 이미 타겟 분포의 표본으로 포함된 x(t−1)에 비해 candidate x∗가 타겟 분포에 포함될 확률이 얼마나 되는지 그 비율을 나타낸 것으로 이해하면 좋다.
때문에, 여기서의 Acceptance rate은 1보다 큰 값을 가질 수 있다.
여기서 발전된 M-H 방법 같은 경우, 제안 분포가 대칭 분포가 아니어도 적용이 가능하도록 응용되었기 때문에, 그에 따른 보정이 acceptance rate α에 포함되어 step 3에 나타난 α의 형태가 되었다고 이해하면 쉽다.
Step 4. α∗=min 그리고 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.
'Statistics > Bayesian With Python' 카테고리의 다른 글
표본 샘플링 방법(3) - 깁스 샘플링(Gibbs Sampling) (0) | 2023.06.01 |
---|---|
표본 샘플링 방법(1) - 기각 샘플링(Rejection Sampling) (0) | 2022.03.14 |
다양한 사전 분포 _ Prior distribution (0) | 2021.01.29 |
베이지안 모델링의 기본 (0) | 2021.01.13 |
댓글