§ 11.8 — 연습문제 + Ch.11 결산

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

Ch.11 § 11.8 의 7 문제는 이론·구현·실증 세 축을 반복.

유형	문제
이론 증명	1 (MH 정상성), 5 (var+ 불편성), 6 (ESS 유도)
알고리즘 구현	2 (Metropolis bioassay), 3 (Gibbs 3 모형), 4 (분산 계층), 7 (survey)

이 포스트는 각 문제의 수식 유도 + 구현 골자 + 직관. 단순 답안이 아닌 “이 연습을 풀고 나면 무엇을 얻는가” 에 집중한다.

Overview (02-11-0), § 11.1~11.3 (02-11-1), § 11.4~11.6 (02-11-2) 와 함께 Ch.11 심화 시리즈 4 포스트로 완결된다.

2 문제 1 — Metropolis-Hastings 정상성 증명

문제. MH 의 정상 분포가 실제로 \(p(\theta \mid y)\) 임을 보이시오.

2.1 풀이 — Detailed Balance 재확인

§ 11.1~11.3 심화 (02-11-1) 에서 이미 수행한 증명의 요지.

전이 밀도 (proposal \(J\), acceptance \(\alpha\)):

\[ T(\theta' \mid \theta) = J(\theta' \mid \theta) \alpha(\theta, \theta') + \delta(\theta - \theta') \left[1 - \int J(\tilde{\theta} \mid \theta) \alpha(\theta, \tilde{\theta}) d\tilde{\theta}\right] \]

수용 확률:

\[ \alpha(\theta, \theta') = \min\!\left(1, \frac{p(\theta' \mid y) J(\theta \mid \theta')}{p(\theta \mid y) J(\theta' \mid \theta)}\right) \]

2.2 Detailed Balance

\(\theta \ne \theta'\) 인 경우 (\(\delta\) 항 0). WLOG \(p(\theta' \mid y) J(\theta \mid \theta') \ge p(\theta \mid y) J(\theta' \mid \theta)\) 로 라벨. 그러면:

• \(\theta \to \theta'\): \(\alpha = 1\). \(T(\theta' \mid \theta) = J(\theta' \mid \theta)\).
• \(\theta' \to \theta\): \(\alpha = \frac{p(\theta \mid y) J(\theta' \mid \theta)}{p(\theta' \mid y) J(\theta \mid \theta')}\). \(T(\theta \mid \theta') = J(\theta \mid \theta') \cdot \alpha\).

양변 계산:

\[ p(\theta \mid y) T(\theta' \mid \theta) = p(\theta \mid y) J(\theta' \mid \theta) \]

\[ p(\theta' \mid y) T(\theta \mid \theta') = p(\theta' \mid y) J(\theta \mid \theta') \cdot \frac{p(\theta \mid y) J(\theta' \mid \theta)}{p(\theta' \mid y) J(\theta \mid \theta')} = p(\theta \mid y) J(\theta' \mid \theta) \]

같음. Detailed balance 성립 → \(p(\theta \mid y)\) 가 정상 분포. \(\square\)

2.3 Gibbs 의 특수 케이스

Gibbs 에서 제안 \(J_j(\theta_j^* \mid \theta_{-j}^{t-1}) = p(\theta_j^* \mid \theta_{-j}^{t-1}, y)\) 대입 시 \(\alpha = 1\) — 항상 수용. 역시 정상 분포 조건 만족.

직관 — 이 증명이 가르치는 것

증명의 핵심은 라벨링 자유. \(p(\theta') J(\theta \mid \theta') \ge p(\theta) J(\theta' \mid \theta)\) 라고 가정해도 일반성 손실 없음 — 대칭 상황에서는 라벨만 바꾸면 됨.

MH 의 “\(r\) = 비율의 비율” 구조가 왜 정확한지: 어느 쪽이 더 확률 높은가에 상관없이 detailed balance 를 유지하도록 설계됨.

이 증명을 이해하면 임의의 MH 변형 (slice sampling, Hamiltonian MC) 의 정합성도 같은 방식으로 검증 가능.

3 문제 2 — Bioassay Metropolis 재현

문제. § 3.7 의 bioassay 예제를 Metropolis 로 재현. 시작점, 점프 규칙, 로그 밀도 계산 명시.

3.1 데이터

\(x_i\) = 로그 용량, \(n_i\) = 처리 동물 수, \(y_i\) = 사망 수.

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

3.2 모형

\[ y_i \sim \mathrm{Bin}(n_i, \theta_i), \quad \theta_i = \mathrm{logit}^{-1}(\alpha + \beta x_i) \]

Prior: \(p(\alpha, \beta) \propto 1\) (비정보).

3.3 로그 사후

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

\(\theta_i\) 가 \(\alpha + \beta x_i\) 에 의존.

3.4 Metropolis 구현

• 시작점: MLE 근처 또는 과분산. 예: \((\alpha^0, \beta^0) = (1, 10)\).
• 점프 규칙: \((\alpha^*, \beta^*) \sim \mathrm{N}((\alpha^{t-1}, \beta^{t-1}), c^2 \Sigma)\). \(\Sigma\) 는 사후 공분산 추정 (Laplace 근사).
• \(c\) 선택: Roberts-Gelman-Gilks 권고 \(c = 2.4/\sqrt{2} \approx 1.7\).
• 로그 스케일: \(\log r = \log p^* - \log p^{t-1}\).

3.5 Python 구현

import numpy as np
from scipy.stats import norm
from numpy.linalg import cholesky

rng = np.random.default_rng(112)

x_bio = np.array([-0.86, -0.30, -0.05, 0.73])
n_bio = np.array([5, 5, 5, 5])
y_bio = np.array([0, 1, 3, 5])

def log_posterior(params):
    a, b = params
    eta = a + b * x_bio
    theta = 1 / (1 + np.exp(-eta))
    eps = 1e-12
    return (y_bio * np.log(theta + eps) + (n_bio - y_bio) * np.log(1 - theta + eps)).sum()

# Laplace 추정 대용 수작업 공분산
Sigma_prop = np.array([[1.0, 0.5], [0.5, 3.0]])
L = cholesky(Sigma_prop)
c = 2.4 / np.sqrt(2)

n_iter = 10_000
n_chain = 4
chains = []
for ch in range(n_chain):
    params = np.zeros((n_iter, 2))
    params[0] = rng.normal([1.0, 10.0], [0.5, 2.0])  # 과분산 시작
    accept = 0
    for t in range(1, n_iter):
        prop = params[t-1] + c * (L @ rng.normal(0, 1, size=2))
        log_r = log_posterior(prop) - log_posterior(params[t-1])
        if np.log(rng.uniform()) < log_r:
            params[t] = prop
            accept += 1
        else:
            params[t] = params[t-1]
    print(f"Chain {ch+1}: acceptance = {accept/n_iter:.3f}")
    chains.append(params[n_iter//2:])

all_samples = np.concatenate(chains)
print(f"alpha mean: {all_samples[:, 0].mean():.2f}, sd: {all_samples[:, 0].std():.2f}")
print(f"beta mean:  {all_samples[:, 1].mean():.2f}, sd: {all_samples[:, 1].std():.2f}")
print(f"LD50 mean:  {(-all_samples[:, 0] / all_samples[:, 1]).mean():.3f}")

기대 결과: \(\alpha \approx 1\), \(\beta \approx 10\), LD50 \(\approx -0.1\), 수용률 30~35%. Ch.3.7 의 정확 격자 결과와 일치.

3.6 핵심 체크

• 과분산 시작점에서 여러 체인 수렴 확인 → \(\widehat{R}\).
• 수용률이 25~40% 구간이면 점프 크기 적절.
• \(\log r\) 로 계산해 overflow 방지.

4 문제 3 — 6 Machine Quality Control

문제. Table 11.4 의 6 머신 × 5 측정 데이터에 separate / pooled / hierarchical 3 모형 적합. 각 모형에서 (i) 6 번째 머신 평균 사후, (ii) 6 번째 머신 다음 측정 예측, (iii) 관측 없는 7 번째 머신 사후 보고.

4.1 데이터 (Table 11.4)

Machine	Measurements
1	83, 92, 92, 46, 67
2	117, 109, 114, 104, 87
3	101, 93, 92, 86, 67
4	105, 119, 116, 102, 116
5	79, 97, 103, 79, 92
6	57, 92, 104, 77, 100

4.2 세 모형

4.2.1 Separate

각 머신 독립: \(y_{ij} \sim \mathrm{N}(\theta_j, \sigma^2)\), 모든 \(\theta_j\) 는 완전 독립 (계층 없음). 사실상 풀링 없음.

\(\theta_6 \mid y \sim \mathrm{N}(\bar{y}_{\cdot 6}, \sigma^2/n_6)\).

7 번째 머신: 예측 불가 (사전 정보 없음).

4.2.2 Pooled

모든 머신 같은 평균: \(y_{ij} \sim \mathrm{N}(\mu, \sigma^2)\). 즉 \(\theta_1 = \ldots = \theta_6 = \mu\).

\(\mu \mid y \sim \mathrm{N}(\bar{y}, \sigma^2/n)\). 6 번째 머신 평균 = \(\mu\). 7 번째 머신 평균도 = \(\mu\).

4.2.3 Hierarchical (§ 11.6 와 동일 구조)

\[ y_{ij} \mid \theta_j, \sigma \sim \mathrm{N}(\theta_j, \sigma^2), \quad \theta_j \mid \mu, \tau \sim \mathrm{N}(\mu, \tau^2) \]

Prior: \(p(\mu, \log \sigma, \tau) \propto 1\).

Gibbs: § 11.6 의 식 (11.9)~(11.17) 그대로 적용 (단 \(J = 6\)).

4.3 세 질문의 답

(i) 6 번째 머신 평균 \(\theta_6\):

• Separate: \(\bar{y}_{\cdot 6} \approx 86\), 분산 \(\sigma^2/5 \approx 80\). 넓은 구간.
• Pooled: \(\mu \approx 90\), 분산 \(\sigma^2/30 \approx 13\). 좁지만 그룹 6 의 특성 반영 안 됨.
• Hierarchical: shrinkage 형태 \(\widehat{\theta}_6 = \frac{\mu/\tau^2 + 5 \bar{y}_{\cdot 6}/\sigma^2}{1/\tau^2 + 5/\sigma^2}\). 양 극단 사이 절충.

(ii) 6 번째 머신 다음 측정 \(\tilde{y}_6\) 예측:

• Separate: \(\mathrm{N}(\theta_6, \sigma^2)\) on \(\theta_6\) 의 사후.
• Pooled: \(\mathrm{N}(\mu, \sigma^2)\) on \(\mu\) 의 사후.
• Hierarchical: \(\mathrm{N}(\theta_6, \sigma^2)\) on \(\theta_6\) (shrinkage 적용) 의 사후.

(iii) 관측 없는 7 번째 머신 평균 \(\theta_7\):

• Separate: 예측 불가. 사전 없음.
• Pooled: \(\mu\) (모든 머신 같음).
• Hierarchical: \(\theta_7 \sim \mathrm{N}(\mu, \tau^2)\), \((\mu, \tau)\) 사후에서 추출.

직관 — 세 모형이 주는 세 다른 세계관

Separate: “모든 머신은 완전 다르다” — 관측 있는 머신만 추론 가능. Pooled: “모든 머신은 같다” — 1 개 평균만. 새 머신 예측 가능하지만 그룹별 차이 무시. Hierarchical: “머신들은 분포에서 추출된 관련 있는 단위” — 새 머신 예측에 불확실성 계층으로 반영.

새 머신의 예측 분산:

• Pooled: \(\sigma^2\) (측정 오차만).
• Hierarchical: \(\tau^2 + \sigma^2\) (머신 간 분산 + 측정 오차).

실무에서는 Hierarchical 이 정직 — 새 머신이 분포의 어디에 있을지 모르는 불확실성을 분산에 포함.

4.4 수치 기대

Separate vs Pooled vs Hierarchical 비교:

• \(\tau\) 가 작으면 Hierarchical \(\to\) Pooled (모든 머신 유사).
• \(\tau\) 가 크면 Hierarchical \(\to\) Separate (머신 간 큰 차이).
• Table 11.4 의 데이터는 중간: \(\tau \approx 10{-}15\), \(\sigma \approx 12{-}15\). 부분 shrinkage.

7 번째 머신 예측 구간:

• Pooled: \([\mu - 2\sigma, \mu + 2\sigma]\) 대략 \([60, 120]\).
• Hierarchical: \([\mu - 2\sqrt{\tau^2+\sigma^2}, \ldots]\) 더 넓음.

5 문제 4 — 머신별 분산 계층화

문제. 문제 3 의 모형을 확장 — 머신별 분산 \(\sigma_j^2\) 를 Inv-\(\chi^2(\nu_0, \sigma_0^2)\) 로 계층화.

5.1 확장 모형

\[ y_{ij} \mid \theta_j, \sigma_j \sim \mathrm{N}(\theta_j, \sigma_j^2) \]

\[ \theta_j \mid \mu, \tau \sim \mathrm{N}(\mu, \tau^2), \quad \sigma_j^2 \sim \mathrm{Inv}\text{-}\chi^2(\nu_0, \sigma_0^2) \]

\(\nu_0\) 은 고정 (예: \(\nu_0 = 4\)), \(\sigma_0^2\) 는 미지 hyperparameter. Prior: \(p(\sigma_0^2) \propto 1\).

5.2 도전 — \(\sigma_0^2\) 의 조건부

\(\sigma_j^2 \sim \mathrm{Inv}\text{-}\chi^2(\nu_0, \sigma_0^2)\) 의 pdf:

\[ p(\sigma_j^2 \mid \nu_0, \sigma_0^2) \propto (\sigma_j^2)^{-\nu_0/2 - 1} \exp\!\left(-\frac{\nu_0 \sigma_0^2}{2 \sigma_j^2}\right) \]

\(\sigma_0^2\) 의 조건부 사후:

\[ p(\sigma_0^2 \mid \sigma_{1:J}^2, \nu_0) \propto \prod_j (\sigma_0^2)^{\nu_0/2} \exp\!\left(-\frac{\nu_0 \sigma_0^2}{2 \sigma_j^2}\right) = (\sigma_0^2)^{J \nu_0/2} \exp\!\left(-\frac{\nu_0 \sigma_0^2}{2} \sum_j \frac{1}{\sigma_j^2}\right) \]

Gamma 형태: \(\sigma_0^2 \mid \sigma_{1:J}^2 \sim \mathrm{Gamma}(J\nu_0/2 + 1, \frac{\nu_0}{2} \sum_j \sigma_j^{-2})\).

직접 추출 가능. Gibbs 적용.

5.3 Gibbs 갱신

각 반복:

1. \(\theta_j\): 식 (11.9)(11.10)(11.11) 에서 \(\sigma\) 를 \(\sigma_j\) 로 대체.
2. \(\mu\): 식 (11.12)(11.13).
3. \(\tau^2\): 식 (11.16)(11.17).
4. \(\sigma_j^2\): Inv-\(\chi^2(\nu_0 + n_j, \cdot)\) 공액 업데이트.
5. \(\sigma_0^2\): 위 유도한 Gamma 분포.

5.4 \(\nu_0\) 고정 이유

“데이터가 \(\nu_0\) 를 식별하기에 부족함”. \(\nu_0\) 와 \(\sigma_0^2\) 가 Inv-\(\chi^2\) 에서 교환 가능 한 구조 → prior 지배. 고정 (\(\nu_0 = 4, 8\) 등) 후 민감도 분석.

직관 — 왜 \(\sigma_0\) 가 아닌 \(\nu_0\) 를 고정하는가

Inv-\(\chi^2(\nu_0, \sigma_0^2)\) 에서 \(\nu_0\) 는 얼마나 많은 “가상 관측” 이 prior 에 있는가. 자유도가 10 이면 “10 개 관측에 해당하는 강한 prior”.

6 머신 × 5 측정 = 30 관측의 데이터로 \(\nu_0\) (prior 강도 자체) 와 \(\sigma_0\) (prior 스케일) 를 둘 다 추정하려 하면 두 parameters 가 confounded.

실무 전략: \(\nu_0\) 를 “약간 정보적” 수준 (4~8) 에 고정, \(\sigma_0\) 만 추정. \(\nu_0\) 값을 바꿔 가며 민감도 분석.

6 문제 5 — \(\widehat{\mathrm{var}}^+\) 불편성 증명

문제. 식 (11.3) \(\widehat{\mathrm{var}}^+\) 가 시작 분포 = 목표 분포 하에서 주변 사후 분산의 불편 추정량임을 증명. 힌트: 다른 체인의 시뮬레이션 간 제곱 차이의 반의 평균으로 표현.

6.1 풀이 — 쌍별 차이 표현

핵심 힌트: \(\widehat{\mathrm{var}}^+\) 를 서로 다른 체인의 관측 간 제곱 차이 의 평균으로 표현.

설정: \(m\) 체인, 각 \(n\) 관측, \(\psi_{ij}\). 체인 평균 \(\bar{\psi}_{\cdot j}\).

과분산 시작 = 목표 분포 가정: 모든 \(\psi_{ij}\) 가 independent draws from \(p(\psi \mid y)\).

이 가정 하에서:

\[ \mathbb{E}[\psi_{ij}] = \mu_\psi, \quad \mathrm{Var}[\psi_{ij}] = \mathrm{var}(\psi \mid y) \]

다른 체인 (\(j \ne j'\)) 독립:

\[ \mathbb{E}\!\left[\tfrac{1}{2}(\psi_{ij} - \psi_{ij'})^2\right] = \tfrac{1}{2} \cdot 2 \mathrm{var}(\psi \mid y) = \mathrm{var}(\psi \mid y) \]

6.2 \(B\) 와 \(W\) 의 기댓값

\(\bar{\psi}_{\cdot j} \sim \mathrm{N}(\mu_\psi, \mathrm{var}(\psi)/n)\) 이므로:

\[ \mathbb{E}[B] = \mathbb{E}\!\left[\frac{n}{m-1}\sum_j (\bar{\psi}_{\cdot j} - \bar{\psi}_{\cdot \cdot})^2\right] = n \cdot \frac{\mathrm{var}(\psi \mid y)}{n} = \mathrm{var}(\psi \mid y) \]

(표준 분산 추정의 불편성.)

\(s_j^2\) 의 불편성 (\(n - 1\) 자유도) 때문에:

\[ \mathbb{E}[W] = \mathrm{var}(\psi \mid y) \]

6.3 \(\widehat{\mathrm{var}}^+\) 계산

\[ \mathbb{E}[\widehat{\mathrm{var}}^+] = \frac{n-1}{n} \mathbb{E}[W] + \frac{1}{n} \mathbb{E}[B] = \frac{n-1}{n} \mathrm{var}(\psi) + \frac{1}{n} \mathrm{var}(\psi) = \mathrm{var}(\psi \mid y) \]

정확히 사후 분산. \(\square\)

6.4 \(n \to \infty\) 극한 — 무작위 시작

시작 분포가 목표 분포와 다르더라도 \(n \to \infty\) 에서 Markov chain 수렴 → \(\mathbb{E}[W] \to \mathrm{var}(\psi)\).

\(B\) 의 행동: 체인 평균 \(\bar{\psi}_{\cdot j} \to \mu_\psi\) a.s. → \(B \to 0\). 따라서:

\[ \widehat{\mathrm{var}}^+ = \frac{n-1}{n} W + \frac{1}{n} B \to W \to \mathrm{var}(\psi) \]

두 경우 모두 점근 불편. 유한 \(n\) + 과분산 시작 → 과대 추정 (upper bound).

직관 — \(\widehat{\mathrm{var}}^+\) 의 세 얼굴

1. 이상적 (시작 = 목표): 정확한 사후 분산.
2. 수렴 후: 점근 불편.
3. 수렴 전 (과분산): 과대 추정 (upper bound).

\(\widehat{R} = \sqrt{\widehat{\mathrm{var}}^+/W}\) 의 의미 재해석: \(W\) 가 과소, \(\widehat{\mathrm{var}}^+\) 가 과대 → 비율 \(> 1\). 수렴 시 둘이 같아져 1 에 수렴.

이 증명이 보여주는 것: \(\widehat{R}\) 은 유한 \(n\) 에서 “현재 추정된 분산이 얼마나 참값보다 큰가” 의 보수적 척도.

7 문제 6 — Effective Sample Size 유도

문제. (a) 식 (11.5) 점근 분산 공식 유도. (b) 어떤 예제로 \(\widehat{n}_{\mathrm{eff}}\) 시간 추적.

7.1 (a) 점근 분산 유도

설정: 정상 Markov chain \(\{\psi^t\}\). 평균 \(\bar{\psi}_n = n^{-1} \sum_{t=1}^n \psi^t\). 자기상관 \(\rho_t = \mathrm{Cor}(\psi^s, \psi^{s+t})\).

분산:

\[ \mathrm{Var}[\bar{\psi}_n] = \frac{1}{n^2} \sum_{s, t} \mathrm{Cov}(\psi^s, \psi^t) = \frac{\mathrm{Var}[\psi]}{n^2} \sum_{s, t} \rho_{|s-t|} \]

대각 (s = t): \(n\) 개, 각 \(1\). 대각 외: 각 \((s, t)\), \(s \ne t\) — \(\rho_{|s-t|}\).

\[ \sum_{s, t} \rho_{|s-t|} = n + 2 \sum_{t=1}^{n-1} (n - t) \rho_t \]

\(n\) 으로 나눔:

\[ \frac{1}{n} \sum_{s, t} \rho_{|s-t|} = 1 + 2 \sum_{t=1}^{n-1} \left(1 - \frac{t}{n}\right) \rho_t \]

\(n \to \infty\) (ergodic): \((1 - t/n) \to 1\), 합이 \(\sum_{t=1}^\infty \rho_t\) 로 수렴:

\[ \lim_{n \to \infty} n \cdot \mathrm{Var}[\bar{\psi}_n] = \mathrm{Var}[\psi] \left(1 + 2 \sum_{t=1}^\infty \rho_t\right) \]

\(m\) 체인 결합:

\[ \lim_{n \to \infty} mn \cdot \mathrm{Var}[\bar{\psi}_{\cdot \cdot}] = \mathrm{Var}[\psi] \left(1 + 2 \sum_{t=1}^\infty \rho_t\right) \tag{11.5} \]

(각 체인 독립이므로 \(m\) 배 효율.)

7.2 (b) \(\widehat{n}_{\mathrm{eff}}\) 시간 추적

이상적 행동: \(\widehat{n}_{\mathrm{eff}}\) 가 시간 \(t\) (누적 반복 수) 에 비례 증가.

\[ \widehat{n}_{\mathrm{eff}} \approx \frac{mn}{\tau_{\mathrm{ac}}} \]

\(\tau_{\mathrm{ac}}\) 가 수렴 후 상수. 따라서 \(n\) 증가 시 \(\widehat{n}_{\mathrm{eff}}\) 선형 증가.

실무 확인: \(\widehat{n}_{\mathrm{eff}}(t)\) 를 여러 \(t\) 에 대해 계산해 plot.

# 가정: chains 배열 (m, n_total, D) 있음
n_steps = [100, 200, 500, 1000, 2000, 5000]
for n in n_steps:
    subset = chains[:, :n, :]
    n_eff = ess_from_chains(subset)  # 문제 3 함수
    print(f"n = {n:5d}: n_eff for mu = {n_eff['mu']:.0f}")

기대 결과: \(\widehat{n}_{\mathrm{eff}}\) 가 \(n\) 에 비례 증가. 비례 깨짐 시:

• 수렴 전이면 \(\widehat{n}_{\mathrm{eff}}\) 가 흔들림.
• 수렴 후에도 깨지면 \(\tau_{\mathrm{ac}}\) 추정 불안정 — 체인 더 길게.

직관 — \(\widehat{n}_{\mathrm{eff}}\) 의 안정성 지표

수렴 후 \(\widehat{n}_{\mathrm{eff}}\) 의 기울기 가 일정해야. 기울기가 변하면:

• 증가 기울기 감소: \(\tau_{\mathrm{ac}}\) 가 실제보다 작게 추정되다가 교정 — 초기 부족.
• 증가 기울기 증가: 체인이 아직 탐색 중 — 수렴 미완료.
• 기울기 음수: 버그 또는 정지 규칙 오류.

이 시간 추적이 \(\widehat{R}\) 이 잡지 못하는 미세한 문제까지 감지. Stan 의 ess_bulk diagnostic 이 이 원리.

8 문제 7 — Stratified Survey 재현

문제. § 8.3 의 16 strata 계층 다항 모형을 재현.

8.1 모형 (§ 8.3 복습)

\(y_{ij}\) = strata \(j\) 의 Bush/Dukakis/무의견 카운트.

\[ y_{ij} \sim \mathrm{Multinomial}(n_j; \theta_{1j}, \theta_{2j}, \theta_{3j}) \]

Logit 변환:

\[ \alpha_{1j} = \frac{\theta_{1j}}{\theta_{1j} + \theta_{2j}}, \quad \alpha_{2j} = 1 - \theta_{3j} \]

\[ \beta_{1j} = \mathrm{logit}(\alpha_{1j}), \quad \beta_{2j} = \mathrm{logit}(\alpha_{2j}) \]

계층:

\[ (\beta_{1j}, \beta_{2j}) \sim \mathrm{N}_2(\mu, \Sigma) \]

\(\mu = (\mu_1, \mu_2)\), \(\Sigma\) 는 \(2 \times 2\) 공분산 (분산 \(\tau_1^2, \tau_2^2\), 상관 \(\rho\)).

8.2 (a) Nonhierarchical

각 strata 독립 Dirichlet(1,1,1) prior → 독립 Dirichlet 사후. 격자 또는 직접 추출. 식 (8.7) 의 가중 차이 \(\sum_j (N_j/N)(\theta_{1j} - \theta_{2j})\) 계산.

8.3 (b) Hierarchical + Metropolis

37 차원 모수 (\(\beta\) 32 + hyperparameter 5). 전 결합 Metropolis 는 고차원으로 비효율.

전략: Metropolis-within-Gibbs.

1. \(\beta_j \mid \mu, \Sigma, y\): Multinomial-logit — 비공액. Metropolis.
2. \(\mu \mid \beta, \Sigma\): 정규 공액 (32 관측, 2 차원 정규). 직접 추출.
3. \(\Sigma \mid \beta, \mu\): Wishart 공액 — 직접 추출.

\(\beta_j\) 블록 Metropolis: 2 차원, \(J\) 블록 (16 strata). 각 블록 수용률 25~40% 목표.

8.4 Figure 8.1b 와 비교

예상 결과: hierarchical 이 shrinkage 효과로 nonhierarchical 보다 좁은 격차 분포. 중앙 \(\approx 0.10\), 더 집중.

9 Ch.11 심화 시리즈 — 4 포스트 논리 지도

Ch.11 의 완결된 경로.

포스트	역할	핵심 메시지
02-11-0 Overview	지도	MCMC 혁명, Gibbs/MH/진단 축
02-11-1 §11.1~11.3	이론	Detailed balance 증명, Gibbs = MH 특수
02-11-2 §11.4~11.6	실무	\(\widehat{R}\), ESS, coagulation Gibbs
02-11-3 §11.8 (본편)	연습·결산	7 문제 풀이 + 4 포스트 지도

Ch.11 의 한 문장:

MCMC 는 Markov 전이 \(T\) 를 설계해 목표 분포 \(p(\theta \mid y)\) 에서 샘플링하는 일반 전략이며, Gibbs (조건부 공액) 과 Metropolis-Hastings (임의 사후) 의 두 축으로 구현된다. 수렴은 detailed balance 로 보장되고, 실무 수렴은 \(\widehat{R}\), \(n_{\mathrm{eff}}\) 로 확인한다. 계층 정규 모형은 Gibbs 의 전형이며, 복잡 모형은 Ch.12 의 HMC 로 연결된다.

10 전체 코드 — 문제 3 의 Hierarchical Gibbs 구현

문제 3 의 세 모형 중 hierarchical 모형을 Gibbs 로 구현.

10.1 데이터

import numpy as np
rng = np.random.default_rng(113)

machines = [
    [83, 92, 92, 46, 67],
    [117, 109, 114, 104, 87],
    [101, 93, 92, 86, 67],
    [105, 119, 116, 102, 116],
    [79, 97, 103, 79, 92],
    [57, 92, 104, 77, 100],
]
J = len(machines)
y = [np.array(m, dtype=float) for m in machines]
n_j = np.array([len(m) for m in y])
n = int(n_j.sum())
ybar = np.array([m.mean() for m in y])

10.2 Gibbs (§ 11.6 와 같은 구조)

def gibbs_machines(n_iter, seed):
    rng = np.random.default_rng(seed)
    theta = np.array([rng.choice(m) for m in y], dtype=float)
    mu = theta.mean()
    sigma2 = 100.0
    tau2 = 100.0

    chain = np.zeros((n_iter, J + 3))

    for t in range(n_iter):
        V_theta = 1.0 / (1.0/tau2 + n_j/sigma2)
        theta_hat = V_theta * (mu/tau2 + n_j*ybar/sigma2)
        theta = rng.normal(theta_hat, np.sqrt(V_theta))

        mu = rng.normal(theta.mean(), np.sqrt(tau2/J))

        ss = sum(((g - theta[j])**2).sum() for j, g in enumerate(y))
        sigma2 = n * (ss/n) / rng.chisquare(n)

        tau_hat2 = ((theta - mu)**2).sum() / (J - 1)
        tau2 = (J - 1) * tau_hat2 / rng.chisquare(J - 1)

        chain[t] = np.concatenate([theta, [mu, np.sqrt(sigma2), np.sqrt(tau2)]])

    return chain

chains = np.array([gibbs_machines(1000, seed=i) for i in range(10)])
saved = chains[:, 500:, :].reshape(-1, J + 3)

10.3 세 질문 답변

# (i) 6 번째 머신 평균
theta_6 = saved[:, 5]
print(f"theta_6 사후 평균: {theta_6.mean():.1f}, 95% 구간: [{np.quantile(theta_6, 0.025):.1f}, {np.quantile(theta_6, 0.975):.1f}]")

# (ii) 6 번째 머신 다음 측정
sigma_samples = saved[:, J+1]
y_tilde_6 = saved[:, 5] + sigma_samples * rng.normal(0, 1, size=len(saved))
print(f"y_tilde_6 사후 평균: {y_tilde_6.mean():.1f}, 95% 구간: [{np.quantile(y_tilde_6, 0.025):.1f}, {np.quantile(y_tilde_6, 0.975):.1f}]")

# (iii) 7 번째 머신 평균
mu_samples = saved[:, J]
tau_samples = saved[:, J+2]
theta_7 = mu_samples + tau_samples * rng.normal(0, 1, size=len(saved))
print(f"theta_7 사후 평균: {theta_7.mean():.1f}, 95% 구간: [{np.quantile(theta_7, 0.025):.1f}, {np.quantile(theta_7, 0.975):.1f}]")

기대 결과: \(\theta_6 \approx 86\), \(\tilde{y}_6\) 구간 넓음 (\(\sigma^2\) 추가), \(\theta_7\) 구간 더 넓음 (\(\tau^2\) 추가).

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

MCMC 실무의 핵심 15 항목.

설계

1. 조건부 공액 확인 → Gibbs 우선.
2. 비공액 → Metropolis 또는 Metropolis-within-Gibbs.
3. 블록 크기 조정 — 상관 변수 묶기.
4. 과분산 시작점 — 조잡한 추정 + 큰 SD.

구현

5. 로그 밀도로 계산 — overflow 방지.
6. 제안 공분산 \(\Sigma\) — Laplace 근사 또는 이전 run.
7. \(c = 2.4/\sqrt{d}\) — Roberts-Gelman-Gilks 최적.
8. 수용률 20~40% 목표 — 아니면 \(c\) 조정.

진단

9. \(m \ge 4\) 체인 — split 후 \(m \ge 8\).
10. Warm-up 50% — 보수적 기본.
11. \(\widehat{R} < 1.1\) — 모든 추정량에서.
12. \(\widehat{n}_{\mathrm{eff}} \ge 5m\) — 기본 기준.
13. \(\widehat{n}_{\mathrm{eff}}(t)\) 선형 증가 — 기울기 일정성.

문제 해결

14. 느리면 재매개변수화 — 좌표 축 정렬 또는 non-centered.
15. 현대 도구 우선 — Stan, PyMC, NumPyro 가 진단 자동화.

12 관련 주제

선행 지식

• Ch.11 Overview (02-11-0) — MCMC 지도
• § 11.1~11.3 (02-11-1) — Gibbs·MH 이론
• § 11.4~11.6 (02-11-2) — 진단 + 계층 정규

후속 주제

• Ch.12 Efficient MCMC — HMC, NUTS, Stan
• Ch.13 Variational Inference — MCMC 대안
• Ch.15 Hierarchical Linear Models — MCMC 본격 응용

관련 개념

• Gelman & Rubin (1992b) — \(\widehat{R}\) 원저
• Brooks & Gelman (1998) — 다변량 \(\widehat{R}\)
• Geyer (1992) — ESS 부분합 정지
• Vehtari et al. (2021) — Rank-normalized split-\(\widehat{R}\) 현대 표준
• Cook, Gelman, Rubin (2006) — Fake-data simulation 검증
• Brooks, Gelman, Jones, Meng (2011), Handbook of MCMC — 종합 참고서