ROC 곡선과 평가지표 구현¶

import numpy as np

# 양성 120, 음성 180 인 자료를 만든다. 두 집단의 점수 분포가 겹치므로
# 어느 문턱값을 잡아도 오류가 생긴다 — 측도들이 갈리는 자리다.
np.random.seed(42)
n = 300
y_true = np.concatenate([np.ones(120), np.zeros(180)])
scores = np.concatenate([
    np.random.normal(0.65, 0.25, 120),
    np.random.normal(0.35, 0.25, 180)
])
scores = np.clip(scores, 0, 1)

from scipy import stats

q = stats.norm.cdf(0.6)
print(f"q = Phi(0.6) = {q:.6f}")
print(f"E[TP] = {120 * q:.2f}  E[FN] = {120 * (1 - q):.2f}  "
      f"E[FP] = {180 * (1 - q):.2f}  E[TN] = {180 * q:.2f}")
print(f"E[FP]/E[FN] = {180 * (1 - q) / (120 * (1 - q)):.4f}")
print(f"이론 AUC = Phi(0.30/(0.25*sqrt2)) "
      f"= Phi({0.3 / (0.25 * np.sqrt(2)):.6f}) = "
      f"{stats.norm.cdf(0.3 / (0.25 * np.sqrt(2))):.6f}")
print(f"\nclip 으로 0 에 몰린 수 = {(scores == 0).sum()}, "
      f"1 에 몰린 수 = {(scores == 1).sum()}")
print(f"서로 다른 점수의 개수 = {len(np.unique(scores))}"
      f"   (= 300 - 13 - 11 + 2)")

q = Phi(0.6) = 0.725747
E[TP] = 87.09  E[FN] = 32.91  E[FP] = 49.37  E[TN] = 130.63
E[FP]/E[FN] = 1.5000
이론 AUC = Phi(0.30/(0.25*sqrt2)) = Phi(0.848528) = 0.801928

clip 으로 0 에 몰린 수 = 13, 1 에 몰린 수 = 11
서로 다른 점수의 개수 = 278   (= 300 - 13 - 11 + 2)

def confusion_matrix(y_true, y_pred):
    """2x2 혼동행렬 [TN, FP; FN, TP] 을 구한다.

    행이 실제, 열이 예측이다. 아래의 모든 측도가 이 네 칸에서 나온다.
    """
    tp = np.sum((y_true == 1) & (y_pred == 1))
    tn = np.sum((y_true == 0) & (y_pred == 0))
    fp = np.sum((y_true == 0) & (y_pred == 1))
    fn = np.sum((y_true == 1) & (y_pred == 0))
    return np.array([[tn, fp], [fn, tp]])

y_pred = (scores >= 0.5).astype(int)
cm = confusion_matrix(y_true, y_pred)
print("[[TN FP]\n [FN TP]] =")
print(cm)

TN, FP, FN, TP = cm[0, 0], cm[0, 1], cm[1, 0], cm[1, 1]
sd_tp = np.sqrt(120 * q * (1 - q))
sd_fp = np.sqrt(180 * q * (1 - q))
print(f"\nTP {TP} 대 기댓값 {120 * q:.2f}  (SD {sd_tp:.3f}, "
      f"z = {(TP - 120 * q) / sd_tp:.3f})")
print(f"FP {FP} 대 기댓값 {180 * (1 - q):.2f}  (SD {sd_fp:.3f}, "
      f"z = {(FP - 180 * (1 - q)) / sd_fp:.3f})")
print(f"FP/FN = {FP}/{FN} = {FP / FN:.4f}   (기댓값의 비 1.5)")

[[TN FP]
 [FN TP]] =
[[129  51]
 [ 30  90]]

TP 90 대 기댓값 87.09  (SD 4.887, z = 0.596)
FP 51 대 기댓값 49.37  (SD 5.986, z = 0.273)
FP/FN = 51/30 = 1.7000   (기댓값의 비 1.5)

def precision_recall_f1(y_true, y_pred):
    """정밀도, 재현율, F1 을 구한다.

    정밀도는 "양성이라 한 것 중 맞은 비율", 재현율은 "실제 양성 중 잡아낸
    비율"이다. 둘은 서로 맞바꿈 관계라, 문턱을 낮추면 재현율이 오르고
    정밀도가 내린다. F1 은 그 둘의 조화평균이다.
    """
    cm = confusion_matrix(y_true, y_pred)
    tn, fp, fn, tp = cm[0, 0], cm[0, 1], cm[1, 0], cm[1, 1]
    precision = tp / (tp + fp) if (tp + fp) > 0 else 0.0
    recall = tp / (tp + fn) if (tp + fn) > 0 else 0.0
    f1 = (2 * precision * recall / (precision + recall)
          if (precision + recall) > 0 else 0.0)
    return precision, recall, f1

prec, rec, f1 = precision_recall_f1(y_true, y_pred)
accuracy = np.mean(y_true == y_pred)
print(f"정확도 {accuracy:.6f}  정밀도 {prec:.6f}  "
      f"재현율 {rec:.6f}  F1 {f1:.6f}")
print(f"분수로   219/300={219 / 300:.6f}  90/141={90 / 141:.6f}  "
      f"90/120={90 / 120:.6f}  180/261={180 / 261:.6f}")
print(f"FP {FP} > FN {FN} 이므로 정밀도 < 재현율")
print(f"정밀도의 이론값 = 120q/(120q+180(1-q)) = "
      f"{120 * q / (120 * q + 180 * (1 - q)):.6f}")

정확도 0.730000  정밀도 0.638298  재현율 0.750000  F1 0.689655
분수로   219/300=0.730000  90/141=0.638298  90/120=0.750000  180/261=0.689655
FP 51 > FN 30 이므로 정밀도 < 재현율
정밀도의 이론값 = 120q/(120q+180(1-q)) = 0.638229

def roc_curve(y_true, scores):
    """문턱값을 바꿔 가며 ROC 곡선을 그린다.

    문턱을 하나 정하는 대신 모든 문턱에서의 성적을 한 곡선에 담는다.
    그래서 ROC 는 문턱 선택과 무관하게 모형의 순위 매기는 능력만 잰다.
    """
    thresholds = np.sort(np.unique(scores))[::-1]
    fpr_list, tpr_list = [0.0], [0.0]
    n_pos = np.sum(y_true == 1)
    n_neg = np.sum(y_true == 0)

    for thresh in thresholds:
        y_pred = (scores >= thresh).astype(int)
        tp = np.sum((y_true == 1) & (y_pred == 1))
        fp = np.sum((y_true == 0) & (y_pred == 1))
        tpr_list.append(tp / n_pos if n_pos > 0 else 0)
        fpr_list.append(fp / n_neg if n_neg > 0 else 0)

    return np.array(fpr_list), np.array(tpr_list), thresholds

fpr, tpr, thresholds = roc_curve(y_true, scores)

print(f"len(fpr) = {len(fpr)}, len(tpr) = {len(tpr)}, "
      f"len(thresholds) = {len(thresholds)}")

optimal_idx = int(np.argmax(tpr - fpr))
print(f"\nargmax(tpr - fpr) = {optimal_idx}")
print(f"그 곡선 점:        TPR {tpr[optimal_idx]:.6f}  "
      f"FPR {fpr[optimal_idx]:.6f}  J {tpr[optimal_idx] - fpr[optimal_idx]:.6f}")
print(f"옳은 문턱  thresholds[{optimal_idx - 1}] = "
      f"{thresholds[optimal_idx - 1]:.7f}")
print(f"어긋난 문턱 thresholds[{optimal_idx}]  = "
      f"{thresholds[optimal_idx]:.7f}")

for name, t in [("옳은", thresholds[optimal_idx - 1]),
                ("어긋난", thresholds[optimal_idx])]:
    p = (scores >= t)
    tp_r = np.sum((y_true == 1) & p) / 120
    fp_r = np.sum((y_true == 0) & p) / 180
    print(f"{name} 문턱에서  TPR {tp_r:.6f}  FPR {fp_r:.6f}  "
          f"J {tp_r - fp_r:.6f}")

len(fpr) = 279, len(tpr) = 279, len(thresholds) = 278

argmax(tpr - fpr) = 133
그 곡선 점:        TPR 0.766667  FPR 0.283333  J 0.483333
옳은 문턱  thresholds[132] = 0.4995733
어긋난 문턱 thresholds[133]  = 0.4970793
옳은 문턱에서  TPR 0.766667  FPR 0.283333  J 0.483333
어긋난 문턱에서  TPR 0.766667  FPR 0.288889  J 0.477778

def auc_trapezoid(fpr, tpr):
    """사다리꼴 공식으로 곡선 아래 넓이를 구한다.

    AUC 는 "무작위로 고른 양성의 점수가 무작위로 고른 음성의 점수보다
    높을 확률"과 정확히 같다. 0.5 면 동전 던지기, 1 이면 완벽하다.
    """
    order = np.argsort(fpr)
    fpr_sorted = fpr[order]
    tpr_sorted = tpr[order]
    return np.trapz(tpr_sorted, fpr_sorted)

area = auc_trapezoid(fpr, tpr)
print(f"AUC = {area:.4f}")

# (1) 짝을 전부 세어 본다.
pos, neg = scores[y_true == 1], scores[y_true == 0]
d = pos[:, None] - neg[None, :]
wins, ties = int((d > 0).sum()), int((d == 0).sum())
print(f"\n동점 쌍 = 1*12 + 9*2 = "
      f"{((scores == 0) & (y_true == 1)).sum() * ((scores == 0) & (y_true == 0)).sum()}"
      f" + {((scores == 1) & (y_true == 1)).sum() * ((scores == 1) & (y_true == 0)).sum()}"
      f" = {ties}")
mw = (wins + 0.5 * ties) / (len(pos) * len(neg))
print(f"만-휘트니 = ({wins} + 0.5*{ties})/21600 = {mw:.10f}")

# (2) fpr 은 이미 오름차순이다. 정렬이 무엇을 하는지 본다.
print(f"\nfpr 이 이미 오름차순인가? {bool(np.all(np.diff(fpr) >= 0))}")
print(f"정렬 없이          : {np.trapz(tpr, fpr):.10f}")
for kind in ["quicksort", "stable", "heapsort"]:
    o = np.argsort(fpr, kind=kind)
    print(f"argsort({kind:9s}): {np.trapz(tpr[o], fpr[o]):.10f}")

AUC = 0.7935

동점 쌍 = 1*12 + 9*2 = 12 + 18 = 30
만-휘트니 = (17116 + 0.5*30)/21600 = 0.7931018519

fpr 이 이미 오름차순인가? True
정렬 없이          : 0.7931018519
argsort(quicksort): 0.7935416667
argsort(stable   ): 0.7931018519
argsort(heapsort ): 0.7929166667

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))

# 왼쪽: 혼동행렬
im = axes[0].imshow(cm, cmap='Blues', aspect='equal')
axes[0].set_xticks([0, 1]); axes[0].set_yticks([0, 1])
axes[0].set_xticklabels(['Neg', 'Pos'])
axes[0].set_yticklabels(['Neg', 'Pos'])
axes[0].set_xlabel('Predicted'); axes[0].set_ylabel('Actual')
axes[0].set_title('Confusion Matrix')

# 가운데: 두 집단의 점수 분포. 겹치는 만큼이 오류가 되고, 붉은 세로선을
# 좌우로 옮기는 것이 곧 문턱을 바꾸는 일이다.
axes[1].hist(scores[y_true == 0], bins=25, alpha=0.6,
             label='Negative', edgecolor='k')
axes[1].hist(scores[y_true == 1], bins=25, alpha=0.6,
             label='Positive', edgecolor='k')
axes[1].axvline(0.5, color='red', linestyle='--', label='Threshold')
axes[1].set_xlabel('Predicted score'); axes[1].set_ylabel('Frequency')
axes[1].set_title('Score Distributions'); axes[1].legend(fontsize=8)

# 오른쪽: ROC 곡선. 왼쪽 위 모서리에 가까울수록 좋다.
axes[2].plot(fpr, tpr, linewidth=2, label=f'AUC = {area:.3f}')
axes[2].plot([0, 1], [0, 1], 'k--', alpha=0.4, label='Random')
axes[2].set_xlabel('False Positive Rate')
axes[2].set_ylabel('True Positive Rate')
axes[2].set_title('ROC Curve'); axes[2].legend(fontsize=9)

plt.tight_layout()
plt.show()

# 그림에서 읽히지 않는 수를 따로 찍어 둔다.
print("왼쪽 칸의 네 값:", cm.ravel(), " (열지도에 숫자가 없다)")
print(f"가운데 칸 양 끝 막대:  점수 0 에 {(scores == 0).sum()}개"
      f"(음성 {((scores == 0) & (y_true == 0)).sum()}, "
      f"양성 {((scores == 0) & (y_true == 1)).sum()}),  "
      f"점수 1 에 {(scores == 1).sum()}개"
      f"(음성 {((scores == 1) & (y_true == 0)).sum()}, "
      f"양성 {((scores == 1) & (y_true == 1)).sum()})")
print(f"두 정규분포의 겹침 = 2*Phi(-0.6) = "
      f"{2 * stats.norm.cdf(-0.6):.4f}")
k = int(np.argmax(tpr >= 119 / 120 - 1e-12))
print(f"TPR 이 119/120 에 닿는 지점의 FPR = {fpr[k]:.4f}  "
      f"-> 거기서 오른쪽 끝까지 평평하다")
print(f"tau=0.5 인 곡선 위의 점: FPR {51 / 180:.4f}, TPR {90 / 120:.4f}"
      f"  (그림에는 표시가 없다)")

왼쪽 칸의 네 값: [129  51  30  90]  (열지도에 숫자가 없다)
가운데 칸 양 끝 막대:  점수 0 에 13개(음성 12, 양성 1),  점수 1 에 11개(음성 2, 양성 9)
두 정규분포의 겹침 = 2*Phi(-0.6) = 0.5485
TPR 이 119/120 에 닿는 지점의 FPR = 0.7722  -> 거기서 오른쪽 끝까지 평평하다
tau=0.5 인 곡선 위의 점: FPR 0.2833, TPR 0.7500  (그림에는 표시가 없다)

def pr_curve(y_true, scores):
    """문턱값을 바꿔 가며 정밀도-재현율 곡선을 구한다.

        양성이 드문 자료에서는 ROC 보다 이 곡선이 낫다. ROC 의 FPR 은 분모가
        음성 수라, 음성이 압도적으로 많으면 거의 움직이지 않기 때문이다.
        """
    thresholds = np.sort(np.unique(scores))[::-1]
    precisions, recalls = [1.0], [0.0]
    n_pos = np.sum(y_true == 1)

    for thresh in thresholds:
        y_pred = (scores >= thresh).astype(int)
        tp = np.sum((y_true == 1) & (y_pred == 1))
        fp = np.sum((y_true == 0) & (y_pred == 1))
        fn = np.sum((y_true == 1) & (y_pred == 0))
        prec = tp / (tp + fp) if (tp + fp) > 0 else 1.0
        rec = tp / (tp + fn) if (tp + fn) > 0 else 0.0
        precisions.append(prec)
        recalls.append(rec)

    return np.array(recalls), np.array(precisions), thresholds