Conditional Likelihoods — Further Results and Exercises

정규 축약 우도 · REML 유도 · Fieller-Creasy · 사영 행렬 · 초기하 반복 (McCullagh §7.7)

Ch.7 조건부 우도의 심화 결과와 연습문제를 재구성한다. 정규 모형의 축약 우도와 Bartlett (1936) 결과, 공간 공분산의 REML 유도, Fieller-Creasy 문제, 사영 행렬과 일반화역행렬, 초기하 반복 알고리즘, 순서형 분산 근사, Ille-et-Vilaine 확장 모형을 다룬다.

Statistics
GLM
저자

Kwangmin Kim

공개

2026년 04월 18일

1 이 장의 위치

McCullagh §7.7은 Ch.7 전체를 관통하는 19개의 연습문제를 제시한다. 단순 계산 문제가 아니라, 조건부 우도 이론의 핵심 확장과 연결을 담고 있다:

  • Exercises 7.1-7.4: 정규 모형에서 축약 우도(reduced likelihood) \(l^*\) 와 조건부 우도의 관계 — Bartlett (1936)의 고전적 결과
  • Exercise 7.5: 공간 공분산 추정에서 REML의 정확한 유도
  • Exercises 7.6-7.7, 7.19: Fieller-Creasy 문제 — 정규 평균의 비(ratio)에 대한 조건부 추론
  • Exercises 7.8-7.13: 사영 행렬, 일반화역행렬, 잔차 벡터의 동치성 — REML의 행렬 대수적 기초
  • Exercise 7.14: 비중심 초기하의 평균·분산 반복 알고리즘 (식 7.11, 7.12)
  • Exercises 7.15-7.17: 순서형 반응의 분산 근사와 추정 방정식의 성질
  • Exercise 7.18: Ille-et-Vilaine 식도암 자료의 확장 모형

이 포스트에서는 연습문제를 주제별로 묶어 풀이하고, 각 결과의 직관적 의미와 실무적 연결을 설명한다.

2 정규 모형의 축약 우도 (Exercises 7.1–7.4)

2.1 Exercise 7.1 — 완전충분통계량

\(Y_1, \ldots, Y_n \stackrel{\text{iid}}{\sim} N(\mu, \sigma^2)\) 에서 \(\mu = \mu_0\) 가 알려져 있을 때

\[ S_0 \equiv S(\mu_0) = \sum_{j=1}^{n} (Y_j - \mu_0)^2 \]

는 \(\sigma^2\) 에 대한 완전충분통계량(complete sufficient statistic)이다.

풀이: 로그우도를

\[ l(\sigma^2) = -\frac{S_0}{2\sigma^2} - \frac{n}{2}\log\sigma^2 + \text{const} \]

로 쓸 수 있다. \(S_0\) 가 \(\sigma^2\) 에 대한 지수족의 자연 통계량이므로 충분성은 인수분해 정리에서, 완전성은 \(S_0 / \sigma^2 \sim \chi^2_n\) 이 자연 지수족이므로 자동으로 따른다.

직관: \(\mu_0\) 가 알려져 있으면 각 관측값에서 \(\mu_0\) 를 빼고 제곱한 합이 \(\sigma^2\) 에 대한 모든 정보를 담는다. \(\mu\) 를 모르면 \(S(\bar{Y}) = \sum(Y_j - \bar{Y})^2\) 가 충분통계량이 되지만, 이때 자유도가 \(n-1\) 로 줄어든다 — 이것이 바로 REML의 직관적 기원이다.

2.2 Exercise 7.2 — 비중심 카이제곱과 조건부 로그우도

\(\mu\) 가 미지일 때 \(S(\mu_0)\) 의 분포는 비중심 카이제곱이다:

\[ S_0 / \sigma^2 \sim \chi^2_n(\lambda / \sigma^2), \qquad \lambda = n(\mu - \mu_0)^2. \]

조건부 로그우도 \(l_c(\mu, \sigma^2; \mu_0)\) 는 \(S_0\) 로 조건부한 것으로

\[ l_c(\mu, \sigma^2; \mu_0) = -\frac{1}{2\sigma^2}\{S(\mu) - S(\mu_0)\} - \left(\frac{n}{2} - 1\right)\log S(\mu_0) + \frac{\lambda}{2\sigma^2} - \log \sum_{r=0}^{\infty} \left(\frac{\lambda}{2\sigma^2}\right)^r \frac{1}{r!} B\!\left(\frac{n-1}{2}, r + \frac{1}{2}\right). \]

축약 함수(reduced function)는 \(\mu_0 = \mu\) 로 설정하여 \(\lambda = 0\) 으로 만든 것이다:

\[ l^*(\mu, \sigma^2) = l_c(\mu, \sigma^2;\, \mu) = -\frac{1}{2}(n-2)\log S(\mu). \]

직관 — 축약 우도의 의미

\(l^*\) 는 \(\sigma^2\) 가 완전히 사라진 함수이다. \(S(\mu)\) 는 \(\mu\) 에서의 잔차제곱합이므로, \(l^*\) 는 “잔차를 가장 작게 만드는 \(\mu\)”를 찾는다 — 이는 최소제곱과 같은 답을 주지만, 정당화 경로가 다르다(조건부 우도를 경유). 지수 \((n-2)\) 가 \((n-1)\) 이 아닌 이유는 조건부 취함으로써 자유도가 하나 더 줄어들기 때문이다.

2.3 Exercise 7.3 — 두 미분의 비교

\(l^*\) 의 \(\mu\) 에 대한 미분:

\[ U^* = \frac{\partial l^*}{\partial \mu} = \frac{(n-2)}{S(\mu)} \sum(y_j - \mu) = \frac{n(n-2)(\bar{y} - \mu)}{(n-1)s^2 + n(\bar{y} - \mu)^2}. \]

한편 \(l_c\) 의 미분을 \(\mu_0 = \mu\) 에서 평가하면

\[ \left.\frac{\partial l_c}{\partial \mu}\right|_{\mu_0 = \mu} = \frac{1}{\sigma^2} \sum(y_j - \mu), \]

이것은 \(\bar{y}\) 의 단조증가 함수이다.

직관: \(l_c\) 의 미분은 \(\sigma^2\) 에 의존하지만 \(\bar{y}\) 에 대해 단조이므로 \(\mu\) 의 MLE는 항상 \(\hat{\mu} = \bar{y}\) 이다. 반면 \(U^*\) 는 \(\sigma^2\) 에 의존하지 않지만, 분모에 \((\bar{y} - \mu)^2\) 항이 있어 \(\bar{y}\) 에 대해 비단조이다.

특수한 경우:

  • \(n = 1\): \(l^* = -\frac{1}{2} \cdot (-1) \cdot \log S(\mu)\) — \(S(\mu) = (y - \mu)^2\) 이므로 \(l^*\) 가 \(\mu = y\) 에서 최대가 아니라 최소가 된다. \(l^*\) 는 의미 있는 우도가 아니다.
  • \(n = 2\): \(l^* = 0\) (상수). 데이터에 무관하게 \(\mu\) 에 대한 정보가 전혀 없다.

이 기이한 행동은 \(l^*\) 가 엄밀한 로그우도가 아님을 보여준다 — Exercise 7.5에서 이 점이 공식적으로 확인된다.

2.4 Exercise 7.4 — 조건부 모멘트와 Bartlett (1936)

\(S(\mu)\) 가 주어졌을 때 \(Y_i\) 와 \(Y_j\) 는 조건부 비상관(conditionally uncorrelated)이며 \(\operatorname{var}(Y_i \mid S(\mu)) = S(\mu)/n\) 이다. 따라서

\[ \operatorname{var}(\bar{Y} \mid S(\mu)) = S(\mu) / n^2. \]

\(l^*\) 의 미분에 대한 조건부 모멘트:

\[ E(U^* \mid S(\mu)) = 0, \qquad \operatorname{var}(U^* \mid S(\mu)) = (n-2)^2 / S(\mu), \]

\[ -E\!\left(\frac{\partial^2 l^*}{\partial \mu^2} \;\middle|\; S(\mu)\right) = (n-2)^2 / S(\mu). \]

직관: 기대 정보 = 분산 등식이 성립한다. 이것은 \(l^*\) 가 우도 함수의 핵심 성질 — Bartlett 항등식(information identity) — 을 만족함을 의미한다. Bartlett (1936)가 보인 이 결과는 조건부 우도 이론의 초기 성과 중 하나이다.

3 공간 공분산과 REML (Exercise 7.5)

3.1 설정

공간 위치 \(s_1, \ldots, s_n\) 에서의 관측 \(Y\) 가

\[ E(Y) = X\beta, \qquad \operatorname{cov}(Y) = \Sigma(\theta; s) \]

를 따른다. \(\beta\) 는 장해 모수, \(\theta\) (공분산 모수)가 관심 대상이다.

3.2 \(\theta_0\) 에서의 충분통계량

\(\theta = \theta_0\) 에서 \(\beta\) 의 충분통계량은

\[ S_0 = X^T W_0 Y, \qquad W_0 = \Sigma(\theta_0; s)^{-1}. \]

직관: \(W_0\) 는 \(\theta_0\) 에서의 가중 행렬이고, \(S_0 = X^T W_0 Y\) 는 가중 최소제곱의 정규방정식 우변이다. \(\beta\) 에 대한 모든 정보가 이 통계량에 집약된다.

3.3 축약 함수 = REML

\(\theta_0 = \theta\) 로 설정하면 축약 함수가 된다:

\[ l^*(\beta, \theta) = -\frac{1}{2} Y^T \Sigma^{-1}(I - X(X^T \Sigma^{-1} X)^{-1} X^T \Sigma^{-1}) Y - \frac{1}{2}\log\det\Sigma + \frac{1}{2}\log\det(X^T \Sigma^{-1} X). \]

직관 — REML이 여기서 나온다

이 식은 REML(Restricted Maximum Likelihood) 의 정확한 형태이다.

  • 첫째 항: 잔차의 가중 제곱합 — \((I - P_W)Y\) 가 “회귀 잔차”이고 이를 가중 내적한 것
  • 둘째 항: \(-\frac{1}{2}\log\det\Sigma\) — 관측의 공분산 구조에 대한 벌점
  • 셋째 항: \(+\frac{1}{2}\log\det(X^T \Sigma^{-1} X)\) — 이 항이 REML을 일반 ML과 구분짓는다. \(\beta\) 를 적분(marginalize)한 대가로 생기는 보정항이다.

패턴 인식: Exercises 7.1-7.4에서 \((n-2)\) 지수가 등장한 것은 스칼라 버전의 REML 보정이었다. 여기서는 행렬 버전으로 확장된 것이다.

\(l^*\) 는 \(\beta\) 에 의존하지 않는다 — 장해 모수가 완전히 소거되었다.

그러나 \(X = I\) (즉 관측 수 = 모수 수)를 대입하면 \(l^*\) 가 로그우도 함수가 아님을 보일 수 있다. 이는 축약 함수가 일반적으로 진정한 우도가 아니라 근사적 우도(approximate likelihood)라는 사실을 확인한다.

4 Fieller-Creasy 문제 (Exercises 7.6–7.7, 7.19)

4.1 문제 설정

\(Y_1, Y_2\) 가 독립이고 \(Y_i \sim N(\mu_i, 1)\) 일 때 관심 모수는 평균의 비:

\[ \psi = \mu_2 / \mu_1. \]

이 문제는 §7.2.2 끝에서 언급된 Fieller-Creasy 문제의 구체적 전개이다.

4.2 Exercise 7.6 — 조건부 분포의 비유일성

\(Y_1 + \psi Y_2 = C\) 를 조건으로 하면 \(\psi Y_1 - Y_2\) 의 조건부 분포는 \(\psi\) 에만 의존하고 \(\mu_1\) 에는 의존하지 않는다. 그러나 \(Y_1\) 단독의 조건부 분포나 \(Y_2\) 단독의 조건부 분포는 서로 다른 “우도”를 유도한다.

직관: 장해 모수를 소거하는 충분통계량이 \(\psi\) 에 의존하므로 (\(S_\lambda(\psi) = Y_1 + \psi Y_2\)), 서로 다른 \(\psi\) 값에서 서로 다른 통계량으로 조건부를 취하게 된다. 이것이 조건부 우도의 근본적 한계 — 조건부 통계량의 \(\psi\)-의존성 — 이다.

4.3 Exercise 7.7 — 조건부 vs 주변 MLE

식 (7.2)를 사용한 조건부 MLE와 §7.2.1의 잔차 기반 주변 MLE를 비교하면:

  • 조건부 MLE: \(\hat{\psi}_c\) 는 \((Y_1 + \psi Y_2)\) 의 충분성에 기반
  • 주변 MLE: 잔차 \(\psi Y_1 - Y_2\) 의 분포에 기반

두 방법은 일반적으로 다른 추정값을 준다.

4.4 Exercise 7.19 — 편향 보정된 추정 방정식

\(l^*(\psi)\) 의 미분에서 조건부 기댓값을 빼면 편향 보정된 추정 방정식이 된다:

\[ \frac{\partial l^*(\psi)}{\partial \psi} - E\!\left(\frac{\partial l^*(\psi)}{\partial \psi} \;\middle|\; S_\lambda(\psi), \psi\right) = \frac{y_2 - \psi y_1}{1 + \psi^2} \cdot \frac{s_\lambda(\psi)}{1 + \psi^2}. \]

이 비편향 추정 방정식은 \(S_\lambda\) 가 \(\psi\) 에 의존하지 않을 때 원래의 (7.3)과 동일해진다.

직관: 조건부 통계량이 \(\psi\) 에 의존하면 \(l^*\) 의 미분이 조건부로 기댓값 0이 아닐 수 있다 — 즉 Bartlett 항등식이 깨진다. 편향 보정은 이 문제를 교정하지만, 결과적으로 추정 방정식이 더 복잡해진다.

5 사영 행렬과 잔차 동치성 (Exercises 7.8–7.13)

5.1 Exercise 7.8 — 행렬식 분해

직교 행렬 \(H = [H_1 : H_2]\) (\(H_1\): \(n \times (n-p)\), \(H_2\): \(n \times p\))에 대해

\[ \det \Sigma_0 = \frac{\det(H_1^T \Sigma_0 H_1)}{\det(H_2^T \Sigma_0^{-1} H_2)}. \]

직관: 양의 정부호 행렬의 행렬식을 직교 부분공간에서의 “부분 행렬식”으로 분해하는 공식이다. REML 로그우도에서 \(\log\det\Sigma\) 와 \(\log\det(X^T \Sigma^{-1} X)\) 항이 이 분해를 통해 연결된다.

5.2 Exercise 7.9 — DET 연산자

\(P_X = X(X^T X)^{-1} X^T\) 가 랭크 \(p\) 의 사영 행렬이고 \(\Sigma^* = (I - P_X)\Sigma_0(I - P_X)\) 일 때

\[ \log\operatorname{DET}(\Sigma^*) = \log\det\Sigma_0 + \log\det(X^T \Sigma_0^{-1} X) - 2\log\operatorname{DET}(X), \]

여기서 \(\operatorname{DET}(A)\) 는 \(A\) 의 비영 특이값의 곱이다.

직관: \(\Sigma^*\) 는 “회귀 잔차 공간에서의 공분산”이고 랭크가 \(n - p\) 이므로 일반적인 행렬식 대신 \(\operatorname{DET}\) 를 사용한다. 이 공식은 REML 로그우도를 잔차의 언어로 재표현하는 데 핵심적이다.

5.3 Exercises 7.10–7.13 — 주변 = 조건부 동치

\(Y = X\beta + Z\gamma + \epsilon\) 에서 \(R = (I - P_X)Y\) 의 분포는 \(\beta\) 에 의존하지 않는다. 정규 오차일 때

\[ \text{(주변 우도 기반 } R\text{)} = \text{(조건부 우도 기반 } P_X Y\text{)} \]

이 성립한다. 또한 가중 잔차 \(R_W = (I - P_W)Y\) (\(P_W = X(X^T W X)^{-1} X^T W\))에 대해

\[ \operatorname{cov}(R_W) = (I - P_W)\Sigma, \qquad \Sigma_W^{-} = \Sigma^{-1}(I - P_W), \]

그리고

\[ R^T \Sigma_I^{-} R = R_W^T \Sigma_W^{-} R_W = Y^T W Y - Y^T \Sigma^{-1} P_W Y. \]

직관 — 왜 이 동치가 중요한가

실무에서 REML을 구현할 때, “OLS 잔차 \(R\) 에 기반한 우도”와 “GLS 잔차 \(R_W\) 에 기반한 우도”가 동일하다는 사실은 구현의 자유도를 준다. 어느 잔차를 사용하든 같은 REML 추정값을 얻으므로, 수치적으로 안정적인 쪽을 선택하면 된다. lme4 등의 혼합모형 패키지는 이 동치를 활용하여 효율적인 REML 계산을 수행한다.

6 초기하 반복 알고리즘 (Exercise 7.14)

6.1 식 (7.11)-(7.12)의 반복 풀이

§7.3에서 비중심 초기하 \(H(m, s; \psi)\) 의 평균 \(\mu_{11}\) 과 분산 \(\kappa_2\) 가 동시에 다음 연립방정식을 만족한다:

\[ \psi = \frac{\mu_{11}\mu_{22} + \kappa_2}{\mu_{12}\mu_{21} + \kappa_2}, \]

\[ \kappa_2 = \frac{m_\cdot}{m_\cdot - 1} \left(\frac{1}{\mu_{11}} + \frac{1}{\mu_{12}} + \frac{1}{\mu_{21}} + \frac{1}{\mu_{22}}\right)^{-1}. \]

Exercise 7.14는 이 방정식을 반복적으로 푸는 알고리즘을 제시한다.

6.2 알고리즘

초기 추정값 \(\mu_{ij}^{(1)}\) 에서 시작하여:

(i) 분산 갱신:

\[ \kappa_2^{(r)} = \frac{m_\cdot}{m_\cdot - 1} \left(\frac{1}{\mu_{11}^{(r)}} + \frac{1}{\mu_{12}^{(r)}} + \frac{1}{\mu_{21}^{(r)}} + \frac{1}{\mu_{22}^{(r)}}\right)^{-1}. \]

(ii) 오즈비 갱신:

\[ \psi^{(r)} = \frac{\mu_{11}^{(r)}\mu_{22}^{(r)} + \kappa_2^{(r)}}{\mu_{12}^{(r)}\mu_{21}^{(r)} + \kappa_2^{(r)}}. \]

(iii) 평균 갱신 (지수족 성질 \(\partial \mu / \partial \theta = \kappa_2\) 활용):

\[ \mu_{11}^{(r+1)} = \mu_{11}^{(r)} + \kappa_2^{(r)} \cdot \{\log \psi - \log \psi^{(r)}\}. \]

직관: 단계 (iii)는 Newton-Raphson의 한 스텝이다. 비중심 초기하가 지수족이므로 \(\log \psi\) 가 자연 모수이고, \(\mu_{11}\) 이 평균 모수이며, \(\kappa_2 = \partial \mu_{11} / \partial(\log \psi)\) 가 정보(information)이다. 따라서 “현재 오즈비에서 목표 오즈비까지의 자연 모수 차이 \(\times\) 정보 = 평균의 보정량”이 된다.

초기값이 나쁘면 \(\mu_{11}\) 이 허용 범위 밖으로 나갈 수 있으므로 범위 제한이 필요하다. 수렴은 일반적으로 매우 빠르다 (Liao, 1988).

7 순서형 반응의 분산 근사 (Exercises 7.15–7.17)

7.1 Exercise 7.15 — 두 분산 공식의 비교

§7.5.2의 점근 분산 (7.27)에는 두 가지 표현이 있다:

\[ \operatorname{var}(\hat{\Delta}_c) \approx \zeta_\cdot \left\{\sum_{j=1}^{k-1} (\zeta_j + \zeta_{j+1}) d_j\right\}^{-1}, \]

\[ \operatorname{var}(\hat{\Delta}_c) \approx \frac{3(m_\cdot - 1)}{m_\cdot \zeta_\cdot} \left\{1 - \sum \left(\zeta_j / \zeta_\cdot\right)^3\right\}^{-1}. \]

두 공식은 동일하지 않지만 수치적으로 매우 유사하다. Table 7.3 데이터에서 차이는 약 0.5% 이다.

직관: 둘째 공식은 \(d_j\) 를 \(\widetilde{V}\) 의 대각 원소로 근사한 것이다. 도수가 극단적으로 불균등하지 않은 한 이 근사는 매우 정확하다. 실무에서는 어느 공식을 사용하든 결론이 동일하다.

7.2 Exercise 7.16 — 미사용 범주의 삭제

특정 범주에 응답이 전혀 없으면 (\(s_j = 0\)), 그 범주를 삭제해도 추정 방정식 (7.26)에 영향이 없다.

직관: \(s_j = 0\) 이면 \(S_{j-1} = S_j\) 이므로 연속된 두 절단점에서의 조건부 분포가 동일해지고, 해당 항이 텔레스코핑(telescoping)으로 상쇄된다. 이는 데이터가 결정하는 “유효 범주 수”만이 추정에 기여함을 의미한다.

7.3 Exercise 7.17 — 가중치의 전역 의존성

추정 방정식 (7.26)의 가중치 \(w_j^* = (\zeta_j + \zeta_{j+1}) / \zeta_\cdot\) 는 \(S_j\) 뿐만 아니라 전체 벡터 \(S\) 에 의존한다. 이것은 (7.25)의 “각 \(w_j^*\) 가 \((\psi, S_j)\) 에만 의존” 조건을 위반한다.

직관: 이론적 최적 가중치를 정확히 사용하려면 전체 주변 합계 정보가 필요하다. 그러나 가중치 선택이 효율에 미치는 영향은 작으므로 (§7.5.2 본문에서 “at worst a small loss of efficiency”), 이 위반은 실무적으로 무해하다. 단순 가중치 \(w_j^* = 1\) 또는 \(w_j^* = S_j(m_\cdot - S_j)\) 도 일치추정량을 주며, 효율 손실은 미미하다.

8 Ille-et-Vilaine 확장 모형 (Exercise 7.18)

8.1 구간별 상수 오즈비 모형

§7.4.3에서 세 모형을 비교했다. Exercise 7.18은 네 번째 모형을 추가한다: 75세 미만에서는 오즈비가 상수이고, 75세 이상에서는 다른 값을 가지는 모형이다.

\[ \log \psi_i = \begin{cases} \beta_0 & i = 1, \ldots, 5 \\ \beta_0 + \beta_2 & i = 6 \end{cases} \]

이 모형은 자유도 4 (모수 2개, 관측 6개)를 가진다.

직관: Table 7.2에서 75세 이상 층은 비암-고알코올 칸이 0이어서 오즈비가 \(\infty\) 이다. 또한 25-34세 층은 암-저알코올이 0이어서 역시 오즈비가 \(\infty\) 이다. 이 두 극단 층이 모형 (ii)의 이탈도를 주로 부풀리고 있다.

75세 이상을 분리하면 이탈도가 모형 (iii)의 선형 경향보다 더 많이 감소할 가능성이 있다. 그러나 이 모형은 데이터를 보고 선택한 것이므로 유의수준에 선택 효과(selection effect) 보정이 필요하다 — 이탈도 감소량 자체가 아니라, “가능한 모든 구간 분할 중 최선을 선택했다”는 사실을 고려해야 한다.

9 Python 구현

9.1 축약 우도와 REML (Exercises 7.1–7.5)

코드 
import numpy as np
from scipy.optimize import minimize_scalar

# ── Exercise 7.1-7.3: 정규 모형의 축약 우도 ──────
np.random.seed(42)
n = 10
mu_true, sigma_true = 3.0, 2.0
Y = np.random.normal(mu_true, sigma_true, n)

def S(mu, Y):
    return np.sum((Y - mu)**2)

def l_star(mu, Y, n):
    """축약 우도 l*(mu)"""
    return -0.5 * (n - 2) * np.log(S(mu, Y))

# l*의 미분 U*
def U_star(mu, Y, n):
    ybar = np.mean(Y)
    s2 = np.var(Y, ddof=1)
    return n * (n - 2) * (ybar - mu) / ((n - 1) * s2 + n * (ybar - mu)**2)

# MLE: l*를 최대화
res = minimize_scalar(lambda mu: -l_star(mu, Y, n), bounds=(0, 6), method='bounded')
mu_star = res.x

print(f"표본 평균: {np.mean(Y):.4f}")
print(f"l* MLE:    {mu_star:.4f}")
print(f"(두 값은 동일해야 한다)")

# ── Exercise 7.4: 조건부 모멘트 확인 (시뮬레이션) ─
n_sim = 100000
S_target = S(mu_true, Y)  # 고정된 S(mu) 값

# 조건부 분포에서 샘플링: Y | S(mu) = S_target
# chi^2_{n}(lambda/sigma^2) 에서 S_0 = S_target 조건부
ybar_samples = []
for _ in range(n_sim):
    Y_sim = np.random.normal(mu_true, sigma_true, n)
    S_sim = S(mu_true, Y_sim)
    if abs(S_sim - S_target) < S_target * 0.01:  # 근사 조건부
        ybar_samples.append(np.mean(Y_sim))

if len(ybar_samples) > 100:
    ybar_arr = np.array(ybar_samples)
    print(f"\n조건부 var(Y_bar | S): {np.var(ybar_arr):.4f}")
    print(f"이론값 S/(n^2):        {S_target / n**2:.4f}")

# ── Exercise 7.5: REML 비교 ──────────────────────
from scipy.optimize import minimize

# 단순 예: Y ~ N(X*beta, sigma^2 * I), X = [1]
# REML vs ML for sigma^2
n = 20
X = np.ones((n, 1))
beta_true = 5.0
sigma2_true = 4.0
Y = np.random.normal(beta_true, np.sqrt(sigma2_true), n)

def neg_ml(log_sigma2):
    sigma2 = np.exp(log_sigma2)
    beta_hat = np.mean(Y)
    resid = Y - beta_hat
    return 0.5 * n * np.log(sigma2) + 0.5 * np.sum(resid**2) / sigma2

def neg_reml(log_sigma2):
    sigma2 = np.exp(log_sigma2)
    p = X.shape[1]
    beta_hat = np.mean(Y)
    resid = Y - beta_hat
    # REML = ML + 0.5 * log det(X^T Sigma^{-1} X)
    return (0.5 * n * np.log(sigma2) + 0.5 * np.sum(resid**2) / sigma2
            + 0.5 * np.log(n / sigma2))  # = 0.5*log(n) - 0.5*log(sigma2)

res_ml = minimize(neg_ml, x0=0.0)
res_reml = minimize(neg_reml, x0=0.0)

sigma2_ml = np.exp(res_ml.x[0])
sigma2_reml = np.exp(res_reml.x[0])

print(f"\nML  sigma^2 = {sigma2_ml:.4f}  (n 나눗셈)")
print(f"REML sigma^2 = {sigma2_reml:.4f}  (n-1 나눗셈)")
print(f"표본분산 s^2 = {np.var(Y, ddof=1):.4f}")

9.2 초기하 반복 알고리즘 (Exercise 7.14)

코드 
import numpy as np
from math import comb, exp, log

def hyper_exact(m1, m2, s, psi):
    """비중심 초기하의 정확한 평균과 분산"""
    lo, hi = max(0, s - m2), min(m1, s)
    log_terms = []
    for t in range(lo, hi + 1):
        lt = np.log(comb(m1, t)) + np.log(comb(m2, s - t)) + t * np.log(psi)
        log_terms.append(lt)
    log_terms = np.array(log_terms)
    max_lt = np.max(log_terms)
    log_P0 = max_lt + np.log(np.sum(np.exp(log_terms - max_lt)))
    vals = np.arange(lo, hi + 1, dtype=float)
    probs = np.exp(log_terms - log_P0)
    mu = np.sum(vals * probs)
    var = np.sum(vals**2 * probs) - mu**2
    return mu, var

# ── Exercise 7.14: 반복 알고리즘 ──────────────────
m1, m2, s = 3, 4, 3
psi_target = 6.0  # log odds ratio ~ 1.79

# 초기값: 독립 가정 (psi=1)
mu11 = m1 * s / (m1 + m2)
m_dot = m1 + m2

print("=== 반복 알고리즘 (Exercise 7.14) ===")
print(f"목표 psi = {psi_target:.3f}")

for r in range(10):
    mu12 = m1 - mu11
    mu21 = s - mu11
    mu22 = m2 - mu21

    # (i) kappa2
    inv_sum = 1/mu11 + 1/max(mu12, 1e-10) + 1/max(mu21, 1e-10) + 1/max(mu22, 1e-10)
    kappa2 = (m_dot / (m_dot - 1)) / inv_sum

    # (ii) psi^{(r)}
    psi_r = (mu11 * mu22 + kappa2) / (mu12 * mu21 + kappa2)

    # (iii) mu11 갱신
    delta_theta = log(psi_target) - log(max(psi_r, 1e-10))
    mu11_new = mu11 + kappa2 * delta_theta

    # 범위 제한
    mu11_new = max(max(0, s - m2) + 0.01, min(min(m1, s) - 0.01, mu11_new))

    print(f"  iter {r}: mu11 = {mu11:.4f}, psi = {psi_r:.4f}, kappa2 = {kappa2:.4f}")

    if abs(mu11_new - mu11) < 1e-8:
        mu11 = mu11_new
        break
    mu11 = mu11_new

# 정확한 값과 비교
mu_exact, var_exact = hyper_exact(m1, m2, s, psi_target)
print(f"\n반복 결과: mu11 = {mu11:.6f}")
print(f"정확한 값: mu11 = {mu_exact:.6f}")
print(f"반복 결과: kappa2 = {kappa2:.6f}")
print(f"정확한 값: var   = {var_exact:.6f}")

10 R 구현

코드 
# ── Exercise 7.1-7.4: 축약 우도 ──────────────────
set.seed(42)
n <- 10
Y <- rnorm(n, mean = 3, sd = 2)

S_mu <- function(mu) sum((Y - mu)^2)
l_star <- function(mu) -0.5 * (n - 2) * log(S_mu(mu))

# 최적화
opt <- optimize(l_star, interval = c(0, 6), maximum = TRUE)
cat("l* MLE:", opt$maximum, "\n")
cat("표본평균:", mean(Y), "\n")

# ── Exercise 7.5: ML vs REML ─────────────────────
library(lme4)

# 단순 예: 랜덤 절편 모형
set.seed(123)
n_groups <- 20
n_per <- 5
group <- rep(1:n_groups, each = n_per)
u <- rnorm(n_groups, sd = 2)
y <- 5 + u[group] + rnorm(n_groups * n_per, sd = 1)

# ML vs REML
fit_ml   <- lmer(y ~ 1 + (1 | group), REML = FALSE)
fit_reml <- lmer(y ~ 1 + (1 | group), REML = TRUE)

cat("\nML   sigma_u:", sqrt(VarCorr(fit_ml)$group[1]), "\n")
cat("REML sigma_u:", sqrt(VarCorr(fit_reml)$group[1]), "\n")
cat("ML   sigma_e:", sigma(fit_ml), "\n")
cat("REML sigma_e:", sigma(fit_reml), "\n")

# ── Exercise 7.14: 초기하 반복 ───────────────────
hyper_iter <- function(m1, m2, s, psi_target, tol = 1e-8, max_iter = 50) {
  m_dot <- m1 + m2
  mu11 <- m1 * s / m_dot  # 초기값

  for (r in seq_len(max_iter)) {
    mu12 <- m1 - mu11
    mu21 <- s - mu11
    mu22 <- m2 - mu21

    inv_sum <- 1/mu11 + 1/mu12 + 1/mu21 + 1/mu22
    kappa2 <- (m_dot / (m_dot - 1)) / inv_sum

    psi_r <- (mu11 * mu22 + kappa2) / (mu12 * mu21 + kappa2)
    delta <- log(psi_target) - log(psi_r)
    mu11_new <- mu11 + kappa2 * delta

    # 범위 제한
    lo <- max(0, s - m2) + 0.001
    hi <- min(m1, s) - 0.001
    mu11_new <- max(lo, min(hi, mu11_new))

    if (abs(mu11_new - mu11) < tol) break
    mu11 <- mu11_new
  }

  list(mu11 = mu11, kappa2 = kappa2, psi = psi_r, iterations = r)
}

result <- hyper_iter(3, 4, 3, psi_target = 6.0)
cat("\nExercise 7.14 결과:\n")
cat("  mu11 =", result$mu11, "\n")
cat("  kappa2 =", result$kappa2, "\n")

# 정확한 값 (BiasedUrn)
if (requireNamespace("BiasedUrn", quietly = TRUE)) {
  exact_mean <- BiasedUrn::meanFNCHypergeo(3, 4, 3, 6.0)
  exact_var  <- BiasedUrn::varFNCHypergeo(3, 4, 3, 6.0)
  cat("  정확 mu11 =", exact_mean, "\n")
  cat("  정확 var  =", exact_var, "\n")
}

11 핵심 요약

Ch.7 심화 결과의 세 가지 주제

  • 축약 우도와 REML: 정규 모형에서 \(l^*(\mu) = -\frac{1}{2}(n-2)\log S(\mu)\) 는 \(\sigma^2\) 를 소거한 축약 우도이다. 다변량으로 확장하면 REML의 정확한 형태 \(l^* = -\frac{1}{2}Y^T\Sigma^{-1}(I - P_W)Y - \frac{1}{2}\log\det\Sigma + \frac{1}{2}\log\det(X^T\Sigma^{-1}X)\) 가 된다. 이 함수는 \(\beta\) 에 의존하지 않지만 엄밀한 우도는 아니다.
  • Fieller-Creasy: 충분통계량이 관심 모수 \(\psi\) 에 의존할 때 조건부 우도의 한계가 드러난다. 편향 보정이 가능하지만 추정 방정식이 복잡해진다.
  • 초기하 반복: Exercise 7.14의 알고리즘은 지수족의 \(\partial\mu/\partial\theta = \kappa_2\) 성질을 활용한 Newton-Raphson 스텝이며, 비중심 초기하의 평균과 분산을 수치적으로 빠르게 구한다.

Subscribe

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