가중최소제곱¶

(1) 해석적으로. 가중제곱합을 그대로 다시 적는다. \(w_i > 0\)이므로 \(w_i = (\sqrt{w_i})^2\)이고

이다. 목적함수가 글자 그대로 같으므로 최소점도 같다. 정규방정식으로 봐도 \(\tilde{\mathbf X}^\top \tilde{\mathbf X} = \mathbf X^\top \mathbf W \mathbf X\), \(\tilde{\mathbf X}^\top \tilde{\mathbf y} = \mathbf X^\top \mathbf W \mathbf y\)라 바로 확인된다. 왜 이 변환이 쓸모 있는가 하면, \(w_i = 1/\sigma_i^2\)일 때 변환된 오차가

로 등분산이 되기 때문이다. WLS는 새로운 이론이 아니라 자료를 등분산으로 만든 뒤의 OLS다.

단순회귀의 닫힌 꼴은 OLS의 것에서 모든 합과 평균에 \(w\)를 끼워 넣은 모양이다. 가중평균을

로 두면

이다. 중심을 산술평균이 아니라 가중평균으로 잡는다는 것이 핵심이다.

(2) 척도 불변성. \(\mathbf W \to c\mathbf W\)를 넣으면

로 \(c\)가 깨끗이 상쇄된다. 추정값은 가중치의 절대 크기가 아니라 비율만 본다. 그런데 분산 공식은 \((c\mathbf X^\top \mathbf W \mathbf X)^{-1} = \tfrac1c (\mathbf X^\top \mathbf W \mathbf X)^{-1}\)로 \(1/c\)배가 된다.

실무에서 이것은 덫이다. \(\sigma_i\)를 몰라 "\(x\)에 반비례하게" 같은 모양만 맞는 가중치를 쓰면 추정값은 제대로 나오지만 표준오차는 임의의 상수만큼 틀린다. \((\mathbf X^\top \mathbf W \mathbf X)^{-1}\)이 분산이려면 \(w_i\)가 \(1/\sigma_i^2\)에 비례하는 것으로는 모자라고 정확히 같아야 한다. 비율만 아는 경우에는 가중잔차에서 척도를 다시 추정해 \(\hat s^2 (\mathbf X^\top \mathbf W \mathbf X)^{-1}\)을 써야 하며, statsmodels.WLS가 하는 일이 그것이다.

import numpy as np

def ols_fit(X, y):
    """보통최소제곱. 모든 관측값에 같은 무게를 준다."""
    beta = np.linalg.lstsq(X, y, rcond=None)[0]
    return beta

def wls_fit(X, y, w):
    """가중최소제곱. 관측값마다 다른 무게 w 를 준다.

    정규방정식의 X'X 와 X'y 사이에 가중행렬 W 가 끼어드는 것이 전부다.
    분산이 큰 관측값에 작은 무게를 주면 추정의 분산이 줄어든다.
    """
    W = np.diag(w)
    XtW = X.T @ W
    beta = np.linalg.solve(XtW @ X, XtW @ y)
    return beta

아래 보기 2가 만드는 자료로 세 가지를 확인한다. 변환 동등성, 닫힌 꼴, 척도 불변성이다.

import numpy as np

# 보기 2 가 쓸 자료를 먼저 만들어 둔다 (씨앗이 같으므로 같은 자료다).
np.random.seed(42)
n = 120
x = np.random.uniform(1, 10, n)
sigma = 0.5 + 1.5 * x
y = 3.0 + 2.0 * x + np.random.normal(0, sigma)
X = np.column_stack([np.ones(n), x])
w = 1.0 / sigma ** 2

# 1) sqrt(w) 변환 뒤의 OLS 와 같은가
Xt = X * np.sqrt(w)[:, None]
yt = y * np.sqrt(w)
print("변환 뒤 OLS =", np.linalg.lstsq(Xt, yt, rcond=None)[0])
print("wls_fit     =", wls_fit(X, y, w))

# 2) 단순회귀의 닫힌 꼴
xw = (w * x).sum() / w.sum()
yw = (w * y).sum() / w.sum()
b1 = (w * (x - xw) * (y - yw)).sum() / (w * (x - xw) ** 2).sum()
print(f"\n닫힌 꼴 기울기 = {b1:.12f}")
print(f"wls_fit 기울기 = {wls_fit(X, y, w)[1]:.12f}")

# 3) 가중치를 1000 배 해도 추정값은 그대로, 분산 공식은 1/1000 배
se = np.sqrt(np.diag(np.linalg.inv(X.T @ np.diag(w) @ X)))
se_1000 = np.sqrt(np.diag(np.linalg.inv(X.T @ np.diag(1000 * w) @ X)))
print(f"\nw 를 1000 배 한 추정값 = {wls_fit(X, y, 1000 * w)}")
print(f"SE(기울기):  w 그대로 {se[1]:.6f},  1000w {se_1000[1]:.6f}")
print(f"비 = {se[1] / se_1000[1]:.4f},  sqrt(1000) = {np.sqrt(1000):.4f}")

출력:

변환 뒤 OLS = [3.34999074 2.01054088]
wls_fit     = [3.34999074 2.01054088]

닫힌 꼴 기울기 = 2.010540878598
wls_fit 기울기 = 2.010540878598

w 를 1000 배 한 추정값 = [3.34999074 2.01054088]
SE(기울기):  w 그대로 0.258867,  1000w 0.008186
비 = 31.6228,  sqrt(1000) = 31.6228

세 가지가 모두 맞는다. 변환 뒤의 OLS와 wls_fit이 같은 값을 주고, 닫힌 꼴이 소수점 열두째 자리까지 일치하며, 가중치를 \(1000\)배 해도 추정값은 한 자리도 바뀌지 않는다. 그런데 공식이 주는 표준오차는 \(0.258867\)에서 \(0.008186\)으로 \(\sqrt{1000} = 31.62\)배 작아진다. 자료도 추정값도 그대로인데 표준오차만 서른 배 작아졌다면 그 수가 거짓이라는 뜻이고, 거짓인 쪽은 가중치의 척도가 틀린 쪽이다.

(1) 해석적으로. \(\tilde\beta_1 = \sum_i a_i y_i\)가 불편이려면

이 모든 \((\beta_0, \beta_1)\)에 대해 성립해야 하므로 두 제약 \(\sum_i a_i = 0\)과 \(\sum_i a_i x_i = 1\)이 나온다. 오차가 독립이므로 분산은 \(\sum_i a_i^2 \sigma_i^2\)이다.

이제 둘째 제약을 첫째 제약으로 다듬는다. 어떤 상수 \(\bar x_w\)를 빼도 되므로

이다. 이 합을 두 조각으로 갈라 코시–슈바르츠를 적용한다.

오른쪽 둘째 근호 안이 \(w_i = 1/\sigma_i^2\)로 \(\sum_i w_i (x_i - \bar x_w)^2 = S_w\)다. 양변을 제곱해 옮기면

를 얻는다. 여기서 \(\bar x_w\)는 아직 아무 상수여도 좋았다. 부등식은 어느 \(c\)에 대해서도 \(\operatorname{Var} \ge 1/S_w(c)\)로 성립하므로 가장 센 하한은 \(S_w(c) = \sum_i w_i (x_i - c)^2\)을 가장 작게 만드는 \(c\)에서 나온다. 그 \(c\)가 바로 가중평균 \(\bar x_w = \sum_i w_i x_i / \sum_i w_i\)다(\(S_w(c)\)를 \(c\)로 미분해 0으로 두면 나온다).

등호는 두 벡터가 평행할 때, 곧 어떤 \(\lambda\)에 대해

일 때다. 제약 \(\sum_i a_i x_i = 1\)에 넣으면 \(\lambda \sum_i w_i (x_i - \bar x_w) x_i = \lambda S_w = 1\)이므로 \(\lambda = 1/S_w\)이고

이다. 이것이 바로 보기 1의 닫힌 꼴 \(\hat\beta_1 = \sum_i w_i(x_i - \bar x_w)(y_i - \bar y_w)/S_w\)의 가중치다. WLS가 최소분산 선형불편추정량이고, 그 분산이 정확히 \(1/S_w\)다.

등호 조건이 유일하다는 점에 주목하라. 가중치가 \(1/\sigma_i^2\)에 비례하지 않으면 등호가 깨지고 분산이 반드시 커진다. "분산이 큰 쪽에 작은 가중치를 주면 좋다"가 아니라 "정확히 \(1/\sigma_i^2\)에 비례할 때만 최적" 이다.

(2) 수치적으로.

np.random.seed(42)
n = 120
x = np.random.uniform(1, 10, n)

# 잡음의 크기가 x 에 따라 커진다. 전형적인 이분산 자료다.
sigma = 0.5 + 1.5 * x
y = 3.0 + 2.0 * x + np.random.normal(0, sigma)

X = np.column_stack([np.ones(n), x])

# 두 방법 모두 불편추정량이라 계수 자체는 비슷하게 나온다.
beta_ols = ols_fit(X, y)

# 분산의 역수를 무게로 쓰는 것이 최적이다. 여기서는 참 sigma 를 알고 있어
# 그대로 썼지만, 실제로는 sigma 도 자료에서 추정해야 한다.
w = 1.0 / sigma ** 2
beta_wls = wls_fit(X, y, w)

가중치를 바꿔 가며 참 분산을 계산한다. \(\sigma_i\)를 알고 있으므로 모의실험 없이 공식으로 끝난다.

import numpy as np

def slope_var(v):
    """가중치 v 로 적합한 기울기의 참 분산과 그 선형가중치 a."""
    xv = (v * x).sum() / v.sum()
    S = (v * (x - xv) ** 2).sum()
    a = v * (x - xv) / S
    return (a ** 2 * sigma ** 2).sum(), a

print(f"{'가중치':>22}{'Var':>11}{'SE':>10}{'sum a':>10}{'sum a*x':>10}")
for lab, v in [("OLS (w = 1)", np.ones(n)),
               ("WLS 최적 (1/sigma^2)", w),
               ("1/sigma (덜 줌)", 1 / sigma),
               ("1/sigma^4 (지나침)", 1 / sigma ** 4),
               ("1/x (모양만 비슷)", 1 / x)]:
    V, a = slope_var(v)
    print(f"{lab:>22}{V:>11.6f}{np.sqrt(V):>10.6f}"
          f"{a.sum():>10.1e}{(a * x).sum():>10.4f}")

xw = (w * x).sum() / w.sum()
Sw = (w * (x - xw) ** 2).sum()
print(f"\n하한 1/S_w = {1 / Sw:.6f}   (SE = {1 / np.sqrt(Sw):.6f})")

출력:

                   가중치        Var        SE     sum a   sum a*x
           OLS (w = 1)   0.121722  0.348887   2.4e-17    1.0000
    WLS 최적 (1/sigma^2)   0.067012  0.258867  -4.5e-17    1.0000
         1/sigma (덜 줌)   0.080300  0.283372  -2.1e-16    1.0000
       1/sigma^4 (지나침)   0.195189  0.441802   3.3e-16    1.0000
          1/x (모양만 비슷)   0.077782  0.278894   5.2e-17    1.0000

하한 1/S_w = 0.067012   (SE = 0.258867)

부등식이 다섯 줄 모두에서 지켜진다. 어느 가중치든 두 제약 \(\sum a_i = 0\)(\(10^{-16}\) 수준)과 \(\sum a_i x_i = 1\)을 만족하므로 모두 불편이고, 분산은 전부 하한 \(0.067012\) 이상이다. 그리고 등호는 \(w = 1/\sigma^2\) 한 줄에서만 성립한다. 소수점 여섯째 자리까지 하한과 같다.

가중치를 지나치게 주면 OLS보다 나빠진다. \(1/\sigma^4\)의 분산이 \(0.1952\)로 OLS의 \(0.1217\)보다 \(1.6\)배 크다. 정밀한 관측 몇 개에 발언권을 몰아주는 바람에, 사실상 그 몇 개로만 직선을 긋는 꼴이 된 것이다. "가중치를 주면 좋아진다"가 아니라 "올바른 가중치를 주어야 좋아진다"는 (1)의 등호 조건이 여기서 값을 치른다.

반대로 \(1/\sigma\)나 \(1/x\)처럼 방향만 맞는 가중치는 \(0.0803\), \(0.0778\)로 OLS의 \(0.1217\)보다는 좋고 최적의 \(0.0670\)보다는 나쁘다. 실무에서 \(\sigma_i\)를 정확히 알 길이 없으므로 대개 이 중간 어딘가에 머물게 되며, 그 사실이 WLS 대신 로버스트 표준오차를 쓰는 쪽을 택하게 만드는 이유 가운데 하나다.

(1) 해석적으로. \(\hat{\boldsymbol\beta}_{\text{OLS}} = (\mathbf X^\top \mathbf X)^{-1}\mathbf X^\top \mathbf y\)에 \(\mathbf y = \mathbf X \boldsymbol\beta + \boldsymbol\varepsilon\)을 넣으면

이다. 선형변환의 공분산 법칙 \(\operatorname{Var}(\mathbf A \boldsymbol\varepsilon) = \mathbf A \boldsymbol\Sigma \mathbf A^\top\)를 그대로 적용하면

을 얻는다. 가운데에 \(\boldsymbol\Sigma\)가 끼인 모양이라 샌드위치 공식이라 부른다. \(\boldsymbol\Sigma = \sigma^2 \mathbf I\)이면 \(\sigma^2\)이 밖으로 나오고 양쪽의 \((\mathbf X^\top \mathbf X)^{-1}\mathbf X^\top \mathbf X\)가 상쇄되어 익숙한 \(\sigma^2 (\mathbf X^\top \mathbf X)^{-1}\)로 돌아간다.

OLS가 보고하는 것은 \(s^2(\mathbf X^\top\mathbf X)^{-1}\)인데, 이것은 빵만 있고 속이 없는 셈이다. \(\boldsymbol\Sigma\)가 \(\sigma_i^2\)마다 다른데 그것을 하나의 \(s^2\)으로 뭉갰기 때문이다.

어느 쪽으로 틀리는지는 미리 말할 수 없다. 큰 \(\sigma_i^2\)이 지렛대가 큰 자리(설계행렬의 끝)에 몰리면 샌드위치의 속이 두꺼워져 참값이 보고값보다 커지고, 가운데에 몰리면 반대가 된다. 같은 자료에서도 계수마다 방향이 다를 수 있다는 것을 (2)에서 보게 된다. 이것이 이 페이지가 "OLS 표준오차는 너무 클 수도 작을 수도 있다"고 적는 이유다.

WLS의 분산 \((\mathbf X^\top \mathbf W \mathbf X)^{-1}\)은 사정이 다르다. \(\mathbf W = \boldsymbol\Sigma^{-1}\)이면 샌드위치가 \((\mathbf X^\top \mathbf W \mathbf X)^{-1}\mathbf X^\top \mathbf W \boldsymbol\Sigma \mathbf W \mathbf X (\mathbf X^\top \mathbf W \mathbf X)^{-1} = (\mathbf X^\top \mathbf W \mathbf X)^{-1}\)로 저절로 접힌다. 올바른 가중치를 쓰면 보고값이 곧 참값이다.

(2) 수치적으로.

# 여기서 갈린다. OLS 의 표준오차 공식은 등분산을 전제하므로, 이분산
# 자료에서는 그 값 자체를 믿을 수 없다.
resid_ols = y - X @ beta_ols
s2_ols = np.sum(resid_ols ** 2) / (n - 2)
se_ols = np.sqrt(np.diag(s2_ols * np.linalg.inv(X.T @ X)))

# WLS 의 표준오차는 이분산을 제대로 셈에 넣은 값이고, OLS 의 것보다 작다.
# 같은 자료에서 더 정확한 결론을 얻는다는 뜻이다.
W = np.diag(w)
XtWX_inv = np.linalg.inv(X.T @ W @ X)
se_wls = np.sqrt(np.diag(XtWX_inv))

print(f"OLS:  intercept={beta_ols[0]:.3f} (SE={se_ols[0]:.3f}), "
      f"slope={beta_ols[1]:.3f} (SE={se_ols[1]:.3f})")
print(f"WLS:  intercept={beta_wls[0]:.3f} (SE={se_wls[0]:.3f}), "
      f"slope={beta_wls[1]:.3f} (SE={se_wls[1]:.3f})")

출력(참값은 절편 3.0, 기울기 2.0):

OLS:  intercept=4.083 (SE=1.851), slope=1.842 (SE=0.312)
WLS:  intercept=3.350 (SE=0.841), slope=2.011 (SE=0.259)

이제 샌드위치 공식으로 참값을 구하고 모의실험으로 판정한다.

import numpy as np

XtX_inv = np.linalg.inv(X.T @ X)
Sigma = np.diag(sigma ** 2)
V_true = XtX_inv @ X.T @ Sigma @ X @ XtX_inv      # 샌드위치 공식
se_true = np.sqrt(np.diag(V_true))

print(f"{'':>12}{'절편':>12}{'기울기':>12}")
print(f"{'OLS 보고':>12}{se_ols[0]:>12.4f}{se_ols[1]:>12.4f}")
print(f"{'OLS 참값':>12}{se_true[0]:>12.4f}{se_true[1]:>12.4f}")
print(f"{'WLS':>12}{se_wls[0]:>12.4f}{se_wls[1]:>12.4f}")

# 같은 X 를 고정한 채 y 만 20000번 다시 뽑아 확인한다.
r = np.random.default_rng(7)
E = r.normal(0, 1, (20_000, n)) * sigma
B_ols = (np.linalg.inv(X.T @ X) @ X.T @ E.T).T
W = np.diag(w)
B_wls = (np.linalg.inv(X.T @ W @ X) @ X.T @ W @ E.T).T
print(f"{'모의 OLS':>12}{B_ols[:, 0].std(ddof=1):>12.4f}{B_ols[:, 1].std(ddof=1):>12.4f}")
print(f"{'모의 WLS':>12}{B_wls[:, 0].std(ddof=1):>12.4f}{B_wls[:, 1].std(ddof=1):>12.4f}")

출력:

                      절편         기울기
      OLS 보고      1.8509      0.3122
      OLS 참값      1.3636      0.3489
         WLS      0.8406      0.2589
      모의 OLS      1.3581      0.3478
      모의 WLS      0.8366      0.2574

샌드위치 공식이 맞는다. 참값 \(1.3636\)과 \(0.3489\)를 모의실험이 \(1.3581\)과 \(0.3478\)로 재현한다(차이 \(0.4\%\)). WLS도 \(0.8406\) 대 \(0.8366\)으로 맞는다. 보고된 OLS 표준오차만 어느 쪽과도 맞지 않는다.

같은 자료에서 방향이 반대다. 절편에서는 보고값 \(1.8509\)가 참값 \(1.3636\)보다 \(36\%\) 크고, 기울기에서는 보고값 \(0.3122\)가 참값 \(0.3489\)보다 \(11\%\) 작다. (1)에서 "미리 말할 수 없다"고 한 것의 가장 선명한 예다. 기울기만 보고 "이분산이면 표준오차가 과소평가된다"고 외운 사람은 절편에서 틀리게 된다.

실질적으로 위험한 쪽은 기울기다. 참 표준오차가 \(0.3489\)인데 \(0.3122\)로 보고하면 \(t\)값이 \(12\%\) 부풀고 신뢰구간이 그만큼 좁아진다. 절편의 과대평가는 보수적인 쪽이라 덜 위험하다.

WLS와 견주는 올바른 비교. 이 페이지의 요지는 WLS가 더 효율적이라는 것인데, 그것을 보려면 WLS의 \(0.2589\)를 OLS 보고값 \(0.3122\)가 아니라 OLS 참값 \(0.3489\)와 견주어야 한다. 비는 \(0.3489/0.2589 = 1.35\)이고 분산으로는 \(1.82\)배다. 보기 2의 표가 준 \(0.121722/0.067012 = 1.82\)와 같은 수다. OLS 보고값은 틀린 수이므로 어떤 비교에도 쓰면 안 된다.

sigma_rev = 10 - 0.8 * x
y_rev = 3.0 + 2.0 * x + np.random.normal(0, sigma_rev)
w_rev = 1.0 / sigma_rev ** 2
beta_ols_rev = ols_fit(X, y_rev)
beta_wls_rev = wls_fit(X, y_rev, w_rev)

그렇다. 이분산이 존재하는 한 방향과 무관하게 WLS가 OLS보다 낫다. WLS는 올바른 가중치를 쓰므로 이분산 오차로 일반화된 Gauss-Markov 정리 아래에서 최소분산 선형불편추정량(BLUE)이 된다. \(\square\)

모든 \(i\)에 대해 \(w_i = c\)이면 \(\mathbf{W} = c\mathbf{I}\)이므로

OLS는 모든 관측값에 같은 가중치를 주는 WLS의 특수한 경우이며, 등분산 가정이 성립할 때 적절하다. \(\square\)

resid_ols = y - X @ beta_ols
log_resid_sq = np.log(resid_ols ** 2 + 1e-10)
gamma = np.linalg.lstsq(X, log_resid_sq, rcond=None)[0]
sigma_hat = np.sqrt(np.exp(X @ gamma))
w_feas = 1.0 / sigma_hat ** 2
beta_fwls = wls_fit(X, y, w_feas)

실행가능 WLS는 알려진 가중치 대신 추정된 가중치를 쓴다. 추정량은 일치성과 점근적 효율성을 갖지만, 가중치를 아는 WLS에 비해 작은 표본에서는 효율이 떨어질 수 있다. \(\square\)

일반화 Gauss-Markov 정리에 따르면 \(\mathrm{Var}(\boldsymbol{\varepsilon}) = \boldsymbol{\Sigma}\)일 때 \(\boldsymbol{\beta}\)의 BLUE는

\(\boldsymbol{\Sigma}\)가 대각행렬이면 \(\boldsymbol{\Sigma}^{-1} = \mathrm{diag}(1/\sigma_1^2, \ldots, 1/\sigma_n^2) = \mathbf{W}\)이다. 따라서 \(\hat{\boldsymbol{\beta}}_{\text{GLS}} = \hat{\boldsymbol{\beta}}_{\text{WLS}}\)이다. "최선"이란 모든 선형불편추정량 가운데 분산이 가장 작다는 뜻으로, 임의의 선형불편 \(\tilde{\boldsymbol{\beta}}\)와 임의의 방향 \(\mathbf{a}\)에 대해 \(\mathrm{Var}(\mathbf{a}^\top\hat{\boldsymbol{\beta}}_{\text{WLS}}) \leq \mathrm{Var}(\mathbf{a}^\top\tilde{\boldsymbol{\beta}})\)이다. \(\square\)

HC 표준오차는 계수 추정값은 그대로 두고 OLS의 표준오차만 교정한다.

여기서 \(\hat{\boldsymbol{\Omega}} = \mathrm{diag}(e_1^2, \ldots, e_n^2)\)이다(HC0 형태). 절충 관계는 이렇다. (1) HC 표준오차는 분산함수를 지정하지 않고도 임의의 이분산 아래에서 타당하지만 OLS 계수 자체는 비효율적이다. (2) WLS는 효율적인 추정을 주지만 올바른 가중함수를 지정해야 한다. (3) 분산 모형이 잘못 설정되면 WLS는 표준오차에 편향을 들여올 수 있는 반면, HC 표준오차는 여전히 로버스트하다. 실무에서 분산함수를 모를 때는 HC 표준오차가 선호된다. \(\square\)