""" Tutorial 95 - Quantile regression from scratch: the whole margin, not the average. Least squares fits one line through the middle of a cloud and calls it the answer. That line is a conditional MEAN, and a mean says nothing about how wide the cloud is or whether it gets wider as you move along it. Quantile regression fits a line to any percentile you name, so five fits describe the whole conditional distribution instead of its centre. The test case: NFL point spreads. A spread is a forecast of one number, and the usual conversion to a win probability assumes the margin's spread-around-the-line never changes - only its centre moves. This script checks that assumption instead of making it, on every played game in the bundled file. Written by hand: the check (pinball) loss, the iteratively reweighted least squares solver, a brute-force exact solver used to prove the fast one right, and a pairs bootstrap for the slopes. No scipy, no statsmodels. Run: python downloads/95_quantile_regression_from_scratch.py Data: nflverse games table (github.com/nflverse/nfldata), June 2026 snapshot, bundled as nfl_games_lines.csv. """ import math import os import matplotlib.pyplot as plt import numpy as np import pandas as pd import sdt_common as sdt sdt.init("quantile-regression-from-scratch") HERE = os.path.dirname(os.path.abspath(__file__)) CSV = os.path.join(HERE, "nfl_games_lines.csv") TAUS = (0.10, 0.25, 0.50, 0.75, 0.90) # --- the method, written out in full ------------------------------------------------- def check_loss(r, tau): """Koenker-Bassett check loss (the 'pinball'): tau*r above the line, (tau-1)*r below. Under-predicting by a point costs tau; over-predicting by a point costs 1-tau. At tau = 0.5 the two are equal and this is absolute error, whose minimiser is the median. Tilt the penalty and the minimiser slides to another quantile. """ r = np.asarray(r, dtype=float) return float(np.sum(np.where(r > 0, tau * r, (tau - 1) * r))) def quantile_of(v, tau): """The constant that minimises the check loss: an order statistic, by hand.""" s = np.sort(np.asarray(v, dtype=float)) k = min(max(int(math.ceil(tau * len(s))) - 1, 0), len(s) - 1) return float(s[k]) def qreg(X, y, tau, iters=500, eps=1e-9): """Quantile regression by iteratively reweighted least squares. The check loss is |r| scaled by tau or 1-tau, and |r| = r^2 / |r|, so a weighted least-squares step with weight (tau or 1-tau)/|r| solves it. Iterate: the weights are recomputed from the residuals until the coefficients stop moving. eps floors the denominator so a point sitting exactly on the line cannot divide by zero. """ X, y = np.asarray(X, dtype=float), np.asarray(y, dtype=float) b = np.linalg.lstsq(X, y, rcond=None)[0] for _ in range(iters): r = y - X @ b w = np.where(r > 0, tau, 1 - tau) / np.maximum(np.abs(r), eps) XW = X * w[:, None] nb = np.linalg.solve(X.T @ XW, XW.T @ y) if np.max(np.abs(nb - b)) < 1e-12: return nb b = nb return b def qreg_brute(x, y, tau): """Exact simple quantile regression: the optimum interpolates two data points. Slow (every pair), but it needs no numerical argument at all, so it is the thing to check the fast solver against. """ best = (np.inf, None, None) for i in range(len(x)): for j in range(i + 1, len(x)): if x[i] == x[j]: continue b = (y[i] - y[j]) / (x[i] - x[j]) a = y[i] - b * x[i] loss = check_loss(y - (a + b * x), tau) if loss < best[0] - 1e-12: best = (loss, a, b) return best # --- 1. the mean line, and the question it cannot answer ------------------------------- games = pd.read_csv(CSV) played = games.dropna(subset=["result", "spread_line"]).copy() x = played.spread_line.to_numpy(float) y = played.result.to_numpy(float) X = np.column_stack([np.ones(len(y)), x]) b_ols = np.linalg.lstsq(X, y, rcond=None)[0] ols_resid = y - X @ b_ols r2 = 1 - ols_resid.var() / y.var() with sdt.snippet("ols"): print(f"{len(games)} rows in the file, {len(played)} played games with a line " f"({played.season.min()}-{played.season.max()})") print(f" ties kept: {int((y == 0).sum())} pick'ems (spread 0): {int((x == 0).sum())}") print(f"\nleast squares: margin = {b_ols[0]:+.4f} {b_ols[1]:+.4f} x spread") print(f" R-squared {r2:.4f}, residual sd {ols_resid.std(ddof=2):.4f} points") print(f"\nmargin overall: mean {y.mean():.4f}, median {np.median(y):.1f}") print("the mean line is one line. It says nothing about how far from it a game lands.") assert len(games) == 7548 and len(played) == 7276 assert int((y == 0).sum()) == 15 and int((x == 0).sum()) == 31 assert played.season.min() == 1999 and played.season.max() == 2025 assert played.result.notna().all() and played.spread_line.notna().all() assert abs(b_ols[0] - (-0.0031)) < 5e-4 and abs(b_ols[1] - 1.0430) < 5e-4 assert 0.1810 < r2 < 0.1818 assert 13.200 < ols_resid.std(ddof=2) < 13.205 assert 2.344 < y.mean() < 2.345 and np.median(y) == 3.0 assert abs(np.median(y) - y.mean()) > 0.65 # the mean and the median disagree assert (y == np.round(y)).all() # margins are whole points assert set(np.unique(x * 2)) <= set(np.arange(-40, 56).astype(float)) # lines move by halves # --- 2. the check loss, and why its minimiser is a quantile --------------------------- wc = played[(played.season == 2025) & (played.game_type == "WC")] slate = np.sort(wc.result.to_numpy(float)) TAU_EX = 0.75 cands = list(slate) losses = [check_loss(slate - c, TAU_EX) for c in cands] k_rule = int(math.ceil(TAU_EX * len(slate))) best_c = cands[int(np.argmin(losses))] with sdt.snippet("check"): print(f"the {len(slate)} wild-card games of the 2025 season, home margins sorted:") print(" ", [int(v) for v in slate]) print(f"\ncheck loss at tau = {TAU_EX} for each candidate:") for c, L in zip(cands, losses): print(f" a = {c:+6.1f} loss = {L:8.2f}{' <- smallest' if L == min(losses) else ''}") print(f"\nrule: the minimiser is order statistic ceil(tau*n) = " f"ceil({TAU_EX} x {len(slate)}) = {k_rule}, which is {slate[k_rule - 1]:+.0f}") print("\nsame thing on all 7,276 games - IRLS with only an intercept vs the order statistic:") for t in TAUS: a_irls = float(qreg(np.ones((len(y), 1)), y, t)[0]) a_exact = quantile_of(y, t) print(f" tau {t:.2f} IRLS {a_irls:8.4f} order statistic {a_exact:+6.1f} " f"loss {check_loss(y - a_irls, t):10.1f} vs {check_loss(y - a_exact, t):10.1f}") assert len(slate) == 6 and list(slate.astype(int)) == [-24, -4, -3, -3, 4, 13] assert (wc.game_type == "WC").all() and len(wc) == 6 assert best_c == 4.0 and slate[k_rule - 1] == 4.0 and k_rule == 5 assert abs(min(losses) - 19.25) < 1e-9 # hand-checkable: 0.75*9 + 0.25*50 assert abs(check_loss(slate - 4.0, 0.75) - (0.75 * 9 + 0.25 * 50)) < 1e-9 assert abs(check_loss(slate - (-3.0), 0.75) - 22.75) < 1e-9 assert abs(check_loss(slate - 13.0, 0.75) - 23.75) < 1e-9 assert check_loss(slate - 4.5, 0.75) > min(losses) # between data points is worse assert abs(check_loss(slate - np.median(slate), 0.5) - 22.5) < 1e-9 for _t in TAUS: # IRLS never beats the order statistic, and never trails it much _a = float(qreg(np.ones((len(y), 1)), y, _t)[0]) _best = check_loss(y - quantile_of(y, _t), _t) _gap = check_loss(y - _a, _t) - _best assert -1e-9 <= _gap < 0.25 and _gap / _best < 1e-4 assert quantile_of(y, 0.5) == 3.0 and quantile_of(y, 0.9) == 21.0 assert quantile_of(y, 0.1) == -17.0 and quantile_of(y, 0.25) == -7.0 assert quantile_of(y, 0.75) == 11.0 # tau = 0.10 is the one where IRLS stops short of the order statistic: the loss has a # flat stretch there, so both answers are optimal and neither is "the" minimiser. assert abs(float(qreg(np.ones((len(y), 1)), y, 0.10)[0]) - (-17.0)) > 1e-3 # --- 3. five lines, and a proof the solver is right ------------------------------------ FIT = {t: qreg(X, y, t) for t in TAUS} on_line = {t: int((np.abs(y - X @ FIT[t]) < 1e-9).sum()) for t in TAUS} R1 = {t: 1 - check_loss(y - X @ FIT[t], t) / check_loss(y - quantile_of(y, t), t) for t in TAUS} s25 = played[played.season == 2025] x25, y25 = s25.spread_line.to_numpy(float), s25.result.to_numpy(float) X25 = np.column_stack([np.ones(len(x25)), x25]) with sdt.snippet("fit"): print("quantile regression on 7,276 games: margin = a + b x spread\n") print(" tau a b games exactly on the line R1") for t in TAUS: a, b = FIT[t] print(f" {t:.2f} {a:8.4f} {b:7.4f} {on_line[t]:>10} {R1[t]:.4f}") med_hits = played[np.abs(y - X @ FIT[0.50]) < 1e-9] pairs = sorted({(float(r.spread_line), int(r.result)) for r in med_hits.itertuples()}) print(f"\nthe median line runs through {on_line[0.50]} real games, at just " f"{len(pairs)} distinct (spread, margin) pairs:") print(" ", pairs) print(f"\nsame five fits on the {len(x25)} games of 2025 alone, IRLS vs brute force " f"over all {len(x25) * (len(x25) - 1) // 2} pairs of games:") for t in TAUS: bi = qreg(X25, y25, t) Lb, ab, bb = qreg_brute(x25, y25, t) Li = check_loss(y25 - (bi[0] + bi[1] * x25), t) print(f" tau {t:.2f} IRLS a={bi[0]:8.3f} b={bi[1]:6.4f} " f"brute a={ab:8.3f} b={bb:6.4f} loss {Li:9.4f} vs {Lb:9.4f}" f"{'' if abs(bi[1] - bb) < 1e-6 else ' <- a vertex short'}") assert abs(FIT[0.10][0] - (-559 / 33)) < 1e-6 and abs(FIT[0.10][1] - 38 / 33) < 1e-6 assert abs(FIT[0.25][0] - (-227 / 27)) < 1e-6 and abs(FIT[0.25][1] - 28 / 27) < 1e-6 assert abs(FIT[0.50][0] - (-1 / 7)) < 1e-6 and abs(FIT[0.50][1] - 20 / 21) < 1e-6 assert abs(FIT[0.75][0] - (121 / 15)) < 1e-6 and abs(FIT[0.75][1] - 16 / 15) < 1e-6 assert abs(FIT[0.90][0] - (343 / 20)) < 1e-6 and abs(FIT[0.90][1] - 11 / 10) < 1e-6 assert all(on_line[t] >= 2 for t in TAUS) # the optimum interpolates data points assert on_line[0.50] == 51 and on_line[0.10] == 6 and on_line[0.90] == 9 assert pairs == [(-3.0, -3), (7.5, 7)] # 51 games, two distinct coordinates assert 0.090 < min(R1.values()) and max(R1.values()) < 0.108 assert abs(R1[0.50] - 0.0908) < 5e-4 and abs(R1[0.10] - 0.1077) < 5e-4 for t in TAUS: # the fit beats the flat line it nests assert check_loss(y - X @ FIT[t], t) < check_loss(y - quantile_of(y, t), t) _exact = 0 for t in TAUS: # brute force is optimal by construction; IRLS tracks it _b, (_L, _a, _bb) = qreg(X25, y25, t), qreg_brute(x25, y25, t) _Li = check_loss(y25 - (_b[0] + _b[1] * x25), t) assert _Li >= _L - 1e-9 # nothing beats the exact optimum assert (_Li - _L) / _L < 1e-4 # and IRLS is never materially off it _exact += int(abs(_b[1] - _bb) < 1e-6) # the intercept can sit in a flat stretch assert _exact == 4 # at tau 0.90 it stops one vertex short assert len(x25) == 285 # --- 4. read the fan ------------------------------------------------------------------- def band(s, lo=0.10, hi=0.90): return (FIT[hi][0] + FIT[hi][1] * s) - (FIT[lo][0] + FIT[lo][1] * s) grid = np.linspace(x.min(), x.max(), 400) fan = np.array([FIT[t][0] + FIT[t][1] * grid for t in TAUS]) crossings = int((np.diff(fan, axis=0) <= 0).sum()) rng = np.random.default_rng(95) BOOT = 200 ci = {} for t in TAUS: draws = np.empty(BOOT) for i in range(BOOT): idx = rng.integers(0, len(y), len(y)) draws[i] = qreg(np.column_stack([np.ones(len(idx)), x[idx]]), y[idx], t)[1] ci[t] = (float(np.percentile(draws, 2.5)), float(np.percentile(draws, 97.5))) played["resid"] = played.result - played.spread_line played["band"] = pd.cut(played.spread_line.abs(), [-0.01, 2.5, 6.5, 10.5, 30], labels=["0-2.5", "3-6.5", "7-10.5", "11+"]) buckets = {} for lab, d in played.groupby("band", observed=True): q = np.percentile(d.resid, [10, 50, 90]) buckets[lab] = (len(d), q[0], q[1], q[2], q[2] - q[0]) # every whole-number line posted at least 150 times - the hits AND the misses, so the # reader sees the selection rather than a cherry-picked five. KEY = [-7.0, -3.0, -1.0, 1.0, 3.0, 4.0, 6.0, 7.0, 10.0] key_rows = [] for s in KEY: d = played[played.spread_line == s].result key_rows.append((s, len(d), float(d.median()), float(d.mean()))) with sdt.snippet("fan"): print("where the five lines sit, by spread:\n") print(" spread q10 q25 q50 q75 q90 80% band") for s in (-10, -3, 0, 3, 7, 14): vals = [FIT[t][0] + FIT[t][1] * s for t in TAUS] print(f" {s:+5.1f} {vals[0]:6.2f} {vals[1]:6.2f} {vals[2]:6.2f} " f"{vals[3]:6.2f} {vals[4]:6.2f} {band(s):6.2f}") print(f"\nthe 80% band changes by {FIT[0.90][1] - FIT[0.10][1]:+.4f} points per point of " f"spread: {band(-10):.2f} wide at a 10-point road favourite, {band(14):.2f} at a " f"14-point home favourite.") print(f"lines crossing anywhere in the observed range ({x.min():.0f} to {x.max():.0f}): " f"{crossings}") print(f"\nslopes with a 95% pairs-bootstrap interval ({BOOT} resamples):") for t in TAUS: lo, hi = ci[t] # an endpoint landing exactly on 1 is the discreteness of these fits, not a verdict reaches_one = lo <= 1 + 1e-9 and hi >= 1 - 1e-4 print(f" tau {t:.2f} b = {FIT[t][1]:.4f} [{lo:.4f}, {hi:.4f}]" f"{'' if reaches_one else ' <- clear of 1'}") with sdt.snippet("nomodel"): print("the same question with no model at all - margin minus spread, by size of line:\n") print(" |line| n p10 p50 p90 width") for lab, (n_b, p10, p50, p90, w) in buckets.items(): print(f" {lab:<8} {n_b:>5} {p10:6.1f} {p50:6.1f} {p90:6.1f} {w:6.1f}") print("\nwhere the median lands at every whole-number line posted 150+ times:\n") print(" spread n median mean median - spread") for s, n_s, med, mean in key_rows: hit = "exact" if abs(med - s) < 1e-9 else f"{med - s:+.0f}" print(f" {s:+5.1f} {n_s:>5} {med:+6.1f} {mean:+6.2f} {hit}") n_hit = sum(1 for s, _n, med, _m in key_rows if abs(med - s) < 1e-9) print(f"\n{n_hit} of {len(key_rows)} land exactly on the line - and they are the key " f"numbers football scores in.") assert crossings == 0 # no quantile crossing in range assert abs(FIT[0.90][1] - FIT[0.10][1] - (-0.0515)) < 5e-4 assert 33.3 < band(14) < 33.45 and 34.55 < band(-10) < 34.7 assert abs(band(-10) - band(14)) < 1.3 # 1.2 points of 34: flat assert all(abs(FIT[t][1] - 1) < 0.16 for t in TAUS) # every slope within 0.16 of one assert max(FIT[t][1] for t in TAUS) == FIT[0.10][1] assert 0.94 < FIT[0.50][1] < 0.96 assert ci[0.10][0] > 1 and ci[0.25][0] < 1 < ci[0.25][1] assert ci[0.75][1] > 1 and ci[0.90][1] > 1 assert all(ci[t][0] <= 1 + 1e-9 and ci[t][1] >= 1 - 1e-4 for t in (0.25, 0.50, 0.75, 0.90)) # only tau 0.10 is clear of slope 1 assert abs((FIT[0.75][1] - FIT[0.25][1]) - 0.0296) < 5e-4 # the IQR widens as the band narrows assert (FIT[0.75][1] - FIT[0.25][1]) * (FIT[0.90][1] - FIT[0.10][1]) < 0 # opposite signs assert all(hi - lo < 0.22 for lo, hi in ci.values()) # intervals this wide cannot resolve 0.05 assert buckets["0-2.5"][0] == 1481 and buckets["3-6.5"][0] == 3593 assert buckets["7-10.5"][0] == 1648 and buckets["11+"][0] == 554 assert sum(b[0] for b in buckets.values()) == len(played) assert all(33.4 <= b[4] <= 34.1 for b in buckets.values()) # every bucket's 80% band assert abs(buckets["0-2.5"][4] - buckets["11+"][4]) <= 0.5 assert [r[1] for r in key_rows] == [158, 498, 261, 285, 655, 227, 237, 333, 166] hits = [s for s, _n, med, _m in key_rows if abs(med - s) < 1e-9] assert hits == [-7.0, -3.0, 3.0, 7.0, 10.0] # the key numbers, and only those assert [med for s, _n, med, _m in key_rows] == [-7, -3, -2, -1, 3, 5, 5, 7, 10] assert abs(key_rows[7][3] - 8.90) < 5e-3 and key_rows[7][2] == 7.0 # at +7 the mean is not assert abs(key_rows[0][3] - (-8.55)) < 5e-3 assert all(abs(med - s) <= 2 for s, _n, med, _m in key_rows) # misses are small # --- 5. does the fan hold where it was not fitted? -------------------------------------- train = played[played.season <= 2019] test = played[played.season >= 2020] xt, yt = train.spread_line.to_numpy(float), train.result.to_numpy(float) xe, ye = test.spread_line.to_numpy(float), test.result.to_numpy(float) Xt = np.column_stack([np.ones(len(xt)), xt]) FIT_TR = {t: qreg(Xt, yt, t) for t in TAUS} sd_old = (train.result - train.spread_line).std() sd_new = (test.result - test.spread_line).std() below = {t: float((ye < FIT_TR[t][0] + FIT_TR[t][1] * xe).mean()) for t in TAUS} inside80 = float(((ye >= FIT_TR[0.10][0] + FIT_TR[0.10][1] * xe) & (ye <= FIT_TR[0.90][0] + FIT_TR[0.90][1] * xe)).mean()) inside50 = float(((ye >= FIT_TR[0.25][0] + FIT_TR[0.25][1] * xe) & (ye <= FIT_TR[0.75][0] + FIT_TR[0.75][1] * xe)).mean()) with sdt.snippet("holdout"): print(f"fit on {len(xt)} games (1999-2019), scored on {len(xe)} games (2020-2025):\n") print(" tau fitted line share of held-out games below it") for t in TAUS: a, b = FIT_TR[t] print(f" {t:.2f} margin = {a:+8.4f} {b:+.4f} x spread {below[t]:.4f}") print(f"\n inside the 80% band: {inside80:.4f} inside the 50% band: {inside50:.4f}") print(f" worst miss on a nominal level: " f"{max(abs(below[t] - t) for t in TAUS):.4f}") print(f"\n sd of (margin - spread): {sd_old:.4f} in 1999-2019, {sd_new:.4f} in 2020-2025") assert len(xt) == 5583 and len(xe) == 1693 and len(xt) + len(xe) == len(played) assert abs(below[0.50] - 0.4903) < 5e-4 and abs(below[0.25] - 0.2310) < 5e-4 assert abs(below[0.10] - 0.0780) < 5e-4 and abs(below[0.90] - 0.9067) < 5e-4 assert abs(below[0.75] - 0.7513) < 5e-4 assert max(abs(below[t] - t) for t in TAUS) < 0.025 # every level within 2.5 points assert 0.828 < inside80 < 0.829 and 0.520 < inside50 < 0.521 assert inside80 > 0.80 and inside50 > 0.50 # the tails are over-covered, not under assert below[0.10] < 0.10 and below[0.90] > 0.90 # both tails pulled in since 2020 assert sd_new < sd_old # the modern era really is tighter assert 13.363 < sd_old < 13.365 and 12.662 < sd_new < 12.663 # --- 6. the exhibit --------------------------------------------------------------------- BROWN, GREY = sdt.sport_color("football"), "#8A8577" SHADE = ["#C9BFAC", "#9C8F76", BROWN, "#9C8F76", "#C9BFAC"] fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10.6, 4.9), gridspec_kw={"width_ratios": [1.45, 1]}) ax1.scatter(x + np.random.default_rng(0).uniform(-0.18, 0.18, len(x)), y, s=4, color="#6C7079", alpha=0.10, linewidths=0, zorder=1) xs = np.linspace(-19, 27, 100) for t, col in zip(TAUS, SHADE): ax1.plot(xs, FIT[t][0] + FIT[t][1] * xs, color=col, lw=2.0 if t == 0.50 else 1.5, zorder=3) ax1.annotate(f" {t:.2f}", (27, FIT[t][0] + FIT[t][1] * 27), fontsize=8.5, color=col, va="center", annotation_clip=False) ax1.plot(xs, b_ols[0] + b_ols[1] * xs, color="#2C5E8A", lw=1.4, ls="--", zorder=4, label="least squares (the mean)") ax1.set_xlim(-21, 30) ax1.set_ylim(-45, 50) ax1.set_xlabel("point spread (positive = home favoured)") ax1.set_ylabel("final margin, home minus away") ax1.set_title("Five quantiles of the margin", fontsize=11.5) ax1.legend(loc="upper left", fontsize=8.2, frameon=False) ax1.grid(True, axis="both") lo_err = [FIT[t][1] - ci[t][0] for t in TAUS] hi_err = [ci[t][1] - FIT[t][1] for t in TAUS] ax2.errorbar([t for t in TAUS], [FIT[t][1] for t in TAUS], yerr=[lo_err, hi_err], fmt="o", color=BROWN, ecolor=GREY, capsize=4, markersize=6, lw=1.4) ax2.axhline(1.0, color="#2C5E8A", lw=1.2, ls="--") # park the note in the empty quadrant: every interval sits above 1 out here ax2.text(0.59, 0.93, "slope 1 = the spread\nmoves the whole\ndistribution, one for one", fontsize=8.2, color="#2C5E8A", ha="left", va="center") ax2.set_xlim(0.03, 0.97) ax2.set_ylim(0.86, 1.26) ax2.set_xticks(list(TAUS)) ax2.set_xlabel("quantile (tau)") ax2.set_ylabel("fitted slope on the spread") ax2.set_title("Near-parallel, with wide intervals", fontsize=11.5) ax2.grid(True, axis="both") fig.tight_layout() sdt.save_fig(fig, "margin_fan", source="nflverse games table (github.com/nflverse/nfldata)", asof="June 2026") print("\nall asserts passed")