1 개요 — §9.4 가 필요한 이유
§9.2, §9.3 에서 다룬 준스코어 함수 \(U(\beta; y) = D^\top V^{-1}(Y - \mu)\) 는 실은 더 큰 이론 — 추정함수 이론(theory of estimating functions) — 의 특수한 사례다. McCullagh & Nelder §9.4 는 이 관계를 뒤집어 묻는다.
근본 질문: “우도가 없을 때, 또는 우도 조건이 너무 강할 때, \(\theta\) 의 영평균 함수 \(g(Y;\theta)\) 를 어떻게 구성하고 결합해야 가장 효율적인 추정량을 얻는가?”
즉 준스코어는 \((Y-\mu)\) 라는 기본 영평균 함수로부터 \(D^\top V^{-1}\) 가중치로 결합된 선형 추정함수 클래스 내의 최적 결합임을 §9.4 는 증명한다.
1.1 왜 “추정함수”를 별도 이론으로 다루는가
실무에서 우리는 우도보다 훨씬 약한 정보만 가지고 추정을 해야 할 때가 많다.
- 시계열·종단 데이터: 완전한 결합 우도는 복잡하나, 마팅게일성(조건부 평균이 0)은 쉽게 확인된다.
- 부수 모수(nuisance parameter) 제거: 관심 모수에 대한 직접적 추정식을 세울 수 있다면 부수모수를 조건화로 제거 가능 (예: 매칭 짝 데이터의 부수 구간 모수).
- 기하학적/구조적 제약: “원 위의 점”, “타원 궤도” 같은 구조적 모형은 확률모형을 완전히 명시하지 않아도 영평균 관계식을 쉽게 적을 수 있다.
추정함수 이론은 이런 상황에서 “어떤 함수를 0 으로 놓고 풀지” 를 원리적으로 정해주는 프레임워크다.
2 추정함수의 정의
관측 데이터 \(Y\) 와 모수 \(\theta\) 의 함수 \(g(Y;\theta)\) 가 모든 \(\theta\) 값에서 기댓값 0 을 가지면, 즉 \[ E_\theta[g(Y;\theta)] = 0, \quad \forall \theta \in \Theta \] 를 만족하면 \(g\) 를 \(\theta\) 에 대한 추정함수(estimating function) 라 한다. 모수 개수 \(p\) 개만큼의 방정식 \(g(\widehat{\theta};Y)=0\) 을 풀어 추정량 \(\widehat{\theta}\) 을 얻는다.
2.1 스코어 함수와의 차이
스코어 함수 \(U(\theta) = \partial \ell/\partial \theta\) 는 Bartlett 항등식에 의해 \(E[U]=0\) 이므로 추정함수의 특수한 예다. 그러나 추정함수는 Bartlett 두 번째 항등식(\(\operatorname{var}(U) = -E[\partial U/\partial\theta]\)) 을 요구하지 않는다. 더 나아가 \(g\) 의 고차 누적률이 \(\theta\) 에 의존해도 된다. 즉 \(g(Y;\theta)\) 는 pivotal statistic 일 필요가 없다.
2.2 직관: “영평균 함수” 의 자유도
왜 “영평균”만 요구하고 그 외는 자유롭게 두는가. 핵심 아이디어는 다음 동치관계다.
\[ g(\widehat{\theta}; Y) = 0 \; \Longleftrightarrow \; \widehat{\theta} \text{ 는 } g \text{ 의 영점} \]
만일 참값 \(\theta_0\) 에서 \(E[g]=0\) 이고, 표본이 증가하면서 \(g\) 가 자신의 기댓값에 수렴하면(대수의 법칙), \(g\) 의 영점은 점근적으로 참값 \(\theta_0\) 에 수렴한다. 이것이 일치성(consistency) 의 본질이며, 우도의 형태가 맞든 틀리든 상관없이 성립한다.
반사실 — \(E[g] \ne 0\) 이면? \(g\) 의 평균이 상수 \(c \ne 0\) 으로 수렴하면, 영점 방정식 \(g(\widehat\theta; Y) = 0\) 의 해는 “\(\widehat\theta\) 에서 \(g\) 가 \(-c\) 를 상쇄하는 값” 을 가리키므로 참값 \(\theta_0\) 에서 체계적으로 벗어난다. 즉 영평균 조건은 일치성의 필수 전제 이며, 이 조건 없이는 어떤 \(g\) 로도 \(\theta_0\) 을 목표할 수 없다.
2.3 Y − μ(β) 는 가장 단순한 추정함수
§9.2 의 GLM 맥락에서 \(E[Y] = \mu(\beta)\) 이므로 \[ g(Y;\beta) = Y - \mu(\beta) \] 는 \(n\)-차원 벡터값 추정함수다. 문제는 이 \(n\)-벡터를 \(p\)-벡터로 압축하면서 정보 손실을 최소화해야 한다는 것이다. §9.4.2 가 답하는 질문이 바로 이것이다.
3 기본 추정함수와 조건부 무편향성
3.1 시간 순서가 있을 때 — 마팅게일 구조
관측이 시간 순서대로 \(Y_1, Y_2, \ldots, Y_n\) 로 들어온다고 하자. 각 시점 \(t\) 에서 과거까지의 정보 \(Y_{(t)} = (Y_1, \ldots, Y_t)\) 에 의존하는 기본 추정함수(elementary estimating function) \(g_t(Y_{(t)}; \theta)\) 가 \[ E\{g_t(Y_{(t)}; \theta) \mid Y_{(t-1)}\} = 0 \tag{마팅게일 조건} \] 을 만족하면, 누적합 \(\sum_{s \leq t} g_s\) 는 마팅게일(martingale) 을 이룬다.
마팅게일은 “다음 번 이득의 기대값이 지금까지의 이력을 알 때 0” 인 확률과정이다. \(g_t\) 를 이 ‘이득’ 으로 보면:
- \(g_t\) 자체는 양의 값일 수도 음의 값일 수도 있다 (랜덤).
- 그러나 과거를 모두 아는 상태에서의 예측값은 0 이다.
- 즉, \(g_t\) 는 “이미 알려진 정보로는 예측 불가능한 순수 놀라움(innovation)” 이다.
추정함수가 마팅게일이 되면, 여러 시점의 \(g_t\) 를 더해도 편향이 누적되지 않는다. 이것이 종단 데이터·시계열에서 추정함수 이론이 강력한 이유다 (Godambe & Heyde, 1987).
3.2 회귀 맥락 — 조건화 집합 \(A_i\)
회귀·시계열을 통합하기 위해 일반적으로 조건화 집합 \(A_i \equiv A_i(y; \theta)\) 를 도입하여 \[ E\{g_i(Y;\theta) \mid A_i\} = 0 \tag{9.13'} \] 을 요구한다. 이때 \(\{A_i\}\) 는 중첩조건(nested) \(A_{i-1} \subseteq A_i\) 을 만족해야 한다.
| 맥락 | \(A_i\) 의 의미 |
|---|---|
| 선형회귀 | \(A_i = A\) = 공변량 집합 (모든 \(i\) 에 공통) |
| 시계열 | \(A_i = Y_{(i-1)}\) = 시점 \(i-1\) 까지의 이력 |
| 일반 규칙 | \(\dim A_i\) 를 최대로 취한다 — 조건부 정보를 많이 쓸수록 효율이 올라간다 |
3.3 왜 “조건부” 영평균이 중요한가
무조건부 \(E[g_i]=0\) 만으로는 공변량 의존성이 숨어버릴 수 있다. 예컨대 시계열에서 \(g_t = Y_t - \theta Y_{t-1}\) 은 \(Y_{t-1}\) 에 조건부일 때만 영평균이며, 이 조건부성을 통해 \(Y_{t-1}\) 이 자연스럽게 가중치 함수 로 들어온다. §9.4.3 Avebury 예제에서도 조건부 기댓값과 무조건부 기댓값을 혼동하면 계수행렬 \(D\) 의 rank 가 떨어져서 식별불가가 발생한다.
4 추정함수의 최적 결합 — 핵심 공식
4.1 설정
\(n\) 개의 기본 추정함수 \(g_1(Y;\theta), \ldots, g_n(Y;\theta)\) 가 모두 조건부 영평균이고, 조건부로 독립 이라고 하자 (즉 \(V = \operatorname{cov}(g \mid A)\) 는 대각행렬).
이들을 선형결합하여 \(p\)-차원 추정방정식 하나로 줄이려면 \(n \times p\) 가중 행렬이 필요하다. McCullagh & Nelder 는 최적 가중 행렬 을 다음과 같이 정의한다.
\(n \times p\) 행렬 \(D\) 의 \((i, r)\)-성분은 \[ D_{ir} = -E\left\{\frac{\partial g_i(Y;\theta)}{\partial \theta_r} \,\Big|\, A_i\right\} \tag{9.14} \] 즉 기본 추정함수 \(g_i\) 의 \(\theta_r\) 에 대한 도함수의 조건부 음의 기댓값 이다.
조건부 공분산 \(V = \operatorname{cov}(g \mid A)\) 와 계수행렬 \(D\) 를 이용하여 \[ U(\theta; y) = D^\top V^{-1} g(Y; \theta) \tag{9.15} \] 를 \(\theta\) 에 대한 최적 추정함수 로 정의한다. 추정량 \(\widehat{\theta}\) 는 \(U(\widehat{\theta}; y)=0\) 의 해다.
\(D_{ir}\) 과 \(V\) 가 \(A_i\) 의 함수이고 \(g_i\) 가 조건부 영평균이므로, 반복기댓값 법칙 에 의해 \(U\) 도 무조건부 영평균이다.
4.2 왜 하필 \(D^\top V^{-1} g\) 인가 — 가중최소제곱의 일반화
\(g\) 자체가 잡음이 큰 영평균 벡터, \(D\) 가 각 관측이 모수에 반응하는 민감도(신호), \(V\) 가 잡음 분산 이라고 보면:
\[ D^\top V^{-1} g = \underbrace{D^\top}_{\text{감도(신호)}} \cdot \underbrace{V^{-1}}_{\text{잡음 역가중}} \cdot \underbrace{g}_{\text{영평균 잔차}} \]
이것은 정확히 가중최소제곱의 정규방정식 형태다. 분산이 큰 \(g_i\) 의 기여를 작게 하고, 모수 민감도가 큰 \(g_i\) 의 기여를 크게 한다 — 정보량에 비례한 결합.
4.3 선형변환 불변성
기본 추정함수를 \(g^* = B g\) (단, \(B\) 는 \(n \times n\) 가역행렬, \(A\) 의 함수) 로 바꾸면 \[ V^* = B V B^\top, \quad D^* = B D \] 가 되므로, \[ D^{*\top}(V^*)^{-1} g^* = (BD)^\top (B V B^\top)^{-1} (Bg) = D^\top B^\top B^{-\top} V^{-1} B^{-1} B g = D^\top V^{-1} g. \]
이 불변성은 실무적으로 매우 중요하다. 기본 추정함수를 \(Y_t - \theta Y_{t-1}\) 로 쓰든 \(Y_t/\theta - Y_{t-1}\) 로 쓰든(둘 다 영평균), 최적 결합 \(U = D^\top V^{-1} g\) 는 동일한 추정량을 준다. 즉 “어떤 단위로 영평균 함수를 쓰든 최종 추정량은 같다” — 수학적 우아함과 실무적 안정성을 동시에 보장한다.
4.4 AR(1) 예제 — 선형 \(g\) 에서 이차 \(U\) 가 나오는 구조
자기회귀 모형 \[ Y_t = \theta Y_{t-1} + \epsilon_t, \quad \epsilon_t \overset{\text{iid}}{\sim} N(0, 1) \] 에서 기본 추정함수는 \[ g_t = Y_t - \theta Y_{t-1}, \qquad E\{g_t \mid Y_{t-1}\} = 0. \]
- \(V_{tt} = \operatorname{var}(g_t \mid Y_{t-1}) = 1\) (대각 = 1).
- \(D_{t} = -E\{\partial g_t/\partial\theta \mid Y_{t-1}\} = Y_{t-1}\).
따라서 \[ U(\theta; y) = \sum_t Y_{t-1} \cdot (Y_t - \theta Y_{t-1}). \]
\(g_t\) 는 \(Y\) 에 대해 선형 이지만 최적 결합 \(U\) 는 \(Y_{t-1} Y_t\) 항 때문에 \(Y\) 에 대해 이차식 이다. 이는 로그우도 도함수와도 일치한다. 대체 선택 \(g_t^* = Y_t/\theta - Y_{t-1}\) 을 써도(선형변환 \(B_{tt} = 1/\theta\)) 같은 \(U\) 를 얻는다 — 불변성의 구체적 확인.
4.5 점근 분산 — \((D^\top V^{-1} D)^{-1}\)
§9.2.3 와 같은 Taylor 전개에 의해 \[ \operatorname{cov}(\widehat{\theta}) \approx (D^\top V^{-1} D)^{-1}. \]
자기회귀 모형처럼 \(D\) 가 데이터 \(Y_{t-1}\) 를 포함하는 경우, \(D\) 와 \(V\) 의 기댓값이 아닌 관측값(observed) 을 대입하는 쪽이 수치적으로 훨씬 안정적이다. Avebury 예제에서도 마찬가지 — 관측된 각도 \(\cos\epsilon_i\) 의 실제값이 식별에 필수적이다.
5 예제 1 — Fieller-Creasy 문제
5.1 문제 설정
관측이 쌍으로 들어온다: \(Y_i = (Y_{i1}, Y_{i2}), \; i=1,\ldots,n\). 독립이고 \[ E[Y_{i1}] = \mu_i, \quad E[Y_{i2}] = \mu_i/\theta, \quad \operatorname{var}(Y_{i1}) = \operatorname{var}(Y_{i2}) = \sigma^2. \] 관심 모수는 두 평균의 비 \(\theta = E[Y_{i1}]/E[Y_{i2}]\) 이며, \(\mu_i\) 는 부수 모수(nuisance) 다.
5.2 왜 이 문제가 어려운가
\(\mu_i\) 가 \(n\) 개 있으므로 표본이 늘어도 부수 모수 수가 함께 늘어난다 (Neyman-Scott 문제). 우도 기반 접근은 각 \(\mu_i\) 를 소거해야 하며 이는 profile likelihood 를 계산하게 만든다. 추정함수 접근은 \(\mu_i\) 를 건드리지 않고 \(\theta\) 만 추정 할 수 있음을 보여준다.
5.3 기본 추정함수
비율 관계 \(E[Y_{i2}] \theta = E[Y_{i1}]\) 를 이용하여 \[ g_i(Y_i; \theta) = Y_{i1} - \theta Y_{i2}. \] \(E[g_i] = \mu_i - \theta \cdot \mu_i/\theta = 0\) 이므로 \(\mu_i\) 와 무관하게 영평균이다.
\(\operatorname{var}(g_i) = \sigma^2(1+\theta^2)\), \(\;\partial g_i/\partial\theta = -Y_{i2}\).
5.4 계수행렬 D 구하기 — 조건화 선택이 핵심
조건화 집합을 \(A_i(\theta) = Y_{i2} + \theta Y_{i1}\) 로 선택한다 (아래에서 이 선택의 의미를 설명). 그러면 \[ -E\{\partial g_i/\partial\theta \mid A_i\} = E\{Y_{i2} \mid Y_{i2} + \theta Y_{i1}\}. \] 계산하면 (다변량 정규 조건부 평균 공식): \[ D_i = \frac{Y_{i2} + \theta Y_{i1}}{1 + \theta^2}. \]
5.5 최적 추정함수
\[ U(\theta) = \sum_i \frac{D_i \cdot g_i}{V_i} = \sum_{i=1}^n \frac{(Y_{i2} + \theta Y_{i1})(Y_{i1} - \theta Y_{i2})}{\sigma^2 (1+\theta^2)^2}. \]
이것은 §7.3 의 조건부 로그우도 스코어와 완전히 동일하다. 즉 조건화 집합을 적절히 고르면 추정함수 이론이 조건부 우도(conditional likelihood) 를 자동으로 재현 한다.
5.6 왜 두 개의 해가 있는가
\(U(\theta) = 0\) 은 일반적으로 두 해를 갖는다: 하나는 \(\widehat{\theta}\), 다른 하나는 \(-1/\widehat{\theta}\).
- 첫 번째: 로그우도의 최대 에 대응.
- 두 번째: 로그우도의 최소 에 대응 (허해).
정보량이 작으면 \(\widehat{\theta}\) 의 정규근사는 신뢰할 수 없다. 그 대신 스코어 기반 신뢰집합 \[ \{\theta : |U(\theta)/i^{1/2}(\theta)| \leq k^*_{\alpha/2}\}, \quad i(\theta) = \sum \frac{(y_{i2}+\theta y_{i1})^2}{\sigma^2(1+\theta^2)^3} \] 을 직접 사용하는 쪽이 추천된다. 이는 비율 모수의 표준적인 Fieller 구간 해석과 연결된다.
- 부수 모수 \(\mu_i\) 는 추정식에서 완전히 소거 되었다 — 추정함수의 가장 큰 장점.
- 조건화 선택 \(A_i = Y_{i2} + \theta Y_{i1}\) 는 우연이 아니다: \(Y_{i2} + \theta Y_{i1}\) 와 \(Y_{i1} - \theta Y_{i2}\) 는 직교 선형결합 이며, 조건부 우도에서 자연스럽게 등장하는 “부수 통계량 + 충분 통계량” 분해에 해당.
- 이 예제는 추정함수 이론이 조건부 우도를 구성적으로 유도 하는 도구임을 보여준다.
6 예제 2 — 거석환 (Megalithic Stone Rings)
6.1 문제 설정
평면 위의 점 \(Y_i = (Y_{i1}, Y_{i2})\) 이 중심 \((\omega_1, \omega_2)\), 반지름 \(\rho\) 의 원 둘레 근처에 위치한다. 고고학적으로 이들은 거석(standing stones) 이며, 실제로는 원호의 일부만 남아 있다.
데이터 생성 모형: \[ \begin{aligned} Y_{i1} &= \omega_1 + R_i \cos \epsilon_i, \\ Y_{i2} &= \omega_2 + R_i \sin \epsilon_i, \end{aligned} \tag{9.16} \] 여기서:
- \(R_1, \ldots, R_n\) 은 독립·동일분포, \(E[R_i] = \rho\), \(\operatorname{var}(R_i) = \sigma^2 \rho^2\) (변동계수 \(\sigma\) 일정).
- \(\epsilon_1, \ldots, \epsilon_n\) 은 각도이며, 독립이라고 가정할 필요 없다 (돌들은 원주를 따라 대체로 균등 간격이므로 각도들이 상호 상관됨).
6.2 왜 완전한 확률모형이 아닌가
\(\epsilon_i\) 의 결합분포를 명시하지 않는다. 이는 고고학적 현실성 을 반영한다 — 거석이 정확히 어떤 분포로 배치되었는지는 알 수 없다. 완전 우도 접근이 불가능한 전형적 상황이다.
6.3 기본 추정함수
각 관측점에서 관측 반경 − 모수 반경 을 취한다: \[ g_i = \sqrt{(Y_{i1}-\omega_1)^2 + (Y_{i2}-\omega_2)^2} - \rho = R_i(\omega_1, \omega_2) - \rho. \]
\(E[R_i] = \rho\) 이므로 조건화 \(A = (\epsilon_1, \ldots, \epsilon_n)\) 에서 영평균이다. 각도가 상관되어도 상관없다 — 조건부 로 처리하기 때문.
6.4 도함수 및 계수행렬
편미분 \(\partial g_i/\partial(\omega_1, \omega_2, \rho)\) 을 계산하면:
- \(\partial g_i/\partial\omega_1 = -\cos\epsilon_i\) (즉 \(-(Y_{i1}-\omega_1)/R_i\)).
- \(\partial g_i/\partial\omega_2 = -\sin\epsilon_i\).
- \(\partial g_i/\partial\rho = -1\).
따라서 도함수 벡터는 \((\cos\epsilon_i, \sin\epsilon_i, 1)\) 의 음수이며, 이는 \(g_i\) 와 조건부 독립 이다(가정).
6.5 추정방정식 (9.17)
\(V_i = \operatorname{var}(g_i) = \operatorname{var}(R_i) = \sigma^2 \rho^2\) 을 대입하면: \[ \begin{aligned} \sum_i \frac{Y_{i1}-\omega_1}{\rho^2 R_i}(R_i - \rho) &= \sum_i \frac{\cos\epsilon_i \cdot (R_i - \rho)}{\rho^2} = 0, \\ \sum_i \frac{Y_{i2}-\omega_2}{\rho^2 R_i}(R_i - \rho) &= \sum_i \frac{\sin\epsilon_i \cdot (R_i - \rho)}{\rho^2} = 0, \\ \sum_i \frac{R_i - \rho}{\rho^2} &= 0. \end{aligned} \tag{9.17} \]
이 세 식은 반지름 오차의 제곱합을 최소화하는 최소제곱 방정식 과 동일하다 (각도 오차는 무시).
6.6 왜 무조건부 기댓값을 쓰면 안 되는가 — 식별성의 함정
도함수 \((\cos\epsilon_i, \sin\epsilon_i, 1)\) 에서 \(\epsilon_i\) 가 (균등분포는 아니지만) 동일분포라면 무조건부 기댓값 \[ E[\cos\epsilon] = c_1 \text{ (상수)}, \quad E[\sin\epsilon] = c_2 \text{ (상수)} \] 은 모든 \(i\) 에 대해 같은 값이다. 이것을 계수행렬 \(D\) 에 넣으면 \[ D \propto \begin{pmatrix} c_1 & c_2 & 1 \\ c_1 & c_2 & 1 \\ \vdots \\ c_1 & c_2 & 1 \end{pmatrix} \] 이 되어 rank 가 1 이다. 따라서 세 모수 \((\omega_1, \omega_2, \rho)\) 는 식별 불가능하다.
이에 반해 \(D\) 에 관측된 \(\cos\epsilon_i, \sin\epsilon_i\) 를 그대로 넣으면 (조건부 유지) 각 행이 서로 다르고 rank 가 3 이 되어 식별 가능.
→ 조건화 집합 \(A\) 의 적절한 선택이 식별성 자체를 좌우한다.
6.7 산포(분산) 추정
\(\rho^2 \sigma^2\) 의 불편 추정량: \[ \widehat{\rho}^2 \widetilde{\sigma}^2 = \frac{\sum_i (\widehat{R}_i - \widehat{\rho})^2}{n-3}. \] 자유도는 \(n-3\) (모수 3 개).
6.8 Avebury 거석환 수치 결과
Angell & Barber (1977) 데이터를 네 원호 (A, B, C, W) 에 개별 적합:
| 원호 | \(\widehat{\omega}_1\) | \(\widehat{\omega}_2\) | \(\widehat{\rho}\) | \(\widehat{\rho}^2\widetilde{\sigma}^2\) | 표본수 |
|---|---|---|---|---|---|
| C | 530.8 | 651.0 | 638.8 | 5.60 | 7 |
| W | 1472.0 | 1553.4 | 1840.4 | 3.78 | 16 |
| A | 795.0 | 516.5 | 782.8 | 9.00 | 10 |
| B | 512.7 | 533.1 | 545.4 | 0.72 | 7 |
6.9 동질성 검정 — 세 원호를 하나의 원으로 적합
A, B, C 를 공통 원 으로 적합하면 잔차제곱합 \(= 878.8\) (자유도 21). 개별 적합의 풀링 잔차제곱합: \[ 4 \times 5.60 + 7 \times 9.00 + 4 \times 0.72 = 88.3 \quad (\text{자유도 15}). \] \(F\)-비로 검정하면 \[ F = \frac{(878.8 - 88.3)/(21-15)}{88.3/15} = \frac{790.5/6}{5.89} \approx 22.4 \] 자유도 \((6, 15)\) 에서 극단적으로 유의 (\(p < 0.001\)) — 세 원호가 하나의 원이 아니라는 강력한 증거.
6.10 얕은 원호에서의 모수 상관 — 기하학적 직관
원호 C 의 표준오차와 상관행렬:
\[ \begin{aligned} \text{s.e.}(\widehat{\omega}_1) &= 14.57, \quad \text{s.e.}(\widehat{\omega}_2) = 108.6, \quad \text{s.e.}(\widehat{\rho}) = 108.4, \\ \operatorname{corr} &= \begin{pmatrix} 1 & -0.878 & 1.000 \\ \cdot & 1 & -0.876 \\ \cdot & \cdot & 1 \end{pmatrix}. \end{aligned} \]
\(\widehat{\omega}_2\) 와 \(\widehat{\rho}\) 의 상관계수가 0.99995 로 거의 완전상관이다. 왜 그럴까.
- 원호가 얕다 = 데이터가 원의 작은 부분만 차지한다.
- 이 작은 부분만 관측하면 “중심이 약간 위에 있고 반지름이 조금 큰 원” 과 “중심이 조금 아래에 있고 반지름이 조금 작은 원” 이 거의 같은 모양 을 만든다.
- 따라서 \(\omega_2\) 를 키우면 \(\rho\) 도 같이 키워서 보정하는 식으로, 두 모수가 함께 움직인다.
- 완전한 원이라면 상관은 거의 0 이 된다 — 다양한 곡률 방향의 정보가 모수를 분리해준다.
이 현상은 과적합(overfitting) 의 통계적 표현 이며, 회귀에서 설명변수 간 다중공선성과 같은 본질이다.
99% 신뢰집합은 거의 수직축으로 길쭉한 타원이 되며, 원호 C 의 중심 위치 추정의 불확실성이 수직 방향으로 집중된다.
7 최적성 주장의 정확한 의미 (§9.5 미리보기)
“왜 \(U = D^\top V^{-1} g\) 가 최적 인가” 를 엄밀히 정당화하려면 비교 클래스를 명시해야 한다.
\(h(y;\beta) = H^\top (y - \mu(\beta))\) 꼴의 선형 추정함수 클래스 를 생각하자 (\(H\) 는 \(\beta\) 에 의존하나 \(y\) 에는 의존하지 않음). \(h(y; \widetilde{\beta}) = 0\) 의 해 \(\widetilde{\beta}\) 와 §9.2 준스코어 해 \(\widehat{\beta}\) 를 비교하면: \[ \operatorname{var}(a^\top \widetilde{\beta}) \geq \operatorname{var}(a^\top \widehat{\beta}) \quad \text{(모든 선형결합 } a \text{ 에 대해, 점근적으로)} \] 즉 \(\widehat{\beta}\) 가 Löwner 순서의 의미에서 최소 분산 이다.
이 최적성은 대역적(global) 이 아니다. §9.3.3 의 Voter transition 예제에서 보았듯이, 비선형 추정함수까지 포함하는 더 큰 클래스에서는 최대우도가 준스코어보다 더 효율적일 수 있다. 즉 §9.4 가 말하는 “최적” 은 “선형 추정함수 중 최적” 이라는 조건부 최적성이다.
점근 전개 \[ \widetilde{\beta} - \beta \approx (H^\top D)^{-1} h(y; \beta) \] 에서 \(\operatorname{cov}(\widetilde{\beta}) = (H^\top D)^{-1} H^\top V H (H^\top D)^{-\top}\) 과 \(\operatorname{cov}(\widehat{\beta}) = (D^\top V^{-1} D)^{-1}\) 의 차이가 양반정치(positive semi-definite) 임을 증명하는 것이 §9.5 의 주 내용이다.
8 코드 예시 — 추정함수 이론의 구현
8.1 Python Low-level — Fieller-Creasy 문제
import numpy as np
from scipy.optimize import brentq
# ---- 데이터 시뮬레이션 ----
np.random.seed(0)
n = 50
mu_true = np.random.uniform(1.0, 5.0, size=n) # 부수모수 μ_i
theta_true = 1.5 # 관심 모수
sigma = 0.3
Y1 = mu_true + sigma * np.random.randn(n)
Y2 = mu_true / theta_true + sigma * np.random.randn(n)
# ---- 최적 추정함수 U(θ) ----
def U(theta, Y1, Y2, sigma):
num = (Y2 + theta*Y1) * (Y1 - theta*Y2)
den = sigma**2 * (1 + theta**2)**2
return (num / den).sum()
# ---- 관측 정보 i(θ) ----
def info_theta(theta, Y1, Y2, sigma):
num = (Y2 + theta*Y1)**2
den = sigma**2 * (1 + theta**2)**3
return (num / den).sum()
# ---- 뿌리 찾기: 두 해 중 양수 해 ----
theta_hat = brentq(U, 0.1, 10.0, args=(Y1, Y2, sigma))
se_theta = 1.0 / np.sqrt(info_theta(theta_hat, Y1, Y2, sigma))
print(f"theta_hat = {theta_hat:.4f}")
print(f"true theta = {theta_true}")
print(f"s.e.(theta_hat) = {se_theta:.4f}")
# ---- 비허해(minimum root) 확인 ----
theta_minus = -1/theta_hat
print(f"허해 -1/theta_hat = {theta_minus:.4f} (로그우도의 minimum)")8.2 Python — Avebury 원 적합 (minimum working example)
import numpy as np
from scipy.optimize import least_squares
# 원호 C 데이터 (McCullagh Table 9.3)
arc_C = np.array([
[733.7, 44.0], [659.7, 28.0], [624.2, 19.3],
[588.4, 13.9], [551.6, 12.3], [515.1, 9.5],
[478.0, 16.6]
])
# 추정함수 (radial residual)
def residuals(params, data):
w1, w2, rho = params
r = np.sqrt((data[:,0]-w1)**2 + (data[:,1]-w2)**2)
return r - rho
# 초기값: 데이터 중심 + 평균 반지름
x0 = [arc_C[:,0].mean(), arc_C[:,1].mean(), 500.0]
res = least_squares(residuals, x0, args=(arc_C,))
w1_hat, w2_hat, rho_hat = res.x
n, p = len(arc_C), 3
# 잔차 분산 rho^2 * sigma^2 = RSS / (n-p)
rss = (res.fun**2).sum()
rho2_sigma2 = rss / (n - p)
print(f"omega_1_hat = {w1_hat:.1f} (McCullagh: 530.8)")
print(f"omega_2_hat = {w2_hat:.1f} (McCullagh: 651.0)")
print(f"rho_hat = {rho_hat:.1f} (McCullagh: 638.8)")
print(f"rho^2 * sigma^2 = {rho2_sigma2:.3f} (McCullagh: 5.60)")
# ---- 관측 정보 행렬로부터 표준오차 ----
# J_ir = ∂g_i/∂theta_r = -(cos ε_i, sin ε_i, 1)
R = np.sqrt((arc_C[:,0]-w1_hat)**2 + (arc_C[:,1]-w2_hat)**2)
cos_eps = (arc_C[:,0] - w1_hat) / R
sin_eps = (arc_C[:,1] - w2_hat) / R
J = -np.column_stack([cos_eps, sin_eps, np.ones(n)])
# cov ≈ (J'J)^{-1} * rho^2 sigma^2
cov = np.linalg.inv(J.T @ J) * rho2_sigma2
se = np.sqrt(np.diag(cov))
corr = cov / np.outer(se, se)
print("\n표준오차:", np.round(se, 2))
print("상관행렬:\n", np.round(corr, 4))8.3 R — geepack 을 활용한 일반 추정방정식 스타일
library(geepack)
# 자기회귀 시계열 예제: Y_t = theta * Y_{t-1} + eps_t
set.seed(0)
n <- 200
theta_true <- 0.7
y <- numeric(n); y[1] <- rnorm(1)
for (t in 2:n) y[t] <- theta_true * y[t-1] + rnorm(1)
# 데이터프레임 구성: 각 시점을 한 관측으로 보고 lag 를 공변량으로
df <- data.frame(
id = 1:(n-1),
y = y[-1],
ylag = y[-n]
)
# GEE: 시점 간 독립 가정으로 최적 추정함수 sum(y_{t-1} * (y_t - theta * y_{t-1}))
fit <- geeglm(y ~ ylag - 1, data = df, id = id,
family = gaussian, corstr = "independence")
summary(fit)
# Estimate: theta_hat ≈ 0.7 (true value)
# 추정함수 이론의 일치 확인:
# theta_hat = sum(y_{t-1} * y_t) / sum(y_{t-1}^2)
theta_eq <- sum(y[-n] * y[-1]) / sum(y[-n]^2)
cat("수동 계산 theta_hat =", round(theta_eq, 4), "\n")
cat("geeglm theta_hat =", round(coef(fit), 4), "\n")8.4 Python — AR(1) 최적 추정함수의 직접 구현
import numpy as np
# AR(1) 시뮬레이션
np.random.seed(0)
n, theta_true = 200, 0.7
y = np.zeros(n)
for t in range(1, n):
y[t] = theta_true * y[t-1] + np.random.randn()
# 기본 추정함수 g_t = y_t - theta * y_{t-1}
# V_tt = 1, D_t = y_{t-1}
# U(theta) = sum y_{t-1} (y_t - theta y_{t-1}) = 0
# => theta_hat = sum(y_{t-1} y_t) / sum(y_{t-1}^2)
y_lag, y_now = y[:-1], y[1:]
theta_hat = (y_lag * y_now).sum() / (y_lag**2).sum()
# 점근 분산: (D^T V^{-1} D)^{-1} = 1 / sum(y_{t-1}^2)
var_theta = 1.0 / (y_lag**2).sum()
print(f"theta_hat = {theta_hat:.4f} (true = {theta_true})")
print(f"s.e. = {np.sqrt(var_theta):.4f}")9 §9.4 의 통합 정리
9.1 추정함수 프레임워크의 논리 흐름
기본 추정함수 g_i(Y; θ)
E[g_i | A_i] = 0
↓
계수행렬 D_ir = -E[∂g_i/∂θ_r | A_i]
↓
조건부 공분산 V = cov(g | A)
↓
최적 결합 U = D^T V^{-1} g
↓
U(θ_hat; y) = 0
↓
cov(θ_hat) ≈ (D^T V^{-1} D)^{-1}9.2 §9.2, §9.3, §9.4 의 연결
| 레벨 | \(g\) | \(V\) | \(D\) | 적용 맥락 |
|---|---|---|---|---|
| §9.2 (독립) | \(Y-\mu\) | diag(\(\sigma^2 V(\mu)\)) | \(\partial\mu/\partial\beta\) | 독립 GLM |
| §9.3 (의존) | \(Y-\mu\) | 블록 \(\sigma^2 V(\mu)\) | \(\partial\mu/\partial\beta\) | 종단·군집 |
| §9.4 | 임의 영평균 \(g_i\) | 조건부 cov | \(-E[\partial g/\partial\theta \mid A]\) | 시계열·부수모수·구조모형 |
즉 §9.4 는 §9.2, §9.3 을 포함하는 최상위 추상화 다. \(g = Y - \mu\) 를 선택하면 §9.2 가 나오고, 조건화 구조를 활용하면 §9.3 이 나오며, 임의의 영평균 함수를 허용하면 §9.4 가 된다.
9.3 실무에서 추정함수 이론을 언제 쓰는가
| 상황 | 왜 추정함수가 유리한가 |
|---|---|
| 부수 모수 제거 | 조건화 선택으로 \(\mu_i\) 같은 무한차원 부수모수를 우회 |
| 구조적 기하 모형 | 확률모형 완전 명시가 어려우나 영평균 관계식은 자연스러움 (원 적합, 궤도 추정) |
| 시계열·마팅게일 | 조건부 영평균 성질이 자동 보장, 우도 계산이 복잡해도 추정 가능 |
| 과산포·이분산 | \(V\) 지정만 하면 분포 없이 효율적 추정 (Wedderburn, 1974) |
| GEE 틀에서의 종단 데이터 | Liang-Zeger (1986) 는 §9.4 의 직접적 응용 |