The Gaussian is the workhorse distribution of ML — and the multivariate Gaussian is just "many Gaussians that know about each other." Stable Diffusion samples from it. VAEs squeeze data into it. Gaussian Processes ARE it. If you understand the shape of an ellipse in 2D, you already understand multivariate Gaussians — every formula in this lesson is just generalizing that ellipse to 768 dimensions.
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 derive linear regression, Gaussian-process predictions, and DDPM reverse-step posteriors from a single identity
Apply the change-of-variables formula and the reparameterization trick — the two tools that make normalizing flows and VAEs differentiable
The previous lesson taught you how one number wiggles. This one teaches you how many numbers wiggle together. That is the entire story. Every interesting ML system — every classifier, every recommender, every generative model — 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.
The interactive below is the workhorse for this lesson. Drag the off-diagonal of Σ to introduce correlation; watch the ellipse tilt and the marginals on the axes change shape. Come back to it after each section to test your intuition against your eyes.
Loading visualization...
Try it! Open the Python REPL (bottom-right of the screen: click Quick Actions, then Python) and follow along.
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)≥0andx∑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)=y∑P(X=x,Y=y)fX(x)=∫−∞∞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?
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∣y)=f(y)f(x,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. 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∣Y=y]=x∑x⋅P(X=x∣Y=y)E[X]=E[E[X∣Y]]
The Law of Total Expectation — E[X] = E[E[X | Y]] — is a workhorse identity. It 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). This decomposition is everywhere in ML:
Bellman expectation in RL: V(s) = E[r + γ V(s') | s] — the value of a state is the conditional expected reward-plus-future-value, given that state. The whole RL machinery is built on conditioning expectations on states and actions.
DDPM training target: the optimal denoiser at noise level t is E[x_0 | x_t] — the conditional mean of the clean image given the noisy one. The neural network is trained to approximate this conditional expectation.
TD-learning targets: r_t + γ V(s_{t+1}) is a one-sample estimate of E[r + γ V(s') | s_t]. Bootstrapping in TD methods is just an unbiased estimate of the same conditional expectation.
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."
When you stack n random variables into a vector X = (X_1, X_2, ..., X_n), the analog of variance is the covariance matrix Σ:
Σij=Cov(Xi,Xj)=E[(Xi−μi)(Xj−μ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 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.
Now we can write down the star of the show. The multivariate normal distribution N(μ, Σ) over x ∈ ℝⁿ has density:
f(x)=(2π)n/2∣Σ∣1/21exp(−21(x−μ)TΣ−1(x−μ))
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.
dM(x,μ)=(x−μ)TΣ−1(x−μ)
Mahalanobis distance is what classical anomaly detection uses (sklearn.covariance.MahalanobisDistance). It is also the geometry the multivariate Gaussian uses internally — every level set of the PDF is a level set of the Mahalanobis distance.
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.
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?
#Conditioning a Multivariate Gaussian: The Magic Identity
Here is the single most useful fact about the multivariate Gaussian, and the reason it shows up in so many places. Conditioning a jointly-Gaussian distribution on some of its components gives you another Gaussian — and the new mean and covariance have closed-form expressions.
Partition x = (x_a, x_b) into two blocks, with the joint distribution:
[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:
Linear regression with Gaussian noise. If y = X β + ε with ε ~ N(0, σ² I), then (β, y) is jointly Gaussian. Conditioning on the observed y gives the posterior β | y — and its mean is the closed-form least-squares solution β̂ = (X^T X)^{-1} X^T y. Linear regression IS Gaussian conditioning in disguise. Bayesian linear regression is just the same identity with a prior on β.
Gaussian Processes (Optuna's TPE under the hood, gpytorch, sklearn.gaussian_process). A GP says: "function values at any finite set of input points form a multivariate Gaussian." Predicting at a new point is conditioning the joint Gaussian over (training outputs, new output) on the training outputs — the formula above gives you the predictive mean and variance for free.
Kalman filters. The "predict / update" cycle of a Kalman filter is two operations: (1) push the prior Gaussian belief forward through linear-Gaussian dynamics — still Gaussian, easy; (2) condition that prior on the new noisy sensor measurement — exactly the formula above. Sensor fusion in self-driving cars and aerospace is one big stack of these conditioning steps.
Diffusion model reverse posteriors. In DDPM, the forward chain x_t = √(ᾱ_t) x_0 + √(1 − ᾱ_t) ε is multivariate Gaussian conditional on x_0. The Bayes-flipped reverse posterior q(x_{t-1} | x_t, x_0) is also a Gaussian, with mean and covariance that come from this identity. The neural network learns to approximate that conditional mean.
Drop back into the interactive — but this time, turn on the "conditional slice" overlay. Pick a value for Y; watch the conditional distribution of X | Y appear as a 1D Gaussian whose mean tracks the regression line and whose variance is strictly less than Σ_XX.
Loading visualization...
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
Suppose x ~ f_X and you apply an invertible transformation y = g(x). What is the density of y?
fY(y)=fX(g−1(y))detJg−1(y)
This formula is the foundation of normalizing flows — generative models that build complex distributions by chaining a sequence of invertible transformations on top of a simple base (usually a Gaussian). At each step, you apply the change-of-variables formula, accumulating log-determinants of Jacobians. The flow's training objective is exactly log p_X(x) = log p_Z(g^{-1}(x)) + log |det J_{g^{-1}}(x)|, and architectures like RealNVP, Glow, and FFJORD are designed so the Jacobian determinant has a tractable closed form.
It is also the formula that makes diffusion-model variational bounds tractable, because the forward noising chain is built from a sequence of invertible Gaussian shifts.
Here is one final identity that is small but transformative. Sampling x ~ N(μ, Σ) is equivalent to:
x=μ+Lε,LLT=Σ,ε∼N(0,I)
Why is this useful? Because gradients can flow through the transformation. In a VAE, the 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 trick — sample fixed noise, then push it through a learned deterministic map — is the engine behind score-based generative models, normalizing flows, and Gumbel-Softmax for discrete latents. Whenever you see + ε or + noise in a deep-learning architecture, there is a good chance the reparameterization trick is at work.
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.
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 universal "well-behaved noise" in ML. It is the prior over weights in Bayesian deep learning. It is the latent prior in VAEs. It is the noise schedule in diffusion. It is the building block of GP regression. It is the error covariance in Kalman filters. It is the confidence-ellipse generator behind every published 2D scatter plot. When you see one Gaussian, you are usually looking at a slice of a much bigger one.
The next lesson 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.
Two random variables X and Y are jointly Gaussian with Cov(X, Y) = 0. What can you conclude?
Next up: Statistical Inference — hypothesis testing, p-values, and confidence intervals. Then we move into MLE, MAP, and Fisher information, where every loss function in ML is rederived from probabilistic first principles.