정칙화 로지스틱 회귀 실습¶

from sklearn.linear_model import LogisticRegression
import numpy as np

# 변수 20개 중 참으로 쓰이는 것은 앞의 다섯뿐이다. 세 벌점이 이 다섯을
# 어떻게 다루는지 견준다.
np.random.seed(42)
n, p = 200, 20
X = np.random.randn(n, p)
true_beta = np.zeros(p)
true_beta[:5] = [1.5, -1.0, 0.8, -0.5, 0.3]
logit = X @ true_beta
prob = 1 / (1 + np.exp(-logit))
y = np.random.binomial(1, prob)

# sklearn 에서 C 는 벌점의 역수다. C 가 작을수록 벌점이 세다 —
# 다른 책의 lambda 와 방향이 반대이니 헷갈리기 쉽다.
ridge_model = LogisticRegression(penalty='l2', C=1.0, solver='lbfgs',
                                  max_iter=1000)
ridge_model.fit(X, y)
print("Ridge coefficients:", np.round(ridge_model.coef_[0], 3))

# 정상조건 X^T r = beta/C 를 직접 확인한다.
r = y - ridge_model.predict_proba(X)[:, 1]
print(f"\n잔차의 합 = {r.sum():.3e}   (절편은 벌점을 받지 않는다)")
print(f"X^T r  첫 5개 = {np.round(X.T @ r, 4)[:5]}")
print(f"beta   첫 5개 = {np.round(ridge_model.coef_[0], 4)[:5]}")
print(f"max |X^T r - beta| = {np.abs(X.T @ r - ridge_model.coef_[0]).max():.3e}")

# 수렴 기준을 조이면 등식이 더 정확해진다 — 어긋남은 수학이 아니라 풀이기다.
tight = LogisticRegression(penalty='l2', C=1.0, solver='lbfgs',
                           max_iter=20000, tol=1e-8).fit(X, y)
r2 = y - tight.predict_proba(X)[:, 1]
print(f"tol=1e-8 로 다시: max |X^T r - beta| = "
      f"{np.abs(X.T @ r2 - tight.coef_[0]).max():.3e}")

# 벌점을 완전히 끈 최대가능도와 견준다.
mle = LogisticRegression(penalty=None, solver='lbfgs', max_iter=5000).fit(X, y)
print(f"\nMLE  첫 5개 = {np.round(mle.coef_[0], 3)[:5]}")
print(f"능형 노름 {np.linalg.norm(ridge_model.coef_[0]):.4f}, "
      f"MLE 노름 {np.linalg.norm(mle.coef_[0]):.4f}, "
      f"비 {np.linalg.norm(ridge_model.coef_[0]) / np.linalg.norm(mle.coef_[0]):.4f}")

Ridge coefficients: [ 1.346 -1.162  0.955 -0.593 -0.031 -0.187 -0.009 -0.243  0.359  0.041
 -0.176  0.021  0.069  0.145  0.398  0.04   0.047  0.095 -0.046 -0.166]

잔차의 합 = -3.266e-03   (절편은 벌점을 받지 않는다)
X^T r  첫 5개 = [ 1.3531 -1.1535  0.9554 -0.5991 -0.0282]
beta   첫 5개 = [ 1.3458 -1.1624  0.9554 -0.5931 -0.031 ]
max |X^T r - beta| = 8.861e-03
tol=1e-8 로 다시: max |X^T r - beta| = 6.020e-06

MLE  첫 5개 = [ 1.519 -1.318  1.071 -0.673 -0.054]
능형 노름 2.2160, MLE 노름 2.5020, 비 0.8857

# L1 벌점은 lbfgs 로 풀 수 없다. 0 에서 미분이 되지 않기 때문이며,
# 그래서 saga 같은 다른 풀이기를 써야 한다.
lasso_model = LogisticRegression(penalty='l1', C=1.0, solver='saga',
                                  max_iter=5000)
lasso_model.fit(X, y)
print("Lasso coefficients:", np.round(lasso_model.coef_[0], 3))
print(f"Non-zero coefficients: {np.sum(lasso_model.coef_[0] != 0)} / {p}")

r1 = y - lasso_model.predict_proba(X)[:, 1]
s = X.T @ r1
nz = lasso_model.coef_[0] != 0
print(f"\n0 이 된 계수의 번호: {np.where(~nz)[0]}")
print(f"그 셋의 |s_j| = {np.round(np.abs(s[~nz]), 4)}   (모두 1 이하여야 한다)")
print(f"0 이 아닌 {nz.sum()}개의 |s_j| 범위 = "
      f"[{np.abs(s[nz]).min():.4f}, {np.abs(s[nz]).max():.4f}]")
print(f"그 중 1 에서 가장 많이 벗어난 값 = "
      f"{np.abs(np.abs(s[nz]) - 1.0).max():.4f}")
print(f"부호가 계수와 일치하는가: "
      f"{bool(np.all(np.sign(s[nz]) == np.sign(lasso_model.coef_[0][nz])))}")

Lasso coefficients: [ 1.34  -1.146  0.94  -0.579  0.    -0.158  0.    -0.211  0.344  0.028
 -0.134  0.     0.023  0.092  0.37   0.006  0.007  0.055 -0.003 -0.127]
Non-zero coefficients: 17 / 20

0 이 된 계수의 번호: [ 4  6 11]
그 셋의 |s_j| = [0.8042 0.2508 0.6809]   (모두 1 이하여야 한다)
0 이 아닌 17개의 |s_j| 범위 = [0.9985, 1.0023]
그 중 1 에서 가장 많이 벗어난 값 = 0.0023
부호가 계수와 일치하는가: True

# 엘라스틱넷은 l1_ratio 로 두 벌점의 배합비를 정한다. 0.5 면 절반씩이다.
enet_model = LogisticRegression(penalty='elasticnet', C=1.0,
                                 solver='saga', l1_ratio=0.5,
                                 max_iter=5000)
enet_model.fit(X, y)
print("Elastic Net coefficients:", np.round(enet_model.coef_[0], 3))

be = enet_model.coef_[0]
re = y - enet_model.predict_proba(X)[:, 1]
se = X.T @ re
nze = be != 0
pred = 0.5 * be[nze] + 0.5 * np.sign(be[nze])
print(f"\n0 이 된 계수의 번호: {np.where(~nze)[0]}, "
      f"|s_j| = {np.round(np.abs(se[~nze]), 4)}  (문턱 0.5 이하)")
print(f"0 이 아닌 {nze.sum()}개:  "
      f"max |s_j - (0.5 b + 0.5 sign b)| = {np.abs(se[nze] - pred).max():.4f}")

bl, br = lasso_model.coef_[0], ridge_model.coef_[0]
print(f"\n0 이 아닌 개수   능형 {int((br != 0).sum())}  "
      f"엘넷 {int(nze.sum())}  라쏘 {int((bl != 0).sum())}")
print(f"L1 노름         능형 {np.abs(br).sum():.4f}  "
      f"엘넷 {np.abs(be).sum():.4f}  라쏘 {np.abs(bl).sum():.4f}")
print(f"|라쏘| <= |엘넷| 이 20개 모두 성립: "
      f"{bool(np.all(np.abs(bl) <= np.abs(be) + 1e-9))}")
print(f"|엘넷| <= |능형| 이 20개 모두 성립: "
      f"{bool(np.all(np.abs(be) <= np.abs(br) + 1e-9))}")

Elastic Net coefficients: [ 1.342 -1.153  0.947 -0.585 -0.012 -0.172  0.    -0.229  0.351  0.035
 -0.154  0.006  0.046  0.118  0.385  0.022  0.027  0.074 -0.026 -0.146]

0 이 된 계수의 번호: [6], |s_j| = [0.0184]  (문턱 0.5 이하)
0 이 아닌 19개:  max |s_j - (0.5 b + 0.5 sign b)| = 0.0014

0 이 아닌 개수   능형 20  엘넷 19  라쏘 17
L1 노름         능형 6.1297  엘넷 5.8315  라쏘 5.5625
|라쏘| <= |엘넷| 이 20개 모두 성립: True
|엘넷| <= |능형| 이 20개 모두 성립: True

import matplotlib.pyplot as plt

# C 를 키우며 계수 경로를 그린다. 참으로 쓰인 다섯 변수는 굵게, 나머지는
# 흐리게 그려 어느 쪽이 먼저 살아나는지 보이게 한다.
C_values = np.logspace(-3, 3, 50)
coefs = []

for C in C_values:
    model = LogisticRegression(penalty='l2', C=C, solver='lbfgs',
                                max_iter=2000)
    model.fit(X, y)
    coefs.append(model.coef_[0])

coefs = np.array(coefs)

plt.figure(figsize=(10, 5))
for j in range(p):
    plt.plot(np.log10(C_values), coefs[:, j],
             linewidth=2 if j < 5 else 0.8,
             alpha=1.0 if j < 5 else 0.3)
plt.xlabel('log10(C)')
plt.ylabel('Coefficient value')
plt.title('Ridge Logistic Regression: Coefficient Paths')
plt.axhline(0, color='black', linestyle='--', linewidth=0.5)
plt.tight_layout()
plt.show()

# 그림의 양 끝에서 읽어야 할 수.
print(f"왼쪽 끝 log10(C) = -3:  max |beta| = {np.abs(coefs[0]).max():.5f}")
print(f"오른쪽 끝 log10(C) = 3: 첫 5개 = {np.round(coefs[-1][:5], 3)}")
print(f"                        MLE 첫 5개 = {np.round(mle.coef_[0][:5], 3)}")
print(f"잡음 15개의 끝값 범위 = [{coefs[-1][5:].min():.3f}, "
      f"{coefs[-1][5:].max():.3f}]")
print(f"그 가운데 절댓값이 가장 큰 것 = "
      f"{np.abs(coefs[-1][5:]).max():.3f} (번호 "
      f"{5 + int(np.argmax(np.abs(coefs[-1][5:])))})")
print(f"참 신호 5번(beta5=0.3)의 끝값 = {coefs[-1][4]:.3f}")

왼쪽 끝 log10(C) = -3:  max |beta| = 0.03406
오른쪽 끝 log10(C) = 3: 첫 5개 = [ 1.519 -1.318  1.07  -0.673 -0.054]
                        MLE 첫 5개 = [ 1.519 -1.318  1.071 -0.673 -0.054]
잡음 15개의 끝값 범위 = [-0.271, 0.436]
그 가운데 절댓값이 가장 큰 것 = 0.436 (번호 14)
참 신호 5번(beta5=0.3)의 끝값 = -0.054

from sklearn.linear_model import LogisticRegressionCV

# C 는 교차검증으로 고른다. Cs=20 은 격자 점의 개수이며, sklearn 이
# 알아서 로그 눈금으로 펼친다.
model_cv = LogisticRegressionCV(
    Cs=20, penalty='l2', cv=5, scoring='accuracy',
    solver='lbfgs', max_iter=2000
)
model_cv.fit(X, y)
print(f"Best C: {model_cv.C_[0]:.4f}")
print(f"Best CV accuracy: {model_cv.scores_[1].mean(axis=0).max():.4f}")

grid = np.logspace(-4, 4, 20)
print(f"\nCs_ 가 logspace(-4, 4, 20) 인가: "
      f"{bool(np.allclose(model_cv.Cs_, grid))}")
k = int(np.argmin(np.abs(grid - model_cv.C_[0])))
print(f"고른 점은 k = {k},  10^(-4 + 8*{k}/19) = {10 ** (-4 + 8 * k / 19):.6f}")
print(f"이웃 격자점 C_9 = {grid[9]:.4f}, C_11 = {grid[11]:.4f}  "
      f"(한 칸 {10 ** (8 / 19):.3f}배)")

mean_scores = model_cv.scores_[1].mean(axis=0)
print(f"\n격자별 CV 정확도 = {np.round(mean_scores, 3)}")
print(f"최댓값 {mean_scores.max():.4f} 를 주는 격자점 = "
      f"{np.where(mean_scores == mean_scores.max())[0]}")
print(f"맞힌 건수 = {mean_scores.max():.4f} * {len(y)} = "
      f"{mean_scores.max() * len(y):.0f}")

refit = LogisticRegression(penalty='l2', C=model_cv.C_[0],
                           solver='lbfgs', max_iter=2000).fit(X, y)
print(f"\n훈련 정확도   {refit.score(X, y):.4f}")
print(f"교차검증 정확도 {mean_scores.max():.4f}")
print(f"다수범주 기준선 {max(y.mean(), 1 - y.mean()):.4f}")

Best C: 1.6238
Best CV accuracy: 0.7450

Cs_ 가 logspace(-4, 4, 20) 인가: True
고른 점은 k = 10,  10^(-4 + 8*10/19) = 1.623777
이웃 격자점 C_9 = 0.6158, C_11 = 4.2813  (한 칸 2.637배)

격자별 CV 정확도 = [0.535 0.555 0.645 0.71  0.72  0.72  0.72  0.73  0.74  0.74  0.745 0.745
 0.74  0.735 0.735 0.735 0.735 0.735 0.735 0.735]
최댓값 0.7450 를 주는 격자점 = [10 11]
맞힌 건수 = 0.7450 * 200 = 149

훈련 정확도   0.8000
교차검증 정확도 0.7450
다수범주 기준선 0.5050

import numpy as np
from sklearn.linear_model import LogisticRegression

np.random.seed(42)
n, p = 200, 50
X = np.random.randn(n, p)
true_beta = np.zeros(p)
true_beta[:5] = [2.0, -1.5, 1.0, -0.8, 0.5]
logit = X @ true_beta
prob = 1 / (1 + np.exp(-logit))
y = np.random.binomial(1, prob)

for C in [0.01, 0.1, 1.0, 10.0]:
    model = LogisticRegression(penalty='l1', C=C, solver='saga',
                                max_iter=5000)
    model.fit(X, y)
    nnz = np.sum(model.coef_[0] != 0)
    print(f"C = {C:5.2f}: {nnz} non-zero coefficients out of {p}")

C =  0.01: 0 non-zero coefficients out of 50
C =  0.10: 4 non-zero coefficients out of 50
C =  1.00: 36 non-zero coefficients out of 50
C = 10.00: 48 non-zero coefficients out of 50

from sklearn.linear_model import LogisticRegressionCV

model_cv = LogisticRegressionCV(
    Cs=20, penalty='l1', cv=5, scoring='accuracy',
    solver='saga', max_iter=5000
)
model_cv.fit(X, y)
print(f"Best C: {model_cv.C_[0]:.4f}")

best_idx = np.argmax(model_cv.scores_[1].mean(axis=0))
best_acc = model_cv.scores_[1].mean(axis=0)[best_idx]
print(f"Best CV accuracy: {best_acc:.4f}")
print(f"Non-zero coefficients: "
      f"{np.sum(model_cv.coef_[0] != 0)} / {p}")

Best C: 0.0886
Best CV accuracy: 0.8000
Non-zero coefficients: 4 / 50