§ 12.8 — 연습문제 + Ch.12 결산

1 개요 — Ch.12 연습의 역할

Ch.12 § 12.8 의 4 문제는 § 12.1~12.7 의 각 축을 실증.

유형	문제	검증 대상
Adaptive Metropolis	1	§ 12.2 의 자동 튜닝
Simulated Tempering	2	§ 12.3 의 다봉 극복
HMC	3	§ 12.4~12.5 의 leapfrog·gradient
Stan + Rejection	4	§ 12.6 의 자동화 + § 10 의 rejection

네 문제가 서로 다른 알고리즘을 같은 bioassay 또는 유사 모형에 적용 — Ch.12 도구들의 성능 비교 기회.

Overview (02-12-0), § 12.1~12.3 (02-12-1), § 12.4~12.6 (02-12-2) 의 마무리편.

2 문제 1 — Bioassay Adaptive Metropolis

문제. Exercise 11.2 의 bioassay Metropolis 를 § 12.2 의 adaptive 알고리즘으로 재실행.

2.1 배경 재확인

Ch.3.7 의 bioassay. \(y_i \sim \mathrm{Bin}(n_i, \theta_i)\), \(\theta_i = \mathrm{logit}^{-1}(\alpha + \beta x_i)\).

\(x = (-0.86, -0.30, -0.05, 0.73)\), \(n_i = 5\), \(y = (0, 1, 3, 5)\).

2.2 단순 Metropolis 의 문제

§ 11.1~11.3 심화 의 구현:

• 시작점 \((\alpha, \beta) = (1, 10)\).
• 점프 \(J = \mathrm{N}(\cdot, c^2 \Sigma)\). \(\Sigma\) 는 Laplace 근사.

단점: \(\Sigma\) 와 \(c\) 를 수작업으로 튜닝. 사후 모양이 다양한 문제에서 매번 반복.

2.3 Adaptive 알고리즘 — § 12.2

Gelman 의 two-phase protocol.

Phase I — Warm-up adaptation:

1. 초기: \(\Sigma^{(0)} = I\), \(c^{(0)} = 2.4/\sqrt{d}\) (\(d = 2\) → \(c = 1.7\)).
2. 매 \(k = 100\) 반복마다:
   • \(\Sigma^{(k)} = \text{cov}(\theta^{(1:k)})\) (현재까지 샘플의 공분산).
   • 수용률 \(\alpha^{(k)}\) 측정.
   • \(c^{(k+1)} = c^{(k)} \cdot \exp\!\left(\frac{\alpha^{(k)} - 0.234}{\sqrt{k}}\right)\) (Robbins-Monro).

Phase II — Fixed:

Warm-up 후 \(\Sigma, c\) 고정. 본 샘플링.

2.4 기대 결과

Bioassay 의 사후 공분산:

\[ \Sigma_{\text{Laplace}} \approx \begin{pmatrix} 1.2 & 3.4 \\ 3.4 & 14.0 \end{pmatrix} \]

(수치는 대략). \(\alpha, \beta\) 강한 양 상관 (\(\approx 0.8\)).

Adaptive 알고리즘이 몇 수백 반복 후 이 공분산 추정 + \(c \approx 1.5{-}1.7\) 로 수용률 0.23 수렴.

단순 Metropolis (isotropic \(\Sigma\)) 대비 자기상관 시간 2~3 배 단축.

직관 — Adaptive 의 가치

수동 튜닝: “문제마다 \(\Sigma\) 를 추정하고 수용률 보고 \(c\) 조정” — 반복 작업.

Adaptive: 기계가 이 과정을 warm-up 에서 자동화. 본 run 은 fixed 로 정합성 유지.

현대 Stan/PyMC 는 HMC 기반이라 공분산 adaptation 이 mass matrix 조정으로 대체. 같은 원리, 다른 대상.

Adaptive Metropolis 는 오래된 기술이지만 gradient 계산 어려운 문제 (블랙박스 모형, Approximate Bayesian Computation 등) 에서 여전히 유용.

3 문제 2 — Cauchy 다봉 + Simulated Tempering

문제. \(y_i \sim \mathrm{Cauchy}(\theta, 1)\), 균등 prior, \(y_1 = 1.3, y_2 = 15.0\).

3.1 (a) 사후 밀도

\[ p(\theta \mid y) \propto \frac{1}{1 + (y_1 - \theta)^2} \cdot \frac{1}{1 + (y_2 - \theta)^2} \]

두 봉: \(\theta = 1.3\) 근처 하나, \(\theta = 15.0\) 근처 하나. Cauchy 의 꼬리가 꼬리 정규 분포처럼 빨리 줄지 않아 두 관측의 영향이 독립적으로 지속.

두 봉 사이 (\(\theta \approx 8\)): 양쪽 관측 모두 꼬리 → 밀도 매우 낮음.

import numpy as np
import matplotlib.pyplot as plt

theta_grid = np.linspace(-10, 25, 1000)
y1, y2 = 1.3, 15.0
log_post = -np.log(1 + (y1 - theta_grid)**2) - np.log(1 + (y2 - theta_grid)**2)
post = np.exp(log_post - log_post.max())
plt.plot(theta_grid, post)
plt.xlabel(r"$\theta$"); plt.ylabel("posterior (unnormalized)")

3.2 (b) Metropolis — 실패

대칭 Cauchy 점프:

\[ J(\theta^* \mid \theta) \propto \frac{1}{1 + ((\theta^* - \theta)/s)^2} \]

Scale \(s\) 조정. \(s\) 너무 작으면 한 봉에 갇힘. \(s\) 너무 크면 수용률 낮음.

\(s = 1\) 실험: 체인이 시작점 근처 한 봉만 탐색. 수천 반복 후에도 두 봉 중 하나만 방문 — 사후 편향.

\(s = 5\): 가끔 건너뜀. 그러나 두 봉 사이 경로 가 너무 얇아 대부분 기각.

결론: 단순 Metropolis 로는 어떤 \(s\) 를 써도 사후를 정확히 샘플하기 어려움.

3.3 (c) Simulated Tempering

사다리: inverse-temperature \(1/T_k = 0.1, 0.2, \ldots, 1.0\) (10 단).

\[ q_k(\theta) = p(\theta \mid y)^{1/T_k} \]

\(1/T = 0.1\) (고온): 사후가 거의 평평 → 봉 사이 이동 쉬움. \(1/T = 1.0\) (원사후): 두 봉 분리.

알고리즘:

1. 현재 \((\theta, s)\) 에서 \(q_s\) 로 Metropolis step.
2. 인접 온도 \(s \pm 1\) 로 jump 제안.

\(s = 0\) (\(T = 1\), 원사후) 인 표본만 최종 추론에 사용.

3.4 (d) 세 방법 비교

방법	왼쪽 봉 (\(\theta \approx 1.3\)) 표본 비율	오른쪽 봉 (\(\theta \approx 15\)) 비율	정확성
Metropolis s=1	100% 또는 0% (한 쪽에 갇힘)	반대	심각 편향
Metropolis s=5	80~100%	0~20%	약간 편향
Simulated tempering	50%	50%	정확

교훈: 다봉 사후는 근본적으로 다른 알고리즘 필요. Tempering, parallel chains, 또는 HMC 의 다봉 확장 (e.g., NUTS with multiple trajectories).

직관 — 왜 Cauchy 모델이 다봉인가

정규 likelihood 라면 두 관측 \((1.3, 15.0)\) 의 평균 \(\approx 8\) 근처 단일 봉. 정규는 이상치에 민감 — 두 값을 평균.

Cauchy 는 두꺼운 꼬리 로 이상치에 둔감. 따라서 \(\theta = 8\) (두 관측 모두 중간) 이 두 관측 각각의 꼬리에서 거리가 멀기보다, 한 관측 근처에 붙는 게 더 그럴듯. 두 봉이 각 관측에 자리 잡음.

일반 원리: Heavy-tail likelihood → 이상치 존중 → 다봉 가능성 높음. 실무에서 robust regression 은 이 성질 활용 (이상치를 별도 mode 로 취급).

Simulated tempering 이 이 구조에 필수적 — 알고리즘 교체 없이 단순 Metropolis 만으로는 한계.

4 문제 3 — Bioassay HMC

문제. Bioassay 에 HMC 구현. (a) gradient 해석+수치 비교, (b) 초기 튜닝, (c) 65% 수용률, (d) ESS ≥ 100 × 4 체인, (e) Ch.3 와 일관성.

4.1 (a) Gradient 유도

Log-posterior (prior = uniform):

\[ \log p(\alpha, \beta \mid y) = \sum_i \left[y_i \log \theta_i + (n_i - y_i) \log(1 - \theta_i)\right] \]

\(\theta_i = \sigma(\eta_i)\), \(\eta_i = \alpha + \beta x_i\), \(\sigma\) = sigmoid.

\(\alpha\) 편미분:

\[ \frac{\partial \log p}{\partial \alpha} = \sum_i (y_i - n_i \theta_i) \]

유도: chain rule. \(d\sigma / d\eta = \sigma(1-\sigma)\), \(d\log\theta / d\alpha = (1-\theta)\) (since \(d\eta/d\alpha = 1\)).

\(\beta\) 편미분:

\[ \frac{\partial \log p}{\partial \beta} = \sum_i x_i (y_i - n_i \theta_i) \]

\(x_i\) 가중 합.

4.2 수치 미분 검증

def num_grad(fn, theta, h=1e-4):
    g = np.zeros_like(theta)
    for i in range(len(theta)):
        t_p = theta.copy(); t_p[i] += h
        t_m = theta.copy(); t_m[i] -= h
        g[i] = (fn(t_p) - fn(t_m)) / (2 * h)
    return g

# 해석 gradient vs 수치 gradient
theta_test = np.array([1.0, 10.0])
g_analytic = analytical_grad(theta_test)
g_numeric = num_grad(log_p_bioassay, theta_test)
print(f"max diff: {np.abs(g_analytic - g_numeric).max():.2e}")

\(10^{-6}\) 이하면 해석 gradient 정확.

4.3 (b) 튜닝 시작값

Mass matrix \(M = \mathrm{diag}(1, 10^{-2})\) — bioassay 에서 \(\alpha\) scale \(\approx 1\), \(\beta\) scale \(\approx 10\).

\(\epsilon = 0.1\), \(L = 10\) 으로 시작.

4.4 (c) 수용률 튜닝

초기 실험: 수용률 \(\approx 0.3\) (65% 보다 낮음). \(\epsilon\) 너무 큼.

시도: \(\epsilon = 0.05\), \(L = 20\) (\(\epsilon L\) 유지). 수용률 \(\approx 0.70\). 통과.

4.5 (d) ESS

4 체인 × 2000 반복, warm-up 1000 후 1000 저장.

\(n_{\mathrm{eff}}\) 측정:

# 수도코드
n_eff_alpha = ess_from_chains(chains_alpha)
n_eff_beta = ess_from_chains(chains_beta)
print(f"ESS alpha: {n_eff_alpha}, ESS beta: {n_eff_beta}")

기대: \(n_{\mathrm{eff}} \approx 500{-}800\) — 4 체인 × 250 독립 equivalent.

체인당 100 목표 달성.

4.6 (e) Ch.3 결과 비교

Ch.3.7 의 격자 계산:

• \(\alpha\) 사후 평균 \(\approx 0.85\).
• \(\beta\) 사후 평균 \(\approx 7.75\).
• LD50 = \(-\alpha/\beta \approx -0.11\).

HMC 결과가 이 값들의 Monte Carlo 오차 \(\pm 0.05\) 이내면 일관.

직관 — HMC 가 이 문제에 과잉인가?

2 차원 사후 (Bioassay) 에서는 격자 계산 또는 Laplace 가 충분. HMC 는 과잉.

그러나 교육적으로 중요:

1. Gradient 유도 연습 — chain rule 응용.
2. Tuning 경험 — 수용률 조정 감각.
3. ESS 개념 — 상관된 표본의 실질 정보.

고차원 복잡 모형에서 HMC 가 본격 유용 — Ch.15 계층 회귀, Ch.18 결측 imputation 등.

Bioassay 는 HMC 의 “hello world” — 작은 예제로 원리 이해.

5 문제 4 — Stan + Rejection Sampling Coverage

문제. \(y_j \sim \mathrm{Bin}(n_j, \theta_j)\), \(\theta_j = \mathrm{logit}^{-1}(\alpha + \beta x_j)\). Prior: \(\alpha \sim t_4(0, 2^2), \beta \sim t_4(0, 1)\). \(J = 10, x_j \sim \mathrm{U}(-1, 1), n_j \sim \mathrm{Poisson}^+(5)\).

5.1 (a) 데이터 + Stan 적합

import numpy as np
rng = np.random.default_rng(128)
J = 10
x = rng.uniform(-1, 1, J)
n = np.maximum(rng.poisson(5, J), 1)  # Poisson+
alpha_true = 1.0
beta_true = 2.0
theta_true = 1 / (1 + np.exp(-(alpha_true + beta_true * x)))
y = rng.binomial(n, theta_true)

Stan 모델:

data {
  int<lower=0> J;
  array[J] int<lower=0> y;
  array[J] int<lower=1> n;
  vector[J] x;
}
parameters {
  real alpha;
  real beta;
}
model {
  alpha ~ student_t(4, 0, 2);
  beta ~ student_t(4, 0, 1);
  y ~ binomial_logit(n, alpha + beta * x);
}

CmdStanPy:

from cmdstanpy import CmdStanModel
model = CmdStanModel(stan_file="logistic.stan")
fit = model.sample(
    data={"J": J, "y": y.tolist(), "n": n.tolist(), "x": x.tolist()},
    chains=4, iter_warmup=1000, iter_sampling=2000, show_progress=False
)
summary = fit.summary()
print(summary[["Mean", "5%", "50%", "95%", "R_hat"]].loc[["alpha", "beta"]])

시각화: 데이터 + 사후 적합 곡선 + 불확실성.

import matplotlib.pyplot as plt
draws = fit.draws_pd()
x_plot = np.linspace(-1, 1, 200)
fig, ax = plt.subplots(figsize=(8, 5))
# 불확실성: 100 posterior samples
for idx in rng.choice(len(draws), 100):
    a, b = draws["alpha"].iloc[idx], draws["beta"].iloc[idx]
    theta_plot = 1 / (1 + np.exp(-(a + b * x_plot)))
    ax.plot(x_plot, theta_plot, color="gray", alpha=0.2, lw=0.5)
# 참 모형
ax.plot(x_plot, 1 / (1 + np.exp(-(alpha_true + beta_true * x_plot))), 'r-', lw=2, label="True")
# 데이터
ax.scatter(x, y / n, s=50 * n / n.max(), label="Data (y/n)")
ax.set_xlabel("x"); ax.set_ylabel("theta"); ax.legend()

5.2 (b) 50% 구간 Coverage

alpha_25, alpha_75 = np.percentile(draws["alpha"], [25, 75])
beta_25, beta_75 = np.percentile(draws["beta"], [25, 75])
alpha_covered = alpha_25 <= alpha_true <= alpha_75
beta_covered = beta_25 <= beta_true <= beta_75
print(f"Alpha 50% interval: [{alpha_25:.2f}, {alpha_75:.2f}], truth = {alpha_true}, covered: {alpha_covered}")
print(f"Beta 50% interval: [{beta_25:.2f}, {beta_75:.2f}], truth = {beta_true}, covered: {beta_covered}")

한 번 실험: 50% 구간이 참값 포함 여부 2 개 중 1 개 정도 예상 (확률 1/2).

정식 SBC (Simulation-Based Calibration): 여러 데이터셋 반복해 coverage 분포 확인. Talts et al. (2018).

5.3 (c) Rejection Sampling

\((\alpha, \beta)\) 사후에서 독립 표본 1000 개.

Envelope: Stan 사후 근사 분포를 이용. 예를 들어 MCMC 로 얻은 사후 평균·공분산으로 multivariate \(t_4\) 제안.

alpha_mean = draws["alpha"].mean()
beta_mean = draws["beta"].mean()
cov_ab = np.cov(draws[["alpha", "beta"]].T)
L_chol = np.linalg.cholesky(cov_ab)

# Envelope: t_4 distribution
def log_envelope(theta):
    from scipy.stats import multivariate_t
    return multivariate_t.logpdf(theta, loc=[alpha_mean, beta_mean], shape=cov_ab, df=4)

def log_target(theta):
    a, b = theta
    lp = -0.5 * np.log(1 + a**2 / (4 * 4))  # t_4(0, 2^2) log
    lp += -0.5 * np.log(1 + b**2 / 4)        # t_4(0, 1) log
    eta = a + b * x
    lp += (y * eta - n * np.log1p(np.exp(eta))).sum()
    return lp

# M 추정: 격자에서 log_target - log_envelope 최대값
test_points = np.column_stack([
    draws["alpha"].sample(1000, random_state=128),
    draws["beta"].sample(1000, random_state=128)
])
log_ratios = np.array([log_target(t) - log_envelope(t) for t in test_points])
log_M = log_ratios.max() + 0.5  # 여유 마진
print(f"log M = {log_M:.2f}")

# Rejection sampling
n_accepted = 0
samples = []
n_trials = 0
while n_accepted < 1000:
    u = rng.chisquare(4, 1)[0] / 4
    z = rng.normal(0, 1, 2)
    prop = np.array([alpha_mean, beta_mean]) + L_chol @ z / np.sqrt(u)
    log_accept = log_target(prop) - log_envelope(prop) - log_M
    if np.log(rng.uniform()) < log_accept:
        samples.append(prop)
        n_accepted += 1
    n_trials += 1
samples = np.array(samples)
print(f"수용률: {n_accepted / n_trials:.3f}")
print(f"Rejection alpha mean: {samples[:, 0].mean():.2f}")
print(f"Stan alpha mean:      {alpha_mean:.2f}")

두 평균이 Monte Carlo 오차 이내 일치 → 두 방법 교차 검증.

5.4 Rejection 의 한계

고차원 (\(d \ge 10\)) 에서는 \(M\) 이 급격히 커짐 → 수용률 0 에 수렴. 이 문제가 독립 표본 원할 때 rejection 우선 의 한계.

실무: Stan 으로 상관 있는 표본 \(\times\) PSIS 또는 thinning 으로 근사적 독립.

직관 — 두 방법의 역할 분담

Stan (HMC): 범용·고차원. 모든 모형에 기본. Rejection: 교차 검증·독립 표본 필요 시. Stan 과 같은 답 나오면 신뢰.

실무: Stan 이 주, rejection 은 검증 도구.

또 다른 용도: posterior 의 variational approximation 을 envelope 으로 사용한 rejection — 대규모 모형에서 독립 표본 얻기. VI + rejection 하이브리드.

6 Ch.12 심화 시리즈 — 4 포스트 논리 지도

Ch.12 의 완결된 경로.

포스트	역할	핵심 메시지
02-12-0 Overview	지도	무작위 보행 극복의 세 축
02-12-1 §12.1~12.3	변형	재매개변수화·튜닝·slice·reversible
02-12-2 §12.4~12.6	HMC	Hamilton 역학·NUTS·Stan
02-12-3 §12.8 (본편)	연습·결산	4 문제 + 4 포스트 지도

Ch.12 의 한 문장:

MCMC 의 무작위 보행 을 극복하려면 (1) 좌표 변환 + 자동 튜닝 으로 Gibbs/Metropolis 를 개선, (2) slice·reversible·tempering 으로 특수 구조 대응, (3) Hamiltonian Monte Carlo 로 물리 역학을 빌려 근본 가속, 그리고 Stan 같은 자동화 환경이 이 복잡도를 사용자로부터 숨긴다.

7 실전 체크리스트 — Ch.12 전체 결산

MCMC 효율 15 항목.

알고리즘 선택

1. 조건부 공액 → Gibbs.
2. 미분 가능 사후 → HMC/NUTS (Stan/PyMC).
3. 다봉 → Tempering / SMC.
4. 차원 변화 → Reversible jump.

모델 설계

5. Non-centered parameterization — 계층 모형 기본.
6. 제약은 변환으로 — Stan 자동 또는 수동 Jacobian.
7. 이산·연속 분리 — HMC 는 연속만.
8. 보조 변수 — \(t\), mixture 등 비공액 구조 해소.

튜닝

9. Mass matrix \(M\) — 사후 공분산 역수 근사.
10. Step size \(\epsilon\) — 수용률 목표치 (MH 23%, HMC 65%).
11. Warm-up 1000 — adaptation 후 fix.
12. 여러 체인 — 병렬 진단.

검증

13. \(\widehat{R}, n_{\mathrm{eff}}\) — 기본 진단.
14. Divergent transitions — HMC funnel 경고.
15. 사후 예측 점검 (Ch.6) — 모형 타당성.

8 관련 주제

선행 지식

• Ch.12 Overview (02-12-0) — 전체 지도
• § 12.1~12.3 (02-12-1) — 재매개변수화·튜닝
• § 12.4~12.6 (02-12-2) — HMC·NUTS·Stan
• Ch.11 시리즈 — MCMC 기초

후속 주제

• Ch.13 Modal and Distributional Approximations — MCMC 대안
• Ch.15 Hierarchical Linear Models — HMC 본격 활용
• Ch.18 Missing Data — Gibbs 의 최대 응용

관련 개념

• Haario, Saksman, Tamminen (2001) — Adaptive Metropolis 표준 논문
• Andrieu & Thoms (2008) — Adaptive MCMC 종합
• Geyer & Thompson (1995) — Simulated tempering
• Neal (2011) — HMC 종합 리뷰
• Hoffman & Gelman (2014) — NUTS
• Carpenter et al. (2017) — Stan
• Talts, Betancourt, Simpson, Vehtari, Gelman (2018) — Simulation-Based Calibration
• Betancourt (2017) — A Conceptual Introduction to HMC