교차검증 방법 비교 (코드)¶

import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline

def generate_data(n=200, seed=42):
    """sin 곡선에 선형 추세와 잡음을 얹은 자료를 만든다.

    참 관계가 다항식이 아니므로 "옳은 차수"가 따로 없다. 세 교차검증
    방법이 어느 차수를 고르는지 견주는 것이 목적이다.
    """
    rng = np.random.default_rng(seed)
    x = rng.uniform(-3, 3, n)
    y = np.sin(x) + 0.3 * x + rng.normal(0, 0.5, n)
    return x.reshape(-1, 1), y

def poly_pipeline(degree):
    """주어진 차수의 다항회귀 파이프라인."""
    return Pipeline([
        ("poly", PolynomialFeatures(degree=degree, include_bias=False)),
        ("lr", LinearRegression()),
    ])

X, y = generate_data(n=200)
degrees = range(1, 11)

from scipy import integrate

VarX = 3.0
Exsin = integrate.quad(lambda x: x * np.sin(x) / 6, -3, 3)[0]
Esin2 = integrate.quad(lambda x: np.sin(x) ** 2 / 6, -3, 3)[0]
Varg = Esin2 + 0.09 * VarX + 0.6 * Exsin
Cov = Exsin + 0.3 * VarX
print(f"E[X sin X] = {Exsin:.6f},  E[sin^2 X] = {Esin2:.6f}")
print(f"Var(g) = {Varg:.6f},  Cov(g, X) = {Cov:.6f},  최적 기울기 = {Cov/VarX:.6f}")
print(f"1차의 근사오차 = {Varg - Cov**2/VarX:.6f}")
print(f"바닥 sigma^2 = {0.25:.6f},  1차의 모집단 MSE = {0.25 + Varg - Cov**2/VarX:.6f}")

E[X sin X] = 1.037032,  E[sin^2 X] = 0.523285
Var(g) = 1.415504,  Cov(g, X) = 1.937032,  최적 기울기 = 0.645677
1차의 근사오차 = 0.164806
바닥 sigma^2 = 0.250000,  1차의 모집단 MSE = 0.414806

def validation_set_mse(X, y, degrees, n_splits=10, rng=None):
    """검증집합 방법. 자료를 절반으로 갈라 한쪽으로 훈련하고 한쪽으로 잰다.

    가장 간단하지만 두 가지 약점이 있다. 어떻게 가르느냐에 따라 결과가
    들쭉날쭉하고(그래서 여기서는 10번 다르게 갈라 본다), 훈련에 절반밖에
    쓰지 못해 성능을 실제보다 나쁘게 잡는다.
    """
    rng = rng or np.random.default_rng(7)
    n = len(y); n_train = n // 2
    all_mses = {d: [] for d in degrees}
    for _ in range(n_splits):
        perm = rng.permutation(n)
        tr, te = perm[:n_train], perm[n_train:]
        for d in degrees:
            model = poly_pipeline(d).fit(X[tr], y[tr])
            all_mses[d].append(np.mean((y[te] - model.predict(X[te])) ** 2))
    return all_mses

am = validation_set_mse(X, y, degrees)
M = np.array([am[d] for d in degrees])
best = M.argmin(axis=0) + 1
print(f"각 분할이 고른 차수 = {list(best)},  표준편차 = {best.std():.4f}")
print(f"차수 6: 열 번의 평균 {M[5].mean():.6f},  분할에 따른 SD {M[5].std(ddof=1):.6f}")

각 분할이 고른 차수 = [10, 6, 9, 5, 5, 6, 6, 3, 3, 9],  표준편차 = 2.3152
차수 6: 열 번의 평균 0.277278,  분할에 따른 SD 0.035057

from sklearn.model_selection import cross_val_score, LeaveOneOut
import time

# 하나빼기 교차검증(LOOCV). 관측값 하나씩만 남기므로 훈련자료를 거의
# 다 쓰고, 어떻게 가르느냐에 따른 흔들림도 없다. 대신 n 번 적합해야 한다.
t0 = time.time()
loocv_mses = [-cross_val_score(poly_pipeline(d), X, y, cv=LeaveOneOut(),
                               scoring="neg_mean_squared_error").mean()
              for d in degrees]
elapsed = time.time() - t0        # 아래 주의 참조: 출력에는 싣지 않는다

print(f"Best degree (LOOCV): {int(np.argmin(loocv_mses)) + 1}")

Best degree (LOOCV): 6

from sklearn.preprocessing import PolynomialFeatures

def loocv_hat(X, y, d):
    """모자행렬로 한 번에 구하는 LOOCV."""
    Z = PolynomialFeatures(degree=d, include_bias=True).fit_transform(X)
    H = Z @ np.linalg.pinv(Z.T @ Z) @ Z.T
    e = y - H @ y
    return np.mean((e / (1 - np.diag(H))) ** 2)

hat = [loocv_hat(X, y, d) for d in degrees]
print(f"{'d':>3}{'sklearn':>12}{'모자행렬':>12}{'차이':>12}")
for d, a, b in zip(degrees, loocv_mses, hat):
    print(f"{d:>3}{a:>12.6f}{b:>12.6f}{a-b:>12.2e}")
print(f"적합 횟수: sklearn {200*10} 번,  모자행렬 {10} 번")

  d     sklearn        모자행렬          차이
  1    0.442177    0.442177    1.11e-16
  2    0.447928    0.447928    0.00e+00
  3    0.263086    0.263086   -5.55e-17
  4    0.262563    0.262563    3.33e-16
  5    0.259166    0.259166   -5.83e-15
  6    0.257722    0.257722   -2.22e-16
  7    0.261009    0.261009    1.14e-12
  8    0.264039    0.264039    4.01e-14
  9    0.262481    0.262481    1.95e-11
 10    0.267918    0.267918    5.43e-11
적합 횟수: sklearn 2000 번,  모자행렬 10 번

from sklearn.model_selection import KFold

# k겹 교차검증. 앞 두 방법의 절충이다. k 번만 적합하면 되고, 겹마다의
# 훈련자료가 조금씩 겹쳐 LOOCV 보다 추정의 분산이 오히려 작다.
# 실무에서 5 나 10 을 쓰는 까닭이 이것이다.
for k in (5, 10):
    kf = KFold(n_splits=k, shuffle=True, random_state=42)
    mses = [-cross_val_score(poly_pipeline(d), X, y, cv=kf,
                             scoring="neg_mean_squared_error").mean()
            for d in degrees]
    print(f"Best degree ({k}-fold): {int(np.argmin(mses)) + 1}")

Best degree (5-fold): 6
Best degree (10-fold): 6

# (가) 자료는 그대로 두고 겹 배정만 바꾼다.
for k in (5, 10):
    vals = np.array([
        -cross_val_score(poly_pipeline(6), X, y,
                         cv=KFold(n_splits=k, shuffle=True, random_state=s),
                         scoring="neg_mean_squared_error").mean()
        for s in range(20)])
    print(f"k={k:2d}: 겹 배정 20 번의 평균 {vals.mean():.6f},  SD {vals.std(ddof=1):.6f}")
print(f"LOOCV: {loocv_hat(X, y, 6):.6f},  SD 0 (결정적)")

# (나) 자료를 새로 뽑는다.
res = {5: [], 10: [], 'loo': []}
for s in range(100):
    rng = np.random.default_rng(1000 + s)
    xx = rng.uniform(-3, 3, 200)
    yy = np.sin(xx) + 0.3 * xx + rng.normal(0, 0.5, 200)
    XX = xx.reshape(-1, 1)
    for k in (5, 10):
        res[k].append(-cross_val_score(
            poly_pipeline(6), XX, yy,
            cv=KFold(n_splits=k, shuffle=True, random_state=s),
            scoring="neg_mean_squared_error").mean())
    res['loo'].append(loocv_hat(XX, yy, 6))
print(f"\n자료 100 개에서 (d=6)")
for k, lab in ((5, ' 5겹'), (10, '10겹'), ('loo', 'LOOCV')):
    v = np.array(res[k])
    print(f"  {lab}: 평균 {v.mean():.6f},  자료에 따른 SD {v.std(ddof=1):.6f}"
          f"  (SD 추정의 오차 {v.std(ddof=1)/np.sqrt(2*100):.6f})")

k= 5: 겹 배정 20 번의 평균 0.260405,  SD 0.007121
k=10: 겹 배정 20 번의 평균 0.258382,  SD 0.003266
LOOCV: 0.257722,  SD 0 (결정적)

자료 100 개에서 (d=6)
   5겹: 평균 0.265180,  자료에 따른 SD 0.026311  (SD 추정의 오차 0.001860)
  10겹: 평균 0.262959,  자료에 따른 SD 0.026537  (SD 추정의 오차 0.001876)
  LOOCV: 평균 0.261968,  자료에 따른 SD 0.025890  (SD 추정의 오차 0.001831)

import numpy as np
from collections import Counter
from sklearn.model_selection import cross_val_score, LeaveOneOut, KFold

degrees = list(range(1, 11))
res = {"validation": [], "5-fold": [], "10-fold": [], "LOOCV": []}

for trial in range(100):
    rng = np.random.default_rng(1000 + trial)
    x = rng.uniform(-3, 3, 200).reshape(-1, 1)
    y = np.sin(x.ravel()) + 0.3 * x.ravel() + rng.normal(0, 0.5, 200)

    perm = rng.permutation(200); tr, te = perm[:100], perm[100:]
    vm = [np.mean((y[te] - poly_pipeline(d).fit(x[tr], y[tr]).predict(x[te]))**2)
          for d in degrees]
    res["validation"].append(degrees[int(np.argmin(vm))])

    for name, cv in [("5-fold",  KFold(5,  shuffle=True, random_state=0)),
                     ("10-fold", KFold(10, shuffle=True, random_state=0)),
                     ("LOOCV",   LeaveOneOut())]:
        m = [-cross_val_score(poly_pipeline(d), x, y, cv=cv,
                              scoring="neg_mean_squared_error").mean()
             for d in degrees]
        res[name].append(degrees[int(np.argmin(m))])

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import cross_val_score, KFold

X, y = generate_data(n=200, seed=42)
degrees = list(range(1, 11))
all_mses = {d: [] for d in degrees}

for rep in range(5):
    kf = KFold(10, shuffle=True, random_state=rep)
    for d in degrees:
        all_mses[d].append(
            -cross_val_score(poly_pipeline(d), X, y, cv=kf,
                             scoring="neg_mean_squared_error").mean())

means = [np.mean(all_mses[d]) for d in degrees]
stds  = [np.std(all_mses[d]) for d in degrees]
print("Best degree (repeated 10-fold):", degrees[int(np.argmin(means))])

plt.errorbar(degrees, means, yerr=stds, fmt='o-', capsize=4)
plt.xlabel("Polynomial Degree"); plt.ylabel("CV MSE")
plt.title("Repeated 10-fold CV")
plt.show()

Best degree (repeated 10-fold): 5