Week 6: Multivariate Gaussians

From 1D Normal to nD#

The univariate normal (Gaussian) $X \sim \mathcal{N}(\mu, \sigma^2)$ has density:

A multivariate Gaussian (also called a multivariate normal) in $\mathbb{R}^d$ is parameterised by a mean vector $\mu \in \mathbb{R}^d$ and a covariance matrix $\Sigma \in \mathbb{R}^{d \times d}$:

Notation: $X \sim \mathcal{N}(\mu, \Sigma)$. The term $(x-\mu)^T \Sigma^{-1}(x-\mu)$ is the Mahalanobis distance squared — it generalises the squared z-score $(x-\mu)^2 / \sigma^2$ to multiple dimensions.

Level sets are ellipsoids. The set of points where the density is constant satisfies $(x-\mu)^T \Sigma^{-1}(x-\mu) = c^2$ — an ellipsoid centred at $\mu$. The shape and orientation of the ellipsoid are determined by $\Sigma$:

• Diagonal $\Sigma$ with equal entries: spherical level sets (isotropic).
• Diagonal $\Sigma$ with different entries: axis-aligned ellipses.
• Non-diagonal $\Sigma$: tilted ellipses — the off-diagonal entries encode linear dependence between dimensions.

For $d=1$, $\Sigma$ is just $\sigma^2$, and we recover the familiar bell curve.

Mean and Covariance#

The parameters $\mu$ and $\Sigma$ are exactly the mean and covariance:

• $\mathbb{E}[X] = \mu$ (an element-wise expectation: $\mu_i = \mathbb{E}[X_i]$).
• $\mathrm{Cov}(X) = \mathbb{E}[(X-\mu)(X-\mu)^T] = \Sigma$, where $\Sigma_{ij} = \mathrm{Cov}(X_i, X_j)$.

The covariance matrix is symmetric ($\Sigma = \Sigma^T$) and positive semidefinite (PSD): for any vector $v$, $v^T \Sigma v \ge 0$.

Why PSD? Consider the variance of the scalar $v^T X = \sum_i v_i X_i$. This variance must be $\ge 0$:

If $\Sigma$ were not PSD, there would exist some linear combination of the $X_i$ with negative variance — impossible. This connects to the Linear Algebra course: $\Sigma$ is a symmetric PSD matrix, so its eigenvalues are all $\ge 0$, and the condition number governs how "stretched" the ellipsoid is along each eigen-direction.

Correlation from covariance: The correlation matrix $R$ has entries $\rho_{ij} = \Sigma_{ij} / \sqrt{\Sigma_{ii} \Sigma_{jj}}$. Diagonal entries of $R$ are always 1; off-diagonals are in $[-1, 1]$.

Worked example — 2D Gaussian. Let

Then $\mathrm{Var}(X_1) = 4$, $\mathrm{Var}(X_2) = 1$, $\mathrm{Cov}(X_1, X_2) = 1$, and the correlation is $\rho = 1 / \sqrt{4 \cdot 1} = 0.5$. The ellipse is centred at $(3, 1)$, stretched more horizontally (variance 4) than vertically (variance 1), and tilted because of the positive off-diagonal covariance.

Marginals and Conditionals#

Gaussians are closed under marginalisation and conditioning — if the joint is Gaussian, every marginal and conditional is also Gaussian. This is one reason Gaussians dominate statistical modeling.

Marginals#

If $X \sim \mathcal{N}(\mu, \Sigma)$, then any subset of components is Gaussian. For a partition $X = (X_a, X_b)$:

You simply extract the corresponding sub-vector of $\mu$ and sub-matrix of $\Sigma$. No integration required. The other variables are "marginalised out" automatically.

Conditionals#

Given $X \sim \mathcal{N}(\mu, \Sigma)$ partitioned as $(X_a, X_b)$, the conditional $X_a \mid X_b = x_b$ is:

Key observations:

• The conditional mean is a linear function of the observed $x_b$ — the "correction" $\Sigma_{ab}\Sigma_{bb}^{-1}(x_b - \mu_b)$ adjusts $\mu_a$ based on how far $x_b$ is from its mean, scaled by covariance structure.
• The conditional covariance $\Sigma_{aa} - \Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}$ does not depend on the value $x_b$. Conditioning always reduces (or leaves unchanged) variance — information never increases uncertainty.
• The Schur complement $\Sigma_{aa} - \Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}$ is always PSD (smaller than $\Sigma_{aa}$ in the Loewner order).

Special case — bivariate. For $d=2$ with correlation $\rho$:

If $|\rho| = 1$, the conditional variance shrinks to zero — knowing $X_2$ determines $X_1$ exactly. If $\rho = 0$, the conditional mean and variance are the marginal values (no information gained).

Whitening and Mahalanobis#

Whitening transforms correlated data into uncorrelated, unit-variance data. Given $X \sim \mathcal{N}(\mu, \Sigma)$, the whitened variable is:

$\Sigma^{-1/2}$ is the inverse matrix square root — it exists because $\Sigma$ is PSD. After whitening, each component of $Z$ is independent standard normal.

The Mahalanobis distance between $x$ and the distribution is:

This is the Euclidean distance in the whitened space. It accounts for correlation and differing variances — a point far from the mean in a high-variance direction gets less "surprise penalty" than the same Euclidean distance in a low-variance direction. This is used in anomaly detection, clustering, and as a building block in Kalman filters.

Why Gaussians Dominate ML Modeling#

• Central Limit Theorem (Week 7): sums of many independent random variables become approximately Gaussian — errors, noise, and aggregate quantities are often near-normal.
• Closed under linear operations: if $X \sim \mathcal{N}(\mu, \Sigma)$ and $Y = AX + b$, then $Y \sim \mathcal{N}(A\mu + b, A\Sigma A^T)$. Linear transformations preserve Gaussianity.
• Conjugacy (Week 9): a Gaussian prior with a Gaussian likelihood yields a Gaussian posterior — Bayesian updates stay in the same family.
• Maximum entropy: among all distributions with a given mean and covariance, the Gaussian has the largest entropy — it is the "least informative" assumption consistent with those moments.
• Limitations: real data is often heavier-tailed, multimodal, or bounded. But Gaussians remain the default building block because of their mathematical tractability.

Knowledge check#

Browser lab#

Sample a 2D Gaussian, plot scatter with ellipse level curves, and show a conditional slice histogram.

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

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

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]
cond_val = 3.0
near = samples[np.abs(samples[:, 0] - cond_val) < 0.2]
ax.hist(near[:, 1], bins=30, density=True, alpha=0.7, label=f"X₂ | X₁≈ {cond_val}")

# True conditional: X₂ | X₁ = x₁
rho = Sigma[0, 1] / np.sqrt(Sigma[0, 0] * Sigma[1, 1])
cond_mean = mu[1] + rho * np.sqrt(Sigma[1, 1] / Sigma[0, 0]) * (cond_val - mu[0])
cond_var = Sigma[1, 1] * (1 - rho**2)
x2_grid = np.linspace(mu[1] - 3, mu[1] + 3, 200)
cond_pdf = (1 / np.sqrt(2 * np.pi * cond_var)) * np.exp(-0.5 * (x2_grid - cond_mean)**2 / cond_var)
ax.plot(x2_grid, cond_pdf, "tab:red", linewidth=2, label="True conditional")
ax.set_title(f"Conditional X₂ | X₁ = {cond_val}")
ax.set_xlabel("X₂"); ax.legend()

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(f"Conditional X₂ | X₁=3: mean={cond_mean:.3f}, var={cond_var:.3f}")