A GPS fix isn't a dot on the map, it's a blob, and the blob has a shape. The multivariate Gaussian is the math that describes that shape, fuses two of them into one, and tells you exactly how many “sigmas” away a reading has to be before you call it impossible.
Open the raw data from your phone's GPS chip and you'll find two numbers sitting next to your location: a latitude/longitude, and a field called accuracy, reported in meters. Most people read that as “you're somewhere inside a circle of this radius.” That's already useful. But it's also wrong in a way that matters, and untangling why is the entire point of this lesson.
A single accuracy number describes a circle, equal uncertainty in every direction. Real GPS error isn't circular. The receiver estimates your position by triangulating from several satellites, and the satellites overhead are never arranged symmetrically. When the satellites you can see are clustered more toward the north-south sky than east-west, your east-west position gets pinned down tightly while your north-south position stays sloppy. The error isn't a circle; it's a tilted oval, wide in one direction and narrow in another, and the tilt itself carries information, it's telling you exactly which direction the measurement doesn't trust.
Here's the phrase that will anchor this whole lesson: a measurement's uncertainty isn't a number, it's a shape. A single sigma (a standard deviation) is enough to describe uncertainty along one axis, like “my speedometer is accurate to plus or minus 2 km/h.” The moment you have two or more correlated quantities, x and y position, temperature and humidity, this stock's return and that stock's return, the uncertainty stops being a single number and becomes a matrix. That matrix is called the covariance matrix, written Σ (capital Greek sigma), and the shape it draws is an ellipse. This lesson is about learning to read that matrix the way you'd read a speedometer: at a glance, telling direction, and size, and whether two blobs actually agree with each other.
Why should an engineer care about the tilt and not just the size? Because the tilt is where the free information lives. Picture a delivery drone that only trusts its GPS's east-west reading and ignores the correlation with north-south. If the true error cloud is a 45° diagonal cigar, large error, but almost entirely in one specific diagonal direction, a filter that ignores the tilt will systematically either panic over safe readings or accept dangerous ones, because it's using the wrong shape to decide what's plausible. Every optimal sensor fusion algorithm you'll ever use, including the Kalman filter, is doing nothing more exotic than correctly tracking this ellipse as it changes shape over time.
This lesson builds that machinery from the ground up. We start with the familiar single-variable bell curve (Chapter 1), then add a second correlated variable and watch the bell curve become a tilted dome sitting over an ellipse (Chapter 2). From there we learn to measure “how many sigmas away” a point is once the world is tilted, the Mahalanobis distance (Chapter 3), and how learning one variable narrows your uncertainty about a correlated one, which is secretly what every Kalman update is doing (Chapter 4). We'll transform ellipses through matrices (Chapter 5), generate our own correlated random samples from scratch (Chapter 6), and combine two independent uncertain estimates into one sharper belief using nothing but algebra (Chapter 7). Chapter 8 is the payoff: a full two-sensor fusion lab you can drag around with your mouse. Chapter 9 is your map for choosing the right tool and where to go next.
Let's make the blob concrete before touching any formula. Below, imagine 200 independent GPS fixes taken at the exact same spot, the true position never moves, only the receiver's noise does. Drag the sliders for the spread in the east-west direction, the spread in the north-south direction, and the correlation between them, and watch the cloud of dots change shape. At correlation zero with equal spreads you get a fuzzy circle. Push the correlation up and the cloud leans into a diagonal cigar, even though nothing about the individual x-spread or y-spread changed.
The orange dot is the true position. Each blue dot is one noisy fix. Drag the sliders; click “Resample” to draw a fresh 200 points from the current shape.
Watch what happens at the extremes. Set both spreads equal and correlation to zero: the cloud is a fuzzy circle, because error in the east-west direction tells you nothing about error in the north-south direction, they're independent. Now push correlation to +0.9 while keeping the spreads equal: the circle collapses into a tight diagonal line running from bottom-left to top-right. A large positive error in x now almost always comes bundled with a large positive error in y. Flip the correlation to −0.9 and the cigar rotates to run from top-left to bottom-right instead. The spreads never changed, only the correlation did, and yet the shape of “where the true position plausibly is” changed completely.
Put real numbers on it. Suppose your phone's error is tight in the east-west direction, σx = 5 m, but sloppy in the north-south direction, σy = 15 m, with no correlation between them. A reading that jumps 12 meters east of your last fix is 12 / 5 = 2.4σx away, a big, surprising jump on the tight axis, the kind of thing worth flagging. A reading that jumps the same 12 meters to the north is only 12 / 15 = 0.8σy away, unremarkable noise on the sloppy axis. Same 12-meter jump, wildly different level of surprise, purely because the two axes don't trust each other equally. That asymmetry is exactly what a single “accuracy: 12 m” number destroys, and exactly what the covariance matrix preserves.
Before we can tilt an ellipse, we need to be completely fluent with the un-tilted, single-axis case: the ordinary bell curve. You've almost certainly seen it before, but this chapter nails down the exact vocabulary we'll reuse for the rest of the lesson, so don't skip it even if it looks familiar.
A single measurement, a speedometer reading, a resistor's resistance, a delivery time, that's corrupted by lots of small, independent sources of noise tends to land close to some central value most of the time, and far from it rarely, in a very specific bell-shaped pattern. That pattern is the Gaussian distribution (also called the normal distribution), and it's completely described by just two numbers: the mean μ (mu), which is where the peak sits, and the variance σ2 (sigma-squared), which controls how wide the bell is. The square root of the variance, σ (sigma), is the standard deviation, it's in the same units as the original measurement, which is why people quote σ (“accurate to ±2 km/h”) far more often than σ2.
Read the formula in two pieces. The exponent, −(x−μ)2/2σ2, is the important part: it measures how far x is from the mean, in units of σ, then squares that distance and flips its sign. Far from the mean, the squared distance is large and positive, so the exponent is a large negative number, and exp() of a large negative number is close to zero, the curve is low out there. Near the mean, the exponent is close to zero, and exp(0) = 1, the curve peaks there. The front piece, 1/(σ√(2π)), is just a constant that scales the whole curve so the total area underneath sums to exactly 1, because f(x) is a probability density, not a probability by itself.
The single most useful fact about the Gaussian, the one you should be able to recite instantly, is the 68/95/99.7 rule: about 68% of the area under the curve lies within one σ of the mean, about 95% lies within two σ, and about 99.7% lies within three σ. This isn't a coincidence you need to derive from scratch every time, it's a fixed geometric property of this exact bell shape, and it's worth memorizing because it turns “how many sigma away” into an instant gut-check on how surprising something is.
Say a car's cruise control holds speed with mean μ = 60 km/h and standard deviation σ = 2 km/h. First question: what fraction of the time is the actual speed between 58 and 62 km/h? Convert the boundaries into σ units first, this conversion is called a z-score: z = (x − μ) / σ. For x = 58: z = (58 − 60) / 2 = −1. For x = 62: z = (62 − 60) / 2 = +1. So “between 58 and 62” is exactly “within one σ of the mean,” and the 68/95/99.7 rule says that's about 68% of the time.
Second question: how surprising is a reading of 65 km/h? z = (65 − 60) / 2 = 2.5. That's two-and-a-half standard deviations out, beyond the 95% band (which stops at z = 2) but not yet at the 99.7% band (z = 3). A quick lookup (or a calculator's normal CDF) says only about 0.6% of readings should land this far out or farther. A single 65 km/h reading isn't proof the cruise control is broken, but if you see several readings out here, something changed.
The shaded bands mark 1σ and 2σ from the mean. Click on the axis to drop a marker and see its z-score and which band it falls in.
Notice what happens as you drag σ up: the peak drops and the curve flattens, because the total area underneath must always stay exactly 1, spreading the same unit of “probability mass” over a wider range necessarily lowers the peak. This trade-off, height for width, is the entire geometry of a 1-D Gaussian, and it's about to become far richer once we let a second variable in.
python import math def gaussian_pdf(x, mu, sigma): # the two-piece formula, written out explicitly coeff = 1.0 / (sigma * math.sqrt(2 * math.pi)) expo = -((x - mu) ** 2) / (2 * sigma ** 2) return coeff * math.exp(expo) def z_score(x, mu, sigma): return (x - mu) / sigma mu, sigma = 60.0, 2.0 print(z_score(65, mu, sigma)) # 2.5 -- two-and-a-half sigma out print(gaussian_pdf(60, mu, sigma)) # peak density, about 0.1995 # the library equivalent -- same numbers, one line each from scipy.stats import norm print(norm.pdf(60, loc=mu, scale=sigma)) # matches gaussian_pdf print((65 - mu) / sigma) # matches z_score
Now let a second variable in. Instead of tracking only your east-west GPS error, track east-west and north-south together. Each still has its own mean and variance, but there's a new question we couldn't even ask with one variable: when x drifts high, does y tend to drift high too, drift low, or not care at all? The number that answers this is the covariance.
Read it as: for every sample, multiply how far x is above its mean by how far y is above its mean, then average that product over all samples. If x and y tend to be above their means together (or below together), the products are mostly positive and the average, the covariance, comes out positive. If one tends to be above its mean exactly when the other is below, the products are mostly negative and the covariance is negative. If the two have no relationship, positive and negative products cancel out and the covariance lands near zero.
Covariance has one annoying property: its size depends on the units of x and y, so “a covariance of 8” means nothing on its own until you know whether x and y are in millimeters or kilometers. The fix is to normalize by each variable's own spread, giving the correlation coefficient:
ρ (rho) is always between −1 and +1 no matter what units you started in. ρ = 0 means no linear relationship; ρ = +1 means x and y move in perfect lockstep; ρ = −1 means they move in perfect opposition. This is the same ρ you dragged in Chapter 0's slider, now you know exactly what number it is.
Stack the variance of x, the variance of y, and their covariance into a single 2×2 table and you get the covariance matrix Σ, the object this entire lesson is about:
Σ is always symmetric (the top-right entry equals the bottom-left entry, because Cov(X,Y) = Cov(Y,X)) and its diagonal entries are always non-negative (a variance can't be negative). The joint distribution over x and y is now a surface sitting above the plane, a dome, taller where probability is denser, and if you slice that dome horizontally at any height, the slice you get is an ellipse. Σ is exactly what tells you the size, shape, and tilt of every one of those ellipses.
Let's use a concrete covariance matrix and find its ellipse by hand, the same way Chapter 2 of the Eigenvalues lesson found stretch directions for a general matrix, because that's exactly what we're doing here.
Read off the plain-English numbers first: σx2 = 4 so σx = 2; σy2 = 4 so σy = 2; the covariance is 2, so ρ = 2 / (2·2) = 0.5. Two equally-spread variables with moderate positive correlation. If Σ were diagonal (covariance zero), the ellipse would be a circle aligned with the axes. The off-diagonal 2 is exactly what tilts it.
To find the tilt precisely, we need Σ's eigenvalues and eigenvectors, the special directions Σ only stretches, never rotates (exactly the definition from the Eigenvalues lesson). The characteristic equation is det(Σ − λI) = 0. For a 2×2 matrix, the two shortcuts are: eigenvalues sum to the trace (4 + 4 = 8) and multiply to the determinant (4·4 − 2·2 = 16 − 4 = 12). So we need two numbers that sum to 8 and multiply to 12: solving λ2 − 8λ + 12 = 0 with the quadratic formula gives λ = (8 ± √(64−48)) / 2 = (8 ± √16) / 2 = (8 ± 4) / 2, so λ1 = 6 and λ2 = 2. Check: 6 + 2 = 8 ✓ and 6 × 2 = 12 ✓.
Now the directions. For λ1 = 6, solve (Σ − 6I)v = 0: subtracting 6 from the diagonal gives [[−2, 2], [2, −2]], and this matrix sends v = (1, 1) to (−2+2, 2−2) = (0, 0), exactly what “equals zero” means, so v1 = (1, 1), the 45° diagonal. For λ2 = 2, solve (Σ − 2I)v = 0: this gives [[2, 2], [2, 2]], which sends v = (1, −1) to (2−2, 2−2) = (0, 0), so v2 = (1, −1), the 135° diagonal. The two eigenvectors are perpendicular, as they always are for a symmetric matrix like Σ.
Translate the eigen-decomposition into the picture: the ellipse's long axis points along (1, 1), the 45° line, with a semi-axis length proportional to √λ1 = √6 ≈ 2.449. Its short axis points along (1, −1), the 135° line, with semi-axis length proportional to √λ2 = √2 ≈ 1.414. Notice the semi-axis uses √λ, not λ itself, because λ is a variance along that direction, and a length needs a standard deviation, which is variance's square root. This is the single most important geometric fact in the whole lesson: the eigenvectors of Σ are the ellipse's axes, and the square roots of the eigenvalues are how long those axes are.
Drag the sliders. The orange and teal lines are the eigenvector axes, scaled by √λ. Watch them stay perpendicular no matter how you tilt the ellipse.
python import numpy as np Sigma = np.array([[4, 2], [2, 4]]) eigvals, eigvecs = np.linalg.eigh(Sigma) # eigh: Sigma is symmetric print(eigvals) # [2. 6.] -- ascending order print(eigvecs) # columns are the unit eigenvectors, (1,-1)/sqrt2 then (1,1)/sqrt2 axis_lengths = np.sqrt(eigvals) # [1.414, 2.449] -- semi-axis lengths
Chapter 1 gave us a clean way to ask “how surprising is this value” along a single axis: compute the z-score, (x−μ)/σ. Chapter 2 built a tilted ellipse out of two correlated variables. Now we need the two-variable version of “how many sigma away,” and the honest answer is that ordinary straight-line distance gets this badly wrong once the world is tilted.
Here's the trap. Take the ellipse from Chapter 2, Σ = [[4, 2], [2, 4]], centered at the origin. Consider two points, both exactly √10 ≈ 3.16 units from the center in ordinary Euclidean distance: P = (3, 1) and Q = (3, −1). Same straight-line distance from the mean. But P sits almost exactly along the ellipse's long axis (the 45° direction, where the data is spread wide and such a jump is unremarkable), while Q sits closer to the short axis (the 135° direction, where the data is tight and the same jump is a much bigger deal). Treating P and Q as equally surprising, just because a ruler says they're equidistant, throws away everything Chapter 2 taught us about the shape.
The fix is the Mahalanobis distance, named after statistician P. C. Mahalanobis, and it's a beautifully small edit to the ordinary distance formula:
Compare this to ordinary squared distance, which is (x−μ)T(x−μ), the same expression but with an identity matrix silently standing in for Σ−1. Inserting Σ−1 (sigma-inverse) does the job of “dividing out” the scale and correlation before measuring distance: it stretches the cheap (high-variance) directions back down and shrinks the expensive (low-variance) directions back up, until the ellipse becomes a perfect circle again. Once the space is un-tilted like that, plain distance on the corrected coordinates is exactly the right notion of “how many sigma away,” which is precisely what Mahalanobis distance measures.
First we need Σ−1 for Σ = [[4, 2], [2, 4]]. For a 2×2 matrix [[a,b],[c,d]], the inverse is (1/det)·[[d,−b],[−c,a]], where det = ad−bc. Here det = 4·4 − 2·2 = 16−4 = 12 (the same determinant from Chapter 2, no surprise, since det(Σ) is the product of the eigenvalues, 6×2 = 12). So:
Now point P = (3, 1). First multiply Σ−1 by P: top entry = ⅓×3 + (−⅝)×1 = 1 − ⅝ = ⅝×5 = 5⁄6; bottom entry = (−⅝)×3 + ⅓×1 = −½ + ⅓ = −1⁄6. So Σ−1P = (5⁄6, −1⁄6). Then D2 = P · (Σ−1P) = 3×5⁄6 + 1×(−1⁄6) = 15⁄6 − 1⁄6 = 14⁄6 = 7⁄3 ≈ 2.333. So DP = √2.333 ≈ 1.53σ.
Now point Q = (3, −1). Top entry of Σ−1Q: ⅓×3 + (−⅝)×(−1) = 1 + ⅝ = 7⁄6. Bottom entry: (−⅝)×3 + ⅓×(−1) = −½ − ⅓ = −5⁄6. So Σ−1Q = (7⁄6, −5⁄6). Then D2 = 3×7⁄6 + (−1)×(−5⁄6) = 21⁄6 + 5⁄6 = 26⁄6 = 13⁄3 ≈ 4.333. So DQ = √4.333 ≈ 2.08σ.
There it is, in real numbers: P and Q are the exact same Euclidean distance (√10 ≈ 3.16) from the mean, but P is only 1.53 sigma away while Q is 2.08 sigma away, a 36% larger surprise, purely because Q's direction fights the correlation instead of riding along with it. This is the number every anomaly detector, every outlier gate, and every consistency check in this lesson is built on.
How far is “too far”? For an n-dimensional Gaussian, D2 follows a chi-squared distribution with n degrees of freedom. In 2-D (n = 2), the value that 95% of genuine samples fall under is χ20.95,2 ≈ 5.991. Both P (2.33) and Q (4.33) are comfortably under that gate, both are plausible draws from this Gaussian. But take a third point, R = (4, −2): repeating the same recipe gives Σ−1R = (5⁄3, −4⁄3), so D2 = 4×5⁄3 + (−2)×(−4⁄3) = 20⁄3 + 8⁄3 = 28⁄3 ≈ 9.33, which clears the 5.991 gate, R is a 95%-confidence outlier, even though its Euclidean distance from the origin (√20 ≈ 4.47) isn't dramatically larger than P and Q's.
The dashed rings are Mahalanobis-distance contours for Σ = [[4,2],[2,4]] (1σ, 2σ, and the 95% chi-squared gate). Drag the dot; the readout shows both distances live.
python, plain Python, no libraries def mahalanobis_2d(x, y, Sigma): a, b = Sigma[0][0], Sigma[0][1] c, d = Sigma[1][0], Sigma[1][1] det = a * d - b * c if abs(det) < 1e-12: raise ValueError("Sigma is singular") # Sigma^-1 = (1/det) * [[d,-b],[-c,a]], applied directly to (x,y) inv_x = (d * x - b * y) / det inv_y = (-c * x + a * y) / det return x * inv_x + y * inv_y # this is D^2 Sigma = [[4, 2], [2, 4]] print(mahalanobis_2d(3, 1, Sigma)) # 2.333... print(mahalanobis_2d(3, -1, Sigma)) # 4.333...
numpy import numpy as np def mahalanobis(x, mu, Sigma): diff = x - mu return float(diff @ np.linalg.inv(Sigma) @ diff) Sigma = np.array([[4., 2.], [2., 4.]]) mu = np.zeros(2) print(mahalanobis(np.array([3., 1.]), mu, Sigma)) # 2.333... print(mahalanobis(np.array([3., -1.]), mu, Sigma)) # 4.333...
pytorch import torch def mahalanobis_t(x, mu, Sigma): diff = (x - mu).unsqueeze(-1) # column vector, shape (2,1) Sigma_inv = torch.linalg.inv(Sigma) return (diff.T @ Sigma_inv @ diff).squeeze() Sigma = torch.tensor([[4., 2.], [2., 4.]]) mu = torch.zeros(2) print(mahalanobis_t(torch.tensor([3., 1.]), mu, Sigma)) # tensor(2.333...) print(mahalanobis_t(torch.tensor([3., -1.]), mu, Sigma)) # tensor(4.333...) -- autograd-ready
Suppose you have the correlated x-y blob from Chapter 2 and you get to observe x exactly, a wheel encoder tells you the robot's forward distance precisely, and you want to know what that implies about a correlated quantity y, say lateral drift, that you haven't measured. This chapter is about the two different ways of asking “what happens to y,” and why they give very different answers.
The first way is the marginal distribution of y: ignore x entirely and ask what y looks like on its own, averaged over every possible value of x. Geometrically, this is squashing the 2-D dome flat onto the y-axis, like shining a flashlight through the blob from the side and looking at the shadow it casts. The marginal of y is just a plain 1-D Gaussian with mean μy and variance σy2, read straight off the diagonal of Σ. Correlation plays no role here at all, the shadow is the same shape no matter how tilted the blob is, because squashing along one axis erases the tilt.
The second way is the conditional distribution of y given a specific, known value of x: not the whole shadow, but a single slice through the dome at exactly x = c. Picture cutting a loaf of bread at one point and looking at the cross-section, rather than flattening the whole loaf. That cross-section is also a Gaussian, a wonderful, non-obvious fact about the multivariate Gaussian family, but it has a different mean and, crucially, a smaller variance than the marginal, because knowing x removes exactly the part of y's uncertainty that the correlation could explain.
The formulas for the conditional mean and variance follow directly from Σ, and they are worth deriving in words before symbols. If ρ is positive and you observe x above its mean, that pulls the conditional mean of y upward too, by an amount that scales with how strongly correlated they are (ρ) and how the two variables' natural scales compare (σy/σx). And however you observe x, the conditional variance shrinks by a factor of (1 − ρ2), regardless of which specific value x took:
Notice the variance formula doesn't depend on which value of x you observed at all, only on how correlated x and y are. That's a genuinely surprising property of the Gaussian family: how much you learn about y depends only on ρ, while what you learn (the shifted mean) depends on the actual observed value of x.
Use the familiar Σ = [[4, 2], [2, 4]] with μ = (0, 0), so σx = σy = 2 and ρ = 0.5. Suppose we observe x = 3. Plugging in: μy|x=3 = 0 + 0.5 × (2/2) × (3 − 0) = 0.5 × 1 × 3 = 1.5. So learning x = 3 shifts our best guess for y from 0 (the unconditional mean) up to 1.5, the correlation drags y along partway, but not all the way to 3, because ρ = 0.5 is a partial relationship, not a perfect one.
Now the variance: σ2y|x = 4 × (1 − 0.52) = 4 × (1 − 0.25) = 4 × 0.75 = 3, so the conditional standard deviation is √3 ≈ 1.73, down from the marginal σy = 2. Observing x shrank our uncertainty about y from 2 to 1.73, a real but modest 13.5% reduction, exactly what you'd expect from a moderate correlation of 0.5. Crank ρ up toward 1 and the shrink factor (1−ρ2) races to zero, a near-perfect correlation means knowing x tells you almost everything about y. At ρ = 0, the factor is exactly 1: an uncorrelated x tells you nothing, and the conditional equals the marginal, which is exactly the geometric intuition from Chapter 2's circle case.
Top: the 2-D blob with a draggable conditioning line at x = c. Bottom: the marginal (shadow, unconditioned) distribution of y versus the conditional (slice) distribution of y given x = c.
Every ellipse we've drawn so far has been sitting still. Real systems move: a robot predicts where it'll be one time-step from now given where it is now and how fast it's going; a camera projects a 3-D point into 2-D pixels; a unit conversion turns meters into feet. Each of these is a linear transform, multiply by a matrix, maybe add a constant, and this chapter answers the question every one of those systems needs answered: if I know the uncertainty before the transform, what's the uncertainty after it?
Say x is Gaussian with mean μ and covariance Σ, and you compute y = Ax + b for some matrix A and constant vector b. The mean transforms exactly as you'd hope, because expectation is linear: E[Ax+b] = A·E[x] + b = Aμ + b. The covariance is less obvious, but derivable from the definition of covariance and the algebra of expectations, and the result is a clean sandwich:
Why a matrix on both sides rather than, say, just AΣ or an elementwise A·A applied to each variance? Because covariance measures how pairs of components move together, and a linear transform mixes every input component into every output component. The A on the left carries the transform's effect on the “first copy” of the deviation, and the AT on the right carries its effect on the “second copy”, remember covariance is built from a product of two deviations, (x−μ) times itself, so both factors need to pass through the same transform.
Take a state x = (position, velocity) with Σ = [[4, 2], [2, 4]], position variance 4, velocity variance 4, and a covariance of 2 between them (perhaps because both were estimated from the same noisy trajectory fit). A constant-velocity motion model over one time unit says: new position = old position + old velocity, new velocity = old velocity, unchanged. In matrix form that's y = Ax with A = [[1, 1], [0, 1]].
Compute AΣ first, one entry at a time. Row 1 of A is [1, 1]; dotted with column 1 of Σ, (4, 2), gives 1×4 + 1×2 = 6; dotted with column 2, (2, 4), gives 1×2 + 1×4 = 6. Row 2 of A is [0, 1]; dotted with column 1 gives 0×4 + 1×2 = 2; dotted with column 2 gives 0×2 + 1×4 = 4. So AΣ = [[6, 6], [2, 4]].
Now finish with (AΣ)AT, where AT = [[1, 0], [1, 1]]. Row 1 of AΣ is [6, 6]; dotted with column 1 of AT, (1, 1), gives 6+6 = 12; dotted with column 2, (0, 1), gives 0+6 = 6. Row 2 of AΣ is [2, 4]; dotted with column 1 gives 2+4 = 6; dotted with column 2 gives 0+4 = 4. The result:
Read the story in these numbers. Position's variance tripled, from 4 to 12, because the new position inherited both its own old uncertainty and velocity's uncertainty and the cross-talk between them (that's exactly what “new position = old position + old velocity” means for uncertainty: three sources stack up). Velocity's variance stayed exactly 4, because the model says velocity doesn't change, nothing new was added or removed from what we already believed about it. This is not a toy exercise: it is, entry for entry, the predict step of a Kalman filter with no process noise. Compare this to Chapter 4's conditioning formula, which was the update step, together, predict-then-update is the entire Kalman filter, and both halves are this lesson's algebra wearing a different hat.
One more useful fact, needed later: if x1 and x2 are independent Gaussians, their sum x1 + x2 is Gaussian too, with mean μ1+μ2 and covariance Σ1+Σ2, covariances simply add when independent sources of randomness combine. This is exactly how a Kalman filter injects process noise Q into a prediction: it adds Q's covariance to the already-transformed AΣAT, because the true motion model has its own independent randomness layered on top of the deterministic A.
The blue ellipse is the input Σ = [[4,2],[2,4]]. Adjust A's four entries; the orange ellipse shows AΣAT live.
numpy import numpy as np Sigma = np.array([[4., 2.], [2., 4.]]) A = np.array([[1., 1.], [0., 1.]]) # constant-velocity model, dt=1 Sigma_new = A @ Sigma @ A.T print(Sigma_new) # [[12. 6.] [6. 4.]] -- matches the hand calc # independent sums add covariance directly Sigma_1 = np.array([[1., 0.], [0., 1.]]) Sigma_2 = np.array([[0.5, 0.], [0., 0.5]]) print(Sigma_1 + Sigma_2) # [[1.5 0.] [0. 1.5]] -- process noise added to a prediction
Every chapter so far has worked with the analytic ellipse. Sometimes you need the opposite: actual random samples with a given covariance, for example to feed synthetic sensor noise into a filter you're testing before it ever sees a real robot, or to run a Monte Carlo check that your formulas actually behave the way the math says they should. This chapter builds correlated samples from nothing but a source of plain, independent, standard-normal randomness.
“Standard normal” means mean 0, variance 1, no correlation with anything, the easiest Gaussian to generate (the Box-Muller transform, used elsewhere on this site, turns two uniform random numbers into one standard normal draw). Call a pair of independent standard normal numbers z = (z1, z2). The question is: what linear transform turns z into a sample x with our target mean μ and covariance Σ?
Chapter 5 already answered exactly this question in reverse: if z has covariance I (the identity, no correlation, unit variance in every direction) and we compute x = μ + Lz for some matrix L, then x's covariance is L·I·LT = LLT. So we need a matrix L with LLT = Σ. That's called a Cholesky decomposition, and for a symmetric positive-definite matrix (which every valid covariance matrix is) it always exists and is unique if we insist L is lower-triangular, zero above the diagonal.
Solve for L's three entries against Σ = [[4, 2], [2, 4]] by multiplying L by LT symbolically and matching entries. The (1,1) entry of LLT is l112, which must equal 4, so l11 = 2. The (2,1) entry is l21·l11, which must equal 2, so l21 = 2/2 = 1. The (2,2) entry is l212 + l222, which must equal 4, so 1 + l222 = 4, l222 = 3, l22 = √3 ≈ 1.732. So L = [[2, 0], [1, 1.732]]. Quick sanity check on the (1,2) entry (should be 0 since L is lower-triangular and LT is upper-triangular): l11·0 + 0·l22 = 0. ✓
Now draw two samples. For z = (1, 0): x = μ + Lz = (0,0) + (2×1 + 0×0, 1×1 + 1.732×0) = (2, 1). For z = (0, 1): x = (0,0) + (2×0 + 0×1, 1×0 + 1.732×1) = (0, 1.732). Notice the second sample got a nonzero y even though z's y-component (1) never touched a nonzero entry of the first column of L, the l21 entry is exactly what leaks a bit of z1's influence into y when z1 is nonzero, and it's precisely this leak that creates correlation between the two output components. Two independent inputs; one correlated pair out. That's what the Cholesky factor buys you.
Each click draws 20 fresh Cholesky-transformed samples. The dashed ellipse is the true Σ = [[4,2],[2,4]]; the readout compares it to the empirical covariance measured from every sample drawn so far.
numpy import numpy as np def cholesky_2x2(Sigma): a, b, d = Sigma[0,0], Sigma[1,0], Sigma[1,1] l11 = np.sqrt(a) l21 = b / l11 l22 = np.sqrt(d - l21**2) return np.array([[l11, 0.], [l21, l22]]) Sigma = np.array([[4., 2.], [2., 4.]]) L = cholesky_2x2(Sigma) # [[2. 0.] [1. 1.732]] -- matches the hand calc print(np.allclose(L, np.linalg.cholesky(Sigma))) # True -- NumPy agrees rng = np.random.default_rng(0) z = rng.standard_normal((10000, 2)) # 10,000 independent standard normals mu = np.array([0., 0.]) x = mu + z @ L.T # correlated samples, vectorized print(np.cov(x.T)) # converges to [[4,2],[2,4]] as N grows
Here's the question every sensor-fusion system has to answer: you have two independent, uncertain estimates of the same quantity, a GPS fix and an IMU-predicted position, say, each shaped like its own Gaussian blob. How do you combine them into a single, better estimate, and can you prove it's actually better rather than just a reasonable-sounding average?
The answer comes from a genuinely beautiful fact: if you multiply two Gaussian probability densities together as functions of x (not caring yet that the product isn't itself a normalized probability), the result is proportional to another Gaussian. Multiplying two bell-shaped curves gives you back a bell-shaped curve, a closure property most probability distributions don't have. This is precisely what Bayes' Rule is doing when the prior and the likelihood are both Gaussian: posterior ∝ prior × likelihood, and Gaussian × Gaussian ∝ Gaussian, so the whole update stays inside the Gaussian family forever. No numerical integration required, ever.
The cleanest way to see how the two combine is to stop working with covariance and switch to its reciprocal, the precision matrix, Λ = Σ−1 (this is exactly the matrix Chapter 3 needed for Mahalanobis distance). Precision measures confidence directly: a small covariance (tight, confident belief) means a large precision, and vice versa. When you multiply two independent Gaussians together, their precisions simply add, and their precision-weighted means add too:
Read the precision-addition rule as: confidence is additive. Two independent sources of information are always at least as informative together as either one alone, you can never become less certain by adding an independent measurement, no matter how noisy that measurement is. This single fact is the mathematical reason sensor fusion works at all.
Take a prior belief about position from dead-reckoning, vague and isotropic: μ1 = (0, 0), Σ1 = [[4, 0], [0, 4]] (standard deviation 2 in both directions, no correlation). Now a GPS fix comes in, more confident: μ2 = (3, 1), Σ2 = [[1, 0], [0, 1]] (standard deviation 1, four times more precise). Because both are diagonal, inverting is just reciprocating each diagonal entry: Λ1 = [[0.25, 0], [0, 0.25]] and Λ2 = [[1, 0], [0, 1]].
Add the precisions: Λpost = [[0.25+1, 0], [0, 0.25+1]] = [[1.25, 0], [0, 1.25]]. Invert back to covariance (reciprocate again, since it's still diagonal): Σpost = [[1/1.25, 0], [0, 1/1.25]] = [[0.8, 0], [0, 0.8]]. Check the headline claim: 0.8 is smaller than both inputs' variances (4 and 1), fusing two independent estimates always leaves you strictly more confident than your best individual estimate, not merely as confident as the better one.
For the mean: Λ1μ1 = [[0.25,0],[0,0.25]]·(0,0) = (0, 0). Λ2μ2 = [[1,0],[0,1]]·(3,1) = (3, 1). Sum: (0,0) + (3,1) = (3, 1). Finally μpost = Σpost·(3,1) = (0.8×3, 0.8×1) = (2.4, 0.8). Sanity check the weighting: the posterior mean sits exactly 80% of the way from the prior's mean (0,0) toward the measurement's mean (3,1), and 80% is exactly the fraction of total precision that the four-times-sharper measurement contributes (1 / 1.25 = 0.8). The vague prior only gets a 20% vote, precisely because it earned only 20% of the combined confidence.
python, plain Python, no libraries (diagonal case) def fuse_diag(mu1, var1, mu2, var2): # works per-axis since both covariances are diagonal here prec1, prec2 = 1.0 / var1, 1.0 / var2 prec_post = prec1 + prec2 var_post = 1.0 / prec_post mu_post = var_post * (prec1 * mu1 + prec2 * mu2) return mu_post, var_post for axis in [(0.0, 4.0, 3.0, 1.0), (0.0, 4.0, 1.0, 1.0)]: mu1, var1, mu2, var2 = axis print(fuse_diag(mu1, var1, mu2, var2)) # (2.4, 0.8) then (0.8, 0.8) -- matches the hand calc, axis by axis
numpy (general, non-diagonal case) import numpy as np def fuse_gaussians(mu1, Sigma1, mu2, Sigma2): Lam1, Lam2 = np.linalg.inv(Sigma1), np.linalg.inv(Sigma2) Sigma_post = np.linalg.inv(Lam1 + Lam2) mu_post = Sigma_post @ (Lam1 @ mu1 + Lam2 @ mu2) return mu_post, Sigma_post mu1, Sigma1 = np.array([0.,0.]), np.eye(2) * 4. mu2, Sigma2 = np.array([3.,1.]), np.eye(2) * 1. mu_post, Sigma_post = fuse_gaussians(mu1, Sigma1, mu2, Sigma2) print(mu_post) # [2.4 0.8] print(Sigma_post) # [[0.8 0.] [0. 0.8]]
pytorch import torch def fuse_gaussians_t(mu1, Sigma1, mu2, Sigma2): Lam1, Lam2 = torch.linalg.inv(Sigma1), torch.linalg.inv(Sigma2) Sigma_post = torch.linalg.inv(Lam1 + Lam2) mu_post = Sigma_post @ (Lam1 @ mu1 + Lam2 @ mu2) return mu_post, Sigma_post mu1, Sigma1 = torch.zeros(2), torch.eye(2) * 4. mu2, Sigma2 = torch.tensor([3., 1.]), torch.eye(2) mu_post, Sigma_post = fuse_gaussians_t(mu1, Sigma1, mu2, Sigma2) print(mu_post, Sigma_post) # tensor([2.4, 0.8]) -- gradients flow through if mu1/Sigma1 require_grad
The blue ellipse is the prior, the orange ellipse is the measurement, and the green ellipse is the precision-weighted posterior.
Every idea in this lesson comes together here. Below are two independent sensor estimates of the same 2-D quantity, call them Sensor A and Sensor B, each a full ellipse with its own center, spread, and tilt, exactly the objects Chapter 2 taught you to read. Drag either ellipse's center to move it. Use the sliders to stretch, squash, and rotate each one. Watch three things happen automatically, live, every time you touch a control.
First, the lab computes the Mahalanobis distance between the two sensors' means, not a point and a mean this time, but two whole estimates asking each other “are we even talking about the same place?” It does this exactly the way Chapter 3 taught, except the covariance used for the gate is ΣA + ΣB, the combined spread of “how far apart these two readings could plausibly be by chance alone” (this follows from Chapter 5's rule that independent covariances add). Second, it checks that distance against the 95% chi-squared gate from Chapter 3 and shows a clear accept/reject verdict. Third, whenever the gate passes, it computes the precision-weighted posterior from Chapter 7, the fused estimate, and draws it as a third ellipse, always smaller than either input.
Blue = Sensor A. Orange = Sensor B. Green = fused posterior (only drawn when the gate passes, unless you force it).
Play the “break it” scenario. Click Reset to a consistent pair first: Sensor A sits near (0,0) with ΣA = [[4,0],[0,4]] (in its own un-rotated form), Sensor B sits nearby with a smaller, rotated ellipse. Compute the gate by hand: with A's mean near the origin and B's mean only a couple of units away, the Mahalanobis distance between the means, measured against ΣA+ΣB, lands comfortably under 5.991, the gate passes, the green posterior ellipse appears, and it's visibly smaller than both inputs, exactly as Chapter 7 predicts.
Now click Break it. Sensor B's center jumps to somewhere the two ellipses barely overlap. Redo the same hand-calc shape: say A stays at μA = (0,0) with ΣA = [[4,0],[0,4]], and B jumps to μB = (6,6) with ΣB = [[1,0],[0,1]]. The difference of means is (6,6); the combined covariance is ΣA+ΣB = [[5,0],[0,5]], whose inverse is [[0.2,0],[0,0.2]]. The Mahalanobis distance squared is 6×(6×0.2) + 6×(6×0.2) = 7.2 + 7.2 = 14.4, more than double the 5.991 gate. The lab flags this pair as inconsistent and refuses to fuse by default.
One more thing worth dragging around: rotate Sensor A and Sensor B to very different angles, so their long axes point in different directions, and watch the fused posterior's shape. It's neither input's shape, it's typically rounder and smaller than both, because each sensor is confident along a different direction, and fusing them fills in each other's blind spot. This is the deepest reason multi-sensor systems out-perform any single sensor: two ellipses tilted differently cover more of the plane with high confidence than either one alone, and the precision-addition math from Chapter 7 captures that automatically, you never have to reason about it by hand.
Ten chapters, one covariance matrix Σ running through every single one of them. This closing chapter is a map: which tool answers which question, and where each of them shows up again once you leave this lesson.
| Tool | Core question it answers | Cares about correlation? | Reach for it when… |
|---|---|---|---|
| Single Gaussian (1-D) | How surprising is one number? | No, only one axis exists | You have one noisy scalar and need a z-score or a confidence interval. |
| Covariance Σ + ellipse | What shape is my uncertainty? | Yes, ρ is the tilt | Two or more correlated quantities need one confidence region, not separate error bars. |
| Mahalanobis distance | How many sigmas away, given the tilt? | Yes | Deciding whether a point is a plausible member of a cluster or an outlier. |
| Marginal | What do I know about y alone? | Ignores it by construction | You only ever observe y and never get to see the correlated x. |
| Conditional | What does knowing x tell me about y? | Yes, this is exactly what shrinks σ | You just observed one correlated variable and need to update belief in another. |
| Linear transform (AΣAT) | How does uncertainty move through a model? | Yes, preserves and reshapes it | Propagating uncertainty through a motion model, unit conversion, or projection. |
| Cholesky sampling | How do I generate realistic correlated noise? | Yes, injected via the off-diagonal of L | Stress-testing a filter or simulator against a specific synthetic covariance. |
| Product of Gaussians (fusion) | How do I combine two independent estimates? | Yes, both inputs' tilts carry through | Two independent sensors (or a prior and a measurement) estimate the same quantity. |
Every entry in that table has a name you've now derived by hand, not just memorized: 6 and 2 as the eigenvalues of [[4,2],[2,4]]; 1.53σ and 2.08σ as two points at the identical Euclidean distance; 1.5 and 3 as a conditional mean and variance; 12 and 4 as a covariance transformed through a motion model; 2 and 1.732 as a Cholesky factor's entries; 2.4, 0.8, and 0.8 as a fused mean and variance. None of it was magic, every number came from the same handful of 2×2 operations: determinant, inverse, eigen-decomposition, matrix multiply.
This lesson was deliberately built as the missing middle chapter between two things you may already know. Bayes' Rule tells you, in general, how to update a belief with evidence: posterior ∝ prior × likelihood. This lesson showed you exactly what happens to that rule when everything involved is Gaussian, the update collapses into closed-form algebra on means and covariances, no integrals required, which is Chapter 7's entire content. The Kalman Filter, meanwhile, is nothing but this lesson's predict step (Chapter 5's AΣAT) alternating forever with this lesson's update step (Chapter 4's conditioning, equivalently Chapter 7's fusion), if the Kalman Filter ever felt like a black box of Greek letters, it's because you hadn't yet built the covariance-ellipse intuition this lesson exists to give you.
Two lessons worth visiting next, in order. First The Kalman Filter, you'll recognize the predict step as Chapter 5 and the update step as Chapter 4/7 immediately, and the “Kalman gain” will read as exactly the ρ·(σy/σx) term from Chapter 4, generalized to matrices. Then The Extended Kalman Filter, for when the motion or measurement model isn't linear, it linearizes around the current estimate first (a local A or H matrix, called a Jacobian) and then runs the exact same Gaussian algebra this lesson built. If Bayes' Rule itself still feels shaky, Bayes' Rule is the more general foundation this whole lesson specialized.