모형선택 비교¶

import numpy as np

def aic(n, rss, k):
    """AIC. 모수 하나당 벌점이 2 다."""
    return n * np.log(rss / n) + 2 * k

def bic(n, rss, k):
    """BIC. 벌점이 log(n) 이라 n>7 부터 AIC 보다 무겁다.

    그래서 BIC 는 대개 더 작은 모형을 고른다. AIC 는 예측을 잘하는 모형을,
    BIC 는 참 모형을 찾는 것을 겨냥한다고 흔히 말한다.
    """
    return n * np.log(rss / n) + k * np.log(n)

import statsmodels.api as sm

rng = np.random.RandomState(0)          # 전역 난수 상태를 건드리지 않는다
n_demo = 200
X_demo = np.column_stack([np.ones(n_demo), rng.randn(n_demo, 3)])
y_demo = X_demo @ np.array([1.0, 2.0, -1.0, 0.5]) + rng.randn(n_demo)
res = sm.OLS(y_demo, X_demo).fit()
rss_demo = (res.resid ** 2).sum()
k_demo = X_demo.shape[1]

print(f"쪽의 aic()      = {aic(n_demo, rss_demo, k_demo):.4f}")
print(f"statsmodels aic = {res.aic:.4f}")
print(f"차이            = {res.aic - aic(n_demo, rss_demo, k_demo):.4f}")
print(f"n(1 + ln 2pi)   = {n_demo * (1 + np.log(2 * np.pi)):.4f}")
print(f"bic 쪽의 차이   = {res.bic - bic(n_demo, rss_demo, k_demo):.4f}")
print(f"ln n = {np.log(n_demo):.4f},  2 와의 비 = {np.log(n_demo) / 2:.4f},  e^2 = {np.e ** 2:.4f}")

쪽의 aic()      = -8.9125
statsmodels aic = 558.6629
차이            = 567.5754
n(1 + ln 2pi)   = 567.5754
bic 쪽의 차이   = 567.5754
ln n = 5.2983,  2 와의 비 = 2.6492,  e^2 = 7.3891

def cv_mse(X, y, folds=5):
    """k겹 교차검증으로 예측오차를 추정한다.

    AIC·BIC 가 공식으로 벌점을 매긴다면, 교차검증은 실제로 떼어 놓은
    자료에서 재 본다. 가정이 적은 대신 계산이 많이 든다.
    """
    n = len(y)
    indices = np.arange(n)
    np.random.shuffle(indices)
    fold_size = n // folds
    mses = []
    for k in range(folds):
        val_idx = indices[k * fold_size:(k + 1) * fold_size]
        train_idx = np.setdiff1d(indices, val_idx)
        X_tr, y_tr = X[train_idx], y[train_idx]
        X_va, y_va = X[val_idx], y[val_idx]
        beta = np.linalg.lstsq(X_tr, y_tr, rcond=None)[0]
        pred = X_va @ beta
        mses.append(np.mean((y_va - pred) ** 2))
    return np.mean(mses)

n_demo = 200
for f in [5, 6, 7, 9, 30]:
    used = f * (n_demo // f)
    print(f"  folds={f:2d}: fold_size={n_demo // f:3d},  검증에 쓰이는 관측값 {used}개,  한 번도 안 쓰이는 것 {n_demo - used}개")

# 아래 보기 3 과 같은 자료를 RandomState 로 만든다(전역 상태를 건드리지 않는다).
rng2 = np.random.RandomState(42)
p_demo = 8
Xd = rng2.randn(n_demo, p_demo)
yd = Xd @ np.array([3.0, -2.0, 1.5, 0, 0, 0, 0, 0]) + rng2.randn(n_demo) * 2
order = [0, 1, 2, 4, 7, 3, 5, 6]          # 보기 3 이 고르게 되는 순서
designs = [np.column_stack([np.ones(n_demo), Xd[:, order[:s]]]) for s in range(1, 9)]

np.random.seed(0)
col3 = np.array([cv_mse(designs[2], yd) for _ in range(20)])
print(f"크기 3 모형에 cv_mse 를 20번: 평균 {col3.mean():.4f}  표준편차 {col3.std(ddof=1):.4f}  "
      f"최소 {col3.min():.4f}  최대 {col3.max():.4f}")

# 여덟 크기의 CV 열을 500 번 다시 만들어, 최솟값이 어디에 떨어지는지 센다.
import collections
np.random.seed(0)
cnt = collections.Counter(); diff = []
for _ in range(500):
    col = [cv_mse(D, yd) for D in designs]
    cnt[int(np.argmin(col)) + 1] += 1
    diff.append(col[2] - col[3])
diff = np.array(diff)
print("최솟값이 나온 모형 크기의 빈도(500회):", dict(sorted(cnt.items())))
print(f"크기3 - 크기4: 평균 {diff.mean():+.4f}  표준편차 {diff.std(ddof=1):.4f}  "
      f"(평균의 몬테카를로 오차 {diff.std(ddof=1) / np.sqrt(500):.4f})")
print(f"크기 3 이 크기 4 보다 작게 나온 비율 = {(diff < 0).mean():.3f}")

  folds= 5: fold_size= 40,  검증에 쓰이는 관측값 200개,  한 번도 안 쓰이는 것 0개
  folds= 6: fold_size= 33,  검증에 쓰이는 관측값 198개,  한 번도 안 쓰이는 것 2개
  folds= 7: fold_size= 28,  검증에 쓰이는 관측값 196개,  한 번도 안 쓰이는 것 4개
  folds= 9: fold_size= 22,  검증에 쓰이는 관측값 198개,  한 번도 안 쓰이는 것 2개
  folds=30: fold_size=  6,  검증에 쓰이는 관측값 180개,  한 번도 안 쓰이는 것 20개
크기 3 모형에 cv_mse 를 20번: 평균 3.7151  표준편차 0.0527  최소 3.6409  최대 3.8332
최솟값이 나온 모형 크기의 빈도(500회): {3: 217, 4: 154, 5: 85, 6: 35, 7: 7, 8: 2}
크기3 - 크기4: 평균 -0.0164  표준편차 0.0771  (평균의 몬테카를로 오차 0.0034)
크기 3 이 크기 4 보다 작게 나온 비율 = 0.570

# 여덟 변수 중 앞의 셋만 실제로 쓰이고 나머지 다섯은 계수가 0 이다.
# 전진선택이 그 셋을 먼저 집어내는지 보는 것이 이 실험의 목적이다.
np.random.seed(42)
n, p_total = 200, 8
X_raw = np.random.randn(n, p_total)
beta_true = np.array([3.0, -2.0, 1.5, 0, 0, 0, 0, 0])
y = X_raw @ beta_true + np.random.randn(n) * 2

remaining = list(range(p_total))
selected = []
aic_history, bic_history = [], []

# 전진선택: 매 단계에서 AIC 를 가장 많이 낮추는 변수를 하나씩 더한다.
# 모든 부분집합을 따지면 2^8 개지만, 이렇게 하면 훨씬 적게 본다.
# 대신 최적 부분집합을 놓칠 수 있다.
for step in range(p_total):
    best_score, best_j = np.inf, None
    for j in remaining:
        cols = selected + [j]
        X_cand = np.column_stack([np.ones(n), X_raw[:, cols]])
        beta = np.linalg.lstsq(X_cand, y, rcond=None)[0]
        rss = np.sum((y - X_cand @ beta) ** 2)
        score = aic(n, rss, len(cols) + 1)
        if score < best_score:
            best_score, best_j = score, j
    selected.append(best_j)
    remaining.remove(best_j)

    X_sel = np.column_stack([np.ones(n), X_raw[:, selected]])
    beta = np.linalg.lstsq(X_sel, y, rcond=None)[0]
    rss = np.sum((y - X_sel @ beta) ** 2)
    k = len(selected) + 1
    aic_history.append(aic(n, rss, k))
    bic_history.append(bic(n, rss, k))

from scipy import stats

def rss_of(s):
    Xs = np.column_stack([np.ones(n), X_raw[:, selected[:s]]])
    return np.sum((y - Xs @ np.linalg.lstsq(Xs, y, rcond=None)[0]) ** 2)

print("선택 순서:", selected)
print("크기   RSS        AIC       BIC")
for s in range(1, p_total + 1):
    print(f" {s}   {rss_of(s):9.4f}  {aic_history[s - 1]:8.2f}  {bic_history[s - 1]:8.2f}")

rss3, rss4 = rss_of(3), rss_of(4)
print(f"n ln(RSS4/RSS3) = {n * np.log(rss4 / rss3):.4f}")
print(f"AIC 차 = 그것 + 2    = {n * np.log(rss4 / rss3) + 2:.4f}   (표의 차 {aic_history[3] - aic_history[2]:.4f})")
print(f"BIC 차 = 그것 + ln n = {n * np.log(rss4 / rss3) + np.log(n):.4f}   (표의 차 {bic_history[3] - bic_history[2]:.4f})")
print(f"네 번째 변수의 부분 F = {(rss3 - rss4) / (rss4 / (n - 5)):.4f}")
print(f"P(chi2_1 > 2) = {stats.chi2.sf(2, 1):.4f},  P(chi2_1 > ln 200) = {stats.chi2.sf(np.log(n), 1):.4f}")
print(f"RSS3/n = {rss3 / n:.4f},  RSS3/(n-4) = {rss3 / (n - 4):.4f}  (참 sigma^2 = 4)")

선택 순서: [0, 1, 2, 4, 7, 3, 5, 6]
크기   RSS        AIC       BIC
 1   2002.7003    464.79    471.38
 2   1243.9264    371.54    381.44
 3    712.3828    262.06    275.25
 4    707.4385    262.67    279.16
 5    705.3433    264.07    283.86
 6    705.1412    266.02    289.10
 7    705.0570    267.99    294.38
 8    705.0384    269.99    299.67
n ln(RSS4/RSS3) = -1.3929
AIC 차 = 그것 + 2    = 0.6071   (표의 차 0.6071)
BIC 차 = 그것 + ln n = 3.9054   (표의 차 3.9054)
네 번째 변수의 부분 F = 1.3628
P(chi2_1 > 2) = 0.1573,  P(chi2_1 > ln 200) = 0.0213
RSS3/n = 3.5619,  RSS3/(n-4) = 3.6346  (참 sigma^2 = 4)

cv5 = [cv_mse(np.column_stack([np.ones(n), X_raw[:, selected[:s]]]), y, folds=5)
       for s in range(1, p_total + 1)]
cv10 = [cv_mse(np.column_stack([np.ones(n), X_raw[:, selected[:s]]]), y, folds=10)
        for s in range(1, p_total + 1)]