Chapter 59: Threshold Optimization

import numpy as np
from sklearn.linear_model import LogisticRegression

from kaggle_companion.activities._common import clean

SEED = 59
PREVALENCE = 0.05            # share of positive rows when the model and the threshold are fitted
DRAWS = 300                  # independent development sets per size; results are means over draws
AUDIT_ROWS = 50000           # fresh rows used only for scoring
SHIFTS = [0.01, 0.05, 0.20]  # audit prevalences used to test the frozen threshold


def generate(n, prevalence, rng):
    """Five noisy features; each shifts a little for positive rows."""
    y = (rng.random(n) < prevalence).astype(int)
    X = rng.normal(size=(n, 5)) + y[:, None] * np.array([1.0, 0.8, 0.6, 0.4, 0.2])
    return X, y


class Scored:
    """Rows sorted by score, so F1 at any threshold costs one search instead of one pass."""
    def __init__(self, p, y):
        order = np.argsort(p)
        self.p, self.y = p[order], y[order]
        self.positives_from = np.r_[np.cumsum(self.y[::-1])[::-1], 0]   # positives among rows at or above each position
        self.total = int(y.sum())

    def f1(self, threshold):
        """F1 when every row with score >= threshold is called positive."""
        i = np.searchsorted(self.p, threshold, side="left")
        called, hits = len(self.p) - i, self.positives_from[i]
        return 2 * hits / (called + self.total) if called + self.total else 0.0

    def best_threshold(self):
        """The threshold that maximizes F1 on these rows: a midpoint between two neighbouring scores."""
        k = np.arange(len(self.p), 0, -1)                       # rows called positive when cutting at each position
        f1 = 2 * self.positives_from[:-1] / (k + self.total)
        i = int(np.argmax(f1))
        return (self.p[i - 1] + self.p[i]) / 2 if i > 0 else self.p[0] - 1e-9, float(f1[i])


def averaged_fold_threshold(p, y, rng, folds=5):
    """A variant of the chapter's fold aggregation: tune on each fold's training part (four fifths of the rows), then average."""
    order = rng.permutation(len(p))
    found = []
    for k in range(folds):
        keep = np.setdiff1d(order, order[k::folds])
        if y[keep].sum() > 0:
            found.append(Scored(p[keep], y[keep]).best_threshold()[0])
    return float(np.mean(found))


def run(development_rows):
    rng = np.random.default_rng(SEED)
    X, y = generate(5000, PREVALENCE, rng)
    model = LogisticRegression(max_iter=500).fit(X, y)             # fitted once; never sees development or audit rows
    audits = {}
    for prevalence in SHIFTS:
        Xa, ya = generate(AUDIT_ROWS, prevalence, rng)
        audits[prevalence] = Scored(model.predict_proba(Xa)[:, 1], ya)
    audit = audits[PREVALENCE]

    tuned_dev, tuned_audit, fold_audit, thresholds = [], [], [], []
    shifted = {s: [] for s in SHIFTS}
    for _ in range(DRAWS):
        Xd, yd = generate(development_rows, PREVALENCE, rng)
        if yd.sum() == 0:                                           # a draw with no positive row cannot be tuned: give it one
            yd[rng.integers(development_rows)] = 1
        dev = Scored(model.predict_proba(Xd)[:, 1], yd)
        threshold, dev_f1 = dev.best_threshold()                    # fitted on the development rows' own labels
        tuned_dev.append(dev_f1)
        tuned_audit.append(audit.f1(threshold))                     # frozen, then scored on fresh labels
        thresholds.append(threshold)
        fold_audit.append(audit.f1(averaged_fold_threshold(dev.p, dev.y, rng)))
        for s in SHIFTS:
            shifted[s].append(audits[s].f1(threshold))
    return clean({
        "development_rows": development_rows,
        "development_f1": float(np.mean(tuned_dev)), "development_f1_sd": float(np.std(tuned_dev)),
        "audit_f1": float(np.mean(tuned_audit)), "audit_f1_sd": float(np.std(tuned_audit)),
        "fold_average_audit_f1": float(np.mean(fold_audit)),
        "default_audit_f1": audit.f1(0.5), "best_possible_audit_f1": audit.best_threshold()[1],
        "best_possible_threshold": audit.best_threshold()[0],
        "threshold_mean": float(np.mean(thresholds)), "threshold_sd": float(np.std(thresholds)),
        "tuned_beats_default": float(np.mean(np.array(tuned_audit) > audit.f1(0.5))),
        "shift": {str(s): {"frozen": float(np.mean(shifted[s])), "default": audits[s].f1(0.5),
                           "best_possible": audits[s].best_threshold()[1]} for s in SHIFTS},
        "draws": DRAWS, "audit_rows": AUDIT_ROWS,
    })

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

VALUES = [100, 300, 1000, 5000]
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.ch59 import draw, NCOLS, HEIGHT

DEFAULT = 100
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()