E3 · 오차의 크기와 상관을 반영하는 회귀#

1. 마지막 측정의 잡음이 더 크다면 같은 직선을 쓸까#

E1·E2의 자료 \(t=(0,1,2,3)^T\), \(y=(1,2,2,5)^T\), \(X=[\mathbf1,t]\)를 다시 봅니다. 활동량과 산출량은 같은 기준 단위입니다. 이번에는 측정 장치의 교정 정보로 오차 공분산이

\[ \Omega=\operatorname{diag}(1,1,1,4) \]

라고 알려졌다고 가정합니다. 오차 평균은 0이고 \(X\)는 고정입니다. 마지막 관측의 표준편차는 다른 관측의 두 배입니다. 따라서 오차를 각 표준편차로 나눈 뒤 제곱하면 목적함수는

\[ \sum_{i=1}^3(y_i-a-bt_i)^2+\frac14(y_4-a-3b)^2 \]

입니다. 같은 크기의 마지막 잔차에 더 작은 벌점을 줍니다. 가중 정상방정식은

\[\begin{split} X^T\Omega^{-1}X=\frac14\begin{pmatrix}13&15\\15&29\end{pmatrix},\qquad X^T\Omega^{-1}y=\frac14\begin{pmatrix}25\\39\end{pmatrix}. \end{split}\]

첫 행의 13은 \(4(1+1+1+1/4)\), 15는 \(4(0+1+2+3/4)\)입니다. 두 식 \(13a+15b=25\), \(15a+29b=39\)를 풀면

\[ \widehat\beta_{GLS}=(35/38,33/38)^T. \]

OLS의 \((7/10,6/5)\)보다 기울기가 작습니다. 큰 마지막 산출량이 기울기를 끌어올리는 정도가 줄었습니다. 이 결과는 알려진 공분산이 맞다는 가정에 의존합니다. 원하는 직선을 얻도록 관측의 가중치를 임의로 정한 것이 아닙니다.

동일한 네 관측에 대한 OLS와 마지막 관측의 분산을 4로 둔 GLS 직선 비교

그림 137 오차막대는 모형에서 주어진 오차 표준편차 1,1,1,2이며 추정된 평균의 신뢰구간이 아니다. 마지막 관측의 낮은 정밀도를 반영한 GLS는 더 작은 기울기를 갖는다.#

2. 역공분산 내적과 백색화는 같은 계산이다#

\(\Omega\succ0\)일 때 이 장의 내적은

\[ \langle a,b\rangle_{\Omega^{-1}}=a^T\Omega^{-1}b \]

입니다. 공분산 자체와 내적의 가중행렬인 역공분산을 구별합니다. GLS는 이 내적에서 \(y\)\(\operatorname{col}X\)에 사영합니다.

\(\Omega=LL^T\)인 Cholesky 인수를 사용하면 \(\widetilde y=L^{-1}y\), \(\widetilde X=L^{-1}X\), \(\widetilde u=L^{-1}u\)이고 \(\operatorname{Var}\widetilde u=I\). 목적함수는 \(\|\widetilde y-\widetilde X\beta\|^2\)가 됩니다. 앞 예에서는 \(L=\operatorname{diag}(1,1,1,2)\)라 마지막 행 전체, 곧 반응과 설명변수를 모두 2로 나눕니다.

\[ P_\Omega=X(X^T\Omega^{-1}X)^{-1}X^T\Omega^{-1},\quad P_\Omega^2=P_\Omega,\quad P_\Omega^T\Omega^{-1}=\Omega^{-1}P_\Omega. \]

일반적으로 \(P_\Omega^T\ne P_\Omega\)입니다. 유클리드 대칭성이 아니라 선택한 내적에서의 자기수반성이 맞는 조건입니다. 가중 잔차에는 \(X^T\Omega^{-1}(y-X\widehat\beta)=0\)이 성립합니다.

상관도 같은 방식으로 처리합니다. 두 측정의 공분산이 \(\sigma^2\begin{pmatrix}1&\rho\\\rho&1\end{pmatrix}\), \(|\rho|<1\)이면

\[\begin{split} L=\sigma\begin{pmatrix}1&0\\\rho&\sqrt{1-\rho^2}\end{pmatrix},\qquad L^{-1}y=\begin{pmatrix}y_1/\sigma\\(y_2-\rho y_1)/(\sigma\sqrt{1-\rho^2})\end{pmatrix}. \end{split}\]

둘째 좌표는 첫 관측과 함께 움직이는 부분을 뺀 뒤 남은 표준편차로 나눕니다. \(\rho\to1\)에서는 차이 방향의 분산이 0에 가까워져 나눗셈이 민감해집니다. \(\rho=1\)에서는 이 Cholesky 백색화가 정의되지 않고 지지공간을 따로 처리해야 합니다.

3. 계산에서는 삼각풀이 다음 QR을 사용한다#

공식을 설명할 때 역행렬을 쓰더라도 구현에서 전체 역행렬을 만들 필요는 없습니다. \(L\widetilde X=X\), \(L\widetilde y=y\)를 삼각풀이로 해결한 뒤 QR 또는 SVD 최소제곱을 적용합니다. 이 방식은 \(\widetilde X^T\widetilde X\)를 만들 때 생기는 조건수 제곱을 피합니다. 정확산술에서는 \(\kappa_2(\widetilde X^T\widetilde X)=\kappa_2(\widetilde X)^2\)입니다. 이것을 \(\kappa(\Omega)^2\)라는 보편 전진오차 법칙으로 바꾸지는 않습니다.

import numpy as np
import scipy.linalg as la
t=np.arange(4.);X=np.column_stack([np.ones(4),t]);y=np.array([1.,2.,2.,5.])
Omega=np.diag([1.,1.,1.,4.])
def gls_whiten(X,y,Omega):
    L=la.cholesky(Omega,lower=True)
    Xw=la.solve_triangular(L,X,lower=True)
    yw=la.solve_triangular(L,y,lower=True)
    b=la.lstsq(Xw,yw)[0]
    return b,Xw,yw
b,Xw,yw=gls_whiten(X,y,Omega)
assert np.allclose(b,[35/38,33/38])
assert np.allclose(Xw.T@(yw-Xw@b),0.,atol=1e-12)
Sinv=la.solve(X.T@X,np.eye(2),assume_a='pos')
C=Sinv@X.T
V_ols=C@Omega@C.T
V_gls=la.solve(Xw.T@Xw,np.eye(2),assume_a='pos')
expected=27/1900*np.array([[4.,-6.],[-6.,9.]])
assert np.allclose(V_ols-V_gls,expected)
print('GLS 계수:',b,'분산 차이의 고윳값:',la.eigvalsh(V_ols-V_gls))
GLS 계수: [0.92105263 0.86842105] 분산 차이의 고윳값: [-1.04083409e-16  1.84736842e-01]

고유분해 \(\Omega=U\Lambda U^T\)\(\Lambda^{-1/2}U^T\)를 써도 백색화할 수 있습니다. Cholesky는 삼각구조를 제공하고 고유분해는 작은 고윳값과 지지를 드러냅니다. 어느 좌표를 쓰든 정확한 가역 백색화면 같은 GLS가 됩니다. 입력 공분산의 추정오차가 크면 두 수치 방법 모두 그 모형 오차를 고칠 수 없습니다.

4. 알려진 공분산에서 GLS가 최소분산인 이유#

\(C_*=(X^T\Omega^{-1}X)^{-1}X^T\Omega^{-1}\)라 쓰면 GLS는 \(C_*y\)입니다. 다른 선형불편추정량 \(Cy\)\(CX=I\)를 만족해야 합니다. \(D=C-C_*\)라 두면 \(DX=0\)이고

\[ \operatorname{Var}(Cy)-\operatorname{Var}(C_*y)=D\Omega D^T\succeq0. \]

교차항이 0이 되는 과정을 마지막 절에서 확인합니다. 이는 각 계수의 분산뿐 아니라 모든 선형 결합의 분산에 대한 비교입니다. 정규성은 필요하지 않으며 평균·공분산과 선형불편성 조건만 씁니다.

처음 예의 분산행렬은

\[\begin{split} V_{OLS}=\frac1{100}\begin{pmatrix}82&-48\\-48&47\end{pmatrix},\qquad V_{GLS}=\frac1{38}\begin{pmatrix}29&-15\\-15&13\end{pmatrix}, \end{split}\]
\[\begin{split} V_{OLS}-V_{GLS}=\frac{27}{1900} \begin{pmatrix}4&-6\\-6&9\end{pmatrix} =\frac{27}{1900}(2,-3)^T(2,-3). \end{split}\]

따라서 차이는 PSD입니다. 관측한 자료에서 GLS의 잔차제곱합이나 실제 계수오차가 항상 더 작다는 주장은 아닙니다. 비교한 것은 같은 참모형 아래 반복표본 분산입니다.

OLS와 GLS가 모든 \(y\)에 대해 같을 조건은 \(\operatorname{col}X\)\(\Omega\)의 불변부분공간인 것입니다. 대칭성 때문에 이는 \(\Omega P_X=P_X\Omega\)와 동치입니다. 예컨대 절편만 있는 두 관측의 등분산·등상관 공분산은 상수 방향을 보존하므로 GLS도 표본평균입니다. 한 특정 \(y\)에서 우연히 계수가 같다는 것만으로 이 조건을 추론할 수는 없습니다.

5. 공분산을 추정하면 효율성 보장은 별도 문제다#

실제로는 \(\widehat\Omega\)를 넣은 FGLS를 사용합니다. 가중치가 자료에 의존하면 추정량은 더 이상 고정 \(C\)를 곱한 선형추정량이 아닐 수 있고, 위 Aitken 증명을 그대로 적용할 수 없습니다. 잘못 추정한 정밀도는 분산을 늘릴 수 있습니다.

반례를 완전히 계산해 봅시다. 독립인 두 관측 \(Y_i=\mu+u_i\)에 실제 \(\operatorname{Var}u_i=1\)이면 평균의 분산은 \(1/2\)입니다. 독립 예비조사에서 두 공분산 후보 \(\operatorname{diag}(1,99)\)\(\operatorname{diag}(99,1)\) 중 하나를 선택하는 부정확한 절차를 가정합니다. 어느 후보에서도 FGLS 가중치는 \((99/100,1/100)\) 또는 그 반대입니다. 평균에는 불편이지만 분산은

\[ (99/100)^2+(1/100)^2=4901/5000>1/2. \]

이 예는 좋은 일치적 공분산 추정의 점근 효율성을 부정하지 않습니다. 유한표본에서 추정 가중치가 자동으로 개선을 보장하지 않는다는 반례입니다. 특정 AR(1) 실험에서 반드시 역전이 나와야 한다고 시드나 결과를 조정하지 않습니다.

6. 강건 공분산은 계수 추정과 다른 단계다#

고정 설계 OLS의 참분산은 \(S^{-1}X^T\Omega XS^{-1}\), \(S=X^TX\)입니다. GLS로 계수를 다시 추정하는 대신 이 분산만 추정할 수도 있습니다. 이분산, 클러스터, 시계열 상관은 가운데 행렬을 어떻게 추정하는지의 차이로 나타납니다.

서로 독립이라고 가정한 클러스터 \(g=1,\ldots,G\)에 점수합 \(s_g=X_g^Te_g\)를 만들면 통상적인 보정 전 클러스터 공분산은

\[ \widehat V_{cl}=S^{-1}\left(\sum_gs_gs_g^T\right)S^{-1}. \]

OLS 정규방정식 때문에 \(\sum_gs_g=X^Te=0\)입니다. 따라서 가운데 행렬의 계수는 \(\min(k,G-1)\) 이하입니다. 이 결과는 전체 표본 수가 크더라도 클러스터 수가 적으면 공분산의 독립 방향이 제한됨을 보여 줍니다.

처음 자료를 앞 두 관측과 뒤 두 관측의 두 집단으로 나누면 \(s_1=(2/5,1/10)^T\), \(s_2=-s_1\)입니다. \(S^{-1}s_1=(1/4,-1/10)^T\)이므로

\[\begin{split} \widehat V_{cl}=\begin{pmatrix}1/8&-1/20\\-1/20&1/50\end{pmatrix},\qquad \operatorname{rank}\widehat V_{cl}=1. \end{split}\]

두 독립 제약을 이 행렬의 보통 역행렬로 동시에 표준화할 수 없습니다. 유사역을 넣고 자유도만 줄이면 원래 두 제약의 유효한 검정이 자동으로 되는 것도 아닙니다. E2의 알려진 특이 정규분포와 달리 추정 과정이 만든 계수 부족일 수 있기 때문입니다. wild cluster 재표집·무작위화 검정은 별도 설계와 타당성 조건이 필요하며 계수 부족을 보편적으로 해결하지 않습니다.

7. 시계열 가중치가 공분산을 음수로 만들지 않으려면#

시점별 점수 \(g_t=x_te_t\in\mathbb R^k\)를 행으로 쌓은 행렬을 \(G_0\)라 합시다. 시차 \(j\) 가중 \(w_j\)로 만든 HAC 가운데 행렬은

\[ \widehat B=G_0^TKG_0,\qquad K_{ts}=w_{|t-s|}. \]

Toeplitz 구조는 시점 가중행렬 \(K\)에 있습니다. 결과인 변수별 공분산 \(\widehat B\) 전체가 Toeplitz일 필요는 없습니다. Bartlett 가중 \(w_j=\max(1-j/(m+1),0)\)\(K\succeq0\)를 보장합니다.

반대로 일정 시차까지 모두 1을 주고 자르면 항상 PSD가 되지 않습니다. 세 시점, 대역폭 1이면

\[\begin{split} K_{flat}=\begin{pmatrix}1&1&0\\1&1&1\\0&1&1\end{pmatrix},\qquad v=(1,-2,1)^T,\quad v^TK_{flat}v=-2. \end{split}\]

Bartlett은 인접 가중치가 \(1/2\)여서 같은 벡터에 \(v^TK_{Bartlett}v=2\)입니다. 잔차합이 0인 \(v\)를 절편 회귀의 점수로 해석할 수도 있으므로 단순히 불가능한 점수를 만든 반례가 아닙니다.

b_ols=la.lstsq(X,y)[0];e=y-X@b_ols
scores=np.column_stack([X[:2].T@e[:2],X[2:].T@e[2:]])
Vcl=Sinv@scores@scores.T@Sinv
assert np.allclose(scores.sum(axis=1),0.,atol=1e-12)
assert np.linalg.matrix_rank(Vcl,tol=1e-12)==1
assert np.allclose(Vcl,[[1/8,-1/20],[-1/20,1/50]])
lag=np.abs(np.arange(3)[:,None]-np.arange(3)[None,:])
Kflat=(lag<=1).astype(float);Kbart=np.maximum(1-lag/2,0.)
v=np.array([1.,-2.,1.])
assert np.isclose(v@Kflat@v,-2.) and np.isclose(v@Kbart@v,2.)
assert la.eigvalsh(Kbart).min()>0
print('클러스터 계수와 HAC 가중치의 부호 확인')
클러스터 계수와 HAC 가중치의 부호 확인
세 시점의 Bartlett 및 균등 절단 가중행렬과 두 행렬의 고윳값 비교

그림 138 두 가중행렬은 같은 색 범위로 표시한다. 오른쪽 고윳값에서 균등 절단은 \(1-\sqrt2<0\)인 방향을 갖고 Bartlett은 모두 양수이다. 이는 특정한 세 시점 예이며 Bartlett의 일반 PSD 증명은 마지막 절에 있다.#

PSD 추정량에 \(\lambda I\)를 더하거나 항등행렬 쪽으로 축소하면 작은 고윳값을 올릴 수 있습니다. 그러나 수치 가역성을 얻는 일과 원래 검정의 분포·자유도를 정당화하는 일은 다릅니다. 대역폭, 클러스터 독립성, 시계열 혼합과 적률, 축소 편향의 통계 이론은 사용 목적에 맞춰 추가해야 합니다.

8. 반복 가중최소제곱으로 로지스틱 회귀를 푼다#

이번에는 산출량 대신 이진 성공 여부 \(y_i\in\{0,1\}\)를 관측합니다. 독립 Bernoulli 모형에서 성공확률 \(p_i=1/(1+e^{-x_i^T\beta})\)라 하면 로그우도와 점수는

\[ \ell(\beta)=\sum_i[y_i\eta_i-\log(1+e^{\eta_i})],\quad s=X^T(y-p),\quad -H=X^TWX, \quad W_{ii}=p_i(1-p_i). \]

Newton 보정은 \(X^TWX\,\Delta=X^T(y-p)\)입니다. 현재 \(0<p_i<1\)일 때 작업반응 \(z=\eta+W^{-1}(y-p)\)를 만들면 같은 식을 \(X^TWX\beta_{new}=X^TWz\)로 쓸 수 있습니다. 현재 계수로 가중치를 정하고 가중최소제곱을 다시 풀기 때문에 IRLS라 부릅니다.

두 활동 수준 0과 1에서 각각 네 번 관측하여 성공 수가 1과 3이라고 합시다. 절편과 활동 효과가 있으면 각 수준의 확률을 따로 맞출 수 있어 최적 확률은 \(1/4,3/4\). 따라서 \(\beta_0=-\log3\), \(\beta_1=2\log3\)입니다. 초기 \(\beta=0\)에서는 \(p_i=1/2\)이고

\[\begin{split} X^TWX=\begin{pmatrix}2&1\\1&1\end{pmatrix},\quad s=(0,1)^T,\quad \Delta=(-1,2)^T. \end{split}\]

첫 단계가 두 로그오즈 \(-1,1\)을 만들고 다음 단계들이 \(-\log3,\log3\)로 보정합니다.

from scipy.special import expit
def logistic_fit(X,y,limit=80,tol=1e-11):
    beta=np.zeros(X.shape[1]);history=[]
    loss=lambda b:np.sum(np.logaddexp(0,X@b)-y*(X@b))
    for iteration in range(limit):
        eta=X@beta;p=expit(eta);w=p*(1-p);score=X.T@(y-p)
        info=X.T@(w[:,None]*X)
        history.append((loss(beta),la.norm(score),la.eigvalsh(info).min()))
        if la.norm(score)<tol:return beta,np.array(history)
        step=la.solve(info,score,assume_a='pos')
        alpha=1.
        while loss(beta+alpha*step)>loss(beta)-1e-4*alpha*(score@step):
            alpha*=.5
            if alpha<2**-30:raise RuntimeError('보폭 축소 실패')
        beta+=alpha*step
    raise RuntimeError('정지 기준에 도달하지 못함')
xl=np.repeat([0.,1.],4);Xl=np.column_stack([np.ones(8),xl])
yl=np.array([0.,0.,0.,1.,0.,1.,1.,1.])
bl,history=logistic_fit(Xl,yl)
assert np.allclose(bl,[-np.log(3),2*np.log(3)],atol=1e-9)
assert np.all(np.diff(history[:,0])<=1e-12)
print('IRLS/Newton과 닫힌 해:',bl)
IRLS/Newton과 닫힌 해: [-1.09861229  2.19722458]

일반 지수족 GLM에서는 Fisher scoring이 IRLS 형태를 갖습니다. 정준연결인 로지스틱에서는 관측 음의 헤시안과 기대정보가 같아 Newton까지 일치합니다. 비정준연결에서는 잔차가 포함된 헤시안 항 때문에 두 방법이 다를 수 있습니다. Gauss–Newton은 비선형 잔차를 선형화한 별도의 방법이며, 같은 가중최소제곱 계산을 공유한다는 이유만으로 모든 모형에서 세 방법이 동일하다고 할 수 없습니다.

완전분리도 계산해야 합니다. 두 관측 \(x=-1,1\), 반응 \(y=0,1\), 절편 0과 기울기 \(b>0\)에서는 로그우도가 \(b\to\infty\)로 갈수록 증가하여 유한 최대점이 없습니다. 정보행렬은 \(2p(1-p)I_2\)이고 \(p=\operatorname{logit}^{-1}(b)\). 두 고윳값이 함께 0으로 가므로 조건수는 계속 1입니다. 조건수만 점검하면 정보의 절대 크기 소멸을 놓칩니다.

유한 해를 가진 로지스틱 자료의 우도 개선과 완전분리에서 정보 고윳값이 사라지지만 조건수가 1인 예

그림 139 왼쪽은 코드의 음의 로그우도 감소이다. 오른쪽은 별도의 완전분리 두 관측 모형에서 기울기를 증가시킨 결과이며 고윳값 축은 로그이다. 서로 다른 두 자료의 현상을 구별하여 표시했다.#

9. 연습과 전체 풀이#

1. 처음 GLS 예의 마지막 적합값을 구하세요.

풀이. \(35/38+3(33/38)=134/38=67/19\). OLS의 \(43/10\)보다 작습니다. 마지막 관측 자체를 없앤 것이 아니라 그 정밀도를 낮게 반영한 결과입니다.

2. 공분산이 \(\sigma^2I\)이면 GLS와 OLS가 왜 같은가요?

풀이. 목적함수가 OLS 목적함수의 양의 상수 \(1/\sigma^2\)배이므로 최소점이 같습니다. Aitken의 분산 비교는 이 경우 Gauss–Markov 정리가 됩니다.

3. 두 관측의 상관 \(\rho=1/2\), \(\sigma=1\)에서 백색화 둘째 좌표를 구하세요.

풀이. \((y_2-y_1/2)/(\sqrt3/2)=(2y_2-y_1)/\sqrt3\). 원래 공분산의 공유 성분을 제거한 뒤 분산 1로 만든 좌표입니다.

4. 클러스터가 12개일 때 12개의 독립 선형제약에 보통 Wald 역행렬을 쓸 수 있나요?

풀이. 통상적인 OLS 클러스터 공분산 계수는 11 이하라 \(R\widehat VR^T\)가 12차원에서 가역일 수 없습니다. 제약의 정당한 축소나 별도의 검정 설계가 필요합니다. 유사역·재표집을 넣기만 하면 해결된다는 결론은 없습니다.

5. 균등 절단 가중의 음수 이차형식을 직접 계산하세요.

풀이. \(v=(1,-2,1)\)에 대각 기여는 \(1+4+1=6\), 인접 교차 기여는 \(2(1\cdot(-2)+(-2)\cdot1)=-8\)입니다. 합은 \(-2\). Bartlett에서는 교차 기여가 절반인 \(-4\)여서 합 2입니다.

6. 로지스틱 예의 첫 Newton 계수를 구하세요.

풀이. \(2\Delta_0+\Delta_1=0\), \(\Delta_0+\Delta_1=1\)을 풀면 \((-1,2)\). 초기 계수가 0이므로 그 자체가 새 계수입니다. 두 수준의 예측확률은 \(\operatorname{logit}^{-1}(-1)\)\(\operatorname{logit}^{-1}(1)\)입니다.

7. 완전분리 예에서 조건수가 작으면 표준오차도 안정적인가요?

풀이. 아닙니다. 정보의 두 고윳값이 모두 \(2p(1-p)\to0\)이라 역행렬 대각은 무한대로 갑니다. 상대적인 방향 불균형은 없지만 정보의 절대량이 사라지고 유한 MLE도 없습니다.

8. \(\widehat\Omega+\lambda I\)를 사용하면 원래 모형의 정확 F 분포를 되찾나요?

풀이. 양의 이동은 가역성과 수치 안정성을 개선할 수 있지만 공분산의 진실성·정규성·독립성 가정을 만들지는 않습니다. 가중치 추정과 정칙화가 통계량 분포에 미치는 영향은 따로 정당화해야 합니다.

10. 지금까지의 내용을 수학의 언어로 정리해 봅시다#

GLS의 대수·유한표본 분산 비교와 강건 추론의 점근 타당성은 다른 명제입니다. 이 절에서는 사용한 대수를 완전히 증명하고 일반 HAC·클러스터 극한 이론은 전제하지 않습니다.

정리 1. 백색화 사영과 Aitken 하한#

\(X\)가 완전 열계수이고 \(E[u]=0\), \(\operatorname{Var}u=\Omega\succ0\)이면 GLS는 백색화 후 OLS이며 모든 선형불편추정량 중 공분산이 Loewner 의미에서 최소입니다.

증명. \(\Omega=LL^T\)에서 \(r^T\Omega^{-1}r=\|L^{-1}r\|^2\)이므로 두 최소화 문제가 같습니다. \(\widetilde X\)도 완전 열계수여서 해는 유일하고 정규방정식으로 GLS 공식을 얻습니다. \(W=\Omega^{-1}\), \(S_W=X^TWX\)라 쓰면 \(P_\Omega=XS_W^{-1}X^TW\). 직접 곱하여 멱등성을, 전치하여 \(P_\Omega^TW=WP_\Omega\)를 얻습니다. 상은 \(\operatorname{col}X\)이므로 이 내적의 사영입니다.

\(C_*=S_W^{-1}X^TW\)이면 \(C_*X=I\), \(\operatorname{Var}(C_*y)=S_W^{-1}\). 다른 \(CX=I\)\(D=C-C_*\)를 두면 \(DX=0\). 교차항은 \(C_*\Omega D^T=S_W^{-1}X^TD^T=0\)이고 전치 교차항도 0입니다. 따라서 \((C_*+D)\Omega(C_*+D)^T-C_*\Omega C_*^T=D\Omega D^T\succeq0\). 정규성은 어느 단계에서도 사용하지 않았습니다. ∎

정리 2. OLS와 GLS의 일치 조건#

두 계수 추정량이 모든 \(y\)에서 같을 필요충분조건은 \(\Omega\operatorname{col}X\subseteq\operatorname{col}X\)이며, 이는 \(\Omega P_X=P_X\Omega\)와 동치입니다.

증명. \(S=\operatorname{col}X\)가 대칭 \(\Omega\)의 불변공간이면 \(v\in S^\perp\), \(s\in S\)\(\langle\Omega v,s\rangle=\langle v,\Omega s\rangle=0\)이어서 \(S^\perp\)도 불변입니다. 가역성으로 두 공간 모두 \(\Omega^{-1}\)에도 불변입니다. OLS 잔차 \(r\in S^\perp\)\(X^T\Omega^{-1}r=0\)이므로 GLS 정규방정식도 만족하여 해가 같습니다.

역으로 모든 \(y\)의 해가 같다고 하자. 특히 \(y\in S^\perp\)이면 OLS 계수는 0이므로 GLS 정규방정식에서 \(X^T\Omega^{-1}y=0\). 따라서 \(S^\perp\)\(\Omega^{-1}\)에 불변이고 대칭성으로 \(S\)도 불변입니다. 두 제한사상이 가역이므로 \(\Omega\)에도 불변입니다. 마지막으로 두 직교공간을 각각 보존하는 것과 \(P_X\)와 가환하는 것은 각 성분에 작용시켜 양방향으로 확인합니다. ∎

정리 3. 클러스터 계수 상한#

OLS 잔차로 만든 통상적인 클러스터 가운데 행렬은 계수가 \(\min(k,G-1)\) 이하입니다.

증명. \(B=[s_1,\ldots,s_G]\)라 두면 \(\widehat B_{cl}=BB^T\)이고 \(\operatorname{rank}(BB^T)=\operatorname{rank}B\). 정규방정식에서 \(B\mathbf1=0\)이므로 \(G\)개 열은 종속이고 계수가 \(G-1\) 이하입니다. 행 수가 \(k\)라는 상한도 있습니다. 가역 \(S^{-1}\)로 양쪽을 곱해도 계수는 그대로입니다. 전체에 같은 양의 소표본 보정 상수를 곱해도 결론은 유지되지만, 클러스터별 잔차 보정을 바꾼 모든 변형에 이 증명을 자동 적용하지 않습니다. ∎

정리 4. Bartlett HAC의 PSD성#

정수 \(m\ge0\)에서 \(K_{ts}=\max(1-|t-s|/(m+1),0)\)는 모든 표본 길이에 PSD입니다.

증명. 정수 시작점 \(j\)에 길이 \(m+1\) 창의 지시벡터 \(b_j\)\((b_j)_t=1\{j\le t\le j+m\}\)로 정의합니다. 두 시점 \(t,s\)를 동시에 포함하는 창의 개수는 \(\max(m+1-|t-s|,0)\)입니다. 따라서

\[ K=\frac1{m+1}\sum_{j=1-m}^{n}b_jb_j^T. \]

임의 \(v\)\(v^TKv=(m+1)^{-1}\sum_j(b_j^Tv)^2\ge0\). 그러므로 임의 점수행렬에 \(G_0^TKG_0\)도 PSD입니다. 이 증명은 창의 겹침만 사용하며 Fourier 이론이나 무한 시계열 극한을 요구하지 않습니다. ∎

정리 5. GLM의 Fisher scoring과 IRLS#

독립 지수족 관측의 평균 \(\mu_i\), 분산 \(V_i>0\), 선형예측자 \(\eta_i=x_i^T\beta\), 미분 \(d_i=d\mu_i/d\eta_i\ne0\)라 합시다. 알려진 분산모수와 정칙 미분 조건 아래 점수와 기대정보는

\[ s=X^T\operatorname{diag}(d_i/V_i)(y-\mu),\qquad I=X^T\operatorname{diag}(d_i^2/V_i)X. \]

따라서 scoring은 \(W_{ii}=d_i^2/V_i\), \(z_i=\eta_i+(y_i-\mu_i)/d_i\)인 가중최소제곱 갱신입니다. 정준연결에서는 Newton과도 같습니다.

증명. 지수족 로그밀도를 \([y_i\theta_i-b(\theta_i)]/\phi+c(y_i,\phi)\)로 쓰면 \(\mu_i=b'(\theta_i)\), \(V_i=\phi b''(\theta_i)\). 연쇄법칙에서 \(d\theta_i/d\eta_i=d_i/b''(\theta_i)\)이므로 \(\partial\ell_i/\partial\eta_i=(y_i-\mu_i)d_i/V_i\). 독립성으로 점수 공분산을 합하면 기대정보 식을 얻습니다. 또는 한 번 더 미분할 때 \((y_i-\mu_i)\)가 곱해진 항의 기대가 0이어서 같은 식입니다.

scoring 식 \(I\Delta=s\)에서 \(\beta_{new}=\beta+\Delta\)라 쓰면 \(X^TWX\beta_{new}=X^TWX\beta+X^TW(z-X\beta)=X^TWz\). \(X\)가 완전 열계수이면 \(W\succ0\)라 유일한 가중최소제곱 해입니다.

정준연결 \(\theta_i=\eta_i\)에서는 점수가 \(X^T(y-\mu)/\phi\)이고 음의 헤시안이 \(X^T\operatorname{diag}(b''(\eta_i)/\phi)X\)라 반응에 의존하지 않아 기대정보와 일치합니다. 비정준연결에서는 생략했던 잔차 항이 관측 헤시안에 남으므로 일반적 동일성은 없습니다. 비선형 최소제곱의 Gauss–Newton은 잔차의 일차 Taylor 근사에서 별도로 유도되며 O4의 잔차 헤시안 항을 생략합니다. ∎

계산 선택

필요한 조건과 뜻

GLS로 계수 자체를 바꿈

알려진 양의 정부호 공분산

OLS 계수에 강건 공분산을 붙임

오차 구조에 맞는 가운데 행렬 추정

클러스터 점수를 외적합함

정규방정식이 주는 계수 상한

Bartlett 시차 가중을 사용함

창의 외적합으로 PSD 보존

IRLS를 반복함

현재 평균·분산에서 만든 작업 최소제곱

백색화는 공분산에 맞는 내적을 표준 최소제곱 계산으로 옮깁니다. 다음 장에서는 이 곡률을 우도의 정보행렬로 읽고, 관심 모수와 나머지 모수를 분리합니다. E4로 이어 읽기.