Statistical inference

Estimation, fitting and testing.

In this chapter 8

Tick a unit to add it to your PDF basket.

Estimation and testing

Estimation and testing

Definition

Sampling distributions

Adapted

Add to PDF

Defines a statistic and its sampling distribution, fixes the two rival conventions for the sample variance, and proves that for a normal sample the sample mean and the sample variance are independent with laws $N(\mu,\sigma^2/n)$ and $\sigma^2\chi^2(n-1)/(n-1)$ , from which the Student $t$ and Fisher $F$ densities follow; worked on an eight-point sample giving $t=1.063$ and $\chi^2=5.731$ .

Prerequisites: The multivariate normal distribution, Expectation, variance and covariance, Joint distributions and transformations, The central limit theorem, Unitary and orthogonal matrices, Multiple integrals and change of variables, The Gamma function, The Beta function

Note(What this unit is for)Original

Every number computed from data - a mean, a variance, a fitted slope - is a function of random variables and is therefore itself random. It has a distribution, and that distribution, not the number, is what licences any statement about the quantity being measured. The distribution of a statistic is called its sampling distribution, and the whole of estimation, interval construction and testing is bookkeeping on sampling distributions.

For one parent population - the normal - the sampling distributions of the two statistics a physicist computes first, the sample mean and the sample variance, are known exactly, for every sample size, in closed form. That is the content of this unit. Three named distributions appear: chi-squared for the sample variance, Student's $t$ for the mean divided by its own estimated error, and Fisher's $F$ for the ratio of two sample variances. Their appearance is not a coincidence of calculus; it comes from one geometric fact, that the standard normal law on $\mathbb{R}^n$ is invariant under rotations.

Definition(Sample, statistic, sampling distribution)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.3, the set-up on pp. 1229-1230; Grimmett & Stirzaker (2020), §4.10, the opening paragraph, p. 135

A random sample of size $n$ from a distribution $F$ is a list $X_1, \dots, X_n$ of independent random variables each with distribution $F$ . Riley, Hobson and Bence call $F$ the population and write it $P(x\mid a)$ when it is known up to a list $a = (a_1,\dots,a_m)$ of unknown parameters; independence then makes the joint density of the sample

P(x \mid a) \;=\; \prod_{i=1}^{n} P(x_i \mid a) , \qquad x = (x_1,\dots,x_n) \in \mathbb{R}^{n} .

A statistic is any function $T(X_1,\dots,X_n)$ of the sample alone - it may not involve the unknown $a$ , because it has to be computable from the data. Being a function of random variables, $T$ is a random variable, and its distribution is the sampling distribution of $T$ . Riley, Hobson and Bence specify it by the transformation rule $P(\hat a \mid a)\,d\hat a = P(x\mid a)\,d^{n}x$ , where $d^{n}x$ is the volume between the surfaces $T(x) = \hat a$ and $T(x) = \hat a + d\hat a$ ; this is the change-of-variables formula of Multiple integrals and change of variables written for one coordinate of the image. The sampling distribution depends on the parent $F$ , on the function $T$ and on $n$ , and in general on the true values $a$ as well - which is exactly why a statistic can carry information about $a$ .

Definition(Sample mean and sample variance, and the two conventions)AdaptedAdapted from Grimmett & Stirzaker (2020), §4.10, the definitions of $\overline X$ and $S^2$ with the footnote, p. 135; Riley, Hobson & Bence (2006), §31.2.2, eqs. (31.6)-(31.8), pp. 1224-1225; Arfken, Weber & Harris (2013), §23.7, eqs. (23.94)-(23.95), p. 1168

For a sample $X_1,\dots,X_n$ the sample mean is

\overline X \;=\; \frac1n \sum_{i=1}^{n} X_i ,

on which all our sources agree. For the spread there are two conventions in the library, and they differ by a factor $n/(n-1)$ :

S^{2} \;=\; \frac{1}{n-1}\sum_{i=1}^{n}\bigl(X_i - \overline X\bigr)^{2} \qquad\text{(Grimmett and Stirzaker)} ,

s^{2} \;=\; \frac{1}{n}\sum_{i=1}^{n}\bigl(X_i - \overline X\bigr)^{2} \qquad\text{(Riley, Hobson and Bence)} .

This unit uses $S^{2}$ , with the divisor $n-1$ . Grimmett and Stirzaker leave the reasons for this choice to statistics texts, note that with it $\mathbb{E}(S^{2}) = \sigma^{2}$ exactly, and record the rival convention in a footnote. Arfken, Weber and Harris reach the same divisor by a first-order argument: estimating the variance of $\bar x$ by $S^{2}/n$ , the deviations $e_j = x_j - \bar x$ satisfy $\overline{e^{2}} = S^{2}(1 - 1/n)$ (their sample standard deviation $s$ is our $S$ ). The $n-1$ that results, they say, allows for the probable error in $\bar x$ , and they call it Bessel's correction. In the language of the proof below, one degree of freedom has been spent on estimating the mean. Note that $S^{2}$ is undefined at $n = 1$ , as it should be: one reading says nothing about spread. Both multiply out to the same sum of squared deviations, so the conversion between them is

(n-1)\,S^{2} \;=\; n\,s^{2} \;=\; \sum_{i=1}^{n}\bigl(X_i-\overline X\bigr)^{2} ,

and every formula below can be rewritten in the other convention by that identity. Riley, Hobson and Bence also give the computing form $s^{2} = \overline{x^{2}} - \bar x^{2}$ , which needs one pass through the data instead of two.

Two facts hold for any parent population with finite mean $\mu$ and variance $\sigma^{2}$ , normal or not, and are proved in Expectation, variance and covariance from linearity and the additivity of variances over independent summands:

\mathbb{E}\bigl(\overline X\bigr) = \mu , \qquad \operatorname{var}\bigl(\overline X\bigr) = \frac{\sigma^{2}}{n} .

Everything else below needs normality.

Definition(The chi-squared,

t

and

F

distributions)AdaptedAdapted from Grimmett & Stirzaker (2020), §4.4, Examples (4.4.6) and (4.4.8), p. 108, §4.7, Example (4.7.1), p. 122, and §4.10, the definitions of $t(r)$ and $F(r,s)$ with their densities, pp. 136-137; Riley, Hobson & Bence (2006), §30.9.4, eqs. (30.122)-(30.123) and the moments, p. 1193

Let $d, r, s$ be positive integers.

Chi-squared. $X$ has the chi-squared distribution with $d$ degrees of freedom, written $\chi^{2}(d)$ , if it is a gamma variate with parameters $\lambda = \tfrac12$ and $t = \tfrac12 d$ , that is if it has density

f(x) \;=\; \frac{1}{2^{d/2}\,\Gamma(d/2)}\, x^{d/2 - 1} e^{-x/2} , \qquad x > 0 ,

with $\Gamma$ the gamma function of The Gamma function. Its mean is $d$ and its variance $2d$ . Equivalently, and this is the form used throughout statistics, $\chi^{2}(d)$ is the law of $Z_1^{2} + \cdots + Z_d^{2}$ for independent standard normal $Z_i$ : the gamma moment generating function is multiplicative in the shape parameter, so a sum of independent $\chi^{2}(d_1)$ and $\chi^{2}(d_2)$ variates is $\chi^{2}(d_1+d_2)$ , and one standard normal square is $\chi^{2}(1)$ : for $y > 0$ , $\Pr(Z^{2} \le y) = \Pr(-\sqrt y \le Z \le \sqrt y)$ , whose derivative is $(2\pi y)^{-1/2}e^{-y/2}$ , the density above with $d = 1$ (Grimmett and Stirzaker, Example (4.7.1), p. 122). Riley, Hobson and Bence write the same object as $\chi_n^{2} = \sum_i (X_i - \mu_i)^{2}/\sigma_i^{2}$ for independent $X_i \sim N(\mu_i, \sigma_i^{2})$ , and state that if the $X_i$ are subject to $m$ independent linear constraints then the order drops to $n - m$ . The theorem below is the leading instance of that rule.

Student's $t$ . If $V \sim N(0,1)$ and $U \sim \chi^{2}(r)$ are independent, then

T \;=\; \frac{V}{\sqrt{U/r}}

has the $t$ distribution with $r$ degrees of freedom, written $t(r)$ , with density

f_T(t) \;=\; \frac{\Gamma\bigl(\tfrac12(r+1)\bigr)}{\sqrt{\pi r}\;\Gamma\bigl(\tfrac12 r\bigr)} \left(1 + \frac{t^{2}}{r}\right)^{-\frac12 (r+1)} , \qquad -\infty < t < \infty .

It is symmetric about $0$ ; $t(1)$ is the Cauchy distribution, which has no mean, and for $r \ge 2$ the mean is $0$ . As $r \to \infty$ the density tends to $(2\pi)^{-1/2}e^{-t^{2}/2}$ , since $(1+t^{2}/r)^{-r/2} \to e^{-t^{2}/2}$ .

Fisher's $F$ . If $U \sim \chi^{2}(r)$ and $V \sim \chi^{2}(s)$ are independent, then

F \;=\; \frac{U/r}{V/s}

has the $F$ distribution with $r$ and $s$ degrees of freedom, written $F(r,s)$ , with density

f_F(x) \;=\; \frac{r\,\Gamma\bigl(\tfrac12(r+s)\bigr)}{s\,\Gamma\bigl(\tfrac12 r\bigr)\Gamma\bigl(\tfrac12 s\bigr)} \cdot \frac{(rx/s)^{\frac12 r - 1}}{\bigl[1 + (rx/s)\bigr]^{\frac12 (r+s)}} , \qquad x > 0 .

Grimmett and Stirzaker note two consequences that are immediate from the definition: $F^{-1} \sim F(s,r)$ , by exchanging the roles of $U$ and $V$ ; and $T^{2} \sim F(1,r)$ when $T \sim t(r)$ , because $V^{2} \sim \chi^{2}(1)$ . The mean of $F(r,s)$ is $s/(s-2)$ when $s > 2$ . Riley, Hobson and Bence write the normalising factor through the beta function $B\bigl(\tfrac12 r, \tfrac12 s\bigr) = \Gamma(\tfrac12 r)\Gamma(\tfrac12 s)/\Gamma(\tfrac12(r+s))$ of The Beta function; the two forms are the same.

Definition(The sampling distributions of a normal sample)AdaptedAdapted from Grimmett & Stirzaker (2020), §4.10, Theorem (4.10.1) and the discussion of $U$, $V$, $T$ and $F$ that follows it, pp. 135-137

Let $X_1, \dots, X_n$ be independent, each with the $N(\mu, \sigma^{2})$ distribution, $n \ge 2$ , and let $\overline X$ and $S^{2}$ be the sample mean and the sample variance with divisor $n-1$ . Then:

$\overline X$ and $S^{2}$ are independent random variables;
$\overline X \sim N\bigl(\mu,\ \sigma^{2}/n\bigr)$ ;
$(n-1)S^{2}/\sigma^{2} \sim \chi^{2}(n-1)$ .

Consequently, writing

U = \frac{(n-1)S^{2}}{\sigma^{2}} \sim \chi^{2}(n-1) , \qquad V = \frac{\sqrt n\,(\overline X - \mu)}{\sigma} \sim N(0,1) ,

neither of which depends on $\sigma$ once it is formed, the ratio

T \;=\; \frac{V}{\sqrt{U/(n-1)}} \;=\; \frac{\overline X - \mu}{S/\sqrt n}

has the $t(n-1)$ distribution. The unknown $\sigma$ has cancelled: $T$ is computable from the data and $\mu$ alone, and its distribution involves no unknown parameter.

If, further, $Y_1,\dots,Y_m$ is a second sample, independent of the first, drawn from $N(\nu, \sigma^{2})$ with the same variance, and $S_Y^{2}$ is its sample variance, then $(m-1)S_Y^{2}/\sigma^{2} \sim \chi^{2}(m-1)$ independently of $(n-1)S^{2}/\sigma^{2}$ , so the variance ratio

\frac{S_Y^{2}}{S^{2}} \;=\; \frac{\bigl[(m-1)S_Y^{2}/\sigma^{2}\bigr]/(m-1)}{\bigl[(n-1)S^{2}/\sigma^{2}\bigr]/(n-1)} \ \sim\ F(m-1,\, n-1) ,

again free of $\sigma$ .

Idea of the proofOriginal

Only one idea is needed, and it is geometric. Write $Y_i = (X_i-\mu)/\sigma$ ; then the vector $Y = (Y_1,\dots,Y_n)$ has density proportional to $\exp\bigl(-\tfrac12 \lVert y\rVert^{2}\bigr)$ , which depends on $y$ only through its length. So the law of $Y$ is unchanged by any rotation of $\mathbb{R}^{n}$ , and the coordinates of $Y$ in any orthonormal frame are again independent standard normals.

Now choose the frame to suit the question. Let the first axis be the diagonal direction $n^{-1/2}(1,1,\dots,1)$ . The first coordinate in that frame is $Z_1 = \sqrt n\,\overline Y$ , which carries the sample mean; the remaining $n-1$ coordinates are orthogonal to the diagonal, and because the length is preserved, the sum of their squares is $\lVert Y\rVert^{2} - Z_1^{2} = \sum_i (Y_i - \overline Y)^{2}$ , which is the sample variance up to a factor. So the mean lives on one axis and the variance on the orthogonal complement of that axis: independence of $\overline X$ and $S^{2}$ is the independence of coordinates in a rotated frame, and the $n-1$ degrees of freedom are literally the dimension of the complement.

The $t$ and $F$ laws then cost nothing beyond a change of variables: $t$ is the law of a normal divided by the root of an independent chi-squared, $F$ the law of a ratio of two independent chi-squareds, and both densities come out of the standard two-variable transformation formula followed by one gamma integral.

Proof(Independence, and the three laws)AdaptedAdapted from Grimmett & Stirzaker (2020), §4.10, proof of Theorem (4.10.1), eqs. (2)-(4), p. 136, and the two density calculations, pp. 136-137

Step 1: standardise. Put $Y_i = (X_i - \mu)/\sigma$ . Each $Y_i$ is $N(0,1)$ and they are independent, so the joint density of $Y = (Y_1,\dots,Y_n)$ is

f(y) \;=\; (2\pi)^{-n/2}\exp\Bigl(-\tfrac12 \textstyle\sum_{i=1}^{n} y_i^{2}\Bigr) \;=\; (2\pi)^{-n/2} e^{-\lVert y\rVert^{2}/2} .

Also $\overline Y = (\overline X - \mu)/\sigma$ and

\sum_{i=1}^{n}\bigl(Y_i - \overline Y\bigr)^{2} \;=\; \frac{1}{\sigma^{2}}\sum_{i=1}^{n}\bigl(X_i-\overline X\bigr)^{2} \;=\; \frac{(n-1)S^{2}}{\sigma^{2}} ,

because the shift by $\mu$ cancels inside the bracket and the factor $\sigma$ comes out squared. So the theorem is a statement about a standard normal vector.

Step 2: rotate. Let $A$ be any real orthogonal $n\times n$ matrix and let $Z$ be the coordinates of $Y$ in the rotated frame, so that $Y_i = \sum_{j} Z_j a_{ji}$ and, because $A$ preserves lengths (Unitary and orthogonal matrices), $\sum_i Y_i^{2} = \sum_i Z_i^{2}$ . The Jacobian of an orthogonal map has modulus $1$ , so by the change-of-variables formula the density of $Z$ is $(2\pi)^{-n/2}e^{-\lVert z\rVert^{2}/2}$ : the same function. It factorises into $n$ standard normal densities, hence $Z_1,\dots,Z_n$ are independent $N(0,1)$ variables. This is the only place normality is used, and it is used essentially: no other law on $\mathbb{R}^{n}$ with independent identically distributed coordinates has a rotation-invariant density.

Step 3: choose the rotation. Take the first row of $A$ to be $n^{-1/2}(1,1,\dots,1)$ , a unit vector, and complete it to an orthonormal basis of $\mathbb{R}^{n}$ (possible by Gram-Schmidt; the rows $2,\dots,n$ are then any orthonormal basis of the hyperplane $\sum_i y_i = 0$ ). With this choice

Z_1 \;=\; \frac{1}{\sqrt n}\sum_{i=1}^{n} Y_i \;=\; \sqrt n\, \overline Y ,

and $Z_1 \sim N(0,1)$ by Step 2 - which says exactly that $\overline Y \sim N(0,1/n)$ , that is $\overline X \sim N(\mu, \sigma^{2}/n)$ . That is conclusion 2.

Step 4: identify the rest of the length. Subtracting the first square from the total,

\sum_{j=2}^{n} Z_j^{2} \;=\; \sum_{i=1}^{n} Y_i^{2} - Z_1^{2} \;=\; \sum_{i=1}^{n} Y_i^{2} - n\,\overline Y^{2} \;=\; \sum_{i=1}^{n}\bigl(Y_i - \overline Y\bigr)^{2} \;=\; \frac{(n-1)S^{2}}{\sigma^{2}} ,

the third equality being the identity $\sum_i (Y_i-\overline Y)^{2} = \sum_i Y_i^{2} - n\overline Y^{2}$ , got by expanding the square and using $\sum_i (Y_i - \overline Y) = 0$ .

Step 5: read off independence and the chi-squared law. By Step 2 the variable $Z_1$ is independent of the block $(Z_2,\dots,Z_n)$ . By Step 3, $\overline X$ is a function of $Z_1$ alone; by Step 4, $S^{2}$ is a function of $Z_2,\dots,Z_n$ alone. Functions of independent variables are independent, so $\overline X$ and $S^{2}$ are independent: conclusion 1. And by Step 4, $(n-1)S^{2}/\sigma^{2}$ is a sum of the squares of $n-1$ independent standard normal variables, hence $\chi^{2}(n-1)$ by the additivity of chi-squared recorded in the definition: conclusion 3. The order is $n-1$ , not $n$ , because one linear constraint has been imposed - the residuals $Y_i - \overline Y$ sum to zero and so lie in a hyperplane.

Step 6: the $t$ density. With $U$ and $V$ as in the statement and $r = n-1$ , $U$ and $V$ are independent by conclusion 1, so their joint density is the product

f(u,v) \;=\; \frac{(\tfrac12)^{r/2}}{\Gamma(\tfrac12 r)}\,u^{\frac12 r - 1}e^{-u/2} \cdot \frac{1}{\sqrt{2\pi}}\, e^{-v^{2}/2} , \qquad u > 0 .

Change variables to $(w,t)$ with $w = u$ and $t = v\sqrt{r/u}$ , so $v = t\sqrt{w/r}$ and, at fixed $w$ , $dv = \sqrt{w/r}\,dt$ . The Jacobian is therefore $\sqrt{w/r}$ and

f_{U,T}(w,t) \;=\; \sqrt{w/r}\; f\bigl(w,\, t\sqrt{w/r}\bigr) \;=\; \frac{(\tfrac12)^{r/2}}{\Gamma(\tfrac12 r)\sqrt{2\pi r}}\; w^{\frac12(r-1)}\exp\Bigl[-\frac{w}{2}\Bigl(1+\frac{t^{2}}{r}\Bigr)\Bigr] .

Integrate out $w$ over $(0,\infty)$ . Substituting $\xi = \tfrac{w}{2}(1+t^{2}/r)$ turns the integral into $\Gamma\bigl(\tfrac12(r+1)\bigr)$ times $\bigl[\tfrac12(1+t^{2}/r)\bigr]^{-\frac12(r+1)}$ , by the definition of the gamma function. The powers of $\tfrac12$ do not quite cancel: $(\tfrac12)^{r/2}\,(\tfrac12)^{-\frac12(r+1)} = \sqrt2$ , and this factor turns $\sqrt{2\pi r}$ into $\sqrt{\pi r}$ , leaving

f_T(t) \;=\; \frac{\Gamma\bigl(\tfrac12(r+1)\bigr)}{\sqrt{\pi r}\,\Gamma\bigl(\tfrac12 r\bigr)}\Bigl(1+\frac{t^{2}}{r}\Bigr)^{-\frac12(r+1)} ,

the $t(r)$ density. Finally $T = V/\sqrt{U/(n-1)}$ equals $\sqrt n(\overline X-\mu)/\sigma$ divided by $S/\sigma$ , and the $\sigma$ 's cancel, giving $T = (\overline X - \mu)/(S/\sqrt n)$ as claimed.

Step 7: the $F$ density. Let $U\sim\chi^{2}(r)$ and $V\sim\chi^{2}(s)$ be independent and set $x = (U/r)/(V/s)$ . The same two-variable device works: keep $w = V$ as the second coordinate, so $U = (r/s)\,x\,w$ with Jacobian $\partial U/\partial x = (r/s)w$ , multiply the two gamma densities, and integrate $w$ out over $(0,\infty)$ . The $w$ -integral is again a gamma integral, of total order $\tfrac12(r+s)$ , and after collecting constants one obtains

f_F(x) \;=\; \frac{r\,\Gamma\bigl(\tfrac12(r+s)\bigr)}{s\,\Gamma\bigl(\tfrac12 r\bigr)\Gamma\bigl(\tfrac12 s\bigr)} \cdot\frac{(rx/s)^{\frac12 r-1}}{\bigl[1+(rx/s)\bigr]^{\frac12(r+s)}} , \qquad x>0 .

Applying this with $U = (m-1)S_Y^{2}/\sigma^{2}$ and $V = (n-1)S^{2}/\sigma^{2}$ , which are independent because the two samples are, gives the last claim of the statement. $\square$

Remark(Where the hypotheses are used)Original

Normality. Used once, in Step 2, and it is what the whole theorem rests on. Rotation invariance alone does not give independent coordinates (the uniform law on a disc is rotation-invariant, yet its two coordinates are dependent); the fact Step 2 rests on is narrower. A vector of independent standard normals has a rotation-invariant density, and among laws on $\mathbb{R}^{n}$ with independent, identically distributed coordinates and a density, only a centred normal has one. For any other parent law, $\overline X$ and $S^{2}$ are in general dependent. A two-point example makes this concrete: let $n=2$ and let each $X_i$ be $\pm 1$ with probability $\tfrac12$ . Then $\overline X \in \{-1,0,1\}$ and $S^{2} = (X_1-X_2)^{2}/2 \in \{0,2\}$ , and $S^{2} = 0$ holds exactly when $\overline X = \pm 1$ . Knowing one determines the other, so they are as far from independent as possible, although $\mathbb{E}(S^{2}) = \sigma^{2} = 1$ still holds.

Independence of the $X_i$ . Used in Step 1 to write the joint density as a product. Without it there is no rotation-invariant density to exploit. Correlated normal data can be brought back into the scheme by whitening - multiply by $\Sigma^{-1/2}$ , with $\Sigma$ the covariance matrix, as in The multivariate normal distribution - but then the statistic being computed is not the ordinary $S^{2}$ .

Common variance $\sigma^{2}$ for both samples, in the $F$ statement. If the two populations had variances $\sigma_1^{2} \ne \sigma_2^{2}$ , the ratio $S_Y^{2}/S^{2}$ would be $(\sigma_2^{2}/\sigma_1^{2})$ times an $F(m-1,n-1)$ variable, so its distribution would involve the unknown ratio and could not be used without further assumption. This is precisely why the ratio is informative about that ratio.

$n \ge 2$ . For $n = 1$ there is no $S^{2}$ , and the complement of the diagonal axis in $\mathbb{R}^{1}$ is $\{0\}$ ; the theorem degenerates correctly rather than failing.

What survives without normality. Only $\mathbb{E}(\overline X) = \mu$ and $\operatorname{var}(\overline X) = \sigma^{2}/n$ , and - as a limit, not an identity - the approximate normality of $\overline X$ for large $n$ given by the central limit theorem. The exact $t$ and $F$ laws are properties of the normal family alone, though they are mildly robust in practice.

A convention to watch. Riley, Hobson and Bence define the sample variance with divisor $N$ and accordingly write Student's statistic as $t = (\bar x - \mu_0)/\bigl(s/\sqrt{N-1}\bigr)$ , with $N-1$ under the root. Since $(N-1)S^{2} = Ns^{2}$ , one has $s/\sqrt{N-1} = S/\sqrt N$ , so their $t$ and ours are the same number computed from the same data; the degrees of freedom are $N-1$ in both. Only the intermediate symbol $s$ differs. The same substitution turns their $u = Ns^{2}/\sigma_0^{2} \sim \chi^{2}_{N-1}$ into our conclusion 3.

Method(Using the three distributions)Original

Given a sample $x_1,\dots,x_n$ believed to come from a normal population:

Compute $\bar x = n^{-1}\sum_i x_i$ and the sum of squared deviations $Q = \sum_i (x_i - \bar x)^{2}$ ; then $S^{2} = Q/(n-1)$ and $S = \sqrt{S^{2}}$ .
If $\sigma$ is known, use $\sqrt n\,(\bar x - \mu)/\sigma \sim N(0,1)$ .
If $\sigma$ is unknown, use $t = (\bar x - \mu)/(S/\sqrt n) \sim t(n-1)$ . This is the only change: the standard normal is replaced by $t$ with $n-1$ degrees of freedom, whose tails are heavier, the price of not knowing $\sigma$ .
For a statement about $\sigma^{2}$ itself, use $U = (n-1)S^{2}/\sigma^{2} \sim \chi^{2}(n-1)$ .
To compare the spreads of two independent normal samples, of sizes $n_1$ and $n_2$ , use $F = S_2^{2}/S_1^{2} \sim F(n_2-1,\,n_1-1)$ under the hypothesis of equal variances. Since $F^{-1} \sim F(n_1-1,\,n_2-1)$ , the ratio may be taken either way up, but the orientation must be fixed before the data are seen. For a one-sided question ("is sample 2 the noisier?") read the single tail on that side. For a two-sided question ("do the spreads differ?") the tail probability is twice the smaller of $\Pr(F > f)$ and $\Pr(F < f)$ at the observed value $f$ . Putting whichever sample variance happens to be larger on top and then reading one upper tail undoubled roughly halves this figure, because the orientation has been chosen from the data.
Convert the resulting tail probability into an interval or a test with Confidence intervals or Hypothesis testing; this unit supplies the distributions, not the decisions.

Example(A calibration run: mean, variance and a variance ratio)Original

A platinum resistance thermometer is held in a stirred water reference bath and read eight times. The readings, in degrees Celsius, are

20.12,\quad 19.87,\quad 20.31,\quad 20.05,\quad 19.76,\quad 20.44,\quad 19.95,\quad 20.18 .

Treat them as a random sample from $N(\mu,\sigma^{2})$ .

Sample mean. The readings sum to $160.68$ , so with $n=8$

\bar x = \frac{160.68}{8} = 20.085\ ^{\circ}\mathrm{C}.

Sample variance. The deviations from $20.085$ are $0.035,\,-0.215,\,0.225,\,-0.035,\,-0.325,\,0.355,\,-0.135,\,0.095$ , and they sum to zero, as they must. Their squares sum to

Q = \sum_{i=1}^{8}(x_i-\bar x)^{2} = 0.3582 ,

S^{2} = \frac{Q}{n-1} = \frac{0.3582}{7} = 0.05117\ (^{\circ}\mathrm{C})^{2} , \qquad S = 0.22621\ ^{\circ}\mathrm{C} .

For comparison, the Riley-Hobson-Bence convention would give $s^{2} = Q/8 = 0.044775$ and $s = 0.21160$ , smaller by the factor $\sqrt{7/8}$ .

Distance of the mean from a reference value. The bath is believed to sit at $\mu_0 = 20.000\,^{\circ}\mathrm{C}$ . Since $\sigma$ is unknown we use the $t$ statistic of step 3:

t = \frac{\bar x - \mu_0}{S/\sqrt n} = \frac{20.085 - 20.000}{0.22621/\sqrt 8} = \frac{0.085}{0.079977} = 1.0628 ,

to be read against $t(7)$ . The two-sided tail probability is $\Pr\bigl(\lvert T\rvert > 1.0628\bigr) = 0.3232$ for $T \sim t(7)$ , so a departure at least this large happens in about a third of samples of this size even when the bath is exactly at $20\,^{\circ}\mathrm{C}$ . Had $\sigma$ been treated as known and equal to $S$ , the same ratio would have been read against $N(0,1)$ , giving the smaller tail $0.2879$ ; the difference is the penalty for estimating $\sigma$ from seven degrees of freedom.

A statement about the spread. The manufacturer quotes a repeatability of $\sigma_0 = 0.25\,^{\circ}\mathrm{C}$ . Under that value, step 4 gives

U = \frac{(n-1)S^{2}}{\sigma_0^{2}} = \frac{0.3582}{0.0625} = 5.7312 ,

a draw from $\chi^{2}(7)$ , whose mean is $7$ and standard deviation $\sqrt{14} = 3.742$ . The upper tail is $\Pr\bigl(\chi^{2}(7) > 5.7312\bigr) = 0.5715$ : the observed scatter is entirely ordinary for the quoted repeatability.

A second sample and the variance ratio. The bath is re-read six times the next day, after the probe has been re-seated:

20.31,\quad 19.62,\quad 20.55,\quad 19.88,\quad 20.44,\quad 19.70 ,

with $\bar y = 120.50/6 = 20.0833$ and $\sum_j (y_j-\bar y)^{2} = 0.79933$ , hence $S_Y^{2} = 0.79933/5 = 0.159867$ . The question is two-sided - did the spread change between the days? - and we fix the orientation as day 2 over day 1:

F = \frac{S_Y^{2}}{S^{2}} = \frac{0.159867}{0.051171} = 3.1241 ,

a draw from $F(5,7)$ if the two days share one $\sigma$ . The upper tail is $\Pr\bigl(F(5,7) > 3.1241\bigr) = 0.0848$ and the lower tail is $0.9152$ , so by step 5 the two-sided tail probability is $2 \times 0.0848 = 0.1696$ ; for a two-sided test at the $5\%$ level the ratio would have to exceed the $97.5\%$ point of $F(5,7)$ , which is $5.2852$ . The second day is noisier, but not by more than the ratio of two sample variances of these sizes routinely produces. The single tail $0.0848$ would be the right figure only for the one-sided question "did re-seating make the probe noisier?", posed before the second day's readings were taken.

Every number above was computed to the digits shown; the three tail probabilities are $0.3232$ , $0.5715$ and $0.1696$ .

Application(Error bars, instrument qualification, and the analysis of variance)Original

Quoting an error bar. The statement " $\bar x \pm S/\sqrt n$ with $n-1$ degrees of freedom" is the normal-sample result above, and it is what a laboratory writes down after repeating a measurement. The $\sqrt n$ is conclusion 2; the $n-1$ is conclusion 3; the fact that one may divide by the estimated $\sigma$ and still know the exact distribution is conclusion 1 together with the $t$ law. For $n$ of order ten the difference between the $t$ and normal tails is tens of percent, which is why small-sample work quotes degrees of freedom and not just a standard error.

Qualifying an instrument. A detector specification of the form "resolution better than $\sigma_0$ " is checked by reading $(n-1)S^{2}/\sigma_0^{2}$ against $\chi^{2}(n-1)$ , exactly as in the example. The same chi-squared appears, with $n-M$ in place of $n-1$ , as the goodness-of-fit statistic of a least-squares fit with $M$ fitted parameters: each fitted parameter is one more linear constraint on the residual vector, one more dimension removed from the orthogonal complement, in exactly the sense of Step 4.

Comparing two populations. The $F$ ratio is the engine of the analysis of variance, in which a total sum of squares is split into orthogonal pieces - between groups and within groups - whose expectations differ if and only if the group means differ. The splitting is the same rotation argument applied to a coarser partition of $\mathbb{R}^{n}$ , and the independence of the pieces, which is what makes the ratio an $F$ variable, is again the independence of coordinates in a rotated frame.

Spectral estimation. A periodogram ordinate formed from Gaussian noise is a sum of two squared normals, hence $\sigma^{2}\chi^{2}(2)/2$ in distribution, which is why a raw periodogram has a $100\%$ relative error however long the record, and why averaging independent segments - adding degrees of freedom - is the only cure. The same counting gives the degrees of freedom of a smoothed spectral estimate.

Definition

Estimators and their properties

Adapted

Add to PDF

Defines bias, consistency, efficiency and the standard error of an estimator, introduces the score and the Fisher information of a sample, and proves the Cramer-Rao (Fisher) inequality bounding the variance of any estimator from below; worked on twelve Poisson counts, where the sample mean attains the bound $\hat\lambda = 6.833 \pm 0.755$ while the sample variance reaches only $6.3\%$ of it.

Prerequisites: Sampling distributions, Expectation, variance and covariance, Markov, Chebyshev, Jensen and Chernoff inequalities, Convergence of random variables, The central limit theorem, Inner products and norms, Quadratic forms and the Rayleigh quotient, Differentiation under the integral sign

Note(What this unit is for)Original

A measurement does not hand us a parameter; it hands us a sample, from which we compute a number and call that number our estimate. Several different numbers will usually be available for the same parameter - the mean and the median both estimate the centre of a symmetric law - and they are not equally good. This unit sets out the four things one asks of such a rule: does it converge to the truth, is it right on average, how much does it scatter, and how much scatter was unavoidable.

The last question has a sharp answer. The data carry a finite amount of information about the parameter, measurable as a single number, the Fisher information, and no estimator whatever can scatter less than that number allows. That is the Cramer-Rao inequality, and it converts "this estimator looks good" into "this estimator is within a known factor of the best possible". The distributions that the calculations below run on are those of Sampling distributions; the two standard recipes for building estimators are Maximum likelihood and Least-squares fitting, and turning an estimate into a statement with a stated probability is the business of Confidence intervals.

Definition(Estimator, bias, consistency)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.3, pp. 1229-1230, and §31.3.1, eqs. (31.14)-(31.15), pp. 1230-1231

Let the sample $x = (x_1,\dots,x_n)$ be drawn from a population whose joint density $P(x\mid a)$ is known except for the value of a parameter $a$ . (Riley, Hobson and Bence write $N$ for the sample size and allow $a$ to stand for a whole list $a_1,\dots,a_m$ of unknown parameters; we write $n$ for the sample size and treat one parameter at a time until the Fisher matrix below.)

An estimator of $a$ is a statistic $\hat a(x)$ - a function of the sample alone, in the sense of Sampling distributions - whose value is offered as an estimate of $a$ . Being a function of random variables it is itself random, and its law is its sampling distribution $P(\hat a \mid a)$ .

Its expectation and variance may be computed either over the sampling distribution or over the sample space,

\mathbb{E}[\hat a] = \int \hat a\, P(\hat a\mid a)\,d\hat a = \int \hat a(x)\, P(x\mid a)\,d^{n}x ,

\operatorname{var}[\hat a] = \int \bigl(\hat a(x) - \mathbb{E}[\hat a]\bigr)^{2} P(x\mid a)\,d^{n}x ,

the two integrals agreeing by the law of the unconscious statistician of Expectation, variance and covariance.

The bias of $\hat a$ is

b(a) \;=\; \mathbb{E}[\hat a] - a .

It is a function of the true $a$ and of $n$ , never of the observed data. If $b \equiv 0$ the estimator is unbiased.

The estimator is consistent when, as the sample is made larger and larger, it settles down on the parameter it is estimating. Riley, Hobson and Bence write the requirement as $\hat a \to a$ in the limit $n \to \infty$ , and give alongside it the equivalent demand that the sampling distribution should collapse onto a point mass there, $P(\hat a \mid a) \to \delta(\hat a - a)$ . In the vocabulary of Convergence of random variables the statement being made is convergence in probability: for every $\varepsilon > 0$ , $\Pr(\lvert \hat a - a\rvert > \varepsilon) \to 0$ , which is the sense in which the law of large numbers holds.

Definition(Standard error, mean squared error, efficiency)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.3.1, eq. (31.16) and the discussion of efficiency, pp. 1231-1232, and §31.3.3, eq. (31.24) and the root mean square error, pp. 1234-1235

The standard error on an estimate is the standard deviation of its sampling distribution,

\sigma_{\hat a} \;=\; \bigl(\operatorname{var}[\hat a]\bigr)^{1/2} .

It usually depends on the unknown $a$ ; one then substitutes the estimate into the formula and quotes the estimated standard error $\hat\sigma_{\hat a}$ , writing the result as $a = \hat a \pm \hat\sigma_{\hat a}$ .

For a biased estimator the spread about the true value is measured instead by the root mean square error $\epsilon_{\hat a}$ , where

\epsilon_{\hat a}^{2} \;=\; \mathbb{E}\bigl[(\hat a - a)^{2}\bigr] \;=\; \operatorname{var}[\hat a] + b(a)^{2} ,

the decomposition following by adding and subtracting $\mathbb{E}[\hat a]$ inside the square, since the cross term has expectation zero. The two pieces are the statistical and the systematic part of the error.

Of two estimators of the same $a$ , the one with the smaller variance is the more efficient. Writing $V_{\min}$ for the lower bound established below, the efficiency of $\hat a$ is

e \;=\; \frac{V_{\min}}{\operatorname{var}[\hat a]} \;\in\; (0,1] ,

and $\hat a$ is efficient, or minimum-variance, when $e = 1$ . Among unbiased estimators the one of smallest variance is the best unbiased estimator. Riley, Hobson and Bence stress that no single estimator is optimal in all respects: there is a trade-off between bias and variance, and a slightly biased estimator with a small variance may well be preferred to an unbiased one with a large variance. The root mean square error is the quantity that weighs the two against each other.

Definition(Score and Fisher information)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.3.2, eqs. (31.18)-(31.23), pp. 1233-1234, and §31.3.5, eq. (31.36), p. 1241

Suppose $P(x \mid a) > 0$ on a region of sample space that does not depend on $a$ , and that $P$ is twice differentiable in $a$ there. The score of the sample is the logarithmic derivative

U(x,a) \;=\; \frac{\partial \ln P(x\mid a)}{\partial a} ,

and the Fisher information carried by the sample about $a$ is

I(a) \;=\; \mathbb{E}\bigl[U^{2}\bigr] \;=\; \int \left(\frac{\partial \ln P}{\partial a}\right)^{2} P(x\mid a)\, d^{n}x .

Step 2 of the proof below shows that $\mathbb{E}[U] = 0$ , so $I(a)$ is also the variance of the score, and Step 6 shows that under one further differentiation

I(a) \;=\; \mathbb{E}\!\left[-\frac{\partial^{2}\ln P}{\partial a^{2}}\right] ,

which is the form Riley, Hobson and Bence use: the information is the average curvature of the log-density at the true parameter. For several parameters $a_1,\dots,a_M$ the same construction gives the Fisher matrix

F_{ij} \;=\; \mathbb{E}\!\left[-\frac{\partial^{2}\ln P}{\partial a_i\,\partial a_j}\right] ,

a symmetric positive semi-definite matrix (a quadratic form in the sense of Quadratic forms and the Rayleigh quotient, since $v^{\mathsf T} F v = \mathbb{E}[(\sum_i v_i \partial_i \ln P)^{2}] \ge 0$ ). Riley, Hobson and Bence state that for unbiased efficient estimators the Fisher matrix is exactly the inverse of their covariance matrix, $(V^{-1})_{ij} = F_{ij}$ , which is the many-parameter form of the bound proved below.

Information adds over independent observations. If the sample is independent with common one-observation density $p(x_i\mid a)$ , then $\ln P = \sum_i \ln p(x_i\mid a)$ , so the score is a sum of $n$ independent terms each of mean zero, and variances of independent summands add (Expectation, variance and covariance). Hence

I(a) \;=\; n\, i(a) , \qquad i(a) = \mathbb{E}\!\left[\left(\frac{\partial \ln p(x_1\mid a)}{\partial a}\right)^{2}\right] ,

with $i(a)$ the information in a single observation. Every $1/n$ in the error bars below comes from this one line.

Definition(The Cramer-Rao inequality)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.3.1, eq. (31.17), p. 1232, and §31.3.2, eqs. (31.18) and (31.23), pp. 1233-1234

Let $\hat a(x)$ be any estimator of $a$ with bias $b(a) = \mathbb{E}[\hat a] - a$ , and suppose that

the set where $P(x\mid a) > 0$ does not depend on $a$ , and $P$ is differentiable in $a$ there;
the integrals $\int P\,d^{n}x$ and $\int \hat a\, P \,d^{n}x$ may be differentiated with respect to $a$ under the integral sign (Differentiation under the integral sign);
$0 < I(a) < \infty$ and $\operatorname{var}[\hat a] < \infty$ .

Then

\operatorname{var}[\hat a] \;\ge\; V_{\min} \;=\; \frac{\bigl(1 + b'(a)\bigr)^{2}}{I(a)} \;=\; \frac{\bigl(1 + \partial b/\partial a\bigr)^{2}}{\mathbb{E}\bigl[(\partial \ln P/\partial a)^{2}\bigr]} ,

and, if in addition $P$ is twice differentiable in $a$ and $\int P\,d^{n}x$ may be differentiated twice under the integral sign,

\operatorname{var}[\hat a] \;\ge\; \frac{\bigl(1 + \partial b/\partial a\bigr)^{2}}{\mathbb{E}\bigl[-\partial^{2}\ln P/\partial a^{2}\bigr]} .

In particular an unbiased estimator satisfies $\operatorname{var}[\hat a] \ge 1/I(a)$ . Equality holds if and only if the score is an affine function of the estimator,

\frac{\partial \ln P}{\partial a} \;=\; c(a)\,\bigl[\hat a(x) - \mathbb{E}[\hat a]\bigr]

for almost every $x$ , with $c(a)$ depending on $a$ but not on $x$ .

Riley, Hobson and Bence call this Fisher's inequality and also give the name Cramer-Rao inequality; they write it with $P$ the population density of the whole sample, as here.

Idea of the proofOriginal

Think of random variables with finite variance as vectors, with $(F,G) = \mathbb{E}[FG]$ as the inner product. Two vectors matter: the centred estimator $\hat a - \mathbb{E}[\hat a]$ , whose length squared is the variance we want to bound below, and the score $U$ , whose length squared is the information.

Differentiating the normalisation $\int P \,d^{n}x = 1$ says that the score has mean zero, so it is already centred. Differentiating $\int \hat a P\,d^{n}x = a + b(a)$ says that the inner product of the two vectors is fixed at $1 + b'(a)$ , whatever the estimator. So we have two vectors with a prescribed inner product, and the Cauchy-Schwarz inequality of Inner products and norms immediately forbids either of them from being too short: $\lvert(\hat a - \mathbb{E}[\hat a], U)\rvert \le \lVert \hat a - \mathbb{E}[\hat a]\rVert\,\lVert U \rVert$ rearranges into exactly the bound. Equality in Cauchy-Schwarz means the two vectors are parallel, which is the condition stated for an efficient estimator.

So the content of the theorem is a triangle: the estimator must "overlap" the score by a fixed amount, the score has a fixed length, and therefore the estimator cannot be short. Nothing about the estimator enters except that it is an unbiased-up-to- $b$ function of the data.

Proof(The bound, and both forms of the information)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.3.2, the derivation of (31.19)-(31.23) and the identity between the two forms of the information, pp. 1233-1234

Write $\partial$ for $\partial/\partial a$ and suppress the arguments of $P$ .

Step 1: the two identities we shall need. Everything comes from differentiating two integrals in $a$ . Hypothesis 2 licenses taking $\partial$ inside, and hypothesis 1 guarantees that no boundary term appears, because the region of integration where $P>0$ is the same for every $a$ .

Step 2: the score has mean zero. Differentiating $\int P \,d^{n}x = 1$ ,

0 \;=\; \partial \!\int P\,d^{n}x \;=\; \int \partial P \, d^{n}x \;=\; \int (\partial \ln P)\, P\, d^{n}x \;=\; \mathbb{E}[U] ,

where the third equality is the chain rule $\partial P = (\partial \ln P) P$ , valid where $P>0$ . Consequently $I(a) = \mathbb{E}[U^{2}] = \operatorname{var}[U]$ .

Step 3: the score overlaps the estimator by $1+b'$ . The definition of the bias gives $\int \hat a(x) P \,d^{n}x = a + b(a)$ . Differentiating both sides, and using the chain rule again on the left ( $\hat a$ does not depend on $a$ ),

\int \hat a\, (\partial \ln P)\, P\, d^{n}x \;=\; 1 + b'(a) , \qquad\text{that is}\qquad \mathbb{E}[\hat a\,U] \;=\; 1 + b'(a) .

Step 4: centre the estimator. Multiply the identity of Step 2 by the constant $\mathbb{E}[\hat a]$ and subtract it from Step 3. Since $\mathbb{E}[\hat a]$ does not depend on $x$ , this gives

\mathbb{E}\bigl[(\hat a - \mathbb{E}[\hat a])\,U\bigr] \;=\; 1 + b'(a) .

Subtracting any other constant would have produced a larger left-hand factor in the next step; subtracting the mean is what makes the bound sharp, and this is the point at which Riley, Hobson and Bence choose their free function $\alpha(a)$ to be $\mathbb{E}[\hat a]$ .

Step 5: Cauchy-Schwarz. On the space of functions of $x$ with $\int (\cdot)^{2} P \,d^{n}x < \infty$ , the pairing $(F,G) = \int FG\,P\,d^{n}x$ is an inner product, so by Inner products and norms,

\Bigl(\mathbb{E}\bigl[(\hat a - \mathbb{E}[\hat a])\,U\bigr]\Bigr)^{2} \;\le\; \mathbb{E}\bigl[(\hat a - \mathbb{E}[\hat a])^{2}\bigr]\;\mathbb{E}\bigl[U^{2}\bigr] \;=\; \operatorname{var}[\hat a]\; I(a) .

The left-hand side is $(1+b')^{2}$ by Step 4, and $I(a) > 0$ by hypothesis 3, so dividing gives

\operatorname{var}[\hat a] \;\ge\; \frac{(1+b'(a))^{2}}{I(a)} ,

which is the first form. Equality in Cauchy-Schwarz holds exactly when the two entries are proportional, $U = c(a)(\hat a - \mathbb{E}[\hat a])$ , which is the stated equality condition.

Step 6: the curvature form of the information. Differentiate the identity $\int (\partial \ln P)\,P \,d^{n}x = 0$ of Step 2 once more in $a$ , again taking the derivative inside:

\int \bigl(\partial^{2}\ln P\bigr) P \,d^{n}x \;+\; \int \bigl(\partial \ln P\bigr)\,\partial P \,d^{n}x \;=\; 0 .

In the second integral write $\partial P = (\partial \ln P)P$ once more; it becomes $\int (\partial\ln P)^{2} P\,d^{n}x = I(a)$ . Hence

I(a) \;=\; -\int \bigl(\partial^{2}\ln P\bigr)P\,d^{n}x \;=\; \mathbb{E}\!\left[-\frac{\partial^{2}\ln P}{\partial a^{2}}\right] ,

and substituting this into the first form gives the second. $\square$

Remark(Where the hypotheses are used)Original

Fixed support (hypothesis 1). Used in Steps 2, 3 and 6, where the derivative is taken inside the integral with no boundary term. It is not a formality. Let $x_1,\dots,x_n$ be uniform on $(0,a)$ . Then $\ln p = -\ln a$ wherever $p>0$ , so $\partial \ln p = -1/a$ and $i(a) = 1/a^{2}$ , giving the apparent bound $a^{2}/n$ for an unbiased estimator. But $M = \max_i x_i$ has $\mathbb{E}[M] = na/(n+1)$ , so $\hat a = (n+1)M/n$ is unbiased with

\operatorname{var}[\hat a] = \frac{a^{2}}{n(n+2)} ,

smaller than $a^{2}/n$ by a factor $n+2$ . The inequality fails outright, because the region $\{0<x_i<a\}$ moves with $a$ and the differentiations of Steps 2 and 3 pick up boundary terms that were discarded.

Differentiation under the integral (hypothesis 2). Used at the same three places. Differentiation under the integral sign gives the standard sufficient condition - continuity of $P$ and $\partial P$ - which holds for the exponential families met in practice.

Finite, non-zero information (hypothesis 3). $I(a) > 0$ is needed to divide in Step 5; if $I(a) = 0$ the data say nothing about $a$ and the bound reads $\operatorname{var} \ge \infty$ unless $b' = -1$ . $I(a)<\infty$ and $\operatorname{var}[\hat a]<\infty$ are needed for Cauchy-Schwarz to be meaningful.

Twice differentiable (for the second form only). The curvature form is the one usually evaluated, because the second derivative of a log-likelihood is often free of the data or easy to average, but it is a corollary of the first form, not an independent statement.

What the bound does not say. It bounds the variance, not the mean squared error. A biased estimator may have variance below $1/I(a)$ - the factor $(1+b')^{2}$ is then smaller than one - and the example below exhibits one. It also does not assert that an estimator attaining the bound exists; it exists only when the score happens to be affine in some statistic, which is the exponential-family case.

The limit in "consistent". The definition reported above writes (Riley, Hobson & Bence (2006), §31.3.1, p. 1231) $\lim_{n\to\infty}\hat a = a$ . Since $\hat a$ is a random variable this is not an ordinary limit, and the statement intended is the one made precise in Convergence of random variables: convergence in probability, equivalently the collapse of $P(\hat a\mid a)$ onto a point mass, which is the second form they give. The accompanying argument of Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §31.3.1, p. 1232) is sound in that reading: if $\hat a$ is unbiased and $\operatorname{var}[\hat a]\to 0$ , then for any $\varepsilon>0$ Chebyshev's inequality of Markov, Chebyshev, Jensen and Chernoff inequalities gives $\Pr(\lvert\hat a - a\rvert>\varepsilon) \le \operatorname{var}[\hat a]/\varepsilon^{2}\to 0$ . The converse fails: a consistent estimator need not be unbiased at any finite $n$ .

Method(Assessing an estimator)Original

Given a candidate $\hat a(x)$ and a parametric model $P(x\mid a)$ :

Bias. Compute $\mathbb{E}[\hat a]$ and subtract $a$ . For a sum or an average of sample values this needs only linearity of expectation; for a ratio or a square root it usually does not close in elementary form and one expands to order $1/n$ .
Variance. Compute $\operatorname{var}[\hat a]$ , if necessary in terms of the central moments $\nu_r$ of the parent. For the corrected sample variance $\hat\sigma^{2} = (n-1)^{-1}\sum_i (x_i-\bar x)^{2}$ , Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §31.4.2, eqs. (31.47)-(31.49), pp. 1247-1248) give

\operatorname{var}\bigl[\hat\sigma^{2}\bigr] \;=\; \frac{1}{n}\left(\nu_4 - \frac{n-3}{n-1}\,\nu_2^{2}\right) ,

which needs the fourth central moment of the population, not just its variance. 3. Consistency. If the estimator is unbiased and its variance tends to zero, it is consistent, by Chebyshev. Otherwise check that both the bias and the variance tend to zero, which gives $\epsilon_{\hat a}\to 0$ and hence consistency by the same inequality. 4. Information. Write down $\ln P(x\mid a)$ , differentiate twice in $a$ , change sign and average: that is $I(a)$ . For an independent sample use $I = n\,i$ and do the work on one observation. 5. Efficiency. Form $V_{\min} = (1+b')^{2}/I(a)$ and divide by the variance from step 2. If the score is proportional to $\hat a - \mathbb{E}[\hat a]$ the estimator is efficient and the division is unnecessary. 6. Quote. For an unbiased estimator quote $\hat a \pm \hat\sigma_{\hat a}$ , substituting the estimate wherever the unknown $a$ appears in the variance formula. For a biased one quote $\hat a \pm \hat\epsilon_{\hat a}$ , or state the bias separately as a systematic error.

Example(Twelve minutes of muon counts: three estimators of one rate)Original

A scintillator telescope records cosmic-ray muons in twelve consecutive one-minute gates. The counts are

7,\quad 4,\quad 9,\quad 6,\quad 5,\quad 11,\quad 8,\quad 6,\quad 3,\quad 7,\quad 10,\quad 6 ,

with total $\sum_i x_i = 82$ and $n = 12$ . Model them as independent Poisson variables of unknown mean $\lambda$ , so that

P(x\mid\lambda) \;=\; \prod_{i=1}^{12} \frac{\lambda^{x_i} e^{-\lambda}}{x_i!} , \qquad \ln P \;=\; -n\lambda + \Bigl(\textstyle\sum_i x_i\Bigr)\ln\lambda - \textstyle\sum_i \ln (x_i!) .

The information. Differentiating twice, $\partial^{2}\ln P/\partial\lambda^{2} = -\lambda^{-2}\sum_i x_i$ , whose expectation is $-n\lambda/\lambda^{2}$ . Hence

I(\lambda) \;=\; \frac{n}{\lambda} , \qquad\text{so}\qquad V_{\min} = \frac{\lambda}{n} \ \text{ for any unbiased estimator.}

Estimator 1: the sample mean. $\hat\lambda_1 = \bar x = 82/12 = 6.8333$ . It is unbiased, since $\mathbb{E}[x_i]=\lambda$ , and $\operatorname{var}[\bar x] = \lambda/n$ because a Poisson variable has variance equal to its mean. That is exactly $V_{\min}$ , so $\hat\lambda_1$ is efficient, $e=1$ . One can see this from the equality condition directly:

\frac{\partial \ln P}{\partial\lambda} \;=\; -n + \frac{1}{\lambda}\sum_i x_i \;=\; \frac{n}{\lambda}\,\bigl(\bar x - \lambda\bigr) ,

which is $c(\lambda)$ times $\hat\lambda_1 - \mathbb{E}[\hat\lambda_1]$ with $c(\lambda)=n/\lambda$ . Substituting the estimate into the standard error $\sqrt{\lambda/n}$ gives

\lambda \;=\; 6.833 \pm 0.755\ \text{counts per minute} .

Estimator 2: the sample variance. A Poisson law has variance equal to its mean, so $\hat\lambda_2 = \hat\sigma^{2} = \frac{1}{11}\sum_i (x_i-\bar x)^{2}$ is also an unbiased estimator of $\lambda$ . The deviations from $6.8333$ have squares summing to $\sum_i (x_i-\bar x)^{2} = 185/3 = 61.667$ , so

\hat\lambda_2 \;=\; \frac{61.667}{11} \;=\; 5.606 .

Its variance follows from the formula in step 2 of the method with the Poisson central moments $\nu_2 = \lambda$ and $\nu_4 = \lambda + 3\lambda^{2}$ :

\operatorname{var}[\hat\lambda_2] = \frac{1}{n}\left(\lambda + 3\lambda^{2} - \frac{n-3}{n-1}\lambda^{2}\right) = \frac{\lambda}{n} + \frac{2\lambda^{2}}{n-1} ,

since $3 - (n-3)/(n-1) = 2n/(n-1)$ . With $n=12$ and $\lambda \approx 6.8333$ this is $0.569 + 8.490 = 9.059$ , a standard error of $3.010$ against $0.755$ . The efficiency is

e \;=\; \frac{\lambda/n}{\lambda/n + 2\lambda^{2}/(n-1)} \;=\; \frac{1}{1 + 2n\lambda/(n-1)} \;=\; \frac{1}{1 + 164/11} \;=\; 0.0629 .

Both estimators are unbiased and both are consistent, but the second throws away about $94\%$ of the information. Estimating a Poisson rate from the scatter rather than from the total is a quantitatively bad idea, and the bound says so without any need to simulate.

Estimator 3: a biased one that beats the unbiased bound. Take $\hat\lambda_3 = \frac{1}{n+1}\sum_i x_i = 82/13 = 6.3077$ . Then $\mathbb{E}[\hat\lambda_3] = n\lambda/(n+1)$ , so

b(\lambda) = -\frac{\lambda}{n+1} , \qquad b'(\lambda) = -\frac{1}{n+1} , \qquad 1 + b' = \frac{n}{n+1} ,

and the bound becomes $V_{\min} = (n/(n+1))^{2}\lambda/n = n\lambda/(n+1)^{2}$ . The actual variance is $\operatorname{var}[\hat\lambda_3] = n\lambda/(n+1)^{2}$ as well: this estimator is also efficient, and with $\lambda\approx 6.8333$ its variance is $0.485$ , below the unbiased bound $0.569$ . There is no contradiction, because it is biased. Comparing on the root mean square error,

\epsilon^{2}[\hat\lambda_3] = \frac{n\lambda + \lambda^{2}}{(n+1)^{2}} = \frac{82 + 46.694}{169} = 0.762 , \qquad \epsilon^{2}[\hat\lambda_1] = \frac{\lambda}{n} = 0.569 ,

so here the shrunken estimator is worse overall. Solving $\epsilon^{2}[\hat\lambda_3] < \epsilon^{2}[\hat\lambda_1]$ gives $\lambda < (2n+1)/n = 2.083$ : shrinking towards zero pays only when the rate is low. At $6.8$ counts per gate it does not.

The three numbers to check are $6.833$ , $5.606$ and $6.308$ , with efficiencies $1$ , $0.0629$ and $1$ .

Application(Error budgets, detector design, and Fisher forecasts)Original

Quoting a result. Every " $\pm$ " in an experimental paper is a standard error in the sense defined here, and the distinction between it and the root mean square error is the distinction between the statistical and the systematic parts of an error budget. A calibration offset that does not average away is a bias $b$ ; the inequality is silent about it, which is why systematics have to be controlled rather than beaten down with more data.

Designing the measurement before taking it. Because $I(a) = n\,i(a)$ and $V_{\min} = 1/I$ , the Cramer-Rao bound predicts the achievable error bar from the model alone, before any data exist. Choosing between two experimental configurations then means comparing their one-observation informations $i(a)$ : how many hours of beam time, how many pixels, what exposure. In cosmology this calculation is standard enough to have a name - a Fisher forecast - and it is the Fisher matrix $F_{ij}$ , not the scalar, that is computed, its inverse giving the forecast covariance ellipse for the parameters, as in Confidence intervals.

Resolution limits. In signal processing the same bound fixes how precisely a time delay, a frequency or a direction of arrival can be recovered from noisy samples; for a pulse of bandwidth $B$ observed at signal-to-noise ratio $\rho$ the information about the arrival time scales as $B^{2}\rho^{2}$ , which is why timing precision improves linearly with bandwidth. The bound is quoted as the target that a practical algorithm is measured against.

Where the estimators come from. This unit says how to judge an estimator but not how to find one. The two standard constructions are the subject of Maximum likelihood, whose estimators attain the bound asymptotically and exactly whenever any estimator does, and Least-squares fitting, which is the same construction specialised to Gaussian errors. The propagation of the resulting standard errors through subsequent calculations is Propagation of uncertainty.

Method

Maximum likelihood

Adapted

Add to PDF

Defines the likelihood and the maximum-likelihood estimator, proves that maximum likelihood commutes with reparametrisation and that any efficient estimator must be a function of it, states the large-sample normal form of the likelihood with its curvature error bar, and locates a lighthouse from fifteen Cauchy-distributed flashes at $1.163 \pm 0.184$ km where the sample mean gives $-0.64$ km.

Prerequisites: Estimators and their properties, Sampling distributions, Probability spaces, Joint distributions and transformations, Characteristic functions, Markov, Chebyshev, Jensen and Chernoff inequalities, The central limit theorem, Taylor expansion and critical points in several variables

Note(What this unit is for)Original

Estimators and their properties says how to judge a rule for turning data into a number but not where such a rule should come from. Maximum likelihood is the answer that works almost always: write down the probability of having obtained the data you actually obtained, regard it as a function of the unknown parameters, and take the parameters that make it largest.

The method is worth its name for four reasons, all of them below. It is automatic - no cleverness is needed once the model is written down. It is invariant: estimating a rate and estimating a lifetime give consistent answers. It is asymptotically as good as any method can be, in the precise sense of the Cramer-Rao bound. And it produces its own error bar, from the curvature of the log-likelihood at the peak, with no separate calculation.

Definition(Likelihood, log-likelihood, and the ML estimator)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.5, the definition of the likelihood function, p. 1255, and §31.5.1, eqs. (31.68)-(31.69), pp. 1257-1258

Adopt a hypothesis $H$ that the data $x = (x_1,\dots,x_n)$ were drawn from a joint density of a specified functional form, known except for the values of $M$ parameters $a = (a_1,\dots,a_M)$ . The likelihood function is that density read as a function of the parameters with the data held fixed:

L(x; a) \;=\; P(x \mid a, H) .

Riley, Hobson and Bence make the conditioning on $H$ explicit, as here, because the likelihood is only ever defined relative to an assumed model. When the sample is independent the likelihood factorises,

L(x;a) \;=\; \prod_{i=1}^{n} p(x_i \mid a) ,

but the definition does not require independence, nor that the $x_i$ be drawn from a common law.

The log-likelihood is $\ell(a) = \ln L(x;a)$ . A maximum-likelihood estimator (ML estimator) $\hat a$ is a value of the parameters at which $L(x;a)$ attains its largest value over the allowed parameter range. Since $\ln$ is strictly increasing, $L$ and $\ell$ have the same maximisers and the same stationary points, and for a product likelihood $\ell$ is a sum, so one works with $\ell$ .

If the maximum lies at an interior stationary point, $\hat a$ solves the $M$ likelihood equations

\left.\frac{\partial \ell}{\partial a_i}\right|_{a=\hat a} \;=\; 0 , \qquad i = 1,\dots,M .

The left-hand side is the score of Estimators and their properties, so the likelihood equations say: choose the parameters that make the observed data score zero. Riley, Hobson and Bence emphasise that this is a necessary condition only. They list the shapes that occur in practice - a single interior stationary point, several stationary points of which one is the global maximum, and a maximum sitting at an end-point of the allowed range, with or without stationary points elsewhere - and insist that one verify that a solution found is in fact a maximum.

Method(Three properties of maximum likelihood)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.5.2, the invariance property and eq. (31.76), p. 1261, §31.5.3, the efficiency theorem, p. 1262, and §31.5.6, the large-sample form and its consequences, pp. 1269-1270

Let $\hat a$ be an ML estimator of $a$ under the hypothesis $H$ .

1. Invariance under reparametrisation. For any function $\alpha$ of the parameters, the ML estimator of $\alpha(a)$ is

\widehat{\alpha(a)} \;=\; \alpha(\hat a) .

One therefore never has to decide whether to estimate a rate or its reciprocal. The price is that ML estimators are in general biased: if $\hat a$ is unbiased for $a$ and $\alpha$ is strictly convex (or strictly concave) and non-affine, then $\alpha(\hat a)$ is biased upwards (respectively downwards) for $\alpha(a)$ , by Jensen's inequality (Markov, Chebyshev, Jensen and Chernoff inequalities). Riley, Hobson and Bence give the case of a sample from the density $\tau^{-1}e^{-x/\tau}$ on $(0,\infty)$ : the ML estimators are $\hat\tau = \bar x$ and $\hat\lambda = 1/\bar x$ for $\lambda = 1/\tau$ , with

\mathbb{E}[\hat\tau] = \tau , \qquad \mathbb{E}[\hat\lambda] = \frac{n}{n-1}\lambda ,

so the second is biased although the first is not. They note that the bias of an ML estimator always vanishes as $n\to\infty$ .

2. If an efficient estimator exists, it is a function of the ML estimator. Suppose some estimator $\hat a_{\mathrm{eff}}$ attains the Cramer-Rao bound of Estimators and their properties with $P$ replaced by $L$ . Then $\hat a_{\mathrm{eff}}$ is a function of $\hat a$ . Maximum likelihood therefore never misses an efficient estimator that exists.

3. Large-sample form. In the limit $n \to \infty$ the likelihood tends to a multivariate Gaussian in the parameters, centred on the ML estimate:

L(x;a) \;\longrightarrow\; L_{\max}\exp\left[-\tfrac12 Q(a,\hat a)\right] , \qquad Q(a,\hat a) = (a-\hat a)^{\mathsf T} V^{-1} (a - \hat a) ,

\bigl(V^{-1}\bigr)_{ij} \;=\; \left.-\frac{\partial^{2}\ell}{\partial a_i \partial a_j}\right|_{a=\hat a} ,

and in that limit $V^{-1}$ tends to the Fisher matrix $F$ . Riley, Hobson and Bence draw from this the two statements that ML estimators are asymptotically consistent and asymptotically minimum-variance, attributing the Gaussian shape to the central limit theorem of The central limit theorem; they also note that the corresponding statement for the sampling distribution $P(\hat a \mid a)$ holds, with the roles of $a$ and $\hat a$ exchanged, so that in this limit the two ways of setting an error bar agree.

For one parameter, the last display is the curvature standard error

\hat\sigma_{\hat a} \;=\; \left[\left.-\frac{\partial^{2}\ell}{\partial a^{2}}\right|_{a=\hat a}\right]^{-1/2} ,

which Riley, Hobson and Bence justify from the Bayesian reading of the likelihood (Bayesian inference) and recommend as a measure of the error even when $L$ is not Gaussian.

Idea of the proofOriginal

All three properties are statements about the shape of one function of the parameters.

Invariance is the observation that a change of parameter relabels the horizontal axis of the graph of $L$ but does not move its peak: the highest point of a curve is still the highest point after the axis is stretched. Nothing is assumed about the stretch except that it is a well-defined map.

The efficiency statement is the equality case of the Cauchy-Schwarz step in the Cramer-Rao proof. Equality forces the score to be proportional to $\hat a_{\mathrm{eff}}$ minus its mean. But the ML estimator is where the score vanishes, so wherever the score vanishes that bracket vanishes too - which pins $\hat a_{\mathrm{eff}}$ to the ML estimate.

The large-sample form is Taylor's theorem plus the law of large numbers. For an independent sample $\ell$ is a sum of $n$ terms, so its second derivative is a sum of $n$ terms and grows like $n$ , while the cubic term in the expansion about the peak grows like $n$ as well but is multiplied by a displacement of order $n^{-1/2}$ and so is smaller by $n^{-1/2}$ . Expanding about the maximum, the linear term is absent by stationarity and the quadratic term survives: a parabola in $\ell$ , a Gaussian in $L$ , with width set by the curvature.

Proof(Invariance, and the efficiency theorem)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.5.2, p. 1261, and §31.5.3, p. 1262

Invariance. Let $\alpha$ map the parameter space into some set $A$ , and define the induced likelihood of a value $\alpha_0 \in A$ as the largest likelihood available among parameters carrying that value,

L^{*}(\alpha_0) \;=\; \sup\{\, L(x;a) : \alpha(a) = \alpha_0 \,\} .

This is the only definition available when $\alpha$ is not injective, and it reduces to $L(x;\alpha^{-1}(\alpha_0))$ when it is. Now $\sup_{\alpha_0 \in A} L^{*}(\alpha_0) = \sup_a L(x;a) = L(x;\hat a)$ , because every $a$ is counted exactly once on the left, under $\alpha_0 = \alpha(a)$ . The supremum on the left is therefore attained at $\alpha_0 = \alpha(\hat a)$ , which says $\widehat{\alpha(a)} = \alpha(\hat a)$ . No smoothness, monotonicity or invertibility of $\alpha$ was used.

Bias is not preserved. Suppose $\mathbb{E}[\hat a] = a$ and $\alpha$ is strictly convex and non-affine. Jensen's inequality gives $\mathbb{E}[\alpha(\hat a)] > \alpha(\mathbb{E}[\hat a]) = \alpha(a)$ whenever $\hat a$ is not degenerate, so $\alpha(\hat a)$ is biased upwards; applying the same argument to $-\alpha$ gives a downward bias when $\alpha$ is strictly concave and non-affine. This is why invariance and unbiasedness cannot both be demanded of an estimation rule.

Efficiency. Apply the Cramer-Rao inequality of Estimators and their properties with $P(x\mid a)$ taken to be $L(x;a)$ , which is legitimate under $H$ because they are the same function. Suppose $\hat a_{\mathrm{eff}}$ attains the bound. Then Step 5 of that proof holds with equality, and equality in Cauchy-Schwarz requires the two entries to be proportional:

\frac{\partial \ell}{\partial a} \;=\; c(a)\,\bigl[\hat a_{\mathrm{eff}}(x) - \alpha(a)\bigr] \qquad\text{for almost every } x,

with $\alpha(a) = \mathbb{E}[\hat a_{\mathrm{eff}}]$ and $c(a) \ne 0$ independent of $x$ . Evaluate both sides at $a = \hat a$ , where the left-hand side vanishes by the likelihood equation. Since $c(\hat a)\neq 0$ ,

\hat a_{\mathrm{eff}}(x) \;=\; \alpha(\hat a) ,

so $\hat a_{\mathrm{eff}}$ is the function $\alpha$ of the ML estimator - and, by the invariance property just proved, is itself the ML estimator of $\alpha(a)$ . $\square$

The large-sample statement 3 is quoted from Riley, Hobson and Bence, who derive it from the central limit theorem rather than proving it; the conditions under which it holds are discussed in the remark.

Remark(Where the hypotheses are used, and where they fail)Original

The model $H$ . Everything is conditional on the assumed functional form. If $H$ is wrong, the ML estimator converges to the parameter value whose model law is closest to the truth in an information sense, not to anything one wanted; and the curvature error bar, which measures only how sharply the assumed model discriminates, is then an understatement. Checking $H$ is the business of Hypothesis testing, not of this unit.

Interior maximum, and differentiability. The likelihood equations are necessary only where the maximum is interior and $\ell$ is differentiable. The catalogue of likelihood shapes (Riley, Hobson & Bence (2006), §31.5.1, Fig. 31.6, p. 1257) is a list of the ways this can fail, and each occurs in practice. For a sample from the uniform density on $(0,\theta)$ the likelihood is $\theta^{-n}$ for $\theta \ge \max_i x_i$ and zero below, so it is strictly decreasing wherever it is positive: the ML estimator is $\hat\theta = \max_i x_i$ , at the end-point of the allowed range, and the derivative is never zero. The worked example below has three stationary points, two of them spurious.

The regularity behind statement 3. The large-sample Gaussian form needs the conditions already used for the Cramer-Rao bound - support not depending on the parameters, differentiation under the integral sign, finite non-zero information - together with an interior true value and a third derivative of $\ell$ under control, so that the Taylor expansion can be truncated after the quadratic term. The statement of the large-sample form (Riley, Hobson & Bence (2006), §31.5.6, pp. 1269-1270) does not list these; they are what is being used. The uniform example again shows what is lost: there $n(\theta - \hat\theta)/\theta$ converges to an exponential distribution, not a Gaussian, and the error shrinks like $1/n$ rather than $1/\sqrt n$ , which is faster than the regular rate rather than slower.

Asymptotic, not finite-sample. Statement 3 says nothing at any fixed $n$ . The example in this unit has $n = 15$ , and the curvature error bar $0.184$ and the asymptotic bound $0.183$ agree there to one part in a hundred; that agreement is empirical, not guaranteed.

Two readings of the error bar. The curvature formula for $\hat\sigma_{\hat a}$ uses the observed second derivative at the peak, which depends on the data; the Cramer-Rao bound uses its expectation over samples, the Fisher information. The two coincide in the limit and often differ at finite $n$ . The curvature form (Riley, Hobson & Bence (2006), §31.5.5, eq. (31.85), p. 1266) is reached through the Bayesian identification of the normalised likelihood with a posterior under a flat prior, and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §31.5.5, p. 1266) are explicit that a frequentist regards the resulting interval as an approximation to the one in Confidence intervals.

Method(Maximum likelihood in practice)Original

State the model. Write the density of one observation, $p(x\mid a)$ , with every parameter named and its allowed range given.
Form the log-likelihood. For an independent sample, $\ell(a) = \sum_{i=1}^{n}\ln p(x_i\mid a)$ . Drop any additive term free of $a$ ; it shifts $\ell$ without moving the maximum.
Differentiate and solve. Set $\partial\ell/\partial a_i = 0$ for each $i$ . If the equations are algebraic, solve them. If not, solve numerically by Newton's method on the score, as in Root finding, starting from a crude estimate (a sample mean, a median, a method-of-moments value).
Confirm it is the global maximum. Check that the Hessian $\partial^{2}\ell/\partial a_i \partial a_j$ is negative definite at the solution, and scan $\ell$ across the whole allowed range: multiple stationary points are common, and the maximum may sit at an end-point where no derivative vanishes.
Error bars. Invert the negative Hessian at the peak: $V = [-\partial^{2}\ell/\partial a_i\partial a_j]^{-1}$ , whose diagonal square roots are the standard errors and whose off-diagonal entries are the covariances. Equivalently, for one parameter, find where $\ell$ has fallen by $\tfrac12$ from its peak; for a parabola these agree exactly, and the drop-by- $\tfrac12$ rule is the better of the two when the peak is skew.
Reparametrise freely. If the quantity wanted is $\alpha(a)$ , report $\alpha(\hat a)$ . Its error bar follows by Propagation of uncertainty, $\sigma_{\alpha} = \lvert\alpha'(\hat a)\rvert\,\sigma_{\hat a}$ , to first order - or, better, by repeating step 5 in the new parameter.
Report the bias if it matters. ML estimators are biased at finite $n$ by an amount of order $1/n$ . Quote it when it is comparable with the standard error, which happens when $n$ is small or the parametrisation strongly curved.

Example(Locating a lighthouse from fifteen flashes)Original

A beacon lies somewhere on a straight line parallel to a long, straight shore, at a perpendicular distance $\gamma = 0.50$ km known from a chart, but at an unknown position $\mu$ along it. It emits flashes at instants when its beam points in a direction uniformly distributed over all angles. A detector records the position $x$ at which each flash strikes the shore. If $\phi$ is the angle of the beam from the perpendicular, uniform on $(-\pi/2,\pi/2)$ , then $x = \mu + \gamma\tan\phi$ , and the change of variables of Joint distributions and transformations gives the Cauchy density

p(x\mid\mu) \;=\; \frac{1}{\pi}\,\frac{\gamma}{\gamma^{2} + (x-\mu)^{2}} , \qquad x \in \mathbb{R} .

Fifteen flashes are recorded, at the positions (in km, in the order observed)

1.41,\ \ 2.69,\ \ 1.79,\ \ 0.62,\ \ 0.84,\ \ 2.39,\ \ -29.02,\ \ 1.99,\ \ 1.88,\ \ 1.15,\ \ 0.84,\ \ 0.78,\ \ 0.72,\ \ 1.11,\ \ 1.21 .

The sample mean is useless. Its value is $\bar x = -0.64$ km, dragged far from the cluster by the seventh flash, which went out nearly parallel to the shore. This is not bad luck to be averaged away. Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §4.4, paragraph (4.4.7), p. 108) record that the Cauchy law has no moments, and compute (Grimmett & Stirzaker (2020), §7.4, Example (7.4.10), p. 373) the characteristic function of the average of $n$ independent standard Cauchy variables to be $e^{-|t|}$ , the characteristic function of one of them (Characteristic functions); writing $x_i = \mu + \gamma z_i$ with $z_i$ standard carries that statement to the density here, since $\bar x = \mu + \gamma \bar z$ . The sampling distribution of $\bar x$ is therefore the same for $n=15$ as for $n=1$ : taking more data does not improve it at all. The sample median, $1.15$ km, is far better behaved.

The likelihood. Dropping terms free of $\mu$ ,

\ell(\mu) \;=\; -\sum_{i=1}^{15}\ln\bigl[\gamma^{2} + (x_i-\mu)^{2}\bigr] , \qquad \frac{d\ell}{d\mu} \;=\; \sum_{i=1}^{15}\frac{2\,(x_i-\mu)}{\gamma^{2}+(x_i-\mu)^{2}} .

Each term of the score is a bounded, odd function of $x_i - \mu$ , peaking at $\lvert x_i-\mu\rvert = \gamma$ and decaying beyond: a distant point exerts less pull than a near one, which is exactly why the method is not fooled by the outlier.

Solving. The score is not algebraic in $\mu$ ; solve it numerically. Scanning the range covered by the data, it has three roots, at

\mu = -28.897,\qquad \mu = -27.116,\qquad \mu = 1.1628 ,

with log-likelihoods (in the convention above, with the constant dropped) $-94.168$ , $-95.152$ and $+2.173$ . The first two are the spurious pair that the lone far point manufactures - the case Riley, Hobson and Bence illustrate as several stationary points - and the global maximum is

\hat\mu \;=\; 1.1628\ \text{km} .

Error bar from the curvature. Differentiating the score,

-\frac{d^{2}\ell}{d\mu^{2}} \;=\; \sum_{i=1}^{15} \frac{2\bigl[\gamma^{2}-(x_i-\mu)^{2}\bigr]}{\bigl[\gamma^{2}+(x_i-\mu)^{2}\bigr]^{2}} ,

which evaluated at $\hat\mu$ gives $29.586\ \mathrm{km}^{-2}$ . Hence

\hat\sigma_{\hat\mu} \;=\; (29.586)^{-1/2} \;=\; 0.1838\ \text{km}, \qquad\text{so}\qquad \mu = 1.163 \pm 0.184\ \text{km} .

As a cross-check, the Fisher information of one observation for a Cauchy location parameter is

i(\mu) = \int_{-\infty}^{\infty}\left[\frac{2(x-\mu)}{\gamma^{2}+(x-\mu)^{2}}\right]^{2} p(x\mid\mu)\,dx = \frac{1}{2\gamma^{2}} ,

so the Cramer-Rao bound is $2\gamma^{2}/n = 2(0.25)/15 = 0.03333\ \mathrm{km}^{2}$ , a standard error of $0.1826$ km. Observed curvature and expected information agree to $0.7\%$ here. The interval over which $\ell$ falls by $\tfrac12$ from its peak is $[0.986, 1.359]$ , that is $1.163^{+0.196}_{-0.177}$ : very nearly symmetric, so the Gaussian approximation is doing no harm at $n=15$ .

How good is the median? The sample median of a Cauchy sample is asymptotically normal about $\mu$ with standard deviation $\pi\gamma/(2\sqrt n) = 0.2028$ km, so its efficiency in the sense of Estimators and their properties is

e \;=\; \frac{2\gamma^{2}/n}{\pi^{2}\gamma^{2}/(4n)} \;=\; \frac{8}{\pi^{2}} \;=\; 0.811 .

Robust, simple, and only $19\%$ short of optimal - but the ML estimate is what the data actually support, and it costs one short numerical solve.

Invariance in use. Suppose the quantity wanted is the bearing of the beacon from the detector's origin, $\theta = \arctan(\mu/\gamma)$ . By the invariance property, $\hat\theta = \arctan(1.1628/0.50) = 1.1647$ rad $= 66.73^{\circ}$ , and to first order $\sigma_\theta = \gamma\,\sigma_{\hat\mu}/(\gamma^{2}+\hat\mu^{2}) = 0.50\times 0.1838/1.6021 = 0.0574$ rad $= 3.29^{\circ}$ . No new maximisation was needed.

The numbers to check are $\hat\mu = 1.1628$ , $-d^{2}\ell/d\mu^{2} = 29.586$ and $\hat\sigma_{\hat\mu} = 0.1838$ .

Application(Fitting models to data in the physical sciences)Original

Weighted averages. The commonest use of the method is also the simplest. For $n$ independent measurements $x_i$ of one quantity $\mu$ , each Gaussian with its own known standard deviation $\sigma_i$ , the log-likelihood is $-\tfrac12\sum_i (x_i-\mu)^{2}/\sigma_i^{2}$ up to constants, and the likelihood equation gives the inverse-variance weighted mean

\hat\mu \;=\; \frac{\sum_i x_i/\sigma_i^{2}}{\sum_i 1/\sigma_i^{2}} , \qquad \hat\sigma_{\hat\mu} = \Bigl(\textstyle\sum_i \sigma_i^{-2}\Bigr)^{-1/2} ,

which is how every table of world-average physical constants is constructed ((Riley, Hobson & Bence (2006), §31.5.1, eq. (31.72), p. 1259)). When the measurements are correlated with known covariance matrix $N$ , the same calculation gives $\hat\mu = (\mathbb{1}^{\mathsf T} N^{-1} x)/(\mathbb{1}^{\mathsf T} N^{-1}\mathbb{1})$ , reducing to the above when $N$ is diagonal ((Riley, Hobson & Bence (2006), §31.5.1, pp. 1259-1260)).

Least squares. When the model is a curve $f(x;a)$ observed with Gaussian noise, $\ell = -\tfrac12\chi^{2}(a)$ exactly, so maximising the likelihood is minimising $\chi^{2}$ . That specialisation, together with its matrix solution and its goodness-of-fit test, is Least-squares fitting.

Counting experiments. When the number of events is itself random, the extended likelihood (Riley, Hobson & Bence (2006), §31.5.7, eq. (31.88), p. 1270) multiplies the usual product by a Poisson factor $\lambda^{n}e^{-\lambda}/n!$ for the observed count; this is the standard tool for fitting a signal peak over a background in a particle-physics spectrum, where both the shape and the normalisation carry information.

Beyond the Gaussian world. The example of this unit is the pattern for heavy-tailed data generally: in astronomy, in finance and in geophysics, a model whose likelihood penalises far points only logarithmically gives an estimator that is automatically robust, with no explicit outlier rejection. The same likelihood machinery underlies image deconvolution, survival analysis and phylogenetic inference, and, read through Bayes' theorem as in Bayesian inference, it is the bridge to the posterior distributions used throughout modern data analysis.

Method

Least-squares fitting

Adapted

Add to PDF

Sets up the linear fitting problem with a known noise covariance, solves it in closed form for the parameters and their full covariance matrix, proves the Gauss-Markov theorem that this is the best linear unbiased estimator whatever the noise distribution, shows that for Gaussian noise the minimised chi-squared has $n-M$ degrees of freedom, and fits Planck's constant from six photoelectric points as $h = (6.615 \pm 0.093)\times10^{-34}\,\mathrm{J\,s}$ with $\chi^2 = 2.20$ on $4$ degrees of freedom.

Prerequisites: Orthogonal projection and least squares, Estimators and their properties, Sampling distributions, Maximum likelihood, Singular value decomposition, Quadratic forms and the Rayleigh quotient, Unitary and orthogonal matrices, Spectral theorem for Hermitian matrices, The multivariate normal distribution

Note(What this unit is for)Original

A fit is an estimation problem in disguise. Given points that should lie on a curve of known shape but unknown parameters, and given how noisy each point is, the parameters are to be estimated and their uncertainties quoted. If the noise is Gaussian this is maximum likelihood (Maximum likelihood) and nothing new is needed - except that for a model linear in its parameters the maximisation can be done once and for all, in closed form, giving the parameters and their covariance matrix from one matrix inversion.

Three things make this worth a unit of its own. The closed form is exact, not asymptotic. Its optimality survives the loss of normality: the Gauss-Markov theorem below needs only that the errors have zero mean and the stated covariance. And the minimised value of the objective function is itself a test of whether the model was right, with a distribution known in advance.

Definition(The fitting problem and the chi-squared form)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.6, pp. 1271-1272, and §31.6.1, eqs. (31.92)-(31.94), pp. 1272-1273

Let data be a list of pairs $(x_i, y_i)$ , $i = 1,\dots,n$ , in which the $x_i$ are known exactly and each $y_i$ carries a measurement error. Assume

y_i \;=\; f(x_i; a) + \eta_i ,

where $f(\,\cdot\,;a)$ is a function of known form depending on $M$ unknown parameters $a = (a_1,\dots,a_M)$ , and the noise vector $\eta = (\eta_1,\dots,\eta_n)$ has zero mean and known covariance matrix

\mathsf{N}_{ij} \;=\; \operatorname{Cov}[\eta_i,\eta_j] ,

assumed symmetric positive definite. (We write $n$ for the number of data points and $\mathsf{N}$ for the noise covariance; Riley, Hobson and Bence use the letter $N$ for both.) The $x_i$ may be vectors; nothing below changes.

The chi-squared function of the parameters is the quadratic form

\chi^{2}(a) \;=\; \bigl(y - f(a)\bigr)^{\mathsf T}\mathsf{N}^{-1}\bigl(y - f(a)\bigr) \;=\; \sum_{i,j=1}^{n}\bigl[y_i - f(x_i;a)\bigr]\bigl(\mathsf{N}^{-1}\bigr)_{ij}\bigl[y_j - f(x_j;a)\bigr] ,

with $f(a)$ the column vector of model values $f(x_i;a)$ . When the errors are uncorrelated, $\mathsf{N} = \operatorname{diag}(\sigma_1^{2},\dots,\sigma_n^{2})$ and this collapses to the familiar weighted sum

\chi^{2}(a) \;=\; \sum_{i=1}^{n}\left(\frac{y_i - f(x_i;a)}{\sigma_i}\right)^{2} .

The least-squares estimators $\hat a$ are the parameter values minimising $\chi^{2}$ , found in the interior from $\partial\chi^{2}/\partial a_i = 0$ .

The model is linear in the parameters if

f(x;a) \;=\; \sum_{j=1}^{M} a_j\,h_j(x)

for a fixed, linearly independent set of basis functions $h_1,\dots,h_M$ , which may themselves be arbitrarily non-linear in $x$ . Collecting their values at the data points into the $n\times M$ response matrix

R_{ij} \;=\; h_j(x_i) ,

the model vector is $f(a) = Ra$ and

\chi^{2}(a) \;=\; (y - Ra)^{\mathsf T}\mathsf{N}^{-1}(y - Ra) .

Throughout the statements below, $R$ is assumed to have rank $M$ , which requires $n \ge M$ and says that no parameter combination leaves the model values at the data points unchanged.

Method(The linear least-squares solution)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.6.1, eqs. (31.95)-(31.99), pp. 1273-1274

For a model linear in its parameters, with $R$ of rank $M$ and $\mathsf{N}$ symmetric positive definite:

$\chi^{2}$ has a unique minimiser, the solution of the normal equations

R^{\mathsf T}\mathsf{N}^{-1} R\,\hat a \;=\; R^{\mathsf T}\mathsf{N}^{-1} y , \qquad\text{that is}\qquad \hat a \;=\; \bigl(R^{\mathsf T}\mathsf{N}^{-1}R\bigr)^{-1}R^{\mathsf T}\mathsf{N}^{-1}y \;\equiv\; Sy .

$\hat a$ is a linear function of the data and is unbiased, $\mathbb{E}[\hat a] = a$ , for every $a$ .
Its covariance matrix is

V \;=\; \operatorname{Cov}[\hat a_i,\hat a_j] \;=\; S\,\mathsf{N}\,S^{\mathsf T} \;=\; \bigl(R^{\mathsf T}\mathsf{N}^{-1}R\bigr)^{-1} ,

so the matrix inverted to find the parameters already is their covariance matrix. Equivalently $V^{-1} = \tfrac12\,\partial^{2}\chi^{2}/\partial a_i\partial a_j$ , which is independent of the data. 4. $\chi^{2}$ is exactly quadratic:

\chi^{2}(a) \;=\; \chi^{2}(\hat a) + (a - \hat a)^{\mathsf T}V^{-1}(a-\hat a) ,

with no higher terms, so the surfaces of constant $\chi^{2}$ are the concentric ellipsoids of Quadratic forms and the Rayleigh quotient centred on $\hat a$ .

If in addition the noise is Gaussian, $\eta \sim N(0,\mathsf{N})$ in the sense of The multivariate normal distribution, then $L \propto \exp(-\tfrac12\chi^{2})$ , so the least-squares and maximum-likelihood estimators coincide, $V^{-1}$ is the Fisher matrix of Estimators and their properties, and $\hat a$ attains the Cramer-Rao bound exactly, at every $n$ .

Method(Gauss-Markov)Original

Keep the linear model $y = Ra + \eta$ with $R$ of rank $M$ , and assume only that

\mathbb{E}[\eta] = 0 , \qquad \operatorname{Cov}[\eta] = \mathsf{N} \ \text{ symmetric positive definite and known.}

No distributional assumption is made. Let $\tilde a = Cy$ be any estimator that is linear in the data and unbiased for every value of $a$ . Then

\operatorname{Cov}[\tilde a] - \operatorname{Cov}[\hat a] \ \text{ is positive semi-definite,}

with $\hat a = Sy$ the least-squares estimator. In particular $\operatorname{var}[\tilde a_i] \ge \operatorname{var}[\hat a_i] = V_{ii}$ for each parameter, and the same holds for every linear combination $u^{\mathsf T}a$ of the parameters. Equality for all $u$ forces $C = S$ .

The least-squares estimator is thus the best linear unbiased estimator: among rules of the form "a fixed matrix times the data", it is the one of smallest variance, whatever the law of the noise.

Idea of the proofOriginal

Both statements are one geometric picture.

Whiten the problem. Because $\mathsf{N}$ is symmetric positive definite it has a symmetric positive definite square root, and dividing through by it turns the weighted form $\chi^{2}$ into an ordinary squared length: with $z = \mathsf{N}^{-1/2}y$ and $B = \mathsf{N}^{-1/2}R$ ,

\chi^{2}(a) \;=\; \lVert z - Ba\rVert^{2} .

Minimising a squared distance from a point to the column space of $B$ is exactly the problem solved in Orthogonal projection and least squares: the answer is the orthogonal projection of $z$ onto that space, characterised by the residual being perpendicular to it, which is the normal equations. The covariance then follows because a linear map sends a covariance matrix $\mathsf{N}$ to $S\mathsf{N}S^{\mathsf T}$ , and the whitening makes that collapse.

Gauss-Markov is the Pythagorean theorem in the same picture. Any other linear unbiased rule differs from the projection by a matrix that annihilates the model space, so its error splits into the projection error plus an orthogonal piece; the two are uncorrelated, their covariances add, and adding a covariance matrix can only make the total larger.

Finally, the leftover $\chi^{2}(\hat a)$ is the squared length of the part of $z$ perpendicular to an $M$ -dimensional subspace. For Gaussian noise $z$ is a standard normal vector, whose law is rotation-invariant, so that squared length is a sum of $n-M$ independent squared standard normals: a chi-squared variable with $n-M$ degrees of freedom. The $M$ lost dimensions are the parameters that have been fitted.

Proof(The solution, its covariance, and the Gauss-Markov bound)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.6.1, the gradient (31.95), the solution (31.98) and the covariance (31.99), pp. 1273-1274

Write $W = \mathsf{N}^{-1}$ , symmetric positive definite.

Step 1: the gradient. Expanding the quadratic form,

\chi^{2}(a) \;=\; y^{\mathsf T}Wy - 2a^{\mathsf T}R^{\mathsf T}Wy + a^{\mathsf T}R^{\mathsf T}WR\,a ,

where the two cross terms were combined using $W^{\mathsf T} = W$ . Differentiating in $a_p$ , only the last two terms contribute, and

\frac{\partial\chi^{2}}{\partial a_p} \;=\; -2\bigl(R^{\mathsf T}Wy\bigr)_p + 2\bigl(R^{\mathsf T}WR\,a\bigr)_p , \qquad\text{that is}\qquad \nabla\chi^{2} \;=\; -2R^{\mathsf T}W(y - Ra) .

Step 2: the normal equations and uniqueness. Setting the gradient to zero gives $R^{\mathsf T}WR\,\hat a = R^{\mathsf T}Wy$ . The matrix $A = R^{\mathsf T}WR$ is positive definite: for $v \ne 0$ , $v^{\mathsf T}Av = (Rv)^{\mathsf T}W(Rv) > 0$ , because $Rv \ne 0$ (rank $M$ , so $R$ is injective) and $W$ is positive definite. Both hypotheses are used here, and each separately: rank $M$ to get $Rv\neq0$ , positive definiteness of $\mathsf{N}$ to get a strictly positive value. Hence $A$ is invertible, the stationary point is unique, and since the Hessian $\partial^{2}\chi^{2}/\partial a_i\partial a_j = 2A_{ij}$ is positive definite and constant, $\chi^{2}$ is strictly convex and the stationary point is its global minimum. This gives claim 1 with $S = A^{-1}R^{\mathsf T}W$ , and claim 4 is Taylor's theorem for a quadratic, exact because all third derivatives vanish.

Step 3: unbiasedness. $\hat a = Sy = S(Ra + \eta) = SRa + S\eta$ , and

SR \;=\; A^{-1}R^{\mathsf T}WR \;=\; A^{-1}A \;=\; I_M ,

so $\hat a = a + S\eta$ and $\mathbb{E}[\hat a] = a$ , using $\mathbb{E}[\eta] = 0$ . Claim 2.

Step 4: the covariance. From $\hat a - a = S\eta$ ,

V \;=\; \mathbb{E}\bigl[S\eta\,\eta^{\mathsf T}S^{\mathsf T}\bigr] \;=\; S\,\mathsf{N}\,S^{\mathsf T} \;=\; A^{-1}R^{\mathsf T}W\,\mathsf{N}\,W R A^{-1} ,

using $S^{\mathsf T} = W R A^{-1}$ ( $A$ and $W$ are symmetric). Since $W\mathsf{N}W = W$ , the middle collapses to $R^{\mathsf T}WR = A$ , leaving $V = A^{-1}AA^{-1} = A^{-1}$ . Claim 3. Note that only $\mathbb{E}[\eta] = 0$ and $\operatorname{Cov}[\eta] = \mathsf{N}$ were used in Steps 3 and 4: no normality.

Step 5: Gauss-Markov. Let $\tilde a = Cy$ be linear and unbiased. Unbiasedness means $\mathbb{E}[Cy] = CRa = a$ for every $a$ , which forces

CR \;=\; I_M .

Put $D = C - S$ . Then $DR = CR - SR = I_M - I_M = 0$ by Step 3. Now

\operatorname{Cov}[\tilde a] \;=\; C\mathsf{N}C^{\mathsf T} \;=\; (S+D)\mathsf{N}(S+D)^{\mathsf T} \;=\; S\mathsf{N}S^{\mathsf T} + S\mathsf{N}D^{\mathsf T} + D\mathsf{N}S^{\mathsf T} + D\mathsf{N}D^{\mathsf T} .

The cross terms vanish: $S\mathsf{N}D^{\mathsf T} = A^{-1}R^{\mathsf T}W\mathsf{N}D^{\mathsf T} = A^{-1}R^{\mathsf T}D^{\mathsf T} = A^{-1}(DR)^{\mathsf T} = 0$ , and the other is its transpose. Hence

\operatorname{Cov}[\tilde a] \;=\; V + D\,\mathsf{N}\,D^{\mathsf T} ,

and $D\mathsf{N}D^{\mathsf T}$ is positive semi-definite, since $u^{\mathsf T}D\mathsf{N}D^{\mathsf T}u = (D^{\mathsf T}u)^{\mathsf T}\mathsf{N}(D^{\mathsf T}u)\ge 0$ . That is the claim; taking $u = e_i$ gives the statement for a single parameter and general $u$ the statement for a linear combination. Equality for all $u$ forces $D^{\mathsf T}u = 0$ for all $u$ , since $\mathsf{N}$ is positive definite, hence $D = 0$ and $C = S$ . $\square$

Method(The minimised chi-squared is a goodness-of-fit statistic)AdaptedAdapted from Arfken, Weber & Harris (2013), §23.7, eqs. (23.99)-(23.101) and the degrees-of-freedom rule, pp. 1170-1172

Assume in addition that the noise is Gaussian, $\eta \sim N(0,\mathsf{N})$ , and keep $R$ of rank $M$ . Then the minimised value

\chi^{2}_{\min} \;=\; \chi^{2}(\hat a) \;=\; \bigl(y - R\hat a\bigr)^{\mathsf T}\mathsf{N}^{-1}\bigl(y - R\hat a\bigr)

is distributed as $\chi^{2}(n-M)$ , and is independent of $\hat a$ . Hence

\mathbb{E}\bigl[\chi^{2}_{\min}\bigr] = n-M , \qquad \operatorname{var}\bigl[\chi^{2}_{\min}\bigr] = 2(n-M) .

Arfken, Weber and Harris state the counting rule in this form: fitting $M$ parameters to $n$ points leaves $n - M$ effective degrees of freedom. A fit is judged by comparing $\chi^{2}_{\min}$ with $n-M$ and quoting the upper tail probability $\Pr\bigl(\chi^{2}(n-M) > \chi^{2}_{\min}\bigr)$ ; a value far above $n-M$ indicates a wrong model or understated errors, and one far below indicates overstated errors or too many parameters.

If $\mathsf{N}$ is known only up to an overall scale, $\mathsf{N} = \sigma^{2}\mathsf{W}$ with $\mathsf{W}$ known and $\sigma$ unknown, then $\hat a$ is unchanged (the scale cancels in the normal equations),

\hat\sigma^{2} \;=\; \frac{\chi^{2}_{\min}\bigl[\mathsf{W}\bigr]}{\,n-M\,}

is an unbiased estimator of $\sigma^{2}$ , where $\chi^{2}_{\min}[\mathsf{W}]$ is the minimised form computed with $\mathsf{W}$ in place of $\mathsf{N}$ , and the covariance estimate is $\hat\sigma^{2}(R^{\mathsf T}\mathsf{W}^{-1}R)^{-1}$ . No goodness-of-fit test is then available, because the fit has been used to set the error scale.

Proof(Degrees of freedom)Original

Whiten: let $\mathsf{N}^{1/2}$ be the symmetric positive definite square root of $\mathsf{N}$ (it exists by the spectral theorem, Spectral theorem for Hermitian matrices), and set

z \;=\; \mathsf{N}^{-1/2}y , \qquad B \;=\; \mathsf{N}^{-1/2}R , \qquad \zeta \;=\; \mathsf{N}^{-1/2}\eta .

Then $\chi^{2}(a) = \lVert z - Ba\rVert^{2}$ , $z = Ba + \zeta$ , and $\zeta \sim N(0, I_n)$ : the whitened noise is a standard normal vector, by the affine transformation rule for multivariate normals (The multivariate normal distribution). The matrix $B$ still has rank $M$ , since $\mathsf{N}^{-1/2}$ is invertible.

Let $P$ be the orthogonal projection of $\mathbb{R}^{n}$ onto the $M$ -dimensional column space of $B$ (Orthogonal projection and least squares). The minimiser of $\lVert z - Ba\rVert^{2}$ makes $Ba$ equal to $Pz$ , so

\chi^{2}_{\min} \;=\; \lVert (I-P)z\rVert^{2} \;=\; \lVert (I-P)\zeta\rVert^{2} ,

the second equality because $(I-P)Ba = 0$ for every $a$ , $Ba$ lying in the range of $P$ .

Now choose an orthonormal basis of $\mathbb{R}^{n}$ whose first $M$ vectors span the range of $B$ and whose last $n-M$ span its orthogonal complement, and let $u = Q^{\mathsf T}\zeta$ be the coordinates of $\zeta$ in that basis, $Q$ orthogonal. By the rotation invariance of the standard normal law used in Sampling distributions - orthogonal maps preserve $\lVert\cdot\rVert$ and have unit Jacobian, so the density $(2\pi)^{-n/2}e^{-\lVert u\rVert^{2}/2}$ is unchanged - the coordinates $u_1,\dots,u_n$ are again independent standard normals. In these coordinates

\chi^{2}_{\min} \;=\; \sum_{k=M+1}^{n}u_k^{2} , \qquad \hat a - a \;=\; (B^{\mathsf T}B)^{-1}B^{\mathsf T}\zeta \ \text{ depends only on } u_1,\dots,u_M ,

the second because $B^{\mathsf T}\zeta = B^{\mathsf T}P\zeta$ ( $B^{\mathsf T}(I-P) = 0$ , the columns of $B$ being in the range of $P$ ). So $\chi^{2}_{\min}$ is a sum of $n-M$ independent squared standard normals, that is a $\chi^{2}(n-M)$ variable, and it is a function of coordinates independent of those determining $\hat a$ . The mean $n-M$ and variance $2(n-M)$ are those of the chi-squared law recorded in Sampling distributions.

For the scale-only-unknown case, apply the above with $\mathsf{W}$ in place of $\mathsf{N}$ : then $\chi^{2}_{\min}[\mathsf{W}]/\sigma^{2} \sim \chi^{2}(n-M)$ , whose mean is $n-M$ , so $\hat\sigma^{2}$ is unbiased. $\square$

Remark(Where the hypotheses are used)Original

$R$ of rank $M$ . Used in Step 2 to invert $A = R^{\mathsf T}\mathsf{N}^{-1}R$ . If it fails, some combination of parameters is invisible to the data - fitting $a_1 + a_2 x$ and $a_3 x$ simultaneously, say - and $\chi^{2}$ has a flat valley rather than a minimum point. Near-degeneracy is the practical danger: $A$ is then ill-conditioned, and the fitted parameters are enormous and anticorrelated while the fitted curve is perfectly well determined. The cure is not a better inversion routine but a better basis. The singular value decomposition of Singular value decomposition diagnoses it (the ratio of largest to smallest singular value of $\mathsf{N}^{-1/2}R$ is the condition number of Matrix norms and condition numbers) and supplies the minimum-norm solution when the rank really is deficient. Riley, Hobson and Bence make the same point for polynomial fitting and recommend replacing $1,x,x^{2},\dots$ by polynomials orthogonal over the data points, which makes $A$ diagonal.

$\mathsf{N}$ positive definite and known. Positive definiteness is used in Step 2 and in the equality case of Step 5; it fails if two data points are perfectly correlated, which is a statement that one of them is redundant. That $\mathsf{N}$ is known is the most often violated assumption in practice: error bars are themselves estimated. When only the scale is unknown, the last paragraph of the statement applies; when the shape is unknown too, nothing here is valid and the problem becomes one of variance estimation.

$\mathbb{E}[\eta] = 0$ . Used in Steps 3 and 5. A systematic offset common to all points is absorbed into any parameter that can imitate it (an intercept) and biases the others; no amount of data removes it, and the quoted $V$ does not describe it.

The $x_i$ exact. Assumed throughout. When the abscissae carry the error instead, the roles of the coordinates are exchanged and one fits $x$ against $y$ , which gives a different line. When both carry comparable error, Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §23.7, eq. (23.98) and Fig. 23.11, pp. 1169-1170) minimise the sum of squared perpendicular distances instead, parametrising the line by the angle $\alpha$ it makes with the abscissa and obtaining $\tan 2\alpha = 2\sum_j x_j y_j / \sum_j (x_j^{2}-y_j^{2})$ . That estimator is appropriate when the two coordinates are measured with equal precision; it is not the same as either one-sided fit.

Normality. Used only for the identification with maximum likelihood, for the Cramer-Rao statement, and for the distribution of $\chi^{2}_{\min}$ . It is not used for claims 1-4 nor for Gauss-Markov. This is the practical importance of the Gauss-Markov theorem: a fit remains the best linear unbiased estimate even when the errors are not Gaussian, although its $\chi^{2}$ then has no tabulated distribution and the goodness-of-fit reading is lost.

Linear in the parameters, not in $x$ . The basis functions may be anything - powers, sinusoids, exponentials with fixed decay rates. Fitting an exponential with an unknown rate is not a linear problem, and Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §31.6.2, pp. 1276-1277) treat that case by Newton iteration on $\nabla\chi^{2}$ , starting from a guess and solving the linear system $\sum_j (\partial^{2}\chi^{2}/\partial a_i\partial a_j)\,\delta a_j = -\partial\chi^{2}/\partial a_i$ repeatedly. Everything above then holds only approximately and locally, and they warn that $\chi^{2}$ may have many local minima, so that the quality of the initial guess decides which one is found.

A difference in how the optimality is stated. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §31.6.1, p. 1273) write that for a model linear in the parameters the least-squares estimators are unbiased and of minimum variance, whatever the probability density the measurement errors come from. The theorem proved above is the same assertion with one restriction added: minimum-variance among estimators that are linear functions of the data. The restriction cannot be dropped. For heavy-tailed noise a non-linear estimator can beat least squares badly - the lighthouse example of Maximum likelihood is exactly such a case, where the least-squares estimate of a location parameter is the sample mean and is useless, while the maximum-likelihood estimate is not linear in the data and is excellent. This unit follows the restricted statement, which is the one that is true without further assumptions; under normality the two coincide, because then the best linear unbiased estimator is also the best unbiased estimator by the Cramer-Rao argument of claim 4.

Method(Fitting a linear model)Original

Write the model as a basis expansion. Choose $h_1,\dots,h_M$ with $f(x;a) = \sum_j a_j h_j(x)$ . For a straight line, $h_1 = 1$ and $h_2 = x$ .
Build the response matrix $R_{ij} = h_j(x_i)$ , of size $n\times M$ , and the noise covariance $\mathsf{N}$ . For independent errors $\mathsf{N}^{-1} = \operatorname{diag}(w_1,\dots,w_n)$ with weights $w_i = 1/\sigma_i^{2}$ .
Form and solve the normal equations. Compute $A = R^{\mathsf T}\mathsf{N}^{-1}R$ and $b = R^{\mathsf T}\mathsf{N}^{-1}y$ , then solve $A\hat a = b$ by Cholesky factorisation (Solving linear systems numerically), not by forming $A^{-1}$ , unless $M$ is small. For $M = 2$ the explicit $2\times2$ inverse is quicker and is used in the example.
Read off the errors. $V = A^{-1}$ ; the standard error on $\hat a_i$ is $\sqrt{V_{ii}}$ and the correlation coefficient of $\hat a_i$ with $\hat a_j$ is $V_{ij}/\sqrt{V_{ii}V_{jj}}$ . Quote the whole matrix, not just the diagonal: the off-diagonal entries are needed by anything downstream, as in Propagation of uncertainty.
Decorrelate if the correlation is large. Replacing $x$ by $x - \bar x_w$ , where $\bar x_w = \sum_i w_i x_i/\sum_i w_i$ is the weighted mean abscissa, makes $A$ diagonal for a straight line and the two parameters uncorrelated. The fitted line is unchanged; only the parametrisation is.
Test the fit. Compute $\chi^{2}_{\min}$ , compare with $n-M$ , and quote the tail probability. Only then believe the error bars.
Confidence regions. For a joint statement about several parameters use the ellipsoid $(a-\hat a)^{\mathsf T}V^{-1}(a-\hat a) \le c$ , with $c$ the upper point of a chi-squared distribution with $M$ degrees of freedom at the wanted level: $c = 1$ gives the one-parameter one-sigma interval for $M=1$ , while for $M=2$ a joint $68.3\%$ region needs $c = 2.30$ . Confidence intervals develops this.

Example(Planck's constant from the photoelectric effect)Original

A caesium-coated photocathode is illuminated with six mercury lines in turn and the retarding voltage $V$ that just stops the photocurrent is measured for each. Einstein's relation predicts

V \;=\; \frac{h}{e}\,\nu - \frac{\varphi}{e} ,

a straight line in the frequency $\nu$ whose slope gives $h/e$ and whose intercept gives the work function. Writing $x = \nu/10^{14}\,\mathrm{Hz}$ , the data are

$x$	$V/\mathrm{V}$	$\sigma/\mathrm{V}$	$w = 1/\sigma^{2}$
$5.187$	$0.258$	$0.015$	$4444.444$
$5.490$	$0.361$	$0.020$	$2500.000$
$6.879$	$0.951$	$0.012$	$6944.444$
$7.404$	$1.140$	$0.030$	$1111.111$
$8.219$	$1.515$	$0.018$	$3086.420$
$8.976$	$1.805$	$0.025$	$1600.000$

Here $M = 2$ , with $h_1 = 1$ and $h_2 = x$ , so $R$ has a column of ones and a column of the $x_i$ , and $\mathsf{N}^{-1} = \operatorname{diag}(w_i)$ .

The normal matrix. Its entries are weighted sums:

A = R^{\mathsf T}\mathsf{N}^{-1}R = \begin{pmatrix} S_0 & S_1\\ S_1 & S_2\end{pmatrix}, \qquad b = R^{\mathsf T}\mathsf{N}^{-1}y = \begin{pmatrix} T_0 \\ T_1 \end{pmatrix},

with $S_k = \sum_i w_i x_i^{k}$ and $T_k = \sum_i w_i x_i^{k} y_i$ . From the table,

S_0 = 19686.420,\quad S_1 = 132504.717,\quad S_2 = 921857.121,

T_0 = 17483.926,\quad T_1 = 130065.071,

and the determinant is

\Delta = S_0S_2 - S_1^{2} = 18\,148\,066\,234 - 17\,557\,500\,103 = 5.90566\times10^{8} .

The parameters. Inverting the $2\times2$ matrix,

\hat m = \frac{S_0T_1 - S_1T_0}{\Delta} = \frac{2.43813\times10^{8}}{5.90566\times10^{8}} = 0.412846\ \mathrm{V\,per\,10^{14}\,Hz} ,

\hat c = \frac{S_2T_0 - S_1T_1}{\Delta} = \frac{-1.11655\times10^{9}}{5.90566\times10^{8}} = -1.890650\ \mathrm{V} .

The covariance matrix is $V = A^{-1} = \Delta^{-1}\begin{pmatrix} S_2 & -S_1 \\ -S_1 & S_0\end{pmatrix}$ , so

\sigma_{\hat c} = \sqrt{S_2/\Delta} = 0.039509\ \mathrm{V}, \qquad \sigma_{\hat m} = \sqrt{S_0/\Delta} = 0.005774 ,

\operatorname{Cov}[\hat c,\hat m] = -S_1/\Delta = -2.2437\times10^{-4} , \qquad \rho = \frac{-S_1}{\sqrt{S_0S_2}} = -0.9836 .

The two parameters are almost perfectly anticorrelated, because the data sit far from $x=0$ and a small increase in slope is compensated by a large decrease in intercept. Re-fitting in the variable $x - \bar x_w$ with $\bar x_w = S_1/S_0 = 6.7308$ makes the correlation exactly zero and the height of the line at that point is determined to $\sigma = 1/\sqrt{S_0} = 0.00713\,\mathrm{V}$ - five and a half times better than the intercept at the origin.

Goodness of fit. The fitted values $\hat c + \hat m x_i$ and the pulls $(y_i - \hat c - \hat m x_i)/\sigma_i$ are

0.2508,\ 0.3759,\ 0.9493,\ 1.1661,\ 1.5025,\ 1.8151 ; \qquad 0.481,\ -0.744,\ 0.140,\ -0.869,\ 0.693,\ -0.402 .

Squaring and adding the pulls,

\chi^{2}_{\min} = 0.2315+0.5532+0.0196+0.7547+0.4798+0.1618 = 2.2006

on $n - M = 4$ degrees of freedom. The expected value is $4$ with standard deviation $\sqrt{8}=2.83$ , and the upper tail is $\Pr(\chi^{2}(4) > 2.2006) = 0.699$ : an unremarkable fit, with no sign that the straight line is wrong or that the error bars are mis-stated.

The physics. The slope is $h/e$ in volt-seconds once the $10^{14}$ is restored:

\frac{h}{e} = (4.12846 \pm 0.05774)\times10^{-15}\ \mathrm{V\,s} ,

h = (6.6145 \pm 0.0925)\times10^{-34}\ \mathrm{J\,s} ,

using $e = 1.602177\times10^{-19}\,\mathrm{C}$ . The accepted value $6.62607\times10^{-34}\,\mathrm{J\,s}$ lies $0.12$ standard errors away. The intercept gives the work function $\varphi = -\hat c = 1.891 \pm 0.040\ \mathrm{eV}$ , and the threshold frequency is

\nu_0 = -\hat c/\hat m = 4.5796\times10^{14}\ \mathrm{Hz} \qquad (\lambda_0 = 654.6\ \mathrm{nm}) ,

whose standard error must be computed from the full covariance matrix, since $\nu_0$ depends on both parameters. With $\partial\nu_0/\partial \hat c = -1/\hat m$ and $\partial\nu_0/\partial\hat m = \hat c/\hat m^{2}$ , first-order propagation (Propagation of uncertainty) gives $\sigma_{\nu_0} = 0.0347\times10^{14}\,\mathrm{Hz}$ . Had the correlation been ignored, the two terms would have added instead of largely cancelling and the error would have come out more than three times too large.

The numbers to check are $\hat m = 0.412846$ , $\hat c = -1.890650$ , $\chi^{2}_{\min} = 2.2006$ on $4$ degrees of freedom, and $h = 6.6145\times10^{-34}\,\mathrm{J\,s}$ .

Application(Calibration, trends, and everything downstream)Original

Calibration curves. Every instrument that converts a reading into a physical quantity carries a fitted calibration, and the covariance matrix of its parameters is what propagates into the uncertainty on every later measurement made with it. Quoting only the diagonal of $V$ - "slope $\pm$ , intercept $\pm$ " - loses the anticorrelation and, as the threshold frequency above shows, can inflate a downstream error by a large factor.

Trends in time series. Fitting a straight line to a record of global mean temperature, sea level or station pressure is the same calculation with $x$ the date. The subtlety is that $\mathsf{N}$ is then not diagonal: geophysical noise is autocorrelated, as in Autoregressive models and red noise, and treating it as white understates $\sigma_{\hat m}$ badly, sometimes by a factor of several. The machinery above handles it correctly provided the autocovariance is supplied, since nothing in the derivation required $\mathsf{N}$ to be diagonal.

Orbit and trajectory determination. Fitting six orbital elements to many tracking observations is linear least squares after one linearisation about a reference orbit, iterated as in the non-linear scheme. The same structure, applied recursively as data arrive rather than all at once, is the The Kalman filter, in which $V$ is updated at each step instead of being computed from scratch.

Model selection and testing. Because $\chi^{2}_{\min}$ has a known distribution, the difference of the minimised values of two nested models is itself a test statistic: adding a parameter must reduce $\chi^{2}$ by at least about one unit to be worth anything, and the formal version of that statement is the likelihood-ratio test of Hypothesis testing. The same ellipsoids $(a-\hat a)^{\mathsf T}V^{-1}(a-\hat a) = c$ that bound confidence regions are the contours used to display parameter degeneracies in fields from cosmology to pharmacokinetics.

Method

Propagation of uncertainty

Adapted

Add to PDF

Gives the first-order rule for the variance of a function of uncertain inputs, in the matrix form $V_Y = S V_X S^{\mathsf T}$ that keeps the correlations, says exactly when the rule is exact and when its neglected second-order terms matter, and works two examples ending in $g = 9.811 \pm 0.044\,\mathrm{m\,s^{-2}}$ and a predicted-value error of $0.134$ that the naive uncorrelated formula overstates by a factor of $3.2$ .

Prerequisites: Expectation, variance and covariance, Taylor expansion and critical points in several variables, Joint distributions and transformations

Note(What the rule is and is not)Original

A measurement is reported as $x \pm \sigma_x$ , and what is wanted is almost never $x$ itself but some function of it: a period becomes a frequency, two lengths become an aspect ratio, counts and a live-time become a rate. Propagation of uncertainty answers the question "what is the error bar on the answer?" by one device only: replace the function by its tangent plane at the measured point, and use the exact formula for the variance of a linear combination.

Two consequences are worth stating before any algebra. First, because the device is linearisation, the rule is exact when the function is affine and approximate otherwise, with an error controlled by the curvature of the function over the width of the input distribution. Second, because the exact formula for the variance of a linear combination contains covariances, so does the propagation rule. Dropping them is the commonest error in practice, and this unit's second example shows it inflating an error bar by a factor of $3.2$ .

Definition(Inputs, outputs and the covariance matrix)AdaptedAdapted from Riley, Hobson & Bence (2006), §30.12.3, the definition (30.139) of the covariance matrix and the set-up preceding (30.140), pp. 1203-1205

Let $X = (X_1, \dots, X_n)$ be random variables with finite second moments, means $\mu_i = \mathbb{E}[X_i]$ and covariance matrix $V_X$ , the symmetric $n \times n$ matrix with entries

(V_X)_{ij} \;=\; \mathrm{Cov}[X_i, X_j] \;=\; \mathbb{E}\bigl[(X_i - \mu_i)(X_j - \mu_j)\bigr] .

Its diagonal entries are the variances $\sigma_i^2 = V[X_i]$ , and the standard error (the quoted error bar) on $X_i$ is $\sigma_i$ . The correlation matrix has entries $\rho_{ij} = (V_X)_{ij}/(\sigma_i\sigma_j)$ , with unit diagonal. Covariance, correlation and the variance of a linear combination are the subject of Expectation, variance and covariance; $V_X$ is symmetric and positive semi-definite because $t^{\mathsf T}V_X t = V\bigl[\sum_i t_iX_i\bigr] \ge 0$ for every real vector $t$ .

Let $Y = (Y_1, \dots, Y_m)$ with $Y_k = f_k(X_1, \dots, X_n)$ for functions $f_k$ that are differentiable at $\mu = (\mu_1, \dots, \mu_n)$ . The sensitivity matrix (Riley, Hobson and Bence write $S$ ) is the $m \times n$ Jacobian evaluated at the means,

S_{ki} \;=\; \frac{\partial f_k}{\partial x_i}\bigg|_{x = \mu} .

Throughout, "the error on $Y_k$ " means $\sqrt{V[Y_k]}$ , and "errors are independent" means $V_X$ is diagonal. In an experimental report the mean $\mu_i$ is not known and the measured value $x_i$ is used in its place; the remark below says what that substitution costs.

Method(Linearised propagation of uncertainty)AdaptedAdapted from Riley, Hobson & Bence (2006), §30.12.3, eqs. (30.136)-(30.141), pp. 1203-1205; Arfken, Weber & Harris (2013), §23.7, eqs. (23.88) and (23.91), pp. 1166-1167

Assume:

(P1) each $X_i$ has finite variance;
(P2) each $f_k$ is differentiable at $\mu$ , with sensitivity matrix $S$ as above;
(P3) the distribution of $X$ is concentrated near $\mu$ in the sense made precise in the remark: the linearisation is a good approximation to $f_k$ over the region where $X$ has its probability.

Then, to first order,

\mathbb{E}[Y_k] \;\approx\; f_k(\mu) , \qquad \mathrm{Cov}[Y_k, Y_l] \;\approx\; \sum_{i=1}^{n}\sum_{j=1}^{n} S_{ki}\,S_{lj}\,(V_X)_{ij} ,

which in matrix form is

\boxed{\;V_Y \;=\; S\,V_X\,S^{\mathsf T} .\;}

For a single output $Y = f(X_1, \dots, X_n)$ this reads

V[Y] \;\approx\; \sum_{i=1}^{n}\left(\frac{\partial f}{\partial x_i}\right)^{\!2} V[X_i] \;+\; \sum_{i \ne j} \frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\,\mathrm{Cov}[X_i, X_j] ,

all derivatives at $x = \mu$ ; for $n = 2$ this is the formula Riley, Hobson and Bence give as (30.138), and when the $X_i$ are uncorrelated the double sum vanishes and only the squared-derivative terms remain, which is Arfken, Weber and Harris's (23.91).

Exactness. If every $f_k$ is affine, $f_k(x) = c_k + \sum_i S_{ki}x_i$ , then both relations hold with equality and no hypothesis beyond (P1) is needed. Riley, Hobson and Bence state this for (30.138) and for (30.140).

Idea of the proofOriginal

A smooth function looks like its tangent plane on a small enough neighbourhood, and the variance of an affine function of random variables is a formula we already have exactly. So: expand $f_k$ to first order about the means, discard the remainder, and apply the exact formula. The constant term $f_k(\mu)$ carries no variance and drops out; the linear term contributes $\sum_i S_{ki}(X_i - \mu_i)$ , whose covariance with the corresponding sum for $Y_l$ is a double sum of covariances weighted by sensitivities. That double sum is the matrix product $SV_XS^{\mathsf T}$ . Everything else in this unit is about the discarded remainder: it vanishes identically when $f_k$ is affine, and otherwise its size is set by the second derivatives of $f_k$ times the variances of the inputs.

ProofAdaptedAdapted from Riley, Hobson & Bence (2006), §30.12.3, the derivation of (30.138) from (30.136) and (30.137), and of (30.140) by the same route, pp. 1203-1205

Step 1 (the exact affine case). Let $Z = c + \sum_i a_iX_i$ and $W = d + \sum_j b_jX_j$ with real constants. Bilinearity of the covariance (Expectation, variance and covariance) gives, with no approximation and using only (P1),

\mathbb{E}[Z] = c + \sum_i a_i\mu_i , \qquad \mathrm{Cov}[Z, W] = \sum_{i}\sum_{j} a_i b_j\,\mathrm{Cov}[X_i, X_j] . \tag{1}

Taking $W = Z$ and $n = 2$ recovers $V[aX+bY+c] = a^2V[X] + b^2V[Y] + 2ab\,\mathrm{Cov}[X,Y]$ , which is the identity Riley, Hobson and Bence use as their starting point.

Step 2 (linearise). By (P2), Taylor's theorem in several variables (Taylor expansion and critical points in several variables) gives

f_k(X) \;=\; f_k(\mu) \;+\; \sum_{i=1}^{n} S_{ki}\,(X_i - \mu_i) \;+\; R_k(X) , \tag{2}

where $R_k(X)/\lVert X - \mu\rVert \to 0$ as $X \to \mu$ . Write $L_k$ for the affine part, so $Y_k = L_k + R_k$ .

Step 3 (discard the remainder). Hypothesis (P3) is the assumption that $R_k$ may be neglected: we replace $Y_k$ by $L_k$ . This is the only approximation in the argument, and it is where (P3) is used. If every $f_k$ is affine then $R_k \equiv 0$ identically and the step is exact.

Step 4 (apply the exact formula). Applying (1) to $L_k$ and $L_l$ , whose coefficients are $a_i = S_{ki}$ and $b_j = S_{lj}$ and whose constants are $f_k(\mu) - \sum_i S_{ki}\mu_i$ ,

\mathbb{E}[L_k] = f_k(\mu) , \qquad \mathrm{Cov}[L_k, L_l] = \sum_{i}\sum_{j} S_{ki}S_{lj}\,\mathrm{Cov}[X_i, X_j] = \bigl(S V_X S^{\mathsf T}\bigr)_{kl} ,

the last equality being the definition of the matrix product: $(SV_XS^{\mathsf T})_{kl} = \sum_i\sum_j S_{ki}(V_X)_{ij}(S^{\mathsf T})_{jl}$ and $(S^{\mathsf T})_{jl} = S_{lj}$ .

Step 5 (the single-output form). Put $m = 1$ , write $f$ for $f_1$ and split the double sum into its diagonal and off-diagonal parts, using $(V_X)_{ii} = V[X_i]$ :

V[Y] \;\approx\; \sum_i \left(\frac{\partial f}{\partial x_i}\right)^{\!2}V[X_i] \;+\; \sum_{i\ne j}\frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\mathrm{Cov}[X_i,X_j] .

The off-diagonal sum runs over ordered pairs, so each unordered pair $\{i,j\}$ appears twice and contributes $2\,\partial_if\,\partial_jf\,\mathrm{Cov}[X_i,X_j]$ . $\square$

Positive semi-definiteness is inherited. For any real $t \in \mathbb{R}^m$ , $t^{\mathsf T}V_Yt = (S^{\mathsf T}t)^{\mathsf T}V_X(S^{\mathsf T}t) \ge 0$ , so $V_Y$ is again a legitimate covariance matrix. A propagation scheme that dropped the cross terms would not have this property, and can return a "variance" that is larger than the true one by an arbitrary factor, as the second example shows.

Remark(Where the hypotheses are used, and the size of what was thrown away)Original

(P1) is what makes the question well posed. If an input has no finite variance, the output's error bar is not a number that exists, and the formula returns a finite answer for a quantity that is infinite. The standard trap is a ratio: if $X$ and $Y$ are independent $N(0,1)$ then $X/Y$ is Cauchy and has no variance at all (Random variables and distributions), yet the rule applied at any point with $\mu_Y \ne 0$ hands back a finite number. The practical version of this is a quotient whose denominator has a non-negligible chance of being near zero.

(P2) is used only at the single point $\mu$ . The rule needs no smoothness away from the means — but it needs genuine differentiability there. A function with a kink at $\mu$ , such as $f(x) = |x|$ at $\mu = 0$ , has no tangent plane and the rule is meaningless: for $X \sim N(0, \sigma^2)$ the exact answer is $V[|X|] = (1 - 2/\pi)\sigma^2 \approx 0.363\,\sigma^2$ , which no choice of derivative reproduces.

(P3) is the real hypothesis, and it is quantitative. Keeping the quadratic term of (2) for a single input gives the two corrections that linearisation discards:

\mathbb{E}[Y] \;\approx\; f(\mu) + \tfrac12 f''(\mu)\,\sigma^2 , \qquad V[Y] \;\approx\; f'(\mu)^2\sigma^2 + \tfrac12 f''(\mu)^2\sigma^4 ,

the first of which is the correction $\tfrac12\sigma^2f''(\bar x)$ that Arfken, Weber and Harris exhibit in (23.87); in $n$ variables the bias is $\tfrac12\sum_{i,j}\partial_i\partial_jf\,(V_X)_{ij}$ . The ratio of the neglected bias to the propagated error is $\tfrac12|f''|\sigma^2 \big/ |f'|\sigma = \tfrac12\sigma\,|f''/f'|$ , so the linearisation is safe exactly when $\sigma$ is small compared with the scale $|f'/f''|$ on which the slope of $f$ changes — not merely small compared with $\mu$ .

A case where it fails, with exact numbers. Let $X \sim N(1, 0.5^2)$ and $Y = X^2$ . Exactly, $\mathbb{E}[Y] = \mu^2 + \sigma^2 = 1.25$ and $V[Y] = 4\mu^2\sigma^2 + 2\sigma^4 = 1.125$ , so $Y = 1.25 \pm 1.0607$ . Linearisation gives $f(\mu) = 1$ and $|f'(\mu)|\sigma = 1.0$ : the error bar is only $6\%$ low, but the central value is $25\%$ low, and it is low systematically — averaging more experiments will not remove it. Curvature biases the answer before it spoils the error bar.

Dropping the covariances is not conservative. Correlations can either inflate or deflate the output variance, and which one happens depends on the signs of the sensitivities. For $Y = X_1 - X_2$ with equal variances $\sigma^2$ and correlation $\rho$ , $V[Y] = 2\sigma^2(1-\rho)$ , which for $\rho \to 1$ tends to zero: a difference of two strongly correlated measurements is far better determined than either. A shared calibration constant produces exactly this, which is why a ratio of two intensities from the same instrument can beat the error bar on each.

Using $x$ in place of $\mu$ . The sensitivities must be evaluated at the unknown $\mu$ and in practice are evaluated at the measured $x$ . The resulting error in $S$ is itself of first order in $X - \mu$ , so it perturbs $V_Y$ at second order — the same order as the terms already discarded, and no worse.

Relative errors for products. If $Y = X_1^{p_1}\cdots X_n^{p_n}$ with uncorrelated positive inputs, then $\partial_i f/f = p_i/\mu_i$ , and dividing the propagation formula by $f(\mu)^2$ gives

\left(\frac{\sigma_Y}{Y}\right)^{\!2} \;\approx\; \sum_{i=1}^{n} p_i^2\left(\frac{\sigma_i}{\mu_i}\right)^{\!2} ,

the rule that relative errors add in quadrature with the exponents as weights. The cases $C = AB$ and $C = A/B$ are Arfken, Weber and Harris's Exercise 23.7.1, and the two give the same relative error, since the sign of $p$ is squared away.

Method(Propagating uncertainty)AdaptedAdapted from Riley, Hobson & Bence (2006), §30.12.3, the route from (30.137) to (30.141), pp. 1203-1205

Write the output as an explicit function of the measured inputs, $Y_k = f_k(x_1,\dots,x_n)$ , and collect the inputs' error bars and correlations into the covariance matrix $V_X$ . Any quantity measured by the same instrument, read from the same calibration curve, or fitted from the same data set is correlated with its companions until shown otherwise.
Check for curvature and for zeros in denominators. Compare each $\sigma_i$ with the scale $|f'/f''|$ in that variable (remark above). If they are comparable, or if a denominator can plausibly vanish, stop: linearisation is not applicable and the distribution of $Y$ must be obtained directly, by the change-of-variable formula of Joint distributions and transformations or by Monte Carlo.
Differentiate at the measured point to get $S_{ki} = \partial f_k/\partial x_i$ , evaluated at $x$ .
Form $V_Y = SV_XS^{\mathsf T}$ . For one output this is the single-output formula; do not drop the cross terms. For several outputs keep the whole matrix, because the outputs are now correlated even if the inputs were not, and anything computed downstream from them will need those correlations.
For a product or quotient of powers, use relative errors and add $p_i^2(\sigma_i/\mu_i)^2$ in quadrature; it is faster and less error-prone than differentiating.
Quote $Y_k = f_k(x) \pm \sqrt{(V_Y)_{kk}}$ , rounding the error to one or two significant figures and the value to the same decimal place. Keep the full $V_Y$ if the result will be used further.
Sanity-check. Each term of the sum is a contribution to the error budget; identify the dominant one. If one input dominates by an order of magnitude in variance, the others are irrelevant and improving them is wasted effort.

Example(A pendulum measurement of

g

, and which input limits it)Original

A simple pendulum of length $L$ has small-oscillation period $T = 2\pi\sqrt{L/g}$ , so

g \;=\; \frac{4\pi^2 L}{T^2} .

Measurements: $L = 1.000 \pm 0.002\ \mathrm{m}$ (a metre rule read to the nearest millimetre, with the uncertainty dominated by locating the centre of mass of the bob), and the period obtained by timing $20$ swings with a stopwatch, giving $T = 2.006 \pm 0.004\ \mathrm{s}$ . The two are measured by different instruments and are taken to be uncorrelated.

Step 1: central value.

g \;=\; \frac{4\pi^2 (1.000)}{(2.006)^2} \;=\; \frac{39.4784}{4.024036} \;=\; 9.8107\ \mathrm{m\,s^{-2}} .

Step 2: curvature check. In $L$ the function is linear, so there is nothing to check. In $T$ , $f \propto T^{-2}$ gives $|f'/f''| = T/3 = 0.6687\ \mathrm{s}$ , against $\sigma_T = 0.004\ \mathrm{s}$ . The ratio of the neglected bias to the propagated error is $\tfrac12\sigma_T|f''/f'| = 0.5 \times 0.004/0.6687 = 3.0\times10^{-3}$ , three parts in a thousand of the error bar, so linearisation is amply safe; in absolute terms the bias is $\tfrac12 f''\sigma_T^2 = 3g\sigma_T^2/T^2 = 1.2\times10^{-4}\ \mathrm{m\,s^{-2}}$ .

Step 3: relative errors. $g = 4\pi^2 L^{1}T^{-2}$ is a product of powers with $p_L = 1$ and $p_T = -2$ , so by the rule of the remark

\left(\frac{\sigma_g}{g}\right)^{\!2} \;=\; \left(\frac{\sigma_L}{L}\right)^{\!2} + 4\left(\frac{\sigma_T}{T}\right)^{\!2} \;=\; (2.000\times10^{-3})^2 + 4\,(1.9940\times10^{-3})^2 .

The two contributions are $4.000\times10^{-6}$ and $1.5904\times10^{-5}$ : the timing term is $3.98$ times the length term, so the measurement is limited by the stopwatch: halving $\sigma_L$ would shrink $\sigma_g$ only from $4.461\times10^{-3}g$ to $4.111\times10^{-3}g$ , an improvement of $7.8\%$ , whereas halving $\sigma_T$ would nearly halve it.

Step 4: the answer. $\sigma_g/g = \sqrt{1.9904\times10^{-5}} = 4.4614\times10^{-3}$ , hence $\sigma_g = 9.8107 \times 4.4614\times10^{-3} = 0.0438\ \mathrm{m\,s^{-2}}$ , and

g \;=\; 9.811 \pm 0.044\ \mathrm{m\,s^{-2}} .

The standard value is $9.80665\ \mathrm{m\,s^{-2}}$ , so the discrepancy is $+0.0041\ \mathrm{m\,s^{-2}}$ , or $0.09$ standard errors: the measurement and the accepted value agree.

Step 5: what to improve. To halve $\sigma_g$ one must halve $\sigma_T$ , which means timing $40$ swings instead of $20$ at the same reaction-time uncertainty; no improvement in the metre rule will do it.

Example(A correlated pair: predicting from a straight-line fit)Original

Five measurements $y_i$ are made at $x_i = 1,2,3,4,5$ , each with independent error $\sigma = 0.3$ , and a straight line $y = c + mx$ is fitted by least squares (Least-squares fitting). With design matrix $R$ whose rows are $(1, x_i)$ and $V_X = \sigma^2 I$ , the formula $V = (R^{\mathsf T}N^{-1}R)^{-1}$ gives the covariance matrix of $(\hat c, \hat m)$ :

R^{\mathsf T}R = \begin{pmatrix} 5 & 15 \\ 15 & 55\end{pmatrix} , \qquad (R^{\mathsf T}R)^{-1} = \frac{1}{50}\begin{pmatrix} 55 & -15 \\ -15 & 5\end{pmatrix} ,

V_{(\hat c,\hat m)} \;=\; \sigma^2 (R^{\mathsf T}R)^{-1} \;=\; \begin{pmatrix} 0.0990 & -0.0270 \\ -0.0270 & 0.0090 \end{pmatrix} .

So $\sigma_{\hat c} = 0.3146$ , $\sigma_{\hat m} = 0.0949$ , and the correlation is $-0.0270/(0.3146\times0.0949) = -0.905$ : slope and intercept are strongly anti-correlated, because a fit that is steeper must start lower to keep passing through the middle of the data.

The predicted value at $x = 3$ . Put $Y = f(\hat c, \hat m) = \hat c + 3\hat m$ , an affine function, so the rule is exact. The sensitivity row is $S = (1, \ 3)$ and

V[Y] \;=\; S V S^{\mathsf T} \;=\; 0.0990 + 9(0.0090) + 2(3)(-0.0270) \;=\; 0.0990 + 0.0810 - 0.1620 \;=\; 0.0180 ,

giving $\sigma_Y = \sqrt{0.0180} = 0.1342$ .

The check. $x = 3$ is the mean of the $x_i$ , and at the centroid a least-squares line is pinned by the sample mean of the $y_i$ alone, whose standard error is $\sigma/\sqrt{N} = 0.3/\sqrt{5} = 0.13416$ . The propagated value agrees to all five figures, which it must, since both are exact statements about the same affine function.

What dropping the cross term would have cost. Keeping only the squared-derivative terms gives $0.0990 + 0.0810 = 0.1800$ and $\sigma_Y = 0.4243$ — too large by a factor of $3.16$ . An error bar wrong by a factor of three in the safe direction is still wrong: it would hide a real disagreement with a prediction.

The same calculation away from the centroid. At $x = 6$ , one unit beyond the data, $V[Y] = 0.0990 + 36(0.0090) + 12(-0.0270) = 0.0990$ and $\sigma_Y = 0.3146$ . The cancellation is now only partial, and the extrapolated point is $2.3$ times less certain than the interpolated one — the quantitative form of the warning against extrapolating a fit.

Application(Error budgets, fitted parameters and filters)Original

The error budget of an experiment. Writing out the terms of $\sum_i (\partial_i f)^2\sigma_i^2$ before any data are taken is how an experiment is designed: each term is one line of a budget, and the design effort goes where the budget is largest. The pendulum example is the elementary case; in a precision measurement the same table runs to dozens of lines and separates statistical from systematic entries, the latter being exactly the correlated contributions, since a systematic error is one shared between measurements and therefore carried by an off-diagonal entry of $V_X$ .

Fitted parameters and everything computed from them. Riley, Hobson and Bence derive the least-squares covariance matrix $V = (R^{\mathsf T}N^{-1}R)^{-1}$ by applying $V_Y = SV_XS^{\mathsf T}$ to the estimator $\hat a = Sy$ , a case in which the function is linear and the rule is therefore exact (Riley, Hobson & Bence (2006), §31.6, eqs. (31.98)-(31.99), pp. 1274-1275). Any quantity derived from fitted parameters — a predicted value, a ratio of amplitudes, an extrapolated intercept — must then be propagated through the full matrix, as the second example does. This is also why fitting codes return a covariance matrix rather than a list of error bars, and why discarding the off-diagonal entries silently corrupts everything downstream. See Least-squares fitting and Maximum likelihood, whose asymptotic covariance is the inverse Fisher matrix.

Linear filters and data assimilation. The forecast step of the Kalman filter is $V_Y = SV_XS^{\mathsf T} + Q$ : the linear propagation law of this unit applied to the model's state-transition matrix, plus the model error. The extended Kalman filter is the same formula with $S$ the Jacobian of a non-linear model, which is to say it is this unit's approximation, with this unit's failure mode — in a chaotic system the curvature scale $|f'/f''|$ shrinks as the forecast lengthens, and the linearised covariance becomes untrustworthy long before the forecast does. The same structure appears in Optimal interpolation and variational data assimilation.

Where the rule is abandoned. When the curvature test of step 2 fails, the modern practice is to sample: draw inputs from their joint distribution, push each draw through $f$ , and read off the spread of the outputs. That costs nothing in accuracy and is exact in the limit of many draws, and the formula of this unit survives as the thing one checks the sample against.

Method

Confidence intervals

Adapted

Add to PDF

Builds interval estimates by the pivotal method, gives the exact normal-sample intervals for a mean (known and unknown variance) and for a variance, says precisely what the confidence level asserts and what it does not, and works nine speed-of-sound readings to $\mu \in [342.62, 344.18]\,\mathrm{m\,s^{-1}}$ and $\sigma \in [0.684, 1.940]\,\mathrm{m\,s^{-1}}$ at the 95% level.

Prerequisites: Sampling distributions, Estimators and their properties, Random variables and distributions, Expectation, variance and covariance, The central limit theorem

Note(What an interval estimate is for)Original

A point estimate answers the question "what is the value?" and is almost always wrong: a continuous parameter is hit with probability zero. An interval estimate answers the more honest question "what range is consistent with the data, and at what rate does a procedure of this kind succeed?" The machinery for producing such a range is the sampling distribution of Sampling distributions, used backwards. Forwards, it says where the estimator lands when the parameter is known; backwards, it says which parameter values would have made the observed estimate unremarkable.

Everything in this unit is a frequentist construction: the interval is random and the parameter is fixed. The rival construction, in which the parameter is given a distribution and the interval is read off it, is Bayesian inference; the two agree numerically far more often than they agree in meaning.

Definition(Quantiles, and the notation used below)Original

For a random variable with a continuous, strictly increasing distribution function $F$ and for $0 < p < 1$ , the $p$ -quantile is the unique number $q$ with $F(q) = p$ . We write

z_p = \Phi^{-1}(p), \qquad t_{\nu,p}, \qquad \chi^{2}_{\nu,p}

for the $p$ -quantiles of, respectively, the standard normal distribution, Student's $t$ on $\nu$ degrees of freedom and the chi-squared distribution on $\nu$ degrees of freedom, all three defined in Sampling distributions. Symmetry of the first two gives $z_{1-p} = -z_p$ and $t_{\nu,1-p} = -t_{\nu,p}$ ; the chi-squared distribution is not symmetric, and its two tails must be looked up separately. The numbers used below are $z_{0.975} = 1.95996$ , $t_{8,0.975} = 2.30600$ , $\chi^{2}_{8,0.025} = 2.17973$ and $\chi^{2}_{8,0.975} = 17.53455$ . Riley, Hobson and Bence tabulate the last two families in their tables 31.3 and 31.2.

Definition(Confidence interval, confidence level, confidence limits)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.3.4, eqs. (31.25)-(31.29) and the naming paragraph, pp. 1235-1237

Let the data $x$ come from a population depending on one unknown parameter $a$ , all other parameters being known, and let $\hat a$ be an estimator of $a$ with sampling distribution $P(\hat a \mid a)$ . Fix two tail probabilities $\alpha, \beta \ge 0$ with $\alpha + \beta < 1$ . For each value of $a$ , Riley, Hobson and Bence define the two cut points $\hat a_\alpha(a)$ and $\hat a_\beta(a)$ by

\Pr\bigl(\hat a < \hat a_\alpha(a) \mid a\bigr) = \alpha , \qquad \Pr\bigl(\hat a > \hat a_\beta(a) \mid a\bigr) = \beta ,

so that $\hat a$ falls between them with probability $1 - \alpha - \beta$ . Assume both cut points are increasing functions of $a$ - which is what " $\hat a$ estimates $a$ " ought to mean, and which is needed for the next line. Then, given the observed value $\hat a_{\mathrm{obs}}$ , the two equations

\hat a_\alpha(a_+) = \hat a_{\mathrm{obs}} , \qquad \hat a_\beta(a_-) = \hat a_{\mathrm{obs}}

have unique solutions $a_-$ and $a_+$ , and they satisfy

\Pr(a_- < a < a_+) \;=\; 1 - \alpha - \beta .

The random interval $[a_-, a_+]$ is the confidence interval on $a$ at confidence level $1 - \alpha - \beta$ (usually quoted as a percentage); $a_-$ and $a_+$ are the lower and upper confidence limits. Taking $\alpha = \beta$ gives the central interval; taking $\beta = 0$ , so that only $a_+$ is finite, gives a one-sided interval, which is what an upper or lower bound on a parameter means.

Arfken, Weber and Harris describe the same object less formally, as the range that contains a chosen fraction $p$ of the probability, reported in the form $x = x_0 \pm dx$ "with $(100p)\%$ confidence", and they name the technique used throughout this unit the pivotal method: relate the unknown to a quantity whose distribution is known and free of unknowns, and read the interval off that.

Method(Exact intervals for a normal sample)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.3.5, eqs. (31.30)-(31.33), pp. 1238-1239; Riley, Hobson & Bence (2006), §31.7.5, eqs. (31.113) and (31.115) and the inversion on p. 1288; Riley, Hobson & Bence (2006), §31.7.6, eq. (31.122) and the equal-tail paragraph, p. 1292; Arfken, Weber & Harris (2013), §23.7, eqs. (23.116)-(23.118) and table 23.3, pp. 1176-1177

Let $X_1, \dots, X_n$ be independent $N(\mu, \sigma^{2})$ variables, $n \ge 2$ , with sample mean $\overline X$ and sample variance $S^{2} = (n-1)^{-1}\sum_i (X_i - \overline X)^{2}$ . Fix $0 < \alpha < 1$ .

1. Mean, variance known. The interval

\left[\, \overline X - z_{1-\alpha/2}\frac{\sigma}{\sqrt n} \;,\; \overline X + z_{1-\alpha/2}\frac{\sigma}{\sqrt n} \,\right]

covers $\mu$ with probability exactly $1 - \alpha$ . More generally, with tails $\alpha$ and $\beta$ , the limits are $\overline X - z_{1-\beta}\,\sigma/\sqrt n$ and $\overline X + z_{1-\alpha}\,\sigma/\sqrt n$ .

2. Mean, variance unknown. The interval

\left[\, \overline X - t_{n-1,\,1-\alpha/2}\frac{S}{\sqrt n} \;,\; \overline X + t_{n-1,\,1-\alpha/2}\frac{S}{\sqrt n} \,\right]

covers $\mu$ with probability exactly $1 - \alpha$ , for every $n \ge 2$ and every value of the unknown $\sigma$ .

3. Variance. With $a = \chi^{2}_{n-1,\,\alpha/2}$ and $b = \chi^{2}_{n-1,\,1-\alpha/2}$ , the interval

\left[\, \frac{(n-1)S^{2}}{b} \;,\; \frac{(n-1)S^{2}}{a} \,\right]

covers $\sigma^{2}$ with probability exactly $1 - \alpha$ ; taking square roots gives an interval for $\sigma$ of the same level.

4. Large samples, any population. If $\hat a$ is an estimator whose sampling distribution tends to $N(a, \sigma_{\hat a}^{2})$ as $n \to \infty$ - which the central limit theorem of The central limit theorem makes the usual case - then $\hat a \pm z_{1-\alpha/2}\,\sigma_{\hat a}$ has asymptotic coverage $1 - \alpha$ , and an estimate $\hat\sigma_{\hat a}$ may be substituted for $\sigma_{\hat a}$ when the latter is unknown.

In particular, since $z_{0.8413} = 1$ , quoting a result as $\hat a \pm \sigma_{\hat a}$ - the physicist's one-sigma error bar - is the central confidence interval at level $1 - 2(1 - 0.8413) = 0.683$ .

Idea of the proofOriginal

All three exact intervals are one trick applied three times. Find a pivot: a function of the data and of the single unknown parameter whose distribution involves no unknowns at all. Then a probability statement about the pivot can be written down from a table, and because the pivot is monotone in the parameter the statement can be rearranged, with the parameter alone in the middle. The rearrangement is pure algebra; all the probability was spent in the first step.

For the mean with $\sigma$ known the pivot is $\sqrt n(\overline X - \mu)/\sigma$ , standard normal. For the mean with $\sigma$ unknown, $\sigma$ is cancelled by dividing instead by $S$ , and the pivot becomes Student's $t$ - this is exactly why the $t$ distribution was invented. For the variance the pivot is $(n-1)S^{2}/\sigma^{2}$ , chi-squared; being a decreasing function of $\sigma^{2}$ , it inverts with the two tails exchanged.

Proof(The three exact intervals)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.3.5, the inversion (31.31)-(31.32), p. 1238, and §31.7.5, the inversion of the two-tailed $t$-test, p. 1288; Riley, Hobson & Bence (2006), §31.7.6, eq. (31.122) and the equal-tail paragraph, p. 1292, inverted here by the same route

Throughout we use the two facts proved in Sampling distributions for a normal sample: $\overline X \sim N(\mu, \sigma^{2}/n)$ , and $\overline X$ and $S^{2}$ are independent with $(n-1)S^{2}/\sigma^{2} \sim \chi^{2}(n-1)$ .

1. Mean, variance known. Standardising the sample mean,

Z = \frac{\overline X - \mu}{\sigma/\sqrt n} \sim N(0,1) ,

and $Z$ is a pivot: it contains the unknown $\mu$ and nothing else unknown. By the definition of the quantiles and the symmetry $z_{\alpha/2} = -z_{1-\alpha/2}$ ,

\Pr\bigl(-z_{1-\alpha/2} \le Z \le z_{1-\alpha/2}\bigr) = 1 - \alpha .

Multiplying through by $\sigma/\sqrt n > 0$ and subtracting $\overline X$ reverses the inequalities and gives, after a change of sign, the event $\overline X - z_{1-\alpha/2}\sigma/\sqrt n \le \mu \le \overline X + z_{1-\alpha/2}\sigma/\sqrt n$ . The two events are the same event, so they have the same probability. The unequal-tail version is the same computation with $z_{1-\alpha}$ and $z_{1-\beta}$ in the two places.

2. Mean, variance unknown. Write

T = \frac{\overline X - \mu}{S/\sqrt n} = \frac{(\overline X - \mu)/(\sigma/\sqrt n)}{\sqrt{\bigl[(n-1)S^{2}/\sigma^{2}\bigr]/(n-1)}} ,

where the second form is obtained by dividing numerator and denominator by $\sigma$ . The numerator is $N(0,1)$ and the quantity in the square root is $\chi^{2}(n-1)/(n-1)$ ; the independence of $\overline X$ and $S^{2}$ is what allows the definition of Student's $t$ to be applied, and it gives $T \sim t(n-1)$ . Note that $\sigma$ has cancelled: $T$ is a pivot even though $\sigma$ is unknown. Now repeat step 1 with $t_{n-1,1-\alpha/2}$ in place of $z_{1-\alpha/2}$ and $S$ in place of $\sigma$ .

3. Variance. Put $W = (n-1)S^{2}/\sigma^{2} \sim \chi^{2}(n-1)$ , again a pivot. By the definition of $a$ and $b$ as the $\alpha/2$ and $1 - \alpha/2$ quantiles, the two tails each carry probability $\alpha/2$ , so

\Pr(a \le W \le b) = 1 - \alpha .

Since $S^{2} > 0$ with probability one (the sample values are almost surely not all equal, because the normal distribution is continuous) and $a > 0$ , each of the two inequalities may be solved for $\sigma^{2}$ :

a \le \frac{(n-1)S^{2}}{\sigma^{2}} \iff \sigma^{2} \le \frac{(n-1)S^{2}}{a} , \qquad \frac{(n-1)S^{2}}{\sigma^{2}} \le b \iff \frac{(n-1)S^{2}}{b} \le \sigma^{2} .

Taking the two together converts the event $\{a \le W \le b\}$ into the event that $\sigma^{2}$ lies in the stated interval, and the probability is unchanged. The map $\sigma \mapsto \sigma^{2}$ is increasing on $(0,\infty)$ , so taking square roots of the endpoints gives an interval for $\sigma$ with the same coverage.

4. Large samples. By hypothesis $(\hat a - a)/\sigma_{\hat a}$ converges in distribution to $N(0,1)$ , so the probability of the event in step 1 converges to $1 - \alpha$ ; the conclusion is asymptotic, not exact. $\square$

Remark(Where the hypotheses are used, and what fails without them)Original

Monotonicity of $\hat a_\alpha$ and $\hat a_\beta$ in $a$ is what makes $a_-$ and $a_+$ well defined in the general construction. Without it the set of parameter values consistent with the data need not be an interval at all, and the honest output of the procedure is that set rather than a pair of limits.

Independence of $\overline X$ and $S^{2}$ is used once, in step 2, and it is special to the normal distribution: it is equivalent to normality for a sample of size $n \ge 2$ from a non-degenerate distribution with finite variance. If the parent is not normal, $T$ is not $t$ -distributed and the interval is only approximate.

Normality is cheap for interval 2 and expensive for interval 3. For the mean, the central limit theorem rescues the construction: $\sqrt n(\overline X - \mu)/S$ tends to $N(0,1)$ for any parent with finite variance, and since $t_{\nu,p} \to z_p$ the $t$ interval is asymptotically correct whatever the parent. For the variance there is no such rescue. Applying the central limit theorem to the variables $(X_i - \mu)^{2}$ gives $\sqrt n\,(S^{2} - \sigma^{2}) \to N(0, \mu_4 - \sigma^{4})$ , where $\mu_4$ is the fourth central moment of the parent, whereas the chi-squared calibration presumes the normal value $\mu_4 = 3\sigma^{4}$ , that is an asymptotic variance $2\sigma^{4}$ . For an exponential parent $\mu_4 = 9\sigma^{4}$ , so the true standard deviation of $\sqrt n\,(S^{2} - \sigma^{2})$ is $\sqrt{8\sigma^{4}} = 2\sqrt{2}\,\sigma^{2}$ against the assumed $\sqrt{2}\,\sigma^{2}$ : the interval is too short by a factor of two, and its coverage tends to

2\Phi\!\left(\frac{1.95996}{2}\right) - 1 = 2(0.83645) - 1 = 0.6729 ,

not $0.95$ , however large $n$ is. A nominal 95% interval for a variance is worth little unless normality has been checked.

Independence of the observations is used in every step, through the sampling distributions. Positively correlated data - successive readings of a drifting instrument, or a time series - make $S^{2}$ underestimate the variance of $\overline X$ , and all the intervals above come out too narrow. The repair is to replace $n$ by an effective sample size computed from the autocorrelation; see Stationary processes and spectral densities.

One parameter at a time. Intervals 2 and 3 are each exact for their own parameter; the rectangle formed by taking both at level $1-\alpha$ is not a $1-\alpha$ region for the pair $(\mu, \sigma^{2})$ . For several parameters at once, Riley, Hobson and Bence construct a region bounded by the ellipsoid $Q(a, \hat a_{\mathrm{obs}}) = c$ , where $Q(\hat a, a) = (\hat a - a)^{\mathsf T}V^{-1}(\hat a - a)$ is the quadratic form of the large-sample multivariate Gaussian and $c$ is the $1-\alpha$ quantile of $\chi^{2}_{M}$ for $M$ parameters.

Remark(What the confidence level says, and what it does not)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.3.4, the paragraph emphasising that $a_-$ and $a_+$ are random variables, p. 1237; Arfken, Weber & Harris (2013), §23.7, the closing observations of Example 23.7.3, p. 1178

The probability statement is about the interval, not about the parameter. Riley, Hobson and Bence put the emphasis exactly there: $a_-$ and $a_+$ are random variables, and the content of "95% confidence" is that if a large number of samples of the same size were drawn and the interval computed each time, the true value would lie inside on 95% of those occasions. Once numbers are substituted, $[342.62, 344.18]$ either contains the true speed of sound or does not; there is no 95% left in it. A sentence of the form "there is a 95% probability that $\mu$ lies in this interval" is a statement about a posterior distribution for $\mu$ , and needs the prior of Bayesian inference before it means anything.

Three practical corollaries, the first two of which Arfken, Weber and Harris stress:

A higher level buys a wider interval. The level is chosen, not discovered; quoting the level without the interval, or the interval without the level, is meaningless.
The interval is itself sample-dependent. Another sample from the same population gives different endpoints and a different width, the width varying because $S$ does.
Coverage is coverage under the assumed model. A confidence interval answers for statistical fluctuation only. A miscalibrated instrument shifts every sample the same way, and no amount of data narrows that error; systematic uncertainties are estimated and reported separately.

Remark(Three ways the

t

interval is written, and which we follow)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.7.5, eqs. (31.113) and (31.115) and the confidence interval displayed on p. 1288; Arfken, Weber & Harris (2013), §23.7, eqs. (23.95), (23.116)-(23.118) and Example 23.7.3, pp. 1168, 1176-1178

The half-width of interval 2 is written differently by our two sources, because they use different symbols for the sample spread, and the reader who mixes them will get the wrong number.

This unit uses $S^{2} = (n-1)^{-1}\sum_i(X_i-\overline X)^{2}$ and the half-width $t_{n-1,1-\alpha/2}\,S/\sqrt n$ .
Riley, Hobson and Bence use $s^{2} = n^{-1}\sum_i(X_i-\overline X)^{2}$ , define $t = (\bar x - \mu_0)/\bigl(s/\sqrt{N-1}\bigr)$ in their eq. (31.113), and invert it into $\bar x \pm t_{\mathrm{crit}}\,s/\sqrt{N-1}$ . Since $s/\sqrt{N-1} = S/\sqrt N$ identically, this is the same interval as ours, written in the other convention for the sample variance; their own eq. (31.115) says as much, giving $t = (\bar x - \mu_0)/(\hat\sigma/\sqrt N)$ with $\hat\sigma^{2} = Ns^{2}/(N-1) = S^{2}$ .
Arfken, Weber and Harris take $\sigma$ to be the sample standard deviation with divisor $n-1$ , as in their eq. (23.95), and write the interval through eq. (23.118) as $\mu = \bar x - T\sigma/\sqrt n$ . In their Example 23.7.3 the six readings are treated with $n = 5$ , the number of degrees of freedom, in both the divisor of $\sigma^{2}$ and the $\sqrt n$ of the half-width, which yields $5.885 \pm 1.720$ at the 90% level. Our rule applied to their numbers divides by $\sqrt 6$ instead and gives $5.885 \pm 1.570$ .

We follow Riley, Hobson and Bence: the degrees of freedom $n-1$ fix which $t$ distribution is used, and the sample size $n$ fixes the standard error $S/\sqrt n$ . The difference is recorded in the handbook's list of source discrepancies and is not a defect of either book's statement of the $t$ density, on which they agree.

Method(Constructing an interval by the pivotal method)Original

Write down the model and identify the single unknown parameter $\theta$ the interval is for; everything else must be known or eliminated.
Find a pivot $G(\text{data}, \theta)$ whose distribution is completely known - free of $\theta$ and of every nuisance parameter. For normal samples use the table in the statement above; otherwise the sampling distribution of Sampling distributions or the large-sample normal approximation of item 4 is the place to look.
Choose the level $1-\alpha$ and the split of the tails. Equal tails, $\alpha/2$ each, is the default; put the whole of $\alpha$ in one tail if what is wanted is a bound rather than a range.
Look up the two quantiles of the pivot's distribution and write the bracketing statement $\Pr(q_{\alpha/2} \le G \le q_{1-\alpha/2}) = 1-\alpha$ .
Solve both inequalities for $\theta$ , reversing them if the pivot decreases in $\theta$ . The result is the interval.
Substitute the observed numbers, and report the level, the interval and the model assumed. Check the arithmetic by confirming that the point estimate lies inside.

Example(Nine readings of the speed of sound)Original

A student measures the speed of sound in air at $20\,^\circ\mathrm{C}$ by timing a pulse over a fixed baseline, nine times, and obtains (in $\mathrm{m\,s^{-1}}$ )

342.6,\; 344.1,\; 343.2,\; 341.8,\; 345.0,\; 342.9,\; 343.7,\; 344.5,\; 342.8 .

Take the readings to be independent and normally distributed about the true value $\mu$ with unknown standard deviation $\sigma$ . Find 95% central confidence intervals for $\mu$ and for $\sigma$ .

Step 1: the summary statistics. The nine values sum to $3090.6$ , so

\overline x = \frac{3090.6}{9} = 343.400 .

The deviations $x_i - \overline x$ are

-0.8,\; 0.7,\; -0.2,\; -1.6,\; 1.6,\; -0.5,\; 0.3,\; 1.1,\; -0.6 ,

which sum to zero, as they must. Their squares sum to

\sum_i (x_i - \overline x)^{2} = 0.64 + 0.49 + 0.04 + 2.56 + 2.56 + 0.25 + 0.09 + 1.21 + 0.36 = 8.20 ,

so $s^{2} = 8.20/8 = 1.025$ and $s = 1.01242$ . The standard error of the mean is $s/\sqrt 9 = 1.01242/3 = 0.337474$ .

Step 2: the interval for $\mu$ . Here $\sigma$ is unknown, so the pivot is $T = (\overline X - \mu)/(S/\sqrt 9) \sim t(8)$ . With $\alpha = 0.05$ the quantile is $t_{8,0.975} = 2.30600$ , and the half-width is

2.30600 \times 0.337474 = 0.778217 .

Hence the 95% central confidence interval for $\mu$ is

343.400 \pm 0.778 = [\,342.622,\; 344.178\,]\ \mathrm{m\,s^{-1}} .

Step 3: the interval for $\sigma$ . The pivot is $W = 8S^{2}/\sigma^{2} \sim \chi^{2}(8)$ , and the quantiles are $a = \chi^{2}_{8,0.025} = 2.17973$ and $b = \chi^{2}_{8,0.975} = 17.53455$ . With $\sum_i (x_i - \overline x)^{2} = 8.20$ ,

\frac{8.20}{17.53455} = 0.467648 , \qquad \frac{8.20}{2.17973} = 3.761933 ,

so the 95% interval for $\sigma^{2}$ is $[0.4676,\; 3.7619]\ \mathrm{m^{2}\,s^{-2}}$ and, on taking square roots,

\sigma \in [\,0.684,\; 1.940\,]\ \mathrm{m\,s^{-1}} .

Step 4: read the answers. The interval for $\mu$ is nearly symmetric about the estimate, as the $t$ distribution is symmetric; the interval for $\sigma$ is strongly asymmetric about $s = 1.012$ , running from two thirds of it to nearly twice it. Nine readings pin a mean to a fraction of a metre per second and pin a standard deviation hardly at all - the familiar fact that variances need far more data than means. For comparison, had $\sigma = 1.000$ been known in advance, interval 1 would have given $343.400 \pm 1.95996/3 = 343.400 \pm 0.653$ , narrower than the $t$ interval because nothing had to be spent estimating the spread. The accepted value $343.2\ \mathrm{m\,s^{-1}}$ at $20\,^\circ\mathrm{C}$ lies inside the interval for $\mu$ .

Application(Error bars, tolerances and detection limits)Original

The error bar on a plot. A point drawn as $y \pm \sigma_y$ is a 68.3% central confidence interval, by the calculation in the statement above. This is why error bars from different fields can be compared at all, and why a paper that draws two-sigma bars must say so: the same picture then means 95.4%.

Tolerance and qualification. A supplier of resistors claims a spread no worse than some $\sigma_{\max}$ . The quantity under test is a variance, so the relevant object is interval 3, or rather its one-sided version: the upper confidence limit $\sqrt{(n-1)S^{2}/\chi^{2}_{n-1,\alpha}}$ is the number that may honestly be quoted as a bound on $\sigma$ at level $1-\alpha$ . The strong asymmetry seen in the worked example is why qualification tests need large batches.

Detection limits and upper bounds. When an experiment sees no signal, the result is reported as a one-sided interval: the 90% or 95% upper confidence limit on the rate or cross-section. Half a century of particle-physics and astroparticle limits are exactly this construction, with the Poisson rather than the normal sampling distribution.

Trends in geophysical records. A fitted trend in a temperature or sea-level series is reported with a confidence interval obtained from the fit of Least-squares fitting. Because the residuals of such a series are autocorrelated rather than independent, the effective sample size is much smaller than the number of months in the record, and an interval computed as if the data were independent is the standard way to overstate the significance of a trend.

Method

Hypothesis testing

Adapted

Add to PDF

Sets up null and alternative hypotheses, rejection regions, significance, power and the $p$ -value, proves the Neyman-Pearson lemma, derives Student's $t$ -test as a generalised likelihood-ratio test and gives the chi-squared goodness-of-fit test, with eight dark-current readings giving $t = 3.000$ and $p = 0.0199$ and five count bins giving $\chi^{2} = 4.625$ on four degrees of freedom.

Prerequisites: Confidence intervals, Sampling distributions, Maximum likelihood, Least-squares fitting, The central limit theorem

Note(What a test is, and what it is not)Original

An interval estimate reports a range. A test answers a question with two allowed answers. The machinery is the same - a statistic, its sampling distribution, a tail probability - but the logic is deliberately asymmetric, and reading a test as though it were symmetric is the commonest misuse of statistics in the physical sciences.

The asymmetry is this. One hypothesis, the null, is specific enough to pin down the distribution of the data, so the behaviour of any statistic under the null can be computed. The test then asks a single question: is what we saw unremarkable under the null? If it is not, the null is rejected. If it is, nothing whatever has been proved - the data are merely consistent with the null, as they may also be with a hundred other hypotheses that were never considered. Arfken, Weber and Harris define a statistical hypothesis as any assumption about an unknown probability distribution, and that is all a null hypothesis ever is: an assumption put up to be knocked down.

Definition(Hypotheses, test statistic, rejection region, significance and power)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.7.1, p. 1278, and §31.7.2, pp. 1278-1280

A hypothesis about data $x$ specifies, in whole or in part, the distribution from which $x$ was drawn. It is simple if it determines that distribution uniquely, and composite if it fixes the functional form but leaves one or more parameters free. The hypothesis to be tested is the null hypothesis $H_0$ , and the alternative $H_1$ is its complement within a hypothesis space fixed in advance, so that rejecting one is accepting the other. For a normal sample of known variance, $H_0 : \mu = 0$ is simple and the alternative $H_1 : \mu \neq 0$ is composite; if it were known beforehand that $\mu$ is either $0$ or $1$ , then $H_1 : \mu = 1$ would be simple as well.

Take $H_0$ simple. A test statistic is a scalar function $t(x)$ of the data; since $H_0$ fixes the distribution of $x$ , the sampling distribution $P(t \mid H_0)$ can in principle be computed. A rejection region (or critical region) is a set of values of $t$ carrying probability $\alpha$ under $H_0$ :

\Pr\bigl(t \in \mathcal{R} \mid H_0\bigr) = \alpha .

It is one-tailed if it is of the form $t > t_{\mathrm{crit}}$ or $t < t_{\mathrm{crit}}$ , and two-tailed if it consists of both extremes. If the observed $t$ falls in $\mathcal{R}$ , then $H_0$ is rejected at significance level $\alpha$ ; otherwise it is accepted at that level. The complement of $\mathcal{R}$ is the acceptance region.

Two things can go wrong.

An error of the first kind is rejecting $H_0$ when it is true. Its probability is $\alpha$ , by construction; the significance level is the first-kind error rate.
An error of the second kind is accepting $H_0$ when it is false. Its probability is $\beta = \Pr(t \notin \mathcal{R} \mid H_1)$ , which can only be computed when $H_1$ is simple enough to determine $P(t \mid H_1)$ . The quantity $1 - \beta$ is the power of the test: the probability of rejecting $H_0$ when the particular alternative is true.

Choosing $\alpha$ small makes false alarms rare and misses common. The two rates trade off against each other at fixed sample size, and only more data improves both.

Definition(The

p

-value)Original

None of our sources uses the term, so we define it here in terms of theirs. Suppose the rejection region is the upper tail $t > t_{\mathrm{crit}}$ , with $t_{\mathrm{crit}}$ a decreasing function of $\alpha$ . Let $t_{\mathrm{obs}}$ be the observed value. The $p$ -value is

p \;=\; \Pr\bigl(t \ge t_{\mathrm{obs}} \mid H_0\bigr) ,

the probability, computed under the null, of a result at least as extreme as the one obtained. For a two-tailed region symmetric about zero the definition reads $p = \Pr(|t| \ge |t_{\mathrm{obs}}| \mid H_0)$ .

Its meaning is a matter of bookkeeping: $H_0$ is rejected at level $\alpha$ precisely when $p \le \alpha$ , so $p$ is the smallest significance level at which the data reject the null. Reporting $p$ therefore reports the outcome of the test at every level at once, which is why it is preferred to the bare verdict "rejected at 5%". It is emphatically not the probability that $H_0$ is true: that number does not exist until a prior is supplied, as in Bayesian inference.

Method(The Neyman-Pearson lemma)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.7.3, eqs. (31.107)-(31.108), p. 1280

Let $H_0$ and $H_1$ both be simple, with densities $P_0$ and $P_1$ for the data $x$ . Fix a constant $c > 0$ and let

\mathcal{R}_c \;=\; \bigl\{\, x : P_0(x) < c\,P_1(x) \,\bigr\} ,

and suppose $c$ is chosen so that $\Pr(x \in \mathcal{R}_c \mid H_0) = \alpha$ and the boundary set $\{P_0 = cP_1\}$ has probability zero under both hypotheses. Then $\mathcal{R}_c$ has the greatest power among all rejection regions of size at most $\alpha$ : for every region $\mathcal{R}'$ with $\Pr(\mathcal{R}' \mid H_0) \le \alpha$ ,

\Pr\bigl(\mathcal{R}_c \mid H_1\bigr) \;\ge\; \Pr\bigl(\mathcal{R}' \mid H_1\bigr) .

Equivalently, the best scalar test statistic is the likelihood ratio

\Lambda(x) \;=\; \frac{P_0(x)}{P_1(x)} ,

with rejection region $\Lambda < \Lambda_{\mathrm{crit}}$ . Any monotonically increasing function $u = f(\Lambda)$ serves equally well, with region $u < f(\Lambda_{\mathrm{crit}})$ ; any monotonically decreasing $v = g(\Lambda)$ serves with the region reversed, $v > g(\Lambda_{\mathrm{crit}})$ . This is the form in which the lemma is used, because a monotone rewriting usually turns $\Lambda$ into a statistic whose distribution is tabulated.

Riley, Hobson and Bence state the lemma without proof and give the best statistic as $\Lambda$ with rejection region $\Lambda < \Lambda_{\mathrm{crit}}$ , which is the version above and the one used in all their worked examples. Their accompanying inequality (31.107) is printed with the ratio $P(t|H_0)/P(t|H_1)$ greater than a constant in the rejection region; we follow their (31.108), and the difference is recorded in the handbook's list of source discrepancies.

Method(The generalised likelihood-ratio test)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.7.4, eqs. (31.109)-(31.110) and the large-sample statement on p. 1284, pp. 1282-1284; Riley, Hobson & Bence (2006), §31.7.5, eqs. (31.112)-(31.114), pp. 1284-1285

Let the data be drawn from a population $P(x \mid a)$ of known functional form depending on parameters $a = (a_1, \dots, a_M)$ , with likelihood $L(x; a)$ as in Maximum likelihood. Let $H_0$ assert that $a$ lies in a subspace $S$ of the full parameter space $A$ , and $H_1$ that it does not. Write $\hat a$ for the maximum-likelihood estimate over $A$ and $\hat a_S$ for the maximiser over $S$ alone. The generalised likelihood ratio is

\lambda(x) \;=\; \frac{L(x; \hat a_S)}{L(x; \hat a)} .

$\lambda$ depends only on the data, and $0 \le \lambda \le 1$ , since $S \subseteq A$ . The rejection region is $\lambda < \lambda_{\mathrm{crit}}$ : a value near $1$ says the restricted model fits nearly as well as the free one, a small value says it does not. When $H_1$ is also simple, $\lambda$ reduces to the Neyman-Pearson ratio.
Large-sample law (Wilks). If $H_0$ fixes $R$ of the $M$ parameters at given values, $R \le M$ , then for large sample size $N$ the quantity $-2\ln\lambda$ follows approximately a chi-squared distribution of order $R$ . Riley, Hobson and Bence state this without proof.
Student's $t$ -test. For a sample of size $N$ from $N(\mu, \sigma^{2})$ with both parameters unknown, and $H_0 : \mu = \mu_0$ against $H_1 : \mu \neq \mu_0$ with $\sigma^{2}$ free under both,

\lambda \;=\; \left(1 + \frac{t^{2}}{N-1}\right)^{-N/2} , \qquad t \;=\; \frac{\overline x - \mu_0}{S/\sqrt N} , \qquad S^{2} = \frac{1}{N-1}\sum_{i}(x_i - \overline x)^{2} ,

and $t \sim t(N-1)$ under $H_0$ . Since $\lambda$ is a decreasing function of $t^{2}$ , the rejection region $\lambda < \lambda_{\mathrm{crit}}$ is the two-tailed region $|t| > t_{\mathrm{crit}}$ , with $t_{\mathrm{crit}} = t_{N-1,\,1-\alpha/2}$ for significance $\alpha$ .

Idea of the proofOriginal

The Neyman-Pearson lemma is a budgeting argument. A rejection region is a set of outcomes bought with a budget $\alpha$ of null-hypothesis probability; the aim is to buy as much alternative-hypothesis probability as the budget allows. An outcome $x$ costs $P_0(x)$ and pays $P_1(x)$ , so the outcomes worth buying are those with the best pay-to-cost ratio - and buying everything down to a fixed exchange rate $c$ , which is exactly the region $P_0 < cP_1$ , is optimal. The proof below is this argument with the regions compared term by term.

The generalised ratio extends the idea to the usual situation in which the alternative is a whole family: compare the best fit the null allows with the best fit possible, and reject when restricting the model costs too much likelihood. Student's $t$ -test is nothing but this comparison for a normal sample, rearranged until the statistic is one whose distribution is in a table.

Proof(The lemma, and the

t

-test as a special case)Original

1. The Neyman-Pearson lemma. Let $\mathcal{R} = \mathcal{R}_c$ and let $\mathcal{R}'$ be any region with $\Pr(\mathcal{R}' \mid H_0) \le \alpha$ . Write $\Pr_j(\cdot)$ for probability computed with the density $P_j$ . Splitting each region into its common and exclusive parts, $\mathcal{R} = (\mathcal{R}\cap\mathcal{R}') \cup (\mathcal{R}\setminus\mathcal{R}')$ and likewise for $\mathcal{R}'$ , the shared piece cancels from every difference, so

\Pr\nolimits_1(\mathcal{R}) - \Pr\nolimits_1(\mathcal{R}') = \Pr\nolimits_1(\mathcal{R}\setminus\mathcal{R}') - \Pr\nolimits_1(\mathcal{R}'\setminus\mathcal{R}) .

On $\mathcal{R}\setminus\mathcal{R}'$ we are inside $\mathcal{R}$ , where the defining inequality $P_0 < cP_1$ holds, so $P_1 > P_0/c$ pointwise and therefore

\Pr\nolimits_1(\mathcal{R}\setminus\mathcal{R}') \;\ge\; \frac1c \Pr\nolimits_0(\mathcal{R}\setminus\mathcal{R}') .

On $\mathcal{R}'\setminus\mathcal{R}$ we are outside $\mathcal{R}$ , where $P_0 \ge cP_1$ , so $P_1 \le P_0/c$ and

\Pr\nolimits_1(\mathcal{R}'\setminus\mathcal{R}) \;\le\; \frac1c \Pr\nolimits_0(\mathcal{R}'\setminus\mathcal{R}) .

Subtracting the second from the first and reversing the cancellation of the shared piece,

\Pr\nolimits_1(\mathcal{R}) - \Pr\nolimits_1(\mathcal{R}') \;\ge\; \frac1c\Bigl[\Pr\nolimits_0(\mathcal{R}) - \Pr\nolimits_0(\mathcal{R}')\Bigr] \;=\; \frac{1}{c}\bigl[\alpha - \Pr\nolimits_0(\mathcal{R}')\bigr] \;\ge\; 0 ,

the last step because $c > 0$ and $\Pr_0(\mathcal{R}') \le \alpha$ by hypothesis. That is the claim. The reformulation in terms of $\Lambda$ is immediate: $\mathcal{R}_c = \{\Lambda < c\}$ . A monotonically increasing $f$ preserves the ordering of $\Lambda$ and so maps this set to $\{f(\Lambda) < f(c)\}$ , which is the same set of outcomes; a decreasing $g$ reverses it.

2. The $t$ -test is a generalised likelihood-ratio test. For a normal sample the likelihood is

L(x; \mu, \sigma^{2}) = (2\pi\sigma^{2})^{-N/2}\exp\left[-\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}(x_i-\mu)^{2}\right] .

Over the whole parameter space the maximisers are $\hat\mu = \overline x$ and $\hat\sigma^{2} = N^{-1}\sum_i (x_i-\overline x)^{2}$ , as computed in Maximum likelihood. Over the subspace $\mu = \mu_0$ only $\sigma^{2}$ is free, and the same calculation with $\overline x$ replaced by the fixed $\mu_0$ gives $\hat\sigma_S^{2} = N^{-1}\sum_i (x_i - \mu_0)^{2}$ . At either maximiser the exponent evaluates to $-N/2$ , because the fitted variance is exactly the mean squared residual; so the exponentials cancel in the ratio and

\lambda = \left(\frac{\hat\sigma_S^{2}}{\hat\sigma^{2}}\right)^{-N/2} = \left[\frac{\sum_i (x_i-\mu_0)^{2}}{\sum_i (x_i-\overline x)^{2}}\right]^{-N/2} .

Now expand about the sample mean: writing $x_i - \mu_0 = (x_i - \overline x) + (\overline x - \mu_0)$ and squaring, the cross term is $2(\overline x - \mu_0)\sum_i (x_i-\overline x) = 0$ , so

\sum_i (x_i-\mu_0)^{2} = \sum_i (x_i-\overline x)^{2} + N(\overline x - \mu_0)^{2} .

Dividing through by $\sum_i (x_i - \overline x)^{2} = (N-1)S^{2}$ turns the bracket into $1 + N(\overline x-\mu_0)^{2}/[(N-1)S^{2}]$ , and the second term is $t^{2}/(N-1)$ with $t$ as defined. Hence $\lambda = [1+t^{2}/(N-1)]^{-N/2}$ . The map $t^{2} \mapsto \lambda$ is strictly decreasing, so $\lambda < \lambda_{\mathrm{crit}}$ is the same event as $t^{2} > c$ , that is $|t| > t_{\mathrm{crit}}$ . Finally, under $H_0$ the statistic $t$ is exactly the pivot of Confidence intervals, distributed as $t(N-1)$ by Sampling distributions. $\square$

Remark(Where the hypotheses are used, and what fails without them)Original

Both hypotheses simple, in the lemma. The proof compares two fixed densities pointwise. If $H_1$ is composite, $P_1$ is not defined, the power is a function of the alternative rather than a number, and no region is best against all alternatives at once. (A test that is best against every alternative in a family does sometimes exist - for a one-sided normal mean, for instance - but that is a theorem about the family, not a consequence of the lemma.)

The boundary has probability zero. If $\Pr(\Lambda = c) > 0$ - which happens for discrete data, where the attainable sizes form a finite set - no region has size exactly $\alpha$ , and the optimal test must randomise on the boundary. For continuous data the assumption is automatic.

$H_0$ simple, in the definition of the test. Everything rests on being able to compute $P(t \mid H_0)$ . A composite null leaves that distribution undetermined; the usual repairs are to maximise the $p$ -value over the null (conservative) or to eliminate the nuisance parameter, which is precisely what the $t$ -test does with $\sigma$ .

Normality and independence, in the $t$ -test. Normality enters twice: in writing the likelihood, and in the claim $t \sim t(N-1)$ , which needs the independence of $\overline X$ and $S^{2}$ special to the normal law. Independence of the observations enters in factorising the likelihood. For a large sample the central limit theorem makes the test approximately valid for any parent of finite variance; for a small sample from a heavy-tailed parent it is not valid at all.

Large $N$ , in the Wilks approximation. The chi-squared law for $-2\ln\lambda$ is asymptotic. The worked example below shows how rough it can be: at $N = 8$ it returns $p = 0.0101$ where the exact $t$ calculation gives $0.0199$ , too small by a factor of two, and too small is the dangerous direction.

Linearity in the parameters, in the goodness-of-fit test. The count of degrees of freedom as $N-M$ comes from the fitted parameters imposing $M$ linear constraints on the standardised residuals, which is Riley, Hobson and Bence's constraint-counting rule applied to the normal equations. For a model non-linear in its parameters the constraints are non-linear and the chi-squared law is only approximate; those authors note that the test is used anyway when the non-linearity is mild.

The region is fixed before the data are seen. Every probability above is computed over hypothetical repetitions with the rejection region held fixed. Choosing the region, the tail, or the number of data points after looking at the data destroys the meaning of $\alpha$ , and nothing in the mathematics detects that this has happened.

Remark(Reading the verdict honestly)Original

"Accept $H_0$ " means "no evidence against $H_0$ ". A test has power against some alternatives and almost none against others; failing to reject a constant-rate hypothesis with four degrees of freedom says little about a 2% drift. A non-rejection should always be accompanied by the power of the test against an alternative the experiment actually cares about, or by the confidence interval, which displays the whole set of values not rejected.

A $p$ -value is not a posterior probability. $p = 0.02$ does not mean $H_0$ has probability $0.02$ ; it means that data as extreme as these arise two times in a hundred when $H_0$ holds. Small $p$ -values are produced not only by a false null but also by an incorrect error model, by correlated data treated as independent, and by stopping the experiment when the answer looked good.

Many tests, many false alarms. Twenty independent tests at the 5% level produce on average one rejection when every null is true. A search over a hundred bins, a hundred frequencies or a hundred candidate predictors must correct the threshold accordingly; the particle physicists' five-sigma convention, $p \approx 6 \times 10^{-7}$ , is in large part such a correction.

Statistical significance is not size. With enough data, any effect at all - including an instrumental offset of no physical interest - becomes significant. Report the estimate and its interval alongside the test, and say whether the size of the effect matters.

Method(Running a test)Original

State $H_0$ precisely enough to determine the distribution of the data, and say what the alternative is. Decide one-tailed or two-tailed from the alternative, not from the data.
Choose the statistic. If both hypotheses are simple, the likelihood ratio is optimal by the lemma; otherwise use the generalised ratio $\lambda$ , and rewrite it monotonically into a tabulated statistic if one exists.
Obtain $P(t \mid H_0)$ : exactly where possible (normal, $t$ , $\chi^{2}$ , $F$ from Sampling distributions), by the Wilks approximation for large samples, or by simulating the null.
Fix the significance level $\alpha$ and hence the rejection region, before computing the statistic.
Compute the statistic from the data and report the $p$ -value as well as the verdict.
If the null is not rejected, state the power against an alternative of interest, or quote the confidence interval instead.

Method(The chi-squared goodness-of-fit test)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.7.7, the construction of the test and eq. (31.127), pp. 1296-1298; Riley, Hobson & Bence (2006), §30.9.4, eq. (30.124) and the constraint-counting rule, p. 1193; Arfken, Weber & Harris (2013), §23.7, eq. (23.99) and the moments (23.103)-(23.104), pp. 1170-1172

This tests $H_0 :$ "the model is correct" against the vague alternative that it is not, so no likelihood ratio is available; instead the residuals themselves are the statistic.

Fit. With measurements $y_i$ at known $x_i$ , independent and normally distributed about the model value $f(x_i; a)$ with known standard deviations $\sigma_i$ , estimate the $M$ parameters by minimising

\chi^{2}(a) \;=\; \sum_{i=1}^{N}\left[\frac{y_i - f(x_i; a)}{\sigma_i}\right]^{2} ,

as in Least-squares fitting. Arfken, Weber and Harris emphasise that the weights $1/\sigma_i^{2}$ are what make each term a squared standard normal variable, and so make the sum a chi-squared variable. 2. Count the degrees of freedom. If $f$ is linear in $a$ , the $M$ normal equations are $M$ linear constraints on the $N$ standardised residuals, so by the constraint-counting rule the minimum value $\chi^{2}(\hat a)$ is distributed as $\chi^{2}(N-M)$ . Hence $\mathbb{E}[\chi^{2}(\hat a)] = N-M$ and $\operatorname{var}[\chi^{2}(\hat a)] = 2(N-M)$ . 3. Compare. A value in the range $(N-M) \pm \sqrt{2(N-M)}$ is unremarkable. Equivalently the reduced chi-squared $\chi^{2}/(N-M)$ should be near $1$ . 4. Test. For significance $\alpha$ take the one-tailed region $\chi^{2}(\hat a) > k$ , with $k$ the $1-\alpha$ quantile of $\chi^{2}(N-M)$ , tabulated in table 31.2 of Riley, Hobson and Bence and table 23.2 of Arfken, Weber and Harris. A large value says the model is wrong, or the quoted $\sigma_i$ are too small. A conspicuously small value is also a signal: either too many adjustable parameters are soaking up the scatter, or the error bars have been overstated. 5. Counts in bins. When the data are counts $n_i$ in $k$ bins with expected counts $e_i$ predicted by the model, the Poisson law of Random variables and distributions makes the variance of each count equal to its mean, so $\sigma_i^{2} = e_i$ and the statistic becomes

\chi^{2} \;=\; \sum_{i=1}^{k}\frac{(n_i - e_i)^{2}}{e_i} .

Each count must be large enough for the normal approximation to the Poisson distribution to hold - a working rule is $e_i \gtrsim 5$ - and the degrees of freedom are $k - 1 - m$ , one subtracted because the predicted total is matched to the observed total and $m$ more for parameters estimated from the same data. This binned form is our own assembly of the two cited ingredients, not a formula taken from either book.

Example(Eight dark-current readings: a

t

-test, its

p

-value and its power)Original

A photodiode's dark current is measured on eight separate occasions, giving (in pA)

12.9,\; 11.8,\; 13.4,\; 12.2,\; 13.1,\; 12.5,\; 12.0,\; 12.9 .

The manufacturer specifies $12.0$ pA. Test $H_0 : \mu = 12.0$ against $H_1 : \mu \neq 12.0$ at the 5% level, treating the readings as independent and normal with unknown variance.

Step 1: the statistic. The eight values sum to $100.8$ , so $\overline x = 12.600$ . The deviations are

0.3,\; -0.8,\; 0.8,\; -0.4,\; 0.5,\; -0.1,\; -0.6,\; 0.3 ,

summing to zero, with squares summing to

\sum_i (x_i-\overline x)^{2} = 0.09+0.64+0.64+0.16+0.25+0.01+0.36+0.09 = 2.24 .

Hence $S^{2} = 2.24/7 = 0.32$ , $S = 0.565685$ , and the standard error is $S/\sqrt 8 = \sqrt{0.32/8} = \sqrt{0.04} = 0.200$ exactly. Therefore

t = \frac{12.600 - 12.000}{0.200} = 3.000 .

Step 2: the verdict. With $N - 1 = 7$ degrees of freedom and $\alpha = 0.05$ , the two-tailed critical value is $t_{7,0.975} = 2.3646$ . Since $|t| = 3.000 > 2.3646$ , reject $H_0$ at the 5% level. At the 1% level the critical value is $t_{7,0.995} = 3.4995$ , and $H_0$ would not be rejected; the $p$ -value therefore lies between $0.01$ and $0.05$ , and numerically

p = 2\Pr\bigl(t(7) \ge 3.000\bigr) = 0.0199 .

Step 3: the same test as a likelihood ratio. By item 3 of the statement above,

\lambda = \left(1 + \frac{9}{7}\right)^{-4} = \left(\frac{16}{7}\right)^{-4} = \frac{2401}{65536} = 0.036636 ,

so $-2\ln\lambda = 8\ln(16/7) = 6.6134$ . The null fixes $R = 1$ of the two parameters, so the Wilks approximation compares this with $\chi^{2}(1)$ , whose upper tail at $6.6134$ is $0.0101$ . The exact answer is $0.0199$ : at $N = 8$ the asymptotic law overstates the evidence by a factor of two, and the exact $t$ calculation is the one to use.

Step 4: the power. Suppose the quantity of interest is a dark current of $12.5$ pA or more, and treat $\sigma = 0.5657$ as known so that the normal one-tailed test applies. The 5% rejection region is $\overline x > 12.0 + 1.64485 \times 0.200 = 12.3290$ . Under $\mu = 12.5$ the sample mean is $N(12.5, 0.200^{2})$ , so

1 - \beta = \Pr\left(Z > \frac{12.3290 - 12.5}{0.200}\right) = \Pr(Z > -0.8551) = \Phi(0.8551) = 0.804 .

The test misses a true $12.5$ pA diode about one time in five. That is the number to quote when a test fails to reject.

Example(Five count bins: a goodness-of-fit test)Original

A scintillator counts cosmic-ray muons in five consecutive ten-minute intervals, recording

91,\; 72,\; 85,\; 68,\; 84 ,

a total of $400$ counts. Test the hypothesis that the rate is constant over the fifty minutes, at the 5% level.

Step 1: the expected counts. The rate itself is not predicted, only its constancy, so the expectation in each bin is the observed total shared equally: $e_i = 400/5 = 80$ for each $i$ . Each $e_i$ is comfortably above $5$ , so the binned form of the statistic applies.

Step 2: the statistic. The residuals are $11, -8, 5, -12, 4$ , which sum to zero because the total was matched. Their squares are $121, 64, 25, 144, 16$ , summing to $370$ , so

\chi^{2} = \frac{370}{80} = 4.625 .

Step 3: the degrees of freedom and the verdict. There are $k = 5$ bins, one constraint from matching the total and no further fitted parameters, so $k - 1 = 4$ degrees of freedom. The reduced chi-squared is $4.625/4 = 1.156$ , and the expected range is $4 \pm \sqrt 8 = 4 \pm 2.83$ , which contains $4.625$ comfortably. The 5% critical value is $\chi^{2}_{4,0.95} = 9.488$ ; since $4.625 < 9.488$ , the hypothesis of a constant rate is not rejected. The $p$ -value is

p = \Pr\bigl(\chi^{2}(4) \ge 4.625\bigr) = 0.328 .

Step 4: what has and has not been shown. The data are consistent with a constant rate. They are also consistent with a slow drift: the largest count is in the first bin and the smallest in the fourth, and a dedicated test against a linear trend - one parameter, one degree of freedom - would have more power against that specific alternative than this omnibus test, which spreads its four degrees of freedom over every possible departure from flatness.

Application(Where tests are used, and where they are abused)Original

Discovery claims. A new particle, an exoplanet transit or a gravitational-wave event is announced as the rejection of a background-only null hypothesis. The test statistic is almost always a likelihood ratio, justified by the lemma above, and the threshold is set far beyond the conventional 5% precisely because the search ranges over many masses, periods or templates.

Instrument and model validation. The reduced chi-squared of a calibration fit is the routine health check of an experiment: near $1$ , the model and the error bars are both credible; much above $1$ , one of them is wrong; much below $1$ , the error bars have been inflated or the model over-parametrised. Arfken, Weber and Harris read their three-point fit exactly this way, reporting $\chi^{2} = 4.533$ against an expectation of $2$ and a tail probability just above $0.1$ .

Comparing two samples. Asking whether two materials, two detectors or two model runs differ in mean is a $t$ -test; whether they differ in spread is the $F$ -test built on the variance ratio of Sampling distributions. In climate work the same machinery, applied to a detection-and-attribution fingerprint, is how a forced signal is distinguished from internal variability.

Model selection in fitting. Adding a parameter to a fit always lowers $\chi^{2}$ ; the question is whether it lowers it by more than the one degree of freedom it costs. The generalised likelihood-ratio test with $R = 1$ , compared against $\chi^{2}(1)$ , is the standard answer, and the usual reason a spectral line or an extra exponential component is or is not claimed.

Method

Bayesian inference

Adapted

Add to PDF

Turns Bayes' theorem into a method of parameter estimation - prior, likelihood, evidence, posterior - proves that the posterior updates sequentially and that the beta and normal families are conjugate, defines the credible interval and sets it against the confidence interval, and works a detector efficiency to a posterior mean of $9/14 = 0.643$ with a 95% credible interval $[0.386, 0.861]$ .

Prerequisites: Probability spaces, Maximum likelihood, Confidence intervals, Estimators and their properties, Random variables and distributions, The Beta function

Note(The change of question)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.5.5, the remark that the posterior taken as a whole is what a parameter-estimation problem has to deliver, p. 1265

Frequentist inference keeps the parameter fixed and lets the data vary: it asks how a procedure behaves over hypothetical repetitions of the experiment, and reports intervals that cover the truth at a stated rate. Bayesian inference keeps the data fixed - they are what was actually observed - and puts the distribution on the parameter instead. It asks what is now believed about the parameter, given these data and what was believed beforehand.

The price is that "what was believed beforehand" has to be written down as a probability distribution. The return is that every question is answered by one object, the posterior distribution, and answered in the form people instinctively want: a probability that the parameter lies in a range. Riley, Hobson and Bence make the point sharply: on the Bayesian view what a parameter-estimation problem has to deliver is the whole posterior function, and every number one goes on to quote - an estimate, an error bar, an interval - is only a summary of it.

Definition(Prior, likelihood, evidence, posterior)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.5.5, eqs. (31.81)-(31.83) with the naming of the four factors, pp. 1264-1265; Grimmett & Stirzaker (2020), §6.14, Example (1), p. 330

Fix a model $H$ : a family of distributions for the data $x$ indexed by parameters $a = (a_1,\dots,a_m)$ . The sampling density $P(x \mid a, H)$ , read as a function of $a$ with $x$ held at its observed value, is the likelihood $L(x; a)$ of Maximum likelihood. Bayes' theorem, proved for events in Probability spaces and applied here to densities, gives

P(a \mid x, H) \;=\; \frac{P(x \mid a, H)\,P(a \mid H)}{P(x \mid H)} .

The four factors have names:

$P(a \mid H)$ is the prior, the distribution of the parameters before these data were analysed;
$P(x \mid a, H) = L(x;a)$ is the likelihood, through which the data act;
$P(x \mid H) = \int L(x;a)\,P(a\mid H)\,d^{m}a$ is the evidence (or marginal likelihood), which normalises the result;
$P(a \mid x, H)$ is the posterior, the distribution of the parameters after the data.

Because the evidence does not depend on $a$ , it is usually left out and restored at the end by normalisation:

P(a \mid x, H) \;\propto\; L(x; a)\,P(a \mid H) , \qquad \int P(a \mid x, H)\,d^{m}a = 1 .

For a parameter taking values in a countable set $\Theta$ the integral is a sum, and Grimmett and Stirzaker write the posterior mass function with its denominator displayed, $\pi(\theta \mid x) = f(x\mid\theta)\pi(\theta)\big/\sum_{\psi\in\Theta} f(x\mid\psi)\pi(\psi)$ ; what is usually wanted from it is a posterior expectation $\mathbb{E}(g(\theta)\mid x) = \sum_\theta g(\theta)\pi(\theta\mid x)$ .

A prior should encode whatever was known beforehand, including hard physical restrictions - a mass is positive, an efficiency lies in $[0,1]$ . When little is known, Riley, Hobson and Bence indicate it with a uniform or very broad prior; then $P(a\mid H)$ is constant and the posterior is proportional to the likelihood alone.

Definition(Credible interval)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.5.5, eq. (31.84) with the central and equal-ordinate choices, p. 1266

Let $a$ be a single parameter with posterior density $\pi(a \mid x)$ . A credible interval at level $1 - \alpha - \beta$ is an interval $[a_-, a_+]$ with

\int_{a_-}^{a_+} \pi(a \mid x)\,da \;=\; 1 - \alpha - \beta ,

that is, an interval the parameter lies in with the stated posterior probability. Two standard choices pin down the endpoints:

the central interval, with equal tail probabilities $\alpha = \beta$ ;
the interval with equal posterior density at the endpoints, $\pi(a_-\mid x) = \pi(a_+ \mid x)$ , which for a unimodal posterior is the shortest interval of its level and consists of the most probable values.

Riley, Hobson and Bence give both, under the name Bayesian confidence interval; we use the now-standard term credible interval to keep it distinct from the confidence interval of Confidence intervals, which is a different object with a different meaning. In several dimensions the analogue is a region bounded by a surface of constant posterior density.

Definition(Conjugate family)Original

A family $\mathcal{F}$ of prior distributions is conjugate to a likelihood $L(x;\theta)$ if, whenever the prior belongs to $\mathcal{F}$ , the posterior does too: updating moves the parameters of the prior without leaving the family. The practical content is that the integral defining the evidence is done once, in general, rather than numerically for every data set. The term is not used by any of our sources; the two instances below are proved from scratch.

Method(Properties of the posterior)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.5.5, eqs. (31.82)-(31.83), (31.85) and (31.86), pp. 1265-1268

Let $\pi(a)$ be a prior and $L(x;a)$ a likelihood, with posterior $\pi(a\mid x) \propto L(x;a)\pi(a)$ .

Sequential updating. If data arrive in two batches $x_1$ and $x_2$ that are independent given $a$ , then

\pi(a \mid x_1, x_2) \;\propto\; L(x_2; a)\,\pi(a\mid x_1) :

the posterior after the first batch, used as the prior for the second, gives the same answer as analysing both batches at once. The order of the batches is irrelevant.

Flat prior. If $\pi$ is constant over the region where $L$ is appreciable, then $\pi(a\mid x) \propto L(x;a)$ : the posterior is the normalised likelihood, and its maximiser is the maximum-likelihood estimate $\hat a_{\mathrm{ML}}$ . A Bayesian therefore reads the maximum-likelihood estimate as the posterior mode under a uniform prior.
Curvature error. Riley, Hobson and Bence use

\hat\sigma_{\hat a} = \left[-\left.\frac{\partial^{2}\ln L}{\partial a^{2}}\right|_{a = \hat a}\right]^{-1/2}

as an estimate of the standard error. If $L$ is Gaussian in $a$ this is exactly its standard deviation, and the one-sigma limits are then the $68.3\%$ central credible interval; the formula is used as a measure of the error even when $L$ is not Gaussian. For several parameters the same curvature gives $(V^{-1})_{ij} = -\partial^{2}\ln L/\partial a_i \partial a_j$ at $\hat a$ , whose expectation is the Fisher matrix of Estimators and their properties.

Method(Two conjugate families)Original

1. Beta prior, binomial data. Let $\theta \in (0,1)$ be an unknown success probability, let the data be $k$ successes in $n$ independent trials, and let the prior be $\mathrm{Beta}(\alpha_0, \beta_0)$ , with density proportional to $\theta^{\alpha_0-1}(1-\theta)^{\beta_0-1}$ . (The subscripted $\alpha_0, \beta_0 > 0$ are the parameters of the beta distribution; they have nothing to do with the tail probabilities $\alpha, \beta$ of a credible interval.) Then the posterior is

\mathrm{Beta}\bigl(\alpha_0 + k,\; \beta_0 + n - k\bigr) ,

with mean, and hence posterior expectation of $\theta$ ,

\mathbb{E}[\theta \mid k] = \frac{\alpha_0+k}{\alpha_0+\beta_0+n} = w\,\frac{\alpha_0}{\alpha_0+\beta_0} + (1-w)\,\frac{k}{n} , \qquad w = \frac{\alpha_0+\beta_0}{\alpha_0+\beta_0+n} ,

a weighted average of the prior mean and the observed frequency $k/n$ , which is the maximum-likelihood estimate. The prior acts as $\alpha_0$ imagined successes and $\beta_0$ imagined failures, and its weight $w$ falls like $1/n$ .

2. Normal prior, normal data of known variance. Let $x_1,\dots,x_n$ be independent $N(\theta, \sigma^{2})$ with $\sigma^{2}$ known, and let the prior be $\theta \sim N(\mu_0, \tau^{2})$ . Then the posterior is normal,

\theta \mid x \;\sim\; N\!\left(\mu_n, \frac{1}{P_n}\right) , \qquad P_n = \frac{1}{\tau^{2}} + \frac{n}{\sigma^{2}} , \qquad \mu_n = \frac{\tau^{-2}\mu_0 + (n/\sigma^{2})\,\overline x}{P_n} .

In words: precisions add, and the posterior mean is the average of the prior mean and the sample mean weighted by their precisions. As $n \to \infty$ the data term dominates, $\mu_n \to \overline x$ and the posterior width tends to $\sigma/\sqrt n$ , the frequentist standard error.

Idea of the proofOriginal

Everything here follows from one line, $\text{posterior} \propto \text{likelihood} \times \text{prior}$ , plus the discipline of never normalising until the end. Sequential updating is that line applied twice, with the factorisation of the likelihood over independent batches doing all the work. Conjugacy is the observation that for some pairs the product has the same algebraic shape as the prior: a beta density times a binomial likelihood is $\theta^{\text{something}}(1-\theta)^{\text{something}}$ again, and a Gaussian times a Gaussian is a Gaussian, because the exponents are quadratics and quadratics add. In the Gaussian case the addition is cleanest in terms of the coefficient of $\theta^{2}$ , which is the precision; that is why the rule is "precisions add" and not anything about variances.

Proof(Sequential updating and the two conjugacies)Original

1. Sequential updating. Conditional independence of the batches means $L(x_1,x_2; a) = L(x_1;a)L(x_2;a)$ . Hence

\pi(a \mid x_1,x_2) \;\propto\; L(x_1;a)L(x_2;a)\,\pi(a) \;=\; L(x_2;a)\,\bigl[L(x_1;a)\pi(a)\bigr] \;\propto\; L(x_2;a)\,\pi(a\mid x_1) ,

the last step because $\pi(a\mid x_1) \propto L(x_1;a)\pi(a)$ , the omitted factor being a constant in $a$ . Both ends of the chain are probability densities in $a$ , and two densities that are proportional are equal, since each is the other divided by its total mass, which is one. The product $L(x_1;a)L(x_2;a)$ is symmetric in the two batches, so the order of processing cannot matter. $\square$

2. Beta and binomial. The likelihood of $k$ successes in $n$ trials is $\binom{n}{k}\theta^{k}(1-\theta)^{n-k}$ , and the binomial coefficient is a constant in $\theta$ and may be dropped. Multiplying by the prior,

\pi(\theta \mid k) \;\propto\; \theta^{k}(1-\theta)^{n-k}\cdot\theta^{\alpha_0-1}(1-\theta)^{\beta_0-1} \;=\; \theta^{(\alpha_0+k)-1}(1-\theta)^{(\beta_0+n-k)-1} .

This is the kernel of a $\mathrm{Beta}(\alpha_0+k, \beta_0+n-k)$ density. Since a density is determined by its kernel together with the requirement that it integrate to one, and since the beta integral of The Beta function supplies the constant,

\int_0^1 \theta^{p-1}(1-\theta)^{q-1}\,d\theta = B(p,q) = \frac{\Gamma(p)\Gamma(q)}{\Gamma(p+q)} ,

the posterior is exactly that beta distribution - which is also where the hypothesis $\alpha_0, \beta_0 > 0$ is used, since otherwise the integral diverges. Its mean is $p/(p+q)$ with $p = \alpha_0+k$ and $q = \beta_0+n-k$ , as tabulated in Random variables and distributions, giving $(\alpha_0+k)/(\alpha_0+\beta_0+n)$ . The weighted-average form follows by splitting the numerator:

\frac{\alpha_0+k}{\alpha_0+\beta_0+n} = \frac{\alpha_0+\beta_0}{\alpha_0+\beta_0+n}\cdot\frac{\alpha_0}{\alpha_0+\beta_0} + \frac{n}{\alpha_0+\beta_0+n}\cdot\frac{k}{n} . \qquad\square

3. Normal and normal. Keep only the factors depending on $\theta$ . The likelihood is

L \;\propto\; \exp\left[-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}(x_i-\theta)^{2}\right] , \qquad \sum_{i=1}^{n}(x_i-\theta)^{2} = \sum_{i=1}^{n}(x_i-\overline x)^{2} + n(\theta - \overline x)^{2} ,

the identity being the expansion about $\overline x$ whose cross term vanishes. The first sum is free of $\theta$ and joins the constant. Writing $p_0 = 1/\tau^{2}$ and $p_1 = n/\sigma^{2}$ , the posterior is therefore

\pi(\theta\mid x) \;\propto\; \exp\left[-\tfrac12\Bigl(p_0(\theta-\mu_0)^{2} + p_1(\theta-\overline x)^{2}\Bigr)\right] .

The bracket is a quadratic in $\theta$ with leading coefficient $p_0+p_1 = P_n$ . Completing the square,

p_0(\theta-\mu_0)^{2} + p_1(\theta-\overline x)^{2} = P_n\left(\theta - \frac{p_0\mu_0+p_1\overline x}{P_n}\right)^{2} + \frac{p_0p_1}{P_n}(\mu_0-\overline x)^{2} ,

which is checked by expanding both sides: the $\theta^{2}$ coefficients are $P_n$ , the $\theta$ coefficients are $-2(p_0\mu_0+p_1\overline x)$ , and the constants agree because $p_0\mu_0^{2}+p_1\overline x^{2} - (p_0\mu_0+p_1\overline x)^{2}/P_n = p_0p_1(\mu_0-\overline x)^{2}/P_n$ . The second term is free of $\theta$ , so

\pi(\theta\mid x) \;\propto\; \exp\left[-\frac{P_n}{2}(\theta-\mu_n)^{2}\right] , \qquad \mu_n = \frac{p_0\mu_0+p_1\overline x}{P_n} ,

which is the kernel of $N(\mu_n, 1/P_n)$ ; normalisation supplies the constant. $\square$

Remark(Where the hypotheses are used, and what fails without them)Original

Conditional independence, in sequential updating. If the two batches are dependent given $a$ - the same systematic affecting both runs, say - the likelihood does not factorise and feeding the first posterior in as a prior double-counts the shared information, producing a posterior that is too narrow.

A proper posterior. The prior may be improper (a constant on an unbounded range is not a probability density) provided the posterior is normalisable, and that has to be checked. For Riley, Hobson and Bence's exponential model, $L(x;\tau) = \tau^{-N}\exp(-s/\tau)$ with $s = \sum_i x_i$ , a flat prior on $\tau \in (0,\infty)$ gives

\int_0^\infty \tau^{-N}e^{-s/\tau}\,d\tau = s^{1-N}\,\Gamma(N-1) ,

by the substitution $u = s/\tau$ , which is finite only for $N \ge 2$ . With a single observation the posterior does not exist, and any "estimate" read off it is an artefact of truncating the range.

A flat prior is not a neutral prior. Flatness is not preserved by reparametrisation: if $\pi(\theta)$ is constant and $\phi = 1/\theta$ , then $\pi(\phi) = \pi(\theta)|d\theta/d\phi| \propto \phi^{-2}$ , which is far from flat. A uniform prior on a decay time is an informative, decreasing prior on the decay rate. "No prior information" is therefore not a well-defined state, and a paper that uses a flat prior must say flat in what variable.

The model is a hypothesis too. Everything above is conditional on $H$ . The posterior describes the best parameter values within the assumed family and says nothing about whether the family is right; a badly wrong model yields a sharp, confident and useless posterior. Model checking is a separate step, for which the goodness-of-fit test of Hypothesis testing is the frequentist tool and the evidence $P(x\mid H)$ , compared between models, the Bayesian one.

Prior sensitivity. The normal-normal result quantifies how fast the prior stops mattering: its weight is $\tau^{-2}$ against the data's $n/\sigma^{2}$ , so it is negligible once $n \gg \sigma^{2}/\tau^{2}$ . For a parameter the data barely constrain - a nuisance direction in a many-parameter fit - that never happens, and the posterior in that direction is whatever the prior was. Reporting results under two or three different priors is the standard check.

Remark(Credible intervals against confidence intervals)AdaptedAdapted from Riley, Hobson & Bence (2006), §31.5.5, the comparison of the two kinds of interval, pp. 1266-1268, with the worked exponential example on pp. 1263-1264 and 1267

The two constructions produce similar numbers and assert different things.

A confidence interval at level $1-\alpha$ is a random interval; the level is the long-run fraction of repeated samples of the same size whose interval would contain the fixed true value.
A credible interval at level $1-\alpha$ is a fixed interval computed from the one sample obtained; the level is the posterior probability that the parameter lies in it.

Riley, Hobson and Bence make the contrast on a single data set: for ten observed intervals between telephone calls, modelled as exponential with mean $\tau$ , the maximum-likelihood estimate is $\hat\tau = 3.77$ and the curvature error is $\hat\tau/\sqrt{N} = 1.19$ , while the classical $68\%$ central interval computed from the exact sampling distribution of $\hat\tau$ is $[2.86, 5.46]$ and the Bayesian $68\%$ central interval computed from the normalised likelihood is $[3.16, 6.20]$ . Neither matches the symmetric $3.77 \pm 1.19$ , because the relevant distributions are skewed, and the two intervals do not match each other. They record that a frequentist regards the Bayesian interval as an approximation to the classical one, and a Bayesian regards it as the more meaningful of the two.

The disagreement shrinks with sample size. As the data come to dominate the prior, the posterior approaches a normal centred at the maximum-likelihood estimate with variance given by the inverse curvature of $\ln L$ - whose expectation is the inverse Fisher information - and that is also the large-sample sampling distribution of the estimator established in Maximum likelihood. The two intervals then coincide numerically while continuing to mean different things, which is why in large-sample work the choice of school is often a matter of vocabulary. In small samples, in problems with physical boundaries (a mass near zero, an efficiency near one) and in problems with many nuisance parameters, the two can differ a great deal.

Method(A Bayesian analysis)Original

Write the model. Specify $P(x\mid a, H)$ completely; as a function of $a$ at the observed $x$ this is the likelihood.
Choose the prior. Impose hard physical constraints first (positivity, normalisation, a range). Say in which variable any flat prior is flat. If a conjugate family exists and can represent what is actually known, use it; the update is then arithmetic.
Multiply and normalise. Form $L(x;a)\pi(a)$ , recognise it if possible, otherwise normalise numerically. In more than two or three parameters, sample the posterior with Markov chain Monte Carlo rather than evaluating it on a grid.
Summarise. Quote the posterior mean or mode as the estimate, the posterior standard deviation as the error bar, and a credible interval at a stated level - central or highest-density, and say which. For several parameters, integrate out the ones not being reported; marginalising is the Bayesian treatment of nuisance parameters.
Check. Repeat with a different prior to show the conclusion is driven by the data, and test the model itself separately.

Example(A detector efficiency, and a repeated calibration)Original

Part 1: beta-binomial. A new photon counter is tested with ten identical single-photon pulses and registers seven of them. Its efficiency $\theta$ is to be estimated. Previous counters of this design have had efficiencies scattered broadly around one half, which we encode as a $\mathrm{Beta}(2,2)$ prior: a smooth bump on $(0,1)$ with mean $1/2$ , carrying as much weight as two imagined successes and two imagined failures.

By the conjugacy statement with $\alpha_0 = \beta_0 = 2$ , $n = 10$ and $k = 7$ , the posterior is

\mathrm{Beta}(2+7,\; 2+3) = \mathrm{Beta}(9,5) .

Its summaries follow from the beta formulas in Random variables and distributions:

\text{mean} = \frac{9}{14} = 0.6429 , \qquad \text{mode} = \frac{9-1}{9+5-2} = \frac{8}{12} = 0.6667 ,

\text{s.d.} = \sqrt{\frac{9\cdot 5}{14^{2}\cdot 15}} = \sqrt{\frac{45}{2940}} = 0.1237 .

The maximum-likelihood estimate is $k/n = 0.700$ ; the posterior mean is pulled towards the prior mean $0.5$ , by the weight $w = 4/14 = 0.286$ of the statement: $0.286\times 0.5 + 0.714\times 0.7 = 0.643$ , as computed. The $95\%$ central credible interval, from the quantiles of $\mathrm{Beta}(9,5)$ , is

[\,0.3857,\; 0.8614\,] ,

and the posterior probability that the counter beats fifty per cent efficiency is

\Pr(\theta > \tfrac12 \mid k=7) = \frac{7099}{8192} = 0.8666 ,

a statement of exactly the kind a confidence interval is not allowed to make. Ten pulses are few: the interval is half the unit interval wide. Had the prior been flat, $\mathrm{Beta}(1,1)$ , the posterior would be $\mathrm{Beta}(8,4)$ with mean $2/3 = 0.6667$ and interval $[0.3903,\; 0.8907]$ - a visible but not decisive shift, which is the honest report of how much the prior is doing here.

Part 2: normal-normal. A spectral line position is known from an earlier experiment to be $\theta = 5.00 \pm 0.20$ (arbitrary units), which we take as the prior $N(5.00, 0.20^{2})$ . Nine new readings of known instrumental standard deviation $\sigma = 0.30$ have mean $\overline x = 5.30$ . The precisions are

p_0 = \frac{1}{0.20^{2}} = 25 , \qquad p_1 = \frac{9}{0.30^{2}} = 100 , \qquad P_9 = 125 ,

so the posterior standard deviation is $1/\sqrt{125} = 0.0894$ and the posterior mean is

\mu_9 = \frac{25(5.00) + 100(5.30)}{125} = \frac{655}{125} = 5.240 .

The $95\%$ central credible interval is $5.240 \pm 1.95996 \times 0.0894 = [\,5.065,\; 5.415\,]$ . For comparison, the frequentist $95\%$ confidence interval from the new data alone is $5.300 \pm 1.95996 \times 0.30/3 = [\,5.104,\; 5.496\,]$ . The Bayesian interval is shifted towards the earlier measurement and is narrower, $0.351$ against $0.392$ , because it uses both experiments; the whole of that difference is the information in the prior, and nothing else. Note also that by sequential updating the same answer would have come from analysing the two experiments jointly in either order.

Application(Where posteriors are the working object)Original

Cosmology and astrophysics. Parameters of a cosmological model are quoted as posterior means with credible intervals, obtained by sampling a posterior over a dozen or more parameters; positivity and physical ranges enter as priors, and nuisance parameters for instrumental systematics are marginalised away rather than fitted and forgotten. Comparing the evidence $P(x\mid H)$ between models is how an extra parameter is judged to be warranted.

Sequential estimation and filtering. The Kalman filter of The Kalman filter is the normal-normal update of this unit run once per observation, with a forecast step in between; "precisions add" is its measurement update. Operational weather forecasting does the same on a state space of $10^{8}$ variables, where the posterior is represented by an ensemble - see Optimal interpolation and variational data assimilation.

Computation. Beyond two or three parameters the posterior cannot be tabulated, and the integrals wanted from it are posterior expectations. Grimmett and Stirzaker introduce Markov chain Monte Carlo with exactly this example: build a chain whose stationary distribution is the posterior and average along its path. The convergence theory is that of Convergence of Markov chains, and the sampling machinery that of Monte Carlo methods.

Few counts and hard boundaries. When an experiment records three events where two were expected, or measures a quantity known to be non-negative and obtains a negative point estimate, frequentist intervals become awkward and can even be empty. A prior supported on the physical range removes the difficulty by construction, which is why upper limits on masses, rates and cross-sections are so often reported as credible rather than confidence limits - with the prior stated, since the answer depends on it.