Maximum likelihood

Every estimator so far was found by analogy: to estimate a population mean, take a sample mean. That runs out as soon as the parameter is not a simple feature of the distribution, and it gives no guidance at all when two parameters have to be estimated together. This lesson gives the general recipe, which Ronald Fisher developed between 1912 and 1922 and which underlies almost every estimate produced by statistical software today.

The likelihood function

Start with a model that says how probable the data is, given a parameter. For discrete data, $P (data | θ)$ is an ordinary probability; for continuous data it is a density. Now hold the data fixed at what was actually observed and read the expression as a function of $θ$ . That function is the likelihood,

L (θ) = P (x_{1}, dots, x_{n} | θ)

and for independent observations it factorises into $\prod_{i} f (x_{i} | θ)$ .

The idea is then blunt: prefer the value of $θ$ that makes the observed data as probable as it can be. The maximum likelihood estimator is

\hat{θ} = \arg \max_{θ} L (θ)

Two cautions about reading the likelihood, and they are the same caution twice. $L (θ)$ is not the probability that $θ$ is true; $θ$ is a fixed unknown constant and has no probability distribution in this framework. And $L$ is not a density in $θ$ : it does not integrate to one, and its value at a point has no meaning on its own. Only ratios of likelihoods at different $θ$ carry information, which is why the maximum is what gets used.

Why the logarithm

In practice one maximises $ℓ (θ) = \ln L (θ)$ instead. Since the logarithm is strictly increasing, the maximum is at the same place, and the product over observations becomes a sum:

ℓ (θ) = \sum_{i = 1}^{n} \ln f (x_{i} | θ)

Sums differentiate one term at a time, and they do not underflow: the likelihood of 1000 observations each of probability around 0.3 is about $10^{- 523}$ , which is zero in double precision, while its logarithm is a comfortable $- 1204$ . Every serious implementation works on the log scale for both reasons.

The recipe is then calculus. Differentiate $ℓ$ , set the derivative to zero, solve, and check the second derivative is negative. The derivative $ℓ' (θ)$ is called the score, and the equation $ℓ' (θ) = 0$ is the likelihood equation.

The Bernoulli case

Let each observation be 0 or 1 with $P (X = 1) = p$ , and suppose $k$ of $n$ observations are 1. Then

L (p) = p^{k} (1 - p)^{n - k}, ℓ (p) = k \ln p + (n - k) \ln (1 - p)

Differentiating, $ℓ' (p) = \frac{k}{p} - \frac{n - k}{1 - p}$ . Setting this to zero gives $k (1 - p) = (n - k) p$ , so $k = n p$ and

\hat{p} = \frac{k}{n}

The sample proportion, which everybody would have written down anyway, is the maximum likelihood estimator, and now it is derived rather than assumed. The second derivative is $ℓ'' (p) = - k / p^{2} - (n - k) / (1 - p)^{2}$ , negative everywhere, so this is a maximum.

Example. A coin gives 60 heads in 100 tosses. What is the maximum likelihood estimate of $p$ , and how much more likely does it make the data than $p = 0.5$ ?

$\hat{p} = 0.60$ . The likelihood ratio is $\frac{{0.6}^{60} {0.4}^{40}}{{0.5}^{100}}$ , whose logarithm is $60 \ln 0.6 + 40 \ln 0.4 - 100 \ln 0.5 = 2.013$ , so the ratio is $e^{2.013} = 7.5$ . The observed data is 7.5 times more probable under $p = 0.6$ than under a fair coin, which is a preference and nowhere near a proof.

Now you. A different coin gives 7 heads in 20 tosses. What is $\hat{p}$ , and what is the log likelihood ratio against $p = 0.5$ ?

Answer

$\hat{p} = 0.35$ . The log ratio is $7 \ln 0.35 + 13 \ln 0.65 - 20 \ln 0.5 = 0.914$ , so a ratio of 2.5. Weaker evidence from a smaller sample, as expected.

Poisson, exponential and normal

The same three steps handle every standard model, and the answers are worth having.

For Poisson counts with mean $λ$ , $ℓ (λ) = \sum_{i} (x_{i} \ln λ - λ - \ln x_{i}!)$ , so $ℓ' (λ) = \frac{\sum x_{i}}{λ} - n$ , and $\hat{λ} = \overline{x}$ . For exponential waiting times with rate $λ$ , the density is $λ e^{- λ x}$ , so $ℓ (λ) = n \ln λ - λ \sum x_{i}$ , giving $\hat{λ} = 1 / \overline{x}$ .

For the normal with both parameters unknown,

ℓ (μ, σ^{2}) = - \frac{n}{2} \ln (2 π) - \frac{n}{2} \ln σ^{2} - \frac{1}{2 σ^{2}} \sum_{i} (x_{i} - μ)^{2}

Differentiating with respect to $μ$ gives $\frac{1}{σ^{2}} \sum (x_{i} - μ) = 0$ , so $\hat{μ} = \overline{x}$ whatever $σ^{2}$ is. Substituting that and differentiating with respect to $σ^{2}$ gives

{\hat{σ}}^{2} = \frac{1}{n} \sum_{i} (x_{i} - \overline{x})^{2}

which is the divide-by- $n$ version, not $s^{2}$ . Maximum likelihood is biased here, by exactly the factor $(n - 1) / n$ derived in the previous lesson, so at $n = 10$ it is 10 per cent low. That is a real defect and a small one: the bias vanishes as $n$ grows, and no general recipe can be both automatic and unbiased for every parameter at once.

Example. Ladislaus Bortkiewicz recorded deaths from horse kicks in 14 Prussian cavalry corps, and the classic subset covers 200 corps-years with 122 deaths. Estimate the Poisson rate and give its standard error, using $SE = \sqrt{λ / n}$ .

$\hat{λ} = 122 / 200 = 0.61$ deaths per corps-year, and $SE = \sqrt{0.61 / 200} = 0.055$ . The estimate is about eleven standard errors away from zero and about seven from one, so the data pins the rate down well.

Now you. Fifty service times have mean 4.0 minutes. Estimate the exponential rate $λ$ and its standard error, given that $SE (\hat{λ}) = λ / \sqrt{n}$ .

Answer

$\hat{λ} = 1 / 4.0 = 0.25$ per minute, and $SE = 0.25 / \sqrt{50} = 0.0354$ per minute.

Curvature is precision

The recipe so far returns a point. It also returns, for free, an estimate of how precise that point is, and the mechanism is geometric: a sharply peaked log likelihood means that nearby parameter values explain the data much worse, so the data is informative; a flat one means many values do nearly as well.

Sharpness is curvature, so define the observed information as the negative second derivative of the log likelihood at the maximum, $J = - ℓ'' (\hat{θ})$ . The Fisher information is its expected value, $I (θ) = - E [ℓ'' (θ)]$ . The large-sample result is that

SE (\hat{θ}) \approx \frac{1}{\sqrt{I (θ)}}

For the Bernoulli case, $- ℓ''$ evaluated at $\hat{p} = k / n$ works out to $\frac{n}{\hat{p} (1 - \hat{p})}$ , so

SE (\hat{p}) = \sqrt{\frac{\hat{p} (1 - \hat{p})}{n}}

which is the formula used without justification in the first lesson, now derived. For the Poisson, $ℓ'' = - \sum x_{i} / λ^{2}$ , whose negative expectation is $n / λ$ , giving $SE = \sqrt{λ / n}$ . For the normal mean, the same route returns $σ / \sqrt{n}$ , agreeing with the exact derivation two lessons back.

There is a strong theorem behind this. Under regularity conditions, as $n$ grows the maximum likelihood estimator is consistent, its distribution approaches a normal centred at the true $θ$ with variance $1 / I (θ)$ , and that variance attains the Cramér-Rao lower bound. So maximum likelihood is asymptotically the best possible unbiased procedure, which is why it is the default and why almost every standard error printed by software is a square root of an inverse information.

Example. In a sample of 400 voters, 240 support a proposal. Give the maximum likelihood estimate and its standard error.

$\hat{p} = 240 / 400 = 0.60$ and $SE = \sqrt{0.6 \times 0.4 / 400} = 0.0245$ . Doubling the sample from the 100-toss example halved neither the estimate nor the standard error but divided the latter by $\sqrt{4} = 2$ , from 0.049 to 0.0245.

Now you. With 14 successes in 40 trials, what are $\hat{p}$ and its standard error?

Answer

$\hat{p} = 0.35$ and $SE = \sqrt{0.35 \times 0.65 / 40} = 0.0754$ .

Invariance, and what it costs

A convenient property: maximum likelihood is invariant under reparameterisation. If $\hat{θ}$ maximises the likelihood for $θ$ , then $g (\hat{θ})$ is the maximum likelihood estimate of $g (θ)$ for any function $g$ , because relabelling the axis does not move the peak.

So from $\hat{p} = 0.6$ , the maximum likelihood estimate of the odds $p / (1 - p)$ is $0.6 / 0.4 = 1.5$ , and of the log odds is $\ln 1.5 = 0.405$ , with no fresh optimisation. Nothing comparable holds for unbiasedness, which is destroyed by nonlinear transformation as the previous lesson showed. Invariance is one of the reasons the method has spread: parameters get reparameterised constantly, into odds, rates, logs and proportions, and the estimates travel with them.

Standard errors travel too, by the delta method: if $\hat{θ}$ has standard error $σ_{\hat{θ}}$ , then $g (\hat{θ})$ has approximate standard error $| g' (\hat{θ}) | σ_{\hat{θ}}$ , which is the first term of a Taylor expansion. For the log odds, $g' (p) = 1 / (p (1 - p))$ , so the standard error of $\ln (\hat{p} / (1 - \hat{p}))$ at $\hat{p} = 0.6$ with $n = 100$ is $0.049 / (0.6 \times 0.4) = 0.204$ . Every odds ratio reported with an interval in a medical journal comes from that calculation.

Where the method fails

The asymptotic theory rests on regularity conditions, and it is worth knowing the two ways they break.

The first is a parameter on the boundary of the support. For a uniform distribution on $[0, θ]$ , the likelihood is $θ^{- n}$ for $θ \geq \max x_{i}$ and zero below, so it has no stationary point at all: it is maximised at the smallest permitted value, $\hat{θ} = \max x_{i}$ . Differentiating would have found nothing. This estimator is biased low, since the maximum of a sample can never exceed $θ$ , with $E [\hat{θ}] = n θ / (n + 1)$ , and its sampling distribution is nowhere near normal even for large $n$ . The previous lesson's rescaled version repairs the bias, and the normal approximation never arrives.

The second is small samples generally. Everything above is asymptotic, and "asymptotic" is not a sample size. With $n = 8$ , a standard error taken from the curvature of the log likelihood can be badly wrong, and the estimator's distribution can be visibly skewed. Separation in logistic regression is the everyday version: when one group has no failures at all, the likelihood is maximised at an infinite coefficient, software reports an enormous estimate with an enormous standard error, and the honest conclusion is that the data cannot identify that parameter.

Maximum likelihood now supplies both an estimate and a standard error for essentially any model. A single number with a standard error beside it is still a poor way to communicate a result, because it invites the reader to treat the number as the answer. Turning the pair into an interval with a stated and checkable property is the next lesson.