import numpy as np, pandas as pd, warnings
warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier
from sklearn.metrics import (roc_auc_score, average_precision_score, brier_score_loss,
                             precision_recall_curve, roc_curve)
plt.rcParams.update({"figure.dpi":110,"font.size":11,"axes.spines.top":False,"axes.spines.right":False,
    "axes.grid":True,"grid.alpha":0.22,"axes.titleweight":"bold","axes.titlesize":12.5,"axes.titlelocation":"left"})
VI, DK, LT, MUT, GD, RD = "#6d28d9", "#4c1d95", "#a78bfa", "#94a3b8", "#047857", "#dc2626"

BASE = "https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
fn = "capstone-credit-card-fraud.xlsx"
def load(sheet):
    try:    return pd.read_excel("../../data/" + fn, sheet_name=sheet)
    except FileNotFoundError: return pd.read_excel(BASE + fn, sheet_name=sheet)
raw  = load("Transactions")
plan = load("AnalysisPlan")
print("rows in the warehouse export:", f"{len(raw):,}")
raw.head()

rows in the warehouse export: 120,400

for line in plan.analysis_plan.fillna(""):
    print(line)

AGREED WITH FRAUD OPERATIONS BEFORE ANY MODEL WAS FITTED.

THE DECISION. A flagged transaction is held and reviewed by an analyst. An unflagged one settles.
The model produces a queue, and the team can work about 400 reviews a day.

THE COST OF BEING WRONG, from operations, in dollars:
  a missed fraud        = the transaction amount, written off, plus 25 of handling
  a false positive      = 4 of analyst time, plus an estimated 12 of customer annoyance
  These are NOT symmetric and they are not even a fixed ratio: the cost of a miss depends
  on the amount. That rules out any single global threshold chosen from a curve.

M1 SPLIT: by time. The model will score future transactions, so it is validated on later ones.
   A random split would let it learn from transactions that had not happened yet.
M2 METRIC: PRECISION-RECALL AUC as the headline, never ROC-AUC. With 0.28 percent positives
   the false positive rate has a denominator of 120,000 and barely moves, so ROC flatters.
M3 IMBALANCE: three approaches compared on the SAME model and split, and judged on expected
   cost rather than on any ranking metric: (a) random oversampling of the minority, (b) class
   weights, (c) nothing at all, and tune the threshold.
M4 CALIBRATION: checked for each approach. Resampling and class weights both change the base
   rate the model believes in, so a probability from them is not a probability.
M5 THRESHOLD: from the cost matrix and the review capacity, not from Youden or an F1 peak.
M6 EXPLANATION: SHAP on the flagged transactions, used first as a LEAKAGE CHECK and only
   afterwards as an explanation for the analyst working the queue.
M7 EXCLUSIONS: duplicate rows, distance sentinels of -1, and zero-amount authorization holds.

WHAT WE ARE NOT DOING: reporting accuracy. Flagging nothing is 99.72 percent accurate.

d = raw.drop_duplicates(subset="transaction_id").copy()
print(f"raw rows                  {len(raw):,}")
print(f"after the duplicate export {len(d):,}")
d.loc[d.distance_from_home_km == -1, "distance_from_home_km"] = np.nan
d = d.dropna(subset=["distance_from_home_km"])
print(f"with geolocation          {len(d):,}")
d = d[d.amount > 0].copy()
print(f"real purchases            {len(d):,}   (zero-amount rows are authorization holds)")

print(f"\nfraud rate                {d.is_fraud.mean():.4%}   ({d.is_fraud.sum()} of {len(d):,})")
print(f"accuracy of flagging NOTHING: {1-d.is_fraud.mean():.4%}")
print("\nThat number is why accuracy is not mentioned again in this notebook.")

raw rows                  120,400
after the duplicate export 120,000
with geolocation          119,694
real purchases            119,556   (zero-amount rows are authorization holds)

fraud rate                0.2660%   (318 of 119,556)
accuracy of flagging NOTHING: 99.7340%

That number is why accuracy is not mentioned again in this notebook.

d["night"] = d.hour_of_day.between(1, 5).astype(int)
d["log_amt"] = np.log1p(d.amount)
d["dev_online"] = d.new_device * d.online
d["foreign_night"] = d.foreign_merchant * d.night

F  = ["log_amt","night","customer_tenure_days","transactions_last_30d","foreign_merchant",
      "online","new_device","distance_from_home_km","merchant_category_risk"]
FX = F + ["dev_online", "foreign_night"]

tr, te = d[d.transaction_day <= 90], d[d.transaction_day > 90]
ytr, yte = tr.is_fraud.values, te.is_fraud.values
mu, sd = tr[FX].mean(), tr[FX].std()
print(f"train  days 1-90    {len(tr):,} transactions, {ytr.sum()} frauds ({ytr.mean():.3%})")
print(f"test   days 91-120  {len(te):,} transactions, {yte.sum()} frauds ({yte.mean():.4%})")
print("\nA random split would let the model learn from transactions that had not happened yet.")

train  days 1-90    89,973 transactions, 223 frauds (0.248%)
test   days 91-120  29,583 transactions, 95 frauds (0.3211%)

A random split would let the model learn from transactions that had not happened yet.

smd = {c: (tr.loc[tr.is_fraud == 1, c].mean() - tr.loc[tr.is_fraud == 0, c].mean()) /
          np.sqrt((tr.loc[tr.is_fraud == 0, c].var() + tr.loc[tr.is_fraud == 1, c].var())/2)
       for c in F}
top = sorted(smd.items(), key=lambda kv: -abs(kv[1]))

print(f"training window, days 1 to 90: {len(tr):,} transactions, {int(ytr.sum())} frauds, "
      f"one in {len(tr)/max(ytr.sum(),1):,.0f}")
print(f"  amount    median {tr.amount.median():.2f}   "
      f"99th percentile {tr.amount.quantile(.99):.2f}   max {tr.amount.max():,.2f}   "
      f"skew {tr.amount.skew():+.1f}, which is why the model uses log_amt")
print("  strongest separators: " + ", ".join(f"{k} {v:+.2f}" for k, v in top[:3]))
print("  weakest: " + ", ".join(f"{k} {v:+.2f}" for k, v in top[-2:]))
print("  none of them separates the classes on its own, which is the ordinary case "
      "and the reason a model is worth building")

training window, days 1 to 90: 89,973 transactions, 223 frauds, one in 403
  amount    median 21.21   99th percentile 377.34   max 7,890.12   skew +18.0, which is why the model uses log_amt
  strongest separators: new_device +1.69, log_amt +0.86, distance_from_home_km +0.81
  weakest: customer_tenure_days -0.33, transactions_last_30d -0.08
  none of them separates the classes on its own, which is the ordinary case and the reason a model is worth building

fig, axes = plt.subplots(2, 2, figsize=(12.6, 8.2))

ax = axes[0, 0]
cnt = [int((ytr == 0).sum()), int((ytr == 1).sum())]
ax.bar([0, 1], cnt, color=[MUT, RD], width=0.5)
for x, c_ in zip([0, 1], cnt):
    ax.text(x, c_*1.25, f"{c_:,}", ha="center", fontsize=12, fontweight="bold", color=DK)
ax.set_yscale("log")
ax.set_xticks([0, 1]); ax.set_xticklabels(["genuine", "fraud"])
ax.set_ylim(1, cnt[0]*9)
ax.set_ylabel("transactions (log scale)")
ax.set_title(f"One class is {cnt[0]/max(cnt[1],1):,.0f} times the size of the other")
ax.grid(axis="x", alpha=0)

ax = axes[0, 1]
bins = np.linspace(0, 12, 45)
for lab, col, nm in [(0, MUT, "genuine"), (1, RD, "fraud")]:
    ax.hist(tr.loc[tr.is_fraud == lab, "log_amt"], bins=bins, density=True,
            color=col, alpha=0.60, label=nm)
ax.set_xlabel("log of the transaction amount"); ax.set_ylabel("share of that class")
ax.set_title("The classes overlap almost everywhere")
ax.legend(fontsize=9.5, loc="upper left")

ax = axes[1, 0]
vals = sorted(smd.items(), key=lambda kv: kv[1])
vv = [v for k, v in vals]
ax.barh(range(len(vv)), vv, color=[RD if v > 0 else GD for v in vv], height=0.66)
ax.set_yticks(range(len(vv)))
ax.set_yticklabels([k.replace("_", " ") for k, v in vals], fontsize=9.5)
for i, v_ in enumerate(vv):
    ax.text(v_ + (0.03 if v_ > 0 else -0.03), i, f"{v_:+.2f}", va="center",
            ha="left" if v_ > 0 else "right", fontsize=9, fontweight="bold", color=DK)
ax.axvline(0, color=DK, lw=1.6)
ax.set_xlim(min(vv)*2.20, max(vv)*1.30)
ax.set_xlabel("standardized difference, fraud minus genuine")
ax.set_title("What separates them, and none of it cleanly")
ax.grid(axis="y", alpha=0)

ax = axes[1, 1]
hrs = sorted(tr.hour_of_day.unique())
rate = [tr.loc[tr.hour_of_day == h, "is_fraud"].mean()*1000 for h in hrs]
vol  = [int((tr.hour_of_day == h).sum()) for h in hrs]
ax.bar(hrs, rate, color=[RD if 1 <= h <= 5 else VI for h in hrs], width=0.75)
ax.axhline(ytr.mean()*1000, color=DK, lw=1.8, ls="--")
ax.text(23.4, ytr.mean()*1000 + 0.4, f"overall {ytr.mean()*1000:.1f}", ha="right",
        fontsize=9.5, fontweight="bold", color=DK)
ax.set_xlabel("hour of day"); ax.set_ylabel("frauds per 1,000 transactions")
ax.set_ylim(0, max(rate)*1.20)
ax.set_title("Rare everywhere, and least rare overnight (red)")
ax.grid(axis="x", alpha=0)
plt.tight_layout(); plt.show()

def fit_lr(cols, train, y, cols_mu=None, cols_sd=None):
    m_, s_ = (train[cols].mean(), train[cols].std()) if cols_mu is None else (cols_mu, cols_sd)
    mod = LogisticRegression(max_iter=3000).fit((train[cols]-m_)/s_, y)
    return mod, m_, s_

for cols, tag in [(FX, "honest features"), (FX + ["chargeback_code_filed"], "+ chargeback_code_filed")]:
    mod, m_, s_ = fit_lr(cols, tr, ytr)
    p = mod.predict_proba((te[cols]-m_)/s_)[:, 1]
    print(f"  {tag:26s} ROC-AUC {roc_auc_score(yte,p):.4f}   PR-AUC {average_precision_score(yte,p):.4f}")
print(f"\n  a random ranking scores PR-AUC = the base rate = {yte.mean():.4f}")
print("\nA chargeback code is filed WEEKS LATER, once a human has worked the case.")
print("At the moment the model must decide whether to hold the transaction, it does not exist.")

  honest features            ROC-AUC 0.9604   PR-AUC 0.2984
  + chargeback_code_filed    ROC-AUC 0.9997   PR-AUC 0.9818

  a random ranking scores PR-AUC = the base rate = 0.0032

A chargeback code is filed WEEKS LATER, once a human has worked the case.
At the moment the model must decide whether to hold the transaction, it does not exist.

results = {}
mod, m_, s_ = fit_lr(F, tr, ytr)
results["logistic, additive"] = mod.predict_proba((te[F]-m_)/s_)[:, 1]
mod, m_, s_ = fit_lr(FX, tr, ytr)
results["logistic + 2 interactions"] = mod.predict_proba((te[FX]-m_)/s_)[:, 1]
results["gradient boosting"] = HistGradientBoostingClassifier(
    max_iter=200, learning_rate=0.05, max_leaf_nodes=8, min_samples_leaf=60,
    l2_regularization=1.0, random_state=1).fit(tr[F], ytr).predict_proba(te[F])[:, 1]
results["random forest"] = RandomForestClassifier(
    n_estimators=400, min_samples_leaf=25, random_state=1, n_jobs=-1
    ).fit(tr[F], ytr).predict_proba(te[F])[:, 1]

print(f"{'model':30s}{'ROC-AUC':>9}{'PR-AUC':>9}{'lift':>7}{'Brier':>10}")
for nm, p in results.items():
    ap = average_precision_score(yte, p)
    print(f"  {nm:28s}{roc_auc_score(yte,p):>9.4f}{ap:>9.4f}{ap/yte.mean():>7.0f}{brier_score_loss(yte,p):>10.5f}")

model                           ROC-AUC   PR-AUC   lift     Brier
  logistic, additive             0.9554   0.2752     86   0.00269
  logistic + 2 interactions      0.9604   0.2984     93   0.00264
  gradient boosting              0.9565   0.2584     80   0.00272
  random forest                  0.9364   0.2496     78   0.00285

best = FX
mod, m_, s_ = fit_lr(best, tr, ytr)
Xtr = ((tr[best]-m_)/s_).values
Xte = ((te[best]-m_)/s_).values
p_plain = mod.predict_proba(Xte)[:, 1]

rng = np.random.default_rng(3)
extra = rng.choice(np.where(ytr == 1)[0], (ytr == 0).sum() - (ytr == 1).sum(), replace=True)
idx = np.concatenate([np.arange(len(ytr)), extra])
p_over = LogisticRegression(max_iter=3000).fit(Xtr[idx], ytr[idx]).predict_proba(Xte)[:, 1]
p_wt = LogisticRegression(max_iter=3000, class_weight="balanced").fit(Xtr, ytr).predict_proba(Xte)[:, 1]

print(f"{'':32s}{'PR-AUC':>9}{'Brier':>10}{'mean predicted':>16}{'observed':>10}{'ratio':>8}")
for nm, p in [("(a) oversample to 1:1", p_over), ("(b) class weights", p_wt),
              ("(c) untouched, tune threshold", p_plain)]:
    print(f"  {nm:30s}{average_precision_score(yte,p):>9.4f}{brier_score_loss(yte,p):>10.5f}"
          f"{p.mean():>16.4%}{yte.mean():>10.4%}{p.mean()/yte.mean():>7.1f}x")

                                   PR-AUC     Brier  mean predicted  observed   ratio
  (a) oversample to 1:1            0.2985   0.05136        11.5708%   0.3211%   36.0x
  (b) class weights                0.2981   0.05109        11.5307%   0.3211%   35.9x
  (c) untouched, tune threshold    0.2984   0.00264         0.2334%   0.3211%    0.7x

C_FP = 4 + 12          # analyst time plus customer annoyance
amt = te.amount.values

def total_cost(p, t):
    flag = p >= t
    fp = int(((flag) & (yte == 0)).sum()) * C_FP
    fn = float((amt[(~flag) & (yte == 1)] + 25).sum())
    return fp + fn, int(flag.sum()), int(((flag) & (yte == 1)).sum())

print(f"false positive = ${C_FP}    missed fraud = the amount + $25\n")
print(f"{'threshold':>10}{'flagged':>9}{'caught':>8}{'missed':>8}{'cost':>12}")
for t in (0.002, 0.01, 0.05, 0.0883, 0.20, 0.40):
    c, fl, ca = total_cost(p_plain, t)
    print(f"{t:>10.4f}{fl:>9}{ca:>8}{yte.sum()-ca:>8}{c:>12,.0f}")

grid = np.unique(np.round(p_plain, 5))
bc, bt = min((total_cost(p_plain, t)[0], t) for t in grid)
c_none, _, _ = total_cost(p_plain, 2.0)
print(f"\nflag nothing            ${c_none:,.0f}")
print(f"cost-optimal threshold  {bt:.4f}  ->  ${bc:,.0f}   ({1-bc/c_none:.0%} lower)")
fl = int((p_plain >= bt).sum())
print(f"which flags {fl} transactions over 30 days, about {fl/30:.0f} a day")

false positive = $16    missed fraud = the amount + $25

 threshold  flagged  caught  missed        cost
    0.0020     2196      85      10      35,011
    0.0100      887      67      28      16,976
    0.0500      251      44      51      10,402
    0.0883      146      37      58       9,975
    0.2000       58      24      71      11,195
    0.4000       26      17      78      11,609

flag nothing            $14,476
cost-optimal threshold  0.0883  ->  $9,959   (31% lower)
which flags 145 transactions over 30 days, about 5 a day

K = 400 * 30                       # the review team's stated capacity over the test window
thr_cap = np.sort(p_plain)[::-1][K-1]
c_cap, fl_cap, ca_cap = total_cost(p_plain, thr_cap)
print(f"using the full review capacity of {K:,} cases:")
print(f"   threshold {thr_cap:.5f}   flags {fl_cap:,}   catches {ca_cap} of {yte.sum()} frauds")
print(f"   cost ${c_cap:,.0f}   against ${c_none:,.0f} for doing nothing at all")
print(f"\nFilling the queue costs {c_cap/c_none:.0f} times more than having no model.")

using the full review capacity of 12,000 cases:
   threshold 0.00012   flags 12,000   catches 93 of 95 frauds
   cost $190,667   against $14,476 for doing nothing at all

Filling the queue costs 13 times more than having no model.

fig, axes = plt.subplots(1, 2, figsize=(12.6, 4.3))

ax = axes[0]
for nm, p, col, st in [("logistic + interactions", results["logistic + 2 interactions"], VI, "-"),
                       ("gradient boosting", results["gradient boosting"], MUT, "-"),
                       ("random forest", results["random forest"], LT, "--")]:
    pr, rc, _ = precision_recall_curve(yte, p)
    ax.plot(rc, pr, color=col, lw=2.4, ls=st, label=f"{nm}  ({average_precision_score(yte,p):.3f})")
ax.axhline(yte.mean(), color=RD, lw=1.8, ls=":")
ax.text(0.97, yte.mean()+0.012, "a random ranking", ha="right", color=RD, fontsize=9, fontweight="bold")
ax.set_xlabel("recall"); ax.set_ylabel("precision"); ax.set_ylim(0, 0.72)
ax.set_title("Precision-recall, where the base rate is visible"); ax.legend(fontsize=8.5)

ax = axes[1]
for nm, p, col in [("(a) oversampled", p_over, MUT), ("(b) class weights", p_wt, LT),
                   ("(c) untouched", p_plain, VI)]:
    q = pd.qcut(p, 10, labels=False, duplicates="drop")
    t = pd.DataFrame({"p": p, "y": yte, "q": q}).groupby("q").agg(pr=("p","mean"), ob=("y","mean"))
    t = t[t.ob > 0]
    ax.plot(t.pr*100, t.ob*100, "o-", color=col, lw=2.2, ms=6, label=nm)
lims = [0.02, 100]
ax.plot(lims, lims, color=DK, lw=1.6, ls="--")
ax.text(5.0, 9.0, "perfect calibration", color=DK, fontsize=9, fontweight="bold", rotation=40)
ax.set_xscale("log"); ax.set_yscale("log")
ax.set_xlim(*lims); ax.set_ylim(*lims)
ax.set_xlabel("predicted fraud rate (%, log)"); ax.set_ylabel("observed fraud rate (%, log)")
ax.text(1.2, 0.145, "rebalanced: predicts 1 to 60 percent,\nobserves under a tenth of one",
        fontsize=8.8, fontweight="bold", color=RD)
ax.set_title("Same ranking, three different beliefs"); ax.legend(fontsize=9, loc="lower right")
plt.tight_layout(); plt.show()

fig, axes = plt.subplots(1, 2, figsize=(12.6, 4.3))

ax = axes[0]
ts = np.linspace(0.001, 0.5, 400)
cs = [total_cost(p_plain, t)[0] for t in ts]
ax.plot(ts*100, cs, color=VI, lw=2.8)
ax.axhline(c_none, color=RD, lw=2, ls="--")
ax.text(49, c_none+320, "flag nothing", ha="right", color=RD, fontweight="bold", fontsize=9.5)
ax.plot([bt*100], [bc], "o", ms=13, color=GD, mec=DK, mew=1.2, zorder=5)
ax.annotate(f"cost-optimal\n{bt*100:.1f}%, ${bc:,.0f}", xy=(bt*100, bc), xytext=(bt*100+9, bc+2600),
            fontsize=9.5, fontweight="bold", color=DK,
            arrowprops=dict(arrowstyle="->", color=DK, lw=1.5))
ax.set_xlabel("threshold (%)"); ax.set_ylabel("total cost over 30 days ($)")
ax.set_title("The threshold is an accounting question")

ax = axes[1]
labs = ["cost-optimal\nthreshold", "flag\nnothing", "fill the review\ncapacity"]
vals = [bc, c_none, c_cap]
cols = [GD, MUT, RD]
bars = ax.bar(labs, vals, color=cols, edgecolor=DK, lw=0.7, width=0.6)
for r_, v in zip(bars, vals):
    ax.text(r_.get_x()+r_.get_width()/2, v*1.06, f"${v:,.0f}", ha="center", fontweight="bold", fontsize=10)
ax.set_yscale("log"); ax.set_ylabel("total cost over 30 days ($, log scale)")
ax.set_ylim(5e3, 5e5); ax.grid(axis="x", alpha=0)
ax.set_title("Filling the queue is the worst")
plt.tight_layout(); plt.show()

Step 1 · The decision, and what each error costs¶

Step 2 · Cleaning, and the accuracy trap¶

Step 3 · Split by time, because that is how it will be used¶

Step 4 · First look¶

Step 5 · The leak, found before anything else¶

Step 6 · The bake-off, and the baseline that wins¶

Step 7 · Three answers to imbalance¶

Step 8 · The threshold, from the cost matrix¶

Step 9 · The capacity that is not the constraint¶

Step 10 · The pictures¶

Step 11 · Ethics, bias, and limits¶

	transaction_id	transaction_day	amount	hour_of_day	customer_tenure_days	transactions_last_30d	online	new_device	distance_from_home_km	merchant_category_risk
0	T2023984	24	2.31	6	299	25	0	0	12.8	0.041
1	T2017726	18	14.96	20	1354	25	0	1	16.9	0.457
2	T2064569	65	21.39	14	487	33	0	0	596.4	0.227
3	T2002284	3	29.58	5	445	28	1	0	1.6	0.230
4	T2011587	12	44.38	13	2023	16	0	0	4.3	0.158