04 - Multivariate Distributions and Covariance

You roll two dice and write down two numbers: the larger and the smaller. They are two random variables defined on the same experiment, and they are not independent: if the maximum is 1, the minimum can only be 1.

How much do they "move together"? How do you describe a phenomenon with several quantities, such as height and weight or the returns of two stocks, where the link between the quantities matters besides the distribution of each? The answer has three ingredients: the joint law, which gives the probabilities of the pairs; the covariance, which measures the linear link; and, in the continuous case, the multivariate normal, in which everything is computed with matrix algebra.

Here you build them on the two dice and on a Gaussian example, and you test the most misunderstood idea: that zero covariance means independence. You need the discrete variables and the continuous ones, and a minimum of matrix products.

Joint law, marginals, conditionals

The joint law of two variables X and Y assigns a probability to every pair of values: in the discrete case it is the table P(X = x, Y = y). For the dice the thirty-six outcomes are equally likely and the table is built by counting in how many ways each (maximum, minimum) pair is obtained: once if the two dice coincide, twice if they differ, because the order of the rolls does not matter.

The table has everything: summing a row gives you the marginal of the maximum, summing a column that of the minimum, because the marginal of X is P(X = x) = Σ_y P(X = x, Y = y).

Dividing the joint by the marginal gives you the conditional, P(X = x | Y = y) = P(X = x, Y = y)/P(Y = y), that is the probability of X when you know the value of Y: it is the same definition as in article 01, applied to a pair.

Two variables are independent if the joint is the product of the marginals for every pair of values. In the continuous case the sums become integrals and the table a density f(x, y), with probabilities computed as volumes, but the definitions are the same.

Covariance, correlation, covariance matrix

The covariance Cov(X, Y) = E[(X − E[X])(Y − E[Y])] = E[XY] − E[X]E[Y] measures whether the two variables move away from their own mean in the same direction (positive value), in opposite directions (negative) or without a linear link (zero).

The covariance of a variable with itself is the variance. It depends on the unit of measurement, so it is normalized into the correlation coefficient ρ = Cov(X, Y)/(σ_Xσ_Y), which always lies between −1 and 1 (a consequence of the Cauchy-Schwarz inequality) and equals ±1 only if one variable is a linear function of the other.

The covariance enters the variance of a linear combination: Var(aX + bY) = a²Var(X) + b²Var(Y) + 2ab Cov(X, Y); the last term is the one that did not appear in the previous article, because for independent variables it is zero.

With more variables the covariances are collected in the covariance matrix Σ, with the variances on the diagonal and Cov(X_i, X_j) in position (i, j); it is symmetric and, since w^TΣw = Var(w^TX) ≥ 0 for every vector w of coefficients, it is positive semi-definite.

For the dice the calculation is finite and is done with exact fractions: marginals, means, E[XY], covariance and correlation, plus the mean of the maximum conditional on the minimum being 3. Then a million simulated pairs estimate the sample covariance, with a standard error computed from the centred products themselves; the assertion requires that the estimate fall within four standard errors of the exact value. A detail that simplifies the calculation: since maximum times minimum equals the product of the two dice, E[XY] is 7/2 · 7/2 = 49/4.

from fractions import Fraction as F
from itertools import product
import numpy as np

# Two dice: X = maximum, Y = minimum. Probability space: 36 equally likely pairs.
joint = {}
for d1, d2 in product(range(1, 7), repeat=2):
    chiave = (max(d1, d2), min(d1, d2))
    joint[chiave] = joint.get(chiave, 0) + F(1, 36)
assert sum(joint.values()) == 1

pX = {x: sum(p for (a, _), p in joint.items() if a == x) for x in range(1, 7)}   # marginal of the maximum
pY = {y: sum(p for (_, b), p in joint.items() if b == y) for y in range(1, 7)}   # marginal of the minimum
EX = sum(x * p for x, p in pX.items()); EY = sum(y * p for y, p in pY.items())
EXY = sum(x * y * p for (x, y), p in joint.items())
cov = EXY - EX * EY
varX = sum(x * x * p for x, p in pX.items()) - EX**2
varY = sum(y * y * p for y, p in pY.items()) - EY**2
rho = float(cov) / (float(varX) * float(varY)) ** 0.5
print("E[max] =", EX, " E[min] =", EY, " E[max*min] =", EXY)
print("Cov =", cov, "=", round(float(cov), 4), " Var(max) = Var(min) =", varX, varY == varX, " rho =", round(rho, 4))
print("P(max=1, min=6) =", joint.get((1, 6), 0), " but P(max=1) P(min=6) =", pX[1] * pY[6], " -> dependent")
cond = sum(x * p for (x, y), p in joint.items() if y == 3) / pY[3]
print("E[max | min=3] =", cond, "=", round(float(cond), 4), "versus E[max] =", round(float(EX), 4))

N = 10**6
rng = np.random.default_rng(31)
d = rng.integers(1, 7, size=(N, 2))
x, y = d.max(axis=1), d.min(axis=1)
prodotto = (x - x.mean()) * (y - y.mean())
cov_s, se = prodotto.mean(), prodotto.std() / N**0.5
print(f"Monte Carlo: Cov = {cov_s:.4f} (exact {float(cov):.4f}, standard error {se:.4f}), rho = {np.corrcoef(x, y)[0, 1]:.4f}")
assert abs(cov_s - float(cov)) < 4 * se

# Output:
# E[max] = 161/36  E[min] = 91/36  E[max*min] = 49/4
# Cov = 1225/1296 = 0.9452  Var(max) = Var(min) = 2555/1296 True  rho = 0.4795
# P(max=1, min=6) = 0  but P(max=1) P(min=6) = 1/1296  -> dependent
# E[max | min=3] = 33/7 = 4.7143 versus E[max] = 4.4722
# Monte Carlo: Cov = 0.9430 (exact 0.9452, standard error 0.0019), rho = 0.4786

The exact covariance is 1225/1296, about 0.9452, and the correlation 0.4795: maximum and minimum grow together moderately, because a high maximum makes a high minimum possible. The line on independence shows the dependence in the sharpest way: the pair (maximum 1, minimum 6) has probability 0, while the product of the marginals would be 1/1296, so they cannot be independent.

The conditional mean confirms the link: knowing that the minimum is 3, the expected maximum rises from 4.4722 to 33/7 = 4.7143. The simulation finds 0.9430 with a standard error of 0.0019, that is a deviation of 1.2 standard errors from the exact value, and a correlation of 0.4786.

Uncorrelated does not mean independent

If X and Y are independent, E[XY] = E[X]E[Y] and the covariance is zero. The converse is false, and the simplest counterexample is a variable and a function of it: X uniform on {−2, −1, 0, 1, 2} and Y = X².

Here Y is determined by X, and yet Cov(X, Y) = E[X³] − E[X]E[X²] = 0, because the values of X balance around zero and the link (a parabola) is not linear. Exercise 23 of section 6.2 of Grinstead and Snell asks exactly this: to prove that independence implies zero covariance and to find an example in which zero covariance and dependence coexist. The covariance sees only the linear component of the link.

The second case concerns the normal, because you often hear that for Gaussians uncorrelated means independent. It is true for jointly normal variables, not for two variables that are each normal. The program builds X standard normal and Y = S·X, with S = ±1 with probability 1/2 independent of X. Since the standard normal is symmetric, Y is also standard normal; but Y always has the same modulus as X. The program prints the correlation, the statistics of Y, the correlation of the squares and two conditional probabilities.

from fractions import Fraction as F
import numpy as np

# X uniform on {-2,-1,0,1,2}, Y = X^2: Y is a function of X, but the covariance is zero.
valori = [-2, -1, 0, 1, 2]
p = F(1, 5)
EX = sum(v * p for v in valori)
EY = sum(v * v * p for v in valori)
EXY = sum(v * v * v * p for v in valori)
print("Cov(X, X^2) =", EXY - EX * EY, " with E[X] =", EX, "E[X^2] =", EY)
P_Y4 = sum(p for v in valori if v * v == 4)
print("P(Y=4) =", P_Y4, " but P(Y=4 | X=2) = 1 -> dependent")

# Continuous, Gaussian case: two standard normals, uncorrelated, but not independent
N = 10**6
rng = np.random.default_rng(37)
x = rng.standard_normal(N)
segno = rng.choice([-1.0, 1.0], size=N)       # +/-1 with probability 1/2, independent of x
y = segno * x                                   # y is also standard normal (the normal is symmetric)
print(f"corr(X, Y) = {np.corrcoef(x, y)[0, 1]:+.4f} (standard error about {1 / N**0.5:.4f})")
print(f"mean and deviation of Y: {y.mean():+.4f}, {y.std():.4f}   ->  standard normal marginal")
print(f"corr(X^2, Y^2) = {np.corrcoef(x**2, y**2)[0, 1]:.4f}   (Y^2 = X^2 always)")
print(f"P(|X|<0.5) = {np.mean(np.abs(x) < .5):.4f}, P(|X|<0.5 | |Y|>2) = {np.mean(np.abs(x[np.abs(y) > 2]) < .5):.4f}")

# Output:
# Cov(X, X^2) = 0  with E[X] = 0 E[X^2] = 2
# P(Y=4) = 2/5  but P(Y=4 | X=2) = 1 -> dependent
# corr(X, Y) = -0.0001 (standard error about 0.0010)
# mean and deviation of Y: -0.0016, 0.9987   ->  standard normal marginal
# corr(X^2, Y^2) = 1.0000   (Y^2 = X^2 always)
# P(|X|<0.5) = 0.3841, P(|X|<0.5 | |Y|>2) = 0.0000

With the five values the covariance is exactly 0, but P(Y = 4) = 2/5 while P(Y = 4 | X = 2) = 1. In the Gaussian case the simulated correlation between X and Y is −0.0001, compatible with zero (the standard error is about 0.001), and Y has mean −0.0016 and deviation 0.9987, so a standard normal marginal; but the correlation between the squares is 1.0000 and the probability that |X| is below 0.5 goes from 0.3841 to 0 if you know that |Y| exceeds 2.

The pair (X, Y) is not bivariate normal: its distribution lives on two lines, not on an ellipse. Practical consequence: a correlation coefficient close to zero does not prove that two quantities are independent, and to look for non-linear dependences you need other tools, such as mutual information or kernel-based tests, which I do not cover here.

The multivariate normal

A vector X with d components is multivariate normal, X ~ N(μ, Σ), if it has a density proportional to exp(−(x − μ)^TΣ⁻¹(x − μ)/2), where μ is the vector of means and Σ the covariance matrix, here assumed positive definite (invertible). The normalization constant is (2π)^−d/2 det(Σ)^−1/2.

The quantity (x − μ)^TΣ⁻¹(x − μ) is the square of the Mahalanobis distance from μ: a distance that takes into account the scales and orientation of the cloud of points. If Z has independent standard normal components and Σ = LL^T (Cholesky factorization, with L lower triangular), then X = μ + LZ has covariance exactly Σ: it is the standard way to generate samples.

The example has two variables with standard deviations 2 and 1 and correlation 0.9 (covariance 1.8), hence a very elongated cloud. The program generates 200,000 points with Cholesky and checks the sample mean and covariance.

Then it diagonalizes Σ with numpy.linalg.eigh, suited to symmetric matrices, to find the principal directions; checks the variance of a linear combination with w^TΣw; checks that the square of the Mahalanobis distance has the chi-square distribution with 2 degrees of freedom (so 95% of the points lie below the 95% quantile of that distribution); and computes the regression of X₁ on X₂, whose coefficient is Σ₁₂/Σ₂₂ and whose residual variance is Σ₁₁ − Σ₁₂²/Σ₂₂.

import numpy as np
from scipy.stats import chi2

mu = np.array([1.0, 2.0])
Sigma = np.array([[4.0, 1.8],
                  [1.8, 1.0]])                  # sigma1 = 2, sigma2 = 1, rho = 1.8/(2*1) = 0.9
L = np.linalg.cholesky(Sigma)                   # Sigma = L L^T
rng = np.random.default_rng(41)
N = 200_000
Z = rng.standard_normal((N, 2))                  # independent standard components
X = mu + Z @ L.T                                 # X = mu + L Z, row by row

print("sample mean:", X.mean(axis=0).round(3), " (true", mu, ")")
S = np.cov(X, rowvar=False)
print("sample covariance:\n", S.round(3))
assert np.allclose(S, Sigma, atol=0.05) and np.allclose(X.mean(axis=0), mu, atol=0.02)

autoval, autovet = np.linalg.eigh(Sigma)         # eigh: symmetric matrices, eigenvalues in increasing order
print("eigenvalues:", autoval.round(4), " variance share of the 1st principal component:", round(autoval[-1] / autoval.sum(), 4))
print("direction of the 1st component:", np.abs(autovet[:, -1]).round(4))

w = np.array([1.0, -1.0])                        # combination X1 - X2
print(f"Var(X1-X2): formula w'Sigma w = {w @ Sigma @ w:.3f}, simulation {np.var(X @ w):.3f}")

# Mahalanobis distance: d^2 = (x-mu)' Sigma^-1 (x-mu) follows a chi-square with 2 degrees of freedom
D = X - mu
d2 = np.einsum('ij,ij->i', D @ np.linalg.inv(Sigma), D)
soglia = chi2.ppf(0.95, df=2)
oss = np.mean(d2 <= soglia)
tol = 3 * (0.95 * 0.05 / N) ** 0.5
print(f"95% quantile of chi-square(2) = {soglia:.4f}; share with d^2 <= threshold: {oss:.4f} (expected 0.95 +/- {tol:.4f})")
assert abs(oss - 0.95) < tol
dist_euclidea = np.sum(D**2, axis=1)
print(f"with the Euclidean distance the same threshold would keep {np.mean(dist_euclidea <= soglia):.4f} inside")

# Regression of X1 on X2: coefficient Sigma12/Sigma22 and residual variance Sigma11 - Sigma12^2/Sigma22
beta = Sigma[0, 1] / Sigma[1, 1]
res = X[:, 0] - mu[0] - beta * (X[:, 1] - mu[1])
print(f"beta = {beta:.2f}; residual variance: formula {Sigma[0,0] - Sigma[0,1]**2 / Sigma[1,1]:.3f}, simulation {res.var():.3f}; corr(residual, X2) = {np.corrcoef(res, X[:, 1])[0, 1]:+.4f}")

# Output:
# sample mean: [0.991 1.996]  (true [1. 2.] )
# sample covariance:
#  [[4.01  1.805]
#  [1.805 1.002]]
# eigenvalues: [0.1569 4.8431]  variance share of the 1st principal component: 0.9686
# direction of the 1st component: [0.9056 0.4242]
# Var(X1-X2): formula w'Sigma w = 1.400, simulation 1.401
# 95% quantile of chi-square(2) = 5.9915; share with d^2 <= threshold: 0.9498 (expected 0.95 +/- 0.0015)
# with the Euclidean distance the same threshold would keep 0.7272 inside
# beta = 1.80; residual variance: formula 0.760, simulation 0.757; corr(residual, X2) = +0.0017

The eigenvalues are 4.8431 and 0.1569, and the first principal component collects 96.86% of the total variance, with direction (0.9056, 0.4242): the cloud is almost a line, and projecting onto that axis loses 3% of the information. The variance of X₁ − X₂ is 1.400 from the formula w^TΣw and 1.401 from the simulation.

The chi-square threshold is 5.9915 and the share of points within that Mahalanobis distance is 0.9498, inside the tolerance 0.95 ± 0.0015; with the Euclidean distance, which ignores the correlation, the same threshold would keep inside only 72.72% of the points, and this is why the Mahalanobis distance is the right choice for spotting the outliers of a correlated cloud.

The regression gives coefficient 1.80 and residual variance 0.757 against 0.760 theoretical, with a correlation between residual and X₂ of +0.0017: in the Gaussian case the residual is uncorrelated (and indeed independent) from the explanatory variable, and it is the property that underpins linear regression with normal errors.

What is left out

This article does not prove the properties of the multivariate normal (Gaussian marginals and conditionals, the characterization by linear combinations), which are obtained with the moment generating function. It does not cover copulas, which separate the marginals from the dependence structure, nor mutual information, nor the estimation of Σ from few data in high dimension, where the sample covariance becomes singular or noisy.

This last point is a matter for inference. To know how much the mean of many observations depends on the link between them, you need the central limit theorem, which in the version for vectors has precisely a multivariate normal as its limit.

Frequently asked questions

Are correlation and causation the same thing?

No. Correlation measures an observed linear link between two quantities and does not say in which direction it acts, nor whether there is a common cause. This article does not go into causal inference: it only establishes what the coefficient measures.

Why must the covariance matrix be positive semi-definite?

Because the variance of any linear combination of the variables, w^TΣw, cannot be negative. If it is zero for some non-zero w, the matrix is not invertible and the vector lives in a subspace, without a density in the usual sense.

When do you use the Mahalanobis distance instead of the Euclidean one?

When the components have different scales or are correlated: it normalizes by the covariance matrix and makes distances comparable. In the program the same threshold contains 95% of the points with Mahalanobis and 72.72% with the Euclidean distance.

What you can do now

You can read a joint law, derive marginals and conditionals from it, compute covariance and correlation with exact fractions and interpret them for what they are, a linear link. You know that zero covariance is not enough to conclude independence, and you can produce the counterexample, even a Gaussian one. You can generate, diagonalize and query a multivariate normal with a few lines of linear algebra, and you know why the Mahalanobis distance is the right yardstick in a correlated cloud. In the next article you move from distributions to the behaviour of the means of many observations: the law of large numbers and the central limit theorem.

The same series also includes Exploratory Data Analysis (EDA) with Python. The full list is on the series page, Probability and Statistics.