부분집합 선택과 단계적 선택¶

import numpy as np
from sklearn.linear_model import LinearRegression

# 앞의 넷만 참 계수가 0 이 아니다. 세 방법이 이 넷을 찾아내는지 견준다.
np.random.seed(42)
n, p = 200, 8
X = np.random.randn(n, p)
true_beta = np.array([3.0, 1.5, -2.0, 0.8, 0, 0, 0, 0])
y = X @ true_beta + np.random.normal(0, 2, n)
names = [f"x{i+1}" for i in range(p)]

import statsmodels.api as sm

print(f"sigma/sqrt(n) = {2 / np.sqrt(n):.4f}")
fit = sm.OLS(y, sm.add_constant(X)).fit()
print("실제 SE  :", fit.bse[1:].round(4))
print("기대 t   :", (true_beta / (2 / np.sqrt(n))).round(2))
print("실제 t   :", fit.tvalues[1:].round(2))
print("추정 계수:", fit.params[1:].round(3))
V_signal = (true_beta ** 2).sum()
print(f"모집단 R^2 = {V_signal} / ({V_signal} + 4) = {V_signal / (V_signal + 4):.4f},  실제 R^2 = {fit.rsquared:.4f}")
print(f"그림의 어려운 설정(n=60, beta4=0.6)에서 기대 t = {0.6 / (2 / np.sqrt(60)):.2f}")

sigma/sqrt(n) = 0.1414
실제 SE  : [0.1547 0.1426 0.1396 0.1381 0.1449 0.1272 0.1364 0.135 ]
기대 t   : [ 21.21  10.61 -14.14   5.66   0.     0.     0.     0.  ]
실제 t   : [ 18.64  10.39 -13.14   6.02  -1.15  -0.15  -0.07   0.74]
추정 계수: [ 2.884  1.481 -1.834  0.832 -0.166 -0.019 -0.01   0.1  ]
모집단 R^2 = 15.89 / (15.89 + 4) = 0.7989,  실제 R^2 = 0.7790
그림의 어려운 설정(n=60, beta4=0.6)에서 기대 t = 2.32

from itertools import combinations

def best_subset(X, y, max_k=None):
    """모든 부분집합을 다 따져 크기별 최선을 찾는다.

    크기 k 마다 가능한 조합을 남김없이 본다. 답은 확실하지만 부분집합이
    2^p 개라 변수가 스물만 넘어도 감당할 수 없다.
    """
    n, p = X.shape
    if max_k is None:
        max_k = p
    results = {}
    for k in range(1, max_k + 1):
        best_rss, best_features = np.inf, None
        for combo in combinations(range(p), k):
            model = LinearRegression().fit(X[:, combo], y)
            rss = np.sum((y - model.predict(X[:, combo])) ** 2)
            if rss < best_rss:
                best_rss, best_features = rss, combo
        results[k] = {"features": best_features, "rss": best_rss}
    return results

from math import comb

print("크기별 후보 개수:", [comb(p, k) for k in range(1, p + 1)])
print(f"합 = {sum(comb(p, k) for k in range(1, p + 1))},  2^p - 1 = {2 ** p - 1}")
print(f"p = 20 이면 {2 ** 20 - 1:,} 개")

best = best_subset(X, y)
print("  k  변수집합                     RSS")
for k in range(1, p + 1):
    print(f"  {k}  {str(best[k]['features']):28s} {best[k]['rss']:9.3f}")
print("RSS 가 k 에 대해 단조감소:",
      all(best[k + 1]['rss'] <= best[k]['rss'] + 1e-9 for k in range(1, p)))

크기별 후보 개수: [8, 28, 56, 70, 56, 28, 8, 1]
합 = 255,  2^p - 1 = 255
p = 20 이면 1,048,575 개
  k  변수집합                     RSS
  1  (0,)                          1972.005
  2  (0, 2)                        1303.749
  3  (0, 1, 2)                      861.547
  4  (0, 1, 2, 3)                   711.610
  5  (0, 1, 2, 3, 4)                707.136
  6  (0, 1, 2, 3, 4, 7)             705.141
  7  (0, 1, 2, 3, 4, 5, 7)          705.057
  8  (0, 1, 2, 3, 4, 5, 6, 7)       705.038
RSS 가 k 에 대해 단조감소: True

def forward_stepwise(X, y):
    """빈 모형에서 시작해 RSS 를 가장 많이 줄이는 변수를 하나씩 더한다.

    따지는 모형이 p(p+1)/2 개로 줄어 훨씬 빠르다. 다만 한 번 들어간 변수는
    빠지지 않으므로 최적 부분집합을 놓칠 수 있다.
    """
    n, p = X.shape
    selected, remaining = [], list(range(p))
    results = {}
    for k in range(1, p + 1):
        best_rss, best_feature = np.inf, None
        for f in remaining:
            trial = selected + [f]
            model = LinearRegression().fit(X[:, trial], y)
            rss = np.sum((y - model.predict(X[:, trial])) ** 2)
            if rss < best_rss:
                best_rss, best_feature = rss, f
        selected.append(best_feature)
        remaining.remove(best_feature)
        results[k] = {"features": tuple(selected), "rss": best_rss}
    return results

print(f"forward 가 따지는 모형 수 = p + (p-1) + ... + 1 = {p * (p + 1) // 2},  "
      f"best_subset 은 {2 ** p - 1} (비 {(2 ** p - 1) / (p * (p + 1) // 2):.1f}배)")
print(f"p = 20 이면 {20 * 21 // 2} 대 {2 ** 20 - 1:,}")

fwd = forward_stepwise(X, y)
print("들어온 순서:", fwd[p]['features'])
print("  k  변수집합(정렬)              fwd RSS     best RSS   차이")
for k in range(1, p + 1):
    print(f"  {k}  {str(tuple(sorted(fwd[k]['features']))):28s} {fwd[k]['rss']:9.3f}  "
          f"{best[k]['rss']:9.3f}  {fwd[k]['rss'] - best[k]['rss']:.2e}")

forward 가 따지는 모형 수 = p + (p-1) + ... + 1 = 36,  best_subset 은 255 (비 7.1배)
p = 20 이면 210 대 1,048,575
들어온 순서: (0, 2, 1, 3, 4, 7, 5, 6)
  k  변수집합(정렬)              fwd RSS     best RSS   차이
  1  (0,)                          1972.005   1972.005  0.00e+00
  2  (0, 2)                        1303.749   1303.749  0.00e+00
  3  (0, 1, 2)                      861.547    861.547  0.00e+00
  4  (0, 1, 2, 3)                   711.610    711.610  1.14e-13
  5  (0, 1, 2, 3, 4)                707.136    707.136  0.00e+00
  6  (0, 1, 2, 3, 4, 7)             705.141    705.141  1.14e-13
  7  (0, 1, 2, 3, 4, 5, 7)          705.057    705.057  -2.27e-13
  8  (0, 1, 2, 3, 4, 5, 6, 7)       705.038    705.038  0.00e+00

def backward_stepwise(X, y):
    """전체 모형에서 시작해 RSS 를 가장 적게 늘리는 변수를 하나씩 뺀다.

    시작점이 전체 모형이므로 n > p 여야 쓸 수 있다. 변수가 관측보다 많으면
    전진선택으로 가야 한다.
    """
    n, p = X.shape
    current = list(range(p))
    results = {}
    model = LinearRegression().fit(X, y)
    results[p] = {"features": tuple(current),
                  "rss": np.sum((y - model.predict(X)) ** 2)}
    for k in range(p - 1, 0, -1):
        best_rss, best_remove = np.inf, None
        for f in current:
            trial = [x for x in current if x != f]
            model = LinearRegression().fit(X[:, trial], y)
            rss = np.sum((y - model.predict(X[:, trial])) ** 2)
            if rss < best_rss:
                best_rss, best_remove = rss, f
        current.remove(best_remove)
        results[k] = {"features": tuple(current), "rss": best_rss}
    return results

print(f"backward 가 따지는 모형 수 = 1 + (p + (p-1) + ... + 2) = {1 + sum(range(2, p + 1))}")

bwd = backward_stepwise(X, y)
print("  k  변수집합                     RSS      fwd 와 같은 집합인가")
for k in range(p, 0, -1):
    print(f"  {k}  {str(bwd[k]['features']):28s} {bwd[k]['rss']:9.3f}   "
          f"{set(bwd[k]['features']) == set(fwd[k]['features'])}")
print("세 방법이 모든 k 에서 같은 집합:",
      all(set(best[k]['features']) == set(fwd[k]['features']) == set(bwd[k]['features'])
          for k in range(1, p + 1)))

# n <= p 면 출발점이 무너진다.
X_small, y_small = X[:6], y[:6]
print(f"n = {len(y_small)}, p = {p} 일 때 전체모형 설계행렬의 rank = "
      f"{np.linalg.matrix_rank(np.column_stack([np.ones(len(y_small)), X_small]))}  (열 {p + 1}개)")
full = LinearRegression().fit(X_small, y_small)
print(f"  그 모형의 훈련 RSS = {np.sum((y_small - full.predict(X_small)) ** 2):.3e}")

backward 가 따지는 모형 수 = 1 + (p + (p-1) + ... + 2) = 36
  k  변수집합                     RSS      fwd 와 같은 집합인가
  8  (0, 1, 2, 3, 4, 5, 6, 7)       705.038   True
  7  (0, 1, 2, 3, 4, 5, 7)          705.057   True
  6  (0, 1, 2, 3, 4, 7)             705.141   True
  5  (0, 1, 2, 3, 4)                707.136   True
  4  (0, 1, 2, 3)                   711.610   True
  3  (0, 1, 2)                      861.547   True
  2  (0, 2)                        1303.749   True
  1  (0,)                          1972.005   True
세 방법이 모든 k 에서 같은 집합: True
n = 6, p = 8 일 때 전체모형 설계행렬의 rank = 6  (열 9개)
  그 모형의 훈련 RSS = 6.306e-29

n_train = 140
X_tr, X_val = X[:n_train], X[n_train:]
y_tr, y_val = y[:n_train], y[n_train:]

best = best_subset(X_tr, y_tr)
fwd = forward_stepwise(X_tr, y_tr)
bwd = backward_stepwise(X_tr, y_tr)

# 방법마다 검증 RSS 가 가장 작은 k 를 찾는다
for method_name, res in [("Best", best), ("Fwd", fwd), ("Bwd", bwd)]:
    val_rss = [np.sum((y_val - LinearRegression().fit(X_tr[:, res[k]["features"]],
               y_tr).predict(X_val[:, res[k]["features"]])) ** 2)
               for k in range(1, 9)]
    opt_k = np.argmin(val_rss) + 1
    print(f"{method_name}: k={opt_k}, features={res[opt_k]['features']}")

Best: k=4, features=(0, 1, 2, 3)
Fwd: k=4, features=(0, 2, 1, 3)
Bwd: k=4, features=(0, 1, 2, 3)

def aic(n, rss, k):
    return n * np.log(rss / n) + 2 * (k + 1)  # +1 for intercept

fwd = forward_stepwise(X_tr, y_tr)
aic_vals = []
for k in range(1, 9):
    aic_vals.append(aic(n_train, fwd[k]['rss'], k))
opt_k_aic = np.argmin(aic_vals) + 1
print("AIC per k: [" + ", ".join(f"{v:.2f}" for v in aic_vals) + "]")
print(f"Optimal k (AIC): {opt_k_aic}")

AIC per k: [310.90, 267.13, 204.61, 178.62, 178.86, 179.83, 181.65, 183.59]
Optimal k (AIC): 4