표본분포 (단순 최소제곱)¶

먼저 추정 분산을 계산한다.

기울기의 표준오차는

이다. \(t_{0.025, 8} = 2.306\)을 쓰면 95% 신뢰구간은

이다.

최소제곱 기울기 추정량은

이다. 마지막 등호는 \(\sum(x_i - \bar{x})\bar{Y} = \bar{Y}\sum(x_i - \bar{x}) = 0\)을 쓴 것이다. \(c_i = (x_i - \bar{x})/S_{xx}\)로 두면 \(\hat{\beta}_1 = \sum c_i Y_i\)이다.

\(Y_i\)가 독립이고 \(\operatorname{Var}(Y_i) = \sigma^2\)이므로

이다. \(\square\)

잔차제곱합은 (\(\sigma^2\)으로 나눈 뒤) 모자 행렬 \(\mathbf{H}\)에 대해 \(\text{SSE} = \mathbf{Y}^\top(\mathbf{I} - \mathbf{H})\mathbf{Y}/\sigma^2\)로 쓸 수 있다. 행렬 \(\mathbf{I} - \mathbf{H}\)는 계수가 \(n - 2\)인 멱등행렬이다(단순회귀에서는 절편과 기울기에 대응하여 \(\text{rank}(\mathbf{H}) = 2\)이므로).

멱등행렬의 계수가 결과로 나오는 카이제곱분포의 자유도와 같다. 직관적으로는, \(\mathbf{Y}\)에 \(n\)개의 독립적인 정보 조각이 있는 상태에서 출발하지만 두 모수 \(\beta_0\)과 \(\beta_1\)을 적합하는 데 자유도 2를 "써버려" \(\sigma^2\)을 추정할 자유도로 \(n - 2\)가 남는다.

\(\bar{Y} = \sum_i \frac{1}{n}Y_i\)이고 \(\hat{\beta}_1 = \sum_i c_i Y_i\)이므로

이다. 절편 추정량도 반응의 선형함수이며, 독립 정규확률변수의 선형결합이므로 정규분포를 따른다.

가중치 \(d_i\)는 두 항등식을 만족한다. \(\sum_i c_i = 0\)에서 \(\sum_i d_i = 1\)이고, \(\sum_i c_i x_i = 1\)에서 \(\sum_i d_i x_i = \bar{x} - \bar{x} = 0\)이다. 이 둘이 곧 불편성이다.

분산은 \(\sum_i d_i^2\)만 계산하면 된다.

이므로 \(\operatorname{Var}(\hat{\beta}_0) = \sigma^2\sum_i d_i^2 = \sigma^2\!\left(\frac{1}{n} + \frac{\bar{x}^2}{S_{xx}}\right)\)이다. 교차항이 \(\sum_i c_i = 0\) 때문에 사라진 것이 정리 2의 증명에서 \(\operatorname{Cov}(\bar{Y}, \hat{\beta}_1) = 0\)을 쓴 것과 같은 사실임에 주의하라.

import numpy as np

rng = np.random.default_rng(4)
x = np.arange(1., 9.)
n, b0, b1, sig = len(x), 2., 3., 1.5
Sxx = ((x - x.mean()) ** 2).sum()

c = (x - x.mean()) / Sxx
d = 1. / n - x.mean() * c
print("sum d_i = 1 인가    :", np.isclose(d.sum(), 1.0))
print("sum d_i x_i = 0 인가:", np.isclose((d * x).sum(), 0.0))

X = np.column_stack([np.ones(n), x])
Y = b0 + b1 * x + rng.normal(0, sig, size=(200_000, n))
beta = Y @ X @ np.linalg.inv(X.T @ X)

print("가중치 d 로 계산한 값이 최소제곱해와 같은가:", np.allclose(Y @ d, beta[:, 0]))
print(f"E[b0]    모의 {beta[:, 0].mean():.4f}     이론 {b0:.4f}")
print(f"Var(b0)  모의 {beta[:, 0].var():.6f}   이론 {sig**2 * (1/n + x.mean()**2/Sxx):.6f}")

출력:

sum d_i = 1 인가    : True
sum d_i x_i = 0 인가: True
가중치 d 로 계산한 값이 최소제곱해와 같은가: True
E[b0]    모의 2.0030     이론 2.0000
Var(b0)  모의 1.356516   이론 1.366071

\(\bar{x} = 4.5\), \(S_{xx} = 42\)인 이 설계에서 절편의 분산은 \(1.366\)인데 기울기의 분산은 \(0.054\)에 불과하다. \(\bar{x}^2/S_{xx} = 0.482\)가 \(1/n = 0.125\)를 압도하기 때문이다. \(x\)가 원점에서 멀면 절편은 자료로부터 멀리 외삽한 값이 되고, 그 대가가 이 항에 그대로 나타난다. \(\square\)

\(\hat{\beta}_0 = \bar{Y} - \hat{\beta}_1\bar{x}\)이고 본문에서 \(\operatorname{Cov}(\bar{Y}, \hat{\beta}_1) = 0\)임을 보였으므로

이다. 따라서 무상관일 필요충분조건은 \(\bar{x} = 0\), 곧 예측변수가 중심화되어 있는 것이다.

직관. \(\bar{x} > 0\)이면 자료가 원점 오른쪽에 몰려 있다. 회귀직선은 \((\bar{x}, \bar{Y})\)를 반드시 지나므로, 기울기를 조금 키우면 그 지렛대 효과로 절편이 내려간다. 음의 상관이 생기는 이유다.

import numpy as np

rng = np.random.default_rng(0)
x = np.arange(1., 9.)
n, b0, b1, sig = len(x), 2., 3., 1.5
Sxx = ((x - x.mean()) ** 2).sum()

X = np.column_stack([np.ones(n), x])
XtXi = np.linalg.inv(X.T @ X)
Y = b0 + b1 * x + rng.normal(0, sig, size=(200_000, n))
beta = Y @ X @ XtXi

print("Cov 모의실험:", round(np.cov(beta.T)[0, 1], 6))
print("Cov 이론값  :", round(-sig**2 * x.mean() / Sxx, 6))

xc = x - x.mean()                      # 중심화하면
Xc = np.column_stack([np.ones(n), xc])
beta_c = (b0 + b1 * x + rng.normal(0, sig, size=(200_000, n))) @ Xc @ np.linalg.inv(Xc.T @ Xc)
print("중심화 후 Cov:", round(np.cov(beta_c.T)[0, 1], 6))

출력:

Cov 모의실험: -0.241028
Cov 이론값  : -0.241071
중심화 후 Cov: 0.000156

중심화하면 공분산이 0에 가까워진다. 예측변수를 중심화하는 실무 관행의 한 가지 근거가 이것이다. 절편과 기울기의 추정이 서로 얽히지 않는다. \(\square\)

\(\hat{\mu}_0 = \bar{Y} + \hat{\beta}_1(x_0 - \bar{x})\)로 다시 쓰는 것이 요령이다(\(\hat{\beta}_0 = \bar{Y} - \hat{\beta}_1\bar{x}\)를 대입하면 된다). \(\bar{Y}\)와 \(\hat{\beta}_1\)이 무상관이므로

이다.

분산은 \((x_0 - \bar{x})^2\)에 대해 증가하므로 \(x_0 = \bar{x}\)에서 최소이고 그 값은 \(\sigma^2/n\)이다. 자료의 중심에서 멀어질수록 추정이 부정확해지며, 그 증가는 이차식이다. 신뢰띠가 가운데가 잘록한 나비 모양이 되는 이유다.

예측구간과의 차이. 새 관측값 \(Y_0\) 자체를 예측할 때는 오차항 \(\varepsilon_0\)의 변동이 더해진다.

괄호 안의 \(1\)이 결정적이다. \(n \to \infty\)이면 신뢰구간의 폭은 0으로 가지만 예측구간의 폭은 \(2z_{\alpha/2}\sigma\)로 남는다. 자료를 아무리 모아도 개별 관측의 무작위성은 사라지지 않는다. \(\square\)

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
x = np.arange(1., 9.)
n, b0, b1, sig = len(x), 2., 3., 1.5
Sxx = ((x - x.mean()) ** 2).sum()

B = 200_000
Y = b0 + b1 * x + rng.normal(0, sig, size=(B, n))
X = np.column_stack([np.ones(n), x])
beta = Y @ X @ np.linalg.inv(X.T @ X)
b1_hat = beta[:, 1]

print(f"평균   모의 {b1_hat.mean():.4f}   이론 {b1:.4f}")
print(f"분산   모의 {b1_hat.var():.6f}   이론 {sig**2/Sxx:.6f}")

# 정규성: 표준화한 뒤 분위수를 비교한다
z = (b1_hat - b1) / np.sqrt(sig**2 / Sxx)
for q in (0.05, 0.25, 0.5, 0.75, 0.95):
    print(f"  q={q:<5} 모의 {np.quantile(z, q):+.4f}   N(0,1) {stats.norm.ppf(q):+.4f}")

출력:

평균   모의 2.9992   이론 3.0000
분산   모의 0.053105   이론 0.053571
  q=0.05  모의 -1.6431   N(0,1) -1.6449
  q=0.25  모의 -0.6756   N(0,1) -0.6745
  q=0.5   모의 -0.0028   N(0,1) +0.0000
  q=0.75  모의 +0.6678   N(0,1) +0.6745
  q=0.95  모의 +1.6347   N(0,1) +1.6449

평균과 분산이 이론값과 맞고 분위수도 표준정규와 일치한다.

왜 정확히 정규인가. \(\hat{\beta}_1 = \sum_i c_i Y_i\)가 독립 정규확률변수의 선형결합이기 때문이다. 중심극한정리에 기댄 근사가 아니라 \(n\)이 작아도 정확히 성립한다. 여기서도 \(n = 8\)뿐이다. 다만 이는 오차가 정규일 때의 이야기이며, 오차가 정규가 아니면 \(\hat{\beta}_1\)의 정규성은 \(n\)이 커질 때의 근사로만 성립한다. \(\square\)

두 제약은 최소제곱의 정규방정식 그 자체다. \(\sum_i (Y_i - \beta_0 - \beta_1 x_i)^2\)을 \(\beta_0\)과 \(\beta_1\)로 각각 편미분해 0으로 두면

를 얻는다.

import numpy as np

rng = np.random.default_rng(2)
x = np.arange(1., 9.)
y = 2. + 3. * x + rng.normal(0, 1.5, size=len(x))

X = np.column_stack([np.ones(len(x)), x])
beta = np.linalg.lstsq(X, y, rcond=None)[0]
e = y - X @ beta

print("sum e_i    =", round(e.sum(), 12))
print("sum x_i e_i =", round((x * e).sum(), 12))
print("잔차 개수 =", len(e), " 자유로운 성분 =", len(e) - 2)

출력:

sum e_i    = 0.0
sum x_i e_i = 0.0
잔차 개수 = 8  자유로운 성분 = 6

기하적 의미. 잔차벡터 \(\mathbf{e}\)는 \(\mathbf{1}\)과 \(\mathbf{x}\) 둘 다에 직교한다. 곧 \(n\)차원 공간에서 두 방향이 막힌 \((n-2)\)차원 부분공간에 놓인다. 잔차가 \(n\)개 있어도 자유롭게 움직일 수 있는 것은 \(n-2\)개뿐이며, 그래서 \(s^2 = \text{SSE}/(n-2)\)가 불편추정량이 된다.

예측변수가 \(p\)개인 다중회귀에서는 제약이 \(p\)개(절편 포함)가 되어 자유도가 \(n - p\)다. \(\operatorname{tr}(\mathbf{I} - \mathbf{H}) = n - p\)가 같은 사실의 행렬 표현이다. \(\square\)

import numpy as np

rng = np.random.default_rng(3)
x = np.arange(1., 9.)
n, b0, b1, sig = len(x), 2., 3., 1.5

B = 200_000
X = np.column_stack([np.ones(n), x])
Y = b0 + b1 * x + rng.normal(0, sig, size=(B, n))
beta = Y @ X @ np.linalg.inv(X.T @ X)
sse = ((Y - beta @ X.T) ** 2).sum(axis=1)
s2 = sse / (n - 2)

print(f"E[s^2]  모의 {s2.mean():.4f}    sigma^2 = {sig**2:.4f}")
print(f"Var((n-2)s^2/sig^2) 모의 {(sse/sig**2).var():.4f}   chi2_{n-2} 이론 {2*(n-2)}")
print(f"corr(beta1_hat, s^2) = {np.corrcoef(beta[:,1], s2)[0,1]:+.5f}  (독립이면 0)")

출력:

E[s^2]  모의 2.2488    sigma^2 = 2.2500
Var((n-2)s^2/sig^2) 모의 11.9591   chi2_6 이론 12
corr(beta1_hat, s^2) = -0.00447  (독립이면 0)

\(s^2\)은 불편이고, \((n-2)s^2/\sigma^2\)의 분산이 \(\chi^2_{6}\)의 분산 \(12\)와 맞으며, \(\hat{\beta}_1\)과의 상관이 사실상 0이다.

왜 필요한가. t-통계량은

의 꼴이다. \(t\) 분포의 정의는 표준정규 \(Z\)와 카이제곱 \(V\)가 독립일 것을 요구한다. 위 모의실험이 확인한 두 가지가 정확히 그 조건이다. 독립성이 없으면 이 비는 \(t\) 분포를 따르지 않는다.

독립성은 기하에서 나온다. \(\hat{\boldsymbol{\beta}}\)은 \(\operatorname{col}(\mathbf{X})\) 안의 사영으로 결정되고 \(s^2\)은 그 직교여공간의 잔차로 결정되는데, 정규분포에서 직교하는 성분은 독립이기 때문이다.

같은 이야기를 이미 한 번 했다. 4.2절 정규분포의 연습문제 7에서 \(\bar{X}\)와 \(S^2\)의 독립성을 증명할 때, 첫 행이 \(\mathbf{1}/\sqrt{n}\)인 직교행렬 \(\mathbf{Q}\)로 회전시켜 \(\bar{X}\)는 \(Y_1\)만의 함수, \(S^2\)은 \(Y_2, \dots, Y_n\)만의 함수임을 보였다. 그 \(Y_1\) 방향이 여기서는 \(\operatorname{col}(\mathbf{X})\)이고 나머지 \(n-1\)개 방향이 그 직교여공간이다. 절편만 있는 모형이 정확히 그 경우이며, 설명변수를 넣으면 \(\operatorname{col}(\mathbf{X})\)의 차원이 1에서 2로 늘어나 잔차 쪽이 \(n-1\)에서 \(n-2\)로 줄어든다. 직교변환으로 쓰던 논법을 사영의 언어로 옮겨 적은 것이 전부다. \(\square\)

\(S_{xx} = \sum_i (x_i - \bar{x})^2\)을 최대화하려면 자료를 양 끝으로 몰아야 한다. \(x = 0\)에 5개, \(x = 10\)에 5개를 두면 \(S_{xx} = 10 \times 5^2 = 250\)으로 최댓값이다.

import numpy as np

designs = {
    "양 끝에 몰기": np.array([0.]*5 + [10.]*5),
    "균등 배치   ": np.linspace(0, 10, 10),
    "가운데 몰기 ": np.array([4., 4.5, 4.5, 5., 5., 5., 5., 5.5, 5.5, 6.]),
}
for name, x in designs.items():
    Sxx = ((x - x.mean()) ** 2).sum()
    print(f"{name}: Sxx = {Sxx:7.2f},  Var(b1) = sigma^2 / {Sxx:.2f}")

출력:

양 끝에 몰기: Sxx =  250.00,  Var(b1) = sigma^2 / 250.00
균등 배치   : Sxx =  101.85,  Var(b1) = sigma^2 / 101.85
가운데 몰기 : Sxx =    3.00,  Var(b1) = sigma^2 / 3.00

양 끝 설계의 \(S_{xx} = 250\)은 균등 배치의 \(101.85\)보다 약 2.5배 크다. 기울기의 분산은 그만큼 줄어든다.

그런데 이 설계는 위험하다. 가운데에 자료가 하나도 없으므로 직선성 가정을 검증할 수 없다. 참 관계가 곡선이어도 두 점만으로는 언제나 완벽한 직선이 그려진다. 곡률을 탐지할 힘이 0이다.

실무의 절충은 이렇다.

• 직선 모형을 확신할 수 있다면(물리 법칙 등) 양 끝 설계가 효율적이다.
• 모형을 검증해야 한다면 가운데에도 점을 남긴다. 예컨대 양 끝에 40%씩, 가운데에 20%를 배치한다.

이것이 실험계획의 기본 긴장이다. 추정의 효율과 모형 검증 능력은 맞바꾸는 관계다. \(\square\)