[MCMC] MCMC 개요

어떤 확률분포를 가진 랜덤변수에서 데이터를 생성하는 과정을 샘플링(sampling)이라고 한다. 파이썬(Python)이나 매트랩(Matlab)과 같은 많은 컴퓨터 언어에서는 기본적으로 가우시안 확률분포나 균등 확률분포와 같은 표준 형태의 확률밀도함수로부터 샘플을 추출할 수 있는 함수를 제공한다. 하지만 그렇지 않은 경우에는 어떤 방법으로 샘플을 추출해야할까.

 

 

먼저 직접적인 방법이 있다 (https://pasus.tistory.com/49). 누적분포함수(cumulative distribution function)의 역함수를 이용하는 방법으로서 정확한 샘플링 방법이다. 하지만 랜덤변수가 다차원(multi-dimension)을 갖거나 복잡한 확률밀도함수를 갖는 경우에는 이 방법을 적용하기가 어렵다.

두번째는 거절 샘플링(rejection sampling) 방법이다. 먼저 샘플링하기 쉬운 확률분포(가우시안 분포나 균등 분포)를 골라 그로부터 샘플링한다. 그 다음 각 샘플에 대해서 수용확률을 계산하여 그 확률로만 샘플을 수용한다. 이때 수용확률은 샘플이 원하는 목표 확률분포에서 나온 샘플일 확률에 비례한다. 이 방법은 간단하고 정확하며 누적함수의 역함수를 계산할 필요가 없는 장점이 있다. 하지만 목표 분포함수와 실제 샘플링을 수행한 분포함수가 근사할 때에만 효율적이고, 고차원 랜덤변수에서는 잘 적용되지 않는다.

세번째는 중요 샘플링(importance sampling) 방법이 있다 (https://pasus.tistory.com/52). 샘플을 추출하여 기댓값(expectation)을 계산하는 경우에 사용할 수 있다. 중요 샘플링은 기댓값을 계산하고자 하는 확률분포함수는 알고 있지만 샘플을 생성하기가 어려울 때 해당 분포함수 대신에 샘플을 생성하기가 쉬운 다른 분포함수(가우시안 분포나 균등 분포)를 이용해 기댓값을 추정하는 방법이다.

네번째는 이 글에서 소개할 MCMC(Markov Chain Monte Carlo) 방법이다. 마르코프 체인의 장기 거동이 특정 확률분포에 수렴한다는 성질을 활용하여, 복잡한 확률분포에서 간접적으로 샘플링할 수 있는 기법이다. 직접 샘플링이 불가능하거나 거절 샘플링을 쓰면 효율이 극도로 떨어지는 분포에서도 잘 작동한다.

MCMC는 두개의 용어가 결합된 단어이므로 우선 각각을 설명하고자 한다. 먼저 두번째 용어인 몬테카를로부터 시작한다.

몬테카를로 방법은 무작위 샘플링을 통해 수치적 결과를 근사하는 알고리즘의 총칭이다. 이 방법의 핵심 아이디어는 플고자 하는 문제를 확률분포와 확률변수로 표현하고, 해당 분포에서 무작위 샘플을 대량으로 생성한 후, 생성된 샘플의 통계값(평균이나 분산 등)을 이용해 솔루션(적분값이나 확률 등)의 근사값을 계산하는 것이다. 여기서 '풀고자 하는 문제'는 해석적으로 풀기 어려운 문제로서 확정적(deterministic) 문제가 될 수도 있고 확률적(stochastic)문제가 될 수도 있다.

 

 

예를 들어보자. 몬테카를로 방법을 설명할 때 가장 많이 등장하는 예제로서 면적을 계산하는 문제가 있다. 면적을 계산하는 문제는 확정적 문제인데 이를 확률적 문제로 변환하는 것이 몬테카를로 방법의 핵심이다. 원점에 중심이 있고 반지름이 \(1\) 인 원의 면적은 \(\pi\) 다. 이 문제를 몬테카를로 방법으로 풀기 위해서 우선 \([-1 \ 1 ] \times [-1 \ 1 ] \) 영역에서 균등하게 \(N\) 개의 샘플을 발생시킨다. 샘플을 특정 영역의 균등 분포에서 무작위로 추출한다는 것은 곧 해당 영역 내부의 어느 위치이든 동일한 확률로 선택한다는 뜻이므로 원 안에 샘플이 들어갈 확률은 그 원의 면적에 비례하게 된다. 만약 원 안에 들어간 샘플 수가 \(M\) 개라면, 한 변 길이가 \(2\) 인 정사각형 면적과 원의 면적의 비는 전체 샘플 갯수와 원안에 있는 샘플 갯수의 비와 같다. 따라서 원의 면적은 \(area=4M/N\) 이 된다.

 

 

매트랩 코드는 다음과 같다.

 

% Monte Carlo example 1 - circle area computation
% st.watermelon

clear all

% sample generation [-1 1] x [-1 1]
N = 1e4;
x = rand(N,1)*2-1;
y = rand(N,1)*2-1;  

% check x^2 + y^2 <= 1
is_inCircle = (x.^2 + y.^2 <= 1);

% plotting
theta = linspace(0, 2*pi, 200);
r = 1;
x_c = r * cos(theta);
y_c = r * sin(theta);
plot(x,y,'g.', x_c, y_c, 'r')

% count samples inside cicle
M = sum(is_inCircle);

area = 4*M/N;

fprintf('number of total samples N = %d\n', N);
fprintf('number of samples inside cicle M = %d\n', M);
fprintf('cicle area = %.6f\n', area);
fprintf('error = %.6f\n', (pi-area));

 

>> mc_example1
number of total samples N = 10000
number of samples inside cicle M = 7862
cicle area = 3.144800
error = -0.003207

 

샘플 \(N=10,000\) 개를 이용했을 때 오차가 약 \(0.003\) 정도가 생겼는데 샘플 수를 증가시키면 이 오차는 줄어든다.

여기서는 단순하게 원의 면적을 구하는 문제를 예시로 들었지만 복잡한 도형의 경우 면적을 계산하는 것보다 샘플이 도형안에 존재하는지 밖에 존재하는지 판단하는게 훨씬 쉽다는 점을 생각하면 이 방법이 효율적임을 알 수 있을 것이다. 특히 복잡하고 고차원의 적분의 경우 해석적으로 접근이 어려울 때가 많지만 몬테카를로는 샘플링을 이용해 이런 적분을 근사할 수 있다.

이와 같이 실제 문제(주가 변동, 입자 물리 실험 등)가 확정적이든 확률적이든, 이론적 공식을 이용하기 보다는 실제와 유사한 무작위 시뮬레이션을 이용하는 것이 더 직관적이고 간편한 경우가 많다. 하지만 몬테카를로 시뮬레이션은 특정한 확률분포(정규분포, 균등분포 등)에서 샘플링할 수 있을 때만 사용할 수 있다. 만약 그 확률분포에서 직접 샘플을 추출하기 어려울 때 사용할 수 있는 방법이 MCMC다. MCMC는 마르코프 체인이 시간이 지남에 따라 우리가 원하는 확률분포로 수렴하도록 만든 뒤 샘플링하는 방법이다.

 

 

다음으로 MCMC의 첫번째 용어인 마르코프 체인에 대해서 알아보자. 연속공간(continuous space)에서 정의된 랜덤변수의 프로세스 또는 시퀀스 \(\mathbf{x}_{0:t} = \{ \mathbf{x}_0, \mathbf{x}_1, ..., \mathbf{x}_t \}\) 가 다음 조건을 만족하면 마르코프 프로세스 또는 마르코프 시퀀스(Markov sequence)라고 정의한다 (https://pasus.tistory.com/145 ).

 

\[ \begin{align} p(\mathbf{x}_t \vert \mathbf{x}_0, \mathbf{x}_1, ... , \mathbf{x}_{t-1})=p( \mathbf{x}_t \vert \mathbf{x}_{t-1}) \end{align} \]

 

여기서 조건부 확률밀도함수 \(p( \mathbf{x}_t \vert \mathbf{x}_{t-1}) \) 를 천이커널(transition kernel)이라고 하고 \(K( \mathbf{x}_t \vert \mathbf{x}_{t-1})\) 로 표시한다. 천이커널은 시불변(time-invariant)이라고 가정한다. 정의에 의하면 마르코프 시퀀스는 현재 상태가 직전 상태에만 의존하는 시퀀스이다.

 

 

시간스텝 \(t\) 에서의 확률밀도함수는 다음과 같이 계산할 수 있다.

 

\[ \begin{align} p(\mathbf{x}_t ) &= \int p(\mathbf{x}_{t-1}, \mathbf{x}_t) \ d \mathbf{x}_{t-1} \\ \\ &= \int K(\mathbf{x}_t \vert \mathbf{x}_{t-1}) \ p(\mathbf{x}_{t-1}) \ d \mathbf{x}_{t-1} \end{align} \]

 

만약 랜덤변수가 이산공간(discrete space)에서 정의되고 랜덤변수의 시퀀스 \(s_{0:t}= \{ s_0, s_1, ..., s_t \}\) 가 다음 조건을 만족하면 마르코프 체인(Markov chain)이라고 한다.

 

\[ \begin{align} p(s_t \vert s_0, s_1, ..., s_{t-1} )=p(s_t \vert s_{t-1} ) \end{align} \]

 

여기서 이산 상태 \(s_t\) 는 이산 상태의 집합 \(s_t \in \mathcal{S}= \{ S^{(1)}, S^{(2)}, ..., S^{(n)} \}\) 에서만 값을 취할 수 있다. 따라서 \(p(s_t \vert s_{t-1} )\) 는 조건부 확률이 되는데, 이를 \(T=T(s_t\vert s_{t-1} )\) 로 표시하고 천이행렬(transition matrix)이라고 한다. 천이행렬이 시불변이라고 가정하면 \(T\) 는 상수행렬이 되는데 각 구성요소는 다음과 같이 정의된다.

 

\[ \begin{align} T_{ij}=T(s_t=S^{(j)} \vert s_{t-1}=S^{(i)} ) \end{align} \]

 

이와 같이 상태 공간 유형에 따라 천이확률(transition probability)을 천이행렬로 나타낼지 또는 천이커널함수로 나타낼지가 달라진다.

문헌에 따라서는 마르코프 시퀀스와 마르코프 체인을 구별하지 않고 모두 마르코프 체인이라고 하기도 한다. MCMC의 개념을 이해할 때는 이산 상태공간이 보다 직관적이어서 이를 위주로 설명하겠지만 MCMC 알고리즘은 연속공간이든 이산공간이든 모두 적용된다.

정의에 의하면 천이행렬 \(T\) 는 \(n \times n\) 행렬이며, 각 행은 특정 상태에서 다음 상태로 이동할 수 있는 확률을 나타내기 때문에 행의 합은 \(1\) 이 되어야 한다.

 

\[ \begin{align} \sum_{j=1}^n T_{ij} =1 \end{align} \]

 

또한 시간스텝 \(t\) 에서 특정 이산 상태의 확률은 다음과 같이 계산할 수 있다.

 

\[ \begin{align} p(s_t=S^{(j) } ) &= \sum_{i=1}^n p(s_t = S^{(j) }, s_{t-1}=S^{(i)} ) \\ \\ &= \sum_{i=1}^n p(s_{t-1}= S^{(i)}) T(s_t=S^{(j)} \vert s_{t-1}=S^{(i)} ) \end{align} \]

 

따라서 시간스텝 \(t\) 에서 이산 상태의 확률분포는 다음과 같이 행렬/벡터의 곱으로 계산할 수 있다.

 

\[ \begin{align} \mathbf{p}(s_t )= \mathbf{p}(s_{t-1}) T \end{align} \]

 

여기서 확률분포 \(\mathbf{p}(s_t )\) 는 행 벡터로서 \(\mathbf{p}(s_t )=[ p(s_t=S^{(1) } ) \ \ \cdots \ \ p(s_t=S^{(n) } ) ]\) 이다. 초기 상태의 확률분포가 \( \mathbf{p}(s_0 )\) 이라면 시간스텝 \(t\) 에서의 확률분포는 다음과 같이 계산할 수 있다.

 

\[ \begin{align} \mathbf{p}(s_t ) &= \mathbf{p}(s_{t-1} ) T = \mathbf{p}(s_{t-2} ) T^2 = \cdots \\ \\ &= \mathbf{p}(s_0 ) T^t \end{align} \]

 

만약 시간이 충분히 흘러서 \(t \to \infty\) 일 때 확률분포가 일정한 값으로 수렴하면 마르코프 체인의 확률분포를 정상분포(stationary distribution)라고 한다.

 

\[ \begin{align} \lim_{t \to \infty} \mathbf{⁡p}(s_t ) = \mathbf{p}(s_{\infty} ) = \mathbf{p}_{ss} \end{align} \]

 

\(\mathbf{p}_{ss} \) 를 정상분포라고 부르는 이유는 \(\mathbf{p}_{ss} \) 가 천이행렬 \(T\) 로 정의된 마르코프 체인의 극한 분포라면 \(\mathbf{p}_{ss} = \mathbf{p}_{ss} T\) 를 만족하므로 시간이 지나도 일정한 값으로 고정되어 있기 때문이다. 또한 \(\mathbf{p}_{ss} =\mathbf{p}_{ss} T\) 라는 사실에서 정상분포 \(\mathbf{p}_{ss} \) 는 마르코프 체인의 천이행렬 \(T\)의 (좌)고유벡터이며 고유값은 \(1\) 이라는 것을 알 수 있다.

예를 들어 \(3\) 개의 상태를 갖는 마르코프 체인의 천이행렬이 다음과 같을 때,

 

\[ \begin{align} T= \begin{bmatrix} 0.4 & 0 & 0.6 \\ 0.5 & 0.3 & 0.2 \\ 0 & 1 & 0 \end{bmatrix} \end{align} \]

 

각 행의 합은 \(1\) 이 되며, 상태 \(S^{(2)}\) 에서 상태 \(S^{(3)}\) 으로 이동할 확률은 \(0.2\) 이고, 상태 \(S^{(3)}\) 에서 상태 \(S^{(2)}\) 로 이동할 확률은 \(1\) 이 된다는 것을 알 수 있다. 만약 초기 상태의 확률분포가 \(\mathbf{p}(s_0 )=[0.5 \ \ 0.2 \ \ 0.3]\) 이라면 다음 상태의 확률분포는 다음과 같다.

 

\[ \begin{align} \mathbf{p}(s_1 )= \mathbf{p}(s_0) T =[0.3 \ \ 0.36 \ \ 0.34] \end{align} \]

 

그 다음 시간스텝에서의 확률분포는 다음과 같다.

 

\[ \begin{align} \mathbf{p}(s_2 ) &= \mathbf{p}(s_1) T =[0.3 \ \ 0.36 \ \ 0.34] \\ \\ \mathbf{p}(s_3 ) &= \mathbf{p}(s_2) T =[0.3 \ \ 0.36 \ \ 0.34] \end{align} \]

 

즉 시간이 충분히 흐르면 확률분포는 \( [0.3 \ \ 0.36 \ \ 0.34]\) 로 수렴한다. 또한 행렬 \(T\) 의 고유값이 \(1\) 일 때 (좌)고유벡터는 \( [0.3 \ \ 0.36 \ \ 0.34]\) 이어서 수렴된 확률분포와 같다는 것을 확인할 수 있다.

이와 같은 수렴 결과는 MCMC 시뮬레이션에서 근본적인 역할을 한다. 만약 마르코프 체인이수렴 특성을 갖는다면 어떤 시작점에서든 체인은 정상분포 \(\mathbf{p}_{ss}\) 로 수렴한다. 연속 상태공간에서 마르코프 시퀀스가 수렴 특성을 갖는다면 천이행렬 \(T\) 는 천이커널 \(K\) 가 되고 \( p_{ss} (\mathbf{x})\) 는 이에 대응하는 고유함수가 된다.

 

\[ \begin{align} p_{ss} (\mathbf{x}^\prime )= \int K(\mathbf{x}^\prime \vert \mathbf{x}) p_{ss} (\mathbf{x}) \ d\mathbf{x} \end{align} \]

 

 

 

그렇다면 마르코프 체인이 정상분포로 수렴하기 위한 조건은 무엇일까.

이산 상태공간에서는 환원불가능성(irreducibility)과 무주기성(aperiodicity), 재귀성(recurrence)의 조건을 만족하면 확률분포가 정상상태로 수렴한다는 것이 증명되었다. 여기서 환원불가능성이란 어느 상태에서 시작하더라도 충분히 긴 시간이 지나면 다른 모든 상태에 도달할 확률이 \(0\) 이 아님을 의미하고, 무주기성이란 특정 상태에서 출발했을 때, 다시 그 상태 로 돌아오는 시간이 일정 주기를 갖지 않아야 한다는 뜻이다. 또한 재귀성이란 체인이 특정 상태(또는 상태 집합)로 되돌아올 수 있어야 한다는 뜻으로 이산 상태가 유한개로 구성되어 있다면 재귀성은 보장된다.

연속 상태공간에서는 위와 같은 환원불가능성, 무주기성, 재귀성의 개념을 연속 공간으로 확장하여 사용한다.

수렴된 확률분포가 원하는 확률분포 또는 확률밀도함수가 되기 위해서는 추가적으로 세부균형(detailed balance) 조건을 만족해야 한다. 세부균형은 각 상태 \(\mathbf{x}\) 와 \(\mathbf{x}^\prime\) 의 쌍에 대해 \(\mathbf{x}\) 에 도달한 다음 \(\mathbf{x}^\prime\) 으로 이동하는 것이 \(\mathbf{x}^\prime\) 에 도달한 다음 \(\mathbf{x}\) 로 이동하는 것과 확률 또는 확률밀도함수가 같다는 것을 의미한다.

그렇다면 수렴된 분포가 우리가 원하는 확률분포가 되도록 마르코프 체인을 어떻게 구성할 수 있을까?

다음에는 MCMC에서 가장 널리 사용되는 알고리즘인 메트로폴리스-헤이스팅스(Metropolis-Hastings) 알고리즘과 깁스 셈플링(Gibbs Sampling), 그리고 해밀토니안 몬테카를로(HMC, Hamiltonian Monte Carlo) 방법에 대해서 알아보기로 한다.