Statistics: trusting your numbers

You already know how to compute a number. The hard part of production is knowing whether the number means anything: whether 0.93 really beats 0.91, whether last week's dip was real, whether the model you shipped is better than the one it replaced. That is a statistics problem, and this module is the working kit for it.

Cold warm-up

1 · Why this module is the one that pays the salary

Ask what a data team actually delivers and the honest answer is not models, it is claims: this ranking is better, this change raised retention, this anomaly is real. Each claim is a comparison between two noisy numbers. Statistics is the tooling that turns "the number went up" into "the number went up because of the change, and here is how sure we are".

It protects you from the two expensive mistakes. Shipping noise: you A/B test, see +0.4%, ship, and the next cohort reverts — because you never computed an interval and the interval would have contained zero. Missing real effects: you run a test on 500 users, see $p = 0.12$, conclude "no difference", and lock in an inferior design for a year — the test had no power to see the effect that was there. Both failures are invisible without this material.

2 · Population vs sample, statistic vs parameter

The population is everything you want to generalise to — all users, all requests, all future inputs. The sample is the finite pile of observations you actually have. Numbers computed from the population are parameters; numbers computed from the sample are statistics. You never see a parameter, so you use a statistic as its estimator.

Statistics are only trustworthy if the sample is a random draw from the population — roughly, every member had a known, non-zero chance of being included, and inclusion did not depend on the answer. Also distinguish the sampling frame (the list you drew from: opted-in users, logged events, US phone numbers) from the population you claim to describe (all users, all customers, all citizens). Frames are rarely the population, and the gap is a bias no sample size fixes.

Without random sampling you get selection bias — the sample is systematically different from the target. Survivorship bias in logged data is the ML-flavoured version: your training data is a sample of the world, and often a biased one. Clicks only exist for items that were shown; labels only exist for cases a human reviewed; purchases only exist for users who were not churned last month. The model learns the policy that generated the log, not the world.

3 · Estimators and their quality

An estimator $\hat\theta$ is a function of the sample that guesses a parameter $\theta$. Many functions could do the job — $\bar x$ estimates $\mu$, but so does the sample median, or the midpoint of the largest and smallest values. So estimators compete on quality, and there are four standard axes.

1. Unbiasedness — the estimator is centred on the truth on average: $\mathbb{E}[\hat\mu] = \mu$. Over many samples it does not systematically err high or low. (Unbiasedness is about the average over samples, not about any single answer.)
2. Variance — how much the estimator jitters between samples. Unbiased but wildly noisy is not useful.
3. Consistency — as $n \to \infty$ the estimator converges to the parameter. The sample mean is consistent; "the largest observed value" generally is not.
4. Efficiency — among unbiased estimators, who has the smallest variance.

They are combined into one score, the mean squared error of the estimator — and it splits exactly like the model decomposition you will meet in section 9:

Why $n-1$. The variance estimator $\frac1n\sum(X_i-\bar X)^2$ is biased because $\bar X$ was fitted to the same data: it is closer to each $X_i$ than the true $\mu$ is, so the squared deviations come out too small. Two lines:

Taking expectations, with $\mathbb{E}[\sum(X_i-\mu)^2] = n\sigma^2$ and $\mathbb{E}[(\bar X-\mu)^2] = \sigma^2/n$, gives the promised answer:

Maximum likelihood — Module 04's estimator — is (asymptotically) efficient: in the large-sample limit the MLE attains the smallest variance any unbiased estimator can have, the Cramér–Rao bound. Stated here, not proved: the practical rule is that MLE is the default you reach for, and small-sample corrections are the exception.

4 · Sampling distributions, CLT, standard error, confidence intervals

$\bar x$ is itself a random variable — redraw your sample and you get a different $\bar x$. The sampling distribution is the distribution of that variable, and the central limit theorem says that for large $n$ it is approximately normal regardless of the shape of the data:

Because $\sigma$ is unknown we estimate it with $s$, and the standard deviation of an estimator is called the standard error:

The 95% confidence interval is then the point estimate plus and minus about two standard errors:

The $\sqrt n$ controls how fast you get precision: quadruple the data and the interval halves; multiply the data by 100 and you gain exactly one decimal place. There is no way around that arithmetic — it is why "we need 100× the traffic to resolve a 0.01 difference" is a scheduling statement, not pessimism.

5 · Hypothesis testing, done properly

A test is a decision procedure in five steps, and it is worth being able to recite them:

For two conversion rates the statistic is the pooled z-score — the observed difference divided by how big a difference pure chance would produce under $H_0$:

The p-value is the probability, assuming $H_0$ is true, of observing data at least this extreme. If $p \lt \alpha = 0.05$, reject $H_0$. That is the machinery. Now the four corrections that separate a professional report from a naive one:

Significance ≠ practical significance. With $n = 1{,}000{,}000$ observations of a rate near 50%, $\mathrm{SE} = \sqrt{0.5\cdot0.5/10^6} = 0.05$ percentage points, so a difference of 0.1 points gives $z = 2$ and $p \approx 0.046$: statistically significant, worth nothing. As $n$ grows the test will find something for every trivial difference — the p-value answers "is it exactly zero?", and the honest answer for almost any real system is yes, it is not exactly zero.

Absence of evidence is not evidence of absence. A non-significant result means "not enough data to distinguish this from noise", which is the power problem: $1-\beta$ is the probability of detecting a real effect of a given size. Power rises with effect size and with $n$, and the sample size you need scales as the inverse square of the effect:

Multiple comparisons / p-hacking. Test 20 features at $\alpha = 0.05$ and you expect one false positive by construction; the chance of at least one is $1 - 0.95^{20} \approx 64\%$. Two standard fixes: Bonferroni divides the threshold by the number of tests ($0.05/20 = 0.0025$), which controls the chance of any false positive — conservative and simple. False discovery rate control (Benjamini–Hochberg) instead limits the share of your discoveries that are false: sort the p-values $p_{(1)} \le \dots \le p_{(m)}$, find the largest $k$ with $p_{(k)} \le k\alpha/m$, and declare the first $k$ real — better when you are screening many candidates and can tolerate a few misses.

6 · A/B testing, as it is actually done

An A/B test is a randomised experiment, and randomisation is the whole trick. It is the one mechanism that makes treatment and control exchangeable: everything else about the users — motivation, device, timezone, latent churn risk — is, on average, split evenly between arms, so any remaining difference in the outcome is attributable to the treatment rather than to who ended up in which arm. Observational comparisons have to model their way to that property; randomisation gets it by construction.

• Unit of randomisation = the unit that experiences the treatment. Randomise by user, not by request: the same person on both arms gets overlapping, self-inconsistent experiences (interference), and their requests are not independent, which quietly invalidates the standard error. One user, one arm, hashed from a stable id.
• Metric hierarchy. One primary metric that decides the ship, defined before launch. Secondary metrics that explain why. Guardrail metrics that can veto — latency, error rate, unsubscribe rate. A win on a secondary metric with a guardrail breached is not a win.
• Sample size from the MDE. Decide the smallest effect worth detecting ($\delta$), then size the test for it at 80% power: for a baseline near 8.5% and $\delta = 1$ percentage point that is about 12,200 users per arm; for $\delta = 0.1$ point, roughly 100× that. If you cannot afford that, your test cannot answer the question you are asking.
• Stop at the pre-set horizon. Peeking daily and stopping at the first $p \lt 0.05$ inflates false positives — in Block 2's simulation, 40 daily looks turn a 5% false-positive rate into 30%. Fixes: a fixed horizon (simplest), or sequential testing — SPRT / always-valid p-values — which spends the error budget as you look and lets you stop early legitimately.
• Variance reduction. CUPED regresses the metric on each user's pre-experiment value and uses the residual, which can cut variance (and therefore sample size) substantially on metrics with strong user-level autocorrelation — same idea as conditioning, one sentence, big win.

Common failure modes, in the order they bite: novelty effects (users react to what is new, the lift decays over days — plot the effect by cohort day, not just pooled); seasonality (a test spanning a weekend or a sale is confounded with the calendar — run whole cycles); Simpson's paradox in segments (the aggregate wins while every segment loses, or the reverse — section 7 builds the numbers, and the cause is always the same: the mixture of segments differs between arms, so the pooled number is a confounded average); and sample ratio mismatch (if the split is not what you asked for, something upstream is broken and the result is void).

7 · Correlation, causation, and the difference that pays

Pearson correlation measures the strength of a linear relationship between two variables:

Its limitations are worth memorising because they are the ones people trip over: it sees only linear structure (a perfect parabola gives $r \approx 0$), it is extremely sensitive to outliers (one wild point can invent or destroy a correlation), and it is symmetric — $r(x,y) = r(y,x)$, so it can never tell you which way an arrow points. The only road from correlation to causation is design, not more correlation.

A confounder, concretely. Ice-cream sales and drowning deaths correlate; ice cream causes neither — summer (temperature) drives both. The same shape appears in product data: premium users convert better, but is it the premium badge or because motivated users buy premium? A confounder is a third variable that influences both the treatment and the outcome, so the treatment and outcome move together even if there is no arrow between them.

Simpson's paradox, worked. A hospital compares a new treatment with the old one, patients arrive at differing severity, and severity is not randomised:

The treatment wins in both subgroups and loses overall. The reason is a confounded mixture: the control arm is 97% mild cases (70 of 72) while the treatment arm is only 71% mild (20 of 28), and mild cases recover far more often. Pooling mixes the subgroups with different weights, and the weights — not the treatment — produce the reversal. The fix is to compare within strata or reweight both arms to a common mixture.

Randomised experiments are the causal gold standard because randomisation breaks every backdoor path at once. Observational studies must assume their way there, and the two assumptions have names: ignorability (conditional on the covariates you measured, there is no unmeasured confounder — every backdoor path from treatment to outcome is blocked) and positivity (every type of unit has a non-zero chance of each treatment; if some subgroup never receives the treatment, its counterfactual is unidentified). Drawn as a DAG in words: $Z \to T$, $Z \to Y$, $T \to Y$ — the path $T \leftarrow Z \to Y$ is a backdoor, and conditioning on $Z$ blocks it; conditioning on a collider (a variable caused by both $T$ and $Y$) does the opposite and opens a path that was closed.

8 · Linear regression as a statistical model, not just an algorithm

Scikit-learn sells linear regression as a predictor. The statistical version is a generative claim about how the data were produced:

Take that model seriously and you discover that ordinary least squares is exactly maximum likelihood under Gaussian noise (Module 04) — minimising $\sum_i (y_i - x_i^\top\beta)^2$ is maximising the log-likelihood. What the statistical framing adds is the ability to say how uncertain the coefficients are, which is what turns a fitted line into a claim.

Residual analysis is the diagnostic. Plot residuals against fitted values: a funnel shape means heteroscedasticity, a U-shape means you left out nonlinearity, runs and clumps mean dependence. A model with a bad $R^2$ can be honest; a model with a beautiful $R^2$ and structured residuals is lying about its form — Block 4 fits a line to a parabola and the residuals correlate $+0.94$ with $x^2$ while the correct quadratic gives $+0.00$.

$R^2$ vs adjusted $R^2$. $R^2$ is the share of variance explained and it can only go up when you add a feature — even a useless one absorbs noise. The adjusted version charges rent per parameter:

Two cautions to keep with it. Extrapolation: a model is only trustworthy inside the support of its training data — fitting a line to housing prices in the 100k–500k range says nothing about 2M, and no amount of fit-quality inside the range changes that. Multicollinearity: highly correlated features do not damage predictions but they inflate the variance of individual coefficients, so signs flip and magnitudes swing between runs — the matrix is nearly rank-deficient, which is Module 02's rank condition showing up as noisy inference.

9 · Bias–variance for models

Ask one question: if I retrained this model on a different sample of the same size, how different would the prediction be — and how far is it from the truth on average? The answer decomposes into exactly three terms:

• $\mathrm{Bias}^2[\hat f(x)] = \big(\mathbb{E}[\hat f(x)] - f(x)\big)^2$ — the squared distance between the average fit over all possible training sets and the truth. High bias means the model's idea of the world is wrong in the same direction every time: too simple, wrong shape, missing features. It does not go away with more data.
• $\mathrm{Var}[\hat f(x)] = \mathbb{E}\big[(\hat f(x) - \mathbb{E}[\hat f(x)])^2\big]$ — how far individual fits scatter around that average. High variance means the model chases the particular sample it happened to see: swap the training set and the prediction moves. It does shrink with more data.
• $\sigma^2$ — the noise in the outcome itself. No model, of any complexity, removes it. This is the floor under your test error.

The picture to keep is a bullseye. Each training set produces one shot. High bias, low variance: the shots land in a tight knot — consistent, and consistently in the wrong place (systematic offset). High variance, low bias: the shots scatter widely around the centre — on average right, but any single shot could be anywhere. Low bias, low variance: a tight knot in the bullseye — that is the goal, and you pay for every step towards it.

Why the levers work the way they do: a more flexible model (higher-degree polynomial, deeper tree, unregularised net) can represent more shapes, so the average fit can get closer to the truth — bias falls — but it also has more freedom to fit each sample's accidents, so variance rises. More data constrains those accidents: every extra point is another vote against the sample-specific wiggle, so variance falls while bias stays exactly where it was. Irreducible noise moves for nobody. Block 3 runs the whole story on noisy sine data: at degree 1 train and test error are both high (bias), by degree 5 test error has hit the floor, and from degree 7 train error keeps falling (0.053) while test error explodes (0.036 → 5.8 → 38 → 8036) — variance, all of it.

Learning curves are the empirical view: plot train and validation error against training-set size. High variance shows as a large, slowly closing gap between the curves (more data would help — go get it). High bias shows as both curves already close together and already high (more data will not help — change the model). Reading that gap correctly tells you which of the three levers to pull, and it is the single most useful diagnostic in applied ML.

10 · Generalisation, validation and regularization

The test set exists to estimate the error of the final model on unseen data. The validation set exists because something has to absorb the damage of all your decisions: hyperparameters, feature choices, stopping time, model family. The test set stops being a test set the first time you look at it to make a decision — it has become validation, and its number is now optimistically biased by exactly as much as you looked.

k-fold cross-validation spends small data better: split into $k$ folds, train on $k-1$, score the held-out fold, repeat so every point is scored once, average. What it estimates is the expected performance of a model trained on $(k-1)/k$ of your data — a slightly pessimistic proxy for the final model trained on all of it — plus a spread across folds that tells you how much the score moves with the sample. It is a model-selection instrument; it is not a substitute for a final held-out number, because every choice you made was made against CV scores.

Leakage is the bug that makes everything look fine: scaling or imputing before splitting (the scaler saw the test data), duplicate rows or the same user appearing in both splits (the model recognises rather than generalises), and temporal leakage in time series — training on last month to predict last month. Time-ordered data must be split by time, always: train on the past, validate on the future, because that is the direction the world arrives in.

Regularization is variance reduction with a bias budget. Ridge and lasso add a penalty that shrinks coefficients towards zero — you accept a little bias in exchange for a large drop in coefficient variance, which is exactly the trade of section 9 moved into parameter space. In the MAP language of Module 04, ridge is a Gaussian prior on the weights and lasso a Laplace one: same objective, Bayesian justification. Early stopping does the same job implicitly — stopping before convergence keeps the model in the low-variance region of its training path.

11 · Metrics and their error bars

Choosing a metric is choosing what you will optimise, and each of the standard ones is blind to something:

ROC-AUC has a probabilistic reading worth internalising: AUC is the probability that a randomly chosen positive example is scored above a randomly chosen negative one. 0.5 is random ordering, 1.0 is perfect. It measures ranking, not thresholding and not probability quality. Log-loss measures the quality of the probabilities themselves and has the property Module 06 builds from scratch: it is a proper scoring rule, meaning the way to minimise your expected score is to report your true probability — which is exactly what you want a training loss to do.

Calibration is the claim that "0.8 means 80%". Check it with a reliability curve (or reliability diagram): bucket predictions by confidence and plot observed frequency against predicted. Perfect calibration is the diagonal; below it you are overconfident, above it underconfident. A model can be superbly ranked (high AUC) and badly calibrated at the same time — which is why ranking metrics and probability metrics answer different questions.

That last point is the module in miniature. If 0.912 → 0.918 looks like a 0.6-point win, the honest question is whether that difference survives resampling the test set. The clean way to answer it is a paired bootstrap: because both models scored the same examples, you resample examples (not predictions), recompute the difference each time, and read the 2.5th and 97.5th percentiles. The pairing matters — it removes the noise shared by both models and gives a much tighter interval than treating the two scores as independent. The exercise below does exactly this on 40,000 predictions.

Block 1 · Bootstrap confidence interval vs the normal approximation

The bootstrap is statistics without formulas: resample the data you have, recompute the statistic, repeat 10,000 times, and read off the percentiles. It works for any statistic — including ones (like a median or an AUC) with no tidy textbook standard error.

import numpy as np
rng = np.random.default_rng(0)

x = rng.lognormal(0.0, 1.0, size=500)     # skewed data: the median's own sampling
n = x.size                                # distribution is skewed too
B = 10000                                 # bootstrap resamples
idx = rng.integers(0, n, size=(B, n))     # B resamples of size n, with replacement

med = np.median(x)
boot = np.median(x[idx], axis=1)          # 10 000 medians of resampled data
lo, hi = np.percentile(boot, [2.5, 97.5])
se = boot.std(ddof=1)                     # bootstrap standard error
print(f"median {med:.3f}  bootstrap [{lo:.3f}, {hi:.3f}]  normal [{med-1.96*se:.3f}, {med+1.96*se:.3f}]")
# -> median 0.951  bootstrap [0.842, 1.079]  normal [0.828, 1.074]
# the percentile interval sits to the right and is asymmetric about the estimate:
# gaps of -0.109 / +0.128, because the sampling distribution inherits the skew.

mu = x.mean()                             # same data, a symmetric statistic
bm = x[idx].mean(axis=1)
lo2, hi2 = np.percentile(bm, [2.5, 97.5])
se2 = x.std(ddof=1) / np.sqrt(n)          # textbook SE for the mean
print(f"mean {mu:.3f}  bootstrap [{lo2:.3f}, {hi2:.3f}]  normal [{mu-1.96*se2:.3f}, {mu+1.96*se2:.3f}]")
# -> mean 1.598  bootstrap [1.441, 1.779]  normal [1.425, 1.771]   (agree to ~1%)

Where they disagree is the lesson: for a symmetric statistic the two agree to within rounding, and for a skewed one the normal approximation is symmetric when the truth is not. Percentile bootstrap costs nothing extra and respects the shape.

Block 2 · A/B test by hand, then calibrate it 500 times

import numpy as np
from math import erfc, sqrt

def ztest(x1, n1, x2, n2):
    """Two-proportion pooled z-test -> (diff, z, two-sided p, SE)."""
    p1, p2 = x1/n1, x2/n2
    pp = (x1 + x2)/(n1 + n2)              # pooled rate: what H0 assumes
    se = sqrt(pp*(1-pp)*(1/n1 + 1/n2))
    z = (p2 - p1)/se                      # difference over chance-level noise
    return p2 - p1, z, erfc(abs(z)/sqrt(2)), se

n = 5000
r = np.random.default_rng(42)
a, b = r.binomial(n, 0.08), r.binomial(n, 0.09)   # a REAL 1pp lift
d, z, p, se = ztest(a, n, b, n)
print(f"arms {a/n:.4f} vs {b/n:.4f}  diff {d:+.4f}  z {z:.2f}  p {p:.4f}  CI [{d-1.96*se:+.4f}, {d+1.96*se:+.4f}]")
# -> arms 0.0780 vs 0.0928  diff +0.0148  z 2.65  p 0.0081  CI [+0.0038, +0.0258]

r2 = np.random.default_rng(7)
hit = sum(ztest(r2.binomial(n,0.08), n, r2.binomial(n,0.09), n)[2] < 0.05 for _ in range(2000))
print("power 8% vs 9% at n=5000:", round(hit/2000, 3))       # -> 0.448 (a real lift, often missed)

r3 = np.random.default_rng(1005)                             # now NO lift at all
fp = sum(ztest(r3.binomial(n,0.08), n, r3.binomial(n,0.08), n)[2] < 0.05 for _ in range(500))
print(f"false positives: {fp}/500 = {fp/500:.3f}")            # -> 25/500 = 0.050  (calibrated!)

def peek(seed, days=40, per_day=500, R=2000):
    """Peek daily and stop at the first p < 0.05, when the null is true."""
    g = np.random.default_rng(seed); fired = 0
    for _ in range(R):
        ca = cb = 0
        for t in range(1, days+1):
            ca += g.binomial(per_day, 0.08); cb += g.binomial(per_day, 0.08)
            if ztest(ca, per_day*t, cb, per_day*t)[2] < 0.05:
                fired += 1; break
    return fired / R
print("40 daily looks:", round(peek(1), 3), " one look at the end:", round(peek(1, days=1, per_day=20000), 3))
# -> 40 daily looks: 0.301   one look at the end: 0.05

Three numbers to carry out of this block. The test is calibrated: with no real difference it still fired 25 times in 500 — that is what $\alpha = 0.05$ means, and Block 2 lets you watch it happen. It is also underpowered at this sample size: a real 1pp lift is caught only 45% of the time, so "the test showed nothing" is mostly a sample-size statement. And peeking destroys the calibration: 40 looks and you are celebrating a false win 30% of the time.

Block 3 · Bias–variance, demonstrated

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(3)
f = lambda t: np.sin(2*np.pi*t)
NOISE, N_TRAIN, R = 0.3, 30, 50           # noise sd, train size, repetitions
grid = np.linspace(0, 1, 100)             # fixed test grid, noise-free
y_true = f(grid)
DEG = np.arange(1, 13)

train = np.zeros((R, DEG.size)); test = np.zeros((R, DEG.size))
preds = np.zeros((R, DEG.size, grid.size))
for r in range(R):                        # a NEW training set every repetition
    xt = rng.uniform(0, 1, N_TRAIN)
    yt = f(xt) + rng.normal(0, NOISE, N_TRAIN)
    for k, d in enumerate(DEG):
        beta = np.polyfit(xt, yt, d)
        train[r, k] = np.mean((np.polyval(beta, xt) - yt)**2)
        preds[r, k] = np.polyval(beta, grid)
        test[r, k] = np.mean((preds[r, k] - y_true)**2)

mean_pred = preds.mean(0)                 # the AVERAGE model over repetitions
bias2 = ((mean_pred - y_true)**2).mean(1) # how far the average fit sits from truth
var = preds.var(0, ddof=1).mean(1)        # how far the fits scatter around it

print("deg   train     test     bias^2      var   bias^2+var")
for k in range(DEG.size):
    if DEG[k] in (1, 3, 5, 7, 9, 11):
        print(f"{DEG[k]:3d} {train[:,k].mean():8.4f} {test[:,k].mean():9.4f}"
              f" {bias2[k]:9.4f} {var[k]:9.4f} {bias2[k]+var[k]:10.4f}")
#  deg   train     test     bias^2      var   bias^2+var
#    1   0.2645   0.2242    0.2035   0.0211     0.2246
#    3   0.0782   0.0239    0.0058   0.0185     0.0243
#    5   0.0704   0.0360    0.0007   0.0360     0.0367
#    7   0.0644   5.8316    0.1563   5.7911     5.9475
#    9   0.0593  38.0266    0.7456  38.0418    38.7874
#   11   0.0534 8035.99  166.4049 8030.19   8196.59
plt.plot(DEG, train.mean(0), label="train MSE")
plt.plot(DEG, test.mean(0), label="test MSE")
plt.plot(DEG, bias2, label="bias^2"); plt.plot(DEG, var, label="variance")
plt.xscale("linear"); plt.yscale("log"); plt.legend(); plt.show()

Read the columns: bias² + variance matches the test error to within a couple of percent — the decomposition is an identity, not a metaphor. Degree 1 is dominated by bias (0.20 of 0.22), degree 5 is the sweet spot (bias² 0.0007, variance 0.0360), and from degree 7 the variance term swallows everything while train error happily keeps dropping. Train error is monotone in complexity; test error is not. That divergence is what overfitting is.

Block 4 · Residual analysis and honest validation

import numpy as np
from sklearn.datasets import make_regression
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold, cross_val_score

X, y = make_regression(n_samples=250, n_features=100, n_informative=8,
                       noise=50, random_state=0)
n, p = X.shape
m = LinearRegression().fit(X, y)
r2 = m.score(X, y)                        # resubstitution score: always flattering
adj = 1 - (1 - r2) * (n - 1) / (n - p - 1)
cv = cross_val_score(LinearRegression(), X, y,
                     cv=KFold(5, shuffle=True, random_state=0), scoring="r2")
print(f"train R2 {r2:.4f} | adjusted R2 {adj:.4f} | 5-fold CV R2 {cv.mean():.4f} +/- {cv.std():.4f}")
# -> train R2 0.9673 | adjusted R2 0.9454 | 5-fold CV R2 0.8662 +/- 0.0320

resid = y - m.predict(X)
print(f"residual mean {resid.mean():+.2e}, sd {resid.std(ddof=p+1):.2f} (true noise sd 50)")
# -> residual mean -5.91e-15, sd 48.98 (true noise sd 50)

# a WRONG model leaves a PATTERN in the residuals, not just a worse score
rng = np.random.default_rng(0)
x = np.linspace(-3, 3, 200); yy = x**2 + rng.normal(0, 1, 200)
r1 = yy - np.polyval(np.polyfit(x, yy, 1), x)      # fit a line to a parabola
r2c = yy - np.polyval(np.polyfit(x, yy, 2), x)     # fit the right shape
print(f"line  : corr(resid, x^2) = {np.corrcoef(r1, x**2)[0,1]:+.3f}, resid sd {r1.std(ddof=2):.3f}")
print(f"curve : corr(resid, x^2) = {np.corrcoef(r2c, x**2)[0,1]:+.3f}, resid sd {r2c.std(ddof=3):.3f}")
# -> line  : corr(resid, x^2) = +0.940, resid sd 2.831
# -> curve : corr(resid, x^2) = +0.000, resid sd 0.965

Two lessons in the numbers. The train score (0.967) flatters the model by a tenth of a point against honest 5-fold CV (0.866), and adjusted $R^2$ (0.945) only partly corrects — it charges a flat rent per feature and cannot know that 92 of those features were useless. And the residual pattern detects a wrong model that the score alone would not: correlation with the missing term is +0.94 for the straight line, 0.00 once you fit the right shape.

import numpy as np
rng = np.random.default_rng(0)
n = 40000
counts = np.array([0.900, 0.012, 0.018, 0.070]) * n   # both right / only A / only B / both wrong
cell = np.repeat(np.arange(4), np.rint(counts).astype(int))
rng.shuffle(cell)
a = (cell == 0) | (cell == 1)                         # model A correct
b = (cell == 0) | (cell == 2)                         # model B correct
print(f"acc_A {a.mean():.4f}  acc_B {b.mean():.4f}  diff {b.mean()-a.mean():+.4f}")
# -> acc_A 0.9120  acc_B 0.9180  diff +0.0060

dvec = b.astype(int) - a.astype(int)                  # -1, 0, +1 per example
diffs = np.empty(5000)
for i in range(5000):                                 # resample EXAMPLES, both models together
    diffs[i] = dvec[rng.integers(0, n, n)].mean()
print("paired 95% CI:", np.round(np.percentile(diffs, [2.5, 97.5]), 4))   # -> [0.0043 0.0077]

diffs2 = np.empty(5000)                               # unpaired: forget they shared examples
for i in range(5000):
    diffs2[i] = b[rng.integers(0, n, n)].mean() - a[rng.integers(0, n, n)].mean()
print("unpaired 95% CI:", np.round(np.percentile(diffs2, [2.5, 97.5]), 4))  # -> [0.0023 0.0098]

Multiple choice — instant feedback

Open questions — answer aloud, no notes

Re-answer all eight in 1 day, then 3 days, then 7 days — from memory, in writing. If question 4 or 6 comes back fuzzy, re-read only that section.