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

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

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

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

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

 

$$\begin{array}{r}p(\mathbf{x})=\frac{\tilde{p}(\mathbf{x})}{Z}\end{array}$$

 

여기서 $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}|{\mathbf{x}}^{(t)})$
여기서 $q(\mathbf{x}|{\mathbf{x}}^{(t)})$ 는 제안(proposal) 확률밀도함수로서, 이를 사용하여 현재 샘플 ${\mathbf{x}}^{(t)}$ 를 조건으로 하는 후보 샘플 $\mathbf{x}$ 를 생성한다.

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

 

$$\begin{array}{rl}\alpha & =min(1,\frac{p(\mathbf{x})q({\mathbf{x}}^{(t)}|\mathbf{x})}{p({\mathbf{x}}^{(t)})q(\mathbf{x}|{\mathbf{x}}^{(t)})})\\ \\ & =min(1,\frac{\tilde{p}(\mathbf{x})q({\mathbf{x}}^{(t)}|\mathbf{x})}{\tilde{p}({\mathbf{x}}^{(t)})q(\mathbf{x}|{\mathbf{x}}^{(t)})})\end{array}$$

 

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

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

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

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

 

$$\begin{array}{r}q(\mathbf{x}|{\mathbf{x}}^{(t)})=q(\mathbf{x})\end{array}$$

 

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

 

$$\begin{array}{r}q(\mathbf{x}|{\mathbf{x}}^{(t)})=q({\mathbf{x}}^{(t)}|\mathbf{x})\end{array}$$

 

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

 

$$\begin{array}{r}\alpha =min(1,\frac{\tilde{p}(\mathbf{x})}{\tilde{p}({\mathbf{x}}^{(t)})})\end{array}$$

 

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

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

 

 

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

 

$$\begin{array}{r}p(x)\propto 0.5\exp (-(x+2{)}^2)+0.5\exp (-(x-2{)}^2)\end{array}$$

 

 

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

 

$$\begin{array}{rl}q(x│x^{(t)})& =\mathcal{N}(x^{(t)},1)\\ \\ & =\frac{1}{\sqrt{2\pi }}\exp (-\frac{(x-x^{(t)}{)}^2}{2})\end{array}$$

 

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)));