Statistical estimation

Maximum likelihood and MAP, what makes an estimator good, the bias-variance decomposition, and interval estimates computed by bootstrap.

Maximum likelihood and MAP

Maximum likelihood picks the parameters under which the observed data is most probable: argmax P(data | theta). Maximum a posteriori adds a prior: argmax P(data | theta) P(theta). Under a Gaussian prior on the weights, MAP is exactly L2-regularised maximum likelihood — this is the source of weight decay.

import numpy as np
from scipy.optimize import minimize_scalar

data = np.array([2.1, 1.8, 2.4, 2.0, 1.9, 5.2])   # one outlier

def neg_log_lik_gaussian(mu, sigma=1.0):
    return -np.sum(-0.5 * ((data - mu) / sigma) ** 2 - np.log(sigma * np.sqrt(2 * np.pi)))

mle = minimize_scalar(neg_log_lik_gaussian, bounds=(-10, 10), method="bounded").x
mle, data.mean()        # identical: the mean is the MLE for a Gaussian mean

# MAP with a Gaussian prior centred at 0 = ridge shrinkage toward 0
def neg_log_post(mu, prior_sigma=2.0):
    return neg_log_lik_gaussian(mu) + 0.5 * (mu / prior_sigma) ** 2

map_est = minimize_scalar(neg_log_post, bounds=(-10, 10), method="bounded").x
round(map_est, 4), round(mle, 4)      # shrunk toward zero

• MLE is consistent and asymptotically efficient, but can be badly biased in small samples (the MLE for a Gaussian variance divides by n, not n-1).
• A prior is not cheating: it is a statement about plausible parameter values, and it is what makes a model identifiable when data is scarce.
• Regularisation strength is prior strength. L2 weight decay is a Gaussian prior; L1 is a Laplace prior and produces sparsity.
• Outliers move the mean because squared error is dominated by them. A heavier-tailed likelihood (Huber, Student-t) is a modelling choice, not a trick.

Bias, variance and the decomposition

rng = np.random.default_rng(0)

def true_fn(x):
    return np.sin(2 * np.pi * x)

def fit_polynomial(x, y, degree):
    return np.polyfit(x, y, degree)

degrees = [1, 3, 15]
grid = np.linspace(0, 1, 60)

for deg in degrees:
    preds = []
    for _ in range(50):
        xs = rng.random(15)
        ys = true_fn(xs) + 0.2 * rng.normal(size=15)
        preds.append(np.polyval(fit_polynomial(xs, ys, deg), grid))
    preds = np.array(preds)
    bias2 = ((preds.mean(axis=0) - true_fn(grid)) ** 2).mean()
    variance = preds.var(axis=0).mean()
    print(f"degree={deg:2d} bias^2={bias2:.4f} var={variance:.4f} total={bias2 + variance:.4f}")

Confidence intervals and bootstrap

def bootstrap_ci(values, statistic=np.mean, n=5000, alpha=0.05):
    rng = np.random.default_rng(0)
    values = np.asarray(values, float)
    stats = np.array([statistic(rng.choice(values, size=len(values), replace=True))
                      for _ in range(n)])
    lo, hi = np.quantile(stats, [alpha / 2, 1 - alpha / 2])
    return float(lo), float(hi)

accuracy = np.array([1, 1, 0, 1, 0] * 40)     # 60% on 200 examples
bootstrap_ci(accuracy, n=2000)                # roughly (0.53, 0.67)

• A confidence interval describes the sampling variability of the estimate, not the probability that the true value lies in this particular interval. The distinction matters when you are deciding whether two models differ.
• For a metric like accuracy the standard error is sqrt(p(1-p)/n); with 200 examples and p = 0.6 that is about 0.035, so a 5-point difference is not significant.
• Bootstrap the whole evaluation, including sampling of examples, to test whether a model comparison is real.
• Paired comparisons are far more sensitive than comparing two independent intervals: evaluate both models on the same examples and bootstrap the difference.

FAQ

Related

Information theory for machine learning Evaluation and overfitting

Last refreshed 2026-09-18.