Covariance and correlation · Probability & Statistics

Independence is all-or-nothing. Most pairs of variables are somewhere in between — related, but not perfectly. Covariance and correlation measure how much, and in which direction.

Covariance

$\operatorname{Cov}(X, Y) = E\big[(X - \mu_X)(Y - \mu_Y)\big]$

Read the product inside. When X is above its mean and Y is above its mean, both factors are positive, so the product is positive. When both are below, both factors are negative — and the product is positive again. But when one is above and the other below, the product is negative.

So covariance averages agreement: positive when the two tend to move together, negative when they move oppositely, near zero when there's no consistent pattern.

The computational form

Expand and use linearity, exactly as with variance:

$\operatorname{Cov}(X, Y) = E[XY] - E[X]E[Y]$

And notice the special case:

$\operatorname{Cov}(X, X) = E[X^2] - (E[X])^2 = \operatorname{Var}(X)$

Variance is just covariance with itself. That's not a coincidence of notation — it's why the two behave so similarly.

Properties

$\operatorname{Cov}(X, Y) = \operatorname{Cov}(Y, X) \qquad \text{(symmetric)}$

$\operatorname{Cov}(aX + b, \; cY + d) = ac\operatorname{Cov}(X, Y)$

$\operatorname{Cov}(X, Y + Z) = \operatorname{Cov}(X, Y) + \operatorname{Cov}(X, Z) \qquad \text{(bilinear)}$

The shift constants b and d drop out — moving a distribution doesn't change how its parts co-vary. And now the variance rule from last lesson has a proper derivation:

$\operatorname{Var}(X + Y) = \operatorname{Cov}(X + Y, X + Y) = \operatorname{Var}(X) + \operatorname{Var}(Y) + 2\operatorname{Cov}(X, Y)$

Independent \Rightarrow covariance 0 \Rightarrow variances add. The general case carries that cross term.

The problem with covariance

Covariance has units — the product of X's units and Y's units.

Measure height in centimetres and weight in kilograms and you get a covariance in cm·kg. Switch to metres and the number shrinks by 100, though nothing about the relationship changed.

So the magnitude of a covariance is uninterpretable on its own. Is 47 a strong relationship? Unanswerable without knowing the scales. Only the sign means anything.

Correlation

Fix it by dividing out both standard deviations:

$\rho = \operatorname{Corr}(X, Y) = \frac{\operatorname{Cov}(X, Y)}{\sigma_X\,\sigma_Y}$

Now the units cancel, and the result is dimensionless and bounded:

$-1 \le \rho \le 1$

The bound follows from the Cauchy–Schwarz inequality, and its extremes are exactly the perfectly linear cases:

\rho	Meaning
+1	perfect increasing linear relationship, Y = aX + b with a > 0
0	no linear relationship
-1	perfect decreasing linear relationship, a < 0

Correlation is invariant to shifts and to positive rescaling — cm or m, °C converted linearly, it doesn't matter. A negative scaling flips the sign.

Three warnings

1. Correlation measures linear association only. The Y = X^2 example from last lesson has \rho = 0 despite Y being determined by X. A curved relationship can be perfectly strong and completely invisible to \rho. Always plot the data.

2. \rho = 0 does not mean independent (except for jointly Normal variables). It means uncorrelated, which is strictly weaker.

3. Correlation is not causation. Ice cream sales correlate with drownings; neither causes the other, and the confounder is summer. Correlation is a statement about co-movement in observed data, nothing more.

To that add a fourth, less often stated: correlation is not robust. A single outlier can drag \rho from near 0 to near 1. When data is messy, Spearman's rank correlation — Pearson's \rho applied to the ranks — is far steadier, and it also detects any monotone relationship, curved or not.

Worked example

Joint PMF of X (rows, values 0/1) and Y (columns, values 0/1):

	Y=0	Y=1
X=0	0.3	0.2
X=1	0.1	0.4

Find \operatorname{Cov}(X,Y) and \rho.

Marginals first. p_X(1) = 0.1 + 0.4 = 0.5, so E[X] = 0.5. And p_Y(1) = 0.2 + 0.4 = 0.6, so E[Y] = 0.6.

Since both are 0/1 valued, X^2 = X and Y^2 = Y, which makes the variances easy:

$\operatorname{Var}(X) = 0.5 - 0.25 = 0.25, \qquad \operatorname{Var}(Y) = 0.6 - 0.36 = 0.24$

Now E[XY]. The product XY is 1 only in the cell where both are 1:

$E[XY] = 1 \times 0.4 = 0.4$

$\operatorname{Cov}(X, Y) = 0.4 - (0.5)(0.6) = 0.4 - 0.3 = 0.1$

$\rho = \frac{0.1}{\sqrt{0.25}\sqrt{0.24}} = \frac{0.1}{0.5 \times 0.4899} \approx 0.408$

Positive and moderate. Interpretation: P(Y = 1 \mid X = 1) = 0.4/0.5 = 0.8 against P(Y = 1 \mid X = 0) = 0.2/0.5 = 0.4. Knowing X = 1 doubles the chance of Y = 1 — a real association, but far from deterministic, which is what \rho \approx 0.41 is telling you.

Doing it in Python

From the definitions, on a joint table:

import numpy as np

joint = np.array([[0.3, 0.2],     # X = 0
                  [0.1, 0.4]])    # X = 1
xs = np.array([0, 1])
ys = np.array([0, 1])

px, py = joint.sum(axis=1), joint.sum(axis=0)
EX, EY = (xs * px).sum(), (ys * py).sum()
EXY = sum(x * y * joint[i, j] for i, x in enumerate(xs) for j, y in enumerate(ys))

cov = EXY - EX * EY
sdx = ((xs**2 * px).sum() - EX**2) ** 0.5
sdy = ((ys**2 * py).sum() - EY**2) ** 0.5

print("E[X], E[Y] =", EX, EY)
print("E[XY]      =", round(EXY, 4))
print("Cov        =", round(cov, 4))
print("rho        =", round(cov / (sdx * sdy), 4))

On data, NumPy gives you both matrices at once:

import numpy as np

rng = np.random.default_rng(0)
n = 5_000
height = rng.normal(170, 10, n)                       # cm
weight = 0.5 * height - 20 + rng.normal(0, 5, n)      # kg, related to height

print("Cov matrix (cm, kg):\n", np.cov(height, weight).round(3))
print("Corr matrix:\n", np.corrcoef(height, weight).round(3))

# Change the units and watch covariance move but correlation hold still
height_m = height / 100
print("\ncov in cm.kg:", round(np.cov(height, weight)[0, 1], 4))
print("cov in m.kg :", round(np.cov(height_m, weight)[0, 1], 4), " (100x smaller)")
print("corr cm     :", round(np.corrcoef(height, weight)[0, 1], 4))
print("corr m      :", round(np.corrcoef(height_m, weight)[0, 1], 4), " (identical)")

The three failure modes, side by side — this is the most useful block in the lesson:

import numpy as np
from scipy.stats import pearsonr, spearmanr

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

x = rng.uniform(-3, 3, n)

cases = {
    "linear":            (x, 2 * x + rng.normal(0, 1, n)),
    "quadratic (y=x^2)": (x, x**2),
    "monotone but bent": (x, np.exp(x)),
}
for name, (a, b) in cases.items():
    print(f"{name:20} pearson {pearsonr(a, b)[0]:+.3f}   spearman {spearmanr(a, b)[0]:+.3f}")

# And the outlier problem
clean_x = rng.normal(0, 1, 100)
clean_y = rng.normal(0, 1, 100)
dirty_x = np.append(clean_x, 30)
dirty_y = np.append(clean_y, 30)
print(f"\n{'no outlier':20} pearson {pearsonr(clean_x, clean_y)[0]:+.3f}")
print(f"{'one outlier added':20} pearson {pearsonr(dirty_x, dirty_y)[0]:+.3f}   "
      f"spearman {spearmanr(dirty_x, dirty_y)[0]:+.3f}")

The quadratic case reports a Pearson correlation near zero for a deterministic relationship. One added point takes two unrelated variables to a correlation above 0.9. Both are exactly the traps described above.

Your turn

1. \operatorname{Var}(X) = 4, \operatorname{Var}(Y) = 9, \operatorname{Cov}(X,Y) = 3. Find \rho and \operatorname{Var}(X + Y).

2. If Y = 3X + 5, what is \rho(X, Y)? What if Y = -3X + 5?

3. Can \operatorname{Cov}(X,Y) = 10 when \sigma_X = 2 and \sigma_Y = 3?

Solutions

$\rho = \frac{3}{\sqrt{4}\sqrt{9}} = \frac{3}{6} = 0.5$

$\operatorname{Var}(X+Y) = 4 + 9 + 2(3) = 19$

Had they been independent it would be 13. The positive covariance adds 6 — when variables move together, their sum swings more widely.

2. For Y = 3X + 5: covariance is bilinear and the shift drops out, so \operatorname{Cov}(X, Y) = 3\operatorname{Var}(X), while \sigma_Y = 3\sigma_X. Therefore

$\rho = \frac{3\operatorname{Var}(X)}{\sigma_X \cdot 3\sigma_X} = \frac{3\sigma_X^2}{3\sigma_X^2} = 1$

For Y = -3X + 5: the covariance becomes -3\operatorname{Var}(X) but \sigma_Y = |-3|\sigma_X = 3\sigma_X — standard deviation takes the absolute value. So \rho = -1.

Any perfect linear relationship gives \rho = \pm 1, with the sign of the slope. The magnitude of the slope is irrelevant — which is the point of standardising.

3. No. Rearranging |\rho| \le 1:

$|\operatorname{Cov}(X,Y)| \le \sigma_X \sigma_Y = 2 \times 3 = 6$

A covariance of 10 would force \rho = 10/6 \approx 1.67, which is impossible. This is the Cauchy–Schwarz bound, and it's a genuinely useful sanity check on any covariance matrix — if an off-diagonal entry exceeds the geometric mean of the corresponding diagonals, the matrix isn't a valid covariance matrix at all.

Check yourself in code

Compute the covariance and correlation for the 2×2 joint table above, and confirm the Cauchy–Schwarz bound holds.

Print exactly this:

Cov = 0.1
rho = 0.4082
|Cov| <= sd_x * sd_y: True

Round the covariance to 4 decimal places and rho to 4.

import numpy as np

joint = np.array([[0.3, 0.2],     # X = 0
                  [0.1, 0.4]])    # X = 1
xs = np.array([0, 1])
ys = np.array([0, 1])

px, py = joint.sum(axis=1), joint.sum(axis=0)
EX, EY = (xs * px).sum(), (ys * py).sum()
EXY = sum(x * y * joint[i, j] for i, x in enumerate(xs) for j, y in enumerate(ys))

print("Cov =", round(EXY - EX * EY, 4))

# Compute the two standard deviations, then rho, then check |Cov| <= sd_x*sd_y.

import numpy as np

joint = np.array([[0.3, 0.2],
                  [0.1, 0.4]])
xs = np.array([0, 1])
ys = np.array([0, 1])

px, py = joint.sum(axis=1), joint.sum(axis=0)
EX, EY = (xs * px).sum(), (ys * py).sum()
EXY = sum(x * y * joint[i, j] for i, x in enumerate(xs) for j, y in enumerate(ys))

cov = EXY - EX * EY
print("Cov =", round(cov, 4))

sdx = ((xs**2 * px).sum() - EX**2) ** 0.5
sdy = ((ys**2 * py).sum() - EY**2) ** 0.5
print("rho =", round(cov / (sdx * sdy), 4))
print("|Cov| <= sd_x * sd_y:", abs(cov) <= sdx * sdy)

Covariance measures whether two variables move together, but its units make its size meaningless. Correlation divides those units out, landing in [-1, 1] — where \pm 1 means a perfect linear relationship and 0 means no linear relationship at all. Neither one sees curves, neither one implies causation, and neither survives a determined outlier.

Next: what happens to a distribution when you push it through a function.

17. Covariance and correlation