Expectation and variance

Intuition

Almost every page on this site quietly assumes these two numbers. A loss is an expectation. The bias-variance decomposition is an identity about them. A confidence interval is a standard deviation in disguise.

The mean says where a quantity tends to sit, and it "tells us essentially where the random variable tends to be located" [1, §22.6.1.4]. That is not enough on its own, and d2l puts the reason better than most textbooks: "making a profit of $10 ± $1 per sale is very different from making $10 ± $15 per sale despite having the same average value. The second one has a much larger degree of fluctuation, and thus represents a much larger risk" [1, §22.6.1.4].

So you need both. This page is about the three places where knowing both is still not enough: variance does not simply add, the standard deviation's guarantee is much weaker than people assume, and for some distributions neither number exists at all.

Formal definition

The mean of a discrete variable is "the weighted average: sum the values times the probability that the random variable takes on that value" [1, §22.6.1.4]. It has two properties worth memorising, both from [1, §22.6.1.4]:

• μ_aX+b = a μ_X + b
• μ_X+Y = μ_X + μ_Y, for any X and Y

That second one has no conditions. Expectations always add.

The variance is "a quantitative measure of how far a random variable deviates from the mean" [1, §22.6.1.5]. Its properties look similar and are not [1, §22.6.1.5]:

• Var(X) ≥ 0, with equality only if X is constant
• Var(aX + b) = a² Var(X) — the square is why
• Var(X + Y) = Var(X) + Var(Y) only if X and Y are independent

The standard deviation exists because the variance is in the wrong units. Measuring a star rating, "we want to look at (X − μ_X)², which is in units of squared stars. This means that the variance itself is not comparable to the original measurements" [1, §22.6.1.5]. Taking the square root puts it back in stars.

Chebyshev's inequality is what makes the standard deviation interpretable for any distribution at all: P(X ∉ [μ − ασ, μ + ασ]) ≤ 1/α². Stated verbally for α = 10, "99% of the samples from any random variable fall within 10 standard deviations of the mean" [1, §22.6.1.6]. No assumption of normality, no assumption of anything.

Think Stats reaches the variance from the sample side and notes a trap that matters in code: "the name 'sample variance' can refer to either S² or S²_n−1, and the symbol S² is used for either or both" [2, §8.2]. NumPy's var divides by n; pandas divides by n−1.

How it works

The continuous case changes nothing. Split the line into intervals of length ε, and the sum becomes an integral: μ_X = ∫ x p(x) dx [1, §22.6.1.7]. Everything above survives.

Except when the integral does not converge. The Cauchy distribution has density 1/(1+x²), which "looks innocent, and indeed consulting a table of integrals will show it has area one under it, and thus it defines a continuous random variable" [1, §22.6.1.7]. Its variance integral is infinite, and worse, the mean integral "is in essence log(∞) = ∞, so there is no well-defined average value either" [1, §22.6.1.7].

That is what the figure is showing. Not a slow convergence, not a big variance — no mean for the sample mean to converge to.

Where it's used

• Every loss function. Mean squared error is an expected squared deviation, measured from the prediction rather than the mean; cross-entropy is an expected negative log-likelihood.
• The bias-variance decomposition, which is an identity in these two quantities and nothing else.
• Standard errors and intervals, where the whole game is estimating a standard deviation of an estimate [2, §8.3].
• Portfolio and risk work, where Var(X + Y) = Var(X) + Var(Y) + 2Cov(X, Y) is the entire subject rather than a footnote.
• Effect sizes, which exist to "compare the difference between groups to the variability within groups" — Cohen's d divides a difference in means by a pooled standard deviation, so the answer is in standard deviations rather than in whatever the data was measured in [2, §2.9].

Why & when

Report both, always. A mean alone is a sentence with the verb missing.

Check independence before adding variances. This is the most common error of the three, and the worked example measures how large it gets: a correlated pair whose variances sum to 20.2 has a sum with variance 36.2, because in general Var(X + Y) = Var(X) + Var(Y) + 2Cov(X, Y) [1, §22.6.1.10].

Use Chebyshev when you know nothing, and expect it to be loose. It is a guarantee that holds for every distribution, so it has to hold for the worst one: d2l's three-point example sits exactly on it [1, §22.6.1.6], and most distributions sit well inside.

Do not use the standard deviation as a summary of a heavy tail. If the tail is heavy enough, the standard deviation does not exist, and the number you computed is an artefact of your sample — a thousand points estimate it no better than one [3, Comments].

Standardise when the units are meaningless, and remember that the variance is in squared units, so z-scoring is division by σ and never by σ² [1, §22.6.1.5].

Worked example

Three questions. Do the rules hold when measured? How loose is Chebyshev in practice? And what does a distribution without a mean actually look like?

import numpy as np

rng = np.random.default_rng(0)
N = 2_000_000

# --- 1. the rules, measured rather than quoted ------------------------------
X = rng.normal(3.0, 2.0, N)
Y = rng.gamma(2.0, 1.5, N)  # independent of X, and not symmetric
Z = 2.0 * X + rng.normal(0, 0.5, N)  # deliberately correlated with X

print("  E and Var over 2,000,000 draws")
print(f"    E[5X + 1]        {(5 * X + 1).mean():>9.4f}   5 E[X] + 1     {5 * X.mean() + 1:>9.4f}")
print(f"    Var(5X + 1)      {(5 * X + 1).var():>9.4f}   25 Var(X)      {25 * X.var():>9.4f}")
print(f"    Var(X + Y)       {(X + Y).var():>9.4f}   Var X + Var Y  {X.var() + Y.var():>9.4f}   independent")
print(f"    Var(X + Z)       {(X + Z).var():>9.4f}   Var X + Var Z  {X.var() + Z.var():>9.4f}   correlated")
print(f"      the shortfall is 2 Cov(X, Z) = {2 * np.cov(X, Z)[0, 1]:.4f}")

# --- 2. how loose is Chebyshev? --------------------------------------------
# The bound is a statement about a distribution, so the true mean and standard
# deviation go in, not the sample's.
cases = [
    ("normal(0,1)", lambda n: rng.normal(0, 1, n), 0.0, 1.0),
    ("exponential(1)", lambda n: rng.exponential(1.0, n), 1.0, 1.0),
    ("lognormal(0,1)", lambda n: rng.lognormal(0, 1, n), np.exp(0.5), np.sqrt(np.e**2 - np.e)),
    ("t, 3 dof", lambda n: rng.standard_t(3, n), 0.0, np.sqrt(3.0)),
]
print("\n  P(|X - mu| >= alpha sigma), measured against the 1/alpha^2 bound")
print(f"{'alpha':>8}{'bound':>10}" + "".join(f"{c[0]:>17}" for c in cases))
for alpha in (1.5, 2.0, 3.0, 5.0):
    row = "".join(f"{np.mean(np.abs(d(N) - mu) >= alpha * sd):>17.4f}" for _, d, mu, sd in cases)
    print(f"{alpha:>8.1f}{1 / alpha**2:>10.4f}" + row)

# --- 3. the bound is not loose for every distribution -----------------------
# d2l's three-point variable: +/-2 with probability p each, 0 otherwise. Its
# variance is 8p, so with p = 1/(2 alpha^2) the two outer points sit exactly
# alpha standard deviations from the mean. Computed from the pmf, not sampled.
print("\n  a family that attains the bound, computed exactly from its pmf")
for alpha in (1.5, 2.0, 3.0, 5.0):
    p = 1 / (2 * alpha**2)
    sigma = np.sqrt(8 * p)
    print(f"    alpha {alpha:>4.1f}   p {p:.5f}   sigma {sigma:.4f}   alpha*sigma {alpha * sigma:.4f}   P = 2p = {2 * p:.4f}   bound {1 / alpha**2:.4f}")

# --- 4. a distribution with no mean at all ---------------------------------
# Eight independent runs of each, so the comparison is a property of the
# distribution rather than of one lucky stream.
rng4 = np.random.default_rng(1)  # its own seed, so this block reproduces alone
RUNS, M = 8, 1_000_000
norm = rng4.normal(0, 1, (RUNS, M)).cumsum(1) / np.arange(1, M + 1)
cauchy = rng4.standard_cauchy((RUNS, M)).cumsum(1) / np.arange(1, M + 1)
print(f"\n  running mean of n draws, {RUNS} independent runs of each")
print(f"{'n':>12}{'normal: worst |mean|':>24}{'Cauchy: worst |mean|':>24}")
for n in (10, 100, 1_000, 10_000, 100_000, 1_000_000):
    print(f"{n:>12,}{np.abs(norm[:, n - 1]).max():>24.4f}{np.abs(cauchy[:, n - 1]).max():>24.4f}")
sums = cauchy * np.arange(1, M + 1)
print(f"    largest single Cauchy draw seen {np.abs(np.diff(sums, axis=1)).max():>16,.0f}")

  E and Var over 2,000,000 draws
    E[5X + 1]          16.0090   5 E[X] + 1       16.0090
    Var(5X + 1)        99.9654   25 Var(X)        99.9654
    Var(X + Y)          8.4965   Var X + Var Y     8.5016   independent
    Var(X + Z)         36.2382   Var X + Var Z    20.2433   correlated
      the shortfall is 2 Cov(X, Z) = 15.9949

  P(|X - mu| >= alpha sigma), measured against the 1/alpha^2 bound
   alpha     bound      normal(0,1)   exponential(1)   lognormal(0,1)         t, 3 dof
     1.5    0.4444           0.1334           0.0818           0.0562           0.0806
     2.0    0.2500           0.0458           0.0500           0.0370           0.0407
     3.0    0.1111           0.0027           0.0183           0.0181           0.0138
     5.0    0.0400           0.0000           0.0025           0.0058           0.0033

  a family that attains the bound, computed exactly from its pmf
    alpha  1.5   p 0.22222   sigma 1.3333   alpha*sigma 2.0000   P = 2p = 0.4444   bound 0.4444
    alpha  2.0   p 0.12500   sigma 1.0000   alpha*sigma 2.0000   P = 2p = 0.2500   bound 0.2500
    alpha  3.0   p 0.05556   sigma 0.6667   alpha*sigma 2.0000   P = 2p = 0.1111   bound 0.1111
    alpha  5.0   p 0.02000   sigma 0.4000   alpha*sigma 2.0000   P = 2p = 0.0400   bound 0.0400

  running mean of n draws, 8 independent runs of each
           n    normal: worst |mean|    Cauchy: worst |mean|
          10                  0.4403                  2.2880
         100                  0.1060                  2.6014
       1,000                  0.0543                  5.9038
      10,000                  0.0214                  5.5972
     100,000                  0.0049                  1.9442
   1,000,000                  0.0017                 10.8063
    largest single Cauchy draw seen       10,244,088

The linear rules hold exactly, to four decimals on two million draws: E[5X+1] = 16.0090 both ways, Var(5X+1) = 99.9654 both ways. Independent variances add: 8.4965 against 8.5016, a 0.06% gap that is sampling noise.

Correlated variances do not. Var(X + Z) is 36.24 where the sum of the variances is 20.24, and the missing 16.00 is exactly 2 Cov(X, Z) = 15.9949. This is not a small correction. The rule that people carry around — variances add — understates the true spread here by 44%, and in the direction that matters: it says the sum is less variable than it is.

Chebyshev is loose by a factor of three to twelve. At α = 2 the bound allows 25% of the mass outside, and the normal puts 4.58% there. At α = 5 the bound allows 4% and the normal has none. This is what "for any distribution" costs: the bound is set by the worst case.

It is not loose for the worst case. The three-point family sits exactly on the bound at every α tried. So Chebyshev cannot be improved without assuming something, and every tighter rule you have seen — the 68-95-99.7 figures, for instance — is a statement about the normal distribution, not about spread.

And then there is the Cauchy. The normal's worst running mean over eight runs falls from 0.44 to 0.0017 as n goes from ten to a million. The Cauchy's goes from 2.29 to 10.81 — it gets worse. The largest single draw in those eight million numbers was 10,244,088, and one value like that can move the mean of a million draws by ten.

The uncomfortable part is that nothing in the sample announces this. Compute a mean and a standard deviation of Cauchy data and you will get numbers; they are summaries of your particular draws and of nothing else, and they will not stabilise no matter how much data you add.

Pitfalls

Adding variances without checking independence. The single most expensive error here, and the worked example puts a number on it.

Assuming the standard deviation bounds the tail tightly. Chebyshev is all you are entitled to without a distributional assumption, and it cannot be tightened [1, §22.6.1.6]; the worked example finds it loose by a factor of three to twelve.

Using a mean on heavy-tailed data — latency, income, claim sizes, file sizes. If it exists at all, it is dominated by observations you probably have not seen yet. Report a median and some quantiles too.

Mixing up the two sample variances. np.var divides by n and pandas.Series.var by n−1. The names are the same [2, §8.2]; the numbers are not.

Treating the variance as if it were in the data's units. It is in squared units [1, §22.6.1.5], which is why you standardise by σ.

Reading "variance" as "risk" without checking symmetry. Variance counts a surprise upward the same as a surprise downward, and most applications do not.