2 분 소요

본 포스팅에서는 베이지안 추론에서 핵심이 되는 MCMC(Markov Chain Monte Carlo) 알고리즘의 원리와 구현 방식을 살펴보겠습니다.


Introduction

많은 베이지안 모델은 사후분포를 구하고 이를 이용해 원하는 값을 계산하는 것을 목표로 합니다.

관측 데이터 $\mathbf{y}$가 주어졌을 때, 파라미터 $\theta$의 사후분포는 다음과 같이이 Bayes’ Rule에 의해 계산됩니다.

\[p(\theta \mid \mathbf{y})=\dfrac{p(\mathbf{y}\mid \theta)\,p(\theta)}{p(\mathbf{y})}\]


하지만 분모에 등장하는 정규화 상수 혹은 Evidence $p(\mathbf{y})$는 다음과 같은 고차원 적분으로 정의됩니다.

\[p(\mathbf{y})=\int p(\mathbf{y}∣\theta)p(\theta)d\theta\]


따라서 대부분의 실제 모델에서는 이 적분을 closed form으로 계산하기가 거의 불가능합니다( = intractable). 특히 파라미터 $\theta$의 차원이 커짐에 따라 차원의 저주로 인해 수치적분 방식은 계산의 비용이 아주 급격하게 증가하게 됩니다.

그래서 직접적 계산이 아닌, 다른 방식의 접근이 필요합니다.


Markov Chain

Markov Chain은 현재 상태가 오직 이전 상태에만 의존하는 확률 과정(stochastic process)입니다.


markov_chain
Markov Chain


Markov Chain은 다음과 같은 구성 요소로 정의됩니다:

  1. Markov Property(무기억성)
    • 현재 상태 $X_t$만이 다음 상태 $X_{t+1}$에 영향을 미치며, 과거 상태 $X_{t-1},X_{t-2},…$에는 직접 의존하지 않습니다.
      수식: $P(X_{t+1}=x’\mid X_t=x, …) = P(X_{t+1}=x’\mid X_t=x)$
  2. State Space
    • 이산 또는 연속 공간이며, 각각의 $x$가 가능한 상태를 나타냅니다.
  3. Transition Kernel(전이 커널)
    • $K(x\to x’) = P(X_{t+1}=x’\mid X_t=x)$.
    • n-step transition: $K^n(x\to x’) = P(X_{t+n}=x’\mid X_t=x)$.
    • Chapman–Kolmogorov: $K^{m+n}(x\to z) = \int K^m(x\to y)\,K^n(y\to z)\,dy$.
  4. Stationary Distribution(정착 분포)
    • 분포 $\pi(x)$가 전이 커널의 고정점이 되려면
      $\pi(x’) = \int \pi(x)\,K(x\to x’)\,dx$
      를 만족해야 합니다.
    • 유한 상태 체인의 경우 $\pi = \pi K$로 표현합니다.
  5. Ergodicity
    • Irreducibility(귀환성), Aperiodicity(주기 없음), Positive recurrence(양의 재귀성)
    • 이 세 조건이 만족되면
      $\lim_{t\to\infty}P(X_t=x)=\pi(x)$가 성립합니다.
  6. Detailed Balance
    • Metropolis–Hastings에서 자주 사용되는 조건으로
      $\pi(x)\,K(x\to y)=\pi(y)\,K(y\to x)$
      가 성립하면 $\pi$는 정착분포가 됩니다.

MCMC에서는 이 전이 커널 $K$를 목표 사후분포 $p(\theta\mid y)$가 Stationary distribution가 되도록 설계합니다. 충분히 오랜시간 동안 반복하면, 상태 분포가 $p(\theta\mid y)$에 가까워지게 됩니다.


MCMC (Markov Chain Monte Carlo)

MCMC는 말 그대로 Markov Chain을 이용하여 Monte Carlo Sampling을 수행하는 알고리즘입니다.

Objective

고차원의 복잡한 사후분포 $p(\theta \mid \mathbf{y})$에서 직접 샘플링하는 것은 어려우므로,
적절한 전이 커널을 구성하여 Markov Chain을 생성하고,
이로부터 반복적으로 샘플링함으로써 근사적으로 사후분포를 재현합니다.


Metropolis-Hastings Algorithm

MH algorithm
Intuitutive figure of Metropolis-Hastings algorithm

Basic idea

  1. 현재 상태 $\theta^{(t)}$에서 시작하여
  2. 제안 분포 $q(\theta’ \mid \theta^{(t)})$로부터 새로운 후보 $\theta’$를 샘플링
  3. 다음 확률로 $\theta’$를 수용(Accept)
\[\alpha = \min\left(1,\ \dfrac{p(\mathbf{y} \mid \theta')\,p(\theta')\,q(\theta^{(t)} \mid \theta')}{p(\mathbf{y} \mid \theta^{(t)})\,p(\theta^{(t)})\,q(\theta' \mid \theta^{(t)})}\right)\]
  1. $\theta^{(t+1)} = \theta’$ with probability $\alpha$, otherwise $\theta^{(t+1)} = \theta^{(t)}$

Note: $q$가 대칭 분포(예: 정규분포)라면 $q$ 항이 약분되며, 이 경우를 Metropolis 알고리즘이라 부름

Pseudo code

for t in range(T):
    theta_prime = propose(theta[t])  # q(·|theta[t])에서 샘플링
    r = min(1, posterior(theta_prime) / posterior(theta[t]))
    if np.random.rand() < r:
        theta[t+1] = theta_prime
    else:
        theta[t+1] = theta[t]

Gibbs Sampling

Gibbs 샘플링은 조건부 분포를 통해 각 변수에 대해 순차적으로 샘플링하는 방법입니다.

조건 • 파라미터 벡터 $\theta = (\theta_1, \theta_2, …, \theta_d)$가 있을 때, • 각 $p(\theta_i \mid \theta_{-i}, \mathbf{y})$를 샘플링할 수 있어야 함

Pseudo code

python initialize θ = (θ₁, θ₂, …, θ_d) for t in range(T): for i in range(d): θ[i] = sample_from_conditional(θ[-i], y)

주의: Full conditional 분포를 정확히 유도할 수 있을 때만 사용 가능하며, conjugate prior가 유리함


Convergence and Diagnostics

MCMC 샘플의 수렴성 여부는 다음을 통해 확인할 수 있습니다.

  • Burn-in: 초기 몇 천 개 샘플은 버림
  • Trace plot: 파라미터 샘플 경로 시각화 $\rightarrow$ 상관이 없는 무작위한 패턴을 보여야 함.
  • Autocorrelation plot: 자기상관 구조 시각화 $\rightarrow$ 자기상관이 사라져야 함.

References

업데이트: