Advanced30 min

The Multivariate Normal Distribution

Mean vectors and covariance matrices, contours of constant density, marginals and conditionals.

Assumes you know

The Multivariate Normal Distribution

Intuition first

The multivariate normal extends the bell curve to several dimensions at once. Instead of a mean and a variance, it has a mean vector saying where the centre is and a covariance matrix saying how the dimensions vary and how they vary together.

Geometrically, the contours of constant density are ellipses (or ellipsoids). The covariance matrix determines their shape and orientation: equal variances and zero covariance give circles, unequal variances give axis-aligned ellipses, and non-zero covariance tilts them. Reading a covariance matrix as an ellipse is the single most useful mental image in multivariate statistics.

It matters out of proportion to its apparent narrowness because it is the only multivariate distribution that is genuinely tractable. Linear combinations stay normal, marginals are normal, conditionals are normal and their means are linear in the conditioning variables — which is exactly where linear regression comes from. PCA, LDA, Gaussian mixtures, Kalman filters and Gaussian processes all rest on this one distribution.

Definition

\mathbf{X} \sim \mathcal{N}_d(\boldsymbol{\mu}, \Sigma), \qquad f(\mathbf{x}) = \frac{1}{(2\pi)^{d/2}\lvert\Sigma\rvert^{1/2}}\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\mathsf{T}}\Sigma^{-1}(\mathbf{x}-\boldsymbol{\mu})\right)

Notation used in this lesson
Symbol	Meaning	Read aloud
μ	Mean vector, d components	mu
Σ	Covariance matrix, d × d, symmetric positive definite	sigma
\|Σ\|	Determinant — the generalised variance	det sigma
Σ⁻¹	Precision matrix	sigma inverse
d²	Squared Mahalanobis distance	d squared

Reading the covariance matrix

For two dimensions:

\Sigma = \begin{pmatrix} \sigma_1^2 & \rho\sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 & \sigma_2^2\end{pmatrix}

Constant-density contours. Equal variances with zero covariance give circles; unequal variances give axis-aligned ellipses; non-zero covariance tilts them, with the tilt direction set by the sign of ρ.

The eigenvectors of $\Sigma$ give the ellipse axes and the eigenvalues their squared lengths. That is precisely what PCA computes: the principal components are the eigenvectors of the covariance matrix, ordered by how much variance each direction carries.

The three properties that make it tractable

Linear combinations stay normal. For any matrix $A$ and vector $b$ :

A\mathbf{X} + b \sim \mathcal{N}\big(A\boldsymbol{\mu} + b,\; A\Sigma A^{\mathsf{T}}\big)

Marginals are normal. Just pick out the relevant entries of $\boldsymbol\mu$ and the corresponding block of $\Sigma$ — no integration required.

Conditionals are normal, with linear means. Partition $\mathbf{X} = (\mathbf{X}_1, \mathbf{X}_2)$ :

\mathbf{X}_1 \mid \mathbf{X}_2 = \mathbf{x}_2 \;\sim\; \mathcal{N}\Big(\boldsymbol{\mu}_1 + \Sigma_{12}\Sigma_{22}^{-1}(\mathbf{x}_2 - \boldsymbol{\mu}_2),\;\; \Sigma_{11} - \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\Big)

What the conditional formula is telling youAdvanced

Two features deserve attention.

The conditional mean is linear in $\mathbf{x}_2$ . The coefficient matrix $\Sigma_{12}\Sigma_{22}^{-1}$ is exactly the matrix of regression coefficients. In the bivariate case it reduces to

\E[X_1 \mid X_2 = x_2] = \mu_1 + \rho\frac{\sigma_1}{\sigma_2}(x_2 - \mu_2)

and $\rho\sigma_1/\sigma_2$ is the familiar least-squares slope. So linear regression is the conditional mean of a multivariate normal. That is why linear regression is the right model when the joint distribution is normal, and why it is only an approximation otherwise.

The conditional variance does not depend on $\mathbf{x}_2$ . The term $\Sigma_{11} - \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}$ — the Schur complement — is a constant. This is homoscedasticity, assumed by ordinary least squares, and here it is a consequence of joint normality rather than an extra assumption.

The conditional variance is also never larger than the marginal variance: conditioning on information can only reduce uncertainty, and the reduction is exactly $\rho^2$ of it in the bivariate case, since

\sigma_{1\mid 2}^2 = \sigma_1^2(1 - \rho^2)

That $1 - \rho^2$ is the residual fraction of variance, so $\rho^2$ is the fraction explained — which is $R^2$ .

Solved problem 1 · Bivariate normal: marginals, conditionals, and a probability

Heights ( $X_1$ , cm) and weights ( $X_2$ , kg) are bivariate normal with

\boldsymbol\mu = \begin{pmatrix}175 \\ 75\end{pmatrix}, \qquad \Sigma = \begin{pmatrix}64 & 60 \\ 60 & 144\end{pmatrix}

(a) Find the correlation. (b) Find the distribution of weight given height 185 cm. (c) Find $\Prob(\text{weight} > 85 \mid \text{height} = 185)$ .

Step 1 — part (a): extract the correlation

\sigma_1 = \sqrt{64} = 8, \qquad \sigma_2 = \sqrt{144} = 12

\rho = \frac{\Cov}{\sigma_1\sigma_2} = \frac{60}{8 \times 12} = \frac{60}{96} = 0.625

Step 2 — part (b): the conditional mean

Using the bivariate form with $X_2$ as the target and $X_1$ as the conditioner:

\E[X_2 \mid X_1 = x_1] = \mu_2 + \rho\frac{\sigma_2}{\sigma_1}(x_1 - \mu_1)

The slope is

\rho\frac{\sigma_2}{\sigma_1} = 0.625 \times \frac{12}{8} = 0.625 \times 1.5 = 0.9375 \text{ kg per cm}

At $x_1 = 185$ :

\E[X_2 \mid X_1 = 185] = 75 + 0.9375(185 - 175) = 75 + 9.375 = 84.375 \text{ kg}

Step 3 — the conditional variance

\sigma_{2\mid 1}^2 = \sigma_2^2(1 - \rho^2) = 144\big(1 - 0.625^2\big) = 144(1 - 0.390625) = 144 \times 0.609375 = 87.75

\sigma_{2\mid1} = \sqrt{87.75} \approx 9.3675

X_2 \mid X_1 = 185 \;\sim\; \mathcal{N}(84.375,\; 87.75)

Step 4 — part (c): standardise within the conditional

z = \frac{85 - 84.375}{9.3675} = \frac{0.625}{9.3675} \approx 0.06672

\Prob(X_2 > 85 \mid X_1 = 185) = 1 - \Phi(0.0667) \approx 1 - 0.5266 = 0.4734

Step 5 — compare with the unconditional answer

Without knowing the height:

z = \frac{85 - 75}{12} = 0.8333 \quad\Longrightarrow\quad \Prob(X_2 > 85) \approx 1 - 0.7977 = 0.2023

Knowing the person is 10 cm above average height more than doubles the probability that they exceed 85 kg, from 20% to 47%. Conditioning also cut the standard deviation from 12 to 9.37, a reduction of exactly $\sqrt{1-\rho^2} = \sqrt{0.609} = 0.781$ .

Answer

(a) $\rho = 0.625$ ; (b) $\mathcal{N}(84.375,\; 87.75)$ with SD $9.37$ ; (c) $\approx 0.473$ , against an unconditional $0.202$ .

Independence and zero correlation coincide here

For general random variables, zero correlation does not imply independence. For jointly normal variables it does:

\rho = 0 \iff X_1 \ind X_2 \qquad \text{(jointly normal only)}

Degenerate cases

If $\Sigma$ is singular the density does not exist, because $\Sigma^{-1}$ and $\lvert\Sigma\rvert^{-1/2}$ are undefined. This happens when one variable is an exact linear combination of others — the distribution then lives on a lower-dimensional subspace.

python

import numpy as np
from scipy import stats

mu = np.array([175.0, 75.0])
Sigma = np.array([[64.0, 60.0],
                  [60.0, 144.0]])

s1, s2 = np.sqrt(Sigma[0,0]), np.sqrt(Sigma[1,1])
rho = Sigma[0,1] / (s1 * s2)
print(f"rho {rho:.4f}   sigma1 {s1:.1f}  sigma2 {s2:.1f}")

# Conditional distribution of weight given height = 185.
x1 = 185.0
slope = Sigma[1,0] / Sigma[0,0]                 # = rho * s2/s1
cond_mu = mu[1] + slope * (x1 - mu[0])
cond_var = Sigma[1,1] - Sigma[1,0]**2 / Sigma[0,0]
print(f"\nslope {slope:.4f} kg/cm   cond mean {cond_mu:.4f}  cond var {cond_var:.4f}")
print(f"P(W > 85 | H = 185) {stats.norm.sf(85, cond_mu, np.sqrt(cond_var)):.4f}")
print(f"P(W > 85)           {stats.norm.sf(85, mu[1], s2):.4f}")

# Eigenstructure: the ellipse axes are the principal components.
vals, vecs = np.linalg.eigh(Sigma)
print(f"\neigenvalues {vals.round(3)}  (axis lengths squared)")
print(f"variance explained by PC1: {vals[-1]/vals.sum():.1%}")

# Verify by simulation, including the conditional slice.
rng = np.random.default_rng(0)
X = rng.multivariate_normal(mu, Sigma, 2_000_000)
band = np.abs(X[:,0] - 185) < 0.25
print(f"\nempirical cond mean {X[band,1].mean():.4f}  sd {X[band,1].std():.4f}")
print(f"theory              {cond_mu:.4f}       {np.sqrt(cond_var):.4f}")

# Mahalanobis distance is chi-squared with d degrees of freedom.
diff = X - mu
d2 = np.einsum('ij,jk,ik->i', diff, np.linalg.inv(Sigma), diff)
print(f"\nMahalanobis²: mean {d2.mean():.4f} (theory 2)  "
      f"95th pct {np.percentile(d2, 95):.4f} (chi2(2) {stats.chi2.ppf(0.95, 2):.4f})")

That final check is worth knowing: $d^2 \sim \chi^2_d$ , which is how you set a principled multivariate outlier threshold rather than eyeballing one per dimension.

Exercise 1

$\mathbf{X} \sim \mathcal{N}_2\!\left(\begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}4 & 2\\ 2 & 9\end{pmatrix}\right)$ . Find the distribution of $Y = 2X_1 - X_2$ .

Show solution

Write $Y = A\mathbf{X}$ with $A = (2, -1)$ , a $1 \times 2$ matrix. Linear combinations of a multivariate normal are normal, so $Y$ is univariate normal.

Mean:

\E[Y] = A\boldsymbol\mu = 2(0) - 1(0) = 0

Variance:

\Var(Y) = A\Sigma A^{\mathsf{T}} = \begin{pmatrix}2 & -1\end{pmatrix}\begin{pmatrix}4 & 2\\2&9\end{pmatrix}\begin{pmatrix}2\\-1\end{pmatrix}

First multiply $\Sigma A^{\mathsf T}$ :

\begin{pmatrix}4 & 2\\2&9\end{pmatrix}\begin{pmatrix}2\\-1\end{pmatrix} = \begin{pmatrix}8 - 2\\ 4 - 9\end{pmatrix} = \begin{pmatrix}6\\-5\end{pmatrix}

Then:

\begin{pmatrix}2 & -1\end{pmatrix}\begin{pmatrix}6\\-5\end{pmatrix} = 12 + 5 = 17

So $Y \sim \mathcal{N}(0, 17)$ .

Check against the scalar formula:

\Var(2X_1 - X_2) = 4\Var(X_1) + \Var(X_2) - 4\Cov(X_1,X_2) = 16 + 9 - 8 = 17 \;\checkmark

Note the positive covariance reduced the variance of the difference — the two variables move together, so subtracting cancels part of the shared fluctuation.

Exercise 2

Explain why Mahalanobis distance is preferred to Euclidean distance for detecting multivariate outliers.

Show solution

Euclidean distance treats all directions as equivalent and ignores both scale and correlation. Mahalanobis distance,

d^2 = (\mathbf{x} - \boldsymbol\mu)^{\mathsf T}\Sigma^{-1}(\mathbf{x} - \boldsymbol\mu)

corrects for both.

Scale. With height in centimetres (SD 8) and weight in kilograms (SD 12), a Euclidean distance adds centimetres to kilograms — meaningless, and dominated by whichever variable happens to have larger units. Change height to metres and the answer changes completely. Mahalanobis divides each direction by its standard deviation, making it unit-free.

Correlation. With $\rho = 0.625$ , a person 10 cm taller and 9 kg heavier than average is entirely typical, since that is precisely the pattern the correlation predicts. A person 10 cm taller and 9 kg lighter is unusual. Both points are the same Euclidean distance from the mean. Mahalanobis distance, via $\Sigma^{-1}$ , knows the ellipse is tilted and gives the second point a much larger value.

It has a calibrated threshold. For multivariate normal data, $d^2 \sim \chi^2_d$ , so the 95th percentile of $\chi^2_d$ is a principled cut-off. Euclidean distance has no such reference distribution, leaving the threshold arbitrary.

The practical caveat: $\Sigma$ must be estimated, and the sample covariance is itself distorted by outliers — the very things being detected. Robust estimators such as the minimum covariance determinant exist for this reason.

That completes the distributions in this module. Next: Statistics and Inference, which uses all of this to reason from samples back to populations.