표본분포 (일반 최소제곱)¶

앞 절의 단순회귀 결과는 행렬대수를 이용해 일반적인 다중회귀 모형으로 확장된다. 일반적인 경우에 최소제곱추정량 \(\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}\)은 다변량 정규벡터이고, 잔차제곱합을 \(\sigma^2\)으로 나눈 것은 카이제곱을 따르며, 이 둘은 독립이다. 정규성, 카이제곱, 독립성이라는 이 세 사실이 결합되어 개별 계수에 대한 t-통계량과 계수 집합을 검정하는 F-통계량이 나온다. 증명은 이 장에서 앞서 전개한 사영행렬과 이차형식 이론에 기댄다.

일반 선형모형¶

행렬 형태의 모형은

이며, 여기서

• \(\mathbf{y} \in \mathbb{R}^n\)은 반응벡터다.
• \(\mathbf{X} \in \mathbb{R}^{n \times p}\)는 \(\operatorname{rank}(\mathbf{X}) = p\)인(완전 열계수, \(p \leq n\)) 고정된 계획행렬이다.
• \(\boldsymbol{\beta} \in \mathbb{R}^p\)는 미지의 모수벡터다.
• \(\boldsymbol{\varepsilon} \sim N(\mathbf{0}, \sigma^2\mathbf{I}_n)\): 오차는 평균 0, 공통 분산 \(\sigma^2\)을 갖는 독립 동일분포 정규확률변수다.

이 가정들로부터 \(\mathbf{y} \sim N(\mathbf{X}\boldsymbol{\beta}, \sigma^2\mathbf{I}_n)\)이다.

최소제곱추정량¶

최소제곱추정량은 \(\lVert\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\rVert^2\)을 최소화한다.

이것은 \(\mathbf{y}\)의 선형함수다. 행렬 \((\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\)가 \(\mathbf{y}\)를 \(\hat{\boldsymbol{\beta}}\)으로 보낸다.

최소제곱추정량의 표본분포¶

그림으로 보는 공분산 구조¶

단순회귀에서 기울기의 분산은 \(\sigma^2/S_{xx}\)라는 숫자 하나였다. 다중회귀에서 그 자리에 들어오는 \(\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}\)은 행렬이고, 대각 성분 바깥에 있는 비대각 성분은 단순회귀에 없던 이야기를 한다. 계수들의 추정이 서로 얽히며, 그 얽힘은 자료의 잡음이 아니라 계획행렬의 기하에서 나온다. 그림은 그 대응을 가능한 한 깨끗한 설정에서 보인다. \(n = 40\), 절편과 설명변수 둘. 두 설명변수 열을 절편 열에 직교시키고 길이를 같게 맞춘 뒤 둘 사이의 상관 \(\rho\)만 조절하면, \(\mathbf{X}^\top\mathbf{X}\)의 해당 \(2 \times 2\) 블록은 대각 성분이 \(n\)이고 비대각 성분이 \(n\rho\)인 행렬이 되어

를 얻는다. \(2 \times 2\) 역행렬이 비대각 성분의 부호를 뒤집는 것이 전부다.

왼쪽 그림은 \(\rho = 0\)인 경우다. \(\mathbf{X}^\top\mathbf{X}\)가 대각행렬이면 역행렬도 대각행렬이고, 표본분포의 95% 등고선은 원이 된다. \(60{,}000\)번의 모의실험에서 \(\operatorname{Var}(\hat{\beta}_1) = 0.02522\)로 이론값 \(\sigma^2/n = 0.025\)와 맞고, 두 계수 추정량의 상관은 \(+0.0052\)로 이론값 0과 맞는다. 한 계수를 이번 표본에서 과대추정했다는 사실이 다른 계수에 대해 아무것도 말해 주지 않는 상태다.

가운데 그림은 \(\rho = 0.9\)다. 타원이 기울면서 동시에 커진다. 측정된 상관은 \(-0.9007\)로 이론값 \(-\rho = -0.9\)와 맞는다. 부호가 음인 것을 읽는 법은 이렇다. 두 열이 거의 같은 방향을 가리키므로 \(\beta_1 x_1 + \beta_2 x_2\)는 사실상 \((\beta_1 + \beta_2)x_1\)에 가깝고, 자료는 합 \(\beta_1 + \beta_2\)는 잘 집어내지만 \(\beta_1\)을 올리고 \(\beta_2\)를 그만큼 내린 조합과 그 반대를 구별하지 못한다. 타원의 긴 축이 바로 그 구별되지 않는 방향이고, 짧은 축은 여전히 잘 결정되는 방향이다. \(\operatorname{Var}(\hat{\beta}_1)\)은 \(0.13294\)로 이론값 \(0.13158\)과 맞아, 왼쪽보다 다섯 배 넘게 부풀었다.

오른쪽 그림은 그 부풂을 \(\rho\)의 함수로 그린 것이다. 분산팽창은 정확히 \(1/(1 - \rho^2)\)이며, \(\rho = 0.9\)에서 모의값 \(5.2552\) 대 이론값 \(5.2632\), \(\rho = 0.95\)에서 \(10.1603\) 대 \(10.2564\)로 맞는다. \(\rho \to 1\)에서 값이 발산하는 것은 두 열이 완전히 겹치면 \(\mathbf{X}^\top\mathbf{X}\)가 특이행렬이 되어 \((\mathbf{X}^\top\mathbf{X})^{-1}\) 자체가 없어지기 때문이다. 모형의 가정에 \(\operatorname{rank}(\mathbf{X}) = p\)가 들어 있는 이유가 여기 있다. 그림이 덧붙이는 것은 그 가정이 있다와 없다로 갈리는 문제가 아니라 정도의 문제라는 사실이다. 동시에, 두 타원의 중심이 모두 참값 \((2, -1)\)에 그대로 있다는 점도 눈여겨볼 만하다. 다중공선성은 추정값을 치우치게 하지 않는다. \(\hat{\boldsymbol{\beta}}\)은 여전히 불편이며, 다만 분산이 부풀 뿐이다.

적합값과 잔차¶

적합값과 잔차는

이며, 여기서 \(\mathbf{H} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\)는 모자 행렬이고 \(\mathbf{M} = \mathbf{I} - \mathbf{H}\)는 잔차생성행렬이다.

이들의 분포는 다음과 같다.

\(\mathbf{M}\mathbf{X} = \mathbf{0}\)이므로 \(\mathbf{e}\)의 평균은 \(\boldsymbol{\beta}\)와 무관하게 0이지만, 공분산행렬 \(\sigma^2\mathbf{M}\)은 특이행렬이므로(계수 \(n - p\)) 잔차들끼리는 서로 독립이 아니다.

잔차제곱합의 분포¶

추정량과 SSE의 독립성¶

개별 계수에 대한 t-통계량¶

\(\hat{\beta}_j\)의 정규분포, \(\text{SSE}/\sigma^2\)의 카이제곱분포, 그리고 이 둘의 독립성을 결합하면 t-분포가 나온다.

여러 계수를 검정하는 F-통계량¶

계수의 부분집합이 동시에 0인지 검정하려면, 계수 \(q\)인 \(q \times p\) 행렬 \(\mathbf{C}\)에 대해 \(H_0: \mathbf{C}\boldsymbol{\beta} = \mathbf{0}\)을 검정하는 것을 생각한다.

특수한 경우: 전체 F-검정¶

(\(\mathbf{X}\)의 첫 열이 절편이라 가정하고) \(H_0: \beta_1 = \beta_2 = \cdots = \beta_{p-1} = 0\)(절편 \(\beta_0\)은 남겨 두고 모든 기울기가 0)을 검정하면 전체 F-통계량을 얻는다. 이는 정리 5에서 \(\mathbf{C}\)를 절편을 제외한 \(p-1\)개 좌표를 뽑는 \((p-1) \times p\) 행렬로 잡은 경우이며, 따라서 \(q = p - 1\)이다.

여기서 \(\text{SSR} = \hat{\mathbf{y}}^\top\hat{\mathbf{y}} - n\bar{Y}^2\)이 회귀제곱합이다. \(H_0\) 아래에서 \(F \sim F_{p-1,\,n-p}\)이다.

사영의 언어로는 \(\text{SSR} = \lVert(\mathbf{H} - \tfrac{1}{n}\mathbf{J})\mathbf{y}\rVert^2\)이다. 절편만 있는 모형의 사영 \(\tfrac{1}{n}\mathbf{J}\)를 \(\mathbf{H}\)에서 걷어낸 것이 계수 \(p-1\)인 직교사영이기 때문에 자유도가 \(p\)가 아니라 \(p-1\)이다.

가우스–마르코프 정리¶

정규성을 가정하지 않아도 최소제곱추정량은 최적성을 갖는다.

핵심 분포 요약¶

이 모든 분포는 세 가지 재료에 의존한다. (1) \(\hat{\boldsymbol{\beta}}\)이 정규벡터 \(\mathbf{y}\)의 선형함수이고, (2) \(\mathbf{M}\)이 대칭 멱등이며, (3) \(\mathbf{X}^\top\mathbf{M} = \mathbf{0}\)이므로 \(\hat{\boldsymbol{\beta}}\)과 \(\text{SSE}\)가 독립이라는 것이다.

연습문제¶

import numpy as np

rng = np.random.default_rng(0)
n, p, sig = 30, 3, 2.0
X = np.column_stack([np.ones(n), rng.normal(size=(n, p - 1))])
beta = np.array([1., 2., -1.])
XtXi = np.linalg.inv(X.T @ X)

B = 200_000
Y = X @ beta + rng.normal(0, sig, size=(B, n))
bhat = Y @ X @ XtXi

print("평균 모의:", bhat.mean(axis=0).round(4), "  참값:", beta)
print("\n공분산 모의:\n", np.cov(bhat.T).round(4))
print("이론 sigma^2 (X'X)^-1:\n", (sig**2 * XtXi).round(4))

평균 모의: [ 0.9998  1.9991 -1.0002]   참값: [ 1.  2. -1.]

공분산 모의:
 [[ 0.164   0.0536 -0.0876]
 [ 0.0536  0.1971 -0.1124]
 [-0.0876 -0.1124  0.2714]]
이론 sigma^2 (X'X)^-1:
 [[ 0.1629  0.0532 -0.0863]
 [ 0.0532  0.1963 -0.1122]
 [-0.0863 -0.1122  0.2705]]

import numpy as np

rng = np.random.default_rng(0)
n, p, sig = 30, 3, 2.0
X = np.column_stack([np.ones(n), rng.normal(size=(n, p - 1))])
beta = np.array([1., 2., -1.])
XtXi = np.linalg.inv(X.T @ X)
x0 = np.array([1., 0.5, -0.3])

Y = X @ beta + rng.normal(0, sig, size=(200_000, n))
bhat = Y @ X @ XtXi
mu0 = bhat @ x0

print(f"Var(mu0) 모의 {mu0.var():.5f}   이론 {sig**2 * x0 @ XtXi @ x0:.5f}")
print(f"새 관측 예측오차 분산 이론 {sig**2 * (1 + x0 @ XtXi @ x0):.5f}")

Var(mu0) 모의 0.37766   이론 0.37503
새 관측 예측오차 분산 이론 4.37503

import numpy as np

rng = np.random.default_rng(2)
n = 30
x1 = rng.normal(size=n)
y = 1. + 2. * x1 + rng.normal(0, 1., n)

def fit(X):
    b = np.linalg.lstsq(X, y, rcond=None)[0]
    e = y - X @ b
    sse = e @ e
    sst = ((y - y.mean()) ** 2).sum()
    r2 = 1 - sse / sst
    adj = 1 - (sse / (len(y) - X.shape[1])) / (sst / (len(y) - 1))
    return r2, adj, np.linalg.inv(X.T @ X)[1, 1]

X1 = np.column_stack([np.ones(n), x1])
for k in (0, 1, 5, 10):
    noise = rng.normal(size=(n, k))
    X = np.column_stack([X1, noise]) if k else X1
    r2, adj, v = fit(X)
    print(f"잡음변수 {k:>2}개:  R^2 = {r2:.4f}   수정 R^2 = {adj:.4f}   "
          f"Var(b1) 인자 = {v:.4f}")

잡음변수  0개:  R^2 = 0.8446   수정 R^2 = 0.8390   Var(b1) 인자 = 0.0331
잡음변수  1개:  R^2 = 0.8451   수정 R^2 = 0.8336   Var(b1) 인자 = 0.0337
잡음변수  5개:  R^2 = 0.8674   수정 R^2 = 0.8328   Var(b1) 인자 = 0.0471
잡음변수 10개:  R^2 = 0.9025   수정 R^2 = 0.8429   Var(b1) 인자 = 0.0726

import numpy as np

rng = np.random.default_rng(0)
n, p, sig = 30, 3, 2.0
X = np.column_stack([np.ones(n), rng.normal(size=(n, p - 1))])
beta = np.array([1., 2., -1.])

w = np.linspace(1, 5, n)                    # 분산이 5배까지 커진다
V = np.diag(w)
Vinv = np.linalg.inv(V)

B = 200_000
Y = X @ beta + rng.normal(size=(B, n)) * np.sqrt(w) * sig

b_ols = Y @ X @ np.linalg.inv(X.T @ X)
b_gls = (Y @ Vinv @ X) @ np.linalg.inv(X.T @ Vinv @ X)

print("OLS 평균:", b_ols.mean(0).round(4), " 분산:", b_ols.var(0).round(5))
print("GLS 평균:", b_gls.mean(0).round(4), " 분산:", b_gls.var(0).round(5))

OLS 평균: [ 0.9997  1.9983 -1.0009]  분산: [0.48438 0.62429 0.89718]
GLS 평균: [ 0.9996  1.9988 -0.9999]  분산: [0.38248 0.51366 0.77013]

정리하며¶

최소제곱의 행렬 정식화는 단순회귀의 표본분포 이론을 하나의 통합된 틀로 압축한다. 추정량 \(\hat{\boldsymbol{\beta}}\)은 공분산이 \(\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}\)인 다변량 정규분포를 따르고, 잔차제곱합 \(\text{SSE}/\sigma^2\)은 자유도 \(n - p\)인 카이제곱을 따르며, 이 둘은 독립이다. 최소제곱의 직교사영 구조와 오차의 정규성에서 따라오는 이 세 사실이 개별 계수에 대한 모든 표준 t-검정과 계수 집합에 대한 F-검정을 만들어낸다. 나아가 가우스–마르코프 정리는 정규성 가정이 없어도 최소제곱이 선형불편추정량 가운데 최적임을 보여준다.