Chapter 39: Imbalanced Data

import numpy as np
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import log_loss, roc_auc_score

from kaggle_companion.activities._common import clean

SEED = 39
REPLICATES = 5          # independent worlds per keep rate; results are means over them
N_TRAIN, N_DEV, N_TEST = 40_000, 30_000, 60_000
EPS = 1e-6              # probabilities are kept inside [1e-6, 1 - 1e-6] so log loss stays finite


def generate(n, rng):
    """Eight features, a bent and interacting score, about 1.3% positives."""
    X = rng.normal(size=(n, 8))
    score = X[:, 0] + 0.8 * np.maximum(X[:, 1], 0) + 0.7 * X[:, 2] * X[:, 3] - 0.6 * np.abs(X[:, 4]) + 0.4 * X[:, 5]
    return X, (rng.random(n) < 1 / (1 + np.exp(-(1.2 * score - 5.6)))).astype(int)


def to_logit(p):
    p = np.clip(p, EPS, 1 - EPS)
    return np.log(p / (1 - p))


def correct(p, keep):
    """The chapter's correction: natural odds equal sampled odds times the negative keep rate r."""
    return p * keep / (p * keep + 1 - p)


def one_world(keep, seed):
    rng = np.random.default_rng(seed)
    X_train, y_train = generate(N_TRAIN, rng)
    X_dev, y_dev = generate(N_DEV, rng)                       # natural prevalence: used only to fit the calibrator
    X_test, y_test = generate(N_TEST, rng)                    # natural prevalence: used only for scoring
    sampled = (y_train == 1) | (rng.random(N_TRAIN) < keep)   # keep every positive, each negative with probability `keep`
    out = {"training_rows": int(sampled.sum()), "prevalence": float(y_test.mean())}
    makers = {"gbm": lambda: HistGradientBoostingClassifier(max_iter=60, max_leaf_nodes=15, random_state=0),
              "logistic": lambda: LogisticRegression(max_iter=500)}
    for name, make in makers.items():
        model = make().fit(X_train[sampled], y_train[sampled])
        p_dev = np.clip(model.predict_proba(X_dev)[:, 1], EPS, 1 - EPS)
        p_test = np.clip(model.predict_proba(X_test)[:, 1], EPS, 1 - EPS)
        calibrator = LogisticRegression(C=1e6).fit(to_logit(p_dev)[:, None], y_dev)      # Platt scaling on development rows
        outputs = {"raw": p_test, "corrected": np.clip(correct(p_test, keep), EPS, 1 - EPS),
                   "calibrated": calibrator.predict_proba(to_logit(p_test)[:, None])[:, 1]}
        out[name] = {k: [log_loss(y_test, v), roc_auc_score(y_test, v), float(v.mean())] for k, v in outputs.items()}
    return out


def run(keep_rate):
    worlds = [one_world(keep_rate, SEED * 100 + r) for r in range(REPLICATES)]
    # Reference: the same models trained on every row (keep rate 1), same worlds.
    base = worlds if keep_rate == 1 else [one_world(1, SEED * 100 + r) for r in range(REPLICATES)]

    def mean(group, model, arm, j):
        return float(np.mean([w[model][arm][j] for w in group]))
    result = {"keep_rate": keep_rate, "training_rows": float(np.mean([w["training_rows"] for w in worlds])),
              "prevalence": float(np.mean([w["prevalence"] for w in worlds])), "replicates": REPLICATES}
    for model in ("gbm", "logistic"):
        result[model] = {arm: {"log_loss": mean(worlds, model, arm, 0), "auc": mean(worlds, model, arm, 1),
                               "mean_prediction": mean(worlds, model, arm, 2)} for arm in ("raw", "corrected", "calibrated")}
        result[model]["no_downsampling_log_loss"] = mean(base, model, "raw", 0)
        result[model]["no_downsampling_auc"] = mean(base, model, "raw", 1)
    result["calibrated_beats_corrected_gbm"] = float(np.mean([w["gbm"]["calibrated"][0] < w["gbm"]["corrected"][0] for w in worlds]))
    return clean(result)

import pandas as pd
from kaggle_companion.activities.ch39 import explain

VALUES = [1, 0.02, 0.005, 0.002]
results = {value: run(value) for value in VALUES}
table = pd.DataFrame({str(v): explain(results[v], v)['metrics'] for v in VALUES})
print(table.to_string())

import matplotlib.pyplot as plt
from kaggle_companion.activities.ch39 import draw, NCOLS, HEIGHT

DEFAULT = 0.02
fig, axes = plt.subplots(1, NCOLS, figsize=(6.6 if NCOLS == 1 else 10.4, HEIGHT), layout='constrained')
draw(axes, results[DEFAULT], DEFAULT)
plt.show()