능형회귀·라쏘·엘라스틱넷 비교¶

import numpy as np
from sklearn.preprocessing import StandardScaler

def generate_data(n=200, p=20, s=5, rho=0.8, noise=1.0):
    """서로 상관된 설명변수를 갖는 회귀자료를 만든다.

    이웃한 변수끼리 rho, 두 칸 떨어지면 rho^2 로 상관이 줄어드는 구조다.
    콜레스키 분해로 독립 정규에 이 상관을 입힌다. 세 벌점회귀가 갈리는
    곳이 바로 이런 상관 구조에서다.

    n: 관측 수, p: 변수 수, s: 참으로 0 이 아닌 계수의 수
    rho: 이웃한 변수 사이의 상관
    """
    Sigma = np.array([[rho**abs(i-j) for j in range(p)] for i in range(p)])
    L = np.linalg.cholesky(Sigma)
    X = np.random.randn(n, p) @ L.T

    beta_true = np.zeros(p)
    beta_true[:s] = np.array([3, -2, 1.5, -1, 0.5])

    y = X @ beta_true + noise * np.random.randn(n)
    return X, y, beta_true

np.random.seed(42)
X, y, beta_true = generate_data()
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# --- 유도한 값이 맞는지 확인한다 ---
rho, p = 0.8, 20
Sigma = np.array([[rho**abs(i-j) for j in range(p)] for i in range(p)])

R = np.corrcoef(X, rowvar=False)
for lag in (1, 2, 3):
    emp = np.mean([R[i, i+lag] for i in range(p - lag)])
    print(f"상관 lag {lag}:  이론 {rho**lag:.4f}   표본평균 {emp:.4f}")

sig2 = beta_true @ Sigma @ beta_true
print(f"신호분산 b'Sb = {sig2:.4f},   Var(y) 이론 = {sig2 + 1:.4f},  표본 = {y.var(ddof=1):.4f}")
print(f"이론 R^2 = {sig2 / (sig2 + 1):.4f}")

상관 lag 1:  이론 0.8000   표본평균 0.7943
상관 lag 2:  이론 0.6400   표본평균 0.6348
상관 lag 3:  이론 0.5120   표본평균 0.5093
신호분산 b'Sb = 5.3128,   Var(y) 이론 = 6.3128,  표본 = 6.1189
이론 R^2 = 0.8416

from sklearn.linear_model import RidgeCV, LassoCV, ElasticNetCV

# 세 방법 모두 교차검증으로 벌점을 고른다. 엘라스틱넷은 벌점의 세기와
# 배합비(l1_ratio) 둘을 함께 골라야 하므로 격자가 2차원이 된다.
alphas = np.logspace(-4, 2, 100)

ridge_cv = RidgeCV(alphas=alphas, cv=5)
ridge_cv.fit(X_scaled, y)

lasso_cv = LassoCV(n_alphas=100, cv=5, max_iter=10000)
lasso_cv.fit(X_scaled, y)

enet_cv = ElasticNetCV(
    l1_ratio=[0.1, 0.5, 0.7, 0.9, 0.95],
    n_alphas=100, cv=5, max_iter=10000
)
enet_cv.fit(X_scaled, y)

for name, m in [("Ridge", ridge_cv), ("Lasso", lasso_cv), ("ElasticNet", enet_cv)]:
    nz = np.sum(np.abs(m.coef_) > 1e-6)
    err = np.sum((m.coef_ - beta_true) ** 2)
    print(f"{name:11s} alpha={m.alpha_:.4f}  0이 아닌 계수 {nz:2d}개  "
          f"||b-b_true||^2 = {err:.4f}")
print(f"엘라스틱넷이 고른 l1_ratio = {enet_cv.l1_ratio_}")
print(f"능형의 alpha 를 라쏘 척도로 = {ridge_cv.alpha_ / (2 * len(y)):.6f}")

# 능형은 닫힌 꼴이 있다. 직접 풀어 sklearn 과 맞춰 본다.
beta_cf = np.linalg.solve(X_scaled.T @ X_scaled + ridge_cv.alpha_ * np.eye(20),
                          X_scaled.T @ (y - y.mean()))
print(f"닫힌 꼴과 sklearn 능형 계수의 최대 차이 = {np.abs(beta_cf - ridge_cv.coef_).max():.2e}")
d = np.linalg.svd(X_scaled, compute_uv=False)
print(f"유효자유도 df(lambda) = {np.sum(d**2 / (d**2 + ridge_cv.alpha_)):.4f}  (p = 20)")

Ridge       alpha=0.7565  0이 아닌 계수 20개  ||b-b_true||^2 = 0.2300
Lasso       alpha=0.0240  0이 아닌 계수  8개  ||b-b_true||^2 = 0.1178
ElasticNet  alpha=0.0253  0이 아닌 계수  8개  ||b-b_true||^2 = 0.1323
엘라스틱넷이 고른 l1_ratio = 0.95
능형의 alpha 를 라쏘 척도로 = 0.001891
닫힌 꼴과 sklearn 능형 계수의 최대 차이 = 2.66e-15
유효자유도 df(lambda) = 19.6479  (p = 20)

import matplotlib.pyplot as plt

# 참 계수와 세 방법의 계수를 나란히 그린다. 0 이 아닌 계수를 붉게 칠해
# 어느 방법이 몇 개를 살렸는지 한눈에 보이게 한다. 능형은 회색 막대가
# 하나도 없을 것이다 — 정확히 0 이 되지 못하기 때문이다.
fig, axes = plt.subplots(1, 4, figsize=(16, 4), sharey=True)
p = X_scaled.shape[1]

for ax, name, coefs in [
    (axes[0], "True", beta_true),
    (axes[1], "Ridge", ridge_cv.coef_),
    (axes[2], "Lasso", lasso_cv.coef_),
    (axes[3], "Elastic Net", enet_cv.coef_),
]:
    colors = ['#d32f2f' if abs(c) > 1e-6 else '#90a4ae' for c in coefs]
    ax.bar(range(p), coefs, color=colors, edgecolor='black', linewidth=0.3)
    ax.set_title(name)
    ax.set_xlabel("Feature index")
    ax.axhline(0, color='black', linewidth=0.5)

axes[0].set_ylabel("Coefficient value")
plt.tight_layout()
plt.show()

for name, coefs in [("True", beta_true), ("Ridge", ridge_cv.coef_),
                    ("Lasso", lasso_cv.coef_), ("ElasticNet", enet_cv.coef_)]:
    nz = np.where(np.abs(coefs) > 1e-6)[0]
    print(f"{name:11s} 0이 아닌 자리 {list(nz)}")
    print(f"{'':11s}  앞 다섯 계수 {np.round(coefs[:5], 3)}")

True        0이 아닌 자리 [0, 1, 2, 3, 4]
             앞 다섯 계수 [ 3.  -2.   1.5 -1.   0.5]
Ridge       0이 아닌 자리 [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]
             앞 다섯 계수 [ 2.791 -1.942  1.628 -1.21   0.705]
Lasso       0이 아닌 자리 [0, 1, 2, 3, 4, 6, 13, 16]
             앞 다섯 계수 [ 2.737 -1.849  1.486 -1.038  0.626]
ElasticNet  0이 아닌 자리 [0, 1, 2, 3, 4, 6, 13, 16]
             앞 다섯 계수 [ 2.722 -1.827  1.468 -1.026  0.621]

from sklearn.linear_model import Ridge, Lasso

alphas_path = np.logspace(-3, 3, 200)

# 능형 경로: 계수가 부드럽게 0 으로 다가가되 닿지는 않는다.
ridge_coefs = []
for a in alphas_path:
    model = Ridge(alpha=a).fit(X_scaled, y)
    ridge_coefs.append(model.coef_.copy())
ridge_coefs = np.array(ridge_coefs)

# 라쏘 경로: 계수가 하나씩 0 에 닿아 그대로 머문다. 두 그림의 이 차이가
# 곧 변수 선택을 하느냐 못 하느냐의 차이다.
lasso_coefs = []
alphas_lasso = np.logspace(-4, 1, 200)
for a in alphas_lasso:
    model = Lasso(alpha=a, max_iter=10000).fit(X_scaled, y)
    lasso_coefs.append(model.coef_.copy())
lasso_coefs = np.array(lasso_coefs)

print(f"능형 경로에서 계수가 정확히 0 인 횟수 = {(ridge_coefs == 0).sum()}")
print(f"능형 경로의 |계수| 최솟값 = {np.abs(ridge_coefs).min():.3e}  (0 이 아니다)")

n_nz = (np.abs(lasso_coefs) > 1e-8).sum(axis=1)
print(f"라쏘 경로의 0 아닌 계수 개수: {n_nz[0]} (alpha=1e-4) -> {n_nz[-1]} (alpha=10)")

lam_max = np.abs(X_scaled.T @ (y - y.mean())).max() / len(y)
print(f"이론 lambda_max = max_j |x_j'y|/n = {lam_max:.4f}")
print(f"격자에서 모든 계수가 0 이 된 가장 작은 alpha = "
      f"{alphas_lasso[np.where(n_nz == 0)[0][0]]:.4f}")

# LARS 는 꺾임점을 격자 없이 정확히 찾아 준다.
from sklearn.linear_model import lars_path
knots, _, _ = lars_path(X_scaled, y - y.mean(), method='lasso')
print(f"LARS 가 찾은 꺾임점 {len(knots)}개, 가장 큰 것 = {knots[0]:.4f}")

능형 경로에서 계수가 정확히 0 인 횟수 = 0
능형 경로의 |계수| 최솟값 = 9.924e-05  (0 이 아니다)
라쏘 경로의 0 아닌 계수 개수: 20 (alpha=1e-4) -> 0 (alpha=10)
이론 lambda_max = max_j |x_j'y|/n = 1.9458
격자에서 모든 계수가 0 이 된 가장 작은 alpha = 1.9792
LARS 가 찾은 꺾임점 23개, 가장 큰 것 = 1.9458

def plot_shrinkage_operators(lam=1.0):
    """정규직교 설계에서 세 축소 연산자가 최소제곱 추정값을 어떻게 바꾸는지 그린다.

    설계행렬이 정규직교이면 세 방법의 해가 최소제곱 추정값의 간단한 함수로
    나온다. 능형은 일정 비율로 줄이고(직선), 라쏘는 일정 크기만큼 깎아
    작은 값은 0 으로 보내며(꺾인 직선), 경성 문턱은 문턱 아래를 통째로
    0 으로 만든다(계단).
    """
    z = np.linspace(-4, 4, 500)
    ridge = z / (1 + lam)
    lasso = np.sign(z) * np.maximum(np.abs(z) - lam, 0)
    hard = z * (np.abs(z) > lam)

    fig, ax = plt.subplots(figsize=(7, 5))
    ax.plot(z, z, 'k--', alpha=0.3, label="OLS (no shrinkage)")
    ax.plot(z, ridge, linewidth=2, label=f"Ridge")
    ax.plot(z, lasso, linewidth=2, label=f"Lasso")
    ax.plot(z, hard, linewidth=2, label=f"Hard threshold")
    ax.set_xlabel("OLS estimate")
    ax.set_ylabel("Regularized estimate")
    ax.set_title("Shrinkage Operators (Orthonormal Design)")
    ax.legend()
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.show()

plot_shrinkage_operators()

# (1) 에서 손으로 계산한 세 값을 그대로 확인한다.
def shrink(z, lam=1.0):
    return (z / (1 + lam),
            np.sign(z) * np.maximum(np.abs(z) - lam, 0),
            z * (np.abs(z) > lam))

for z0 in (2.5, 0.6, -3.0):
    r, l, h = shrink(np.array(z0))
    print(f"OLS {z0:+.1f}  ->  능형 {r:+.3f}   라쏘 {l:+.3f}   경성 {h:+.3f}")

OLS +2.5  ->  능형 +1.250   라쏘 +1.500   경성 +2.500
OLS +0.6  ->  능형 +0.300   라쏘 +0.000   경성 +0.000
OLS -3.0  ->  능형 -1.500   라쏘 -2.000   경성 -3.000

import numpy as np
from sklearn.linear_model import RidgeCV, LassoCV, ElasticNetCV
from sklearn.preprocessing import StandardScaler

np.random.seed(42)
X, y, beta_true = generate_data(n=200, p=20, s=5, rho=0.9)
X_s = StandardScaler().fit_transform(X)

ridge = RidgeCV(alphas=np.logspace(-4, 2, 100), cv=5).fit(X_s, y)
lasso = LassoCV(n_alphas=100, cv=5, max_iter=10000).fit(X_s, y)
enet = ElasticNetCV(
    l1_ratio=[0.1, 0.5, 0.7, 0.9], n_alphas=100, cv=5
).fit(X_s, y)

for name, m in [("Ridge", ridge), ("Lasso", lasso), ("Elastic Net", enet)]:
    nz = np.sum(np.abs(m.coef_) > 1e-6)
    print(f"{name}: {nz} nonzero coefficients")

Ridge: 20 nonzero coefficients
Lasso: 9 nonzero coefficients
Elastic Net: 10 nonzero coefficients

import numpy as np
from sklearn.linear_model import Ridge, Lasso
from sklearn.preprocessing import StandardScaler

np.random.seed(0)
alphas_test = np.logspace(-3, 2, 30)
n_sim = 200
_, _, beta_true_bv = generate_data(n=2, p=20, s=5)

ridge_mse = {a: [] for a in alphas_test}
lasso_mse = {a: [] for a in alphas_test}

for _ in range(n_sim):
    X_sim, y_sim, _ = generate_data(n=100, p=20, s=5)
    X_sim = StandardScaler().fit_transform(X_sim)
    for a in alphas_test:
        r = Ridge(alpha=a).fit(X_sim, y_sim)
        l = Lasso(alpha=a, max_iter=5000).fit(X_sim, y_sim)
        ridge_mse[a].append(np.sum((r.coef_ - beta_true_bv)**2))
        lasso_mse[a].append(np.sum((l.coef_ - beta_true_bv)**2))

ridge_avg = {a: np.mean(v) for a, v in ridge_mse.items()}
lasso_avg = {a: np.mean(v) for a, v in lasso_mse.items()}

best_ridge = min(ridge_avg, key=ridge_avg.get)
best_lasso = min(lasso_avg, key=lasso_avg.get)
print(f"Ridge best lambda: {best_ridge:.4f}, MSE: {ridge_avg[best_ridge]:.4f}")
print(f"Lasso best lambda: {best_lasso:.4f}, MSE: {lasso_avg[best_lasso]:.4f}")

Ridge best lambda: 0.8532, MSE: 1.1883
Lasso best lambda: 0.0161, MSE: 0.8138