5.9.1 Likelihood and the maximum likelihood estimator

In this lesson

The likelihood: the probability of the data you saw, read as a function of the unknown parameter
Why a likelihood is not a probability distribution for the parameter
The log-likelihood: a sum instead of a product, with the same peak
The maximum likelihood estimator by calculus, for a coin, for waiting times and for a normal, with the conditions of lesson 4.4.1
The same estimate by computer: minimising the negative log-likelihood with scipy.optimize.minimize
Where the method needs care: a peak at the edge, several peaks, and parameters that must stay positive

You need independent events (lesson 5.1.3), the binomial, exponential and normal distributions (lessons 5.2.2, 5.4.2 and 5.4.3) and the rules of logarithms (lesson 1.4.3). The calculus is the second-derivative test (lesson 2.4.2) and the gradient, Hessian and optimality conditions of lessons 4.2.1, 4.3.1 and 4.4.1. Lesson 5.8.1 treats an estimator as a random variable with a bias.

The likelihood

A coin lands heads with an unknown probability $p$ . You toss it 10 times and see H H T H H H T H T H: 7 heads and 3 tails. For a given $p$ the tosses are independent, so the chance of this exact sequence is a product of 10 factors (lesson 5.1.3): $p^7(1 - p)^3$ .

Before the tosses, $p$ is fixed and the outcome is random. Afterwards the outcome is fixed and $p$ is what you do not know. Read the same formula as a function of $p$ and it is called the likelihood, $L(p)$ .

$p$	0.3	0.5	0.6	0.7	0.8	0.9
$L(p) = p^7(1 - p)^3$	0.000075	0.000977	0.001792	0.002224	0.001678	0.000478

The data are most likely when $p = 0.7$ . They are 2.28 times as likely under $p = 0.7$ as under a fair coin, and 29.6 times as likely as under $p = 0.3$ .

Rule: Likelihood

\begin{gathered} L(\theta) = f(x_1; \theta)\, f(x_2; \theta) \cdots f(x_n; \theta) \\[4pt] \text{data } x_1, \ldots, x_n \text{ fixed, } \theta \text{ varies} \end{gathered}

$\theta$ (theta) stands for the unknown parameter, and $f(x; \theta)$ is the pmf or density of one observation when the parameter is $\theta$ . The product needs independent observations. For a continuous model the factors are densities (lesson 5.4.1), not probabilities.

It is not a distribution for the parameter

A pmf or density must have total probability 1. The likelihood need not: the area under $L(p)$ from 0 to 1 is $\frac{1}{1320}$ . It says how well each value of $p$ explains the data, not how probable $p$ is.

Only comparisons between values of $\theta$ mean something, for the same data. That is why the table quotes ratios.

L(p)

import numpy as np from scipy.integrate import quad def L(p): return p**7 * (1 - p)**3 for p in [0.3, 0.5, 0.6, 0.7, 0.8, 0.9]: print(p, f"{L(p):.6f}", round(L(p) / L(0.7), 4)) print(round(L(0.7) / L(0.5), 4), round(L(0.7) / L(0.3), 1)) print(round(quad(L, 0, 1)[0], 10), round(1 / 1320, 10)) grid = np.linspace(0.01, 0.99, 99) print(round(float(grid[np.argmax(L(grid))]), 2))

⌘/Ctrl + Enter runs the cell

The log-likelihood

A product of many numbers below 1 shrinks fast. With 1,000 tosses and 700 heads, $L(0.7) = 0.7^{700}\,0.3^{300}$ is about $5.07 \times 10^{-266}$ . With 2,000 tosses Python gives 0, because the number is below the smallest a computer stores.

Logs fix this. The log of a product is the sum of the logs (lesson 1.4.3), so the log-likelihood is a sum of moderate numbers. And $\ln$ is increasing, so $L$ and $\ln L$ peak at the same $\theta$ .

Rule: Log-likelihood

\ell(\theta) = \ln L(\theta) = \sum_{i=1}^{n} \ln f(x_i; \theta)

$\ell$ is the script letter l. It has the same maximiser as $L$ . A difference $\ell(\theta_1) - \ell(\theta_2)$ is the log of the ratio $\frac{L(\theta_1)}{L(\theta_2)}$ .

For the coin, $\ell(p) = 7\ln p + 3\ln(1 - p)$ . Then $\ell(0.7) = -6.1086$ and $\ell(0.5) = -6.9315$ . The difference, 0.8228, is $\ln 2.2769$ : the ratio from the table.

The likelihood and the log-likelihood of 7 heads in 10 tosses. Taking logs changes the heights, not the place of the peak: both are highest at $p = 0.7$ .

\ell(0.7)

import math for n in [100, 1000, 2000]: k = int(0.7 * n) print(n, f"{0.7**k * 0.3**(n - k):.3g}", round(k * math.log(0.7) + (n - k) * math.log(0.3), 2)) def ell(p): return 7 * math.log(p) + 3 * math.log(1 - p) print(round(ell(0.7), 4), round(ell(0.5), 4), round(math.exp(ell(0.7) - ell(0.5)), 4))

⌘/Ctrl + Enter runs the cell

The maximum likelihood estimator

The maximum likelihood estimate $\hat{\theta}$ is the value of $\theta$ that makes the observed data most likely: the maximiser of $L$ , or of $\ell$ . It is computed from the data, so before the data arrive it is a random variable, an estimator in the sense of lesson 5.8.1.

A coin, by calculus

Differentiate $\ell(p) = 7\ln p + 3\ln(1 - p)$ , using $\frac{d}{dp}\ln(1 - p) = -\frac{1}{1 - p}$ (the chain rule, lesson 2.3.3):

\ell'(p) = \frac{7}{p} - \frac{3}{1 - p}

Setting $\ell'(p) = 0$ gives $7(1 - p) = 3p$ , so $\hat{p} = 0.7$ . With $k$ heads in $n$ tosses the same steps give $\hat{p} = \frac{k}{n}$ : the share of heads.

The second derivative is $\ell''(p) = -\frac{7}{p^2} - \frac{3}{(1 - p)^2}$ , negative for every $p$ . By lesson 2.4.2's test the critical point is a maximum. At $\hat{p} = 0.7$ it is $-47.62$ , a number lesson 5.9.2 puts to use.

Rule: Maximum likelihood, one parameter

\begin{gathered} \ell'(\hat{\theta}) = 0 \\[4pt] \ell''(\hat{\theta}) < 0 \end{gathered}

The slope $\ell'(\theta)$ is called the score. A zero score is needed at a peak inside the range of $\theta$ ; a negative second derivative makes it a maximum.

Waiting times

Five waits for a bus are 2, 0.5, 1.5, 4 and 2 minutes. Model each as exponential with rate $\lambda$ per minute (lesson 5.4.2), with density $f(x; \lambda) = \lambda e^{-\lambda x}$ . The log of each factor is $\ln\lambda - \lambda x$ , so

\begin{aligned} \ell(\lambda) &= 5\ln\lambda - \lambda\,(2 + 0.5 + 1.5 + 4 + 2) \\ &= 5\ln\lambda - 10\lambda \end{aligned}

Then $\ell'(\lambda) = \frac{5}{\lambda} - 10 = 0$ gives $\hat{\lambda} = 0.5$ per minute, and $\ell''(\lambda) = -\frac{5}{\lambda^2} < 0$ . In general $\hat{\lambda} = \frac{1}{\bar{x}}$ : one over the mean wait, 2 minutes.

$\lambda$	0.25	0.4	0.5	0.6	1
$\ell(\lambda)$	$-9.4315$	$-8.5815$	$-8.4657$	$-8.5541$	$-10$

\frac{7}{0.5} - \frac{3}{0.5} = 8

import math import numpy as np from scipy.stats import expon def ell(p): return 7 * math.log(p) + 3 * math.log(1 - p) h = 1e-6 print(round((ell(0.7 + h) - ell(0.7 - h)) / (2 * h), 6), round((ell(0.5 + h) - ell(0.5 - h)) / (2 * h), 6)) d = 1e-4 print(round((ell(0.7 + d) - 2 * ell(0.7) + ell(0.7 - d)) / d**2, 2), round(-7 / 0.7**2 - 3 / 0.3**2, 2)) waits = np.array([2, 0.5, 1.5, 4, 2]) print(1 / waits.mean()) for lam in [0.25, 0.4, 0.5, 0.6, 1]: print(lam, round(5 * math.log(lam) - 10 * lam, 4), round(expon(scale=1 / lam).logpdf(waits).sum(), 4))

⌘/Ctrl + Enter runs the cell

Two parameters: the normal

Five measurements, 2, 4, 5, 6 and 8, come from $\mathcal{N}(\mu, \sigma^2)$ with both $\mu$ and $\sigma$ unknown. The log of the normal density (lesson 5.4.3) is $-\frac{1}{2}\ln(2\pi) - \ln\sigma - \frac{(x - \mu)^2}{2\sigma^2}$ . Adding over the $n$ points:

\begin{aligned} \ell(\mu, \sigma) ={}& -\tfrac{n}{2}\ln(2\pi) - n\ln\sigma \\ &- \frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i - \mu)^2 \end{aligned}

The conditions of lesson 4.4.1

The highest point of $\ell$ is the lowest point of $-\ell$ , so lesson 4.4.1's conditions apply to $-\ell$ . First order: the gradient of $\ell$ is zero. Second order: the Hessian of $-\ell$ is positive definite, so the Hessian of $\ell$ is negative definite.

The two partial derivatives (lesson 4.2.1) are

\begin{gathered} \frac{\partial \ell}{\partial \mu} = \frac{1}{\sigma^2}\sum_{i}(x_i - \mu) \\[6pt] \frac{\partial \ell}{\partial \sigma} = -\frac{n}{\sigma} + \frac{1}{\sigma^3}\sum_{i}(x_i - \mu)^2 \end{gathered}

The first is zero when $\mu$ is the mean $\bar{x} = 5$ . The second is then zero when $\sigma^2 = \frac{1}{n}\sum(x_i - \bar{x})^2 = \frac{9 + 1 + 0 + 1 + 9}{5} = 4$ , so $\hat{\sigma} = 2$ .

Rule: Maximum likelihood for a normal

\begin{gathered} \hat{\mu} = \bar{x} \\[4pt] \hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2 \end{gathered}

$\hat{\sigma}^2$ divides by $n$ , not $n - 1$ : here 4, against the sample variance 5. It is the biased estimator of lesson 5.8.1, whose bias shrinks as $n$ grows. Case study 5.3.5's volatility, which divides by $n$ , is this estimate.

Differentiating again gives the Hessian of $\ell$ (lesson 4.3.1). At $(\hat{\mu}, \hat{\sigma})$ the cross term $-\frac{2}{\sigma^3}\sum(x_i - \mu)$ is zero, and

\begin{gathered} H = \begin{pmatrix} -n/\hat{\sigma}^2 & 0 \\ 0 & -2n/\hat{\sigma}^2 \end{pmatrix} \\[6pt] = \begin{pmatrix} -1.25 & 0 \\ 0 & -2.5 \end{pmatrix} \end{gathered}

Both eigenvalues are negative, so $-H$ is positive definite and $(5, 2)$ is a strict local maximum. Far away $\ell$ falls without limit: as $\sigma \to 0$ , as $\sigma \to \infty$ and as $\mu$ moves off. So by lesson 4.4.1 this peak is the global maximum.

A grid of $\ell(\mu, \sigma)$ agrees: every neighbour of $(5, 2)$ is lower.

	$\sigma = 1.5$	$\sigma = 2$	$\sigma = 2.5$
$\mu = 4$	$-12.1776$	$-11.1854$	$-11.1761$
$\mu = 5$	$-11.0665$	$-10.5604$	$-10.7761$
$\mu = 6$	$-12.1776$	$-11.1854$	$-11.1761$

Maximising over $\sigma^2$ instead of $\sigma$ gives the same point: the estimate of $\sigma$ is the square root of the estimate of $\sigma^2$ . In general the maximum likelihood estimate of a function $g(\theta)$ is $g(\hat{\theta})$ .

(5, 2)

import numpy as np from scipy.stats import norm x = np.array([2, 4, 5, 6, 8.0]) print(x.mean(), x.var(), x.std(), x.var(ddof=1)) def ell(m, s): return norm.logpdf(x, m, s).sum() for m in [4, 5, 6]: print(m, [round(float(ell(m, s)), 4) for s in [1.5, 2, 2.5]]) def hessian(f, a, d=1e-4): a = np.array(a, dtype=float) n = len(a) E = np.eye(n) * d H = np.empty((n, n)) for i in range(n): for j in range(n): H[i, j] = (f(a + E[i] + E[j]) - f(a + E[i] - E[j]) - f(a - E[i] + E[j]) + f(a - E[i] - E[j])) / (4 * d * d) return H H = hessian(lambda p: ell(p[0], p[1]), [5, 2]) print(np.round(H, 4) + 0.0, np.round(np.linalg.eigvalsh(H), 4))

⌘/Ctrl + Enter runs the cell

By computer: minimise the negative log-likelihood

Most models have no formula for $\hat{\theta}$ ; case study 5.9.4 fits one, the Student t. Then a computer searches for the peak. Optimisers minimise, so hand them the negative log-likelihood $-\ell(\theta)$ , as lesson 4.4.1 turned a maximum of $f$ into a minimum of $-f$ .

Rule: Numerical maximum likelihood

\hat{\theta} = \text{the minimiser of } -\ell(\theta)

In SciPy: minimize(nll, x0), with nll the negative log-likelihood and x0 a starting guess. The result holds the minimiser res.x, the minimum res.fun $= -\ell(\hat{\theta})$ , and res.success.

A standard deviation must stay positive, but an optimiser does not know that: norm.logpdf with a negative scale returns nan. The cure is to search over $a = \ln\sigma$ , which can be any number, and report $\hat{\sigma} = e^{\hat{a}}$ . By the rule for $g(\hat{\theta})$ above, this is still the maximum likelihood estimate.

(\mu, a)

import math import numpy as np from scipy.optimize import minimize from scipy.stats import norm x = np.array([2, 4, 5, 6, 8.0]) def nll(params, data): m, a = params return -norm.logpdf(data, m, math.exp(a)).sum() res = minimize(nll, [0.0, 0.0], args=(x,)) print(res.success, round(res.x[0], 4), round(math.exp(res.x[1]), 4), round(res.fun, 4)) print(norm.logpdf(5, 5, -2)) y = np.random.default_rng(591).normal(10, 3, size=500) res = minimize(nll, [0.0, 0.0], args=(y,)) print(round(res.x[0], 4), round(math.exp(res.x[1]), 4), round(y.mean(), 4), round(y.std(), 4))

⌘/Ctrl + Enter runs the cell

Where the method needs care

A peak at the edge. Toss a coin 5 times and see no heads. Then $L(p) = (1 - p)^5$ , largest at $p = 0$ , where its slope is not zero.

The score $-\frac{5}{1 - p}$ is never 0, so setting it to zero finds nothing. Check the ends of the range, as lesson 2.4.1 did on a closed interval: $\hat{p} = 0$ , a coin that never lands heads.

Several peaks. Two points, $-3$ and 3, come from a t with $\nu = 1$ (lesson 5.4.4) shifted by an unknown $\theta$ , with density $\frac{1}{\pi(1 + (x - \theta)^2)}$ . Its log-likelihood has two equal peaks, at $\theta = \pm\sqrt{8} \approx \pm 2.8284$ , and a dip at 0.

$\theta$	$-4$	$-2.8284$	$-2$	0	2	2.8284	4
$\ell(\theta)$	$-6.8946$	$-5.873$	$-6.2407$	$-6.8946$	$-6.2407$	$-5.873$	$-6.8946$

An optimiser climbs the hill nearest its start. From 0.5 it finds 2.8284, from $-0.5$ it finds $-2.8284$ , and from 0 it stays at the dip, where the score is also zero. Check the second-order condition and try several starts, as case study 4.4.5 did.

Parameters with limits. A standard deviation, a rate or a probability has a range. Search over a transformed parameter that can take any value, such as $\ln\sigma$ , or tell the optimiser the bounds.

\theta

import math import numpy as np from scipy.optimize import minimize from scipy.stats import t data = np.array([-3.0, 3.0]) def ell(theta): return t.logpdf(data, 1, theta).sum() for theta in [-4, -math.sqrt(8), -2, 0, 2, math.sqrt(8), 4]: print(round(theta, 4), round(ell(theta), 4)) for start in [0.5, -0.5, 0.0]: res = minimize(lambda p: -ell(p[0]), [start]) print(start, round(res.x[0], 4)) d = 1e-4 print(round((ell(d) - 2 * ell(0) + ell(-d)) / d**2, 4))

⌘/Ctrl + Enter runs the cell

Lesson 5.9.2 reads how sharp the peak is, which says how far $\hat{\theta}$ may be from the truth. Case study 5.9.4 fits a Student t to daily returns this way.

Common mistakes

The likelihood is not the probability of $\theta$ . For 7 heads in 10 its area over $p$ is $\frac{1}{1320}$ . Only ratios for the same data mean something.
Multiply only independent factors. The product rule for $L$ needs independent observations.
Use logs in code. The product for 2,000 tosses underflows to 0; the log-likelihood does not.
Hand the optimiser $-\ell$ , not $\ell$ . Minimising $\ell$ itself runs away from the peak.
The normal's $\hat{\sigma}^2$ divides by $n$ . It is biased (lesson 5.8.1), and it is the maximum likelihood estimate.
A zero score is not enough. It misses a peak at the edge, such as $\hat{p} = 0$ , and it also holds at a dip. Check $\ell''$ , or the Hessian, and try several starts.
Keep parameters in range. Search over $\ln\sigma$ , not $\sigma$ .

Likelihood and the maximum likelihood estimator

0 of 12 correct

1.
A coin shows 4 heads in 10 independent tosses. Find the likelihood at $p = 0.5$ (a fraction)
2.
For the same 4 heads in 10 tosses, find the maximum likelihood estimate $\hat{p}$
3.
One head and one tail give $L(p) = p(1 - p)$ . Find the area under it from 0 to 1
4.
For 7 heads in 10, how many times as likely are the data under $p = 0.7$ as under $p = 0.5$ ? (4 significant figures)
5.
$\ell(p) = 7\ln p + 3\ln(1 - p)$ . Find the score $\ell'(p)$ at $p = 0.5$
6.
Waits of 1, 3 and 2 minutes are modelled as exponential with rate $\lambda$ . Find $\hat{\lambda}$ per minute
7.
Four exponential waits add up to 10 minutes, so $\ell(\lambda) = 4\ln\lambda - 10\lambda$ . Find $\ell''(\hat{\lambda})$
8.
Measurements 1, 3, 4, 6, 6 come from a normal. Find the maximum likelihood estimate $\hat{\sigma}^2$
9.
For the same data, find the maximum likelihood estimate $\hat{\sigma}$ (type sqrt)
10.
A normal fit to $n = 5$ points has $\hat{\sigma}^2 = 3.6$ . Find the entry $\frac{\partial^2 \ell}{\partial \mu^2}$ of the Hessian at the peak
11.
A coin shows no heads in 8 tosses. Find $\hat{p}$
12.
An optimiser searches over $a = \ln\sigma$ and returns $\hat{a} = 0.5$ . Find $\hat{\sigma}$ (type e^(…))

How to type answers

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

L(\theta) = \prod_i f(x_i; \theta)

Likelihood and the maximum likelihood estimator

The likelihood

Rule: Likelihood

It is not a distribution for the parameter

The log-likelihood

Rule: Log-likelihood

The maximum likelihood estimator

A coin, by calculus

Rule: Maximum likelihood, one parameter

Waiting times

Two parameters: the normal

The conditions of lesson 4.4.1

Rule: Maximum likelihood for a normal

By computer: minimise the negative log-likelihood

Rule: Numerical maximum likelihood

Where the method needs care

Common mistakes

Likelihood and the maximum likelihood estimator

Summary card

Free further reading