교차검증과 람다 조율¶

import numpy as np
from sklearn.linear_model import Lasso
from sklearn.model_selection import KFold

rng = np.random.default_rng(42)

# 참으로 쓰이는 변수는 앞의 셋뿐이고 나머지 열일곱은 잡음이다.
n, p = 100, 20
X = rng.normal(size=(n, p))
beta_true = np.zeros(p)
beta_true[:3] = [3.0, -2.0, 1.5]
y = X @ beta_true + rng.normal(0, 1, n)

# lambda 격자는 로그 눈금으로 잡는다. 벌점의 효과가 곱셈으로 작동하므로
# 등간격보다 등비간격이 알맞다.
lambdas = np.logspace(0.5, -2.5, 25)

kf = KFold(n_splits=5, shuffle=True, random_state=0)

# 겹마다의 MSE 를 따로 남긴다. 평균만 구하면 1-표준오차 규칙을 쓸 수 없다.
fold_mse = np.zeros((5, len(lambdas)))
for k, (tr, va) in enumerate(kf.split(X)):
    for i, lam in enumerate(lambdas):
        model = Lasso(alpha=lam, max_iter=10000).fit(X[tr], y[tr])
        fold_mse[k, i] = np.mean((y[va] - model.predict(X[va])) ** 2)

cv_mean = fold_mse.mean(axis=0)
cv_se = fold_mse.std(axis=0, ddof=1) / np.sqrt(5)

i_min = int(np.argmin(cv_mean))

# 1-표준오차 규칙: 최솟값에서 1 표준오차 안에 드는 lambda 중 가장 큰 것을
# 고른다. CV 곡선의 최소점은 그 자체가 흔들리는 추정값이므로, 조금 더
# 단순한 모형 쪽으로 물러서는 편이 안전하다는 생각이다.
threshold = cv_mean[i_min] + cv_se[i_min]
i_1se = int(np.where(cv_mean <= threshold)[0][0])

for name, i in [("최소 CV", i_min), ("1-표준오차", i_1se)]:
    beta = Lasso(alpha=lambdas[i], max_iter=10000).fit(X, y).coef_
    print(f"{name:>10}: lambda = {lambdas[i]:.4f}, "
          f"CV MSE = {cv_mean[i]:.3f} (SE {cv_se[i]:.3f}), "
          f"0 이 아닌 계수 = {np.sum(np.abs(beta) > 1e-8)}개")
print("참으로 0 이 아닌 계수: 3개")

# --- (1) 의 바닥과 천장을 확인한다 ---
print(f"잡음분산 sigma^2 = 1 이므로 CV MSE 의 하한은 1.000")
print(f"최소 CV MSE = {cv_mean[i_min]:.4f}  ->  하한 위로 {cv_mean[i_min] - 1:+.4f}")
print(f"1-SE 문턱 = {cv_mean[i_min]:.4f} + {cv_se[i_min]:.4f} = {threshold:.4f}")
print(f"격자 첫 값 lambda = {lambdas[0]:.4f},  이론 lambda_max = |X'y|_inf/n = "
      f"{np.abs(X.T @ (y - y.mean())).max() / n:.4f}")
print(f"가장 큰 lambda 에서의 CV MSE = {cv_mean[0]:.4f}  (영모형: Var(y) = {y.var():.4f})")

     최소 CV: lambda = 0.1000, CV MSE = 1.001 (SE 0.044), 0 이 아닌 계수 = 9개
    1-표준오차: lambda = 0.1778, CV MSE = 1.038 (SE 0.083), 0 이 아닌 계수 = 5개
참으로 0 이 아닌 계수: 3개
잡음분산 sigma^2 = 1 이므로 CV MSE 의 하한은 1.000
최소 CV MSE = 1.0014  ->  하한 위로 +0.0014
1-SE 문턱 = 1.0014 + 0.0441 = 1.0455
격자 첫 값 lambda = 3.1623,  이론 lambda_max = |X'y|_inf/n = 3.1113
가장 큰 lambda 에서의 CV MSE = 17.4719  (영모형: Var(y) = 17.5023)

import numpy as np
from sklearn.model_selection import KFold

rng = np.random.default_rng(0)
X = rng.normal(5, 3, (100, 4))          # 평균과 척도가 0/1이 아닌 자료
y = X @ np.array([1.0, -1.0, 0.5, 0.0]) + rng.normal(0, 1, 100)

for k, (train_idx, val_idx) in enumerate(KFold(5).split(X)):
    X_train, X_val = X[train_idx], X[val_idx]
    mu = X_train.mean(axis=0)
    sigma = X_train.std(axis=0)
    X_train_s = (X_train - mu) / sigma
    X_val_s = (X_val - mu) / sigma
    # 훈련자료로 적합하고 검증자료에서 잰다
    if k == 0:
        print("훈련 겹 평균:", np.round(X_train_s.mean(axis=0), 4))
        print("검증 겹 평균:", np.round(X_val_s.mean(axis=0), 4))

훈련 겹 평균: [-0.  0. -0. -0.]
검증 겹 평균: [-0.1022  0.1113  0.5175  0.243 ]

def one_se_rule(lambdas, cv_mean, cv_se):
    """
    Select the largest lambda whose CV error is within
    one SE of the minimum CV error.
    """
    best_idx = np.argmin(cv_mean)
    threshold = cv_mean[best_idx] + cv_se[best_idx]

    # 교차검증 오차가 문턱 이하인 lambda 중 가장 큰 것을 고른다
    candidates = np.where(cv_mean <= threshold)[0]
    best_1se_idx = candidates[np.argmax(lambdas[candidates])]
    return lambdas[best_1se_idx]

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Lasso
from sklearn.model_selection import KFold
from sklearn.preprocessing import StandardScaler

np.random.seed(42)
n, p = 200, 30
X = np.random.randn(n, p)
beta_true = np.zeros(p)
beta_true[:5] = [3, -2, 4, -1, 2]
y = X @ beta_true + np.random.randn(n)

X_s = StandardScaler().fit_transform(X)
lambdas = np.logspace(1, -3, 50)
kf = KFold(n_splits=5, shuffle=True, random_state=42)

mse_folds = np.zeros((len(lambdas), 5))
for fold_idx, (tr, va) in enumerate(kf.split(X_s)):
    mu, sig = X_s[tr].mean(0), X_s[tr].std(0)
    Xtr = (X_s[tr] - mu) / sig
    Xva = (X_s[va] - mu) / sig
    for i, lam in enumerate(lambdas):
        m = Lasso(alpha=lam, max_iter=10000).fit(Xtr, y[tr])
        mse_folds[i, fold_idx] = np.mean((y[va] - m.predict(Xva))**2)

cv_mean = mse_folds.mean(axis=1)
cv_se = mse_folds.std(axis=1) / np.sqrt(5)

best_idx = np.argmin(cv_mean)
lam_min = lambdas[best_idx]
lam_1se = one_se_rule(lambdas, cv_mean, cv_se)

plt.errorbar(np.log10(lambdas), cv_mean, yerr=cv_se, fmt='o-', ms=4)
plt.axvline(np.log10(lam_min), color='red', ls='--', label='lambda_min')
plt.axvline(np.log10(lam_1se), color='blue', ls='--', label='lambda_1SE')
plt.xlabel('log10(lambda)')
plt.ylabel('CV MSE')
plt.legend()
plt.show()