Ch.5.4~5.6 — 정규 계층 모형·8 학교·메타분석 심화

Gelman BDA Ch.5.4~5.6 상세 — \(p(\tau \mid y)\) 유도·SAT 코칭·베타차단제 22 연구

Gelman et al. Bayesian Data Analysis (3rd ed., 2013) Ch.5 중반 세 절을 교재 원문 수준으로 심화한다. § 5.4 정규 교환가능 모델의 공동 사후 (5.16)· 조건부 \(\theta_j \mid \mu, \tau, y\) 의 정밀도 가중 평균 (5.17)·\(\mu\) 주변화 후 \(p(\tau \mid y)\) 유도 (5.21)·\(p(\mu \mid \tau, y)\) 의 정규성·계산 전략· 빈도주의 \(\hat\tau^2\) 가 음수 되는 결함, § 5.5 8 학교 SAT 코칭 실험 — no pooling vs complete pooling 의 문제·\(p(\tau \mid y)\) 와 shrinkage 패턴· “학교 A 가 최고” 라는 결론이 사라지는 이유·다중 비교의 자동 해결, § 5.6 베타차단제 22 임상 시험 메타분석·\(\mu\) 의 95% 구간 odds ratio [0.69, 0.90]·complete pooling [0.70, 0.85] 의 과신·새 연구 예측 분포의 넓은 꼬리까지 수식·직관·코드로 완결.

Statistics
Bayesian
저자

Kwangmin Kim

공개

2026년 04월 20일

1 이 포스트의 위치 — Ch.5 심화의 두 번째 조각

§ 5.1~5.3 심화 가 Beta-Binomial 계층 모형으로 계층 구조를 해부했다면, 이 포스트는 정규 계층 모형에서 같은 구조를 반복 하고 8 학교 · 메타분석 에서 실무 적용을 본다. 8 학교는 Gelman BDA 전체를 통틀어 가장 자주 인용되는 예제.

§ 5.4~5.6 의 한 줄 요약

“정규 교환가능 모델에서 \(\mu\) 를 해석적으로 적분 소거하면 \(\tau\) 만의 1 차원 주변 사후 \(p(\tau \mid y)\) 가 남고, 이 분포의 모양이 ‘그룹 간 변동이 얼마나 큰가’ 를 데이터가 말한다. 8 학교와 베타차단제 22 연구에서 \(\tau\) 의 사후 분포가 complete pooling (\(\tau = 0\)) 의 과신을 드러낸다.”

Ch.5 의 상징 예제 (8 schools) 가 여기서 완전 전개 (Gelman et al., 2013, Ch.5.4~5.6).

2 § 5.4 정규 교환가능 평균 모델

2.1 데이터 구조

\(J\) 개 독립 실험, 실험 \(j\) 에서 \(n_j\) 개 관측치 \(y_{ij}\) (각각 알려진 오차 분산 \(\sigma^2\)).

\[ y_{ij} \mid \theta_j \sim N(\theta_j, \sigma^2), \quad i = 1, \ldots, n_j; \ j = 1, \ldots, J \tag{5.11} \]

충분통계량 — 각 그룹의 표본 평균.

\[ \bar{y}_{\cdot j} = \frac{1}{n_j}\sum_i y_{ij}, \quad \bar{y}_{\cdot j} \mid \theta_j \sim N(\theta_j, \sigma_j^2), \quad \sigma_j^2 = \sigma^2/n_j \tag{5.12} \]

2.2 계층 사전

\(\theta_j\) 가 교환가능 정규 모집단에서 추출.

\[ \theta_j \mid \mu, \tau \sim N(\mu, \tau^2), \quad j = 1, \ldots, J \]

공동 사전.

\[ p(\theta_1, \ldots, \theta_J) = \int \prod_{j=1}^J N(\theta_j \mid \mu, \tau^2) \, p(\mu, \tau) \, d(\mu, \tau) \tag{5.14} \]

2.3 Hyperprior

\(\mu\) 에 대해 균등 (데이터가 \(\mu\) 에 대해 정보 많음).

\[ p(\mu, \tau) \propto p(\tau) \tag{5.15} \]

\(\tau\) 의 사전 은 § 5.7 에서 상세. 지금은 균등 \(p(\tau) \propto 1\) 사용.

2.4 공동 사후 (5.16)

가능도와 사전의 곱.

\[ p(\theta, \mu, \tau \mid y) \propto p(\mu, \tau) \prod_{j=1}^J N(\theta_j \mid \mu, \tau^2) \prod_{j=1}^J N(\bar{y}_{\cdot j} \mid \theta_j, \sigma_j^2) \tag{5.16} \]

2.5 조건부 사후 \(p(\theta \mid \mu, \tau, y)\)

\((\mu, \tau)\) 고정 하에서 \(\theta_j\) 들이 독립 이고 각각 정규-정규 켤레.

\[ \theta_j \mid \mu, \tau, y \sim N(\hat\theta_j, V_j) \]

\[ \hat\theta_j = \frac{\frac{1}{\sigma_j^2}\bar{y}_{\cdot j} + \frac{1}{\tau^2}\mu}{\frac{1}{\sigma_j^2} + \frac{1}{\tau^2}}, \quad V_j = \frac{1}{\frac{1}{\sigma_j^2} + \frac{1}{\tau^2}} \tag{5.17} \]

정밀도 가중 평균 — Ch.2.5 의 정규-정규 공식. 관측 정밀도 \(1/\sigma_j^2\) 와 모집단 정밀도 \(1/\tau^2\) 의 가중.

2.6 주변 사후 \(p(\mu, \tau \mid y)\)

\(\theta\) 를 적분 소거. 정규 모델에서 \(\bar{y}_{\cdot j}\) 의 주변 분포 가 닫힌 형태.

\[ \bar{y}_{\cdot j} \mid \mu, \tau \sim N(\mu, \sigma_j^2 + \tau^2) \]

유도 — 반복 기댓값·분산 공식 (Ch.1.8).

  • \(E(\bar{y}_{\cdot j} \mid \mu, \tau) = E(E(\bar{y}_{\cdot j} \mid \theta_j) \mid \mu, \tau) = E(\theta_j \mid \mu, \tau) = \mu\)
  • \(\text{Var}(\bar{y}_{\cdot j} \mid \mu, \tau) = E(\text{Var}(\bar{y}_{\cdot j} \mid \theta_j)) + \text{Var}(E(\bar{y}_{\cdot j} \mid \theta_j)) = \sigma_j^2 + \tau^2\)

그룹 내 분산 + 그룹 간 분산. 이것이 분산 분해의 베이즈 버전.

주변 사후.

\[ p(\mu, \tau \mid y) \propto p(\mu, \tau) \prod_{j=1}^J N(\bar{y}_{\cdot j} \mid \mu, \sigma_j^2 + \tau^2) \tag{5.18} \]

2.7 \(\mu\) 를 한 번 더 적분

\(p(\mu, \tau \mid y) = p(\mu \mid \tau, y) \, p(\tau \mid y)\) 로 분해.

\(\mu \mid \tau, y\) — \(\tau\) 가 알려졌다면 (5.18) 의 \(\mu\) 에 대한 함수가 이차 지수 → 정규.

\[ \mu \mid \tau, y \sim N(\hat\mu, V_\mu) \]

\[ \hat\mu = \frac{\sum_j \frac{1}{\sigma_j^2 + \tau^2} \bar{y}_{\cdot j}}{\sum_j \frac{1}{\sigma_j^2 + \tau^2}}, \quad V_\mu^{-1} = \sum_{j=1}^J \frac{1}{\sigma_j^2 + \tau^2} \tag{5.20} \]

\(\mu\) 는 \(\bar{y}_{\cdot j}\) 들의 정밀도 가중 평균. 정밀도는 \(1/(\sigma_j^2 + \tau^2)\).

2.8 \(\tau\) 의 주변 사후

\(p(\mu, \tau \mid y) / p(\mu \mid \tau, y)\) 계산. 이 항등식이 모든 \(\mu\) 에서 성립 → \(\mu = \hat\mu\) 대입이 단순.

\[ p(\tau \mid y) \propto p(\tau) V_\mu^{1/2} \prod_{j=1}^J (\sigma_j^2 + \tau^2)^{-1/2} \exp\left(-\frac{(\bar{y}_{\cdot j} - \hat\mu)^2}{2(\sigma_j^2 + \tau^2)}\right) \tag{5.21} \]

\(\tau\) 의 함수 로 닫힌 형태지만 복잡. 1 차원 수치 계산 (격자).

직관 — \(p(\tau \mid y)\) 의 의미

\(\tau\) 는 “그룹 간 변동의 크기”. 데이터가 \(\tau\) 에 대해 주는 정보 —

  • \(\tau \to 0\) (complete pooling): 모든 \(\bar{y}_{\cdot j}\) 가 \(\mu\) 주변에 분산 \(\sigma_j^2\) 으로 분포. 데이터의 \(\bar{y}_{\cdot j}\) 들이 서로 가까우면 \(\tau\) 작음.
  • \(\tau \to \infty\) (no pooling): 각 \(\bar{y}_{\cdot j}\) 의 주변 분산이 \(\sigma_j^2 + \tau^2 \to \infty\). 데이터의 \(\bar{y}_{\cdot j}\) 들이 매우 달라야 \(\tau\) 큼.
  • 중간값 — 그룹 간 차이와 그룹 내 측정 오차의 균형.

\(p(\tau \mid y)\) 의 모양이 데이터가 “얼마나 shared” 인지 를 직접 말한다. 8 학교에서는 \(\tau\) 가 0 근처가 아닐 확률이 있지만 크지도 않다 — “약간의 공유”.

2.9 계산 전략

3 단계 순차 샘플링.

  1. \(\tau^{(s)} \sim p(\tau \mid y)\) — 격자 + 역 CDF (Ch.1.9)
  2. \(\mu^{(s)} \mid \tau^{(s)}, y \sim N(\hat\mu, V_\mu)\) — (5.20) 에서 정규 샘플
  3. \(\theta_j^{(s)} \mid \mu^{(s)}, \tau^{(s)}, y \sim N(\hat\theta_j, V_j)\) — 각 \(j\) 독립 (5.17)

Ch.3.1 의 조건부-주변 분해 전략의 3 단계 버전.

2.10 사후 예측 분포

시나리오 A — 기존 그룹 \(j\) 의 새 관측.

\(\theta_j\) 사후에서 \(\tilde{y} \sim N(\theta_j^{(s)}, \sigma^2)\) 추출.

시나리오 B — 새 그룹 \(\tilde{j}\) 의 관측.

  1. \((\mu, \tau)\) 사후 추출
  2. \(\tilde\theta \sim N(\mu, \tau^2)\) — 새 그룹 모수
  3. \(\tilde{y} \sim N(\tilde\theta, \sigma^2)\) — 새 관측

두 시나리오의 차이 — B 가 \(\tau\) 변동까지 포함해 예측 분포가 훨씬 넓다.

2.11 빈도주의 점 추정의 결함

ANOVA 기반 비편향 추정.

\[ \hat\mu = \bar{y}, \quad \hat\tau^2 = (\text{MS}_B - \text{MS}_W)/n \tag{5.22} \]

\(\text{MS}_B\) = between-group, \(\text{MS}_W\) = within-group mean squares.

치명적 결함 — \(\hat\tau^2\) 가 음수일 수 있음.

“\(\text{MS}_W > \text{MS}_B\) 이면 \(\hat\tau^2 < 0\). 문제를 피하려 \(\hat\tau^2 = 0\) 으로 두는 관행이 있지만, 이는 너무 강한 주장 — ‘모든 그룹 평균이 정확히 같다’ 는 것. 실제로는 표본 크기가 \(\tau^2\) 를 0 과 구분할 수 없을 만큼 작다 는 뜻인데.” (교재)

베이즈의 우월성 — \(\tau^2\) 의 사후 분포 전체 를 다룸. 0 에 대한 강한 주장 없이 데이터가 말하는 만큼만.

3 § 5.5 8 학교 SAT 코칭 실험

3.1 연구 배경

ETS (Educational Testing Service) 가 8 개 고등학교 의 SAT-V 코칭 프로그램 효과 측정. 각 학교 독립 무작위 실험.

데이터 (표 5.2).

학교 \(y_j\) (추정 효과) \(\sigma_j\) (표준오차)
A 28 15
B 8 10
C -3 16
D 7 11
E -1 9
F 1 11
G 18 10
H 12 18

단위 — SAT-V 점수 증가 (원 점수 200~800, 평균 약 500, SD 약 100). 8 점 ≈ 1 문제 더 맞음.

\(y_j\) 는 ANCOVA 조정된 추정 (PSAT-M, PSAT-V 공변량 사용), \(\sigma_j^2\) 는 빈도주의 표준오차. 각 학교 30 명 이상이라 정규 근사 합리적.

3.2 사전 고려 — 학교 간 선험적 구분 없음

“어떤 프로그램이 다른 것보다 효과적이라거나, 어떤 것이 다른 것보다 유사하다는 선험적 이유는 없었다.” (교재)

→ 교환가능 모델 적합.

3.3 Nonhierarchical Approach 1 — Separate Estimates

각 학교 독립 추정.

문제 — 표 5.2 만 보면 “A (28) 와 G (18) 는 효과 큰 프로그램” 처럼 보인다. 그러나 표준오차 고려 시 모든 학교의 95% 구간이 겹친다.

\[ y_A \pm 1.96 \sigma_A = 28 \pm 29.4 = [-1.4, 57.4] \]

\[ y_C \pm 1.96 \sigma_C = -3 \pm 31.4 = [-34.4, 28.4] \]

두 구간이 크게 중첩 — “A 가 C 보다 낫다” 는 주장이 통계적으로 확립되지 않음.

3.4 Nonhierarchical Approach 2 — Complete Pooling

“모든 학교 효과가 같다” 가정 → pooled estimate.

\[ \bar{y}_{\cdot\cdot} = \frac{\sum_j y_j/\sigma_j^2}{\sum_j 1/\sigma_j^2} = 7.7 \]

\[ \text{Var} = \left(\sum_j 1/\sigma_j^2\right)^{-1} = 16.6 \]

공통 효과 7.7 점, SE = 4.1. 95% 구간 \([-0.2, 15.6]\) — 0 을 포함 (효과 없음 가능).

3.5 두 접근의 문제

Separate — 각 학교 독립 추정의 극단값 (A: 28) 을 그대로 신뢰. 과도한 “세분”.

Complete pooling — 학교 간 변동 완전 무시. 과도한 “통합”.

합리적 답 — 둘 사이 어딘가 부분 풀링. 계층 모형이 정확히 이를 수행.

3.6 계층 모형

\[ y_j \mid \theta_j \sim N(\theta_j, \sigma_j^2), \quad \theta_j \mid \mu, \tau \sim N(\mu, \tau^2) \]

Hyperprior \(p(\mu, \tau) \propto 1\) (균등).

3.7 \(p(\tau \mid y)\) 의 모양

식 (5.21) 에 대입. \(\tau\) 의 격자 계산.

관찰.

  • \(\tau\) 의 주변 사후가 0 근처에 양의 밀도 를 갖지만 \(\tau = 0\) 에서 모드가 아니다
  • 모드는 약 \(\tau \approx 5\) 근처
  • 95% 구간이 약 \([0.1, 20]\) — 넓은 불확실성
  • \(\tau > 0\) 가 많이 지지되지만 \(\tau = 0\) 도 배제 안 됨

3.8 학교별 사후

각 \(\tau\) 에서 조건부 \(\theta_j\) 의 정규 샘플 → 각 학교의 사후 평균과 SD.

학교 \(y_j\) \(\sigma_j\) 사후 평균 \(E(\theta_j \mid y)\) 사후 SD
A 28 15 약 11 약 9
B 8 10 약 8 약 6
C -3 16 약 6 약 8
D 7 11 약 7 약 7
E -1 9 약 5 약 6
F 1 11 약 6 약 7
G 18 10 약 10 약 7
H 12 18 약 8 약 8

극단값 shrinkage.

  • 학교 A (원 28) → 사후 평균 11 — 크게 감소
  • 학교 C (원 -3) → 사후 평균 6 — 크게 증가

중앙값 근처는 거의 변화 없음 (B, D).

3.9 “학교 A 가 최고” 결론의 해체

Separate 관점 — A (28) 가 G (18) 나 H (12) 보다 명백히 큼.

계층 관점 — 사후에서.

  • \(\Pr(\theta_A > \theta_G \mid y) \approx 0.55\) — 거의 50:50
  • \(\Pr(\theta_A > \theta_C \mid y) \approx 0.66\) — 약한 증거

즉 “A 가 최고” 라는 결론이 사라진다. 원 관측 \(y_A = 28\) 이 상당 부분 운 일 가능성.

직관 — 8 학교의 교훈

\(y_A = 28\) 은 거품일 수 있다. 계층 모형이 말하는 진실.

  1. 진짜 효과 \(\theta_A\) 는 약 11 근처 (전체 평균 8 쪽으로 shrunk)
  2. 학교 A 와 C 의 실제 차이는 거의 구분 불가
  3. 모든 학교의 진짜 효과가 약 5~12 범위

정책 함의 — “학교 A 의 코칭을 전국 확산” 은 성급한 결정. 실제 기대 효과는 약 7~10 점 수준. 이것이 빈도주의 표준편차만 보면 놓치는 베이즈의 실용적 가치.

3.10 다중 비교의 자동 해결

28 개 짝 비교 — \(\theta_i - \theta_j\) 각각.

빈도주의 — Bonferroni 보정 필요. 유의 수준 \(0.05/28 \approx 0.002\). 또는 Tukey HSD.

베이즈 — 그냥 사후 분포에서 각 쌍 차이의 95% 구간 계산.

결과 — 28 개 쌍 모두 0 을 포함하는 95% 구간. 즉 “학교 간 차이가 통계적으로 확립되지 않았다”.

\(\tau\) 의 작은 사후 (\(\tau \approx 5\)) 가 자동으로 모든 차이를 shrink — 다중 비교 보정이 내장.

3.11 Sensitivity to Prior

교재의 민감도 분석 — \(p(\tau) \propto 1/\tau\) (스케일 비정보) vs 균등. 결과 — 결론이 사전에 민감, 특히 작은 \(\tau\) 영역. Gelman (2006a) 의 half-Cauchy 권장 (§ 5.7).

4 § 5.6 메타분석 — 베타차단제 22 연구

4.1 메타분석의 목적

메타분석 (meta-analysis) — 여러 독립 연구 결과를 통합하여 공통 효과 추정. 의학 연구에서 특히 중요.

4.2 데이터 — 베타차단제 심근경색 임상 시험

22 개 임상 시험 — 베타차단제 (beta-blocker) 가 심근경색 후 사망률에 미치는 효과.

각 시험 \(j\) 에서.

  • \(r_{0j}\) — 대조군 사망자 수, \(n_{0j}\) — 대조군 크기
  • \(r_{1j}\) — 처치군 사망자 수, \(n_{1j}\) — 처치군 크기

요약 통계 — 각 시험의 로그 오즈비 (log odds ratio).

\[ y_j = \log\left(\frac{r_{1j}(n_{0j} - r_{0j})}{r_{0j}(n_{1j} - r_{1j})}\right) \]

분산.

\[ \sigma_j^2 = \frac{1}{r_{1j}} + \frac{1}{n_{1j} - r_{1j}} + \frac{1}{r_{0j}} + \frac{1}{n_{0j} - r_{0j}} \]

22 시험의 \(y_j, \sigma_j^2\) 가 표 5.4 의 입력.

4.3 계층 모형 적용

정규 근사 (대부분 시험에서 표본 크기 수백 명).

\[ y_j \mid \theta_j \sim N(\theta_j, \sigma_j^2), \quad \theta_j \mid \mu, \tau \sim N(\mu, \tau^2) \]

\(\mu\) — 평균 로그 오즈비, \(\tau\) — 연구 간 이질성.

4.4 \(p(\tau \mid y)\) 결과

“\(\tau\) 의 주변 사후 밀도가 0 이 아닌 값에서 최고 지만, 0 근처 값도 분명히 plausible — 0 의 사후 밀도가 모드의 약 75% 수준.” (교재)

즉 \(\tau \approx 0\) 가능성 무시 못 함 — connective tissue 가 있지만 많지는 않음.

4.5 사후 요약 (표 5.5)

추정량 2.5% 25% 중앙 75% 97.5%
\(\mu\) (평균) -0.37 -0.29 -0.25 -0.20 -0.11
\(\tau\) (SD) 0.02 0.08 0.13 0.18 0.31
\(\tilde\theta_j\) (예측) -0.58 -0.34 -0.25 -0.17 0.11

해석.

  • \(\mu\) 의 95% 구간 \([-0.37, -0.11]\), 오즈비 스케일 [0.69, 0.90] — 베타차단제가 사망 오즈를 10~31% 감소
  • \(\tau\) 의 95% 구간 \([0.02, 0.31]\) — 연구 간 이질성 존재
  • \(\tilde\theta_j\) — 새 연구의 예측 효과, 95% 구간 \([-0.58, 0.11]\) → 오즈비 [0.56, 1.12] — 10% 확률로 새 연구에서 해로울 수도

4.6 Complete Pooling 과의 비교

Complete pooling — \(\tau = 0\) 가정.

\(\mu\) 의 95% 구간 → 오즈비 [0.70, 0.85]. 계층의 [0.69, 0.90] 보다 좁다.

“원 출판 논문은 이 [0.70, 0.85] 가 unusually narrow range of uncertainty 라고 논평했다. 계층 베이즈 분석이 밝혀낸 것 — 부적절한 모델 사용. 모든 연구가 동일하다는 주장이었다.” (교재)

4.7 핵심 교훈

Complete pooling 의 과신 — \(\tau = 0\) 가정이 실제로 존재하는 연구 간 이질성 을 무시. 신뢰 구간이 너무 좁다.

Fixed-effect meta-analysis (= complete pooling) vs Random-effects meta-analysis (= 계층) 의 선택이 결과를 크게 바꾼다.

직관 — 메타분석의 3 가지 질문
  1. \(\mu\) (평균 효과) — “전체 연구에서의 평균 효과는?” 정책 결정용
  2. \(\theta_j\) (개별 연구 효과) — “이 특정 연구의 효과는?” 재분석용
  3. \(\tilde\theta\) (새 연구 예측) — “비슷한 새 환자 집단에서의 예상 효과는?” 미래 적용

세 번째가 가장 보수적 — \(\mu\) 의 불확실성 + \(\tau\) 의 이질성 + 예측 잡음 모두 포함. 베타차단제에서 10% 확률로 해로울 수 있다 는 결과가 여기서 나옴.

빈도주의 메타분석은 보통 1 번만 다룬다. 베이즈가 세 질문 모두 자연스럽게 처리.

4.8 연구 간 이질성의 추가 해석

연구 간 이질성의 이유 — 환자 집단 차이, 약물 용량, 추적 기간, 병원 특성 등. \(\tau > 0\) 은 “연구들이 완벽히 동일한 조건이 아니다” 를 반영.

\(\tau\) 의 사후 의미 — “새 병원에서 이 약을 쓸 때 기대할 수 있는 효과의 변동성”. 임상 적용에서 중요.

5 세 절의 구조적 통합

절 주제 교훈
§ 5.4 정규 계층 모형 \(p(\tau \mid y)\) 유도·3 단계 샘플링
§ 5.5 8 학교 Separate/Pooling 의 두 극단 피하기
§ 5.6 메타분석 Complete pooling 의 과신

공통 패턴 — \(\tau\) 의 사후가 “얼마나 공유” 를 말한다. \(\tau \to 0\) (pool), \(\tau \to \infty\) (separate) 의 두 극단 사이를 데이터가 결정.

6 빈도주의와의 대응

질문 빈도주의 베이즈 (§ 5.4~5.6)
여러 그룹 평균 ANOVA F-test 정규 계층, \(p(\tau \mid y)\)
Mixed-effects 모델 REML, ML Full Bayes
메타분석 Fixed/Random effects 정규 계층
다중 비교 Bonferroni, Tukey 자동 shrinkage
새 연구 예측 해결 어려움 \(\tilde\theta\) 사후 예측

REML (Restricted Maximum Likelihood) 이 빈도주의의 \(\tau^2\) 추정 표준. 베이즈 사후 평균 ≈ REML 추정 (대표본에서), 그러나 사후 전체 분포 가 주는 정보가 훨씬 풍부.

7 코드 예제 — 8 학교 완전 분석

7.1 Step 1: 순수 Python — \(p(\tau \mid y)\) 격자와 3 단계 샘플링

import math
import random

random.seed(42)

# 8 학교 데이터
y = [28, 8, -3, 7, -1, 1, 18, 12]
sigma = [15, 10, 16, 11, 9, 11, 10, 18]
J = 8

def log_post_tau(tau, y, sigma):
    # 식 (5.21): 균등 p(τ)
    # p(τ|y) ∝ V_μ^{1/2} × ∏ (σ_j² + τ²)^{-1/2} × exp(-(y_j-μ̂)²/(2(σ_j²+τ²)))
    V_inv = sum(1 / (s**2 + tau**2) for s in sigma)
    mu_hat = sum(y[j] / (sigma[j]**2 + tau**2) for j in range(J)) / V_inv

    log_p = 0.5 * math.log(1/V_inv)  # V_μ^{1/2}
    for j in range(J):
        log_p += -0.5 * math.log(sigma[j]**2 + tau**2)
        log_p += -(y[j] - mu_hat)**2 / (2 * (sigma[j]**2 + tau**2))
    return log_p

# τ 격자 [0.01, 40]
tau_grid = [0.01 + i * 0.5 for i in range(81)]
log_post_vals = [log_post_tau(tau, y, sigma) for tau in tau_grid]

# 정규화
max_lp = max(log_post_vals)
probs = [math.exp(lp - max_lp) for lp in log_post_vals]
total = sum(probs)
probs = [p / total for p in probs]

# τ 사후 분위수
cum = 0
q_25 = q_50 = q_975 = None
for i, p in enumerate(probs):
    cum += p
    if q_25 is None and cum >= 0.025: q_25 = tau_grid[i]
    if q_50 is None and cum >= 0.5: q_50 = tau_grid[i]
    if q_975 is None and cum >= 0.975: q_975 = tau_grid[i]; break

print(f"τ 사후 분위수: 2.5% = {q_25:.2f}, 중앙 = {q_50:.2f}, 97.5% = {q_975:.2f}")

# 3 단계 샘플링
S = 5000
theta_samples = [[] for _ in range(J)]

for _ in range(S):
    # 1. τ ~ p(τ|y) 격자에서
    u = random.random()
    cum = 0
    tau = tau_grid[-1]
    for i, p in enumerate(probs):
        cum += p
        if u <= cum:
            tau = tau_grid[i]
            break

    # 2. μ | τ, y ~ N(μ̂, V_μ)
    V_inv = sum(1 / (sigma[j]**2 + tau**2) for j in range(J))
    mu_hat = sum(y[j] / (sigma[j]**2 + tau**2) for j in range(J)) / V_inv
    mu = random.gauss(mu_hat, math.sqrt(1/V_inv))

    # 3. θ_j | μ, τ, y ~ N(θ̂_j, V_j) 각 j 독립
    for j in range(J):
        V_j = 1 / (1/sigma[j]**2 + 1/tau**2)
        theta_hat_j = V_j * (y[j]/sigma[j]**2 + mu/tau**2)
        theta_samples[j].append(random.gauss(theta_hat_j, math.sqrt(V_j)))

# 학교별 사후 평균·95% 구간
print(f"\n{'학교':<6} {'y_j':<6} {'σ_j':<6} {'사후 평균':<12} {'95% 구간':<20}")
labels = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H']
for j in range(J):
    s = sorted(theta_samples[j])
    mean = sum(s)/S
    lo = s[int(0.025*S)]
    hi = s[int(0.975*S)]
    print(f"{labels[j]:<6} {y[j]:<6} {sigma[j]:<6} {mean:<12.1f} [{lo:.1f}, {hi:.1f}]")

예상 출력 — 학교 A (원 28) 사후 평균 약 11, 95% 구간 약 \([-5, 30]\). 여전히 불확실하지만 원 관측 28 보다 중심에 가깝다.

7.2 Step 2: NumPy - 8 학교의 다중 비교 자동 해결

import numpy as np

# Step 1 의 theta_samples 를 numpy 로
ts = np.array(theta_samples)  # shape (J, S)

# 28 개 짝 비교 — 모든 쌍의 Pr(θ_i > θ_j | y)
print(f"\n쌍별 Pr(θ_i > θ_j | y):")
print(f"{'':5s}", end='')
for k in range(J):
    print(f"{labels[k]:>6s}", end='')
print()
for i in range(J):
    print(f"{labels[i]:>3s}: ", end='')
    for j in range(J):
        if i == j:
            print(f"{'—':>6s}", end='')
        else:
            p = (ts[i] > ts[j]).mean()
            print(f"{p:>6.2f}", end='')
    print()

# 차이의 95% 구간이 0 을 포함하는 비율
n_pairs = 0
n_contains_zero = 0
for i in range(J):
    for j in range(i+1, J):
        n_pairs += 1
        diffs = ts[i] - ts[j]
        lo, hi = np.percentile(diffs, [2.5, 97.5])
        if lo <= 0 <= hi:
            n_contains_zero += 1
print(f"\n28 개 짝 차이 중 95% 구간이 0 포함: {n_contains_zero}/{n_pairs}")

예상 출력 — 28/28 짝 모두 0 포함. 다중 비교가 자동으로 해결되어 “어느 학교가 더 낫다” 는 결론 없음.

8 관련 주제

Ch.5 의 다른 심화

Ch.1~4 심화 (선행)

Part I~V 전체

빈도주의 대응

9 참고자료

  • Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). CRC Press. Ch.5 (§ 5.4~5.6).
  • Rubin, D. B. (1981). Estimation in parallel randomized experiments. Journal of Educational Statistics, 6(4), 377–401. [8 학교 원전]
  • Yusuf, S., Peto, R., Lewis, J., Collins, R., & Sleight, P. (1985). Beta blockade during and after myocardial infarction. Progress in Cardiovascular Diseases, 27(5), 335–371. [베타차단제 메타분석]
  • DerSimonian, R., & Laird, N. (1986). Meta-analysis in clinical trials. Controlled Clinical Trials, 7(3), 177–188.
  • Gelman, A., & Hill, J. (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.

Subscribe

Enjoy this blog? Get notified of new posts by email: