Week 6: Multivariate Gaussians

A sensor reads $x_1 = 5$. A second quantity $X_2$ moves with the reading, and the pair is jointly Gaussian with mean $\mu = (3, 1)$, standard deviations $\sigma_1 = 2$ and $\sigma_2 = 1$, and correlation $\rho$. What should you now believe about $X_2$? The surface below is the joint density. The amber curve is a vertical slice at the observed $x_1$.

Slide $x_1$. The slice's peak follows the line on the floor, and its height falls as $x_1$ leaves the mean. Its width stays the same. Switch $\rho$ to 0.9 and the surface narrows into a ridge, and every slice gets thinner. At $\rho = 0$ the slice is centred at $x_2 = 1$ wherever you put it.

Mean and covariance matrix#

A multivariate Gaussian in $\mathbb{R}^d$ has a mean vector $\mu$ and a covariance matrix $\Sigma$:

The exponent holds the squared Mahalanobis distance $(x-\mu)^T \Sigma^{-1}(x-\mu)$, the multivariate form of the squared z-score $(x-\mu)^2/\sigma^2$. Its level sets are ellipses centred at $\mu$. The scene's default is the worked example

with $\Sigma_{12} = \rho\,\sigma_1\sigma_2$, so $\rho = 1/\sqrt{4 \cdot 1} = 0.5$. The ellipses stretch along $x_1$ and tilt because the covariance is positive.

The slice is the conditional#

Fix $x_1$ and read $f(x_1, x_2)$ as a function of $x_2$. That is the slice. The conditional density is

and the marginal $f_1(x_1)$ is one number for a fixed $x_1$. So the slice is the conditional scaled by $f_1(x_1)$: the same centre and width, a different height. For two variables,

At $\rho = 0.5$ the mean is $1 + 0.25\,(x_1 - 3)$, which is 1.5 at $x_1 = 5$ and 0.25 at $x_1 = 0$. The variance is $1 \cdot (1 - 0.25) = 0.75$, a standard deviation of 0.87, at every $x_1$. At $\rho = 0.9$ the slope is 0.45 and the standard deviation $\sqrt{0.19} = 0.44$. At $|\rho| = 1$ the variance is 0.

The common mistake. A surprising reading, far from the mean, feels as if it should leave you less certain. It does not. At $x_1 = 0$ the marginal density is 0.065, against 0.121 at $x_1 = 5$, so the slice is about half as tall. Its width is the same. The reading moves only the centre. Not every joint does this: a multivariate t widens its conditional after a surprising reading.

The general partitioned form, for blocks $X = (X_a, X_b)$, is

The mean is linear in $x_b$. The covariance is the Schur complement of $\Sigma_{bb}$ in $\Sigma$. It and $\Sigma_{aa}$ minus it are both PSD, and it does not depend on $x_b$. The Kalman filter's update step uses the same formula.

Marginals#

The marginal of a block is read off the parameters, with no integration: $X_a \sim \mathcal{N}(\mu_a, \Sigma_{aa})$. In the example $X_1 \sim \mathcal{N}(3, 4)$, and its density at $x_1 = 5$ is 0.121. This is the factor that sets the slice's height.

Why $\Sigma$ must be PSD#

The covariance $\Sigma_{ij} = \mathrm{Cov}(X_i, X_j)$ is symmetric. For any $v$, the scalar $v^T X$ has variance

so $\Sigma$ is positive semidefinite. The density needs $\Sigma^{-1}$, so it needs $\Sigma$ positive definite. The eigenvectors of $\Sigma$ are the axes of the ellipses, and the square roots of its eigenvalues set their lengths. For the example these eigenvalues are about 4.30 and 0.70.

Whitening and Mahalanobis distance#

Whitening maps correlated data to independent unit-variance components:

The Mahalanobis distance $D_M(x) = \sqrt{(x-\mu)^T \Sigma^{-1}(x-\mu)}$ is the Euclidean distance after whitening. A step along a high-variance direction costs less than the same step along a low-variance one.

Why Gaussians dominate#

Sums of many independent effects tend to Gaussians (the CLT, Week 7). Linear maps keep them Gaussian: $AX + b \sim \mathcal{N}(A\mu + b, A\Sigma A^T)$. A Gaussian prior and likelihood give a Gaussian posterior (Week 9). Real data is often heavier-tailed or multimodal, but these closed forms keep Gaussians the default.

Knowledge Check#

Browser lab: the slice from samples#

Sample the worked example, draw the level sets, and compare histograms of samples near $x_1 = 0$ and $x_1 = 5$ with the conditionals from the slice.

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import eigh

rng = np.random.default_rng(42)
n = 100000

mu = np.array([3.0, 1.0])
Sigma = np.array([[4.0, 1.0],
                   [1.0, 1.0]])

samples = rng.multivariate_normal(mu, Sigma, size=n)

def plot_ellipse(ax, mu, Sigma, n_std=2, **kwargs):
    vals, vecs = eigh(Sigma)
    theta = np.linspace(0, 2*np.pi, 200)
    circle = np.array([np.cos(theta), np.sin(theta)])
    ellipse = (vecs @ np.diag(np.sqrt(vals) * n_std)) @ circle + mu[:, None]
    ax.plot(ellipse[0], ellipse[1], **kwargs)

fig, axes = plt.subplots(1, 3, figsize=(14, 4.5))

ax = axes[0]
ax.scatter(samples[:800, 0], samples[:800, 1], alpha=0.3, s=5)
plot_ellipse(ax, mu, Sigma, n_std=1, color="tab:red", linewidth=2, label="1σ")
plot_ellipse(ax, mu, Sigma, n_std=2, color="tab:red", linewidth=1.5, linestyle="--", label="2σ")
ax.scatter(*mu, color="black", marker="x", s=80, label="μ")
ax.set_aspect("equal")
ax.set_xlabel("X₁"); ax.set_ylabel("X₂")
ax.set_title("Scatter + level sets")
ax.legend()

ax = axes[1]
rho = Sigma[0, 1] / np.sqrt(Sigma[0, 0] * Sigma[1, 1])
cond_var = Sigma[1, 1] * (1 - rho**2)
x2_grid = np.linspace(mu[1] - 4, mu[1] + 4, 200)
for cond_val, color in [(0.0, "tab:orange"), (5.0, "tab:red")]:
    near = samples[np.abs(samples[:, 0] - cond_val) < 0.2]
    # Conditional mean moves along a line in x1; the variance does not move
    cond_mean = mu[1] + rho * np.sqrt(Sigma[1, 1] / Sigma[0, 0]) * (cond_val - mu[0])
    cond_pdf = np.exp(-0.5 * (x2_grid - cond_mean)**2 / cond_var) / np.sqrt(2 * np.pi * cond_var)
    ax.hist(near[:, 1], bins=30, density=True, alpha=0.4, color=color,
            label=f"X₂ | X₁≈{cond_val} ({len(near)} samples)")
    ax.plot(x2_grid, cond_pdf, color=color, linewidth=2)
    print(f"x₁ = {cond_val}: {len(near):4d} samples, mean {near[:, 1].mean():.3f} "
          f"(exact {cond_mean:.3f}), sd {near[:, 1].std():.3f} (exact {np.sqrt(cond_var):.3f})")
ax.set_title("Two slices, one width")
ax.set_xlabel("X₂"); ax.legend(fontsize=8)

ax = axes[2]
x1_marginal = np.linspace(mu[0] - 4, mu[0] + 4, 200)
marg_pdf = (1 / np.sqrt(2 * np.pi * Sigma[0, 0])) * np.exp(-0.5 * (x1_marginal - mu[0])**2 / Sigma[0, 0])
ax.hist(samples[:, 0], bins=50, density=True, alpha=0.5, label="Empirical marginal")
ax.plot(x1_marginal, marg_pdf, "tab:blue", linewidth=2, label="True marginal")
ax.set_title("Marginal X₁")
ax.set_xlabel("X₁"); ax.legend()

plt.tight_layout()
plt.show()

print(f"Sample mean: {samples.mean(axis=0)}")
print(f"Sample cov:\n{np.cov(samples.T)}")
print("Notice: both slices have sd near the exact 0.87, the scene's width at ρ = 0.5. "
      "x₁ = 0 has about half as many samples (a lower slice), not a wider spread.")

Key takeaways#

• A multivariate Gaussian is a mean and a covariance. Its level sets are ellipses whose axes are the eigenvectors of $\Sigma$.
• A slice of the joint at $x_1$ has the shape of the conditional. Its centre moves along a line in $x_1$, and its width is set by $\Sigma$ alone.
• The general conditional covariance is the Schur complement. It is PSD and never larger than the marginal covariance.
• Marginals are read off $\mu$ and $\Sigma$. $\Sigma$ is PSD because $v^T \Sigma v$ is a variance.
• Whitening turns correlated components into independent standard normals, and the Mahalanobis distance is Euclidean distance after whitening.

Further reading#

• Christopher M. Bishop, Pattern Recognition and Machine Learning (Springer, 2006): §2.3 is the standard treatment of the multivariate Gaussian, marginals, and conditionals.
• Kevin P. Murphy, Machine Learning: A Probabilistic Perspective (MIT Press, 2012): Chapter 4 covers Gaussian models, covariance, and linear-Gaussian systems.