Maximum Likelihood Estimation

Maximum likelihood estimation (MLE) chooses the parameters of a probabilistic model so that the data actually observed are as probable as possible under the model. It is the most common general recipe for turning a model of how images are produced into an estimator: the noise level of a camera, the rate of a photon count, the colour distribution of an object, the translation between two frames, the pose of a camera or the weights of a classifier. Many of the losses used in vision are maximum-likelihood objectives in disguise. Least squares is maximum likelihood under Gaussian noise, the sum of absolute differences under Laplace noise, and the cross-entropy loss of a classifier is the negative log-likelihood of its labels. This article defines the likelihood, works through the standard cases, states the properties that make the method attractive and the conditions they need, and shows where it breaks down. It builds on Probability and Statistics, which introduces estimators, their bias and variance, and Bayes’ rule.

Definition

A statistical model is a family of distributions $p(x \mid \theta)$ indexed by parameters $\theta$. Once data $\mathcal{D}$ have been observed, the probability (or probability density) of those data, read as a function of $\theta$ with $\mathcal{D}$ held fixed, is the likelihood function. The maximum likelihood estimate is the parameter value at which the likelihood is largest. In practice its logarithm, the log-likelihood, is maximized instead, or equivalently the negative log-likelihood (NLL) is minimized.

Intuition

Photograph a flat grey card and look at one row of pixels. Under a Gaussian noise model, each candidate pair of values for the true grey level and the noise standard deviation assigns some probability density to the readings actually recorded. A grey level far from the readings makes them improbable, since every pixel would then have needed a large noise excursion. A noise level that is too small makes the scatter of the readings improbable, and one that is too large spreads the density so thinly that every reading becomes improbable. The maximum likelihood estimate is the pair that makes the recorded row the least surprising: here the average of the readings and their root-mean-square spread about it.

The likelihood scores parameters, not data, and it is not a probability that the parameters are correct. The same numbers that, for a fixed parameter, form a distribution over possible images form, for a fixed image, a function of the parameters that need not sum or integrate to one (Bishop, 2006, §1.2.3). Fisher, who named the method, insisted on the distinction: likelihood is assigned to a point of the parameter range and, unlike a probability density, cannot be integrated over it (Fisher, 1922, p. 327).

Formal Definition

The notation follows Probability and Statistics: $p(x)$ for a probability mass or density function and $P(\cdot)$ for the probability of an event.

Likelihood and log-likelihood

Let $\mathcal{D} = \{x_1, \dots, x_n\}$ be observations assumed independent and identically distributed (i.i.d.) with density or mass function $p(x \mid \theta)$: the pixel values of a flat patch, the colours of an object’s pixels, the photon counts of a sensor region. The likelihood and log-likelihood are

$$L(\theta) = p(\mathcal{D} \mid \theta) = \prod_{i=1}^{n} p(x_i \mid \theta), \qquad \ell(\theta) = \log L(\theta) = \sum_{i=1}^{n} \log p(x_i \mid \theta),$$

and the maximum likelihood estimate is

$$\hat{\theta}_{\text{ML}} = \arg\max_{\theta} \ell(\theta) = \arg\min_{\theta} \Big( -\sum_{i=1}^{n} \log p(x_i \mid \theta) \Big).$$

Because the logarithm is increasing, maximizing $\ell$ and maximizing $L$ give the same $\theta$. The logarithm turns the product into a sum of per-pixel or per-example terms, which is easier to differentiate, and it avoids the numerical underflow that a product of many small probabilities causes (Bishop, 2006, §1.2.4). For supervised models, such as a classifier mapping an image $\mathbf{x}_i$ to a label $y_i$, the same definition applies to the conditional likelihood $\prod_i p(y_i \mid \mathbf{x}_i, \theta)$ (Murphy, 2022, §4.2.1).

Worked cases

When the log-likelihood is differentiable and has a single interior maximum, the estimate solves $\nabla_\theta \ell(\theta) = \mathbf{0}$. Several models used throughout vision have closed-form solutions.

Bernoulli: a binary pixel or outcome. If each $x_i \in \{0, 1\}$ (a pixel is foreground or not, a detection is correct or not) with $P(x_i = 1) = \theta$, then $\ell(\theta) = n_1 \log \theta + n_0 \log(1 - \theta)$, where $n_1$ and $n_0$ count the ones and zeros. Setting the derivative to zero gives the observed fraction (Murphy, 2022, §4.2.3):

$$\hat{\theta} = \frac{n_1}{n_0 + n_1}.$$

Categorical: a histogram. If each pixel’s quantized colour or grey level falls in one of $K$ bins with probabilities $\theta_1, \dots, \theta_K$ summing to one, and bin $k$ receives $N_k$ of the $N$ pixels, then $\ell(\boldsymbol{\theta}) = \sum_k N_k \log \theta_k$. Maximizing under the sum-to-one constraint with a Lagrange multiplier gives the relative frequencies (Murphy, 2022, §4.2.4):

$$\hat{\theta}_k = \frac{N_k}{N}.$$

The normalized image histogram of a region is therefore the maximum likelihood estimate of a categorical model of its pixel values. A bin with no pixels gets probability exactly zero, which causes trouble on new data (see the third example below).

Poisson: photon counts. The number of photoelectrons a pixel collects during an exposure is modelled as Poisson with rate $\lambda$ (Image Noise). For counts $k_1, \dots, k_n$ from pixels that share the same expected value, a flat region of a raw image, the log-likelihood and its derivative are

$$\ell(\lambda) = \sum_{i=1}^{n} \big( k_i \log \lambda - \lambda - \log k_i! \big), \qquad \frac{d\ell}{d\lambda} = \frac{1}{\lambda} \sum_{i=1}^{n} k_i - n,$$

so $\hat{\lambda} = \frac{1}{n} \sum_i k_i$, the mean count (a one-line derivation).

Gaussian: mean and variance. For pixel values $x_i \sim \mathcal{N}(\mu, \sigma^2)$, the Gaussian log-likelihood is

$$\ell(\mu, \sigma^2) = -\frac{1}{2\sigma^2} \sum_{i=1}^{n} (x_i - \mu)^2 - \frac{n}{2} \log \sigma^2 - \frac{n}{2} \log 2\pi,$$

and maximizing it gives the sample mean and the variance with denominator $n$ (Bishop, 2006, §1.2.4; Murphy, 2022, §4.2.5):

$$\hat{\mu} = \frac{1}{n} \sum_{i=1}^{n} x_i, \qquad \hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \hat{\mu})^2 .$$

The variance estimate is biased low by the factor $(n-1)/n$, which is why the unbiased estimator divides by $n - 1$; Probability and Statistics derives the bias and shows on $3 \times 3$ patches that the biased estimator can still have the smaller mean squared error. For a $d$-dimensional Gaussian, such as the RGB colours of a region, the estimates are the mean vector and the covariance matrix with denominator $n$ (Murphy, 2022, §4.2.6). NumPy’s np.var returns the maximum likelihood variance by default, while np.cov divides by $n - 1$ unless called with bias=True (NumPy 2.3.5).

Properties

Noise models become losses

Many estimation problems in vision have the form “measurement = model prediction + noise”: observed pixels are a warped template plus noise, detected keypoints are projected 3D points plus noise, regressed values are a network’s output plus noise. Writing $r_i(\theta) = y_i - f_i(\theta)$ for the residual of measurement $i$, the noise distribution fixes the negative log-likelihood, and with it the loss.

• Gaussian noise, equal variance: $-\ell(\theta) = \frac{1}{2\sigma^2} \sum_i r_i(\theta)^2 + \frac{n}{2} \log 2\pi\sigma^2$, so the maximum likelihood $\theta$ minimizes the sum of squared residuals, whatever $\sigma$ is (Bishop, 2006, §1.2.5). This is the argument that makes least squares and nonlinear least squares maximum likelihood estimators, and makes the minimum of the reprojection error in bundle adjustment the maximum likelihood reconstruction under i.i.d. Gaussian image noise. Maximizing over $\sigma^2$ as well gives $\hat{\sigma}^2 = \frac{1}{n} \sum_i r_i(\hat{\theta})^2$, the mean squared residual (Bishop, 2006, §1.2.5).
• Gaussian noise, known unequal variances: with $r_i \sim \mathcal{N}(0, \sigma_i^2)$, the NLL is $\sum_i r_i^2 / (2\sigma_i^2)$ plus terms that do not depend on $\theta$, which is weighted least squares with weights $1 / \sigma_i^2$ (Murphy, 2022, §11.2.2.4). Under the Poisson–Gaussian model of a raw sensor, whose variance grows with intensity, bright pixels deserve less weight than dark ones.
• Laplace noise: with density proportional to $\exp(-|r| / b)$, the NLL is $\frac{1}{b} \sum_i |r_i|$ plus a constant: the sum of absolute residuals. Its heavier tails make the estimate far less sensitive to outliers (Murphy, 2022, §11.6.1). For estimating one location from repeated measurements, the squared loss gives the mean and the absolute loss the median, as the burst-denoising example in Loss Function shows.
• Categorical labels: if a classifier’s softmax output is the probability of each class, the NLL of the true labels is the cross-entropy loss, summed over images, or over pixels for segmentation. Logistic regression is the linear model trained this way.

More generally, minimizing the NLL averaged over the data is the same as minimizing the Kullback–Leibler divergence from the empirical distribution of the data to the model (Murphy, 2022, §4.2.2): maximum likelihood picks the member of the model family that is closest, in this sense, to the data actually seen. Replacing the log-loss by an arbitrary per-example loss gives empirical risk minimization, of which maximum likelihood is the special case (Murphy, 2022, §4.3).

Large-sample behaviour

The case for maximum likelihood rests on how it behaves as the amount of data grows, provided the model family contains the distribution that generated the data, its parameters are identifiable, meaning that different parameter values give different distributions (Murphy, 2022, §5.3.2), and the model satisfies regularity conditions of smoothness that Murphy leaves unstated, giving the asymptotic results only “for certain models” (Murphy, 2022, §4.7.2).

• Consistency. The estimate converges in probability to the true parameter as $n \to \infty$ (Murphy, 2022, §5.3.2). Fisher’s own 1922 criterion of consistency was a different notion: a statistic computed from the whole population should equal the parameter (Fisher, 1922, p. 309).
• Asymptotic normality. For large $n$, the sampling distribution of the estimate approaches a Gaussian centred on the true value $\theta^*$ (Murphy, 2022, §4.7.2):

$$\hat{\theta}_{\text{ML}} \;\approx\; \mathcal{N}\!\left(\theta^*,\, \big(n F(\theta^*)\big)^{-1}\right), \qquad F(\theta) = \mathbb{E}\!\left[\nabla_\theta \log p(x \mid \theta)\, \nabla_\theta \log p(x \mid \theta)^\top\right],$$

where $F$ is the Fisher information of one observation. Under regularity conditions it also equals the expected Hessian of the per-observation NLL, $-\mathbb{E}[\nabla_\theta^2 \log p(x \mid \theta)]$ (Murphy, 2022, §4.7.2). The curvature of the log-likelihood at its peak therefore measures how well the data determine the parameters: a sharp peak means a small variance. In practice the Hessian of the NLL at the estimate (the observed information) is inverted to give standard errors. Geometric fits such as camera calibration and bundle adjustment obtain their optional parameter covariances in this spirit, from the curvature of the reprojection error at its minimum. The first example below checks such standard errors against repeated experiments.

• Efficiency. The Cramér–Rao bound states that, under smoothness assumptions on the model, any unbiased estimator has variance at least $1 / (n F(\theta^*))$ (Murphy, 2022, §4.7.6.2). The maximum likelihood estimate reaches this bound asymptotically, so it has the smallest asymptotic variance of any unbiased estimator (Murphy, 2022, §4.7.6.2). For the Poisson rate, $F(\lambda) = 1/\lambda$, so the bound is $\lambda / n$, exactly the variance of the mean count: averaging a flat region’s photon counts is already the best unbiased estimate of its brightness (a short calculation from the definitions above).

None of these properties is promised for small $n$, for misspecified models or for unidentifiable parameters, the cases listed under Limitations below.

Invariance under reparameterization

If $\hat{\theta}$ is the maximum likelihood estimate of $\theta$ and $g$ is a one-to-one function, then $g(\hat{\theta})$ is the maximum likelihood estimate of $g(\theta)$: relabelling the parameter does not move the peak of the likelihood. The maximum likelihood estimate of the noise standard deviation is the square root of the maximum likelihood variance. Fisher made this point in 1922 to separate his method from the mode of a posterior distribution, which moves when the scale on which the parameter is measured changes, whereas the position of maximum likelihood is unchanged by any such transformation (Fisher, 1922, p. 327). Unbiasedness, by contrast, is not preserved: by Jensen’s inequality, the square root of an unbiased variance estimate underestimates the standard deviation on average.

Numerical maximization

Closed forms exist only for simple models. When the parameters enter the likelihood nonlinearly, as for the translation of an image, the pose of a camera, a noise variance that depends on the signal or the weights of a network, the NLL is minimized numerically. Gaussian noise with a nonlinear model gives nonlinear least squares, solved by Gauss–Newton or Levenberg–Marquardt iterations; other likelihoods are minimized with gradient descent, stochastic gradient descent for large data sets, or quasi-Newton methods, the main solver families of optimization. A Gaussian mixture model has no closed-form solution because a sum over components sits inside the logarithm (Bishop, 2006, §9.2.1); expectation–maximization is the standard alternative there, alternating between computing soft assignments of the data to components and re-estimating each component in closed form.

Numerical maximization brings the usual caveats: the log-likelihood may have several local maxima, so the starting point matters, and a cheap estimate such as the method of moments, which matches sample moments to model moments, is often used to initialize the iterations (Murphy, 2022, §4.4.1). Constrained parameters such as variances and probabilities are kept valid by bounds or by optimizing their logarithm or softmax logits instead, which invariance makes harmless.

MAP estimation and regularization

Treating $\theta$ as a random variable with a prior $p(\theta)$ and maximizing the posterior from Bayes’ rule gives maximum a posteriori (MAP) estimation (Murphy, 2022, §4.5):

$$\hat{\theta}_{\text{MAP}} = \arg\max_{\theta} \big[ \log p(\mathcal{D} \mid \theta) + \log p(\theta) \big].$$

It is maximum likelihood plus a log-prior, and it equals the maximum likelihood estimate when the prior is uniform (Murphy, 2022, §4.2.2). The log-prior acts as a penalty, which makes MAP estimation the probabilistic counterpart of regularization. A zero-mean Gaussian prior on weights gives the squared $L^2$ penalty, weight decay (Murphy, 2022, §4.5.3). A Laplace prior gives the $L^1$ penalty, which favours sparse solutions (Murphy, 2022, §11.4.1). A beta prior on a Bernoulli parameter, or a Dirichlet prior on histogram probabilities, adds pseudo-counts to the observed counts; with one pseudo-count per bin this is add-one smoothing, the standard remedy for bins that happen to be empty (Murphy, 2022, §4.5.1). Unlike the maximum likelihood estimate, a MAP estimate depends on how the parameter is expressed, for the reason Fisher gave. Keeping the whole posterior instead of its peak is Bayesian inference.

Examples

Estimating a camera’s noise parameters

Raw sensor noise follows, to a good approximation, the Poisson–Gaussian model of Image Noise: a pixel with expected value $\mu$ has variance $a\mu + b$, where $a$ is the gain from photoelectrons to output units (digital numbers, DN) and $b$ the variance of the read noise. The code renders a raw image of a step-wedge chart with eight flat patches whose expected values range from 5 to 200 DN, with $a = 0.5$ and $b = 4$. It estimates $a$, $b$ and the eight patch means by maximizing the likelihood of the Gaussian approximation $y \sim \mathcal{N}(\mu_j, a\mu_j + b)$, which has no closed-form maximum, and compares the result with a common simpler method: a straight line fitted through the (mean, variance) pairs of the patches. Standard errors come from the inverse Hessian of the NLL, and 200 independent images measure the actual spread of both estimators.

import numpy as np
from scipy import optimize

def render(rng, gain=0.5, read_sd=2.0):
    """Raw image of a step-wedge chart: 8 flat patches of 32 x 32 pixels with
    expected values 5 to 200 DN. Each pixel is gain * Poisson(mu / gain), the
    photoelectrons converted to DN, plus Gaussian read noise, so its variance
    is a * mu + b with a = gain = 0.5 and b = read_sd**2 = 4."""
    mu = np.geomspace(5, 200, 8)
    lam = np.repeat(mu / gain, 32 * 32).reshape(8, 32 * 32)
    return gain * rng.poisson(lam) + rng.normal(0, read_sd, lam.shape)

def nll(theta, y):
    """Negative log-likelihood of the Gaussian approximation
    y ~ N(mu_j, a * mu_j + b), theta = (a, b, mu_1, ..., mu_8)."""
    a, b, mu = theta[0], theta[1], theta[2:, None]
    var = a * mu + b
    return 0.5 * np.sum(np.log(2 * np.pi * var) + (y - mu) ** 2 / var)

def fit(y):
    # Moment-based fit: a straight line through (patch mean, patch variance).
    m, v = y.mean(axis=1), y.var(axis=1, ddof=1)
    a_ls, b_ls = np.polyfit(m, v, 1)
    # Maximum likelihood has no closed form here: minimize the NLL
    # numerically, starting from the line fit, keeping a and b positive.
    x0 = np.r_[max(a_ls, 0.01), max(b_ls, 0.01), m]
    bounds = [(1e-6, None), (1e-6, None)] + [(None, None)] * 8
    res = optimize.minimize(nll, x0, args=(y,), method="L-BFGS-B", bounds=bounds)
    return np.array([a_ls, b_ls]), res.x

def hessian(f, x, h=1e-4):
    """Central-difference Hessian of f at x."""
    n = len(x)
    H = np.empty((n, n))
    E = np.eye(n) * h
    for i in range(n):
        for j in range(n):
            H[i, j] = (f(x + E[i] + E[j]) - f(x + E[i] - E[j])
                       - f(x - E[i] + E[j]) + f(x - E[i] - E[j])) / (4 * h * h)
    return H

y = render(np.random.default_rng(0))
ls, theta = fit(y)
# Standard errors from the observed Fisher information (Hessian of the NLL).
se = np.sqrt(np.diag(np.linalg.inv(hessian(lambda t: nll(t, y), theta))))
print("truth       a = 0.500          b = 4.000")
print(f"line fit    a = {ls[0]:.3f}          b = {ls[1]:.3f}")
print(f"ML          a = {theta[0]:.3f} +- {se[0]:.3f}  b = {theta[1]:.3f} +- {se[1]:.3f}")

# Repeat over 200 independent images to measure each estimator's spread.
est = []
for seed in range(1, 201):
    ls_s, theta_s = fit(render(np.random.default_rng(seed)))
    est.append(np.r_[ls_s, theta_s[:2]])
est = np.array(est)
for name, c in [("line fit", 0), ("ML", 2)]:
    m, sd = est[:, c:c + 2].mean(axis=0), est[:, c:c + 2].std(axis=0)
    print(f"{name:9s} 200 images: a mean {m[0]:.3f} sd {sd[0]:.3f}, "
          f"b mean {m[1]:.3f} sd {sd[1]:.3f}")

Output:

truth       a = 0.500          b = 4.000
line fit    a = 0.522          b = 3.249
ML          a = 0.508 +- 0.013  b = 3.901 +- 0.244
line fit  200 images: a mean 0.501 sd 0.020, b mean 3.935 sd 0.716
ML        200 images: a mean 0.499 sd 0.012, b mean 3.984 sd 0.236

Both estimators are centred on the truth over the 200 images, but the maximum likelihood estimates scatter much less: by a factor of about 1.7 for the gain and 3 for the read-noise variance. The likely reason is that the line fit treats the eight patch variances as equally reliable, although the variance measured on a bright patch is itself far noisier than on a dark one, and the read noise $b$, the intercept, is determined by the dark patches that the noisy bright ones outweigh. The likelihood weights every pixel by its own variance automatically. The gap is therefore a cost of the unweighted fit rather than of fitting a line: by our reasoning, the sample variance of Gaussian data has a spread proportional to the variance itself, so weighting each patch by the inverse square of its expected variance would attack the same weakness. The standard errors computed from the Hessian of one image’s NLL, 0.013 and 0.244, predict the spread actually observed over 200 images, 0.012 and 0.236: with 8,192 pixels, the large-sample Gaussian approximation already holds. The likelihood used here is the Gaussian approximation of the Poisson–Gaussian model, which is accurate when every patch collects more than a few photoelectrons per pixel (here at least 10 on average); for very dark images the exact Poisson–Gaussian likelihood would be needed. Once $a$ and $b$ are known, the weights of a weighted least-squares fit to such a camera’s images follow from them.

Image alignment: SSD, SAD and the noise model

Estimating the translation between a known template and a noisy observed image by maximum likelihood minimizes the sum of squared differences (SSD) when the noise is Gaussian and the sum of absolute differences (SAD) when it is Laplace. The code renders a template of 40 smooth blobs that can be evaluated exactly at any sub-pixel shift, observes it shifted by $(1.3, -0.8)$ pixels with three kinds of noise of the same standard deviation, 4 grey levels (Gaussian; Laplace; and Gaussian with 10% of the pixels replaced by random values, as an occluder or dead pixels would), and estimates the shift by minimizing each cost over 20 random scenes.

import numpy as np
from scipy import optimize

yy, xx = np.mgrid[:64, :64].astype(float)

def make_scene(rng, k=40):
    """A known template: 40 Gaussian blobs (sigma 3 px) with random positions
    and brightness, rendered exactly at any sub-pixel translation t."""
    pos = rng.uniform(-8, 72, (k, 2))
    amp = rng.uniform(-60, 60, k)
    def render(t):
        d2 = (yy[..., None] - pos[:, 0] - t[0]) ** 2 + (xx[..., None] - pos[:, 1] - t[1]) ** 2
        return 120 + np.exp(-d2 / (2 * 3.0 ** 2)) @ amp
    return render

def observe(rng, render, true, noise):
    """The template translated by `true`, plus noise of standard deviation 4."""
    img = render(true)
    if noise == "Gaussian":
        return img + rng.normal(0, 4, img.shape)
    if noise == "Laplace":                         # scale 4/sqrt(2): sd 4
        return img + rng.laplace(0, 4 / np.sqrt(2), img.shape)
    out = img + rng.normal(0, 4, img.shape)        # Gaussian, then 10% of
    bad = rng.random(img.shape) < 0.10             # pixels replaced by
    out[bad] = rng.uniform(0, 255, bad.sum())      # random values
    return out

def estimate(obs, render, cost):
    """Translation minimizing SSD (Gaussian NLL) or SAD (Laplace NLL)."""
    def f(t):
        r = obs - render(t)
        return np.sum(r ** 2) if cost == "SSD" else np.sum(np.abs(r))
    return optimize.minimize(f, np.zeros(2), method="Powell",
                             options={"xtol": 1e-4, "ftol": 1e-10}).x

true = np.array([1.3, -0.8])
print(f"{'noise':>18}  {'SSD error (px)':>23}  {'SAD error (px)':>23}")
for noise in ["Gaussian", "Laplace", "Gaussian+outliers"]:
    err = {"SSD": [], "SAD": []}
    for seed in range(20):
        rng = np.random.default_rng(seed)
        render = make_scene(rng)
        obs = observe(rng, render, true, noise)
        for cost in err:
            err[cost].append(np.linalg.norm(estimate(obs, render, cost) - true))
    print(f"{noise:>18}  " + "  ".join(
        f"RMS {np.sqrt(np.mean(np.square(e))):.4f} max {np.max(e):.4f}"
        for e in err.values()))

Output:

             noise           SSD error (px)           SAD error (px)
          Gaussian  RMS 0.0269 max 0.0444  RMS 0.0378 max 0.0670
           Laplace  RMS 0.0257 max 0.0444  RMS 0.0202 max 0.0381
 Gaussian+outliers  RMS 0.2042 max 0.3589  RMS 0.0448 max 0.0836

Each cost wins under the noise it is the likelihood of. With Gaussian noise, SSD, the maximum likelihood estimator, gives the smaller error, and SAD pays for ignoring the magnitude of residuals: its RMS error is 1.4 times larger over these 20 scenes. The median shows the same loss of efficiency against the mean: for Gaussian data its variance is $\pi/2$ times larger (Murphy, 2022, §5.3.1.1). By our reasoning the same factor applies here, because the large-sample variances of the squared-loss and absolute-loss fits depend on the noise only through its variance and its density at zero residual respectively, exactly as those of the mean and the median do. It predicts an RMS ratio near $\sqrt{\pi/2} \approx 1.25$; the 1.4 measured over only 20 scenes is a noisy estimate of it. With Laplace noise of the same standard deviation, the roles swap and SAD is more accurate. With 10% of the pixels replaced, the Gaussian model is badly wrong: every corrupted pixel pulls the SSD estimate in proportion to its residual, and its error grows more than sevenfold, while SAD, whose heavier-tailed noise model gives every pixel the same bounded pull however large its residual, loses little accuracy. Choosing a cost function is choosing a noise model, and when the noise has outliers, a heavier-tailed likelihood or a robust loss (Nonlinear Least Squares) is the better assumption.

A colour histogram from few pixels: maximum likelihood versus MAP

The code renders a textured object region, quantizes its colours into 512 bins, keeps 20,000 pixels aside as test data and estimates the colour distribution from $n$ training pixels. It compares the maximum likelihood histogram with MAP estimates that add $c$ pseudo-counts to every bin (a Dirichlet prior), scoring each by the average log-likelihood it gives the test pixels.

import numpy as np

rng = np.random.default_rng(0)

# A textured object region: 200 x 200 pixels whose RGB colour mixes a red
# and an orange material with per-pixel shading and sensor noise.
h = w = 200
mat = rng.random((h, w)) < 0.7
base = np.where(mat[..., None], [190.0, 60.0, 50.0], [220.0, 140.0, 40.0])
img = base * (1 + 0.12 * rng.standard_normal((h, w, 1))) + rng.normal(0, 10, (h, w, 3))
img = np.clip(img, 0, 255)
q = (img // 32).astype(int)                       # 8 levels per channel
colour = ((q[..., 0] * 8 + q[..., 1]) * 8 + q[..., 2]).ravel()   # 512 bins

test = colour[:20000]                             # held-out pixels
pool = colour[20000:]
p_true = np.bincount(colour, minlength=512) / colour.size
print("bins used by the whole region:", np.count_nonzero(p_true))

print(f"{'n':>6} {'unseen':>7} {'ML':>6} {'MAP c=1':>8} {'MAP c=0.05':>10}")
for n in [50, 200, 1000, 5000, 20000]:
    counts = np.bincount(rng.choice(pool, n, replace=False), minlength=512)
    p_ml = counts / n                              # ML: relative frequencies
    unseen = np.sum(counts[test] == 0)             # test pixels in empty bins
    ll = []
    for p in [p_ml] + [(counts + c) / (n + 512 * c) for c in (1, 0.05)]:
        with np.errstate(divide="ignore"):         # MAP: c pseudo-counts per bin
            ll.append(np.mean(np.log(p[test])))    # mean test log-likelihood
    print(f"{n:>6} {unseen:>7d} {ll[0]:>6.2f} {ll[1]:>8.2f} {ll[2]:>10.2f}")
print(f"histogram of the whole region (test pixels included): "
      f"{np.mean(np.log(p_true[test])):.2f}")

Output:

bins used by the whole region: 70
     n  unseen     ML  MAP c=1 MAP c=0.05
    50    2750   -inf    -5.01      -3.66
   200     831   -inf    -4.09      -3.19
  1000      93   -inf    -3.36      -3.03
  5000      40   -inf    -3.07      -3.00
 20000       6   -inf    -3.01      -2.99
histogram of the whole region (test pixels included): -2.98

The maximum likelihood histogram fits its training pixels as well as any histogram can, yet it assigns probability zero to every colour it did not see, so a single test pixel in an empty bin makes the test likelihood zero. Even with 20,000 training pixels, six test pixels fall in empty bins: the region uses 70 colour bins, some of them rarely. This is overfitting in its simplest form, the zero-count problem. Both MAP estimates keep every probability positive and approach the whole-region reference as $n$ grows, and the prior’s influence fades: at $n = 20{,}000$ the two pseudo-counts differ by 0.02. With few pixels the prior’s strength matters. Add-one smoothing spreads 512 pseudo-counts over 512 bins, ten times more than the 50 real pixels, and flattens the histogram towards uniform, while 0.05 pseudo-counts per bin keep the empty bins possible without drowning the data. The pseudo-count is a hyperparameter, chosen on held-out data like any regularization strength. The skin-colour classifier in Probability and Statistics adds one pseudo-count per bin for the same reason.

Limitations

• Small samples and overfitting. The guarantees above are asymptotic. With few observations, the estimate can be biased (the Gaussian variance), extreme (a coin that came up heads three times gets $\hat{\theta} = 1$; Bishop, 2006, §1.2.3) or degenerate (the empty histogram bins above). A flexible model fitted by maximum likelihood fits the noise in its training data, and its training likelihood rises while its likelihood on new data falls; Bishop traces the overfitting of such models to the bias of maximum likelihood (Bishop, 2006, §1.2.4). Priors, regularization and held-out evaluation are the remedies.
• Unbounded likelihoods. For some models the likelihood has no maximum at all. A Gaussian mixture model can place one component’s mean exactly on a single data point and shrink its variance towards zero; that point’s density then grows without bound while the other components keep the remaining points probable, so the log-likelihood goes to infinity (Bishop, 2006, §9.2.1). A single Gaussian fitted to several distinct points has no such singularity, because the collapsed density makes all the other points improbable. Similarly, the maximum likelihood covariance of $d$-dimensional feature vectors becomes singular in high dimensions (Murphy, 2022, §4.5.2); it always is when there are no more samples than dimensions, since the centred samples span at most $n - 1$ dimensions. Practical fitting code detects collapsing components and resets them (Bishop, 2006, §9.2.1), adds a small constant to variances, or uses a prior.
• Non-identifiability. When different parameter values produce the same distribution of the data, the likelihood has ridges or several equal peaks, and the estimate is not unique. The $K$ components of a mixture can be relabelled in $K!$ ways without changing the model (Bishop, 2006, §9.2.1). A 3D reconstruction from images is determined only up to a similarity transformation (for calibrated cameras), the gauge freedom of bundle adjustment, which must be fixed by convention before the parameters are unique and their Fisher information invertible. Asymptotic normality is stated for identifiable parameters (Murphy, 2022, §4.7.2).
• Misspecification. Maximum likelihood is only as good as the model. When the noise is not what the model assumes, as with outliers under a Gaussian model in the alignment example, the estimate can be badly wrong, and standard errors computed from the Hessian, which assume the model is right, need not describe its real precision. When the true distribution is outside the model family, the estimate still minimizes the Kullback–Leibler divergence from the empirical distribution to the model (Murphy, 2022, §4.2.2), so it settles on the member of the family closest to the data in that sense, which may not be the parameter the user had in mind.
• Computation. Without a closed form, the estimate is only as good as the optimizer that finds it: local maxima, poor starting points and unconstrained parameters (negative variances, probabilities outside $[0, 1]$) are common failure points.

Common Misconceptions

• “The likelihood is the probability that the parameters are right.” It is the probability of the data for given parameters, read as a function of the parameters. It does not integrate to one over them, and turning it into a probability over parameters requires a prior and Bayes’ rule.
• “Maximum likelihood estimates are unbiased.” Many are not; the Gaussian variance with denominator $n$ is the standard example. They are consistent and asymptotically efficient under regularity conditions, which is a different promise.
• “Least squares is always the maximum likelihood estimate.” Only under independent Gaussian noise of equal variance. Unequal variances call for weighted least squares, heavy-tailed noise for absolute or robust losses, and counts at low light for a Poisson likelihood.
• “A higher likelihood means a better model.” Only on data the model was not fitted to. On the training data, a more flexible model always reaches a higher likelihood, which is exactly how overfitting shows up.
• “MAP estimation is just maximum likelihood with an extra term, so it shares its properties.” A MAP estimate depends on how the parameters are expressed and on the prior’s strength; maximum likelihood depends on neither. With enough data the likelihood dominates and the two agree.

Where It Is Used

• Noise estimation and denoising: fitting Poisson–Gaussian or Gaussian image noise models from flat regions or calibration frames, and the weights of subsequent fits.
• Geometric estimation: least squares fits of lines, homographies and camera matrices; camera calibration and bundle adjustment minimize the reprojection error, the maximum likelihood objective under Gaussian image noise.
• Image alignment and tracking: SSD and SAD costs in template matching, optical flow and Lucas–Kanade tracking, and the measurement models of the Kalman filter.
• Classifiers and networks: the cross-entropy loss and other negative log-likelihood losses (Loss Function) that train image classifiers, detectors and segmentation networks; logistic regression on fixed features; weight decay as a MAP prior.
• Density models: colour histograms and Gaussian colour models for skin detection and segmentation; Gaussian mixture models fitted by expectation–maximization for background subtraction and clustering; generative models trained by maximizing the likelihood or a bound on it, such as the variational autoencoder and its evidence lower bound.

History

Gauss justified least squares in 1809 by asking which law of errors would make the arithmetic mean the most probable value, an argument in the spirit of the method (Gaussian Distribution). As Aldrich (1997) traces it, R. A. Fisher first presented the numerical procedure in 1912 as an undergraduate, under the name “absolute criterion” and with a justification from inverse probability. In 1921 he called the most likely value the “optimum” and separated likelihood from probability. The 1922 paper “On the Mathematical Foundations of Theoretical Statistics” introduced the name “maximum likelihood”. It defined likelihood as proportional to the probability that, for the given parameter values, the observations would be those observed (Fisher, 1922, p. 310). It also set out the criteria of consistency, efficiency and sufficiency against which estimators are still judged. Fisher believed the method leads automatically to sufficient statistics, but said he could not yet prove it with full rigour (Fisher, 1922, p. 323). Aldrich (1997) reports Pratt’s (1976) finding that Edgeworth had anticipated a good part of the 1922 version in 1908, which nobody noticed until a decade or so after Fisher had redone it. Fisher also used invariance under a change of parameter to argue against assigning uniform prior probabilities to parameters (Aldrich, 1997, §8), a point on which maximum likelihood and Bayesian estimation still differ.