1 이 포스트의 위치 — Ch.3 심화의 마지막 조각
§ 3.1~3.3 이 정규 공동 사후, § 3.4~3.6 이 다항·다변량 정규 켤레를 다뤘다면, 이 포스트는 “비켤레 다모수 모델의 실전 계산” 과 “Ch.3 전체 기법을 체화하는 연습” 을 담는다. Bioassay 예제는 Part III MCMC 가 없던 시대에 이미 격자 계산만으로 2 모수 비선형 문제를 해결 한 고전.
“Ch.3 의 마지막은 ‘닫힌 형태가 없어도 격자 + 시뮬레이션으로 간다’ 는 교훈이다. 2 모수까지는 격자, 그 이상은 Part III 의 MCMC — Ch.3 이 닫는 Part I 의 실용 매뉴얼.”
§ 3.7 의 bioassay 는 Gelman 이 Ch.16 GLM 의 예고편으로 배치한 사례 (Gelman et al., 2013, Ch.3.7~3.10).
2 § 3.7 Bioassay 실험의 분석
2.1 과학적 배경
약물 · 화학물질 개발에서 급성 독성 시험 (acute toxicity tests) 또는 bioassay 실험 을 동물에 수행. 절차 — 다양한 용량 수준에서 동물 배치에 투여, 이분 반응 (죽음/생존, 종양/비종양) 관측.
데이터 형식.
\[ (x_i, n_i, y_i), \quad i = 1, \ldots, k \]
- \(x_i\): \(i\) 번째 용량 (대개 로그 스케일)
- \(n_i\): 이 용량 받은 동물 수
- \(y_i\): 양성 결과 (예: 사망) 수
2.2 교재의 실제 데이터 (Racine et al., 1986)
| 용량 \(x_i\) (log g/ml) | 동물 수 \(n_i\) | 사망 수 \(y_i\) |
|---|---|---|
| \(-0.86\) | 5 | 0 |
| \(-0.30\) | 5 | 1 |
| \(-0.05\) | 5 | 3 |
| \(0.73\) | 5 | 5 |
관찰 — 용량이 높을수록 사망률 증가. \(x = -0.86\) 에서 0%, \(x = 0.73\) 에서 100%. 중간 용량에서의 전환 구간이 관심.
2.3 모델링 — 이분 반응의 로지스틱 회귀
그룹 내 이항 — 각 용량 \(x_i\) 에서 \(n_i = 5\) 마리 동물의 결과는 교환가능, 조건부 독립, 같은 사망 확률 \(\theta_i\).
\[ y_i \mid \theta_i \sim \text{Bin}(n_i, \theta_i) \]
(전염병 같은 비독립 상황은 제외.)
2.4 용량-반응 관계
단순한 선택: \(\theta_i = \alpha + \beta x_i\) (선형).
문제 — \(x_i \to \pm\infty\) 에서 \(\theta_i\) 가 \([0, 1]\) 을 벗어난다. 확률 모수의 제약 위반.
표준 해결 — 로짓 변환.
\[ \text{logit}(\theta_i) = \alpha + \beta x_i \tag{3.14} \]
\(\text{logit}(\theta) = \log(\theta/(1-\theta))\). 이것이 로지스틱 회귀 (logistic regression) — Part IV Ch.16 의 핵심 모델.
\(\theta \in (0, 1)\) → \(\text{logit}(\theta) \in (-\infty, \infty)\). 선형 관계 \(\alpha + \beta x\) 가 실수 전체에서 정의되므로 자연스럽게 연결.
\(\theta_i = \text{logit}^{-1}(\alpha + \beta x_i) = \frac{e^{\alpha + \beta x_i}}{1 + e^{\alpha + \beta x_i}}\)
이것이 S 자 (sigmoid) 곡선 — 용량이 낮으면 \(\theta \approx 0\), 높으면 \(\theta \approx 1\), 중간에서 급격한 전환. 중간 전환점이 \(\alpha + \beta x = 0\) 즉 \(x = -\alpha/\beta\) — 이것이 LD50 (사망률 50% 용량).
2.5 가능도
\[ p(y_i \mid \alpha, \beta, n_i, x_i) \propto [\text{logit}^{-1}(\alpha + \beta x_i)]^{y_i} [1 - \text{logit}^{-1}(\alpha + \beta x_i)]^{n_i - y_i} \]
공동 가능도 (네 관측 곱).
\[ p(y \mid \alpha, \beta) = \prod_{i=1}^4 p(y_i \mid \alpha, \beta, n_i, x_i) \]
2.6 비정보 균등 사전
\[ p(\alpha, \beta) \propto 1 \]
\((\alpha, \beta) \in \mathbb{R}^2\) 에 균등. Improper 이지만 사후는 proper. 정보적 사전이 필요하면 다른 bioassay 실험 데이터로 구성 가능 (교재 단서).
2.7 공동 사후
\[ p(\alpha, \beta \mid y, n, x) \propto p(\alpha, \beta) \prod_i [\text{logit}^{-1}(\alpha + \beta x_i)]^{y_i} [1 - \text{logit}^{-1}(\alpha + \beta x_i)]^{n_i - y_i} \tag{3.15} \]
비켤레 — 닫힌 형태 없음. 격자 계산 으로 해결.
2.8 Rough 추정 — 격자 범위 결정
격자 계산 전에 대략적 \((\alpha, \beta)\) 위치 파악. 빈도주의 로지스틱 회귀 (MLE).
\[ (\hat{\alpha}, \hat{\beta}) = (0.8, 7.7) \]
표준오차 — \(\text{SE}(\hat{\alpha}) = 1.0, \text{SE}(\hat{\beta}) = 4.9\).
2.9 격자 설정과 계산
격자 범위 — \((\alpha, \beta) \in [-5, 10] \times [-10, 40]\). MLE 중심의 넉넉한 영역.
계산 단계.
- 격자 각 점 \((\alpha, \beta)\) 에서 로그 비정규화 사후 계산
- 수치 오버플로 방지 — 최대값 차감 후 지수화
- 격자에서 합으로 정규화 → 이산 사후 확률
교재 권장.
“로그 비정규화 사후분포 를 계산하고 지수화 전에 최대값을 차감 하는 것이 좋은 관행이다. 이것은 최댓값 1 의 비정규화 이산 근사를 만들고, 그 뒤 (격자 총 확률을 1 로 설정하여) 정규화할 수 있다.” (교재)
오버플로 주의 — 가능도의 곱이 매우 작거나 커질 수 있다. 로그 스케일에서 합 으로 안정화.
2.10 격자에서 표본 추출
원래 연속 사후를 격자 이산 분포 로 근사한 뒤 표본 추출.
알고리즘.
- \(\alpha\) 의 주변 사후 — 격자에서 \(\beta\) 축으로 합 (수치 주변화)
- For \(s = 1, \ldots, 1000\):
- \(\alpha^{(s)}\) 를 이산 \(p(\alpha \mid y)\) 에서 추출 (역 CDF 법, § 1.9)
- \(\beta^{(s)}\) 를 이산 조건부 \(p(\beta \mid \alpha^{(s)}, y)\) 에서 추출
- 두 값에 균등 난수 jitter (격자 간격 폭) 추가 — 이산을 연속으로 부드럽게
결과는 \((\alpha^{(s)}, \beta^{(s)})\) 의 1000 개 샘플. 교재 그림 3.3b 의 산점도.
2.11 그림 3.3a 의 등고선 플롯
0.05·0.15·…·0.95 × (모드 밀도) 의 등고선. 사후 \((\alpha, \beta)\) 의 양의 상관 구조 가 명확 — \(\alpha\) 가 크면 \(\beta\) 도 크다는 경향.
데이터가 \(\hat\alpha, \hat\beta\) 근처의 특정 조합을 선호할 뿐, 정확한 값 을 결정하지 못한다. \(\alpha\) 가 작으면 \(\beta\) 도 작게 조정 (곡선을 같은 위치에 유지). 두 모수가 데이터를 맞추기 위해 trade-off 관계.
이것이 모수 식별성 (identifiability) 의 약한 실패 — 모수들이 개별적으로는 잘 식별되지 않지만 결합으로는 잘 식별되는 경우. Part III 의 재매개변수화 (reparameterization) 로 개선 가능하지만, 2 모수에서는 격자가 여전히 작동.
2.12 LD50 — 사망률 50% 용량
정의.
\[ E\left(\frac{y_i}{n_i}\right) = \text{logit}^{-1}(\alpha + \beta x_i) = 0.5 \]
\(\text{logit}(0.5) = 0\) 이므로 \(\alpha + \beta x_i = 0\), 즉
\[ \text{LD50} = -\frac{\alpha}{\beta} \]
사후 시뮬레이션이 간단 — 각 \((\alpha^{(s)}, \beta^{(s)})\) 에서 \(\text{LD50}^{(s)} = -\alpha^{(s)}/\beta^{(s)}\) 를 계산, 히스토그램.
2.13 \(\beta \leq 0\) 인 경우의 해석 문제
LD50 은 \(\beta > 0\) (용량 증가가 사망률 증가로 이어짐) 에서만 의미. \(\beta < 0\) 이면 약물이 예방 효과 가 있는 셈 — LD50 의 의미 상실.
현재 실험에서는 — 1000 표본 모두 \(\beta > 0\). 따라서
\[ \Pr(\beta > 0 \mid y) > 1 - 1/1000 = 0.999 \]
약물이 해롭다는 확률 실질적 확실.
2.14 LD50 의 사후 요약
LD50 분포 (그림 3.4). \(\beta > 0\) 조건부로 LD50 히스토그램. 교재의 보고.
- Pr(β > 0) > 0.999 — 약물이 사망률 증가
- LD50 | β > 0 의 사후 분포 — 중앙값·95% 구간 제공
교재의 중요한 교훈.
“이 예는 주변 사후 평균이 항상 모수 추론의 좋은 요약은 아니다 라는 것을 보여준다. LD50 의 사후 평균은 용량-반응 관계가 음수인 경우를 포함하기 때문에, 일반적으로 우리는 LD50 의 사후 평균에 관심이 없다.” (교재)
\(\beta \approx 0\) 인 시뮬이 있으면 \(\text{LD50} = -\alpha/\beta\) 가 ±∞ 로 발산. 평균이 극단값에 오염된다. 조건부 \(\beta > 0\) 에서의 분포 가 훨씬 의미 있다.
이 교훈이 Ch.6 의 사후 점검에서도 반복 — “의미 있는 요약” 의 선택은 모수의 물리적 해석에 의존.
2.15 격자 계산의 실무 주의사항
교재의 경고.
“2 차원 격자 근사 적용 시 여러 실무적 고려사항. 너무 작은 영역에 정의된 격자는 사후의 중요한 특징을 놓칠 수 있다. 큰 영역에 넓은 간격 으로 정의된 격자도 점 사이의 중요한 특징을 놓칠 수 있다.”
실전 가이드.
- Rough 추정 으로 위치 파악
- 넉넉한 영역 설정 (MLE ± 3·SE 정도)
- 촘촘한 간격 (격자 100~200 포인트 이상)
- 꼬리 확인 — 격자 끝에서 밀도가 충분히 작은지
- 로그 스케일 계산 — 오버플로 방지
2.16 격자의 한계
\((\alpha, \beta)\) 는 2 차원이라 격자 \(100 \times 100 = 10^4\) 점으로 충분. 3 차원이면 \(10^6\), 5 차원이면 \(10^{10}\) — 실용적 한계.
“더 복잡한 모델은 Part III 의 계산 방법을 써야 한다.”
bioassay 는 “격자가 가능한 마지막 영역” 의 예시.
3 § 3.8 Ch.1~3 종합 — 베이즈 계산 5 단계 전략
교재의 알고리즘 표기.
- 가능도 작성 \(p(y \mid \theta)\) — \(\theta\) 와 무관한 요인은 생략
- 사후 밀도 작성 \(p(\theta \mid y) \propto p(\theta) \, p(y \mid \theta)\)
- 사전 정보가 잘 정립되면 \(p(\theta)\) 에 포함
- 아니면 약정보적 사전 또는 일시적 \(p(\theta) \propto \text{constant}\) — 나중에 추가 정보·구조 반영 가능
- 조잡한 모수 추정 — 시작점과 계산 비교 기준
- 사후 시뮬레이션 \(\theta^{(1)}, \ldots, \theta^{(S)}\) — 관심 함수의 사후 밀도 계산
- 예측량 시뮬레이션 \(\tilde{y}^{(1)}, \ldots, \tilde{y}^{(S)}\) — 각 \(\theta^{(s)}\) 에서 \(\tilde{y}^{(s)} \sim p(\tilde{y} \mid \theta^{(s)})\)
3.1 비켤레 모델에서의 어려움
“비켤레 모델에서 위 단계 4 는 어려울 수 있다. 복잡한 모델에서 사후 시뮬을 뽑는 다양한 방법이 개발 되었다 (Part III 에서 논의). 때때로 고차원 문제는 해석적·수치적 시뮬 방법의 결합 으로 해결 가능. \(\theta\) 가 1~2 성분이면 이전 절 bioassay 처럼 격자 계산 이 유효.” (교재)
3.2 복잡성 해결의 세 이유
교재가 제시하는 “다모수 모델의 닫힌 형태 부재” 가 실무 제약이 아닌 이유.
- 모수 적으면 시뮬이 해결 — bioassay 격자 같은 방법
- 계층적·조건부 구조 — Ch.5, Part III 의 전략
- 정규 근사 — Ch.4 주제, 정규 밖 모델에도 유용
Ch.3 이 Part I 의 실용 “단순 모델 편” 을 닫고 Ch.4 (점근) → Ch.5 (계층) → Part III (MCMC) 로 이어지는 다리.
4 § 3.9 참고문헌 주해
4.1 다변량 정규와 정규 추론
Box & Tiao (1973) — Ch.2 에서 정규 단변량/다변량과 평균 차·분산 비 등 관련 문제를 상세 처리. 당시 시뮬 방법이 편하지 않았기에 켤레 가족의 해석적 주변 사후 밀도 유도 에 큰 노력 투입.
Mardia, Kent, & Bibby (1979) — 다변량 정규의 수학적 특성 상세 (모든 주변·조건부가 정규).
4.2 Newcomb 데이터
Stigler (1977) — Simon Newcomb 광속 데이터와 실험 논의.
4.3 다항 모델과 사전
Good (1965), Fienberg (1977) — 다항 모델과 정보적/비정보적 사전의 체계적 논의. 로그선형 모델 확장은 Ch.16.
4.4 Bioassay
Racine et al. (1986) — Bioassay 데이터와 모델의 원전. 제약 산업에서 유용했던 간단한 베이즈 분석 여러 예.
5 § 3.10 선정 연습문제 풀이
5.1 Exercise 1 — Dirichlet 의 비율 주변
문제. \((y_1, \ldots, y_J) \sim \text{Multinomial}(\theta_1, \ldots, \theta_J)\), \(\theta \sim \text{Dirichlet}\). \(\alpha = \theta_1 / (\theta_1 + \theta_2)\) 의 사후 주변은?
풀이. Dirichlet 의 분할 특성 (aggregation property).
\((\theta_1, \theta_2, \theta_3 + \ldots + \theta_J) \sim \text{Dirichlet}(\alpha_1, \alpha_2, \alpha_3 + \ldots + \alpha_J)\).
\(\alpha_1, \alpha_2\) 의 상대 비율. \((\theta_1, \theta_2)\) 를 \(\theta_1 + \theta_2\) 로 정규화.
\[ \frac{\theta_1}{\theta_1 + \theta_2} \sim \text{Beta}(\alpha_1, \alpha_2) \]
이것이 Dirichlet 의 성질.
(b) \(y_1\) 을 \(\text{Bin}(y_1 + y_2, \alpha)\) 로 취급했을 때와 동일 증명.
사전 \(\alpha \sim \text{Beta}(\alpha_{1,\text{prior}}, \alpha_{2,\text{prior}})\), 가능도 \(p(y_1 \mid y_1 + y_2, \alpha) = \binom{y_1 + y_2}{y_1} \alpha^{y_1}(1-\alpha)^{y_2}\).
사후 \(\alpha \sim \text{Beta}(\alpha_{1,\text{prior}} + y_1, \alpha_{2,\text{prior}} + y_2)\).
Dirichlet 사후에서 \(\alpha\) 의 주변도 Beta(\(\alpha_{1,\text{prior}} + y_1, \alpha_{2,\text{prior}} + y_2\)) — 정확히 일치.
3 이상 범주의 다항 문제에서 특정 두 범주의 상대 비율만 관심이면, 나머지 범주를 무시하고 이항 모델 로 취급해도 된다. Dirichlet 의 수학적 우아함이 정당화.
1988 CBS 예제 (§ 3.4 심화) 에서 “Bush vs Dukakis 중 Bush 의 비율” \(\alpha = \theta_1/(\theta_1+\theta_2)\) 만 관심이면 \(\text{Beta}(728, 584)\) 사후. 기타/무응답 \(\theta_3\) 를 무시.
5.2 Exercise 2 — 1988 ABC 토론 전후 지지 변화
데이터 (표 3.2). 1988년 9월 25일 대선 토론 전후 ABC News 설문 각 639 명.
| 조사 | Bush | Dukakis | 무응답/기타 | 총 |
|---|---|---|---|---|
| 토론 전 | 294 | 307 | 38 | 639 |
| 토론 후 | 288 | 332 | 19 | 639 |
모델. 독립 다항 \(\alpha_j\) = 조사 \(j\) 에서 “선호 후보가 있는” 응답자 중 Bush 지지 비율. 즉 \(\alpha_j = \theta_{j,\text{Bush}} / (\theta_{j,\text{Bush}} + \theta_{j,\text{Dukakis}})\).
Ex 1 의 결과 적용 — 각 조사에서 \(\alpha_j \sim \text{Beta}(\text{Bush}_j + 1, \text{Dukakis}_j + 1)\) (균등 Dirichlet 사전).
- \(\alpha_1 \sim \text{Beta}(295, 308)\)
- \(\alpha_2 \sim \text{Beta}(289, 333)\)
\(\alpha_2 - \alpha_1\) 사후 분포 — 시뮬레이션.
- \(\alpha_1^{(s)} \sim \text{Beta}(295, 308)\)
- \(\alpha_2^{(s)} \sim \text{Beta}(289, 333)\)
- \(\text{diff}^{(s)} = \alpha_2^{(s)} - \alpha_1^{(s)}\)
관찰 — 토론 후 Bush 비율이 낮아진 경향. \(\Pr(\alpha_2 - \alpha_1 < 0 \mid y)\) 즉 “토론 후 Bush 가 상대적으로 감소” 가 얼마나 확실한지 계산.
수동 계산 — \(\alpha_1 \approx 295/603 = 0.489\), \(\alpha_2 \approx 289/622 = 0.464\). 차이 \(\approx -0.025\). 각 사후 표준편차 \(\approx \sqrt{0.5 \cdot 0.5 / 600} \approx 0.020\). 정규 근사로 \(\Pr(\alpha_2 - \alpha_1 < 0) \approx \Phi(-0.025/(0.020\sqrt{2})) \approx \Phi(-0.88) \approx 0.81\).
결론 — 81% 확률로 토론 후 Bush 비율이 감소 — 약한 증거. 결정적이지는 않다.
5.3 Exercise 3 — Behrens-Fisher 문제 (닭 칼슘)
문제. 통제 집단 32 마리 (\(\bar{y}_c = 1.013\), \(s_c = 0.24\)), 처치 집단 36 마리 (\(\bar{y}_t = 1.173\), \(s_t = 0.20\)).
비정보 사전 \(p(\mu_c, \mu_t, \log \sigma_c, \log \sigma_t) \propto 1\).
(a) 각 \(\mu\) 의 주변 사후.
§ 3.2 결과 적용 — 두 독립 정규 분석.
\[ \mu_c \mid y \sim t_{n_c - 1}(\bar{y}_c, s_c^2/n_c) = t_{31}(1.013, 0.24^2/32) \]
\[ \mu_t \mid y \sim t_{n_t - 1}(\bar{y}_t, s_t^2/n_t) = t_{35}(1.173, 0.20^2/36) \]
(b) \(\mu_t - \mu_c\) 의 사후.
두 독립 \(t\) 의 차 — Behrens-Fisher 문제. 닫힌 형태 복잡 → 시뮬레이션이 간단.
- \(\mu_c^{(s)} \sim t_{31}(1.013, 0.24^2/32)\)
- \(\mu_t^{(s)} \sim t_{35}(1.173, 0.20^2/36)\)
- \(\text{diff}^{(s)} = \mu_t^{(s)} - \mu_c^{(s)}\)
근사 95% 구간. 차이의 평균 \(= 1.173 - 1.013 = 0.160\). 차이의 표준편차 \(\approx \sqrt{0.24^2/32 + 0.20^2/36} \approx \sqrt{0.0018 + 0.0011} = \sqrt{0.0029} \approx 0.054\). 정규 근사 95% 구간 \(\approx [0.054, 0.266]\).
결론 — 처치 집단의 칼슘 유출이 더 높음 (약 16% 차이). 0 이 구간 밖 — 효과 존재의 강한 증거.
빈도주의 Behrens-Fisher 문제 — 두 모집단 정규에서 분산이 다른 경우 평균 차이의 정확한 분포가 복잡하다 (Welch 근사 등 사용).
베이즈에서는 — 각 모집단에 독립 사후를 세우고 시뮬레이션 차이로 끝. 수식적 복잡성이 시뮬로 소거 된다. Ch.3 의 큰 실용 이점.
5.4 Exercise 4 — 베타 차단제 2×2 테이블
문제. 통제 674 명 중 39 사망, 처치 680 명 중 22 사망. \(p_0, p_1\) = 각 사망률.
(a) 비정보 사전·사후 시뮬.
Beta(1, 1) 각각 — \(p_0 \sim \text{Beta}(40, 636)\), \(p_1 \sim \text{Beta}(23, 659)\).
사후 평균 — \(\hat{p}_0 = 40/676 \approx 0.059\), \(\hat{p}_1 = 23/682 \approx 0.034\).
(b) 교차비 (odds ratio) 의 사후.
\[ \text{OR} = \frac{p_1/(1-p_1)}{p_0/(1-p_0)} \]
시뮬 기반. \(p_0^{(s)}, p_1^{(s)}\) 추출 후 각 OR 계산.
\(\log(\text{OR})\) 이 더 자연스럽다 (대칭적 분포 기대). \(\log \text{OR} \approx \log(0.034/0.966) - \log(0.059/0.941) = -3.35 - (-2.77) = -0.58\). OR ≈ 0.56.
즉 처치군 사망 오즈가 통제군의 약 56% (약 44% 감소).
(c) 사전 민감도. Beta(1, 1) 대 Beta(0.5, 0.5) 대 Beta(0, 0) — 표본이 커서 (수백) 결과 비슷하지만, 희귀 사건 (40, 22 사망) 이라 사전이 약간 영향. Ch.5 에서 계층 모형으로 개선.
5.5 Exercise 5 — 반올림 측정
문제. 5 회 측정: 10, 10, 12, 11, 9 (파운드, 가장 가까운 정수로 반올림). 비정보 사전.
(a) 반올림 무시한 사후.
\(\bar{y} = 10.4, s^2 = 1.3\). § 3.2 결과 — \(\mu \mid y \sim t_4(10.4, 1.3/5) = t_4(10.4, 0.26)\).
(b) 반올림 고려한 정확한 사후.
진짜 측정 \(z_i \in (y_i - 0.5, y_i + 0.5]\). 가능도.
\[ p(y \mid \mu, \sigma^2) = \prod_{i=1}^5 [\Phi(y_i + 0.5 \mid \mu, \sigma) - \Phi(y_i - 0.5 \mid \mu, \sigma)] \]
비켤레 — 격자 또는 시뮬.
(c) 차이. 표본이 작고 (5) 표준편차가 작으면 반올림 오차 1 단위 가 상대적으로 큼 → 반올림 고려 사후가 더 넓고 분산 약간 더 큼.
(d) \((z_1 - z_2)^2\) 의 사후 평균. 원자료 복원 시뮬.
반올림 모델에서는 \(z_1, z_2\) 가 각각 독립적으로 \((9.5, 10.5], (9.5, 10.5]\) 에 존재. 같은 범위 → 차가 0 에 가까움. 그러나 \(\mu, \sigma^2\) 의 분포까지 반영하면 \((z_1 - z_2)^2\) 의 사후 평균이 0 이 아닌 양수.
반올림은 검열 (censoring) 의 특수 사례 — 정확한 값 대신 구간만 관측. 베이즈 접근 — 정확한 값 \(z_i\) 를 잠재 변수로 도입하고 데이터 증대 (data augmentation) — Ch.18 의 기본 전략.
단순 예에서는 격자 계산 또는 정규-정규 근사로 해결. 복잡한 검열은 MCMC.
6 네 절의 구조적 통합
| 절 | 주제 | 핵심 메시지 |
|---|---|---|
| § 3.7 | Bioassay 격자 | 2 모수 비켤레도 격자로 해결 |
| § 3.8 | 5 단계 전략 | Ch.1~3 실용 매뉴얼 |
| § 3.9 | 참고문헌 | 정규·다항·bioassay 의 지적 계보 |
| § 3.10 | 연습 | Ch.3 전체 기법 훈련 |
7 코드 예제 — Bioassay 격자 계산과 LD50 사후
7.1 Step 1: 순수 Python — 격자 계산과 1000 표본
import math
import random
random.seed(42)
# 데이터
x = [-0.86, -0.30, -0.05, 0.73]
n = [5, 5, 5, 5]
y = [0, 1, 3, 5]
def sigmoid(z):
return 1.0 / (1.0 + math.exp(-z))
def log_posterior(alpha, beta):
# 로그 비정규화 사후 (균등 사전)
lp = 0.0
for xi, ni, yi in zip(x, n, y):
p = sigmoid(alpha + beta * xi)
# 수치 안정을 위해 log(p), log(1-p) 를 조심스럽게
if p <= 0 or p >= 1:
return float("-inf")
lp += yi * math.log(p) + (ni - yi) * math.log(1 - p)
return lp
# 격자 설정
alpha_grid = [-5 + i * 0.15 for i in range(101)] # 101 points
beta_grid = [-10 + i * 0.5 for i in range(101)] # 101 points
# 로그 비정규화 사후 계산
log_post = [[log_posterior(a, b) for b in beta_grid] for a in alpha_grid]
# 최대값 차감 후 지수화
max_lp = max(max(row) for row in log_post)
unnorm = [[math.exp(log_post[i][j] - max_lp) for j in range(len(beta_grid))]
for i in range(len(alpha_grid))]
# 정규화
total = sum(sum(row) for row in unnorm)
post = [[unnorm[i][j] / total for j in range(len(beta_grid))]
for i in range(len(alpha_grid))]
# α 주변 — β 축으로 합
alpha_marg = [sum(post[i]) for i in range(len(alpha_grid))]
# 격자 표본 추출 (역 CDF 법)
def sample_discrete(probs, values):
cum = []
s = 0.0
for p in probs:
s += p
cum.append(s)
u = random.random()
for i, c in enumerate(cum):
if u <= c:
return values[i]
return values[-1]
S = 1000
alpha_samples = []
beta_samples = []
for _ in range(S):
# α 주변에서 추출
a_idx = sample_discrete(alpha_marg, range(len(alpha_grid)))
# β | α 조건부
beta_cond = post[a_idx]
total_cond = sum(beta_cond)
beta_cond = [p / total_cond for p in beta_cond]
b_idx = sample_discrete(beta_cond, range(len(beta_grid)))
# 격자 jitter
a_jit = alpha_grid[a_idx] + (random.random() - 0.5) * 0.15
b_jit = beta_grid[b_idx] + (random.random() - 0.5) * 0.5
alpha_samples.append(a_jit)
beta_samples.append(b_jit)
# LD50 — β > 0 인 표본에서
ld50 = [-a / b for a, b in zip(alpha_samples, beta_samples) if b > 0]
p_beta_pos = sum(1 for b in beta_samples if b > 0) / S
print(f"Pr(β > 0): {p_beta_pos:.4f}")
ld50.sort()
print(f"LD50 사후 중앙값 (β > 0 조건부): {ld50[len(ld50) // 2]:.3f}")
print(f"LD50 95% 구간: [{ld50[int(0.025 * len(ld50))]:.3f}, {ld50[int(0.975 * len(ld50))]:.3f}]")예상 출력 — \(\Pr(\beta > 0) = 1.000\) 에 가깝고 (모든 시뮬이 양수), LD50 중앙값 ≈ \(-0.1\) 근처, 95% 구간 \(\approx [-0.3, 0.1]\).
7.2 Step 2: scipy + Behrens-Fisher 닭 칼슘 예제
import numpy as np
from scipy import stats
np.random.seed(0)
# 통제 vs 처치
n_c, mean_c, s_c = 32, 1.013, 0.24
n_t, mean_t, s_t = 36, 1.173, 0.20
S = 10000
# 독립 t 사후에서 샘플
# μ ~ t_{n-1}(ȳ, s²/n)
mu_c = stats.t.rvs(df=n_c-1, loc=mean_c, scale=s_c/np.sqrt(n_c), size=S)
mu_t = stats.t.rvs(df=n_t-1, loc=mean_t, scale=s_t/np.sqrt(n_t), size=S)
diff = mu_t - mu_c
lo, hi = np.percentile(diff, [2.5, 97.5])
p_positive = (diff > 0).mean()
print(f"μ_t - μ_c 사후 평균: {diff.mean():.4f}")
print(f"μ_t - μ_c 95% 구간: [{lo:.4f}, {hi:.4f}]")
print(f"Pr(μ_t > μ_c | y) = {p_positive:.4f}")예상 출력 — 평균 차 ≈ 0.160, 95% 구간 [약 0.055, 0.265], \(\Pr(\mu_t > \mu_c) > 0.99\). 처치 집단의 칼슘 유출이 높다는 강한 증거.
8 관련 주제
Ch.3 의 다른 심화 포스트
- Ch.3 개요
- § 3.1~3.3 심화 — 주변화·정규 공동 사후
- § 3.4~3.6 심화 — 다항·다변량
Ch.1~2 심화 (선행)
- Ch.1 개요 + § 1.1~1.4 · § 1.5~1.8 · § 1.9·1.10·1.12
- Ch.2 개요 + § 2.1~2.4 · § 2.5~2.7 · § 2.8~2.11
Part I~V 전체
빈도주의 대응
- GLM 이항 모델 — 로지스틱 회귀의 GLM 관점
- GLM 이론 기초 — 지수족·링크 함수
- MLE · 점 추정
- Monte Carlo — 시뮬 기반 추론
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.3 (§ 3.7~3.10).
- Racine, A., Grieve, A. P., Fluhler, H., & Smith, A. F. M. (1986). Bayesian methods in practice: Experiences in the pharmaceutical industry. Applied Statistics, 35(2), 93–150.
- Mardia, K. V., Kent, J. T., & Bibby, J. M. (1979). Multivariate Analysis. Academic Press.
- Raftery, A. E. (1988). Inference for the binomial \(N\) parameter: A hierarchical Bayes approach. Biometrika, 75(2), 223–228.
- Heitjan, D. F. (1989). Inference from grouped continuous data: A review. Statistical Science, 4(2), 164–179.