Ch.5.1~5.3 — 계층 모형의 언어와 쥐 종양 Beta-Binomial 완전 분석

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

Ch.5 개요 가 계층 모형 전체 (§ 5.1~5.7) 를 훑었다면, 이 포스트는 § 5.1~5.3 을 교재 원문 수준으로 풀어 쓴 심화편이다. 세 절이 묶이는 이유 — § 5.1 이 동기, § 5.2 가 수학적 기반, § 5.3 이 구체 계산 의 삼각 구조로 Beta-Binomial 계층 모형 하나를 완전히 해부한다.

§ 5.1~5.3 의 한 줄 요약

“이력 데이터가 있을 때 ‘공통 모집단’ 이라는 교환가능성 가정이 자연스럽고, 이 가정이 계층 모형의 수학적 기초가 된다. 쥐 종양 예제에서 \((\alpha, \beta)\) 의 격자 계산과 각 \(\theta_j\) 의 조건부 Beta 추출로 1000 개 사후 표본을 얻는 전 과정이 Part III MCMC 의 단순 버전이다.”

Ch.5 의 나머지 (정규 교환가능, 8 학교, 메타분석, 약정보적 사전) 가 모두 이 세 절의 패턴을 반복 (Gelman et al., 2013, Ch.5.1~5.3).

---

2 § 5.1 매개변수화된 사전 분포의 구성

2.1 동기 — 한 실험의 한계

한 실험의 소표본 이 주는 추정은 불안정. 현재 실험 외에 유사한 이전 실험들 (historical data) 이 있다면 — 현재와 과거를 공통 모집단에서 추출된 무작위 표본 으로 수학적 모델링 가능.

2.2 쥐 종양 예제 (Tarone, 1982)

배경. F344 암컷 실험용 쥐 — 약물 평가를 위한 표준 모델 동물. 통제 그룹 (zero dose) 에서의 endometrial stromal polyp (자궁내막 기질 용종, 일종의 종양) 발생 확률 \(\theta\) 추정.

현재 실험 — 14 마리 중 4 마리 종양 발생.

단순 추정. \(\hat\theta = 4/14 \approx 0.286\), 95% Wilson 구간 \(\approx [0.12, 0.54]\). 매우 불안정.

2.3 이력 데이터 (표 5.1)

70 개 이전 실험 의 \(y_j/n_j\) 값. 교재 인용 — 범위가 “0/20, 0/19, 0/18, …” 에서 “16/52, 15/47, 15/46” 까지. 샘플 비율 분포가 광범위.

표본 통계량.

• 70 값 \(y_j/n_j\) 의 표본 평균 0.136
• 표본 표준편차 0.103
• 대부분 0~0.3 범위, 몇 개 극단값

2.4 단순 접근 1 — 고정 사전

만약 모집단 \(\theta\) 분포가 알려진 Beta(\(\alpha, \beta\)) 라면.

\[ \theta \mid y \sim \text{Beta}(\alpha + 4, \beta + 10) \]

\((\alpha, \beta)\) 는 이력 데이터의 평균과 분산 으로 매칭.

Beta 평균·분산 공식.

\[ E = \frac{\alpha}{\alpha+\beta} = 0.136, \quad \text{Var} = \frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)} = 0.103^2 \]

이 두 식을 풀면 \((\alpha, \beta) \approx (1.4, 8.6)\). 사후.

\[ \theta \mid y \sim \text{Beta}(1.4 + 4, 8.6 + 10) = \text{Beta}(5.4, 18.6) \]

사후 평균 \(5.4/24 \approx 0.225\).

2.5 고정 사전의 문제

두 가지 결함.

1. \((\alpha, \beta)\) 의 불확실성 무시 — 70 이력 실험도 노이즈가 있다. 추정된 \(\alpha, \beta\) 자체에 불확실성이 있어야 한다
2. 현재 실험의 기여 0 — 현재 데이터가 이력 모집단 분포 \((\alpha, \beta)\) 의 추정에도 기여해야 한다 (모두 같은 모집단에서 추출)

2.6 Full Bayes 접근 — 계층 모형

\((\alpha, \beta)\) 를 미지 모수로 두고 hyperprior 부여. 70 + 1 = 71 실험의 데이터 모두를 사용하여 \((\alpha, \beta, \theta_1, \ldots, \theta_{71})\) 의 공동 사후 계산.

이것이 계층 모형.

2.7 3 단계 계층 다이어그램

\[ \begin{array}{ccc} (\alpha, \beta) & \sim & p(\alpha, \beta) \quad \text{(hyperprior)} \\ \downarrow & & \\ \theta_1, \ldots, \theta_{71} & \mid \alpha, \beta \sim & \text{Beta}(\alpha, \beta) \quad \text{(모집단 분포)} \\ \downarrow & & \\ y_j & \mid \theta_j \sim & \text{Bin}(n_j, \theta_j) \quad \text{(관측 모델)} \end{array} \]

계층 = 정보 흐름의 위계. \((\alpha, \beta)\) 는 모집단 전체의 특성, \(\theta_j\) 는 각 실험의 특성, \(y_j\) 는 관측.

직관 — 계층 모형이 “두 세계” 를 연결한다

“이 실험만” (소표본, 불안정) vs “모든 실험이 같다” (과도한 공유) 의 두 극단 사이에서, 계층 모형이 데이터가 말하는 만큼 공유 하게 한다.

• \(\text{Var}(\theta_j) = 0\) 이면 모든 실험이 같음 (pooling)
• \(\text{Var}(\theta_j) = \infty\) 이면 완전 독립 (no pooling)
• 데이터가 중간 \(\text{Var}(\theta_j) = \tau^2 > 0\) 을 선호하면 partial pooling

\(\tau\) 가 “공유 정도” 를 자동 결정. 이것이 Ch.5 의 핵심 메시지.

---

3 § 5.2 교환가능성과 계층 모형

3.1 교환가능성의 정의 재확인

\((\theta_1, \ldots, \theta_J)\) 의 공동 분포 \(p(\theta_1, \ldots, \theta_J)\) 가 인덱스 순열에 불변 하면 교환가능.

\[ p(\theta_1, \ldots, \theta_J) = p(\theta_{\sigma(1)}, \ldots, \theta_{\sigma(J)}) \quad \forall \sigma \]

3.2 교환가능성 = 정보 부족의 수학

교재의 선언.

“데이터 \(y\) 외에 \(\theta_j\) 들을 구별할 정보가 없고 어떤 순서나 그룹핑도 불가능하면, 모수들의 사전 분포에 대칭을 가정해야 한다. 이 대칭이 확률적으로 교환가능성으로 표현된다.” (교재)

3.3 쥐 종양에서의 교환가능성

71 실험 중 — \(n_j\) (표본 크기) 외에 구별 정보 없음. 그리고 \(n_j\) 가 \(\theta_j\) (종양률) 와 관련 있을 선험적 이유도 없음 → 교환가능 모델 정당.

3.4 iid 혼합 형태 (5.2)

교환가능성의 가장 간단한 형태 — 각 \(\theta_j\) 가 미지 \(\phi\) 조건부로 독립 동일 분포.

\[ p(\theta \mid \phi) = \prod_{j=1}^J p(\theta_j \mid \phi) \tag{5.1} \]

\(\phi\) 가 미지이므로 불확실성을 반영해 주변화.

\[ p(\theta) = \int \left[\prod_{j=1}^J p(\theta_j \mid \phi)\right] p(\phi) \, d\phi \tag{5.2} \]

“독립 동일 분포들의 혼합” — 대부분의 실무에서 교환가능성을 포착하는 충분한 표현.

3.5 de Finetti 정리 — 이론적 기반

\(J \to \infty\) 에서 임의의 적절한 교환가능 분포 가 (5.2) 형태의 iid 혼합 으로 표현 가능 (Ch.1.2 에서 소개).

주의. \(J\) 가 유한이면 정리가 엄밀하게는 성립 안 함 (Ex 5.1, 5.2, 5.4). 실무적으로는 근사적으로 성립.

반례. 주사위 6 면의 확률 \(\theta_1, \ldots, \theta_6\). 교환가능하지만 \(\sum \theta_j = 1\) 제약 때문에 독립 혼합 모델 불가. 하지만 여전히 교환가능으로 모델링 가능 (Dirichlet 같은 제약 있는 분포).

3.6 이혼율 사고실험 — 교환가능성의 스펙트럼

교재의 유명한 사고 실험. 미국 8 주의 1981 년 이혼률 \(y_1, \ldots, y_8\). \(y_8\) 에 대해 무엇을 말할 수 있나?

시나리오 A — 아무 정보 없음.

독자는 완전 교환가능 사전 사용. 8 개 \(y_j\) 에 공통 분포 부여. 예 — logit-normal, 특정 범위에 집중된 Beta 등.

\(y_1, \ldots, y_7 = (5.8, 6.6, 7.8, 5.6, 7.0, 7.1, 5.4)\) 관측 → \(y_8\) 의 사후 예측은 약 6.5 중심, 5.0~8.0 범위. 관측 인덱스 재라벨링해도 결과 같음 (교환가능성).

시나리오 B — “8 개 모두 산악 주”.

주 이름 목록 제공. Arizona, Colorado, Idaho, Montana, Nevada, New Mexico, Utah, Wyoming. 어떤 \(y_j\) 가 어떤 주인지는 알려지지 않음 → 여전히 교환가능.

그러나 사전 정보가 바뀜.

• Utah — 몰몬교 인구 많음, 낮은 이혼률 예상
• Nevada — 관대한 이혼법, 높은 이혼률 예상

사전 분포에 두꺼운 꼬리 (이상치 허용) 를 배치. 관측된 7 값이 매우 가까우니 빠진 \(y_8\) 이 Nevada 또는 Utah 일 가능성 → 이봉 또는 삼봉 사후.

시나리오 C — \(y_8\) = Nevada 라고 알려줌.

이제 \(y_8\) 이 다른 7 개와 구별됨 → 교환가능성 무너짐. \(y_8\) 이 관측된 최댓값 (7.8) 보다 클 확률이 매우 높다.

실제 값. Nevada 1981 년 이혼률 = 13.9 per 1,000 population. 관측 최댓값의 거의 2 배 — 사후 예측이 극단값을 포착한 판단이 옳았다.

직관 — 교환가능성이 정보 수준에 따라 어떻게 변하는가

정보 수준	교환가능성	처리
정보 0	완전	단순 교환가능 사전
공통 범주	부분	사전 모양 조정 (이상치 허용)
일부 인덱스 식별	일부 비교환	그 인덱스를 조건부로 분리
완전 식별	무너짐	공변량 \(x\) 로 회귀

핵심 통찰 — “교환가능성” 은 이분법이 아니라 정보 수준에 따른 연속 스펙트럼. 관련 공변량을 모델에 투입할수록 잔여가 점점 교환가능에 가까워진다. 이것이 회귀 모델의 계층적 정당화.

3.7 공변량 조건부 교환가능성

관측이 \((x_j, y_j)\) 쌍이고 \(x_j\) 가 \(\theta_j\) 와 관련된 공변량일 때. \(y_j\) 자체는 비교환가능해도 \((x_j, y_j)\) 는 교환가능.

\[ p(\theta_1, \ldots, \theta_J \mid x_1, \ldots, x_J) = \int \left[\prod_{j=1}^J p(\theta_j \mid \phi, x_j)\right] p(\phi \mid x) \, d\phi \]

조건부 독립 구조. 이것이 회귀 모델의 계층적 정당화 — Part IV Ch.14 의 수학적 기반.

3.8 쥐 종양 예제에서의 공변량

\(n_j\) 가 유일한 잠재적 공변량. 교재의 추측.

“아마도 일부 연구 \(j\) 는 투자자가 정확히 희귀한 사건을 의심해서 더 큰 \(n_j\) 를 설정했을지도 모른다 — 즉 작은 \(\theta_j\) 와 작은 \(y_j/n_j\) 기댓값. 그러나 실제로 \(y_j/n_j\) vs \(n_j\) 플롯은 명백한 관계를 보이지 않는다.” (교재)

따라서 \(n_j\) 를 공변량으로 쓰지 않고 단순 교환가능 모델.

3.9 교환가능성에 대한 반대

흔한 반대. “71 실험이 실제로 다른 시간, 다른 쥐, 다른 실험실 에서 수행됐으므로 교환가능 모델은 부적절”.

교재의 답변.

“그러나 그런 정보가 교환가능성을 무효화하지는 않는다. 실험들이 다르다는 것은 \(\theta_j\) 가 다르다는 뜻이지만, 공통 분포에서 추출된 것처럼 간주 하는 것은 완벽히 수용 가능. 실제로 구별 정보가 없으면 \(\theta_j\) 를 교환가능하게 모델링하는 것 외에 논리적 선택이 없다.” (교재)

핵심 논리 — “실험들이 다르다” ≠ “모델이 틀렸다”. 교환가능성은 “실험들이 모두 같다” 는 주장이 아니라 “사전 모집단에서 추출된 대칭 구조” 이다.

3.10 Full Bayesian Treatment

계층 모형의 “계층” 부분 — \(\phi\) 가 미지이고 고유 사전 \(p(\phi)\) 를 갖는다.

공동 사전.

\[ p(\phi, \theta) = p(\phi) \, p(\theta \mid \phi) \]

공동 사후.

\[ p(\phi, \theta \mid y) \propto p(\phi, \theta) \, p(y \mid \phi, \theta) = p(\phi) \, p(\theta \mid \phi) \, p(y \mid \theta) \tag{5.3} \]

가능도가 \(\theta\) 만 통해 \(\phi\) 에 의존 — 이것이 “hyperparameter” 라는 이름의 의미. \(\phi\) 가 \(\theta\) 를 통해 간접적으로 데이터에 영향.

3.11 Hyperprior 선택 원리

교재의 지침.

“\(\phi\) 에 대해 알려진 것이 적으면 diffuse 사전 부여 가능, 단 improper 사전 사용 시 사후 proper 인지 반드시 확인. 결론이 이 단순화 가정에 얼마나 민감한지 평가해야 한다.” (교재)

실무 권장 — 단순하고 비교적 비정보적 사전에서 시작, 필요하면 정보 추가.

3.12 사후 예측 분포의 두 종류

계층 모형에서 두 가지 예측 이 가능.

1. 기존 \(\theta_j\) 에 대한 새 관측 — 예: 실험 \(j\) 에 쥐 추가. 기존 \(\theta_j\) 의 사후에서 표본 → \(\tilde{y} \sim p(\tilde{y} \mid \theta_j)\)
2. 새 실험의 새 \(\tilde\theta\) 에 대한 관측 — 먼저 \(\tilde\theta \sim p(\theta \mid \phi)\) (\(\phi\) 사후에서 \(\phi\) 추출) 후 \(\tilde{y} \sim p(\tilde{y} \mid \tilde\theta)\)

두 예측 모두 모델 점검 (Ch.6) 에 사용.

직관 — 두 예측의 차이

예측 1 — “이 실험에서 쥐 한 마리 더 관측하면?”. 해당 실험의 \(\theta_j\) 불확실성만 반영.

예측 2 — “새 실험 수행하면?”. \(\theta\) 변동 + \(\tilde{y}\) 관측 변동 모두 반영. 예측 분포가 더 넓음.

실무 함의 — 신약 개발에서 “현재 실험의 결과” vs “새 실험 예측” 은 서로 다른 질문. 두 분포 모두 필요.

---

4 § 5.3 Beta-Binomial 계층의 Bayesian 분석

4.1 전략 — 공동 → 조건부 × 주변

Ch.3.1 의 기법을 계층 모형에 적용. 3 단계 해석적 작업.

1. 공동 사후 \(p(\theta, \phi \mid y)\) 를 비정규화 형태 로 쓴다
2. 조건부 사후 \(p(\theta \mid \phi, y)\) 를 해석적으로 유도 (켤레 덕분에 단순)
3. 주변 사후 \(p(\phi \mid y)\) 계산 — 적분 또는 조건부 확률 공식

4.2 조건부 확률 공식

\[ p(\phi \mid y) = \frac{p(\theta, \phi \mid y)}{p(\theta \mid \phi, y)} \tag{5.5} \]

임의의 \(\theta\) 값에서 평가 — 분자는 공동 사후, 분모는 조건부 사후. 두 분포 모두 닫힌 형태이면 주변이 해석적.

4.3 경고 — 정규화 상수

교재의 주의.

“분모 \(p(\theta \mid \phi, y)\) 를 \(\theta, \phi\) 모두의 함수로 보면, 정규화 인자가 \(\phi\) 에 의존 할 수 있다. 베이즈 정리에서 비례 상수가 실제로 \(\phi\) 에 대해 상수 인지 확인해야 한다.” (교재)

4.4 쥐 종양 모델 설정

가능도.

\[ y_j \sim \text{Bin}(n_j, \theta_j), \quad j = 1, \ldots, J = 71 \]

모집단.

\[ \theta_j \mid \alpha, \beta \sim \text{Beta}(\alpha, \beta) \]

Hyperprior — 나중에 지정.

4.5 공동 사후 (5.6)

\[ p(\theta, \alpha, \beta \mid y) \propto p(\alpha, \beta) \prod_{j=1}^J \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta_j^{\alpha - 1}(1-\theta_j)^{\beta - 1} \prod_{j=1}^J \theta_j^{y_j}(1-\theta_j)^{n_j - y_j} \]

Beta 밀도 × 이항 가능도의 곱.

4.6 조건부 사후 (5.7)

\((\alpha, \beta)\) 고정 하에서 \(\theta_j\) 들이 독립 이고 각각 Beta 형태.

\[ \theta_j \mid \alpha, \beta, y \sim \text{Beta}(\alpha + y_j, \beta + n_j - y_j) \]

\[ p(\theta \mid \alpha, \beta, y) = \prod_{j=1}^J \frac{\Gamma(\alpha + \beta + n_j)}{\Gamma(\alpha + y_j)\Gamma(\beta + n_j - y_j)} \theta_j^{\alpha + y_j - 1}(1 - \theta_j)^{\beta + n_j - y_j - 1} \tag{5.7} \]

4.7 주변 사후 (5.8)

식 (5.5) 에 대입, \(\theta_j\) 약분 후.

\[ p(\alpha, \beta \mid y) \propto p(\alpha, \beta) \prod_{j=1}^J \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \frac{\Gamma(\alpha + y_j)\Gamma(\beta + n_j - y_j)}{\Gamma(\alpha + \beta + n_j)} \tag{5.8} \]

\((\alpha, \beta)\) 의 2 차원 함수 — 격자 계산 가능.

이 곱은 해석적 단순화 불가 이지만 감마 함수 계산으로 수치 평가 쉽다.

4.8 매개변수화 선택

\(\alpha, \beta > 0\) 제약 공간 → 제약 없는 실수 공간으로 변환.

교재 선택.

\[ \text{logit}\left(\frac{\alpha}{\alpha+\beta}\right) = \log\left(\frac{\alpha}{\beta}\right), \quad \log(\alpha + \beta) \]

해석.

• \(\log(\alpha/\beta)\) — 모집단 평균 \(\alpha/(\alpha+\beta)\) 의 logit. 평균이 0.5 이면 0, 작으면 음수, 크면 양수
• \(\log(\alpha+\beta)\) — “유효 표본 크기” \(\alpha + \beta\) 의 로그. Beta 의 precision 또는 “사전 관측 수”

이 두 값이 대략 독립 이라는 직관 (평균과 precision 은 서로 독립적 정보).

직관 — 왜 이 매개변수화가 “자연스러운가”

Beta 분포의 해석.

• 평균 \(\alpha/(\alpha+\beta)\) — 분포 중심 (확률의 모드)
• \(\alpha + \beta\) — 정보량 (큰 값 = 뾰족한 분포, 작은 값 = 평평한 분포)

이 두 측도가 수학적으로 \(\alpha, \beta\) 의 함수이지만 해석적으로 독립. 평균을 logit 으로 (확률에서 실수 전체), precision 을 log 로 (양수에서 실수 전체) 변환.

통계학의 많은 모델이 같은 패턴 — “중심” 과 “scale/precision” 의 분리. 정규의 \((\mu, \log\sigma)\), 감마의 \((\log E, \log \kappa)\) 등.

4.9 Hyperprior 선택

첫 시도 — \((\log(\alpha/\beta), \log(\alpha+\beta))\) 균등.

문제 — \((\alpha + \beta) \to \infty\) 에서 사후 적분 발산. Improper 사후.

해결 — (5.9) 권장.

\[ p(\alpha, \beta) \propto (\alpha + \beta)^{-5/2} \tag{5.9} \]

이것은 \(\alpha/(\alpha+\beta)\) 에 균등, \((\alpha+\beta)^{-1/2}\) 에 균등인 사전에 해당. 야코비안 변환 후.

\[ p(\log(\alpha/\beta), \log(\alpha+\beta)) \propto \alpha\beta(\alpha+\beta)^{-5/2} \tag{5.10} \]

결과 사후 proper 임을 확인 가능.

4.10 Hyperprior 선택의 교훈

“Proper hyperprior 를 쓰면 수학적 노력을 피할 수 있다. 또 다른 접근은 잠정적으로 flat 사전 (예: \(p(\frac{\alpha}{\alpha+\beta}, \alpha+\beta) \propto 1\) 또는 \(p(\alpha, \beta) \propto 1\)) 사용 후 사후 등고선과 시뮬을 계산. 결과가 무한대로 drift 하면 사후가 integrable 하지 않다는 신호 — 사전 수정 필요.” (교재)

추가 경고 — improper 사전의 “flat on \([-10^{10}, 10^{10}]\)” 같은 대체는 해결이 아니다. 거의 모든 사후 질량이 \(\infty\) 쪽에 몰려 \(\text{Var}(\theta_j) = 0\) (complete pooling) 의 극단으로 귀결.

4.11 격자 계산

범위 설정. 교재의 첫 시도 — \((\log(\alpha/\beta), \log(\alpha+\beta)) \in [-2.5, -1] \times [1.5, 3]\). 점 추정 \((-1.8, 2.3)\) 중심.

그림 5.2 관찰. 모드는 점 추정 근처. 그러나 사후의 중요한 부분이 범위 밖 — 범위 확대 필요.

수정된 범위. \([-2.3, -1.3] \times [1, 5]\). 그림 5.3a 가 전체 사후 덮기 — 모드 약 \((-1.75, 2.8)\), \((α, β) \approx (2.4, 14.0)\) 대응.

4.12 표본 추출 알고리즘

1. 격자에서 \((\log(\alpha/\beta), \log(\alpha+\beta))\) 추출 — 이산 근사, 역 CDF 법 (Ch.3.7 bioassay 와 동일)
2. \((\alpha, \beta)\) 로 역변환
3. 각 \(j\) 에 대해 \(\theta_j \mid \alpha, \beta, y \sim \text{Beta}(\alpha + y_j, \beta + n_j - y_j)\) 독립 추출

\(L = 1000\) 회 반복.

4.13 수치 결과

격자 기반 사후 모멘트.

\[ E(\alpha \mid y) = 2.4, \quad E(\beta \mid y) = 14.3 \]

모집단 평균 \(\alpha/(\alpha+\beta) \approx 0.144\), 모집단 표준편차 \(\sqrt{\alpha\beta/[(\alpha+\beta)^2(\alpha+\beta+1)]} \approx 0.080\).

이력 데이터 요약 (0.136, 0.103) 과 유사 — 베이즈 추정이 데이터 기반 요약을 잘 재현.

4.14 모든 71 실험의 사후

각 \(\theta_j\) 의 사후 중앙값과 95% 구간 (그림 5.4 산점도 vs \(y_j/n_j\)).

관찰.

1. 45° 선 (unpooled) 에서 모집단 평균 (0.14) 쪽으로 shrinkage
2. \(n_j\) 가 작을수록 더 많이 shrunk — 공유 정보 의존도 증가
3. \(n_j\) 가 클수록 거의 변화 없음 — 자체 데이터에 가깝게 유지
4. 95% 구간이 unpooled 구간보다 좁음 — 이력 데이터의 정보 기여

4.15 현재 실험 (71 번째) — 구체 수치

\(y_{71} = 4, n_{71} = 14\). 사후 \(\theta_{71}\) 의 평균.

\[ E(\theta_{71} \mid y) \approx \frac{\alpha + y_{71}}{\alpha + \beta + n_{71}} = \frac{2.4 + 4}{2.4 + 14.3 + 14} \approx \frac{6.4}{30.7} \approx 0.208 \]

비교.

• Unpooled: \(4/14 = 0.286\)
• 이력 평균: 0.136
• 계층 사후: 약 0.21 (두 값의 중간, 현재 쪽으로 치우침)

\(n_{71} = 14\) 가 상대적으로 작아서 이력 쪽으로 약 40% shrinkage.

4.16 Empirical Bayes 와의 미묘한 차이

교재의 언급.

“결과는 표면적으로 hyperparameter 점 추정에 기반한 결과와 유사 하다 — 실험이 많기 때문에. 그러나 중요한 차이 가 남는다: Full Bayesian 분석에서 사후 변동이 더 높다 — hyperparameter 의 사후 불확실성을 반영하기 때문.” (교재)

Empirical Bayes (EB) — \((\hat\alpha, \hat\beta)\) 점 추정 후 고정. Full Bayes — \((\alpha, \beta)\) 의 사후 분포 전체 반영. 두 \(\theta_j\) 의 점 추정은 비슷하지만 FB 의 구간이 EB 보다 넓다.

---

5 세 절의 구조적 통합

절	주제	핵심 결과
§ 5.1	동기	이력 데이터 → 계층 구조
§ 5.2	수학적 기반	교환가능성 → iid 혼합 (5.2)
§ 5.3	계산	공동 → 조건부 × 주변, 격자

Beta-Binomial 하나 로 계층 모형의 전 과정을 체험. § 5.4~5.7 의 정규·8 학교·메타분석이 같은 패턴 반복.

---

6 빈도주의 대응 — 계층 모형의 첫 번째 사례

질문	빈도주의	베이즈 (§ 5.1~5.3)
이력 데이터 활용	메타분석 고정 효과	계층 모형 full Bayes
소표본 실험	큰 신뢰 구간	이력 평균으로 shrinkage
모집단 분포 \((\alpha, \beta)\)	Method of moments	주변 사후 \(p(\alpha, \beta \mid y)\)
개별 \(\theta_j\)	단일 실험 분석	부분 풀링

James-Stein 추정량 이 이 결과의 빈도주의 대응 — 여러 그룹의 추정을 공통 평균으로 shrink 하면 MSE 감소.

---

7 코드 예제 — 쥐 종양 계층 모형 완전 구현

7.1 Step 1: 순수 Python — 격자 계산과 71 실험 shrinkage

import math
import random

random.seed(42)

# 쥐 종양 데이터 — 70 이력 + 1 현재 (교재 표 5.1 일부 재구성)
# 실제로는 71개 (y_j, n_j) 쌍. 예시로 대표적 몇개 + 현재
ys = [0, 0, 0, 0, 1, 1, 2, 2, 5, 5, 9, 10, 4, 4, 6, 6, 16, 15, 4]
ns = [20, 20, 19, 18, 20, 19, 20, 20, 49, 46, 48, 50, 20, 20, 22, 23, 52, 47, 14]  # 마지막이 현재
J = len(ys)

def log_gamma(x):
    # Ramanujan 근사
    return math.lgamma(x)

def log_marginal_alpha_beta(alpha, beta, ys, ns):
    # log p(α, β | y) — 하이퍼사전 (α+β)^{-5/2} 포함
    if alpha <= 0 or beta <= 0:
        return float("-inf")
    log_hyperprior = -2.5 * math.log(alpha + beta)
    log_lik = 0.0
    for y, n in zip(ys, ns):
        log_lik += log_gamma(alpha + beta) - log_gamma(alpha) - log_gamma(beta)
        log_lik += log_gamma(alpha + y) + log_gamma(beta + n - y) - log_gamma(alpha + beta + n)
    return log_hyperprior + log_lik

# 격자 — log(α/β) ∈ [-2.5, -1], log(α+β) ∈ [1, 5]
grid_log_ratio = [-2.5 + i * 0.05 for i in range(31)]  # 31 points
grid_log_sum = [1 + i * 0.1 for i in range(41)]  # 41 points

log_post = []
for lr in grid_log_ratio:
    for ls in grid_log_sum:
        alpha_beta_sum = math.exp(ls)
        alpha = alpha_beta_sum / (1 + math.exp(-lr))
        beta = alpha_beta_sum - alpha
        # 야코비안 포함
        lp = log_marginal_alpha_beta(alpha, beta, ys, ns) + math.log(alpha) + math.log(beta) - 2.5 * math.log(alpha_beta_sum)
        # Note: 야코비안 d(α,β)/d(log(α/β), log(α+β)) ∝ αβ/(α+β)
        log_post.append((lr, ls, lp))

# 정규화
max_lp = max(lp for _, _, lp in log_post)
weights = [(lr, ls, math.exp(lp - max_lp)) for lr, ls, lp in log_post]
total = sum(w for _, _, w in weights)
probs = [(lr, ls, w / total) for lr, ls, w in weights]

# 모드 찾기
mode_idx = max(range(len(probs)), key=lambda i: probs[i][2])
print(f"사후 모드: log(α/β) = {probs[mode_idx][0]:.2f}, log(α+β) = {probs[mode_idx][1]:.2f}")

# α, β 점 추정
alpha_mode = math.exp(probs[mode_idx][1]) / (1 + math.exp(-probs[mode_idx][0]))
beta_mode = math.exp(probs[mode_idx][1]) - alpha_mode
print(f"α ≈ {alpha_mode:.2f}, β ≈ {beta_mode:.2f}")

# 현재 실험 (j = J-1) 의 사후 평균 — 사전 모드 기반
y_cur, n_cur = ys[-1], ns[-1]
theta_cur_mean = (alpha_mode + y_cur) / (alpha_mode + beta_mode + n_cur)
theta_cur_unpooled = y_cur / n_cur

print(f"\n현재 실험 (y=4, n=14):")
print(f"  Unpooled: {theta_cur_unpooled:.4f}")
print(f"  계층 사후 (α, β 모드 사용): {theta_cur_mean:.4f}")

예상 출력 — \((\alpha, \beta) \approx (2.4, 14.0)\), 현재 실험 unpooled 0.286 → 계층 약 0.21.

7.2 Step 2: 완전 사후 표본 추출 — 71 실험의 shrinkage

import random

def sample_beta(a, b):
    x = random.gammavariate(a, 1)
    y = random.gammavariate(b, 1)
    return x / (x + y)

# Step 1 에서 probs 를 이용, (log(α/β), log(α+β)) 추출 후 θ_j 추출
L = 1000
theta_samples = [[] for _ in range(J)]

probs_flat = [(lr, ls, p) for lr, ls, p in probs]
cum_probs = [0.0]
for _, _, p in probs_flat:
    cum_probs.append(cum_probs[-1] + p)

for _ in range(L):
    # (log(α/β), log(α+β)) 추출
    u = random.random()
    idx = 0
    for i in range(1, len(cum_probs)):
        if u <= cum_probs[i]:
            idx = i - 1
            break
    lr, ls, _ = probs_flat[idx]
    alpha_beta = math.exp(ls)
    alpha = alpha_beta / (1 + math.exp(-lr))
    beta = alpha_beta - alpha

    # 각 θ_j 추출
    for j in range(J):
        theta_j = sample_beta(alpha + ys[j], beta + ns[j] - ys[j])
        theta_samples[j].append(theta_j)

# 71 실험의 shrinkage 요약
print(f"\n{'실험':<6} {'y/n':<10} {'unpooled':<10} {'계층 평균':<10} {'shrinkage':<10}")
for j in [0, 5, 10, 15, J-1]:  # 몇 개만 출력
    unpool = ys[j] / ns[j]
    hier_mean = sum(theta_samples[j]) / L
    shrink = unpool - hier_mean
    print(f"{j:<6} {ys[j]}/{ns[j]:<8} {unpool:<10.4f} {hier_mean:<10.4f} {shrink:<+10.4f}")

예상 관찰 — \(n_j\) 가 작고 극단 \(y_j\) (예 \(y/n = 0/20 = 0\)) 일수록 shrinkage 가 크다. 현재 실험 \(4/14\) 는 약 0.08 감소.

---

8 관련 주제

Ch.5 의 다른 심화 (후속)

• Ch.5 개요
• § 5.4~5.5 정규 계층 · 8 학교 (작성 예정)
• § 5.6~5.7 메타분석 · 분산 사전 (작성 예정)

Ch.1~4 심화 (선행)

• Ch.1 개요 · Ch.2 개요 · Ch.3 개요 · Ch.4 개요
• Ch.2.7 암 발생률 shrinkage — Ch.5 의 축소판

Part I~V 전체

• Part I · Part II · Part III · Part IV · Part V

빈도주의 대응

• 혼합 효과 모형
• James-Stein shrinkage

---

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.1~5.3).
• Tarone, R. E. (1982). The use of historical control information in testing for a trend in proportions. Biometrics, 38(1), 215–220.
• de Finetti, B. (1974). Theory of Probability. Wiley.
• Lindley, D. V., & Smith, A. F. M. (1972). Bayes estimates for the linear model. Journal of the Royal Statistical Society. Series B, 34(1), 1–41.
• Morris, C. N. (1983). Parametric empirical Bayes inference: Theory and applications. Journal of the American Statistical Association, 78(381), 47–55.