M4-T3 — The Multivariate Gaussian and Covariance Ellipses

Learning objectives

After this lesson, you should be able to:

• Write the density of a multivariate Gaussian and explain each factor, including the role of the determinant.
• Turn a 2 × 2 covariance matrix into an ellipse using its eigenvalues and eigenvectors, and draw 1, 2 and 3 sigma contours.
• Compute the Mahalanobis distance, and explain why it differs from the Euclidean distance.
• Derive the probability inside a k-sigma ellipse, and choose a gating threshold from the chi-square distribution.

The intuition

The last lesson ended with a covariance matrix [[13, 12], [12, 13]] cm² for a diagonal drive. What does a Gaussian with this covariance look like? Every contour line of its density, the set of points that are equally plausible, is an ellipse. The ellipse is centred on the mean, its long axis points in the direction of greatest uncertainty, and its short axis in the direction of least. Here the long axis lies along the 45° diagonal and is five times longer than the short one.

Two things in this figure deserve a careful look.

• Distance depends on direction. The point (4, 0) is 4.0 cm from the centre, and the point (3, 3) is 4.2 cm away. But (4, 0) lies across the thin direction of the cloud, where errors of this size are rare, and is outside the 2σ ellipse. The point (3, 3) lies along the cloud and is inside the 1σ ellipse. A fair measure of “how far is this point from the mean?” must count distance in units of the local spread, in the direction of the point. That is the Mahalanobis distance: 2.88σ for the first point and 0.85σ for the second.
• The “1σ ellipse” is not a 68% region. In one dimension ±1σ holds 68%. In two dimensions the 1σ ellipse holds only 39%, the 2σ ellipse 86% and the 3σ ellipse 98.9%. To contain 95% in two dimensions you must draw the ellipse at 2.45σ. In three dimensions, 2.80σ. Spreading probability over more directions leaves less of it near the centre.

A covariance ellipse is the set of points at the same Mahalanobis distance from the mean. Its axes point along the eigenvectors of the covariance matrix and their lengths are the square roots of its eigenvalues, times the size k of the ellipse.

Definitions

Multivariate Gaussian $\mathcal{N}(\boldsymbol\mu, \Sigma)$

A distribution over a vector $\mathbf{x}$ of $n$ quantities with mean vector $\boldsymbol\mu$ and covariance matrix $\Sigma$ (symmetric, positive definite).

Determinant $|\Sigma|$

A single number that measures the volume of the cloud: the area of the 2D ellipse is proportional to $\sqrt{|\Sigma|}$. For a diagonal $\Sigma$ it is the product of the variances.

Precision (information) matrix $\Sigma^{-1}$

The inverse of the covariance. A large entry means a tightly known direction.

Mahalanobis distance $d$

The distance from the mean in units of the standard deviation in the direction of the point: $d^2 = (\mathbf{x}-\boldsymbol\mu)^{\top}\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu)$.

Eigenvalue $\lambda_i$ and eigenvector $\mathbf{v}_i$

The directions $\mathbf{v}_i$ along which $\Sigma$ only stretches, and the stretch factors $\lambda_i$ ($\Sigma \mathbf{v}_i = \lambda_i \mathbf{v}_i$). They are the variances along the principal axes of the cloud.

Derivation: the ellipse and the probability inside it

Why the contours are ellipses with those axes

1. $$\mathbf{y} = R^{\top}(\mathbf{x} - \boldsymbol\mu)$$

   Rotate the coordinates to the principal axes of the cloud, centred on the mean. The rotation $R^{\top}$ does not change distances.

2. $$\Sigma^{-1} = R\,\Lambda^{-1}R^{\top} \;\;\Rightarrow\;\; d^2 = \mathbf{y}^{\top}\Lambda^{-1}\mathbf{y} = \sum_{i} \dfrac{y_i^2}{\lambda_i}$$

   Substitute the decomposition. In the principal coordinates the quadratic form has no cross terms: each coordinate is divided by its own variance.

3. $$d^2 = k^2 \;\Longleftrightarrow\; \sum_i \dfrac{y_i^2}{k^2\lambda_i} = 1$$

   This is the equation of an ellipse (ellipsoid) with semi-axes $k\sqrt{\lambda_i}$ along the axes $y_i$, that is, along the eigenvectors $\mathbf{v}_i$.

4. $$\mathrm{Cov}(\mathbf{y}) = R^{\top}\Sigma R = \Lambda$$

   In the principal coordinates the covariance is diagonal: the components $y_i$ are uncorrelated and, being jointly Gaussian, independent, with variances $\lambda_i$.

The normalizing constant

1. $$\int \exp\!\bigl(-\tfrac12 \mathbf{y}^{\top}\Lambda^{-1}\mathbf{y}\bigr)\,d\mathbf{y} = \prod_{i}\int e^{-y_i^2/(2\lambda_i)}\,dy_i = \prod_i \sqrt{2\pi\lambda_i} = (2\pi)^{n/2}\sqrt{|\Lambda|}$$

   The integral factors into one-dimensional Gaussian integrals (the last module), one per principal axis.

2. $$|\Lambda| = |R^{\top}\Sigma R| = |\Sigma|$$

   A rotation does not change the determinant. Dividing by this constant gives the density above.

The probability inside the ellipse, in two dimensions

Standardize each principal coordinate, $z_i = y_i/\sqrt{\lambda_i}$: then $d^2 = z_1^2 + z_2^2$, the squared distance from the origin of a point with two independent standard normal coordinates. Its density is rotation-symmetric, so switch to polar coordinates, as in the derivation of the Gaussian constant:

1. $$P(d \le k) = \int_0^{2\pi}\!\int_0^{k} \dfrac{1}{2\pi}\, e^{-r^2/2}\, r\,dr\,d\theta$$

   The standard bivariate normal density is $\tfrac{1}{2\pi}e^{-(z_1^2+z_2^2)/2}$, and the area element in polar coordinates is $r\,dr\,d\theta$.

2. $$= \int_0^{k} r\,e^{-r^2/2}\,dr = 1 - e^{-k^2/2}$$

   The angle integral cancels the $1/(2\pi)$; the radial integral has the antiderivative $-e^{-r^2/2}$.

3. $$k = 1: 0.393 \qquad k = 2: 0.865 \qquad k = 3: 0.989 \qquad P = 0.95 \;\Rightarrow\; k = \sqrt{-2\ln 0.05} = 2.448$$

   The numbers behind the figure. Solving for the 95% ellipse inverts the formula.

In general $d^2$ is a sum of $n$ squared independent standard normals, which has the chi-square distribution with $n$ degrees of freedom. The threshold $d^2 \le \chi^2_{n,\,p}$ that contains a probability $p$ is available as chi2_ppf(p, n): 3.84, 5.99 and 7.81 for $n = 1, 2, 3$ at $p = 0.95$, which are the squares of 1.96, 2.45 and 2.80.

Use the chi-square quantile, not a fixed number of sigmas, to gate in more than one dimension: the 95% gate is 1.96σ in one dimension, 2.45σ in two and 2.80σ in three.

A worked example

Take the covariance matrix of the previous lesson, $\Sigma = \begin{pmatrix} 13 & 12 \\ 12 & 13 \end{pmatrix}$ cm², with mean at the origin.

1. Eigenvalues and axes. $\det(\Sigma - \lambda I) = (13 - \lambda)^2 - 144 = 0$ gives $\lambda = 13 \pm 12 = 25,\; 1$. The eigenvector of 25 is $(1, 1)/\sqrt{2}$, at 45°; that of 1 is $(-1, 1)/\sqrt{2}$. The 1σ ellipse has semi-axes $\sqrt{25} = 5$ cm and $\sqrt{1} = 1$ cm. The 2σ ellipse is 10 by 2 cm.
2. Density at the mean. $|\Sigma| = 169 - 144 = 25$, so $f(\mathbf{0}) = 1/(2\pi\sqrt{25}) = 0.0318$ per cm².
3. Two points, two distances. The inverse is $\Sigma^{-1} = \tfrac{1}{25}\begin{pmatrix} 13 & -12 \\ -12 & 13 \end{pmatrix}$. For $(4, 0)$: $d^2 = 13 \times 16/25 = 8.32$, $d = 2.88$σ. For $(3, 3)$: $d^2 = (117 - 216 + 117)/25 = 0.72$, $d = 0.85$σ. The Euclidean distances, 4.00 and 4.24 cm, would have ranked them the other way round.
4. Gate. A 95% gate in two dimensions is $d^2 \le 5.99$. The point (4, 0) fails it (8.32) although it is close in the ordinary sense; (3, 3) passes comfortably.
5. Check the shares by simulation. 20 000 samples: 39.6%, 86.6% and 98.9% lie inside the 1σ, 2σ and 3σ ellipses, against 39.35%, 86.47% and 98.89% from $1 - e^{-k^2/2}$.
6. Other dimensions. The 95% radius is 1.96σ in one dimension, 2.448σ in two, 2.795σ in three.

import math
from robotprob import Rng, sample_mvn, mahalanobis2, eig_sym2, covariance_ellipse, chi2_ppf, mvn_pdf

S = [[13.0, 12.0], [12.0, 13.0]]
l1, l2, ang = eig_sym2(S)
print('eigenvalues', round(l1, 3), round(l2, 3), ' long axis at', round(math.degrees(ang), 1), 'degrees')
a, b, _ = covariance_ellipse(S, 1.0)
print('1-sigma semi-axes =', round(a, 2), round(b, 2), 'cm')

mu = [0.0, 0.0]
for p in ([4.0, 0.0], [3.0, 3.0]):
    d2 = mahalanobis2(p, mu, S)
    print(p, 'Euclidean', round(math.hypot(*p), 2), 'cm   Mahalanobis', round(math.sqrt(d2), 2), 'sigma   d^2 =', round(d2, 2))
print('density at the mean =', round(mvn_pdf(mu, mu, S), 4), ' (1 / (2 pi sqrt(det)) =', round(1 / (2 * math.pi * 5), 4), ')')

rng = Rng(12)
samples = [sample_mvn(rng, mu, S) for _ in range(20000)]
d2s = [mahalanobis2(p, mu, S) for p in samples]
for k in (1, 2, 3):
    inside = sum(1 for d in d2s if d <= k * k) / len(d2s)
    print(k, 'sigma ellipse holds', round(inside, 4), ' formula 1 - exp(-k^2/2) =', round(1 - math.exp(-k * k / 2), 4))
print('95% gate: d^2 <=', round(chi2_ppf(0.95, 2), 3), '=', round(math.sqrt(chi2_ppf(0.95, 2)), 3), 'sigma;  in 1D it is 1.96 sigma, in 3D', round(math.sqrt(chi2_ppf(0.95, 3)), 3))

eigenvalues 25.0 1.0  long axis at 45.0 degrees
1-sigma semi-axes = 5.0 1.0 cm
[4.0, 0.0] Euclidean 4.0 cm   Mahalanobis 2.88 sigma   d^2 = 8.32
[3.0, 3.0] Euclidean 4.24 cm   Mahalanobis 0.85 sigma   d^2 = 0.72
density at the mean = 0.0318  (1 / (2 pi sqrt(det)) = 0.0318 )
1 sigma ellipse holds 0.3963  formula 1 - exp(-k^2/2) = 0.3935
2 sigma ellipse holds 0.8661  formula 1 - exp(-k^2/2) = 0.8647
3 sigma ellipse holds 0.9891  formula 1 - exp(-k^2/2) = 0.9889
95% gate: d^2 <= 5.991 = 2.448 sigma;  in 1D it is 1.96 sigma, in 3D 2.795

In robotics

Drawing and reading uncertainty

The ellipse that RViz draws around a robot's pose, the “covariance ellipse” in SLAM visualizations and every plot of a Kalman filter's position estimate is one of the contours above, usually the 1σ or 2σ one. Read the axes and the tilt: a long thin ellipse means one direction is poorly known, a round ellipse means the filter has no preferred direction. Always check which k is drawn, and remember that a 1σ ellipse holds 39% of the probability, not 68%. The SLAM course draws these in Understand Pose Uncertainty.

Gating measurements

Before a filter uses a measurement it asks how plausible it is: is the Mahalanobis distance between the measurement and its prediction below a threshold? The threshold comes from the chi-square distribution with as many degrees of freedom as the measurement has components, so a 95% gate is d² ≤ 3.84 for one number, 5.99 for a range and a bearing, 7.81 for three numbers. Using 1.96σ in every dimension would reject far too many good measurements as the dimension grows. M8-T2: EKF Localization with Landmarks applies this to landmark observations, and the UAV course to innovation monitoring (Normalized Estimation Error Squared).

Data association

When a robot sees an object that might be one of several known landmarks, the one with the smallest Mahalanobis distance is the best match, because the distance measures the observation in the units of its own uncertainty. A match by Euclidean distance would prefer a landmark that is close in the direction in which the observation is least reliable (data association, Data Association).

Comparing two estimators

Two filters that report the same variance of x can have very different covariance matrices. Compare the determinants (the volumes) for an overall measure, the eigenvalues for the worst direction, and the ellipses overlaid on a plot for the geometry. The largest eigenvalue is the worst-case standard deviation, the number that determines whether a robot fits through a door.

Plausibility checks

For a ground-truth experiment, the fraction of true errors that fall inside the filter's own k-sigma ellipse should match the figure above (39%, 86%, 99%). A filter whose 2σ ellipse contains the truth only 50% of the time is overconfident. This is the idea behind the NEES test of module 7.

Common misconceptions

Key takeaways

1. A multivariate Gaussian N(μ, Σ) has contours that are ellipses centred on μ; its density falls with the squared Mahalanobis distance d² = (x − μ)ᵀΣ⁻¹(x − μ).
2. The ellipse axes lie along the eigenvectors of Σ, with semi-axes k√λᵢ. The eigenvalues are the variances along the principal axes, and the determinant measures the volume.
3. Mahalanobis distance counts distance in units of the standard deviation in the direction of the point, and so differs from the Euclidean distance whenever the covariance is not round.
4. In two dimensions the 1σ, 2σ, 3σ ellipses hold 39%, 86% and 98.9%, from 1 − e^(−k²/2). The 95% ellipse is 1.96σ, 2.45σ, 2.80σ in 1, 2, 3 dimensions.
5. Gate with the chi-square quantile for the dimension of the measurement, not a fixed number of sigmas.

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.

• T. D. Barfoot, State Estimation for Robotics, 2nd edition, Cambridge University Press, 2024 , chapter “Primer on Probability Theory” (Gaussian probability density functions). Free to read as the author’s PDF; the book is copyrighted.
• J. Tsitsiklis and P. Jaillet, MIT RES.6-012 Introduction to Probability , lectures 8 to 12 (continuous random variables, derived distributions, covariance). Free videos and notes, licence CC BY-NC-SA 4.0.
• S. Thrun, W. Burgard and D. Fox, Probabilistic Robotics, MIT Press, 2005, section 2.2 (the multivariate normal) and section 3.2. Paid book.
• P. C. Mahalanobis, “On the generalised distance in statistics”, Proceedings of the National Institute of Sciences of India 2 (1), 1936. The distance named after the author. Paid or archive access.
• K. P. Murphy, Probabilistic Machine Learning: An Introduction, MIT Press, 2022 , chapter 3 (the multivariate Gaussian). Free to read as the draft PDF, licence CC BY-NC-ND.

Knowledge check

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

What’s next

We can now describe the uncertainty of a vector. The next lesson asks what happens to that description when the vector is transformed: rotated, scaled, converted from polar to Cartesian coordinates, or moved by the robot.

Related content

Nothing below is required to finish this lesson.

• Before
  M4-T2: Covariance and Correlation

  Compute covariance and correlation from data, read their signs and sizes physically, and separate correlation from dependence.

• Lab
  M4-L1: Robot Uncertainty Visualizer

  The position estimate comes with a covariance that nobody on the team can read. Compute it from samples, draw its ellipse, and push it through a rotation and a range-bearing conversion.

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

• UAV M3-T2: Gaussian distributions and covariance
• SLAM M2-L3: Understand Pose Uncertainty
• Navigation M7-T3: Navigation Under Pose Uncertainty