""" modelsel.py — numerical companion to the glossary entry 'model selection'. Purpose ------- Backs the entry's claims with numbers: the training error falls as the hypothesis space grows while the risk eventually rises again, so the training error cannot select; the validation error can; and a test set that is used to select underestimates the risk of the candidate finally chosen. The weather days are those of the entry's figure (generated exactly as in pythondemos/validation.py), so the entry's numbers 0.47, 2065, 7.06 and 1.83 are reproduced. Self-contained (numpy and matplotlib only), fixed seeds. Setup ----- Six training days and twenty held-back days with feature x = morning minimum temperature and label y = maximum daytime temperature; the twenty are split fourteen to six into validation set and test set. Candidates c = 0, ..., 5 are the hypotheses that ERM with the squared error loss delivers from polynomial regression of degree c, so the hypothesis spaces grow with c (degree 1 is the linear model). Blocks ------ [B-trainerr] The training error falls monotonically with the degree and reaches zero at degree five, where the polynomial passes through all six days; the line's training error is 0.47. [B-valerr] The validation error on the fourteen days falls from degree 0 to degree 1 and rises afterwards: 7.06 for the line, 2065 for the degree-five polynomial. Selecting by validation error picks the line; its average loss on the six test days, 1.83, is reported as the risk estimate. [B-testreuse] Over 2000 draws of the twenty held-back days, the smallest error among the six candidates on a six-day set used for selecting is, on average, far below the error of that same selected candidate on six fresh days: a test set reused for selecting underestimates the risk. Outputs ------- modelsel_degrees.csv : degree, training error, validation error (14 days), test error (6 days) of each candidate. modelsel_testreuse.csv: one row per draw: error of the selected candidate on the selecting set and on fresh days. modelsel.png : matplotlib preview of the two figures (checking only). """ # /// script # requires-python = ">=3.10" # dependencies = ["numpy==2.3.4", "matplotlib==3.10.7"] # /// import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from pathlib import Path OUT_DIR = Path(__file__).parent report = [] def check(name, ok): report.append((name, bool(ok))) print(f" [{'ok' if ok else 'FAIL'}] {name}") def days(n, seed, noise=1.5): """n days: morning minimum x and maximum daytime temperature y (the generator of pythondemos/validation.py, same seeds).""" gen = np.random.default_rng(seed) xs = gen.uniform(-15.0, 5.0, n) return xs, 4.0 + 0.8 * xs + gen.normal(0.0, noise, n) def avg_sqerr(xs, ys, coef): return float(np.mean((ys - np.polyval(coef, xs)) ** 2)) xt, yt = days(6, seed=20260828) # training set xv, yv = days(20, seed=31) # held back xval, yval, xte, yte = xv[:14], yv[:14], xv[14:], yv[14:] DEGREES = range(0, 6) # ------------------------------------------------------------- [B-trainerr] coefs = {c: np.polyfit(xt, yt, c) for c in DEGREES} tr = {c: avg_sqerr(xt, yt, coefs[c]) for c in DEGREES} check("[B-trainerr] training error by degree: " + ", ".join(f"{c}: {tr[c]:.2f}" for c in DEGREES), all(tr[c + 1] <= tr[c] + 1e-12 for c in range(5)) and tr[5] < 1e-8 and abs(tr[1] - 0.47) < 0.005) # --------------------------------------------------------------- [B-valerr] va = {c: avg_sqerr(xval, yval, coefs[c]) for c in DEGREES} te = {c: avg_sqerr(xte, yte, coefs[c]) for c in DEGREES} c_star = min(va, key=va.get) check("[B-valerr] validation error by degree: " + ", ".join(f"{c}: {va[c]:.2f}" for c in DEGREES) + f"; selected degree {c_star}", c_star == 1 and va[0] > va[1] and all(va[c + 1] > va[c] for c in range(1, 5)) and abs(va[1] - 7.06) < 0.005 and abs(round(va[5]) - 2065) < 1) check(f"[B-valerr] average loss of the selected line on the six test " f"days: {te[c_star]:.2f}", abs(te[c_star] - 1.83) < 0.005) # ------------------------------------------------------------ [B-testreuse] rng = np.random.default_rng(0) R = 2000 sel_err, fresh_err = [], [] for _ in range(R): xs = rng.uniform(-15.0, 5.0, 12) ys = 4.0 + 0.8 * xs + rng.normal(0.0, 1.5, 12) x_sel, y_sel, x_new, y_new = xs[:6], ys[:6], xs[6:], ys[6:] errs = {c: avg_sqerr(x_sel, y_sel, coefs[c]) for c in DEGREES} c_sel = min(errs, key=errs.get) sel_err.append(errs[c_sel]) fresh_err.append(avg_sqerr(x_new, y_new, coefs[c_sel])) sel_err, fresh_err = np.array(sel_err), np.array(fresh_err) check(f"[B-testreuse] error of the selected candidate: {sel_err.mean():.2f} " f"on the six days used for selecting, {fresh_err.mean():.2f} on six " f"fresh days (average over {R} draws)", sel_err.mean() < 0.8 * fresh_err.mean()) # ---------------------------------------------------------------- CSV with open(OUT_DIR / "modelsel_degrees.csv", "w") as fh: fh.write("degree,trainerr,valerr,testerr\n") for c in DEGREES: fh.write(f"{c},{tr[c]:.4f},{va[c]:.4f},{te[c]:.4f}\n") with open(OUT_DIR / "modelsel_testreuse.csv", "w") as fh: fh.write("sel_err,fresh_err\n") for a, b in zip(sel_err, fresh_err): fh.write(f"{a:.4f},{b:.4f}\n") # -------------------------------------------------------------- preview fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9.4, 3.9)) ds = list(DEGREES) ax1.semilogy(ds, [tr[c] + 1e-3 for c in ds], "ko-", label="training error (6 days)") ax1.semilogy(ds, [va[c] for c in ds], "ks--", mfc="none", label="validation error (14 days)") ax1.semilogy(ds, [te[c] for c in ds], "k^:", mfc="none", label="test error (6 days)") ax1.set_xlabel("degree of the polynomial (size of the hypothesis space)") ax1.set_ylabel("average squared error loss") ax1.set_title("training error falls, validation error rises again") ax1.legend(frameon=False, fontsize=8) bins = np.linspace(0.0, 8.0, 33) ax2.hist(sel_err, bins=bins, histtype="step", color="black", ls="-", label="on the six days used for selecting") ax2.hist(fresh_err, bins=bins, histtype="step", color="black", ls="--", label="on six fresh days") ax2.set_xlabel("average squared error loss of the selected candidate") ax2.set_ylabel("number of draws") ax2.set_title("a test set reused for selecting reads too low") ax2.legend(frameon=False, fontsize=8) fig.tight_layout() fig.savefig(OUT_DIR / "modelsel.png", dpi=110) n_ok = sum(ok for _, ok in report) print(f"\n{n_ok}/{len(report)} checks pass") print(f"wrote {OUT_DIR / 'modelsel_degrees.csv'}, {OUT_DIR / 'modelsel_testreuse.csv'}, " f"{OUT_DIR / 'modelsel.png'}") if n_ok != len(report): raise SystemExit(1)