M11-T1 — Maximum Likelihood Estimation

Learning objectives

After this lesson, you should be able to:

• Define the likelihood and the maximum likelihood estimate, and derive the MLE of a Gaussian mean and variance.
• Show that maximum likelihood with Gaussian noise is (weighted) least squares, and solve it by Gauss–Newton.
• Compute the Fisher information and the Cramér–Rao bound, and read the uncertainty ellipse of an estimate.
• Explain why geometry (the placement of landmarks) decides accuracy, and why the 1/n variance estimate is biased.

The intuition

A robot measures its distance to three landmarks and gets three noisy numbers. The true position is the one that would have produced those distances; the readings do not quite agree, so there is no point that fits all three. The maximum likelihood estimate is the point that makes the observed readings most probable: for each candidate position, compute the probability of the three readings (the likelihood), and take the best.

If the noise is Gaussian, the probability of a reading falls off with the square of its error, so maximising the likelihood is the same as minimising the sum of squared errors, each divided by its variance: the familiar least-squares fit.

The left panel repeats the experiment 2 000 times (400 shown) for each of two layouts. When the landmarks surround the robot, the estimates form a round cloud of standard deviation 8.5 and 8.1 cm in x and y: three readings of 10 cm each combine into something slightly better than one. When the landmarks are all on one side, the cloud is stretched: 16 cm in x, 6.7 cm in y. The reason is geometry: a range to a distant landmark only tells you the position along the line to it, and when the three lines are nearly parallel nothing constrains the direction across them. The ellipses are not fitted to the clouds. They are computed from the Jacobian alone, before any data, by the Cramér–Rao bound, and they match: 0.084 and 0.082 against 0.085 and 0.081.

The right panel is a warning that the best-known estimator has a flaw. The maximum likelihood estimate of a variance divides the sum of squares by n, and it is biased low: by −σ²/n. With 3 readings and a true variance of 4 it is wrong by −1.31 on average, and with 50 readings by −0.07. Dividing by n − 1 removes the bias (the measured bias of that version is 0.04, −0.01, 0.00 and 0.01).

Maximum likelihood turns estimation into optimisation: write the negative log-likelihood, minimise it. With Gaussian noise it is a weighted sum of squares, and the curvature at the minimum is the information about the state.

Definitions

Likelihood $L(\theta)$

$p(\mathbf{z}\mid\theta)$ read as a function of the parameters $\theta$ with the data $\mathbf{z}$ fixed. It is not a probability distribution over $\theta$.

Maximum likelihood estimate

$\hat\theta_\mathrm{ML}=\arg\max_\theta\,p(\mathbf{z}\mid\theta)=\arg\min_\theta\,[-\log p(\mathbf{z}\mid\theta)]$.

Residual $\mathbf{r}_i(\theta)$

$\mathbf{z}_i-h_i(\theta)$: the reading minus the reading the state would give.

Fisher information $\mathcal{I}(\theta)$

$-\mathbb{E}\bigl[\nabla^2_\theta\log p(\mathbf{z}\mid\theta)\bigr]$: the expected curvature of the log-likelihood; for Gaussian noise $J^\top R^{-1}J$.

Cramér–Rao bound

The covariance of any unbiased estimator is at least $\mathcal{I}(\theta)^{-1}$. The MLE reaches it as the number of readings grows (for Gaussian noise with linear models, exactly).

Derivation: Gaussian MLE, least squares, and the bound

The Gaussian mean and variance

For $n$ independent readings $z_i\sim\mathcal{N}(\mu,\sigma^2)$:

1. $$\ell(\mu,\sigma^2)=\log L=-\tfrac n2\log(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum_i(z_i-\mu)^2$$

   The log of the product of Gaussian densities: independent readings multiply, so their logs add.

2. $$\frac{\partial\ell}{\partial\mu}=\frac1{\sigma^2}\sum_i(z_i-\mu)=0;\Rightarrow;\hat\mu=\frac1n\sum_iz_i$$

   Setting the derivative to zero gives the sample mean.

3. $$\frac{\partial\ell}{\partial\sigma^2}=-\frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}\sum_i(z_i-\mu)^2=0;\Rightarrow;\hat\sigma^2=\frac1n\sum_i(z_i-\hat\mu)^2$$

   The MLE of the variance divides by n.

4. $$\mathbb{E}[\hat\sigma^2]=\frac{n-1}{n}\sigma^2$$

   Because $\hat\mu$ was fitted to the same data, the residuals are smaller than the true deviations. The expectation is $\sigma^2-\sigma^2/n$: a bias of $-\sigma^2/n$.

Least squares and Gauss–Newton

1. $$-\log p(\mathbf{z}\mid\theta)=\tfrac12\sum_i\mathbf{r}_i^\top R_i^{-1}\mathbf{r}_i+\text{const}$$

   With Gaussian noise $\mathbf{z}_i=h_i(\theta)+\mathbf{v}_i$, the log-density is a quadratic in the residual. Minimising it is the least-squares problem.

2. $$\mathbf{r}(\theta+\Delta)\approx\mathbf{r}(\theta)-J\Delta$$

   Linearise the residuals ($J=\partial h/\partial\theta$) around the current estimate: the same first-order step as the EKF.

3. $$\min_\Delta\tfrac12\lVert\mathbf{r}-J\Delta\rVert^2_{R^{-1}};\Rightarrow;J^\top R^{-1}J\,\Delta=J^\top R^{-1}\mathbf{r}$$

   A linear weighted least-squares problem: the normal equations. Solve, update, repeat. Near the solution the convergence is quadratic for small residuals.

The information and the bound

1. $$\mathcal{I}=-\mathbb{E}[\nabla^2\log L]=J^\top R^{-1}J$$

   The Hessian of the quadratic negative log-likelihood is $J^\top R^{-1}J$ (Gauss–Newton ignores the residual-curvature term, which has expectation zero at the truth).

2. $$\hat\theta-\theta\approx\mathcal{I}^{-1}J^\top R^{-1}\mathbf{v};\Rightarrow;\mathrm{Cov}(\hat\theta)=\mathcal{I}^{-1}J^\top R^{-1}R R^{-1}J\,\mathcal{I}^{-1}=\mathcal{I}^{-1}$$

   To first order the estimation error is a linear function of the noise, so it is Gaussian with covariance equal to the inverse information. This is the Cramér–Rao bound, attained.

3. $$\mathcal{I}_\mathrm{total}=\sum_i J_i^\top R_i^{-1}J_i$$

   Information adds across independent readings, the same additivity as the precisions of the fusion lesson of Module 4.

Geometry

For a range to a landmark in direction $\mathbf{u}_i$ from the robot, $J_i=\mathbf{u}_i^\top$ and $\mathcal{I}=\sigma^{-2}\sum_i\mathbf{u}_i\mathbf{u}_i^\top$. Each landmark adds information only along its own direction. Three landmarks in widely different directions give a well-conditioned $\mathcal{I}$; three almost in the same direction give a nearly singular one, with a small eigenvalue along the unconstrained direction: the long axis of the orange ellipse. The ratio of the standard deviations is the same geometric dilution of precision that GPS engineers quote.

MLE is optimisation; its covariance comes for free as the inverse of the Gauss–Newton matrix JᵀR⁻¹J. The weights R⁻¹ matter: they are what makes the estimate the maximum likelihood and the covariance the right one.

A worked example

Two programs.

1. Gaussian MLE. Readings of a quantity with true variance 4.0, 20 000 trials per n. The 1/n variance has an average error of −1.309 for n = 3, −0.809 for n = 5, −0.400 for n = 10 and −0.072 for n = 50: −σ²/n = −1.33, −0.8, −0.4 and −0.08. The 1/(n−1) version has errors of 0.036, −0.011, −0.001 and 0.008, all within the noise of the experiment.
2. Position from three ranges. A robot at (4, 3) with 10 cm of noise on each range, estimated by Gauss–Newton from a start 0.7 m away; 2 000 trials. Landmarks spread around: standard deviations 0.085 and 0.081 against the bound 0.084 and 0.082, correlation −0.14. Landmarks on one side: 0.160 and 0.067 against 0.160 and 0.065, correlation 0.30. The bound predicts the scatter in both cases.
3. Interpretation. The one-sided layout makes the x error twice as large as the y error because the landmarks, all in the +y direction, constrain y directly and x only through the small angular differences. The root-mean-square position error is 0.117 m for the spread layout and 0.173 m for the one-sided one, although the sensor is the same.

import math
from robotprob import Rng, mean, variance, solve, matmul, transpose, matvec, inverse, sample_mvn

# 1. the maximum likelihood estimates of a Gaussian: the sample mean, and a variance with 1/n
print('n    bias of the 1/n variance   bias of the 1/(n-1) variance   (true variance 4.0, 20000 trials)')
for n in (3, 5, 10, 50):
    rng = Rng(n)
    a, b = [], []
    for _ in range(20000):
        xs = [rng.gauss(10.0, 2.0) for _ in range(n)]
        a.append(variance(xs, ddof=0)); b.append(variance(xs, ddof=1))
    print(str(n).rjust(3), str(round(mean(a) - 4.0, 3)).rjust(16), str(round(mean(b) - 4.0, 3)).rjust(28))

# 2. maximum likelihood position from ranges to landmarks, found by Gauss-Newton, against the Cramer-Rao bound
def locate(landmarks, ranges, sigma, x0):
    x = list(x0)
    for _ in range(20):
        J, r = [], []
        for (lx, ly), z in zip(landmarks, ranges):
            d = math.hypot(x[0] - lx, x[1] - ly)
            J.append([(x[0] - lx) / d / sigma, (x[1] - ly) / d / sigma])     # the rows are scaled by 1/sigma: weighted least squares
            r.append((z - d) / sigma)
        JT = transpose(J)
        step = solve(matmul(JT, J), matvec(JT, r))
        x = [x[0] + step[0], x[1] + step[1]]
        if abs(step[0]) + abs(step[1]) < 1e-9: break
    return x, matmul(JT, J)

true_pos, sigma = (4.0, 3.0), 0.1
setups = {'landmarks spread around': [(0.0, 0.0), (10.0, 0.0), (5.0, 9.0)],
          'landmarks all on one side': [(0.0, 9.0), (3.0, 9.5), (6.0, 9.0)]}
print('MLE of a position from 3 range readings (noise 0.1 m), 2000 trials each')
print('geometry                     std x   std y   Cramer-Rao std x, y    rms error   corr(x, y)')
for name, lm in setups.items():
    rng = Rng(21)
    est = []
    for _ in range(2000):
        z = [math.hypot(true_pos[0] - lx, true_pos[1] - ly) + rng.gauss(0.0, sigma) for lx, ly in lm]
        x, info = locate(lm, z, sigma, (true_pos[0] + 0.5, true_pos[1] - 0.5))
        est.append(x)
    ex, ey = [e[0] for e in est], [e[1] for e in est]
    _, info = locate(lm, [math.hypot(true_pos[0] - lx, true_pos[1] - ly) for lx, ly in lm], sigma, true_pos)
    crb = inverse(info)
    sx, sy = math.sqrt(variance(ex)), math.sqrt(variance(ey))
    cxy = sum((a - mean(ex)) * (b - mean(ey)) for a, b in zip(ex, ey)) / (len(ex) - 1) / (sx * sy)
    rms = math.sqrt(mean([(a - true_pos[0]) ** 2 + (b - true_pos[1]) ** 2 for a, b in est]))
    print(name.ljust(26), str(round(sx, 3)).rjust(6), str(round(sy, 3)).rjust(7), ' ', (round(math.sqrt(crb[0][0]), 3), round(math.sqrt(crb[1][1]), 3)), str(round(rms, 3)).rjust(12), str(round(cxy, 2)).rjust(9))

n    bias of the 1/n variance   bias of the 1/(n-1) variance   (true variance 4.0, 20000 trials)
  3           -1.309                        0.036
  5           -0.809                       -0.011
 10             -0.4                       -0.001
 50           -0.072                        0.008
MLE of a position from 3 range readings (noise 0.1 m), 2000 trials each
geometry                     std x   std y   Cramer-Rao std x, y    rms error   corr(x, y)
landmarks spread around     0.085   0.081   (0.084, 0.082)        0.117     -0.14
landmarks all on one side    0.16   0.067   (0.16, 0.065)        0.173       0.3

In robotics

Everywhere a state is fitted to data

Maximum likelihood is the principle behind calibration (camera intrinsics by reprojection error, IMU biases), trilateration with UWB or range beacons, bundle adjustment, scan matching (ICP minimises squared distances, a least-squares fit) and the front end of every SLAM system. The SLAM course builds on it in its optimisation lessons.

Reading the information matrix

• The eigenvectors of $\mathcal{I}^{-1}$ are the axes of the uncertainty ellipse; the eigenvalues are the variances along them.
• A small eigenvalue of $\mathcal{I}$ (a large one of the covariance) marks a direction that the data barely constrain: a weakly observable direction (Module 7). The solver converges slowly and noise is amplified in that direction.
• A zero eigenvalue means a gauge freedom or an unobservable state, as for the sum of two constants in Module 7. It must be fixed by a prior or an anchor.

Design before data

Since the bound depends only on the Jacobian and the noise, it can be computed before building the system: how many beacons, where to put them, which sensors to add. It gives the best that any estimator can do, which is also the benchmark to compare an implementation against: an estimator whose scatter is far above the bound has a bug or a model error.

Practical points about Gauss–Newton

• Start from a reasonable estimate (an odometry guess). A bad start can converge to a local minimum, or diverge.
• Add damping (Levenberg–Marquardt) when the step is too long: $(J^\top R^{-1}J+\lambda I)\Delta=J^\top R^{-1}\mathbf{r}$.
• Use the covariance in the weights: readings with different noise must be weighted by their own R, not equally.
• Check the Jacobian against finite differences, as for the EKF.

Common misconceptions

Key takeaways

1. Maximum likelihood chooses the state that makes the data most probable: minimise the negative log-likelihood.
2. For Gaussian noise this is weighted least squares: Σ rᵀR⁻¹r. Gauss–Newton solves it by repeated linearisation: JᵀR⁻¹J Δ = JᵀR⁻¹ r.
3. The Fisher information JᵀR⁻¹J is the curvature of the objective, and its inverse is the covariance of the estimate: the Cramér–Rao bound, which the experiments met (0.085 against 0.084).
4. Geometry decides accuracy: landmarks in similar directions leave one direction poorly known (16 cm against 6.7 cm). The bound can be computed before any data.
5. The MLE of a variance has bias −σ²/n; divide by n − 1 for small samples.

References and further reading

Each source is marked free or paid. Free to read does not mean free to reuse: the course resources page lists every source with its licence.

• S. Thrun, W. Burgard and D. Fox, Probabilistic Robotics, MIT Press, 2005, section 3.2 for the relation with the information form, and chapter 11 on SLAM as least squares. Paid book.
• T. D. Barfoot, State Estimation for Robotics, 2nd edition, 2024 , chapter 2 (primer on probability: the MLE, Fisher information and the Cramér–Rao bound) and chapter 4. Free to read as the author’s PDF; the book is copyrighted.
• Y. Bar-Shalom, X. R. Li and T. Kirubarajan, Estimation with Applications to Tracking and Navigation, Wiley, 2001, chapter 2 (the Cramér–Rao lower bound) and chapter 3. Paid book.
• C. Stachniss, “Graph-based SLAM using Pose Graphs” (video lecture) , which uses this least-squares view. Free to watch.

Knowledge check

Three questions. Write an answer, then reveal the explanation. Your answers stay in your browser.

What’s next

Maximum likelihood ignores what we knew before the readings. The next lesson adds a prior, which turns the estimate into a MAP estimate, and shows that the Kalman filter is exactly this kind of problem solved recursively.

Related content

Nothing below is required to finish this lesson.

Other courses on this site treat neighbouring topics from their own angle:

• SLAM M4-T2: Pose Graphs
• SLAM M4-L1: Visualize a Pose Graph
• UAV M4-T12: Pose graph optimization
• SLAM: Robust Optimization