The multivariate Normal · Probability & Statistics

The Normal distribution generalises to several dimensions at once, and the result is the single most important joint distribution in statistics. It's what regression assumes, what the multivariate Central Limit Theorem produces, and what almost every classical method has hiding in its derivation.

From one dimension to many

Recall the univariate density:

$f(x) = \frac{1}{\sigma\sqrt{2\pi}}\exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$

Rewrite the exponent as -\frac{1}{2}(x - \mu)(\sigma^2)^{-1}(x - \mu) — a squared distance from the mean, scaled by the inverse variance. That form generalises directly.

For a k-dimensional vector \mathbf{X} = (X_1, \dots, X_k)^\top:

$f(\mathbf{x}) = \frac{1}{(2\pi)^{k/2}|\boldsymbol\Sigma|^{1/2}}\exp\!\left(-\frac{1}{2}(\mathbf{x} - \boldsymbol\mu)^\top\boldsymbol\Sigma^{-1}(\mathbf{x} - \boldsymbol\mu)\right)$

Written \mathbf{X} \sim N_k(\boldsymbol\mu, \boldsymbol\Sigma), with

\boldsymbol\mu the mean vector, \mu_i = E[X_i]
\boldsymbol\Sigma the covariance matrix, \Sigma_{ij} = \operatorname{Cov}(X_i, X_j)

The diagonal of \boldsymbol\Sigma holds the variances (since \operatorname{Cov}(X_i,X_i) = \operatorname{Var}(X_i)), and the off-diagonal entries hold the covariances. \boldsymbol\Sigma is symmetric and must be positive semi-definite — the matrix version of "a variance can't be negative".

The Mahalanobis distance

The quadratic form in the exponent,

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

is the Mahalanobis distance — a distance that accounts for the spread and correlation of the data. Density depends on \mathbf{x} only through it, so surfaces of constant density are surfaces of constant d: ellipses in 2-D, ellipsoids above.

Why not ordinary Euclidean distance? Because "far" should mean "surprising", and that depends on the spread. A point 3 cm from the mean height is unremarkable; 3 cm from the mean pupil diameter is extraordinary. Dividing by \boldsymbol\Sigma makes the distance dimensionless and comparable — it's the multivariate version of standardising to a z-score.

The ellipse orientation comes from the correlation: zero correlation gives axis- aligned ellipses; positive correlation tilts them along the diagonal; and as \rho \to \pm 1 they collapse toward a line.

Usefully, d^2 \sim \chi^2_k when \mathbf{X} really is multivariate Normal — which gives you both a way to draw confidence ellipses and a way to detect outliers.

The properties that make it special

1. Every marginal is Normal. Drop any variables and what remains is multivariate Normal, with the corresponding sub-vector and sub-matrix. No integration required — just cross out rows and columns.

2. Every conditional is Normal. For a partition into \mathbf{X}_1 and \mathbf{X}_2:

$\mathbf{X}_1 \mid \mathbf{X}_2 = \mathbf{x}_2 \;\sim\; N\!\left(\boldsymbol\mu_1 + \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1}(\mathbf{x}_2 - \boldsymbol\mu_2), \;\; \boldsymbol\Sigma_{11} - \boldsymbol\Sigma_{12}\boldsymbol\Sigma_{22}^{-1}\boldsymbol\Sigma_{21}\right)$

Look closely at that mean: it is linear in \mathbf{x}_2. This is where linear regression comes from (§6) — the best predictor of one Normal variable given others is automatically a linear function.

And the conditional variance doesn't depend on \mathbf{x}_2 at all. It's always smaller than the marginal variance \boldsymbol\Sigma_{11} — conditioning on information can only reduce uncertainty. That's the assumption called homoscedasticity in §6, and here it's a theorem rather than a hope.

3. Any linear combination is Normal. For a matrix \mathbf{A} and vector \mathbf{b}:

$\mathbf{A}\mathbf{X} + \mathbf{b} \sim N(\mathbf{A}\boldsymbol\mu + \mathbf{b},\; \mathbf{A}\boldsymbol\Sigma\mathbf{A}^\top)$

The family is closed under every linear operation — sums, differences, projections, rotations.

4. Uncorrelated \Rightarrow independent. This is the exception flagged back in lesson 2 of this module. If \boldsymbol\Sigma is diagonal, then \boldsymbol\Sigma^{-1} is diagonal too, the quadratic form splits into a sum of separate terms, and the exponential of a sum is a product of exponentials — so the joint density factors.

$\rho = 0 \iff \text{independent} \qquad \textbf{(jointly Normal only)}$

The trap. This requires joint normality, not just Normal marginals. Let X \sim N(0,1) and let S = \pm 1 by a fair coin, independent of X; set Y = SX. Then Y is also standard Normal and \operatorname{Cov}(X,Y) = 0 — yet |Y| = |X| always, so they are wildly dependent. Two Normal marginals do not make a multivariate Normal.

The bivariate case in full

For k = 2 it's worth seeing explicitly:

$\boldsymbol\mu = \begin{pmatrix}\mu_1\\\mu_2\end{pmatrix}, \qquad \boldsymbol\Sigma = \begin{pmatrix}\sigma_1^2 & \rho\sigma_1\sigma_2\\ \rho\sigma_1\sigma_2 & \sigma_2^2\end{pmatrix}$

The conditional distribution simplifies to something you can hold in your head:

$Y \mid X = x \;\sim\; N\!\left(\mu_2 + \rho\frac{\sigma_2}{\sigma_1}(x - \mu_1), \;\; \sigma_2^2(1 - \rho^2)\right)$

Two things fall out. The conditional mean is the regression line, with slope \rho\sigma_2/\sigma_1. And the conditional variance is reduced by the factor (1 - \rho^2) — which is exactly the R^2 of §6, appearing here as the fraction of variance that knowing X removes.

Worked example

Heights and weights are bivariate Normal with \mu_H = 170 cm, \sigma_H = 10, \mu_W = 70 kg, \sigma_W = 12, \rho = 0.6. Someone is 185 cm tall. What's the distribution of their weight?

Apply the conditional formula. The slope:

$\rho\frac{\sigma_W}{\sigma_H} = 0.6 \times \frac{12}{10} = 0.72 \text{ kg per cm}$

The conditional mean:

$E[W \mid H = 185] = 70 + 0.72(185 - 170) = 70 + 10.8 = 80.8 \text{ kg}$

The conditional variance:

$\operatorname{Var}(W \mid H) = 144(1 - 0.36) = 144 \times 0.64 = 92.16, \qquad \sigma = 9.6$

So W \mid H = 185 \sim N(80.8,\; 9.6^2).

Two observations. Knowing the height shrank the standard deviation from 12 to 9.6 — a 20% reduction, and 1 - \sqrt{1-\rho^2} = 0.2 says so in advance.

And note the prediction is not proportional: a person 15 cm above average height is predicted to be only 10.8 kg above average weight, not 15 \times \frac{12}{10} = 18. The correlation of 0.6 pulls the prediction back toward the mean. That's regression to the mean — the phenomenon that gave regression its name, visible here as the factor \rho in the slope.

Doing it in Python

scipy.stats.multivariate_normal handles densities and sampling:

import numpy as np
from scipy.stats import multivariate_normal

mu = np.array([170.0, 70.0])
sd = np.array([10.0, 12.0])
rho = 0.6

cov = np.array([[sd[0]**2,          rho*sd[0]*sd[1]],
                [rho*sd[0]*sd[1],   sd[1]**2       ]])
print("covariance matrix:\n", cov)

mvn = multivariate_normal(mean=mu, cov=cov)
print("\ndensity at the mean:", round(mvn.pdf(mu), 8))
print("density at (185, 81):", round(mvn.pdf([185, 81]), 8))

The conditional formula, checked against simulation:

import numpy as np

mu_h, mu_w, sd_h, sd_w, rho = 170.0, 70.0, 10.0, 12.0, 0.6

slope = rho * sd_w / sd_h
cond_mean = mu_w + slope * (185 - mu_h)
cond_sd = sd_w * np.sqrt(1 - rho**2)
print(f"theory:    mean {cond_mean:.2f}  sd {cond_sd:.2f}")

rng = np.random.default_rng(0)
cov = [[sd_h**2, rho*sd_h*sd_w], [rho*sd_h*sd_w, sd_w**2]]
data = rng.multivariate_normal([mu_h, mu_w], cov, size=2_000_000)

near = np.abs(data[:, 0] - 185) < 0.25          # people about 185 cm tall
print(f"simulated: mean {near.sum() and data[near, 1].mean():.2f}  "
      f"sd {data[near, 1].std():.2f}   (n = {near.sum():,})")

Uncorrelated implies independent — but only under joint normality. This block builds the counterexample:

import numpy as np

rng = np.random.default_rng(1)
n = 500_000

X = rng.standard_normal(n)
S = rng.choice([-1.0, 1.0], n)     # fair sign, independent of X
Y = S * X                          # Y is also standard Normal

print("X is normal-ish: mean", round(X.mean(), 3), " sd", round(X.std(), 3))
print("Y is normal-ish: mean", round(Y.mean(), 3), " sd", round(Y.std(), 3))
print("correlation    :", round(np.corrcoef(X, Y)[0, 1], 4), "-> uncorrelated")
print()
print("But |Y| == |X| always:", bool(np.allclose(np.abs(X), np.abs(Y))))
print("-> Normal marginals + zero correlation, yet totally dependent.")
print("   The pair is NOT jointly Normal, so property 4 does not apply.")

And the Mahalanobis distance as an outlier detector, using d^2 \sim \chi^2_k:

import numpy as np
from scipy.stats import chi2

rng = np.random.default_rng(2)
cov = np.array([[100.0, 72.0], [72.0, 144.0]])
mu = np.array([170.0, 70.0])
data = rng.multivariate_normal(mu, cov, size=10_000)

inv = np.linalg.inv(cov)
d = data - mu
d2 = np.einsum("ij,jk,ik->i", d, inv, d)        # Mahalanobis squared

print("mean d^2:", round(d2.mean(), 3), " (should be k = 2)")
for level in (0.90, 0.95, 0.99):
    cut = chi2.ppf(level, df=2)
    print(f"  inside the {level:.0%} ellipse: {(d2 <= cut).mean():.4f}")

Your turn

1. $\mathbf{X} \sim N_2!\left(\begin{pmatrix}0\0\end{pmatrix}, \begin{pmatrix}1 & 0\0 & 1\end{pmatrix}\right)$. Are X_1, X_2 independent? What shape are the contours?

2. For the height/weight model, what's the distribution of H + W?

3. If \rho = 0.9, by what factor does knowing X reduce the variance of Y?

Solutions

1. The covariance matrix is diagonal, so the variables are uncorrelated — and since they're jointly Normal, uncorrelated does imply independent here.

The contours are circles. With \boldsymbol\Sigma = \mathbf{I} the Mahalanobis distance reduces to ordinary Euclidean distance, d^2 = x_1^2 + x_2^2, so constant-density curves are x_1^2 + x_2^2 = c.

Equal variances give circles; unequal variances with zero correlation give axis-aligned ellipses; nonzero correlation tilts them.

2. Use property 3 with \mathbf{A} = (1\;\;1):

$E[H + W] = 170 + 70 = 240$

$\operatorname{Var}(H + W) = \sigma_H^2 + \sigma_W^2 + 2\rho\sigma_H\sigma_W = 100 + 144 + 2(0.6)(10)(12) = 388$

So H + W \sim N(240,\; 388), with \sigma = \sqrt{388} \approx 19.7.

Note the covariance term contributes 144 — more than a third of the total. Had we wrongly assumed independence we'd have got 244 and understated the spread by 21%.

3. The conditional variance is \sigma_Y^2(1 - \rho^2), so with \rho = 0.9:

$1 - 0.81 = 0.19$

Knowing X removes 81% of the variance, leaving 19%. In standard-deviation terms the reduction is milder: \sqrt{0.19} \approx 0.436, so the SD falls to about 44% of its original size.

That \rho^2 is precisely the R^2 of §6. And the gap between "81% of the variance" and "56% of the standard deviation" is why R^2 tends to sound more impressive than the predictions actually are.

Check yourself in code

Build the height/weight covariance matrix, then compute the conditional distribution of weight given a height of 185 cm.

Print exactly this:

slope 0.72
conditional mean 80.8
conditional sd 9.6
variance reduced by 36.0 %

Round the slope to 2 decimal places, the mean to 1, the sd to 1, and the percentage to 1.

import numpy as np

mu_h, mu_w = 170.0, 70.0
sd_h, sd_w = 10.0, 12.0
rho = 0.6

slope = rho * sd_w / sd_h
print("slope", round(slope, 2))

# Print the conditional mean at h = 185, the conditional sd,
# and the percentage of variance removed by conditioning.

import numpy as np

mu_h, mu_w = 170.0, 70.0
sd_h, sd_w = 10.0, 12.0
rho = 0.6

slope = rho * sd_w / sd_h
print("slope", round(slope, 2))
print("conditional mean", round(mu_w + slope * (185 - mu_h), 1))
print("conditional sd", round(sd_w * np.sqrt(1 - rho**2), 1))
print("variance reduced by", round(100 * rho**2, 1), "%")

The multivariate Normal is defined by a mean vector and a covariance matrix, with elliptical contours set by the Mahalanobis distance. Its marginals and conditionals are all Normal, its conditional means are linear — which is where regression comes from — and it is the one family where zero correlation really does mean independence.

Next: copulas, which separate the marginals of a joint distribution from the dependence structure between them.

20. The multivariate Normal