estimate_1 estimate_2
10.571572 9.980585 21 Estimation
You have data and a model with unknown parameters. Estimation is the business of choosing values for them, and of saying how much to trust the choice.
Two questions run through the chapter. How do we pick an estimator? — and there are several principled answers. How do we judge one? — where the honest answer turns out to be more subtle than “is it unbiased”, and leads directly to the tradeoff that governs all of machine learning.
21.1 Estimators
An estimator is any function of the data used to guess a parameter:
\[ \hat{\theta} = T(X_1, \dots, X_n) \]
The hat marks an estimate. Note carefully: \(\theta\) is a fixed unknown number, while \(\hat{\theta}\) is a random variable — it depends on the sample, and a different sample gives a different value. Everything in this chapter is about the distribution of that random variable, called the sampling distribution.
Same estimator, same population, different answers. The spread of those answers is what “uncertainty in the estimate” means.
21.2 Bias, variance, and mean squared error
Two ways an estimator can be wrong:
\[ \operatorname{Bias}(\hat{\theta}) = \mathbb{E}[\hat{\theta}] - \theta \qquad \operatorname{Var}(\hat{\theta}) = \mathbb{E}\bigl[(\hat{\theta} - \mathbb{E}[\hat{\theta}])^2\bigr] \]
Bias is systematic error — aiming at the wrong place. Variance is instability — scattering around wherever you aim. Mean squared error combines them:
\[ \operatorname{MSE}(\hat{\theta}) = \mathbb{E}\bigl[(\hat{\theta}-\theta)^2\bigr] = \operatorname{Bias}(\hat{\theta})^2 + \operatorname{Var}(\hat{\theta}) \tag{21.1}\]
An estimator with zero bias is unbiased. The classic example is the sample variance, where the divisor matters:
\[ \frac{1}{n}\sum_i (x_i - \bar{x})^2 \quad\text{is biased}; \qquad \frac{1}{n-1}\sum_i (x_i - \bar{x})^2 \quad\text{is not} \]
Dividing by \(n\) understates the spread, because \(\bar{x}\) is itself fitted to the data and sits closer to the points than the true mean does. The correction is exact: \(\mathbb{E}[\hat{\sigma}^2_n] = \frac{n-1}{n}\sigma^2\).
over_n over_n1
mean 3.190612 3.988265
theory 3.200000 4.000000
WarningWatch out
Unbiased does not mean better. Compare the two estimators by MSE rather than bias:
over_n over_n1
simulated 5.756259 7.970684
theory 5.760000 8.000000The biased estimator has the lower MSE. Dividing by \(n\) shrinks the estimate slightly, and the variance saved by shrinking outweighs the bias introduced.
This is not a curiosity — it is the whole logic of regularization (Section 6.11). Ridge regression is deliberately biased and often more accurate for exactly this reason.
21.3 The method of moments
The oldest recipe, and the simplest: set sample moments equal to theoretical moments and solve.
To estimate one parameter, match the mean. For an exponential distribution \(\mathbb{E}[X] = 1/\lambda\), so setting \(\bar{x} = 1/\lambda\) gives \(\hat{\lambda} = 1/\bar{x}\).
sample_mean lambda_hat truth
1.8427354 0.5426715 0.5000000 For two parameters, match the first two moments. It is easy, requires no calculus, and usually gives a reasonable answer — but it can produce estimates outside the valid parameter range, and it generally has larger variance than the method in the next section.
21.4 Maximum likelihood
The dominant approach. Ask: which parameter value would make the data we actually saw most probable?
The likelihood is the joint density of the observed data, viewed as a function of the parameter:
\[ L(\theta) = \prod_{i=1}^{n} f(x_i \mid \theta) \tag{21.2}\]
and the maximum likelihood estimator is the \(\theta\) maximizing it.
The reversal of viewpoint is the key. The same expression, read as a function of \(x\) with \(\theta\) fixed, is a density; read as a function of \(\theta\) with \(x\) fixed, it is a likelihood. A likelihood is not a probability distribution over \(\theta\) — it does not integrate to 1, and treating it as one is the mistake Bayesian methods avoid by introducing a prior.
21.5 The log-likelihood
Nobody maximizes Equation 21.2 directly. Take logarithms:
\[ \ell(\theta) = \log L(\theta) = \sum_{i=1}^{n} \log f(x_i \mid \theta) \tag{21.3}\]
Three reasons, all decisive:
- Products become sums (Section 2.4), and sums are far easier to differentiate.
- Numerical survival. A product of 1,000 densities each around \(0.1\) underflows to zero in double precision; its logarithm is a comfortable \(-2303\).
- Same answer. \(\log\) is strictly increasing, so it does not move the maximum.
product log_sum
0.000 -921.034 The product underflows to exactly zero — every parameter value would look equally good. The log-likelihood is perfectly well behaved.
Worked example. For an exponential with rate \(\lambda\):
\[ \ell(\lambda) = n\log\lambda - \lambda\sum_i x_i \]
Differentiating and setting to zero gives \(n/\lambda = \sum_i x_i\), so \(\hat{\lambda} = 1/\bar{x}\) — the same answer the method of moments gave.
argmax_numeric formula
0.5433584 0.5426715 draw_fun(
function(l) loglik_exp(l, xe),
from = 0.25, to = 1.1,
xlab = "lambda",
ylab = "log-likelihood"
)
NoteIn machine learning
Most loss functions are negative log-likelihoods. Squared error is the negative log-likelihood of a normal with constant variance; binary cross-entropy is the Bernoulli one; Poisson deviance is the Poisson one.
So “minimize the loss” and “maximize the likelihood” are usually the same instruction, and choosing a loss is implicitly choosing a noise model — which is worth knowing when your residuals are obviously not normal.
21.6 Fisher information
How much does the data tell you about \(\theta\)? The answer is curvature of the log-likelihood:
\[ I(\theta) = -\mathbb{E}\left[\frac{\partial^2 \ell}{\partial\theta^2}\right] \tag{21.4}\]
A sharply peaked log-likelihood means nearby parameter values fit much worse, so the data are informative. A flat one means many values fit about equally, and the estimate is poorly determined. This is the second derivative of Section 11.7 doing statistical work — and in several parameters it is the Hessian (Section 12.8).
The payoff is the Cramér–Rao bound: no unbiased estimator can do better than
\[ \operatorname{Var}(\hat{\theta}) \geq \frac{1}{n\,I(\theta)} \tag{21.5}\]
Worked example. For a Bernoulli, \(I(p) = 1/\bigl(p(1-p)\bigr)\), so the bound is \(p(1-p)/n\) — which is exactly the variance of \(\bar{x}\). The sample proportion is as good as any unbiased estimator can be.
fisher_info cramer_rao_bound simulated_variance
4.761904762 0.002100000 0.002089889 The simulated variance sits on the bound. Note this also connects to conditioning: a flat log-likelihood is an ill-conditioned Hessian, and the same near-singularity that made Section 6.10 hard to solve makes parameters hard to estimate.
21.7 Maximum a posteriori estimation
Maximum likelihood uses only the data. MAP estimation adds a prior via Bayes’ theorem (Equation 18.6):
\[ \hat{\theta}_{\text{MAP}} = \arg\max_\theta \; \bigl[\log p(x \mid \theta) + \log p(\theta)\bigr] \tag{21.6}\]
The evidence term is constant in \(\theta\), so it drops out of the maximization.
Read Equation 21.6 as a penalized likelihood and the connection to regularization is immediate:
| Prior on \(\boldsymbol{\beta}\) | Log-prior contributes | Equivalent penalty |
|---|---|---|
| normal, variance \(\tau^2\) | \(-\|\boldsymbol{\beta}\|^2/(2\tau^2)\) | ridge, \(\lambda = \sigma^2/\tau^2\) |
| Laplace | \(-\|\boldsymbol{\beta}\|_1/b\) | Lasso |
| flat | nothing | plain maximum likelihood |
Ridge regression is MAP estimation with a normal prior. The regularization strength \(\lambda\) is a statement about how tightly you believe the coefficients cluster around zero, and a stronger prior means a larger \(\lambda\).
21.8 Bayesian estimation
MAP keeps only the peak of the posterior. Full Bayesian estimation keeps the whole distribution:
\[ p(\theta \mid x) \propto p(x \mid \theta)\,p(\theta) \]
and reports a summary — the posterior mean, say — along with a credible interval covering, for instance, 95% of the posterior mass.
The Beta–binomial pair from Section 20.7 makes this concrete, because the posterior stays in the same family:
# unname() throughout: arithmetic on a named vector keeps the
# name, and c(a = <named>) would nest it into "a.a"
prior_a <- 2 # mild belief centered at 0.5
prior_b <- 2
s <- 7 # successes
f <- 3 # failures
post_a <- prior_a + s
post_b <- prior_b + f
c(
prior_mean = prior_a / (prior_a + prior_b),
posterior_mean = post_a / (post_a + post_b),
mle = s / (s + f),
credible_lower = qbeta(0.025, post_a, post_b),
credible_upper = qbeta(0.975, post_a, post_b)
) prior_mean posterior_mean mle credible_lower credible_upper
0.5000000 0.6428571 0.7000000 0.3857383 0.8614207 The posterior mean sits between the prior mean and the MLE — shrinkage again. With more data the likelihood dominates and the prior fades; with little data the prior does more of the work, which is exactly what you want.
WarningWatch out
A credible interval and a confidence interval answer different questions, and the difference is not pedantic.
A 95% credible interval genuinely has 95% posterior probability of containing \(\theta\) — the interpretation people naturally want. A confidence interval (Section 22.5) makes a statement about the procedure across repeated samples, not about this particular interval.
They often coincide numerically, which is why the distinction gets blurred, but the statements are different.
21.9 The bootstrap
What if you have no formula for the standard error? For a mean there is one; for a median, a correlation, or a ratio of quantiles there may not be.
The bootstrap replaces the formula with resampling: treat the sample as if it were the population, draw new samples from it with replacement, and watch how the estimate varies.
estimate se_formula se_bootstrap ci_lower.2.5% ci_upper.97.5%
9.0411848 0.3674496 0.3573466 8.3718011 9.7484856 draw_histogram(boot, xlab = "bootstrap sample mean")The bootstrap standard error matches the analytic one closely — and for the mean we did not need it. The point is that the same three lines work for a statistic with no formula at all.
median se_bootstrap
8.9558377 0.4888124
WarningWatch out
The bootstrap is not magic. It assumes your sample represents the population, so it cannot fix a biased sample — resampling bad data gives a confident answer to the wrong question.
It also struggles with statistics that depend on extremes, such as the maximum: the resampled maximum can never exceed the observed one, so the bootstrap distribution is systematically wrong in the tail.
21.10 The bias-variance tradeoff
Equation 21.1 splits error into bias and variance, and the central fact of predictive modeling is that you generally cannot reduce both at once.
- A simple model cannot capture the truth: high bias, low variance.
- A complex model chases noise in the particular sample: low bias, high variance.
The total is minimized somewhere between, and finding that point is what model selection is for.
set.seed(12)
truth_f <- function(x) sin(2 * x)
one_run <- function(degree) {
xtr <- seq(-1, 1, length.out = 20)
ytr <- truth_f(xtr) + rnorm(20, sd = 0.3)
fit <- lm(ytr ~ poly(xtr, degree, raw = TRUE))
xte <- seq(-1, 1, length.out = 50)
pred <- predict(fit, data.frame(xtr = xte))
mean((pred - truth_f(xte))^2)
}
degrees <- 1:9
err <- sapply(degrees, function(d) {
mean(replicate(200, one_run(d)))
})
round(err, 4)[1] 0.0376 0.0426 0.0173 0.0205 0.0234 0.0298 0.0320 0.0377 0.0447draw_line(
x = degrees,
y = err,
xlab = "polynomial degree",
ylab = "squared error to the truth"
)The U-shape is the tradeoff made visible. On the left the model is too rigid to follow \(\sin(2x)\); on the right it has enough freedom to interpolate the noise, and every new sample sends it somewhere different.
The curve is not perfectly smooth, and the exception is instructive. Degree 2 is worse than degree 1, reliably. \(\sin(2x)\) is an odd function, so its expansion contains only odd powers — the \(x^2\) term has nothing real to fit and spends a degree of freedom on noise. Degree 3 is the first that adds a term the truth actually has, and the error drops sharply there.
So complexity is not one-dimensional. The right question is not how many parameters a model has but whether they match the structure of the problem, which is what choosing a model class is for.
NoteIn machine learning
This curve is why we have validation sets, cross-validation, early stopping and regularization — all of them are ways of locating the bottom without cheating by looking at the test data.
Modern over-parameterized networks complicate the picture: error can fall again past the interpolation point, a phenomenon called double descent. The tradeoff is real but the simple U-shape is not the whole story.
21.11 Summary
| Concept | Statement |
|---|---|
| Estimator | a function of the data; itself random |
| Bias | \(\mathbb{E}[\hat\theta] - \theta\) |
| MSE | \(\text{bias}^2 + \text{variance}\) |
| Method of moments | match sample and theoretical moments |
| Maximum likelihood | maximize \(\prod f(x_i\mid\theta)\) |
| Log-likelihood | \(\sum \log f(x_i\mid\theta)\); sums, and no underflow |
| Fisher information | curvature of \(\ell\); more curvature, more precision |
| Cramér–Rao | \(\operatorname{Var}(\hat\theta) \geq 1/(nI(\theta))\) |
| MAP | likelihood \(\times\) prior; ridge is the normal-prior case |
| Bootstrap | resample to get a standard error without a formula |
21.12 Exercises
1. Show that \(\bar{X}\) is unbiased for \(\mu\), and find its variance.
TipSolution
By linearity of expectation (Equation 19.3), \(\mathbb{E}[\bar{X}] = \frac{1}{n}\sum_i \mathbb{E}[X_i] = \mu\). For independent observations the variances add, so \(\operatorname{Var}(\bar{X}) = \frac{1}{n^2}\cdot n\sigma^2 = \sigma^2/n\).
mean_of_means target var_of_means theory
5.003298 5.000000 0.159497 0.160000 The \(\sigma^2/n\) is why precision improves only as \(\sqrt{n}\): to halve the standard error you need four times the data.
2. Derive the MLE for a Bernoulli sample.
TipSolution
With \(s\) successes in \(n\) trials, \(\ell(p) = s\log p + (n-s)\log(1-p)\). Differentiating:
\[ \frac{s}{p} - \frac{n-s}{1-p} = 0 \;\Longrightarrow\; \hat{p} = \frac{s}{n} \]
numeric_argmax formula
0.3403808 0.3400000 The sample proportion, which is also what the method of moments gives — the two agree here, though they often do not.
3. For the exponential, compare the method of moments and MLE estimates.
TipSolution
method_of_moments mle identical
0.5426715 0.5426715 1.0000000 They are the same estimator here. That is a feature of the exponential, not a general rule — for most families the two differ, and the MLE usually has smaller variance.
4. Bootstrap the standard error of a correlation, which has an awkward analytic formula.
TipSolution
estimate se_bootstrap ci_lower.2.5% ci_upper.97.5%
0.70607712 0.06374402 0.56434461 0.81307809 Resample pairs, not each variable separately — resampling them independently would destroy the very association being measured. That is the most common bootstrap mistake in multivariate settings.
5. Show that a normal prior on a coefficient gives a ridge penalty.
TipSolution
With \(y_i \sim \mathcal{N}(\mathbf{x}_i^\top\boldsymbol{\beta}, \sigma^2)\) and \(\boldsymbol{\beta} \sim \mathcal{N}(\mathbf{0}, \tau^2\mathbf{I})\), the log-posterior is
\[ -\frac{1}{2\sigma^2}\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\|^2 -\frac{1}{2\tau^2}\|\boldsymbol{\beta}\|^2 + \text{const} \]
Maximizing is minimizing \(\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\|^2 + \lambda\|\boldsymbol{\beta}\|^2\) with \(\lambda = \sigma^2/\tau^2\), which is Equation 6.3.
map ols
[1,] 1.062137 1.066999
[2,] 1.981417 1.996320A tight prior (small \(\tau\)) means a large \(\lambda\) and heavy shrinkage; a diffuse prior recovers ordinary least squares. The regularization strength is a prior belief, which is a more satisfying answer than “a hyperparameter to tune”.
6. Using the bias-variance simulation, find the degree minimizing error and explain the shape on each side.
TipSolution
best_degree best_error degree_1 degree_9
3.0000 0.0173 0.0376 0.0447 Degree 3 is best. Below it the model cannot bend enough to follow \(\sin(2x)\) — that is bias, and it would persist however much data you collected. Above it the extra coefficients fit whatever noise the particular sample contains — that is variance, and it would shrink with more data.
The distinction matters practically: high bias calls for a better model, high variance calls for more data or more regularization. Diagnosing which you have is most of what learning curves are for.