[MCMC] 메트로폴리스-헤이스팅스 알고리즘

앞선 게시글 (https://pasus.tistory.com/358)에서, 천이 커널(transition kernel) \(K\) 를 목표 확률밀도함수(pdf) \(p(\mathbf{x})\) 가 정상(stationary) 분포가 되도록 설계할 수 있다면, 마르코프 체인을 통해 목표 확률 분포에서 추출한 것과 동일한 샘플을 생성할 수 있음을 확인했다. 그렇다면, 우리가 원하는 목표 확률밀도함수에 대해 천이 커널을 어떻게 설계할 수 있을까?

메트로폴리스(Metropolis)가 처음 제안하고 이후 헤이스팅스(Hastings)가 일반화한 방법에 따르면, 이산(discrete) 상태 공간에서도 계산이 어려운 천이 커널의 명시적인 함수를 구하는 대신, 이전 샘플이 주어졌을 때 새로운 샘플을 추출하는 절차를 명확히 정의하는 방식으로 이를 해결한다.

메트로폴리스-헤이스팅스(MH, Metropolis-Hastings) 알고리즘은 가장 널리 사용되는 MCMC 방법으로서 한때 20세기 과학과 공학의 발전에 가장 큰 영향을 미친 10대 알고리즘으로 선정된 적이 있다고 한다.

우선 다음과 같이 목표 확률밀도함수가 주어졌다고 가정한다.

 

\[ \begin{align} p(\mathbf{x})= \frac{\tilde{p}(\mathbf{x})}{Z} \end{align} \]

 

여기서 \(Z\) 는 정규화 상수다. 메트로폴리스-헤이스팅스(MH) 알고리즘의 목표는 확률밀도함수 \(p(\mathbf{x})\) 로 부터 샘플 시퀀스 \(\mathbf{x}^{(t)}, \ t=1,...,T \) 를 생성하는 것이다. MH 알고리즘은 다음과 같다.

다음을 \(T\) 번 반복하며 매 시간스텝마다 새로운 샘플 \(\mathbf{x}^{(t+1)}\) 을 계산한다.

0. 첫 샘플 \(\mathbf{x}^{(1)}\) 을 임의로 선정한다.

1. 새로운 후보 샘플을 추출한다. \(\mathbf{x} \sim q(\mathbf{x} \vert \mathbf{x}^{(t)})\)
여기서 \(q(\mathbf{x} \vert \mathbf{x}^{(t)})\) 는 제안(proposal) 확률밀도함수로서, 이를 사용하여 현재 샘플 \(\mathbf{x}^{(t)}\) 를 조건으로 하는 후보 샘플 \(\mathbf{x}\) 를 생성한다.

2. 후보 샘플의 수락 확률(acceptance probability)을 계산한다.

 

\[ \begin{align} \alpha &= \min \left( 1, \ \frac{p (\mathbf{x}) \ q(\mathbf{x}^{(t)} \vert \mathbf{x}) }{ p( \mathbf{x}^{(t)}) \ q( \mathbf{x} \vert \mathbf{x}^{(t)} ) } \right) \\ \\ &= \min \left( 1, \ \frac{ \tilde{p}(\mathbf{x}) \ q(\mathbf{x}^{(t)} \vert \mathbf{x}) }{ \tilde{p}( \mathbf{x}^{(t)}) \ q( \mathbf{x} \vert \mathbf{x}^{(t)} ) } \right) \end{align} \]

 

3. 확률 \(\alpha\) 로 후보 샘플 \(\mathbf{x}\) 를 새로운 샘플로 수락하거나 거절한다.
즉 균등분포 \(\mathcal{U}[0, 1]\) 에서 샘플 \(u\) 를 추출한 후, \(u \sim \mathcal{U}[0, 1]\), 만약 \(u \lt \alpha\) 이면 수락하고 아니면 거절한다.

4. 수락하면 새로운 샘플로 편입하고 (\(\mathbf{x}^{(t+1)}= \mathbf{x}\)), 거절하면 현재의 샘플을 새로운 샘플로 유지한다 (\(\mathbf{x}^{(t+1)}= \mathbf{x}^{(t)}\)).

MH 알고리즘은 매우 간단하다는 장점이 있지만, 제안 확률밀도함수 \(q(\mathbf{x} \vert \mathbf{x}^{(t)})\) 의 선택에 따라 성능이 크게 달라진다. 예를 들어 제안 확률밀도함수의 표준편차가 너무 좁으면 \(p(\mathbf{x})\) 중에서 특정 모드 주변에서만 샘플링이 될 수도 있고 반면에 너무 넓으면 거절율이 매우 높아져 상관관계가 높아질 수 있다.

제안 확률밀도함수의 쉬운 예로서는 독립 샘플러(independent sampler)와 메트로폴리스 알고리즘을 들 수 있다. 독립 샘플러는 제안 확률밀도함수가 다음과 같이 현재의 샘플에 독립적인 것이다.

 

\[ \begin{align} q(\mathbf{x} \vert \mathbf{x}^{(t) } )=q(\mathbf{x}) \end{align} \]

 

메트로폴리스 알고리즘은 다음과 같이 제안 확률밀도함수가 대칭이라고 가정한다.

 

\[ \begin{align} q(\mathbf{x} \vert \mathbf{x}^{(t) } )=q(\mathbf{x}^{(t) } \vert \mathbf{x}) \end{align} \]

 

그러면 수락 확률은 다음과 같이 목표 확률밀도함수만의 함수가 된다.

 

\[ \begin{align} \alpha = \min \left( 1, \ \frac{ \tilde{p}(\mathbf{x}) }{ \tilde{p}( \mathbf{x}^{(t)}) } \right) \end{align} \]

 

이 밖에 제안 확률밀도함수의 설계 방식에 따라 다양한 MCMC 알고리즘이 고안되어 왔다.

한편 MH 알고리즘으로 생성된 샘플 중 초기 '번인(burn-in)' 기간의 샘플은 무시하고, 마르코프 시퀀스가 정상 상태에 도달한 이후의 샘플만 사용해야 한다.

 

 

다음은 MH 알고리즘의 시뮬레이션 예이다. 목표 확률밀도함수는 다음과 같다.

 

\[ \begin{align} p(x) \propto 0.5 \exp ⁡(-(x+2)^2)+0.5 \exp⁡(-(x-2)^2 ) \end{align} \]

 

 

그림에서 보듯이 두 개의 모드가 있는 확률밀도함수이다. 이런 확률밀도함수로 부터 샘플을 추출하는 것은 어렵다. 다음은 제안 확률밀도함수다. 샘플 추출이 용이하도록 다음과 같이 표준편차가 \(1\) 인 가우시안 확률밀도함수로 선택했다.

 

\[ \begin{align} q(x│x^{(t) } ) &= \mathcal{N}(x^{(t) },1) \\ \\ &= \frac{1}{ \sqrt{2 \pi} } \exp ⁡\left( - \frac{(x-x^{(t) } )^2}{2} \right) \end{align} \]

 

MH 알고리즘을 이용하여 30,000 개의 샘플 시퀀스를 생성하였다.

 

 

초기 번인(burn-in) 기간을 감안하여 5,000번째 이후의 샘플만을 이용하여 샘플의 히스토그램을 그려본 결과 목표 확률밀도함수와 거의 일치함을 확인할 수 있다.

 

 

다음은 예제 시뮬레이션에 사용된 매트랩 코드다.

 

%
% MCMC-MH algorithm
%
% (c) st.watermelon

clear;

% target pdf
p_pdf = @(x) 0.5*exp(-(x+2).^2) + 0.5 * exp(-(x-2).^2);

% Gaussian proposal pdf with mu-mean and 1-std
q_pdf = @(x, mu) 1/sqrt(2*pi)*exp(-(x-mu)^2/2);


% 0. initialize
n_step = 30000;
x = randn(1);

for t = 1 : n_step

    % 1. sampling from proposal q(x_new|x_t)
    x_new = randn(1) + x(t);

    % 2. acceptance probability

    alpha = min(1, ...
           (p_pdf(x_new)*q_pdf(x(t),x_new)) / (p_pdf(x(t))*q_pdf(x_new,x(t)))...
           );
    
    % 3. accept or reject
    u = rand(1);
    if u < alpha
        x = [x; x_new];
    else
        x = [x; x(t)];
    end
end


% result
xx = linspace(-5,5,1000);
figure
plot(x), title('MC'), xlabel('time step')
figure
h = histogram(x(5000:end),'BinWidth',0.2, 'Normalization','probability','FaceColor','y');
hold on
plot(xx, p_pdf(xx)/max(p_pdf(xx))*max(h.Values),'Color','b','linewidth',1.5)
grid on
figure
plot(xx, p_pdf(xx)/max(p_pdf(xx)));