어려운 적분은 가라, 뽑아서 푼다! MCMC와 메트로폴리스-헤이스팅스

베이즈 통계를 공부하다 보면 "사후 확률을 계산하기 너무 어렵다"는 벽에 부딪히는데, MCMC는 이 문제를 "계산하지 않고, 뽑아서(Sampling) 해결한다"는 혁명적인 아이디어입니다. 특히 Metropolis-Hastings 알고리즘은 현대 데이터 과학의 10대 알고리즘으로 꼽힐 만큼 중요합니다.

 

베이즈 통계에서 사후 확률(Posterior)을 구하려면 복잡한 적분을 해야 합니다. 공식이 딱 떨어지는(켤레 사전확률) 운 좋은 경우도 있지만, 현실은 그렇지 않습니다.

"수식을 못 풀겠다면, 그냥 컴퓨터로 수만 번 시뮬레이션해서 맞히면 되지 않을까?"

이것이 바로 MCMC(Markov Chain Monte Carlo)의 핵심입니다. 오늘은 MCMC의 원리와 대표적인 알고리즘인 Metropolis-Hastings를 구현해 보겠습니다.

 

1. 배경 지식: 몬테카를로 적분

복잡한 도형의 넓이(적분)를 구할 때, 수학 공식을 쓰는 대신 무작위로 점을 마구 찍어서(Monte Carlo) 점이 들어간 비율로 넓이를 추정하는 방법입니다.

• 원리: $\frac{\text{도형 안의 점 개수}}{\text{전체 점 개수}} \approx \frac{\text{도형의 넓이}}{\text{전체 사각형 넓이}}$
• 의의: 고차원 적분도 컴퓨터만 있으면 쉽게 근사값을 구할 수 있습니다.

2. MCMC의 탄생: 마르코프 연쇄 + 몬테카를로

우리의 목표는 사후 확률 분포($\pi(\theta|x)$)에서 표본을 뽑는 것입니다. 하지만 그 분포가 어떻게 생겼는지(수식)를 모릅니다.

여기서 마르코프 연쇄의 '안정 상태(Steady State)' 성질을 이용합니다.

1. 우리가 찾고 싶은 사후 확률 분포를 안정 상태로 가지는 마르코프 연쇄를 만듭니다.
2. 이 연쇄를 따라 계속 이동(Transition)하다 보면, 어느 순간 도착하는 곳들이 바로 우리가 원하는 사후 확률 분포의 표본이 됩니다.
3. 이렇게 모은 표본들의 평균을 구하면, 그것이 바로 베이즈 기댓값(모수 추정치)이 됩니다.

3. 메트로폴리스-헤이스팅스 (Metropolis-Hastings) 알고리즘

"어떻게 이동해야 사후 확률 분포에 도달할까?"에 대한 구체적인 방법론입니다.

[알고리즘 4단계]

1. 시작: 아무 값이나 초기값($\theta_0$)으로 잡습니다.
2. 제안 (Propose): 현재 위치에서 근처의 새로운 위치($\theta^*$)를 랜덤하게 제안합니다. (보통 정규분포 사용)
3. 평가 (Accept/Reject): 이동할지 말지 결정합니다.
   • 비율 $r = \frac{\pi(\theta^*|x)}{\pi(\theta_t|x)}$ 계산. (새로운 위치의 확률 / 현재 위치의 확률)
   • $r \ge 1$: 새로운 곳이 더 확률이 높으면 무조건 이동(Accept).
   • $r < 1$: 새로운 곳이 확률이 낮더라도, $r$의 확률로 이동. (가끔은 낮은 곳으로도 가봐야 전체를 훑을 수 있으니까요!)
4. 반복: 이 과정을 수천 번 반복하면, 발자취(Samples)가 결국 사후 확률 분포를 그리게 됩니다.

4. [실습] Python으로 MCMC 구현하기

목표: 평균이 불명확한 어떤 분포(Target)를 MCMC로 찾아내서 시각화하기. (Target 분포는 평균 3, 표준편차 1인 정규분포라고 가정하지만, MCMC는 이를 모르는 상태에서 찾아갑니다.)

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

# 1. 목표 분포 (Target Distribution) - 우리가 알고 싶은 사후 확률이라고 가정
# 실제로는 수식을 모르고 우도*사전확률 값만 알 수 있는 상황임
def target_pdf(x):
    # 예시: 평균 3, 표준편차 1인 정규분포의 확률밀도함수 값 반환 (비례 상수 무시 가능)
    return stats.norm.pdf(x, loc=3, scale=1)

# 2. Metropolis-Hastings 알고리즘
def metropolis_hastings(num_samples, start_value=0):
    samples = [start_value]
    current_x = start_value
    
    for _ in range(num_samples):
        # (1) 제안 (Proposal): 현재 위치에서 정규분포로 다음 위치 추천
        # sigma=0.5는 탐색 보폭(Step size)
        proposed_x = np.random.normal(current_x, 0.5)
        
        # (2) 비율 계산 (Acceptance Ratio)
        # 사후확률 비율 = P(New) / P(Current)
        ratio = target_pdf(proposed_x) / target_pdf(current_x)
        
        # (3) 수락 여부 결정 (Accept/Reject)
        # ratio가 1보다 크면 무조건 수락, 작으면 ratio 확률로 수락
        if np.random.rand() < ratio:
            current_x = proposed_x # 이동
            # (만약 거절되면 current_x는 그대로 유지)
            
        samples.append(current_x)
        
    return np.array(samples)

# 3. 시뮬레이션 실행 (10,000번 반복)
np.random.seed(42)
num_samples = 10000
samples = metropolis_hastings(num_samples, start_value=0)

# 4. 결과 시각화
plt.figure(figsize=(12, 5))

# Trace Plot (샘플링 경로)
plt.subplot(1, 2, 1)
plt.plot(samples)
plt.title("Trace Plot (MCMC Steps)")
plt.xlabel("Iteration"); plt.ylabel("Theta Value")
plt.axhline(3, color='r', linestyle='--', label='True Mean (3)')
plt.legend()

# Histogram (분포 추정)
plt.subplot(1, 2, 2)
plt.hist(samples[1000:], bins=50, density=True, alpha=0.6, color='g', label='MCMC Samples') 
# 앞부분 1000개는 안정화 기간(Burn-in)이라 버림
x = np.linspace(0, 6, 100)
plt.plot(x, target_pdf(x), 'r-', label='True Distribution')
plt.title("Estimated Posterior Distribution")
plt.legend()

plt.tight_layout()
plt.show()

print(f"추정된 평균 (MCMC Mean): {np.mean(samples[1000:]):.4f}")

 

코드 실행 결과 해석

1. Trace Plot (왼쪽 그래프):
   • 0에서 시작했지만, 몇 번 움직이더니 금방 3 근처(True Mean)로 이동해서 그 주변을 왔다 갔다 하는 것을 볼 수 있습니다. 이것이 바로 안정 상태(Steady State)에 도달한 모습입니다.
2. Histogram (오른쪽 그래프):
   • MCMC가 지나간 발자취(초록색 막대)를 모아보니, 우리가 찾고 싶었던 정규분포(빨간 선)와 정확하게 일치합니다.
   • 적분을 하지 않고도, 단지 "이동하고 기록하는(Sampling)" 과정만으로 확률 분포를 완벽하게 찾아냈습니다.

Insight: MCMC가 왜 필요할까?

의료 현장의 데이터는 교과서처럼 예쁜 정규분포를 따르지 않습니다.

예를 들어, 신약의 효과가 환자의 유전자 타입에 따라 복잡하게 얽혀 있는 경우, 수식으로 풀 수 없습니다.

이때 MCMC를 사용합니다.

• 모델: "환자 데이터($x$)를 넣으면 복잡한 확률($\pi$)이 계산되는 블랙박스 함수"만 만듭니다.
• MCMC: 컴퓨터가 이 함수 위를 수만 번 뛰어다니며(Sampling) 지도를 그립니다.
• 결론: "계산해 보니 이 환자에게 약이 효과 있을 확률은 95% 구간에서 82%~88%입니다"라는 정교한 구간 추정이 가능해집니다.

"복잡한 현실을 수식으로 억지로 끼워 맞추지 않고, 시뮬레이션으로 있는 그대로 그려내는 도구." 그것이 MCMC입니다.

 

[Analyst's Insight] MCMC를 쓰는 진짜 이유: "평균의 함정"에서 탈출하기

단순히 평균값 하나(점 추정)만 틱 던져주는 것은 의료 현장에서 매우 위험합니다. "이 수술의 성공률은 80%입니다"라고 했지만, 실제로는 10%~99% 사이에서 널뛰기하는 80%일 수도 있으니까요.

이때 MCMC(베이즈 통계)를 사용하면, 우리가 흔히 쓰는 신뢰구간(Confidence Interval)보다 훨씬 강력하고 직관적인 신용구간(Credible Interval)을 얻을 수 있습니다.

 

1. 정규분포의 강박에서 벗어날 수 있다 (Non-Normal)

• 기존 방식 (빈도주의):
  • 대부분의 데이터가 종 모양(정규분포)일 것이라고 '가정'하고 구간을 계산합니다. ($Mean \pm 1.96 \times SE$)
  • 문제점: 의료 데이터(예: 입원 기간, 의료비)는 정규분포가 아닙니다. 꼬리가 긴(Skewed) 분포거나, 봉우리가 두 개(Multimodal)인 경우, 기존 방식의 구간 추정은 완전히 틀린 답을 줍니다.
• MCMC 방식:
  • 가정 따위는 하지 않습니다.
  • 수만 번 점을 찍어본 결과가 낙타 등처럼 봉우리가 2개라면, 그 모양 그대로 구간을 잡아줍니다.
  • "평균은 의미 없고, 환자는 A그룹(중증) 아니면 B그룹(경증) 둘 중 하나에 속합니다"라는 통찰을 줍니다.

2. '신뢰구간' vs '신용구간' (해석의 직관성)

이 부분이 의사 선생님들을 설득할 때 가장 좋습니다.

• 기존 95% 신뢰구간 (Confidence Interval):
  • "같은 실험을 100번 반복하면, 95번 정도는 이 구간 안에 진짜 평균이 포함될 것이다."
  • (너무 어렵고 와닿지 않습니다. 지금 내 환자가 이 구간에 있을 확률이 95%라는 뜻이 아닙니다.)
• MCMC 95% 신용구간 (Credible Interval):
  • "이 환자의 수치가 이 구간 안에 있을 확률이 진짜로 95%입니다."
  • (우리가 원하던 바로 그 직관적인 대답입니다!)

3. 리스크 관리 (Tail Risk)

• 상황: 신약 A의 효과가 평균적으로는 좋지만, 아주 드물게(0.1% 확률) 치명적인 부작용 수치로 튀는 경우가 있다고 칩시다.
• 기존 방식: 평균과 표준편차만 보면 이 극단적인 0.1%의 위험이 묻혀버립니다.
• MCMC: 시뮬레이션을 돌리다 보면 그 0.1%의 구역에도 점이 찍힙니다.
  • "평균 효과는 좋지만, 최악의 시나리오(Worst Case)에서는 이런 결과가 나올 수도 있습니다"라고 리스크를 정확히 경고할 수 있습니다.

요약: 언제 MCMC를 써야 할까?

1. 데이터가 못생겼을 때: 데이터 분포가 정규분포가 아니고 찌그러져 있거나 봉우리가 여러 개일 때.
2. 직관적인 보고가 필요할 때: "이 구간 안에 있을 확률이 95%입니다"라고 명확하게 말하고 싶을 때.
3. 최악의 상황(Risk)을 봐야 할 때: 단순히 평균이 아니라, 발생 가능한 모든 시나리오의 분포를 보고 싶을 때.

결국 MCMC는 "복잡한 현실을 수식으로 억지로 단순화하지 않고, 있는 그대로(Distribution) 보여주기 위해" 사용하는 가장 정교한 도구입니다.