불균형 자료 다루기¶

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(0)
n = 4000

# 신용점수·부채비율·소득이 연체 위험을 결정한다
score = rng.normal(0, 1, n)
dti = rng.normal(0, 1, n)
income = rng.normal(0, 1, n)
X = np.column_stack([score, dti, income])

# 절편 -2.0이 연체율을 약 19%로 맞춘다
logit = -2.0 - 1.1 * score + 0.9 * dti - 0.5 * income
p_default = 1 / (1 + np.exp(-logit))
y = (rng.random(n) < p_default).astype(int)   # 1 = 연체, 0 = 상환

X_tmp, X_test, y_tmp, y_test = train_test_split(
    X, y, test_size=0.25, random_state=0, stratify=y)
X_train, X_val, y_train, y_val = train_test_split(
    X_tmp, y_tmp, test_size=0.25, random_state=0, stratify=y_tmp)

print(f"n = {n}, 연체율 = {y.mean():.3f}")
print(f"train {len(y_train)}, val {len(y_val)}, test {len(y_test)}")

from scipy import integrate, stats

sd = np.sqrt(1.1**2 + 0.9**2 + 0.5**2)
g = lambda w: 1 / (1 + np.exp(-(-2 + sd * w))) * stats.norm.pdf(w)
theo = integrate.quad(g, -10, 10)[0]
print(f"\nVar(Z) = 2.27,  sd = {sd:.6f}")
print(f"이론 연체율 = {theo:.6f}   sigma(-2) = "
      f"{1 / (1 + np.exp(2)):.6f}")
print(f"관측 연체율 = {y.mean():.6f}  "
      f"(표준오차 {np.sqrt(theo * (1 - theo) / n):.5f})")
print(f"\n조각별 연체율  train {y_train.mean():.6f}  "
      f"val {y_val.mean():.6f}  test {y_test.mean():.6f}")
print(f"조각별 연체 건수  train {y_train.sum()}  val {y_val.sum()}  "
      f"test {y_test.sum()}   합 {y.sum()}")

n = 4000, 연체율 = 0.191
train 2250, val 750, test 1000

Var(Z) = 2.27,  sd = 1.506652
이론 연체율 = 0.190528   sigma(-2) = 0.119203
관측 연체율 = 0.190750  (표준오차 0.00621)

조각별 연체율  train 0.190667  val 0.190667  test 0.191000
조각별 연체 건수  train 429  val 143  test 191   합 763

from sklearn.linear_model import LogisticRegression

# 방법 1: 'balanced' 는 각 범주의 빈도에 반비례해 가중값을 자동으로 정한다.
model = LogisticRegression(class_weight='balanced')

# 방법 2: 직접 정한다. 5.3 은 위 자동값과 비슷하게 맞춘 것이다.
# 가중을 주면 모형이 소수 범주를 더 자주 예측하게 되는데, 이는 실은
# 문턱값을 옮기는 것과 비슷한 일을 하는 셈이다.
weights = [5.3 if yi == 1 else 1.0 for yi in y_train]
model.fit(X_train, y_train, sample_weight=weights)

print("가중 없음 예측 연체율:",
      LogisticRegression().fit(X_train, y_train).predict(X_val).mean().round(4))
print("가중 적용 예측 연체율:", model.predict(X_val).mean().round(4))

n1, n0 = int(y_train.sum()), int((1 - y_train).sum())
eff = 5.3 * n0 / n1
print(f"\n연체 {n1}건, 상환 {n0}건  ->  balanced 비 = {n0 / n1:.6f}")
print(f"실효 가중비 = 5.3 * {n0 / n1:.4f} = {eff:.4f},  "
      f"log = {np.log(eff):.4f}")

plain = LogisticRegression().fit(X_train, y_train)
print(f"\n가중 없음  절편 {plain.intercept_[0]:.6f}  "
      f"기울기 {np.round(plain.coef_[0], 4)}")
print(f"가중 적용  절편 {model.intercept_[0]:.6f}  "
      f"기울기 {np.round(model.coef_[0], 4)}")
print(f"절편 이동 = {model.intercept_[0] - plain.intercept_[0]:.4f}"
      f"   (예측 {np.log(eff):.4f})")

tau = 1 / (1 + eff)
p_plain = plain.predict_proba(X_val)[:, 1]
print(f"\n같은 문턱 tau = 1/(1+{eff:.4f}) = {tau:.6f}")
print(f"가중 없음 모형을 tau 로 자른 예측 연체율 = {(p_plain >= tau).mean():.4f}")
print(f"가중 모형을 0.5 로 자른 예측 연체율      = "
      f"{model.predict(X_val).mean():.4f}")
print(f"두 분류가 엇갈린 관측 = "
      f"{int((model.predict(X_val) != (p_plain >= tau)).sum())} / {len(y_val)}")

가중 없음 예측 연체율: 0.096
가중 적용 예측 연체율: 0.776

연체 429건, 상환 1821건  ->  balanced 비 = 4.244755
실효 가중비 = 5.3 * 4.2448 = 22.4972,  log = 3.1134

가중 없음  절편 -1.948054  기울기 [-0.9578  0.9223 -0.5285]
가중 적용  절편 1.139232  기울기 [-1.0132  0.9481 -0.4773]
절편 이동 = 3.0873   (예측 3.1134)

같은 문턱 tau = 1/(1+22.4972) = 0.042558
가중 없음 모형을 tau 로 자른 예측 연체율 = 0.7960
가중 모형을 0.5 로 자른 예측 연체율      = 0.7760
두 분류가 엇갈린 관측 = 19 / 750

from imblearn.over_sampling import SMOTE

# SMOTE 는 소수 범주의 관측값 사이를 이어 새 점을 만들어 채운다. 단순
# 복제와 달리 같은 점이 겹치지 않는다는 것이 장점이다. 다만 만들어 낸
# 점은 실제 관측이 아니므로, 반드시 훈련자료에만 적용해야 한다.
# 검증·시험자료에 쓰면 성능이 부풀려진다.
X_resampled, y_resampled = SMOTE(random_state=0).fit_resample(X_train, y_train)

model = LogisticRegression()
model.fit(X_resampled, y_resampled)

print(f"원자료:   n = {len(y_train)}, 연체율 = {y_train.mean():.3f}")
print(f"SMOTE 후: n = {len(y_resampled)}, 연체율 = {y_resampled.mean():.3f}")

print(f"\n소수 {n1} -> {int(y_resampled.sum())},  "
      f"합성점 {int(y_resampled.sum()) - n1}개  (= {n0} - {n1})")
print(f"다수 {n0} -> {int((1 - y_resampled).sum())} (그대로)")
print(f"n = {n0} + {n0} = {2 * n0} = {len(y_resampled)}")

원자료:   n = 2250, 연체율 = 0.191
SMOTE 후: n = 3642, 연체율 = 0.500

소수 429 -> 1821,  합성점 1392개  (= 1821 - 429)
다수 1821 -> 1821 (그대로)
n = 1821 + 1821 = 3642 = 3642

from imblearn.over_sampling import BorderlineSMOTE, ADASYN

# BorderlineSMOTE 는 두 범주의 경계 가까이에 있는 점만 골라 늘린다.
# 경계에서 먼 점은 어차피 쉽게 맞히므로 늘려도 보탬이 적다는 생각이다.
X_bl, y_bl = BorderlineSMOTE(random_state=0).fit_resample(X_train, y_train)

# ADASYN 은 분류하기 어려운 점일수록 더 많이 늘린다.
X_ad, y_ad = ADASYN(random_state=0).fit_resample(X_train, y_train)

print(f"BorderlineSMOTE: n = {len(y_bl)}, 연체율 = {y_bl.mean():.3f}")
print(f"ADASYN:          n = {len(y_ad)}, 연체율 = {y_ad.mean():.3f}")

target = n0 - n1
for name, yy in [("BorderlineSMOTE", y_bl), ("ADASYN", y_ad)]:
    made = int(yy.sum()) - n1
    print(f"{name:16s} 소수 {n1} -> {int(yy.sum())},  만든 점 {made} "
          f"(목표 {target}, 차이 {made - target}),  "
          f"연체율 {yy.mean():.6f}")

BorderlineSMOTE: n = 3642, 연체율 = 0.500
ADASYN:          n = 3635, 연체율 = 0.499
BorderlineSMOTE  소수 429 -> 1821,  만든 점 1392 (목표 1392, 차이 0),  연체율 0.500000
ADASYN           소수 429 -> 1814,  만든 점 1385 (목표 1392, 차이 -7),  연체율 0.499037

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (roc_auc_score, brier_score_loss, roc_curve,
                             confusion_matrix)

# 여기부터가 이 절의 결론에 해당하는 절차다. 재표집이나 가중 없이
# 원자료 그대로 적합한다. 불균형 자체는 확률 추정을 망가뜨리지 않는다.
model = LogisticRegression().fit(X_train, y_train)

# 2단계: 순위 매기는 능력(AUC)과 확률의 보정 상태(Brier)를 따로 본다.
# 불균형 자료에서 정확도는 뜻이 없다 — 전부 음성이라 해도 81%가 나온다.
p_val = model.predict_proba(X_val)[:, 1]
print("AUC  :", roc_auc_score(y_val, p_val))
print("Brier:", brier_score_loss(y_val, p_val))

# 3단계: 문턱값은 통계가 아니라 비용이 정한다. 거짓양성 50, 거짓음성
# 1000 이면 최적 문턱은 50/(50+1000) = 0.048 이다. 0.5 를 쓰는 관행에는
# 두 오류의 비용이 같다는 가정이 숨어 있다.
c_fp, c_fn = 50.0, 1000.0
threshold = c_fp / (c_fp + c_fn)

# 4단계: 시험자료는 마지막에 딱 한 번만 쓴다. 여기서 결과를 보고 다시
# 손대면 시험자료가 사실상 검증자료가 되어 버린다.
y_test_pred = (model.predict_proba(X_test)[:, 1] >= threshold).astype(int)

print(f"\ntau* = {c_fp}/({c_fp}+{c_fn}) = {threshold:.6f}")
p_test = model.predict_proba(X_test)[:, 1]
print(f"\n{'문턱':>10} {'TN':>5} {'FP':>5} {'FN':>5} {'TP':>5} "
      f"{'총비용':>9} {'정확도':>8} {'재현율':>8}")
for lab, t in [("전부 음성", 1.01), ("0.5", 0.5), ("tau*", threshold)]:
    pr = (p_test >= t).astype(int)
    cm_t = confusion_matrix(y_test, pr, labels=[0, 1])
    TN, FP, FN, TP = cm_t[0, 0], cm_t[0, 1], cm_t[1, 0], cm_t[1, 1]
    cost = c_fp * FP + c_fn * FN
    print(f"{lab:>10} {TN:5d} {FP:5d} {FN:5d} {TP:5d} {cost:9.0f} "
          f"{(TP + TN) / len(y_test):8.4f} {TP / (TP + FN):8.4f}")

best = min((c_fp * ((p_test >= g) & (y_test == 0)).sum()
            + c_fn * ((p_test < g) & (y_test == 1)).sum(), g)
           for g in np.unique(p_test))
print(f"\n검정자료에서 사후적으로 가장 싼 문턱 = {best[1]:.6f}, "
      f"그때 비용 {best[0]:.0f}")

AUC  : 0.8369258418681812
Brier: 0.111258445223284

tau* = 50.0/(50.0+1000.0) = 0.047619

        문턱    TN    FP    FN    TP       총비용      정확도      재현율
     전부 음성   809     0   191     0    191000   0.8090   0.0000
       0.5   788    21   137    54    138050   0.8420   0.2827
      tau*   204   605     7   184     37250   0.3880   0.9634

검정자료에서 사후적으로 가장 싼 문턱 = 0.046219, 그때 비용 35500

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, brier_score_loss

rng = np.random.default_rng(11)
n = 20000
x = rng.normal(size=n)
p = 1 / (1 + np.exp(-(-3.0 + 1.5 * x)))   # true beta0 = -3, beta1 = 1.5
y = (rng.random(n) < p).astype(int)

ntr = n // 2
Xtr, ytr = x[:ntr].reshape(-1, 1), y[:ntr]
Xte, yte = x[ntr:].reshape(-1, 1), y[ntr:]

def report(name, b0, b1):
    z = b0 + b1 * Xte.ravel()
    ph = 1 / (1 + np.exp(-z))
    print(f"{name:10s} b0={b0:7.4f} b1={b1:.4f} "
          f"AUC={roc_auc_score(yte, ph):.4f} "
          f"Brier={brier_score_loss(yte, ph):.5f} "
          f"mean_p={ph.mean():.4f}")

m = LogisticRegression(penalty=None).fit(Xtr, ytr)
report("plain", m.intercept_[0], m.coef_[0][0])

mb = LogisticRegression(penalty=None, class_weight='balanced').fit(Xtr, ytr)
report("balanced", mb.intercept_[0], mb.coef_[0][0])

i1 = np.where(ytr == 1)[0]
i0 = rng.choice(np.where(ytr == 0)[0], size=len(i1), replace=False)
ii = np.concatenate([i0, i1])
mu = LogisticRegression(penalty=None).fit(Xtr[ii], ytr[ii])
report("under", mu.intercept_[0], mu.coef_[0][0])

tau, pi = 0.5, ytr.mean()
b0_corr = mu.intercept_[0] - np.log(((1 - pi) / pi) * (tau / (1 - tau)))
report("corrected", b0_corr, mu.coef_[0][0])

plain      b0=-3.0714 b1=1.5111 AUC=0.8295 Brier=0.07247 mean_p=0.0910
balanced   b0=-0.8018 b1=1.5079 AUC=0.8295 Brier=0.16354 mean_p=0.3589
under      b0=-0.7882 b1=1.4686 AUC=0.8295 Brier=0.16282 mean_p=0.3596
corrected  b0=-3.0598 b1=1.4686 AUC=0.8295 Brier=0.07258 mean_p=0.0892