Real data is rarely one number; it is vectors of many interrelated quantities — pixels, features, tokens. To model them we need joint distributions that capture how variables vary together, the covariance that quantifies their linear relationship, and the multivariate Gaussian, the workhorse distribution over vectors. This chapter is where probability meets the linear algebra of the companion volume.
Joint, marginal, and conditional distributions
The joint distribution of \(X\) and \(Y\) gives the probability of value pairs: \(p(x,y)=P(X=x,Y=y)\) (discrete) or a joint density \(f(x,y)\) (continuous). From it we recover: \[\textbf{marginal: } p(x)=\sum_y p(x,y),\qquad \textbf{conditional: } p(y\given x)=\frac{p(x,y)}{p(x)}.\]
The joint distribution is the complete probabilistic description of several variables at once — a table (or surface) over all combinations. Marginalizing sums/integrates out the variables you don’t care about (collapsing the table along a dimension); conditioning fixes a variable and renormalizes (taking a slice). Every probabilistic model is, in the end, a recipe for a joint distribution that we then marginalize and condition.
Covariance and correlation
The covariance of \(X\) and \(Y\) measures how they move together: \[\Cov(X,Y)=\E\big[(X-\E X)(Y-\E Y)\big]=\E[XY]-\E[X]\E[Y].\] The correlation normalizes it to \([-1,1]\):\(\Corr(X,Y)=\dfrac{\Cov(X,Y)}{\sigma_X\sigma_Y}\).
Positive covariance: they tend to rise together; negative: one rises as the other falls; zero: no linear relationship. For independent variables, \(\Cov=0\) and \(\Var[X+Y]=\Var[X]+\Var[Y]\).
Two famous traps. (1) Correlation is not causation: a strong correlation may reflect a hidden common cause, not influence. (2) Zero correlation is not independence: covariance detects only linear association. \(Y=X^2\) with \(X\) symmetric about \(0\) has \(\Cov(X,Y)=0\) yet \(Y\) is completely determined by \(X\). Independence implies zero covariance, never the reverse (except for jointly Gaussian variables, where they coincide).
The covariance matrix
For a random vector \(\vx=(X_1,\dots,X_n)\), all pairwise covariances assemble into the covariance matrix \[\boldsymbol{\Sigma}=\E\big[(\vx-\boldsymbol{\mu})(\vx-\boldsymbol{\mu})^\top\big],\qquad \Sigma_{ij}=\Cov(X_i,X_j).\] It is symmetric and positive semidefinite; its diagonal holds the variances. This single matrix is the bridge to linear algebra: its eigenvectors are the principal axes of the data cloud, and its eigendecomposition is Principal Component Analysis (PCA).
The multivariate Gaussian
The multivariate normal \(\Normal(\boldsymbol{\mu},\boldsymbol{\Sigma})\) has density \[f(\vx)=\frac{1}{(2\pi)^{n/2}|\boldsymbol{\Sigma}|^{1/2}}\exp\!\left(-\tfrac12(\vx-\boldsymbol{\mu})^\top\boldsymbol{\Sigma}^{-1}(\vx-\boldsymbol{\mu})\right).\] Its contours of constant density are ellipses (ellipsoids) centered at \(\boldsymbol{\mu}\), oriented along the eigenvectors of \(\boldsymbol{\Sigma}\), with widths set by the square roots of its eigenvalues.
The multivariate Gaussian is governed entirely by a mean vector \(\boldsymbol{\mu}\) (the center) and a covariance matrix \(\boldsymbol{\Sigma}\) (the shape). The quadratic form \((\vx-\boldsymbol{\mu})^\top\boldsymbol{\Sigma}^{-1}(\vx-\boldsymbol{\mu})\) in the exponent is the squared Mahalanobis distance — distance measured in units of standard deviations along each principal direction. Spherical \(\boldsymbol{\Sigma}=\sigma^2 I\) gives circular contours and recovers ordinary Euclidean distance; a general \(\boldsymbol{\Sigma}\) stretches and rotates them. It is the one continuous distribution over vectors that is fully analytically tractable, which is why it is everywhere in ML.
Remarkably, the multivariate Gaussian is closed under the operations we care about: marginals are Gaussian, conditionals are Gaussian, and linear transformations of a Gaussian are Gaussian. This closure is what makes Gaussian processes, Kalman filters, and linear-Gaussian models exactly solvable.
The covariance matrix and the multivariate Gaussian sit at the heart of classical and modern ML. PCA — the most-used dimensionality reduction — is the eigendecomposition of \(\boldsymbol{\Sigma}\): keep the top eigenvectors (directions of greatest variance). Whitening transforms data to \(\boldsymbol{\Sigma}=I\). The multivariate Gaussian is the latent prior \(\Normal(\vzero,I)\) in VAEs, the noise model in diffusion, the function prior in Gaussian processes, and the assumption behind Gaussian discriminant analysis. The Mahalanobis distance powers anomaly detection. And whenever a model reports a covariance over its predictions, it is reporting correlated uncertainty in exactly this language.
import numpy as np
# Sample a correlated 2-D Gaussian; recover its covariance and principal axes
mu = np.array([0.0, 0.0])
Sigma = np.array([[3.0, 1.5], [1.5, 1.0]])
X = np.random.multivariate_normal(mu, Sigma, size=200_000)
print(np.cov(X.T).round(2)) # ~ [[3, 1.5], [1.5, 1]]
vals, vecs = np.linalg.eigh(np.cov(X.T)) # eigenvectors = PCA directions
print(vals.round(2)) # variances along principal axes
# correlation != independence: Y = X^2 has zero covariance
x = np.random.randn(100_000); y = x**2
print(round(np.corrcoef(x, y)[0,1], 3)) # ~0, yet y depends fully on xChapter summary
The joint distribution describes several variables together; marginalize (sum out) and condition (slice + renormalize) to query it.
Covariance/correlation measure linear co-movement; zero correlation \(\neq\) independence, and correlation \(\neq\) causation.
The covariance matrix \(\boldsymbol{\Sigma}\) is symmetric PSD; its eigendecomposition is PCA.
The multivariate Gaussian \(\Normal(\boldsymbol\mu,\boldsymbol\Sigma)\) has ellipsoidal contours and is closed under marginals, conditionals, and linear maps.
These objects underlie PCA, whitening, Gaussian processes, VAE/diffusion priors, and Mahalanobis anomaly detection.