5.8.1 Estimators: bias, variance and consistency

An estimator is a random variable

A distribution has fixed numbers that describe it, such as its mean $\mu$ , its variance $\sigma^2$ or a coin's chance of heads $p$ . Such a number is a parameter, written $\theta$ (theta) in general. In practice $\theta$ is unknown, and data are used to guess it.

The data are a sample $X_1, \ldots, X_n$ : $n$ independent draws, each with the same distribution. This is often shortened to i.i.d., independent and identically distributed.

Rule: Estimator and estimate

\hat{\theta} = g(X_1, \ldots, X_n)

An estimator is a rule $g$ applied to the sample; the hat marks it. Before the draws it is a random variable, with a mean, a variance and a distribution of its own. Its value for one observed sample $x_1, \ldots, x_n$ is an estimate, a single number.

Take a fair die, whose mean $\mu = 3.5$ we pretend not to know, and the estimator $\bar{X} = \frac{1}{2}(X_1 + X_2)$ from two rolls. The 36 equally likely pairs give 11 possible values:

$\bar{x}$	1	1.5	2	2.5	3	3.5	4	4.5	5	5.5	6
pairs	1	2	3	4	5	6	5	4	3	2	1

This is the two-dice total of lesson 5.6.3, halved. Its mean is $3.5$ , the true $\mu$ , and its variance is $\frac{35}{24} \approx 1.4583$ .

Five seeded samples of two rolls give the estimates 2.5, 1, 4.5, 3.5 and 3: one rule, five different numbers. The distribution of an estimator over all possible samples is its sampling distribution, which lesson 5.8.2 studies.

\bar{x}

from fractions import Fraction from itertools import product import numpy as np dist = {} for a, b in product(range(1, 7), repeat=2): xbar = Fraction(a + b, 2) dist[xbar] = dist.get(xbar, 0) + Fraction(1, 36) print([int(36 * q) for x, q in sorted(dist.items())]) mean = sum(x * q for x, q in dist.items()) var = sum((x - mean) ** 2 * q for x, q in dist.items()) print(mean, var, round(float(var), 4)) rng = np.random.default_rng(581) print(rng.integers(1, 7, size=(5, 2)).mean(axis=1))

⌘/Ctrl + Enter runs the cell

Bias, variance and mean squared error

A good estimator lands close to $\theta$ . Three numbers measure how close: where it is centred, how widely it scatters, and the two combined.

Rule: Bias, variance and MSE

\begin{gathered} \operatorname{bias}(\hat{\theta}) = \mathbb{E}[\hat{\theta}] - \theta \\[4pt] \operatorname{MSE}(\hat{\theta}) = \mathbb{E}\big[(\hat{\theta} - \theta)^2\big] \\[4pt] = \operatorname{Var}(\hat{\theta}) + \operatorname{bias}(\hat{\theta})^2 \end{gathered}

An estimator is unbiased when its bias is 0 for every value of $\theta$ . The mean squared error (MSE) is in $\theta$ 's units squared; its square root is in $\theta$ 's units. The split is lesson 5.3.2's $\mathbb{E}[Y^2] = \operatorname{Var}(Y) + (\mathbb{E}[Y])^2$ with $Y = \hat{\theta} - \theta$ .

Toss a coin 10 times to estimate its chance of heads $p$ . The number of heads $H$ is $\operatorname{Bin}(10, p)$ , with mean $10p$ and variance $10p(1 - p)$ (lesson 5.3.2). Compare two estimators.

The share of heads, $\hat{p} = \frac{H}{10}$ . Its mean is $p$ , so it is unbiased, and its variance is $\frac{p(1 - p)}{10}$ .
Add one head and one tail, $\tilde{p} = \frac{H + 1}{12}$ . Its mean is $\frac{10p + 1}{12}$ , pulled towards $\frac{1}{2}$ , so its bias is $\frac{1 - 2p}{12}$ . Its variance is $\frac{10p(1 - p)}{144}$ , smaller than $\hat{p}$ 's.

$p$	MSE of $\hat{p}$	bias of $\tilde{p}$	MSE of $\tilde{p}$
0.5	0.025	0	0.017361
0.3	0.021	0.0333	0.015694
0.1	0.009	0.0667	0.010694

Near $\frac{1}{2}$ , the small bias of $\tilde{p}$ buys a large cut in variance, and it has the smaller error. Near 0 its pull towards $\frac{1}{2}$ costs more than it saves, and the unbiased $\hat{p}$ wins. Neither is better for every $p$ ; case study 4.6.5 shrank estimated returns for the same reason.

H

import numpy as np from scipy.stats import binom k = np.arange(11) for p in [0.5, 0.3, 0.1]: pmf = binom.pmf(k, 10, p) for est in [k / 10, (k + 1) / 12]: mean = est @ pmf var = ((est - mean) ** 2) @ pmf mse = ((est - p) ** 2) @ pmf print(p, round(mean - p, 4) + 0.0, round(var, 6), round(mse, 6)) H = np.random.default_rng(581).binomial(10, 0.3, size=100_000) print(round(np.mean((H / 10 - 0.3) ** 2), 4), round(np.mean(((H + 1) / 12 - 0.3) ** 2), 4))

⌘/Ctrl + Enter runs the cell

The sample mean

The usual estimator of $\mu$ is the sample mean $\bar{X} = \frac{1}{n}(X_1 + \cdots + X_n)$ .

Rule: The sample mean

\begin{gathered} \mathbb{E}[\bar{X}] = \mu \\[4pt] \operatorname{Var}(\bar{X}) = \frac{\sigma^2}{n} \end{gathered}

Unbiased for any distribution with a mean, by linearity (lesson 5.3.2). For independent draws the $n$ variances add and the factor $\frac{1}{n}$ comes out squared (lesson 5.7.2). So its MSE is $\frac{\sigma^2}{n}$ , and its standard deviation $\frac{\sigma}{\sqrt{n}}$ is the standard error of lesson 5.8.2.

For the two rolls above, $\frac{\sigma^2}{n} = \frac{35/12}{2} = \frac{35}{24}$ , the variance the table gave.

Being unbiased is not enough to make an estimator good. From three rolls of a die, here are three estimators of $\mu$ , all unbiased:

estimator	mean	variance
the first roll, $X_1$	3.5	$\frac{35}{12} \approx 2.9167$
the middle value of the three	3.5	$\frac{203}{108} \approx 1.8796$
the sample mean $\bar{X}$	3.5	$\frac{35}{36} \approx 0.9722$

The middle value, the sample median, is centred on 3.5 because the die is symmetric. All three are centred correctly; $\bar{X}$ scatters least, so it has the smallest MSE.

\frac{1}{216}

from itertools import product import numpy as np rolls = np.array(list(product(range(1, 7), repeat=3))) estimators = {"first roll": rolls[:, 0], "median": np.median(rolls, axis=1), "mean": rolls.mean(axis=1)} for name, est in estimators.items(): print(name, est.mean(), round(est.var(), 4))

⌘/Ctrl + Enter runs the cell

Why the sample variance divides by $n - 1$

The variance $\sigma^2$ is the mean squared deviation from $\mu$ . With $\mu$ unknown, the deviations are measured from $\bar{X}$ instead. But $\bar{X}$ is the centre of the sample itself, so the sample sits closer to $\bar{X}$ than to $\mu$ . The squared deviations then come out too small on average.

By a table. A fair coin scores 0 or 1, so $\mu = \frac{1}{2}$ and $\sigma^2 = \frac{1}{4}$ . Take samples of $n = 2$ ; the four outcomes are equally likely.

sample	$\bar{x}$	$\sum (x_i - \bar{x})^2$	divide by $n = 2$	divide by $n - 1 = 1$
(0, 0)	0	0	0	0
(0, 1)	0.5	0.5	0.25	0.5
(1, 0)	0.5	0.5	0.25	0.5
(1, 1)	1	0	0	0
average			$\frac{1}{8}$	$\frac{1}{4}$

Dividing by $n$ gives $\frac{1}{8}$ on average, half the true $\frac{1}{4}$ . Dividing by $n - 1$ gives exactly $\frac{1}{4}$ .

By algebra. The sum of squared deviations from $\bar{X}$ has a mean that can be worked out for any distribution with a finite variance:

1
$\sum_i (X_i - \bar{X})^2 = \sum_i (X_i - \mu)^2 - n(\bar{X} - \mu)^2$
2
$\mathbb{E}\big[\sum_i (X_i - \mu)^2\big] = n\sigma^2$
3
$\mathbb{E}\big[n(\bar{X} - \mu)^2\big] = \sigma^2$
4
$\mathbb{E}\big[\sum_i (X_i - \bar{X})^2\big] = (n - 1)\sigma^2$

Measuring from $\bar{X}$ loses exactly one $\sigma^2$ : the part that is $\bar{X}$ 's own error. Dividing the sum by $n - 1$ puts it back.

Rule: Sample variance

\begin{gathered} s^2 = \frac{1}{n - 1}\textstyle\sum_i (X_i - \bar{X})^2 \\[4pt] \mathbb{E}[s^2] = \sigma^2 \\[8pt] \hat{\sigma}^2 = \frac{1}{n}\textstyle\sum_i (X_i - \bar{X})^2 \\[4pt] \mathbb{E}[\hat{\sigma}^2] = \frac{n - 1}{n}\,\sigma^2 \end{gathered}

$s^2$ is unbiased for any distribution with a finite variance; $\hat{\sigma}^2$ has bias $-\frac{\sigma^2}{n}$ . The number $n - 1$ is the degrees of freedom: the $n$ deviations from $\bar{X}$ add to 0, so only $n - 1$ of them are free to vary.

This is why lesson 5.5.2's sample variance and covariance, np.cov, and np.var with ddof=1 divide by $n - 1$ , as case studies 3.1.4 and 3.2.5 did. For case study 3.1.4's ten months of X, np.var gives 0.002276 and ddof=1 gives 0.002528, larger by the factor $\frac{10}{9}$ .

Estimates of $\sigma^2 = 1$ from samples of 5 normal draws. Dividing by $n - 1$ centres the estimates on 1; dividing by $n$ centres them on 0.8. Either way most single samples fall below 1.

\bar{x}

from itertools import product import numpy as np for sample in product([0, 1], repeat=2): x = np.array(sample) ss = ((x - x.mean()) ** 2).sum() print(sample, x.mean(), ss, ss / 2, ss / 1) x, mu = np.array([2.0, 5.0, 3.0, 6.0]), 3.5 print(((x - x.mean()) ** 2).sum(), ((x - mu) ** 2).sum() - len(x) * (x.mean() - mu) ** 2) rng = np.random.default_rng(5811) normal = rng.standard_normal((100_000, 5)) print(round(normal.var(axis=1).mean(), 4), round(normal.var(axis=1, ddof=1).mean(), 4)) dice = rng.integers(1, 7, size=(100_000, 5)) print(round(dice.var(axis=1).mean(), 4), round(dice.var(axis=1, ddof=1).mean(), 4))

⌘/Ctrl + Enter runs the cell

A known mean, the standard deviation, and higher moments

A known mean. If $\mu$ is known, measure from $\mu$ and divide by $n$ : $\frac{1}{n}\sum_i (X_i - \mu)^2$ is unbiased, by the second step above. The $n - 1$ pays for estimating $\mu$ .
The standard deviation. $s = \sqrt{s^2}$ is biased low even though $s^2$ is not, because the square root bends downwards (Jensen's inequality, lesson 5.7.1). For the coin with $n = 2$ , $s$ is 0 or $\sqrt{0.5}$ , so $\mathbb{E}[s] \approx 0.3536$ against $\sigma = 0.5$ . For 5 normal draws, $\mathbb{E}[s] \approx 0.94\sigma$ . The bias fades as $n$ grows, and $s$ is used anyway.
Skewness and kurtosis. Lesson 5.3.3's sample skewness and kurtosis divide by $n$ , and are biased in small samples. bias=False in SciPy multiplies them by a correction factor, which reduces the bias but does not remove it.

Exponential draws have skewness 2 and excess kurtosis 6. In 100,000 seeded samples of 10, the sample skewness averages 0.974, or 1.155 with bias=False; the excess kurtosis averages 0.18, or 1.282. Higher moments need large samples.

n = 5

import numpy as np from scipy import stats normal = np.random.default_rng(5813).standard_normal((100_000, 5)) print(round((normal**2).mean(axis=1).mean(), 4)) print(round(np.sqrt(0.5) / 2, 4), round(normal.std(axis=1, ddof=1).mean(), 4)) rng = np.random.default_rng(449) market = rng.normal(0.008, 0.04, 10) x = np.round(market + rng.normal(0.002, 0.025, 10), 3) print(round(np.var(x), 6), round(np.var(x, ddof=1), 6), round(np.var(x, ddof=1) / np.var(x), 4)) expo = np.random.default_rng(5814).exponential(1.0, size=(100_000, 10)) print(round(stats.skew(expo, axis=1).mean(), 3), round(stats.skew(expo, axis=1, bias=False).mean(), 3)) print(round(stats.kurtosis(expo, axis=1).mean(), 3), round(stats.kurtosis(expo, axis=1, bias=False).mean(), 3))

⌘/Ctrl + Enter runs the cell

Unbiased is not always the most accurate

Divide the sum of squares $\sum_i (X_i - \bar{X})^2$ by $d$ . For normal draws the MSE of each choice of $d$ is known exactly:

divisor $d$	MSE	MSE for $n = 5$ , $\sigma = 1$
$n - 1$	$\frac{2\sigma^4}{n - 1}$	0.5
$n$	$\frac{(2n - 1)\sigma^4}{n^2}$	0.36
$n + 1$	$\frac{2\sigma^4}{n + 1}$	0.3333

The unbiased $s^2$ has the largest error of the three: its spread costs more than the others' bias. It stays the standard all the same, because the intervals and tests of lessons 5.8.2 and 5.8.3 are built on it. For large $n$ the three are almost equal.

n = 5

import numpy as np n = 5 print(2 / (n - 1), (2 * n - 1) / n**2, round(2 / (n + 1), 4)) X = np.random.default_rng(5812).standard_normal((200_000, n)) ss = ((X - X.mean(axis=1, keepdims=True)) ** 2).sum(axis=1) for d in [n - 1, n, n + 1]: print(d, round(np.mean(ss / d), 4), round(np.mean((ss / d - 1) ** 2), 4))

⌘/Ctrl + Enter runs the cell

Consistency

Bias and variance describe an estimator for one sample size. Consistency asks what happens as the sample grows: does the estimator close in on the truth?

Rule: Consistency

\begin{gathered} \lim_{n \to \infty} P\big(\lvert \hat{\theta}_n - \theta \rvert > \varepsilon\big) = 0 \\[4pt] \text{for every } \varepsilon > 0 \end{gathered}

$\hat{\theta}_n$ is the estimator from $n$ draws and $\varepsilon$ (epsilon) any tolerance, however small. If its bias and its variance both tend to 0, it is consistent: Markov's inequality (lesson 5.7.1) applied to $(\hat{\theta}_n - \theta)^2$ gives $P(\lvert \hat{\theta}_n - \theta \rvert \ge \varepsilon) \le \frac{\operatorname{MSE}}{\varepsilon^2}$ .

The sample mean has MSE $\frac{\sigma^2}{n} \to 0$ , so it is consistent: the law of large numbers of lesson 5.7.2.
Both variance estimators are consistent. For normal draws with $\sigma^2 = 1$ , the chance that $s^2$ misses by more than 0.1 is about 0.83, 0.48 and 0.02 for $n$ = 10, 100 and 1,000, in one seeded run. $\hat{\sigma}^2$ is biased for every $n$ , but its bias $-\frac{\sigma^2}{n}$ vanishes.
The first draw alone is unbiased but not consistent. For $\mathcal{N}(\mu, 1)$ draws, $P(\lvert X_1 - \mu \rvert > 0.1) = 0.9203$ whatever $n$ is: the other draws are ignored.
The mean of a Cauchy sample (the t with 1 degree of freedom, lesson 5.4.4) is not consistent. That distribution has no mean, and the average of $n$ Cauchy draws is again Cauchy. So $P(\lvert \bar{X} \rvert > 1) = 0.5$ for every $n$ .

estimator	unbiased	consistent
$\bar{X}$	yes	yes
$s^2$	yes	yes
$\hat{\sigma}^2$	no	yes
$X_1$ alone	yes	no

How the estimates of $\sigma^2 = 1$ spread as the sample grows. Ninety percent of samples land inside the band, which closes in on 1; the divide-by- $n$ estimate's mean $\frac{n - 1}{n}$ rises to meet it.

s^2

import numpy as np from scipy.stats import norm rng = np.random.default_rng(5816) for n in [10, 100, 1000]: s2 = rng.standard_normal((4000, n)).var(axis=1, ddof=1) print(n, round(np.mean(np.abs(s2 - 1) > 0.1), 3)) print(round(2 * norm.sf(0.1), 4)) cauchy = np.random.default_rng(5815) for n in [10, 100, 1000]: means = cauchy.standard_cauchy((4000, n)).mean(axis=1) print(n, round(np.mean(np.abs(means) > 1), 3), round(2 * norm.sf(0.1 * np.sqrt(n)), 4))

⌘/Ctrl + Enter runs the cell

Lesson 5.8.2 turns the spread of an estimator into an interval of plausible values, and case study 5.8.5 does this for a Sharpe ratio.

Common mistakes

An estimate is one number; the estimator is a random variable. A different sample gives a different estimate, and its spread is part of the answer.
Unbiased is not the same as accurate. The first roll alone is unbiased; compare estimators by their MSE.
Divide by $n - 1$ when the mean is estimated. With a known $\mu$ , divide by $n$ .
$s^2$ is unbiased for $\sigma^2$ , but $s$ is not unbiased for $\sigma$ . An unbiased estimate of $\theta$ does not give an unbiased estimate of $\sqrt{\theta}$ .
Check the divisor. np.var and np.std divide by $n$ unless given ddof=1; np.cov divides by $n - 1$ .
Unbiased and consistent are different properties. $X_1$ alone is unbiased but not consistent; $\hat{\sigma}^2$ is consistent but biased.

Estimators: bias, variance and consistency

0 of 12 correct

1.
An estimator has mean 5.2 when the true value is 5. Find its bias
2.
An estimator has bias 0.3 and variance 0.16. Find its mean squared error
3.
An unbiased estimator has MSE 0.09. Find its standard deviation
4.
A coin is tossed 10 times and $\tilde{p} = \frac{H + 1}{12}$ . Find the bias of $\tilde{p}$ when $p = 0.2$
5.
Draws have variance $\sigma^2 = 9$ . Find the variance of the mean of 4 of them
6.
How many draws make the standard deviation of $\bar{X}$ one tenth of $\sigma$ ?
7.
Draws have variance 12. Find the mean of the divide-by- $n$ estimate $\hat{\sigma}^2$ from samples of 4
8.
Find $s^2$ for the sample 1, 3, 5, dividing by $n - 1$
9.
np.var(x) returns 2.4 for a sample of 6 values. What does np.var(x, ddof=1) return?
10.
The mean $\mu = 0$ is known. Find the unbiased variance estimate from the sample $-1, 3, 2$
11.
A coin scores 0 or 1. For samples of 2, what is the mean of $s$ , the square root of $s^2$ ? (type sqrt)
12.
An estimator has bias $\frac{2}{n}$ and variance $\frac{3}{n}$ . Find its MSE when $n = 100$

How to type answers

Powers with ^: x^7, 2^(1/3), a^-2.
Multiply with * or by writing side by side: 3x^2, 2(x+1).
ln(x) is the natural log, log(x, 2) is log base 2, e is Euler's number, sqrt(x) is a square root.
sin(x), cos(x) and tan(x) take radians; write pi for π and sin(x)^2 for sin²x.
Your first answer to each question is saved: it counts towards your score and brings the question back for review.
Exact answers like ln(20)/ln(3) are best. Decimals need about 4 significant figures.

\hat{\theta} = g(X_1, \ldots, X_n)

Estimators: bias, variance and consistency

An estimator is a random variable

Rule: Estimator and estimate

Bias, variance and mean squared error

Rule: Bias, variance and MSE

The sample mean

Rule: The sample mean

Why the sample variance divides by $n - 1$

Rule: Sample variance

A known mean, the standard deviation, and higher moments

Unbiased is not always the most accurate

Consistency

Rule: Consistency

Common mistakes

Estimators: bias, variance and consistency

Summary card

Free further reading

An estimator is a random variable

Rule: Estimator and estimate

Bias, variance and mean squared error

Rule: Bias, variance and MSE

The sample mean

Rule: The sample mean

Why the sample variance divides by n−1n - 1n−1

Rule: Sample variance

A known mean, the standard deviation, and higher moments

Unbiased is not always the most accurate

Consistency

Rule: Consistency

Common mistakes

Summary card

Free further reading

Why the sample variance divides by $n - 1$