Multivariate Gaussian & Joint Distributions

The Gaussian is the workhorse distribution of ML, and the multivariate Gaussian is just "many Gaussians that know about each other." If you understand the shape of an ellipse in 2D, you already understand the idea: every formula in this lesson generalizes that ellipse to many dimensions. You will work one small example by hand (a person's weight and height), then check it in numpy and by simulation.

Learning Objectives

After this lesson, you will be able to:

Reason about multiple random variables at once using joint, marginal, and conditional distributions, and know which one answers which question
Read and write the multivariate Gaussian PDF — including what each piece (mean vector, covariance matrix, determinant, quadratic form) is doing geometrically
Use the conditioning formula for a jointly-Gaussian pair to compute a conditional mean and a shrunk variance by hand, and see that it reproduces a ridge regression
Apply the change-of-variables formula and the reparameterization trick (sampling as a shift-and-stretch of plain noise), the two tools that make generative models trainable

Prerequisites

Before You Start

Review

Why

Why This Matters

How a self-driving car fuses GPS, lidar and an inertial sensor into one position estimate. Treat (true position, sensor reading) as one multivariate Gaussian, fix the reading you actually saw, and your uncertainty about the position shrinks. That operation is called conditioning, and you will do it by hand in this lesson. Engineers call the repeated version a Kalman filter
How generative models draw random numbers. Image generators and VAEs (networks that squeeze data into a few numbers and expand it back) sample from a multivariate Gaussian using a shift-and-stretch trick from this lesson, which lets gradients flow through the sampling step
Why a 95% ellipse on a scatter plot is an ellipse, not a circle. Its tilt and shape come directly from the eigenvectors and eigenvalues of the covariance matrix you will compute below
Build this → Build a "Sensor Fusion Toy" where you set the covariance of your "true position" and of a noisy sensor; the app draws both ellipses, then shows the much tighter ellipse after conditioning on a sensor reading. Watching the uncertainty shrink is the moment the math clicks

#From One Variable to Many

Earlier lessons taught you how one number wiggles. This one teaches you how many numbers wiggle together. That is the entire story. Almost every interesting ML system has dozens or millions of variables that depend on each other. To reason about them you need joint distributions, conditional distributions, and one specific multi-dimensional bell curve that shows up everywhere.

Intuition

Build Your Intuition

Picture a person's height and weight. Each is a random variable on its own — heights have a mean and a variance, so do weights. But they are not independent: tall people tend to weigh more. If you pick someone at random and I tell you they are 6'4", your guess about their weight should change.

That dependence cannot be captured by talking about each variable separately. You need a joint distribution that says how (height, weight) behave together. Once you have it, everything else falls out: the marginal is what you get if you ignore weight and only ask about height; the conditional is what you get if I fix the height to 6'4" and ask about weight; the expected weight given height is the linear-regression line drawn through the cloud.

Multi-variable thinking is harder than single-variable thinking. But it is also the language modern ML actually speaks, because every dataset is a cloud of correlated columns and every neural network produces a vector of correlated outputs.

Later in this lesson an interactive lets you drag the covariance and watch an ellipse tilt, once the ellipse has been explained. First, the tables.

Try it! Open the Python REPL (bottom-right of the screen: click Quick Actions, then Python) and follow along.

#Joint Distributions

For two random variables X and Y, a joint distribution says how likely each pair of outcomes is. For discrete variables, that is a joint PMF P(X = x, Y = y) — a table with one entry per pair. For continuous variables, it is a joint PDF f(x, y) — a 2D surface where probability is the volume under a region.

Here is a concrete joint PMF for "weather × mood" (totally invented numbers, but the structure is what matters):

	mood = grumpy	mood = neutral	mood = happy	row total
sunny	0.05	0.10	0.30	0.45
cloudy	0.10	0.15	0.10	0.35
rainy	0.15	0.04	0.01	0.20
column total	0.30	0.29	0.41	1.00

Every cell is a joint probability — for example P(sunny, happy) = 0.30. All cells together sum to 1, because some (weather, mood) combination happens.

P(X = x,\, Y = y) \geq 0 \quad\text{and}\quad \sum_{x}\sum_{y} P(X = x,\, Y = y) = 1

For continuous variables, the analogous statement is f(x, y) ≥ 0 and the double integral over the entire plane is 1. Probabilities of regions are volumes under the joint PDF surface.

What Do You Think?

Look at the weather × mood joint PMF above. Reading it directly, which is bigger?

#Marginal Distributions: Ignore Things You Do Not Care About

The marginal of X is what you get by collapsing across Y — adding up (or integrating out) the variable you are not asking about. The row totals in the table above are P(X) — the marginal of weather.

P(X = x) = \sum_{y} P(X = x,\, Y = y) \qquad f_X(x) = \int_{-\infty}^{\infty} f(x, y)\, dy

So P(sunny) = 0.45, P(cloudy) = 0.35, P(rainy) = 0.20. Those numbers add to 1 — they are a valid one-variable PMF.

Quick check

A 2D Gaussian over (X, Y) has Σ = [[1, 0.8], [0.8, 1]] (strong positive correlation). What does the MARGINAL distribution of X look like?

#Conditional Distributions: Once You Know Y

If I tell you it is rainy, what is the distribution over moods? You restrict attention to the rainy row, then renormalize so the row sums to 1.

f(x \mid y) = \frac{f(x, y)}{f(y)}

P(grumpy | rainy) = 0.15 / 0.20 = 0.75. P(happy | rainy) = 0.01 / 0.20 = 0.05. Rainy days are bleak in this fictional world.

This is just Bayes' theorem from the probability lesson, in another costume. P(X | Y) = P(X, Y) / P(Y) is the definition of conditional; the version with P(Y | X) P(X) / P(Y) in the numerator is what you get when you also rewrite P(X, Y) = P(Y | X) P(X).

#Conditional Expectation: A Random Variable in Its Own Right

The conditional expectation E[X | Y] is what X averages to once you know Y. Crucially, this is itself a random variable — it depends on the value of Y, which is random.

For each fixed y, E[X | Y = y] is a number. Take the rainy row of the weather table: after renormalizing it is [0.75, 0.20, 0.05]. Score the moods grumpy = -1, neutral = 0, happy = +1, and then E[mood | rainy] = 0.75(-1) + 0.20(0) + 0.05(+1) = -0.70. Vary y and you get a function. That function, viewed as a transformation of the random variable Y, is what we mean by E[X | Y].

E[X \mid Y = y] = \sum_{x} x \cdot P(X = x \mid Y = y) \qquad E[X] = E[\,E[X \mid Y]\,]

The Law of Total Expectation — E[X] = E[E[X | Y]] — says: instead of averaging X directly, you can first average it within each level of Y, then average those conditional averages weighted by P(Y). It is a bookkeeping identity, and it reappears later in the course (reinforcement learning values a situation by exactly such a conditional average). Here you only need the arithmetic: a conditional mean is the average of one cropped, renormalized slice.

#Independence: One More Time

Recap from the probability lesson: X and Y are independent iff f(x, y) = f(x) f(y) for all x and y. Equivalently, knowing Y tells you nothing about X — the conditional f(x | y) is the same as the marginal f(x).

A weaker condition is uncorrelated: Cov(X, Y) = 0. Independence implies uncorrelated, but the converse fails. The classic counterexample: let X ~ Uniform(-1, 1) and Y = X². Then Cov(X, Y) = 0 because E[XY] = E[X³] = 0 by symmetry, but Y is a deterministic function of X — they could not be more dependent. (You met this in the random-variables lesson.)

For the special case of jointly Gaussian variables, coming up in two sections, independence and uncorrelated coincide. That is one of the many reasons the multivariate Gaussian is so beloved: in this one family, "no linear correlation" is the same thing as "no dependence at all."

#The Covariance Matrix

When you stack n random variables into a vector X = (X_1, X_2, ..., X_n), the analog of variance is the covariance matrix Σ:

\Sigma_{ij} = \text{Cov}(X_i, X_j) = E[(X_i - \mu_i)(X_j - \mu_j)]

Two facts you will use constantly:

Σ is symmetric: Σ_ij = Σ_ji. (Covariance is a symmetric operation.)
Σ is positive semi-definite (PSD): for any vector v, the quadratic form v^T Σ v ≥ 0. (Reason: v^T Σ v = Var(v^T X), and variance can never be negative.)

Symmetric + PSD means Σ has a clean eigendecomposition Σ = U Λ U^T, where U has orthonormal columns (the principal directions) and Λ is diagonal with non-negative entries (the variances along those directions). This is exactly the spectral story from the eigenvalues lesson — applied to the cloud of (X_1, ..., X_n) data points.

PCA (principal component analysis, a method for finding the directions along which data varies most) is the eigendecomposition of Σ; principal components are the eigenvectors of Σ; the variance of the data along each component is the corresponding eigenvalue. When sklearn's PCA(n_components=k) returns its top-k components, it is literally computing the top k eigenvectors of the covariance matrix. The whole field of dimensionality reduction is one identity from this lesson.

Let's see all of this in code. The playground below draws 2000 samples from a multivariate Gaussian with correlation, then recovers the covariance from the samples, eigendecomposes it, and overlays the principal axes — exactly what PCA does on real data.

A Pyodide-backed scratchpad for math lessons.

Loading visualization...

#The Multivariate Gaussian PDF

Now we can write down the star of the show. The multivariate normal distribution N(μ, Σ) over x ∈ ℝⁿ has density:

f(\mathbf{x}) = \frac{1}{(2\pi)^{n/2} \, |\Sigma|^{1/2}} \, \exp\!\left( -\tfrac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^{T} \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu}) \right)

Let's unpack each piece:

Mean vector μ ∈ ℝⁿ: where the bell is centered. Each component is the marginal mean of one coordinate.
Covariance matrix Σ: how wide the bell is in each principal direction, and how tilted it is. Bigger eigenvalues → wider in that direction. Off-diagonal entries → tilt away from the axes.
Quadratic form (x − μ)^T Σ^{-1} (x − μ): this number is small near μ and grows as you move away. The matrix Σ^{-1} (the precision matrix) reweights distance so that "1 unit" along a low-variance direction counts more than "1 unit" along a high-variance one.
Normalization (2π)^{n/2} |Σ|^{1/2}: makes the whole thing integrate to 1. The determinant |Σ| is the product of eigenvalues — the volume of the "bell" in n dimensions.

When n = 1, this collapses to the familiar 1D normal: μ becomes a number, Σ becomes σ², and the formula reduces to (2πσ²)^{−1/2} exp(−(x − μ)² / (2σ²)). The multivariate version is the same shape, generalized to a vector.

Quick check

A 2D Gaussian has Σ = [[4, 0], [0, 1]]. What does its level-set ellipse look like?

#Mahalanobis Distance: Distance That Knows About Correlations

The quadratic form in the exponent has its own name: the squared Mahalanobis distance.

d_M(\mathbf{x}, \boldsymbol{\mu}) = \sqrt{(\mathbf{x} - \boldsymbol{\mu})^{T} \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu})}

Mahalanobis distance is what classical anomaly detection uses: a point far from the cloud in this sense is an outlier. In scikit-learn, EmpiricalCovariance().fit(X).mahalanobis(X) returns the SQUARED distance of every row of X (take the square root to get d_M). It is also the geometry the multivariate Gaussian uses internally — every level set of the PDF is a level set of the Mahalanobis distance.

#Ellipses of Constant Density

If you set the exponent of the multivariate Gaussian PDF to a constant — say (x − μ)^T Σ^{-1} (x − μ) = c. You get an ellipse (in 2D) or an ellipsoid (in higher dimensions). Every point on this ellipse has the same density value, so this is exactly a contour line of the bell.

Where do the ellipse axes come from? From the eigendecomposition of Σ:

The principal axes of the ellipse are the eigenvectors of Σ.
The squared semi-axis lengths along each principal direction are c · λ_i, where λ_i is the corresponding eigenvalue.

So the bell of N(μ, Σ) is "wider along eigendirections with bigger eigenvalues." A spherical Gaussian has Σ = σ² I, all eigenvalues equal, and its level sets are circles. A "cigar-shaped" Gaussian has one big eigenvalue and the rest small — its level sets are stretched ellipses.

Intuition

Build Your Intuition

A 95% ellipse drawn around a scatter of data (a data ellipse) is a level set of a multivariate Gaussian fitted to the data, using the sample mean and sample covariance. Where does the radius 2.45 come from? In two dimensions, the fraction of a Gaussian within Mahalanobis distance r of its mean is 1 − e^(−r²/2). Set that to 0.95 and r² = −2 ln(0.05) = 5.991, so r = 2.448. (Statisticians call 5.991 the 95th percentile of a chi-squared distribution with 2 degrees of freedom. You do not need the name, only the formula.) A simulation of 200,000 draws agrees: about 95.0% of them fall inside.

Do not confuse this with a confidence region for the MEAN. The data ellipse says where individual points fall, and it stays the same size however much data you collect. A confidence ellipse for the mean says where the true center could be, and it is built from Σ/n, so its radius shrinks like 1/√n. With n = 100 it is 10 times smaller (Mahalanobis radius 2.448/10 = 0.245). In 2,000 simulated samples of 100 points, the true mean landed inside that small ellipse 94.65% of the time.

When you see two correlated quantities reported as "(x, y) = (5.2, 3.1) plus an ellipse," the ellipse is a Gaussian belief about (x, y). Engineers who track moving objects keep exactly such an ellipse around their position estimate and shrink it as measurements arrive.

What Do You Think?

Two 2D Gaussians both have variance 1 along each axis (Σ_XX = Σ_YY = 1). The first has correlation ρ = 0; the second has ρ = 0.9. The point (1.5, 1.5) is observed. Which Gaussian assigns it HIGHER density?

Watch Out

The multivariate Gaussian needs Σ to be positive-DEFINITE, not just PSD

Fix: If Σ is positive semi-definite but not strictly positive definite, then |Σ| = 0 and Σ^{-1} does not exist — the PDF formula divides by zero. What is happening geometrically: the distribution has collapsed onto a lower-dimensional subspace (the null space of Σ), so it has no density in ℝⁿ. This is called a degenerate Gaussian. It still exists as a probability law — just on a hyperplane, not on all of ℝⁿ. In code, this is what trips up numpy.linalg.cholesky(Sigma) (the matrix square root used to draw samples, explained below) with the cryptic 'matrix is not positive definite' error. The fix is to add a small ε I regularizer (often called 'jitter') or to drop the redundant linear combination of features.

#Conditioning a Multivariate Gaussian: The Magic Identity

Here is the single most useful fact about the multivariate Gaussian. Conditioning a jointly-Gaussian distribution on some of its components gives you another Gaussian, and the new mean and covariance have closed-form expressions. Numbers first, formula second.

#A worked example: weight given height

Take a person's weight X (kg) and height Y (cm), jointly Gaussian with

mean μ = (70, 170),
variances 144 and 100, so the standard deviations are 12 kg and 10 cm,
covariance 72, so the correlation is 72 / (12 × 10) = 0.6.

As a matrix, Σ = [[144, 72], [72, 100]]. Draw the ellipse with weight along the horizontal axis and height up the vertical axis. Someone is 190 cm tall, which is 20 cm above average. What do we now believe about their weight? Draw a horizontal line at height 190 and look only at the slice of the ellipse along that line. Three questions about that slice:

Where is its center (the conditional mean)? Each extra centimetre of height goes with 72 / 100 = 0.72 kg of extra weight, because 72 is the covariance and 100 is the variance of height. The person is 20 cm above average, so the mean weight moves up by 0.72 × 20 = 14.4 kg, from 70 to 84.4 kg.
How wide is it (the conditional variance)? Height explains part of the spread in weight. The explained part is 72² / 100 = 51.84 out of 144. What is left is 144 − 51.84 = 92.16, a standard deviation of 9.6 kg instead of 12. The variance shrank to 92.16 / 144 = 0.64 = 1 − 0.6² of its old size.
Does it depend on the exact height observed? The center does (a taller person means a bigger shift). The width does not: 92.16 is the same at 150 cm, 170 cm or 190 cm.

Now check this with numpy and a simulation. The simulation draws 200,000 people and keeps only those whose height is within 1 cm of 190.

pythonrunnable cell

import numpy as np

mu = np.array([70.0, 170.0])            # (weight kg, height cm)
Sigma = np.array([[144.0, 72.0],
                  [72.0, 100.0]])

# The formula, in block form: a = weight, b = height
Sab, Sbb, Saa = Sigma[:1, 1:], Sigma[1:, 1:], Sigma[:1, :1]
mean_c = mu[:1] + Sab @ np.linalg.inv(Sbb) @ (np.array([190.0]) - mu[1:])
var_c = Saa - Sab @ np.linalg.inv(Sbb) @ Sab.T
print("formula:", mean_c, var_c)

# Simulation: draw many people, keep those whose height is within 1 cm of 190
rng = np.random.default_rng(7)
X = rng.multivariate_normal(mu, Sigma, size=200_000)
near = X[np.abs(X[:, 1] - 190) < 1, 0]
print("simulation:", near.size, near.mean().round(2), near.var(ddof=1).round(2))

The formula prints a mean of 84.4 and a variance of 92.16. The simulation keeps 2,222 people and finds a mean of 84.14 and a variance of 88.95. Those are close to the hand calculation, with the gap you expect from a slice of only 2,222 draws (and a band 2 cm wide). The hand arithmetic, the matrix formula and the simulated people agree.

That calculation is the general rule. Split the variables into a block a (the one we want, here weight) and a block b (the one we observed, here height):

[x_a]      ([μ_a]    [Σ_aa  Σ_ab])
[x_b] ~ N  ([μ_b],   [Σ_ba  Σ_bb])

Then the conditional distribution x_a | x_b = x_b* is also Gaussian, with:

\boldsymbol{\mu}_{a|b} = \boldsymbol{\mu}_a + \Sigma_{ab} \Sigma_{bb}^{-1} (\mathbf{x}_b^{*} - \boldsymbol{\mu}_b) \\[6pt] \Sigma_{a|b} = \Sigma_{aa} - \Sigma_{ab} \Sigma_{bb}^{-1} \Sigma_{ba}

In the example, Σ_ab = 72 is the covariance of weight and height, Σ_bb = 100 is the variance of height, and the regression coefficient Σ_ab Σ_bb^ is the slope 0.72. The conditional mean is the regression line.

In the interactive below, press "Show conditional X | Y = y₀" to switch on the horizontal slice, then drag y₀ up and down. Watch the center of the slice slide along the line while its width stays fixed, and set ρ to zero to see the slice stop moving.

Interactive 2D multivariate Gaussian visualization.

Loading visualization...

Linear regression turns out to be a special case of this identity. The derivation is short and is in the DeepDive below, along with a preview of other places the same identity reappears.

Quick check

In a 2D Gaussian with positive correlation ρ > 0, you observe Y = +2 (above its mean). What happens to the conditional distribution of X | Y compared to the marginal of X?

#Change of Variables: Densities Under Transformations

Start with numbers. Let x be a standard normal, x ~ N(0, 1), whose density at 0 is 0.3989. Now set y = 2x. The values of y are spread twice as wide, but the total probability is still 1, so it is spread over twice the length and the peak must be half as high: the density of y at 0 is 0.3989 / 2 = 0.1995. The factor 1/2 is the "stretch" of the map, and the general rule below is that idea for any invertible map y = g(x).

f_Y(\mathbf{y}) = f_X\!\bigl(g^{-1}(\mathbf{y})\bigr) \, \bigl|\det J_{g^{-1}}(\mathbf{y})\bigr|

The Jacobian is the matrix of slopes of the map (from the Vector & Matrix Calculus lesson); in one variable it is just the derivative, which is the factor 1/2 above, and its determinant is the local volume stretch.

This formula is the foundation of normalizing flows: generative models that reshape a simple Gaussian into a complicated distribution through a chain of invertible maps, adding up the log of each stretch factor along the way. They are covered in later tracks; here you only need the one-variable idea.

#The Reparameterization Trick

Here is one final identity that is small but transformative. First, numbers. Take the weight-height covariance from the conditioning example. Its Cholesky factor (the lower-triangular matrix L with L Lᵀ = Σ, a matrix "square root" of Σ) is L = [[12, 0], [6, 8]]. Check: row 1 times itself is 144, row 1 times row 2 is 72, row 2 times itself is 36 + 64 = 100, which are the entries of Σ. To draw a person, draw two plain standard-normal numbers ε₁ and ε₂, then set weight = 70 + 12 ε₁ and height = 170 + 6 ε₁ + 8 ε₂. In general, sampling x ~ N(μ, Σ) is equivalent to:

\mathbf{x} = \boldsymbol{\mu} + L \boldsymbol{\varepsilon}, \qquad L L^{T} = \Sigma, \qquad \boldsymbol{\varepsilon} \sim N(\mathbf{0}, I)

Why is this useful? Because gradients can flow through the transformation. A VAE (variational autoencoder) is a network that squeezes each input into a small random code and learns to rebuild the input from it. Its encoder outputs (μ_θ(x), σ_θ(x)) for each input, and the latent z is sampled from N(μ_θ(x), diag(σ_θ²(x))). If you sample z directly, gradients of the reconstruction loss with respect to the encoder parameters θ are not defined — you cannot differentiate through np.random.multivariate_normal. With the reparameterization trick, you sample ε from a fixed standard normal and write z = μ_θ(x) + σ_θ(x) ⊙ ε. Now the chain rule sees a smooth path from θ to z to the loss, and end-to-end backprop just works. VAEs are trainable because of this one identity.

The same pattern, sample fixed noise and then push it through a deterministic map, is used by other generative models too. Whenever you see + ε or + noise in a deep-learning architecture, there is a good chance the reparameterization trick is at work.

#Connection to ML: Where This Lesson Bites

Anecdote

Historical Note

The multivariate Gaussian is the joint child of two giants. Carl Friedrich Gauss arrived at the (univariate) normal distribution in 1809 while developing the method of least squares for astronomical observations of the asteroid Ceres — the bell curve was the noise model that made the math work. Pierre-Simon Laplace, working independently in France, derived essentially the same density and proved the first version of the Central Limit Theorem in 1810, explaining why this particular bell shape kept showing up in nature.

The multivariate generalization came later. Francis Galton studied joint distributions of family heights in 1877 and noticed the bivariate Gaussian's elliptical contours; Karl Pearson formalized the n-dimensional version around 1900 with the modern matrix notation. The Mahalanobis distance is named for the Indian statistician P.C. Mahalanobis, who introduced it in 1936 to distinguish populations of skulls in anthropometric studies — long before anyone wrote a Kalman filter or a VAE, the same identity was being used to ask "is this skull from group A or group B?"

pythonplayground.py · Pyodide

# Build, sample, and condition a multivariate Gaussian by hand.
import numpy as np

# A jointly-Gaussian (X, Y) with negative correlation
mu = np.array([0.0, 0.0])
Sigma = np.array([[1.0, -0.7],
                [-0.7, 1.5]])

# 1) Sample 5 vectors via the reparameterization trick
L = np.linalg.cholesky(Sigma)
samples = mu + (L @ np.random.randn(2, 5)).T
print("Samples (rows):\n", samples)

# 2) Verify: the empirical covariance from many samples should match Sigma
big_sample = mu + (L @ np.random.randn(2, 50_000)).T
print("\nEmpirical Sigma:\n", np.cov(big_sample.T))
print("True Sigma:\n", Sigma)

# 3) Condition X | Y = 1.0 — the formula from the lesson
y_obs = 1.0
mu_x_given_y = mu[0] + Sigma[0, 1] / Sigma[1, 1] * (y_obs - mu[1])
sig_x_given_y = Sigma[0, 0] - Sigma[0, 1] * Sigma[0, 1] / Sigma[1, 1]
print(f"\nP(X | Y={y_obs}) is N({mu_x_given_y:.3f}, {sig_x_given_y:.3f})")
print("Notice: variance shrank because Y told us something about X")

# CHALLENGE: redo step 3 with Sigma off-diagonal = 0 (uncorrelated).
# Verify mu_x_given_y == mu[0] and sig_x_given_y == Sigma[0,0]:
# uncorrelated Gaussians have no information transfer.

Tests · Verify the conditional variance is strictly smaller than Sigma[0, 0] when correlation is nonzero. Run the conditioning formula on a 3D joint Gaussian and convince yourself the formula generalizes.

Auto-Encoding Variational Bayes

Diederik Kingma, Max Welling (2013)

#Why This Lesson Was Worth Sitting Through

Joint, marginal, and conditional are not three different objects — they are three views of the same joint distribution, related by simple operations (sum, divide, multiply). The multivariate Gaussian gives those operations a closed-form home: every marginal of a multivariate Gaussian is a Gaussian, every conditional of a multivariate Gaussian is a Gaussian, and every linear transformation of a multivariate Gaussian is a Gaussian. The family is closed under almost everything you would want to do.

That closedness is why the multivariate Gaussian is the standard "well-behaved noise" in ML, and why, when you see one Gaussian, you are often looking at a slice of a much bigger one.

The next lesson, MLE, MAP & Fisher Information, formalizes how to fit distributions like this from data, and you will see that maximum-likelihood estimation under a Gaussian noise model gives you, for free, both least-squares regression and the L2 loss you have been minimizing all along.

#Try It Yourself

Five thousand samples, one point, one conditioning formula. The exercise below uses mu = [1, 2] and Sigma = [[2.0, 1.2], [1.2, 1.0]], so X1 and X2 are positively correlated. Work through the TODOs in order, run your version, then open the solution and compare.

pythonplayground.py · Pyodide

# Sample, measure and condition a 2D Gaussian by hand.
import numpy as np

rng = np.random.default_rng(7)
mu = np.array([1.0, 2.0])
Sigma = np.array([[2.0, 1.2],
                [1.2, 1.0]])

# TODO 1: draw 5,000 points with the reparameterisation trick x = mu + L z,
#         where L = np.linalg.cholesky(Sigma) and z is standard normal.
#         Print the empirical mean and covariance (np.cov on X.T).
# TODO 2: for the point p = [3, 2], compute the Euclidean distance to mu and
#         the Mahalanobis distance sqrt(d^T Sigma^-1 d). Compare the two.
# TODO 3: apply the conditioning formula to get the mean and variance of X1
#         given X2 = 3. Then keep the samples whose X2 is within 0.05 of 3
#         and compare their mean and variance with the formula.
# TODO 4 (stretch): compare the sample variance of X1 with Sigma[0, 0] to
#         check that the marginal of X1 is N(1, 2).
# TODO 5 (stretch): compute the Mahalanobis distance of every sample from mu and
#         find the fraction within sqrt(5.991) = 2.448. It should be near 0.95.

print("Fill in the TODOs, then compare with the solution.")

Tests · Verify the empirical mean is close to [1, 2] and the empirical covariance close to Sigma, Euclidean distance 2.000 against Mahalanobis 2.673, conditional mean 2.2 and variance 0.56, a marginal variance of X1 close to 2, and about 0.95 of the samples inside the 95% ellipse.

With seed 7 the solution prints an empirical mean of [0.996, 1.986] and an empirical covariance of [[1.986, 1.206], [1.206, 1.008]], close to mu and Sigma as it should be at n = 5,000.

The point [3, 2] sits exactly 2.000 from the mean in plain Euclidean terms, but its Mahalanobis distance is 2.673. The offset [2, 0] runs against the grain of the cloud, where X1 is rarely high while X2 stays at its mean, so the point is further out than a ruler suggests.

The conditioning formula gives X1 given X2 = 3 a mean of 1 + 1.2 x (3 - 2) = 2.2 and a variance of 2 - 1.2^2 / 1 = 0.56. Only 114 of the 5,000 samples have X2 within 0.05 of 3, and those give a mean of 2.257 and a variance of 0.473. That is a small slice, so expect noise; rerun with a larger n or a narrower band to watch it tighten. The variance is well below the marginal 2.0, which is the sharpening the lesson described.

For the stretch, the sample variance of X1 comes out at 1.986 against Sigma[0, 0] = 2, and the sample mean at 0.996 against 1. Marginalising a Gaussian just reads off the matching entries of mu and Sigma.

For the last stretch, 0.952 of the 5,000 samples fall within Mahalanobis distance 2.448 of the mean, close to the 0.95 that the formula 1 - e^(-r^2/2) promises. That is a data ellipse: it describes where individual points fall.

#Quick Check

Quick Check1 / 4

Two random variables X and Y are jointly Gaussian with Cov(X, Y) = 0. What can you conclude?

Next up: MLE, MAP & Fisher Information. You will learn how to fit a distribution to data by maximum likelihood, and see that a Gaussian noise model gives you least-squares regression for free.

Multivariate Gaussian & Joint Distributions

Learning Objectives

After this lesson, you will be able to:

Reason about multiple random variables at once using joint, marginal, and conditional distributions, and know which one answers which question
Read and write the multivariate Gaussian PDF — including what each piece (mean vector, covariance matrix, determinant, quadratic form) is doing geometrically
Use the conditioning formula for a jointly-Gaussian pair to compute a conditional mean and a shrunk variance by hand, and see that it reproduces a ridge regression
Apply the change-of-variables formula and the reparameterization trick (sampling as a shift-and-stretch of plain noise), the two tools that make generative models trainable

Prerequisites

Before You Start

Review

Why

Why This Matters

How a self-driving car fuses GPS, lidar and an inertial sensor into one position estimate. Treat (true position, sensor reading) as one multivariate Gaussian, fix the reading you actually saw, and your uncertainty about the position shrinks. That operation is called conditioning, and you will do it by hand in this lesson. Engineers call the repeated version a Kalman filter
How generative models draw random numbers. Image generators and VAEs (networks that squeeze data into a few numbers and expand it back) sample from a multivariate Gaussian using a shift-and-stretch trick from this lesson, which lets gradients flow through the sampling step
Why a 95% ellipse on a scatter plot is an ellipse, not a circle. Its tilt and shape come directly from the eigenvectors and eigenvalues of the covariance matrix you will compute below
Build this → Build a "Sensor Fusion Toy" where you set the covariance of your "true position" and of a noisy sensor; the app draws both ellipses, then shows the much tighter ellipse after conditioning on a sensor reading. Watching the uncertainty shrink is the moment the math clicks

#From One Variable to Many

Earlier lessons taught you how one number wiggles. This one teaches you how many numbers wiggle together. That is the entire story. Almost every interesting ML system has dozens or millions of variables that depend on each other. To reason about them you need joint distributions, conditional distributions, and one specific multi-dimensional bell curve that shows up everywhere.

Intuition

Build Your Intuition

Later in this lesson an interactive lets you drag the covariance and watch an ellipse tilt, once the ellipse has been explained. First, the tables.

Try it! Open the Python REPL (bottom-right of the screen: click Quick Actions, then Python) and follow along.

#Joint Distributions

Here is a concrete joint PMF for "weather × mood" (totally invented numbers, but the structure is what matters):

	mood = grumpy	mood = neutral	mood = happy	row total
sunny	0.05	0.10	0.30	0.45
cloudy	0.10	0.15	0.10	0.35
rainy	0.15	0.04	0.01	0.20
column total	0.30	0.29	0.41	1.00

Every cell is a joint probability — for example P(sunny, happy) = 0.30. All cells together sum to 1, because some (weather, mood) combination happens.

P(X = x,\, Y = y) \geq 0 \quad\text{and}\quad \sum_{x}\sum_{y} P(X = x,\, Y = y) = 1

For continuous variables, the analogous statement is f(x, y) ≥ 0 and the double integral over the entire plane is 1. Probabilities of regions are volumes under the joint PDF surface.

What Do You Think?

Look at the weather × mood joint PMF above. Reading it directly, which is bigger?

#Marginal Distributions: Ignore Things You Do Not Care About

P(X = x) = \sum_{y} P(X = x,\, Y = y) \qquad f_X(x) = \int_{-\infty}^{\infty} f(x, y)\, dy

So P(sunny) = 0.45, P(cloudy) = 0.35, P(rainy) = 0.20. Those numbers add to 1 — they are a valid one-variable PMF.

Quick check

A 2D Gaussian over (X, Y) has Σ = [[1, 0.8], [0.8, 1]] (strong positive correlation). What does the MARGINAL distribution of X look like?

#Conditional Distributions: Once You Know Y

If I tell you it is rainy, what is the distribution over moods? You restrict attention to the rainy row, then renormalize so the row sums to 1.

f(x \mid y) = \frac{f(x, y)}{f(y)}

P(grumpy | rainy) = 0.15 / 0.20 = 0.75. P(happy | rainy) = 0.01 / 0.20 = 0.05. Rainy days are bleak in this fictional world.

#Conditional Expectation: A Random Variable in Its Own Right

The conditional expectation E[X | Y] is what X averages to once you know Y. Crucially, this is itself a random variable — it depends on the value of Y, which is random.

E[X \mid Y = y] = \sum_{x} x \cdot P(X = x \mid Y = y) \qquad E[X] = E[\,E[X \mid Y]\,]

#Independence: One More Time

#The Covariance Matrix

When you stack n random variables into a vector X = (X_1, X_2, ..., X_n), the analog of variance is the covariance matrix Σ:

\Sigma_{ij} = \text{Cov}(X_i, X_j) = E[(X_i - \mu_i)(X_j - \mu_j)]

Two facts you will use constantly:

Σ is symmetric: Σ_ij = Σ_ji. (Covariance is a symmetric operation.)
Σ is positive semi-definite (PSD): for any vector v, the quadratic form v^T Σ v ≥ 0. (Reason: v^T Σ v = Var(v^T X), and variance can never be negative.)

PCA (principal component analysis, a method for finding the directions along which data varies most) is the eigendecomposition of Σ; principal components are the eigenvectors of Σ; the variance of the data along each component is the corresponding eigenvalue. When sklearn's PCA(n_components=k) returns its top-k components, it is literally computing the top k eigenvectors of the covariance matrix. The whole field of dimensionality reduction is one identity from this lesson.

A Pyodide-backed scratchpad for math lessons.

Loading visualization...

#The Multivariate Gaussian PDF

Now we can write down the star of the show. The multivariate normal distribution N(μ, Σ) over x ∈ ℝⁿ has density:

f(\mathbf{x}) = \frac{1}{(2\pi)^{n/2} \, |\Sigma|^{1/2}} \, \exp\!\left( -\tfrac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^{T} \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu}) \right)

Let's unpack each piece:

Mean vector μ ∈ ℝⁿ: where the bell is centered. Each component is the marginal mean of one coordinate.
Covariance matrix Σ: how wide the bell is in each principal direction, and how tilted it is. Bigger eigenvalues → wider in that direction. Off-diagonal entries → tilt away from the axes.
Quadratic form (x − μ)^T Σ^{-1} (x − μ): this number is small near μ and grows as you move away. The matrix Σ^{-1} (the precision matrix) reweights distance so that "1 unit" along a low-variance direction counts more than "1 unit" along a high-variance one.
Normalization (2π)^{n/2} |Σ|^{1/2}: makes the whole thing integrate to 1. The determinant |Σ| is the product of eigenvalues — the volume of the "bell" in n dimensions.

Quick check

A 2D Gaussian has Σ = [[4, 0], [0, 1]]. What does its level-set ellipse look like?

#Mahalanobis Distance: Distance That Knows About Correlations

The quadratic form in the exponent has its own name: the squared Mahalanobis distance.

d_M(\mathbf{x}, \boldsymbol{\mu}) = \sqrt{(\mathbf{x} - \boldsymbol{\mu})^{T} \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu})}

#Ellipses of Constant Density

Where do the ellipse axes come from? From the eigendecomposition of Σ:

The principal axes of the ellipse are the eigenvectors of Σ.
The squared semi-axis lengths along each principal direction are c · λ_i, where λ_i is the corresponding eigenvalue.

Intuition

Build Your Intuition

What Do You Think?

Watch Out

The multivariate Gaussian needs Σ to be positive-DEFINITE, not just PSD

#Conditioning a Multivariate Gaussian: The Magic Identity

#A worked example: weight given height

Take a person's weight X (kg) and height Y (cm), jointly Gaussian with

mean μ = (70, 170),
variances 144 and 100, so the standard deviations are 12 kg and 10 cm,
covariance 72, so the correlation is 72 / (12 × 10) = 0.6.

Where is its center (the conditional mean)? Each extra centimetre of height goes with 72 / 100 = 0.72 kg of extra weight, because 72 is the covariance and 100 is the variance of height. The person is 20 cm above average, so the mean weight moves up by 0.72 × 20 = 14.4 kg, from 70 to 84.4 kg.
How wide is it (the conditional variance)? Height explains part of the spread in weight. The explained part is 72² / 100 = 51.84 out of 144. What is left is 144 − 51.84 = 92.16, a standard deviation of 9.6 kg instead of 12. The variance shrank to 92.16 / 144 = 0.64 = 1 − 0.6² of its old size.
Does it depend on the exact height observed? The center does (a taller person means a bigger shift). The width does not: 92.16 is the same at 150 cm, 170 cm or 190 cm.

Now check this with numpy and a simulation. The simulation draws 200,000 people and keeps only those whose height is within 1 cm of 190.

pythonrunnable cell

import numpy as np

mu = np.array([70.0, 170.0])            # (weight kg, height cm)
Sigma = np.array([[144.0, 72.0],
                  [72.0, 100.0]])

# The formula, in block form: a = weight, b = height
Sab, Sbb, Saa = Sigma[:1, 1:], Sigma[1:, 1:], Sigma[:1, :1]
mean_c = mu[:1] + Sab @ np.linalg.inv(Sbb) @ (np.array([190.0]) - mu[1:])
var_c = Saa - Sab @ np.linalg.inv(Sbb) @ Sab.T
print("formula:", mean_c, var_c)

# Simulation: draw many people, keep those whose height is within 1 cm of 190
rng = np.random.default_rng(7)
X = rng.multivariate_normal(mu, Sigma, size=200_000)
near = X[np.abs(X[:, 1] - 190) < 1, 0]
print("simulation:", near.size, near.mean().round(2), near.var(ddof=1).round(2))

That calculation is the general rule. Split the variables into a block a (the one we want, here weight) and a block b (the one we observed, here height):

[x_a]      ([μ_a]    [Σ_aa  Σ_ab])
[x_b] ~ N  ([μ_b],   [Σ_ba  Σ_bb])

Then the conditional distribution x_a | x_b = x_b* is also Gaussian, with:

\boldsymbol{\mu}_{a|b} = \boldsymbol{\mu}_a + \Sigma_{ab} \Sigma_{bb}^{-1} (\mathbf{x}_b^{*} - \boldsymbol{\mu}_b) \\[6pt] \Sigma_{a|b} = \Sigma_{aa} - \Sigma_{ab} \Sigma_{bb}^{-1} \Sigma_{ba}

Interactive 2D multivariate Gaussian visualization.

Loading visualization...

Linear regression turns out to be a special case of this identity. The derivation is short and is in the DeepDive below, along with a preview of other places the same identity reappears.

Quick check

In a 2D Gaussian with positive correlation ρ > 0, you observe Y = +2 (above its mean). What happens to the conditional distribution of X | Y compared to the marginal of X?

#Change of Variables: Densities Under Transformations

f_Y(\mathbf{y}) = f_X\!\bigl(g^{-1}(\mathbf{y})\bigr) \, \bigl|\det J_{g^{-1}}(\mathbf{y})\bigr|

#The Reparameterization Trick

\mathbf{x} = \boldsymbol{\mu} + L \boldsymbol{\varepsilon}, \qquad L L^{T} = \Sigma, \qquad \boldsymbol{\varepsilon} \sim N(\mathbf{0}, I)

#Connection to ML: Where This Lesson Bites

Anecdote

Historical Note

pythonplayground.py · Pyodide

# Build, sample, and condition a multivariate Gaussian by hand.
import numpy as np

# A jointly-Gaussian (X, Y) with negative correlation
mu = np.array([0.0, 0.0])
Sigma = np.array([[1.0, -0.7],
                [-0.7, 1.5]])

# 1) Sample 5 vectors via the reparameterization trick
L = np.linalg.cholesky(Sigma)
samples = mu + (L @ np.random.randn(2, 5)).T
print("Samples (rows):\n", samples)

# CHALLENGE: redo step 3 with Sigma off-diagonal = 0 (uncorrelated).
# Verify mu_x_given_y == mu[0] and sig_x_given_y == Sigma[0,0]:
# uncorrelated Gaussians have no information transfer.

Auto-Encoding Variational Bayes

Diederik Kingma, Max Welling (2013)

#Why This Lesson Was Worth Sitting Through

That closedness is why the multivariate Gaussian is the standard "well-behaved noise" in ML, and why, when you see one Gaussian, you are often looking at a slice of a much bigger one.

#Try It Yourself

pythonplayground.py · Pyodide

# Sample, measure and condition a 2D Gaussian by hand.
import numpy as np

rng = np.random.default_rng(7)
mu = np.array([1.0, 2.0])
Sigma = np.array([[2.0, 1.2],
                [1.2, 1.0]])

print("Fill in the TODOs, then compare with the solution.")

With seed 7 the solution prints an empirical mean of [0.996, 1.986] and an empirical covariance of [[1.986, 1.206], [1.206, 1.008]], close to mu and Sigma as it should be at n = 5,000.

#Quick Check

Quick Check1 / 4

Two random variables X and Y are jointly Gaussian with Cov(X, Y) = 0. What can you conclude?

Next up: MLE, MAP & Fisher Information. You will learn how to fit a distribution to data by maximum likelihood, and see that a Gaussian noise model gives you least-squares regression for free.

Multivariate Gaussian & Joint Distributions

#From One Variable to Many

#Joint Distributions

#Marginal Distributions: Ignore Things You Do Not Care About

#Conditional Distributions: Once You Know Y

#Conditional Expectation: A Random Variable in Its Own Right

#Independence: One More Time

#The Covariance Matrix

#The Multivariate Gaussian PDF

#Mahalanobis Distance: Distance That Knows About Correlations

#Ellipses of Constant Density

#Conditioning a Multivariate Gaussian: The Magic Identity

#A worked example: weight given height

Why the conditioning formula has that shape (Schur-complement sketch)

Linear regression is Gaussian conditioning (and the answer is ridge, not least squares)

Where this identity reappears later (a preview, one plain line each)

#Change of Variables: Densities Under Transformations

#The Reparameterization Trick

#Connection to ML: Where This Lesson Bites

#Why This Lesson Was Worth Sitting Through

#Try It Yourself

#Quick Check

Deep Dive

Your Reflection

Continue Learning

Discussion

Multivariate Gaussian & Joint Distributions

#From One Variable to Many

#Joint Distributions

#Marginal Distributions: Ignore Things You Do Not Care About

#Conditional Distributions: Once You Know Y

#Conditional Expectation: A Random Variable in Its Own Right

#Independence: One More Time

#The Covariance Matrix

#The Multivariate Gaussian PDF

#Mahalanobis Distance: Distance That Knows About Correlations

#Ellipses of Constant Density

#Conditioning a Multivariate Gaussian: The Magic Identity

#A worked example: weight given height

Why the conditioning formula has that shape (Schur-complement sketch)

Linear regression is Gaussian conditioning (and the answer is ridge, not least squares)

Where this identity reappears later (a preview, one plain line each)

#Change of Variables: Densities Under Transformations

#The Reparameterization Trick

#Connection to ML: Where This Lesson Bites

#Why This Lesson Was Worth Sitting Through

#Try It Yourself

#Quick Check

Deep Dive

Your Reflection

Continue Learning

Discussion