OLS 회귀 출력 재현¶

import numpy as np
from scipy import stats

def fit_ols(X, y):
    """최소제곱 추정값과, 표준오차를 만드는 데 필요한 두 조각을 돌려준다.

    s 는 잔차의 표준편차(자유도 n-k), cov_matrix 는 (X'X)^-1 이다.
    계수 j 의 표준오차는 s * sqrt(cov_matrix[j, j]) 로 만들어진다.
    """
    n, k = X.shape
    beta_hat = np.linalg.inv(X.T @ X) @ X.T @ y
    y_hat = X @ beta_hat
    residuals = y - y_hat
    s = np.sqrt(np.sum(residuals ** 2) / (n - k))
    cov_matrix = np.linalg.inv(X.T @ X)
    return beta_hat, s, cov_matrix

# 보기 3 의 자료를 여기서 미리 읽어 둔다(보기 3 에서 다시 읽어도 결과는 같다).
import pandas as pd

url = ('https://raw.githubusercontent.com/justmarkham/'
       'scikit-learn-videos/8545c74961398def7724501648fd504dbf061b41/data/Advertising.csv')
data = pd.read_csv(url, usecols=[1, 2, 3, 4])
training_data = data.iloc[:int(len(data) * 0.7)]
y = np.array(training_data.Sales).reshape(-1, 1)
n = y.shape[0]
X = np.concatenate((np.ones((n, 1)), np.array(training_data.iloc[:, :-1])), axis=1)
k = X.shape[1]
beta_hat, s, cov_matrix = fit_ols(X, y)

residuals = y - X @ beta_hat
print(f"n = {n},  k = {k},  s = {s:.6f},  s^2 = {s ** 2:.6f}")
print("(X'X)^-1 의 대각 =", np.diag(cov_matrix).round(8))
print("표준오차 = s * sqrt(diag) =", (s * np.sqrt(np.diag(cov_matrix))).round(6))
print(f"X'e 의 최대 절대값 = {np.abs(X.T @ residuals).max():.3e}")
print(f"잔차의 합          = {residuals.sum():.3e}")

beta_lstsq = np.linalg.lstsq(X, y, rcond=None)[0]
print(f"inv 와 lstsq 의 최대 차이 = {np.abs(beta_lstsq - beta_hat).max():.3e}")
print(f"X'X 의 조건수 = {np.linalg.cond(X.T @ X):.4e},  X 의 조건수 = {np.linalg.cond(X):.4e}")

n = 140,  k = 4,  s = 1.736514,  s^2 = 3.015479
(X'X)^-1 의 대각 = [5.077908e-02 9.600000e-07 3.855000e-05 1.632000e-05]
표준오차 = s * sqrt(diag) = [0.39131  0.001701 0.010782 0.007015]
X'e 의 최대 절대값 = 1.247e-10
잔차의 합          = -1.790e-12
inv 와 lstsq 의 최대 차이 = 4.041e-14
X'X 의 조건수 = 2.0924e+05,  X 의 조건수 = 4.5743e+02

def regression_table(beta_hat, s, cov_matrix, n, k, var_names):
    """회귀 출력표를 직접 만든다.

    statsmodels 의 summary() 가 찍어 주는 계수 표와 같은 내용이다.
    계수 → 표준오차 → t → p-값 → 신뢰구간이 어떤 순서로 만들어지는지를
    보이려고 풀어 썼다.
    """
    df = n - k
    t_crit = stats.t(df).ppf(0.975)

    for name, j in zip(var_names, range(k)):
        coef = beta_hat[j, 0]
        v_j = cov_matrix[j, j]
        se = s * np.sqrt(v_j)
        t_stat = coef / se
        p_val = 2 * stats.t(df).sf(np.abs(t_stat))
        ci_lo = coef - t_crit * se
        ci_hi = coef + t_crit * se
        print(f"{name:10}  coef={coef:.4f}  SE={se:.3f}  "
              f"t={t_stat:.3f}  p={p_val:.3f}  "
              f"CI=({ci_lo:.3f}, {ci_hi:.3f})")

df = n - k
t_crit = stats.t(df).ppf(0.975)
print(f"자유도 = {df},  t 임계값 = {t_crit:.6f}")

names = ["Intercept", "TV", "Radio", "Newspaper"]
for j, nm in enumerate(names):
    se = s * np.sqrt(cov_matrix[j, j])
    t_stat = beta_hat[j, 0] / se
    p_val = 2 * stats.t(df).sf(abs(t_stat))
    lo, hi = beta_hat[j, 0] - t_crit * se, beta_hat[j, 0] + t_crit * se
    print(f"  {nm:10s} |t|>t* {str(abs(t_stat) > t_crit):5s}  "
          f"p<0.05 {str(p_val < 0.05):5s}  CI 가 0 을 배제 {str(not (lo < 0 < hi)):5s}  "
          f"(t={t_stat:+.4f}, p={p_val:.3e})")

se_tv = s * np.sqrt(cov_matrix[1, 1])
print(f"TV 의 표준오차: 참값 {se_tv:.6f},  세 자리로 찍으면 {se_tv:.3f}")
print(f"찍힌 수로 셈한 t = 0.0470/0.002 = {0.0470 / 0.002:.2f},  참 t = {beta_hat[1, 0] / se_tv:.4f}")

자유도 = 136,  t 임계값 = 1.977561
  Intercept  |t|>t* True   p<0.05 True   CI 가 0 을 배제 True   (t=+7.7819, p=1.609e-12)
  TV         |t|>t* True   p<0.05 True   CI 가 0 을 배제 True   (t=+27.6535, p=1.092e-57)
  Radio      |t|>t* True   p<0.05 True   CI 가 0 을 배제 True   (t=+16.6646, p=1.161e-34)
  Newspaper  |t|>t* False  p<0.05 False  CI 가 0 을 배제 False  (t=-0.4285, p=6.690e-01)
TV 의 표준오차: 참값 0.001701,  세 자리로 찍으면 0.002
찍힌 수로 셈한 t = 0.0470/0.002 = 23.50,  참 t = 27.6535

import pandas as pd

url = ('https://raw.githubusercontent.com/justmarkham/'
       'scikit-learn-videos/8545c74961398def7724501648fd504dbf061b41/data/Advertising.csv')
data = pd.read_csv(url, usecols=[1, 2, 3, 4])
# 앞 70%만 훈련에 쓴다.
training_data = data.iloc[:int(len(data) * 0.7)]

y = np.array(training_data.Sales).reshape(-1, 1)
n = y.shape[0]
# 1 로 채운 열을 앞에 붙여 절편을 만든다. 마지막 열이 반응(Sales)이므로
# iloc[:, :-1] 로 설명변수만 고른다.
X = np.concatenate(
    (np.ones((n, 1)), np.array(training_data.iloc[:, :-1])), axis=1
)
k = X.shape[1]

beta_hat, s, cov_matrix = fit_ols(X, y)
var_names = ["Intercept", "TV", "Radio", "Newspaper"]
regression_table(beta_hat, s, cov_matrix, n, k, var_names)

Intercept   coef=3.0451  SE=0.391  t=7.782  p=0.000  CI=(2.271, 3.819)
TV          coef=0.0470  SE=0.002  t=27.653  p=0.000  CI=(0.044, 0.050)
Radio       coef=0.1797  SE=0.011  t=16.665  p=0.000  CI=(0.158, 0.201)
Newspaper   coef=-0.0030  SE=0.007  t=-0.428  p=0.669  CI=(-0.017, 0.011)

Intercept   coef=3.0451  SE=0.391  t=7.782   p=0.000  CI=(2.271, 3.819)
TV          coef=0.0470  SE=0.002  t=27.653  p=0.000  CI=(0.044, 0.050)
Radio       coef=0.1797  SE=0.011  t=16.665  p=0.000  CI=(0.158, 0.201)
Newspaper   coef=-0.0030 SE=0.007  t=-0.428  p=0.669  CI=(-0.017, 0.011)

import statsmodels.api as sm

ref = sm.OLS(y, X).fit()
se_hand = s * np.sqrt(np.diag(cov_matrix))
t_hand = beta_hat.ravel() / se_hand
p_hand = 2 * stats.t(n - k).sf(np.abs(t_hand))
print(f"계수   최대 차이 = {np.abs(ref.params - beta_hat.ravel()).max():.3e}")
print(f"표준오차 최대 차이 = {np.abs(ref.bse - se_hand).max():.3e}")
print(f"t      최대 차이 = {np.abs(ref.tvalues - t_hand).max():.3e}")
print(f"p-값   최대 차이 = {np.abs(ref.pvalues - p_hand).max():.3e}")
print(f"X'X 의 조건수 = {np.linalg.cond(X.T @ X):.4e}")

def coef_and_se(d):
    yy = np.array(d.Sales).reshape(-1, 1)
    nn = len(yy)
    XX = np.concatenate((np.ones((nn, 1)), np.array(d.iloc[:, :-1])), axis=1)
    bb, ss, VV = fit_ols(XX, yy)
    return bb.ravel(), ss * np.sqrt(np.diag(VV))

b_tr, se_tr = coef_and_se(training_data)
b_all, se_all = coef_and_se(data)
rng = np.random.RandomState(0)
b_rnd, se_rnd = coef_and_se(data.iloc[rng.permutation(len(data))[:140]])
print("앞 70%    :", b_tr.round(4), " SE", se_tr.round(4))
print("전체 200행 :", b_all.round(4), " SE", se_all.round(4))
print("무작위 70% :", b_rnd.round(4), " SE", se_rnd.round(4))
print(f"Newspaper 의 t: 앞70% {b_tr[3] / se_tr[3]:.3f},  "
      f"전체 {b_all[3] / se_all[3]:.3f},  무작위70% {b_rnd[3] / se_rnd[3]:.3f}")

계수   최대 차이 = 4.396e-14
표준오차 최대 차이 = 6.661e-16
t      최대 차이 = 1.252e-13
p-값   최대 차이 = 1.410e-14
X'X 의 조건수 = 2.0924e+05
앞 70%    : [ 3.0451e+00  4.7000e-02  1.7970e-01 -3.0000e-03]  SE [0.3913 0.0017 0.0108 0.007 ]
전체 200행 : [ 2.9389e+00  4.5800e-02  1.8850e-01 -1.0000e-03]  SE [0.3119 0.0014 0.0086 0.0059]
무작위 70% : [ 2.9924e+00  4.7700e-02  1.7650e-01 -2.7000e-03]  SE [0.3595 0.0017 0.0107 0.0068]
Newspaper 의 t: 앞70% -0.428,  전체 -0.177,  무작위70% -0.399

lhs = X.T @ X @ beta_hat
rhs = X.T @ y
print(np.allclose(lhs, rhs))  # True

True

y_hat = X @ beta_hat
RSS = np.sum((y - y_hat) ** 2)
TSS = np.sum((y - y.mean()) ** 2)
R2 = 1 - RSS / TSS
adj_R2 = 1 - (1 - R2) * (n - 1) / (n - k)
print(f"R^2 = {R2:.4f}, Adjusted R^2 = {adj_R2:.4f}")

R^2 = 0.8937, Adjusted R^2 = 0.8914

beta_lstsq, residuals, rank, sv = np.linalg.lstsq(X, y, rcond=None)