Statistics, probability run backwards.

The previous chapter built up the machinery of distributions and expectations on the assumption that you know the underlying probability law. In real applications you almost never do. You see a finite sample of measurements — gene expression levels, click-through rates, sensor readings, image pixels — and you have to infer the distribution that produced them, or some functional of it: a mean, a variance, a regression coefficient, a class boundary, a treatment effect. That inversion problem is statistics.

The setup looks like this. We posit a statistical model: a family of probability distributions $\{p_\theta : \theta \in \Theta\}$ indexed by an unknown parameter $\theta$ in some parameter space $\Theta$. The data $X_1, \ldots, X_n$ are assumed to be drawn i.i.d. from one distribution in that family — call it $p_{\theta^*}$ — but we do not know which $\theta^*$ is the true one. The job is to use the data to make defensible statements about $\theta^*$: a single best guess (point estimation), a range of plausible values (confidence intervals), or a yes/no decision about a candidate value (hypothesis testing).

The chapter that follows is structured around the four tasks frequentists care about most. Sections 02–04 build the theory of estimation: what makes one estimator better than another. Section 05 turns single-number estimates into intervals. Sections 06–08 cover testing — formal procedures for deciding between hypotheses, and the modern problem of doing so honestly when you run many tests at once. Sections 09–12 develop regression, the most important applied branch of the field. Section 13 covers resampling — the bootstrap and friends — which made many of the older asymptotic formulas obsolete. Section 14 is a brief frequentist–Bayesian reconciliation. Sections 15–17 cover experimental design and other topics that bridge into causal inference and time series. Section 18 closes by mapping every concept onto its modern incarnation in machine learning.

The prerequisites are exactly the chapters that came before: linear algebra (for regression and multivariate distributions), calculus (for likelihoods and score functions), optimization (for estimating $\theta$ by minimising a loss), and probability (for everything). If those chapters felt like a long detour, this one is where the road rejoins the highway.

A random sample of size $n$ from a distribution $F$ is a sequence $X_1, \ldots, X_n \stackrel{\text{i.i.d.}}{\sim} F$. A statistic is any function of the sample,

that does not depend on the unknown parameter $\theta$. The sample mean $\bar{X}_n = \frac{1}{n}\sum_{i=1}^n X_i$, the sample variance $S_n^2 = \frac{1}{n-1}\sum_{i=1}^n (X_i - \bar{X}_n)^2$, the sample median, the sample maximum, the sample correlation — all are statistics. Because each is a function of random inputs, each has its own distribution: the sampling distribution of the statistic.

The most studied sampling distribution is that of the sample mean. If the $X_i$ have mean $\mu$ and variance $\sigma^2$, then by linearity of expectation and the variance of a sum of independents,

The standard deviation of the sample mean — sometimes called the standard error of the mean — is therefore $\sigma / \sqrt{n}$. By the central limit theorem, for large $n$,

This single fact is the engine behind most classical confidence intervals and tests. It also explains the famous "$1/\sqrt{n}$" rate of convergence: doubling your dataset reduces uncertainty about the mean by only $\sqrt{2}$.

Three other classical sampling distributions show up constantly. If $X_1, \ldots, X_n \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(\mu, \sigma^2)$, then

The chi-squared distribution is the law of a sum of squared standard normals; it powers variance estimation and goodness-of-fit tests. The t-distribution is what you get when you replace the unknown $\sigma$ by the data's own $S_n$ — it has heavier tails than a normal, especially for small $n$. The F-distribution is the ratio of two chi-squareds; it powers ANOVA and tests of nested regression models. All three converge to standard reference distributions as $n \to \infty$, but their finite-sample form is what makes small-sample inference rigorous.

A point estimator for a parameter $\theta$ is a statistic $\hat{\theta}_n = \hat{\theta}(X_1, \ldots, X_n)$ that we use as our best single guess. Since $\hat{\theta}_n$ is a function of random data, it is itself a random variable; its quality is measured by properties of its distribution.

Bias. The bias of an estimator is its expected error,

An estimator is unbiased if its bias is zero for every $\theta$. The sample mean is unbiased for $\mu$. The sample variance with $1/(n-1)$ is unbiased for $\sigma^2$; the version with $1/n$ is biased downward. The MLE is famously biased in finite samples but its bias usually shrinks at rate $1/n$.

Variance and mean squared error. Two unbiased estimators can still differ wildly in spread. The mean squared error of $\hat{\theta}$ is

This is the bias–variance decomposition that will reappear under different names in regression (Section 12) and in machine learning (Section 18). It says you can pay for variance reduction with bias and vice versa, and the right trade is the one with lowest total error — not the one with the most reassuring zero in the bias column.

Consistency. An estimator is consistent if it converges in probability to the true parameter as $n \to \infty$:

Consistency is a minimum requirement: an estimator that doesn't get closer to the truth as you collect more data is unfit for purpose. Sufficient conditions include being unbiased with variance going to zero (immediate from Chebyshev) but the property is more general.

The Cramér–Rao lower bound. Among all unbiased estimators, how small can the variance possibly be? The answer involves the Fisher information

a quantity that measures how sharply the likelihood is curved around the true parameter. The Cramér–Rao inequality states that for any unbiased estimator $\hat{\theta}$ in a regular model,

An estimator that achieves this bound is called efficient. The MLE is asymptotically efficient under mild regularity conditions, which is one of the deepest reasons it dominates statistical practice. (And it generalises: the multivariate Cramér–Rao bound replaces $I(\theta)$ with the Fisher information matrix and the inequality with one between covariance matrices.)

Method of moments (MoM). If a model has $k$ parameters $\theta = (\theta_1, \ldots, \theta_k)$, set the first $k$ population moments equal to the corresponding sample moments and solve:

For a normal model $\mathcal{N}(\mu, \sigma^2)$ with two parameters, this gives the sample mean for $\mu$ and the sample variance (with $1/n$) for $\sigma^2$ — clean and intuitive. MoM estimators are usually consistent but not always efficient; their main role today is producing decent starting points for iterative procedures and supplying closed-form alternatives in models where the likelihood is intractable.

Maximum likelihood (MLE). Given data $X_1, \ldots, X_n \stackrel{\text{i.i.d.}}{\sim} p_\theta$, the likelihood is the joint density viewed as a function of $\theta$:

The log-likelihood $\ell(\theta) = \log L(\theta) = \sum_i \log p_\theta(X_i)$ is easier to differentiate and almost always what people actually optimise. The maximum likelihood estimator is

Setting the derivative — the score function $U(\theta) = \partial \ell / \partial \theta$ — to zero and solving yields the MLE in models smooth enough to admit calculus. For a normal model, $\hat{\mu}_{\text{MLE}} = \bar{X}_n$ and $\hat{\sigma}^2_{\text{MLE}} = \frac{1}{n}\sum_i (X_i - \bar{X}_n)^2$ — note the $1/n$, not $1/(n-1)$: the MLE for the variance is biased downward in this model.

Asymptotic normality of the MLE. Under mild regularity conditions (identifiability, smooth densities, an interior maximum), the MLE satisfies

In words: MLEs are asymptotically unbiased, asymptotically Gaussian, and asymptotically efficient — their variance attains the Cramér–Rao bound. This is the theoretical justification for using $\hat{I}(\hat{\theta})^{-1}/n$ as an approximate covariance matrix when reporting standard errors of fitted parameters in everything from logistic regression to deep neural networks (where it shows up as the observed Fisher information or its diagonal approximation in second-order optimisers).

Numerical reality. Closed-form MLEs are rare outside textbook examples. In practice the log-likelihood is maximised by the optimisation algorithms of Chapter 03: Newton's method (using the Hessian, which is minus the observed Fisher information), Fisher scoring (replacing the Hessian with $-\,I(\theta)$), or first-order gradient methods. For very large datasets the gradient is approximated by stochastic mini-batches, and "MLE" silently becomes "stochastic gradient ascent on the log-likelihood" — i.e., training a neural network.

A $(1-\alpha)$ confidence interval for $\theta$ is a pair of statistics $L(X), U(X)$ such that

The interval is random; the parameter is fixed. The probability statement is about the procedure: if you repeated the entire experiment many times, in $1-\alpha$ of those replications the produced interval would cover the true $\theta$. It is not a statement about the probability that this particular interval contains $\theta$ — that probability is either 1 or 0, since both endpoints and $\theta$ are fixed once data are observed. (This is the line that separates the frequentist and Bayesian readings; the Bayesian credible interval, Section 14, has the natural "probability that $\theta$ lies in the interval" interpretation, but pays a different price.)

The Normal interval for a mean. If $X_1, \ldots, X_n \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(\mu, \sigma^2)$ with known $\sigma$, then $\bar{X}_n \sim \mathcal{N}(\mu, \sigma^2/n)$ and

is an exact $(1-\alpha)$ interval, where $z_{1-\alpha/2}$ is the relevant standard-normal quantile (1.96 for $\alpha = 0.05$). When $\sigma$ is unknown, replace it by the sample standard deviation $S_n$ and the normal quantile by a t-quantile $t_{n-1, 1-\alpha/2}$:

For large $n$ the t-quantile is essentially the normal quantile and you can stop worrying about which one to use; for small samples (say $n < 30$) it matters.

Asymptotic intervals via the MLE. When closed-form pivots are unavailable, the asymptotic normality of the MLE gives

a "Wald interval" valid for large $n$. It is the workhorse for parameter uncertainty in logistic regression, GLMs, and generic likelihood models. Two close cousins — the score interval (inverting a test based on the score function) and the likelihood-ratio interval (the set $\{\theta : 2[\ell(\hat\theta) - \ell(\theta)] \leq \chi^2_{1, 1-\alpha}\}$) — have better small-sample behaviour but require more computation.

Bootstrap intervals. When even asymptotics are inconvenient — heavy-tailed data, small samples, complicated functionals — the bootstrap (Section 13) offers a fully numerical alternative: simulate the sampling distribution by resampling from the data itself. Modern intervals for almost any quantity in modern statistics — variances of trained model parameters, complicated functionals of estimated distributions, derived quantities of regression coefficients — are bootstrap intervals.

The setup is two competing models of the world: the null hypothesis $H_0: \theta \in \Theta_0$ and the alternative hypothesis $H_1: \theta \in \Theta_1$. A test is a function $\phi(X) \in \{0, 1\}$ that decides whether to reject $H_0$ (output 1) or fail to reject (output 0). Two kinds of error are possible:

• Type I error — rejecting $H_0$ when it is true. Its probability under $H_0$ is the test's significance level $\alpha = \mathbb{P}_{H_0}(\phi = 1)$.
• Type II error — failing to reject $H_0$ when $H_1$ is true. Its probability is $\beta(\theta) = \mathbb{P}_\theta(\phi = 0)$ for $\theta \in \Theta_1$. The complement $1 - \beta(\theta)$ is the test's power.

The classical paradigm fixes the level $\alpha$ — typically 0.05 or 0.01 — and then designs a test with the largest possible power against alternatives of interest. Note the asymmetry: a test cannot generally control both errors simultaneously, and the convention is to protect against false rejections.

The Neyman–Pearson lemma. For two simple hypotheses $H_0: \theta = \theta_0$ versus $H_1: \theta = \theta_1$, the most powerful test of level $\alpha$ rejects when the likelihood ratio is large:

where $k$ is chosen so the false-positive rate equals $\alpha$. This is the only "optimality theorem" of classical testing, and it underwrites essentially every test in the next section. For composite hypotheses (e.g. $H_1: \theta \neq \theta_0$) the likelihood-ratio idea generalises to the generalised likelihood ratio test (GLRT)

Under regularity, $-2 \log \Lambda \xrightarrow{d} \chi^2_k$ where $k$ is the difference in the dimension of the two parameter spaces — Wilks' theorem — which is what powers the F-test for nested regression models in Section 10.

One-sided vs. two-sided. When the alternative has a direction ($\theta > \theta_0$), use a one-sided test. When it does not ($\theta \neq \theta_0$), use a two-sided test. The two-sided test halves the rejection region on each side, doubling the p-value of the corresponding one-sided test. Reporting one-sided p-values when the direction was chosen after seeing the data is one of the classic and unfortunately common ways to inflate false positives.

z-test for a mean. $H_0: \mu = \mu_0$ versus $H_1: \mu \neq \mu_0$, with $\sigma^2$ known. The statistic

is compared to standard normal quantiles. In practice the variance is rarely known; this test is a stepping stone, not an end product.

One-sample t-test. Same as above but with $\sigma$ unknown, replaced by $S_n$:

The "Student's t-test" — published anonymously by Gosset in 1908 because his employer, Guinness, considered the result a trade secret. It is robust to mild departures from normality and is the standard tool for comparing a sample mean to a hypothesised value.

Two-sample t-test. $H_0: \mu_1 = \mu_2$ for two independent samples. The pooled-variance version assumes $\sigma_1 = \sigma_2$; Welch's version does not and is the safer default. Either way the statistic is approximately t-distributed under $H_0$.

Paired t-test. When the two samples are paired (same subjects measured before and after), reduce to a one-sample t-test on the differences $D_i = X_i - Y_i$. This eliminates between-subject variability and usually yields much higher power than the two-sample version on the same data.

$\chi^2$ tests. Two flavours, both based on the same statistic but applied differently:

where $O_i$ is the observed count in cell $i$ and $E_i$ is the count expected under $H_0$. The goodness-of-fit test compares observed counts to a hypothesised distribution; df = (cells $-$ 1 $-$ estimated parameters). The test of independence in a contingency table compares observed cell counts to those expected under independence of the row and column variables; df = (rows $-$ 1)(cols $-$ 1).

F-test for equality of variances. Ratio of two sample variances; under $H_0: \sigma_1 = \sigma_2$ in normal samples,

Highly sensitive to non-normality; rarely used as a standalone test today and largely replaced by Levene's or Brown–Forsythe tests. The F-distribution itself remains central — it is the reference distribution for ANOVA and for nested-model comparisons in regression.

ANOVA. Analysis of variance compares means across $k \geq 3$ groups by partitioning the total sum of squares into "between-groups" and "within-groups" components. The ratio

is large when the group means differ substantially relative to within-group variation. ANOVA is exactly the F-test for nested regression models with categorical predictors, a connection made explicit in Section 10.

Non-parametric alternatives. When normality fails or the data are ordinal, the Wilcoxon signed-rank test (paired) and Mann–Whitney U test (two-sample) replace t-tests; the Kruskal–Wallis test replaces ANOVA; the permutation test replaces almost everything by computing the test statistic's exact null distribution from random shuffles of group labels.

Run $m$ independent tests at level $\alpha$ each, all with true nulls. The probability that at least one rejects is

With $\alpha = 0.05$ and $m = 20$, that's already 64%; with $m = 1{,}000$ (a typical genome-wide scan, or a typical hyperparameter sweep) it is essentially 1. Any culture that rewards "find a significant result" while ignoring the count of tests performed will produce a steady stream of false positives. This was the diagnosis of the replication crisis in psychology, biology, and economics in the 2010s — many "discoveries" failed to replicate not because the science was sloppy but because the multiple-testing arithmetic was ignored.

Two error rates have to be distinguished:

• Family-wise error rate (FWER). The probability of any false rejection across the family of tests. Conservative; appropriate when even one false positive is a serious problem.
• False discovery rate (FDR). The expected fraction of false positives among the tests you reject. Less conservative; appropriate when you are screening a large population for promising candidates.

Bonferroni correction controls FWER by testing each hypothesis at level $\alpha / m$:

Simple, conservative, and the right tool for confirmatory studies with a handful of hypotheses. Holm's step-down procedure is uniformly more powerful than Bonferroni at the same FWER and is essentially free; there is no reason to prefer Bonferroni in modern software.

Benjamini–Hochberg (BH). Controls FDR rather than FWER. Sort the p-values $p_{(1)} \leq \cdots \leq p_{(m)}$ and reject all hypotheses with rank $\leq k^*$ where

for a target FDR level $q$. BH is the workhorse of high-dimensional screening — used routinely in genomics, neuroimaging, and any setting where thousands of features are tested and a small fraction of false positives is acceptable.

For ML practitioners the analogues are everywhere: hyperparameter sweeps, dataset pre-processing choices, picking the random seed that worked, comparing against many baselines until one loses. The discipline that has emerged in recent years — train/validation/test splits with the test set touched once, results reported with confidence intervals across seeds, ablations rather than cherry-picked best runs — is essentially the multiple-testing discipline of statistics, transferred to a new field.

The model: an outcome $Y$ is a linear function of predictors $\mathbf{x} = (x_1, \ldots, x_p)^\top$ plus noise,

Stack $n$ observations into a vector $\mathbf{y} \in \mathbb{R}^n$ and an $n \times (p+1)$ design matrix $\mathbf{X}$ whose first column is all ones (for the intercept) and whose remaining columns hold the predictors. Then $\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ with $\boldsymbol{\beta} \in \mathbb{R}^{p+1}$ unknown.

The OLS estimator. The ordinary-least-squares estimator minimises the sum of squared residuals:

The closed-form solution comes from setting the gradient of the squared loss to zero. The matrix $(\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top$ is the Moore–Penrose pseudoinverse when $\mathbf{X}$ has full column rank; numerically you should use a QR or SVD factorisation rather than forming and inverting $\mathbf{X}^\top \mathbf{X}$ explicitly (the condition number squares).

The geometric picture. View $\mathbf{y}$ as a vector in $\mathbb{R}^n$ and the columns of $\mathbf{X}$ as spanning a $(p+1)$-dimensional subspace. Then $\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol{\beta}}$ is the orthogonal projection of $\mathbf{y}$ onto that subspace, and the residual vector $\hat{\boldsymbol{\varepsilon}} = \mathbf{y} - \hat{\mathbf{y}}$ is orthogonal to every column of $\mathbf{X}$. The "normal equations" $\mathbf{X}^\top(\mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}}) = \mathbf{0}$ are exactly that orthogonality condition. The "$R^2$" reported with every regression is the cosine-squared of the angle between the projected and total signals.

The likelihood picture. Under the assumed Gaussian noise, the log-likelihood is

Maximising over $\boldsymbol{\beta}$ for any fixed $\sigma^2$ minimises $\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|^2$ — i.e., gives OLS. So least squares is maximum likelihood under Gaussian noise, and that is why mean-squared error is the natural loss for so many models. The corresponding MLE for the variance is $\hat{\sigma}^2 = \|\mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}}\|^2 / n$; a slightly larger denominator $(n - p - 1)$ gives the unbiased version.

The next two sections build on this scaffold: Section 10 turns the fit into inference about which coefficients matter, and Section 11 generalises beyond Gaussian noise to logistic regression, Poisson regression, and the rest of the GLM family.

Under the standard linear model with Gaussian noise,

The diagonal entries of $\sigma^2 (\mathbf{X}^\top \mathbf{X})^{-1}$ are the variances of the individual coefficient estimates. Replacing $\sigma^2$ by its unbiased estimator $\hat{\sigma}^2 = \|\mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}}\|^2 / (n - p - 1)$ gives standard errors $\operatorname{SE}(\hat{\beta}_j)$, and

These are the t-statistics and p-values printed next to every regression coefficient in every statistical software output. The null against which they are tested is almost always $\beta_j = 0$ — "this predictor has no marginal effect after adjusting for the others" — and the corresponding p-value answers exactly that question, no more.

F-test for nested models. To compare a full model with $p$ predictors against a reduced model with $q < p$ of them, fit both and compute

A large $F$ means the additional predictors meaningfully reduce residual variance. The same machinery is what powers ANOVA: each group is a categorical predictor, and the F-test compares the model with all groups to the intercept-only model.

Confidence vs. prediction intervals. Two distinct intervals are routinely confused. For a new test point $\mathbf{x}_0$:

The confidence interval bounds the regression line; the prediction interval bounds an individual future observation and is always wider, because a new observation has its own irreducible noise $\sigma^2$ on top of the uncertainty in $\hat{\boldsymbol{\beta}}$.

Diagnostics. The classical assumptions — linearity, homoscedasticity, normality of residuals, independence — are checked by inspecting residual plots, not just by tests. A scatter of residuals versus fitted values that fans out signals heteroscedasticity. A Q–Q plot of residuals against normal quantiles checks normality. High-leverage points (large diagonal entries of the hat matrix $\mathbf{H} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top$) are observations whose predictors lie far from the rest; they exert outsized influence on the fit and deserve scrutiny. Cook's distance and the various standardised-residual flavours are formalisations of the same idea.

A generalised linear model has three ingredients:

1. A random component: a distribution for $Y$ from the exponential family — Bernoulli, Binomial, Poisson, Gaussian, Gamma, etc.
2. A linear predictor: $\eta = \mathbf{x}^\top \boldsymbol{\beta}$.
3. A link function $g$ relating the mean of $Y$ to the linear predictor: $g(\mathbb{E}[Y]) = \eta$.

Choosing the distribution and link recovers a familiar model. Gaussian + identity link gives ordinary linear regression. Bernoulli + logit link gives logistic regression:

where $\sigma(z) = 1/(1 + e^{-z})$ is the logistic function. Coefficients in logistic regression are log odds ratios: a one-unit increase in $x_j$ multiplies the odds of $Y = 1$ by $e^{\beta_j}$.

Poisson + log link gives Poisson regression for count data:

A coefficient is then a log rate ratio: the multiplicative effect on the expected count. The standard correction for varying exposure (length of follow-up, area sampled, time at risk) is to add an offset $\log E_i$ to the linear predictor.

Fitting GLMs. The MLE for a GLM has no closed form (logistic regression notoriously does not), so it is computed iteratively. The classical algorithm is iteratively reweighted least squares (IRLS): at each iteration, linearise the link function and solve a weighted OLS problem; the weights depend on the current fit. IRLS is just Newton–Raphson on the log-likelihood with a particularly clean update; the asymptotic covariance of the resulting MLE is the inverse Fisher information, exactly as in Section 04, and software returns it as the standard error matrix for the coefficients.

Inference. Wald tests, likelihood-ratio tests (using deviance differences against $\chi^2_{p-q}$ from Wilks' theorem), and score tests all transfer from linear regression to GLMs with minor changes. The deviance $D = 2[\ell_{\text{sat}} - \ell_{\text{model}}]$ — twice the log-likelihood gap to a "saturated" model that interpolates the data — replaces the residual sum of squares as the basic measure of fit, and behaves like a $\chi^2$ for nested model comparisons.

Modern ML inherits this scaffolding wholesale. A neural network with a sigmoid output trained with binary cross-entropy is a logistic regression with a learned, non-linear feature map. A network with a softmax output trained with categorical cross-entropy is a multinomial logistic regression. The optimisers are different — SGD instead of IRLS — but the loss functions, the link functions, and the asymptotic standard errors all come from this one chapter.

Decompose the prediction error of an estimator $\hat{f}$ at a fixed input $\mathbf{x}_0$:

The irreducible noise $\sigma^2$ is what it is. Bias and variance, however, are knobs the modeller controls. Simple models (few predictors, large penalties) have low variance and high bias; complex models (many predictors, no regularisation) have low bias and high variance. The total error is U-shaped in complexity; the goal is to land in the trough.

Ridge regression. Add an $\ell_2$ penalty on the coefficients to the OLS loss:

Two consequences are worth seeing immediately. First, the matrix $\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I}$ is invertible for any $\lambda > 0$ even if $\mathbf{X}^\top \mathbf{X}$ is singular — ridge handles collinearity and high-dimensional ($p > n$) regimes that OLS cannot. Second, every coefficient is shrunk toward zero, with shrinkage strongest in the directions of small singular values of $\mathbf{X}$ (this is most easily seen via SVD). The penalty parameter $\lambda$ controls the bias–variance trade.

Lasso. Replace the $\ell_2$ penalty with $\ell_1$:

The geometry of the $\ell_1$ ball — a polytope with corners on the coordinate axes — means the optimum frequently lands exactly on those axes, setting some coefficients to exactly zero. Lasso therefore performs variable selection as a side effect of fitting. There is no closed form; the standard solver is coordinate descent or proximal gradient methods, both fast in practice.

Elastic net. A convex combination of the two penalties,

keeps the variable-selection behaviour of lasso while smoothing the path through correlated predictors that lasso treats erratically. It is the default in many high-dimensional applications and the loss surface most ML libraries default to when "$\ell_1$ + $\ell_2$" is offered.

Choosing $\lambda$. Cross-validation. Fit the model for a grid of $\lambda$ values; pick the one with smallest out-of-sample error. The whole pipeline — split, fit, score, repeat — is the empirical core of modern ML and is revisited in Section 13.

The bootstrap. The non-parametric bootstrap replaces the unknown population distribution $F$ with its empirical analogue $\hat{F}_n$ — the discrete distribution that puts mass $1/n$ on each observed datum. Sampling from $\hat{F}_n$ means drawing $n$ observations from the original sample with replacement. Repeat $B$ times to get bootstrap replicates $\hat{\theta}^{*1}, \ldots, \hat{\theta}^{*B}$ of any statistic $\hat{\theta}$.

The empirical distribution of the replicates approximates the sampling distribution of $\hat{\theta}$. Standard errors:

Confidence intervals come in several flavours: the simple percentile interval takes the $\alpha/2$ and $1-\alpha/2$ quantiles of the bootstrap replicates; the basic bootstrap reflects them around $\hat{\theta}$; the BCa (bias-corrected accelerated) interval adjusts for skewness and bias and is the recommended default for modern work. All three are implemented one library call away in R, Python, and Julia.

Why it works. The bootstrap is justified by a single deep observation: $\hat{F}_n \to F$ uniformly (Glivenko–Cantelli), so for "smooth enough" statistics, the sampling distribution under $\hat{F}_n$ converges to the sampling distribution under $F$. The technical conditions (Hadamard differentiability of the functional $\hat{\theta}$ in $F$) hold for almost everything practitioners care about; the famous failures (sample maximum, exact estimators on the boundary) are well documented and have specialised remedies.

The jackknife. An older resampling idea: leave out one observation at a time, recompute $\hat{\theta}_{(-i)}$, and use the variability across the $n$ leave-one-out estimates to estimate $\operatorname{Var}(\hat{\theta})$. The jackknife is computationally cheaper than the bootstrap and conceptually clean; for smooth statistics it gives nearly identical answers and is built into most regression diagnostics (DFFITS, Cook's distance, leave-one-out cross-validation).

Cross-validation. The natural extension to predictive performance. Partition the data into $K$ folds, train on $K-1$, evaluate on the held-out fold, average. The result, $\text{CV}_K = \frac{1}{n}\sum_i L(y_i, \hat{f}_{(-k(i))}(\mathbf{x}_i))$, is a nearly unbiased estimate of the model's expected test error. With $K = n$ this is leave-one-out CV — high variance but minimal bias; with $K = 5$ or 10 a sweet spot for most applications. CV is how every hyperparameter in modern ML — regularisation strengths, network depths, learning rates — gets chosen, and is the practical engine behind the bias–variance trade of Section 12.

The frequentist paradigm of this chapter so far treats $\theta$ as a fixed unknown and randomness as coming from sampling. Probability statements are about procedures: "in 95% of repeated experiments, the interval would contain $\theta$." The Bayesian paradigm treats $\theta$ itself as a random variable with a prior distribution $\pi(\theta)$ encoding beliefs before seeing data. After seeing data $X$, beliefs update via Bayes' theorem to the posterior:

The posterior contains everything the data and prior together imply about $\theta$. From it you read off:

• A point estimate — the posterior mean, median, or mode (the last is the MAP estimate).
• A credible interval $[a, b]$ such that $\mathbb{P}(\theta \in [a, b] \mid X) = 1 - \alpha$ — directly a probability statement about $\theta$.
• A predictive distribution $p(X^{\text{new}} \mid X) = \int p(X^{\text{new}} \mid \theta)\, \pi(\theta \mid X)\, d\theta$ for new data.

Where they agree. With a flat prior and a regular likelihood, the posterior mode equals the MLE and the posterior covariance approaches the inverse Fisher information — so the Bayesian credible interval and the frequentist Wald interval coincide asymptotically. For most practical purposes, MLE-with-Wald-intervals and Bayesian-with-flat-prior return the same numbers; a working data scientist switches between the two phrasings without thinking about it.

Where they disagree. When data are scarce, when the prior is informative, or when the parameter is on the boundary of its space, the two paradigms diverge meaningfully. A Bayesian credible interval for a probability parameter shrinks toward the prior mean (often 0.5) and stays inside $[0, 1]$; a frequentist Wald interval for the same parameter can extend beyond $[0, 1]$ and only catches up asymptotically. Bayesian methods also handle hierarchical structure ("partial pooling" across groups) more naturally and quantify uncertainty in derived quantities through straightforward posterior simulation.

Chapter 07 develops the Bayesian machinery in full — conjugate priors that yield closed-form posteriors, hierarchical models, MCMC and variational inference for everything else. For now the takeaway is that the frequentist methods of this chapter and the Bayesian methods of the next one are not in opposition: they are two complementary languages for the same act of statistical inference, and fluency in both is part of being literate in modern data analysis.

"Correlation does not imply causation" is the slogan; the modern framework that rigorises it goes by the name of potential outcomes (or, in Pearl's formulation, the do-calculus). Let $T \in \{0, 1\}$ denote a treatment and $Y$ an outcome. For each subject define two counterfactual outcomes: $Y(1)$, the outcome they would have had under treatment, and $Y(0)$, the outcome they would have had under control. The fundamental causal estimand is the average treatment effect

The fundamental problem of causal inference is that we observe only one of the two potential outcomes for each subject. The other is missing, and no amount of statistical machinery alone can recover it.

The randomised controlled trial. Random assignment of $T$ severs any dependence between treatment and the potential outcomes: $T \perp \{Y(0), Y(1)\}$. Under this independence,

which is exactly the difference in observed sample means. The RCT is the gold standard not because randomisation is magic but because randomisation makes the missing-data problem solvable. Drug trials, A/B tests, and field experiments in economics all rest on this single identification result.

Observational data. When you cannot randomise — most settings outside controlled experiments — treatment and potential outcomes are not independent. Subjects who choose treatment differ from those who do not in ways that also affect $Y$ (this is confounding). Several frameworks address this:

• Conditioning on confounders. If a set of covariates $\mathbf{X}$ blocks every back-door path between $T$ and $Y$, then $T \perp \{Y(0), Y(1)\} \mid \mathbf{X}$ and $\tau$ is identified by $\mathbb{E}_{\mathbf{X}}[\mathbb{E}[Y \mid T=1, \mathbf{X}] - \mathbb{E}[Y \mid T=0, \mathbf{X}]]$. Regression adjustment, propensity-score weighting, matching, and doubly-robust estimators all fall under this family.
• Instrumental variables. A variable $Z$ that affects $T$ but has no other path to $Y$ — e.g. random encouragement, geographic distance to a service, lottery results — can identify a local treatment effect even with unmeasured confounding.
• Regression discontinuity. When treatment assignment depends on a continuous covariate crossing a threshold (e.g. test-score cutoffs), comparing units just above and just below the threshold gives a clean local causal estimate.
• Difference-in-differences. Comparing the change in outcome over time between a treated group and a control group differences out time-invariant confounders. Standard in policy evaluation.

The relevance to ML is enormous. Predictive models trained on observational data answer rung-one questions ("what is $\mathbb{E}[Y \mid X]$?") perfectly well. They will fail spectacularly if deployed to answer rung-two questions ("what happens if we change $X$?") because they cannot tell correlation from causation by themselves. Recommender systems that recommend books based on past sales conflate "people who already wanted this book bought it" with "showing this book to someone causes them to want it" — and an entire subfield (off-policy evaluation, causal inference for policy learning) exists to disentangle them.

The classical principles, due to Fisher, are three:

• Randomisation. Assign units to treatment groups by a random mechanism. This is what makes the analysis of Section 15 valid and protects against unmeasured confounders.
• Replication. Each treatment is applied to multiple units. Without replication you cannot estimate within-treatment variability and so cannot test hypotheses.
• Blocking. Group experimental units by a known source of variability — sex, age, batch, day of week — and randomise within each block. This removes block-to-block variation from the error term, sharpening the comparison of treatments.

Power and sample size. Before running an experiment, decide how large it needs to be. The power of a test is its probability of detecting a true effect of a given size. For a two-sample t-test detecting a mean difference of $\delta$ with significance level $\alpha$ and power $1 - \beta$, the required per-group sample size is approximately

Halving the detectable effect size quadruples the required sample. Underpowered experiments waste resources and produce a body of literature dominated by inflated effect sizes (the winner's curse: only large estimates survive significance filtering, so published estimates are systematically too large). Power calculations are not optional.

Factorial designs. When studying multiple factors, do not vary them one at a time. A $2^k$ factorial design tests all $2^k$ combinations of $k$ binary factors and lets you estimate main effects and interactions from the same data. With limited resources, a fractional factorial sacrifices some interaction estimates to fit more factors in the same number of runs — the foundation of design of experiments (DoE) in industrial settings.

Bandits and adaptive design. When you can adapt the assignment probabilities as data accumulate, you trade some statistical efficiency for faster identification of the best treatment. Multi-armed bandit algorithms — UCB, Thompson sampling, $\varepsilon$-greedy — are exactly this: experimental designs where the next assignment depends on previous outcomes. They are the bridge from classical experimental design to reinforcement learning, where every action is also a measurement.

A time series $\{X_t\}_{t \in \mathbb{Z}}$ is a sequence of observations indexed by time. The fundamental quantity is the autocovariance

For most useful theory we assume stationarity: the joint distribution of $(X_{t_1}, \ldots, X_{t_k})$ depends only on the spacing of the indices, not on $t$ itself. Under stationarity, $\gamma(h)$ is well-defined and the autocorrelation function $\rho(h) = \gamma(h)/\gamma(0)$ measures linear dependence at lag $h$.

The classical model families. Three workhorses cover most stationary series:

• An autoregressive (AR) process of order $p$: $X_t = \phi_1 X_{t-1} + \cdots + \phi_p X_{t-p} + \varepsilon_t$. The current value is a linear function of the past plus white noise. AR(1) is the textbook starting point and has $\rho(h) = \phi_1^h$.
• A moving-average (MA) process of order $q$: $X_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}$. The current value is a linear function of past shocks.
• An ARMA(p, q) combines both. ARIMA(p, d, q) applies it to the $d$-th difference of a non-stationary series — the bedrock model of econometric forecasting since Box and Jenkins.

Effective sample size. The price of dependence is paid in how much information each new observation contributes. For positively correlated data, the variance of the sample mean is

so the effective sample size $n_{\text{eff}} = n / (1 + 2\sum_h \rho(h))$ can be far smaller than $n$. Standard errors computed assuming independence are then too small, and confidence intervals too narrow. Newey–West and other heteroscedasticity-and-autocorrelation-consistent (HAC) standard errors correct for this in regression on time-series data; cluster-robust standard errors do the analogous job for grouped data.

Modern time-series ML — recurrent networks, temporal convolutions, transformers with positional encodings, state-space models — replaces the linear AR/MA structure with arbitrary non-linear function classes. The statistical principles (stationarity, effective sample size, sequential evaluation) remain. Most published headline results on forecasting that compare deep models to ARIMA do so on truly stationary, well-behaved series; the picture is murkier on the messy real-world data where Box and Jenkins still hold up surprisingly well.

A working list of correspondences worth memorising:

• Empirical risk minimisation. The basic setup of supervised learning — minimise $\frac{1}{n}\sum_i L(y_i, f_\theta(\mathbf{x}_i))$ — is exactly the extremum estimator framework of econometrics, with all its asymptotic theory (consistency, asymptotic normality, sandwich variance) directly applicable.
• Cross-entropy = NLL. Categorical cross-entropy is the negative log-likelihood of a multinomial model. Binary cross-entropy is the NLL of a Bernoulli model (i.e., logistic regression). Mean-squared error is the NLL of a Gaussian. Picking a "loss function" is choosing a likelihood.
• Regularisation = MAP. Weight decay is a Gaussian prior. $\ell_1$ regularisation is a Laplace prior. Dropout has a variational interpretation as approximate Bayesian inference. Every trick used to fight overfitting has a Bayesian shadow.
• Train/validation/test = sample partitioning. The discipline of holding out a test set and touching it once is exactly the discipline of avoiding the multiple-testing problem of Section 08 across a hyperparameter sweep.
• Cross-validation = honest test error. The CV estimator is a nearly unbiased estimator of expected test loss; its variance behaves as Section 13 predicts. Reporting "CV ± SE" is reporting a sampling-distribution-based confidence interval.
• Calibration. A binary classifier is calibrated if among instances assigned probability $p$ the empirical base rate is $p$. The reliability diagram is a goodness-of-fit plot; the Brier score and ECE are scoring rules with statistical pedigrees.
• Conformal prediction. A modern, distribution-free framework for producing prediction intervals with guaranteed coverage from any black-box predictor. It is essentially a clever application of exchangeability and empirical quantiles — Section 13 ideas in modern dress.
• Bandits and reinforcement learning. Exploration vs. exploitation is the experimental-design vs. estimation trade-off of Section 16. Off-policy evaluation is exactly the causal-inference identification problem of Section 15.
• A/B testing infrastructure. Every large internet company runs experiments at scale using the framework of Sections 06–08, with extensions for sequential monitoring and heterogeneous effects.
• Statistical learning theory. Rademacher complexity and VC dimension generalise the bias–variance decomposition to data-dependent function classes; PAC bounds turn concentration inequalities (Section 16 of the previous chapter!) into generalisation guarantees. The whole apparatus is statistics applied to the algorithms ML practitioners run every day.
• Causal ML. Doubly robust estimators, double machine learning, instrumental-variable forests — modern econometric methods built directly on the Section 15 framework, run with neural networks and gradient boosting in place of linear regression.
• Diagnostics and ablations. A residual plot in regression, a learning curve in deep learning, an ablation in a paper — all are model-criticism procedures whose statistical lineage is the diagnostic toolkit of Section 10.

The pattern is consistent. Almost every "advance" in ML is a familiar statistical idea applied at a scale or with a function class that earlier statisticians could not access. Knowing the older language does not make you slower; it makes you faster, because you stop reinventing things and start composing them.

Keep these connections in mind as you read the rest of the compendium. The remaining chapters of Part I — information theory, Bayesian reasoning, signal processing — extend this scaffolding rather than replace it; Part II then uses the whole structure to build modern ML from the ground up.