28. Maximum likelihood estimation

The method of moments asks "which parameter reproduces the sample's moments?" Maximum likelihood asks a better question:

Which parameter value makes the data we actually observed most probable?

It is the dominant estimation method in statistics, and almost everything in the rest of this course is built on it.

The likelihood function

Given data x_1, \dots, x_n that are independent with density (or mass function) f(x \mid \theta), the likelihood is

L(\theta) = \prod_{i=1}^n f(x_i \mid \theta)

This is the same formula as the joint density, read differently. The joint density treats \theta as fixed and x as varying. The likelihood fixes the observed x and treats \theta as the variable.

That distinction matters: L(\theta) is not a probability distribution over \theta. It doesn't integrate to 1 and it has no area interpretation. It's a ranking of parameter values by how well each explains the data.

The maximum likelihood estimator is the value that maximises it:

\hat\theta_{\text{MLE}} = \arg\max_\theta L(\theta)

Always use the log

Products are miserable to differentiate and they underflow numerically — 1,000 densities of 0.1 each multiply to 10^{-1000}, which is zero in floating point.

Take logs:

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

Since \log is strictly increasing, the maximiser is unchanged. Sums are easy to differentiate and numerically stable.

The recipe:

1. Write the likelihood.
2. Take the log.
3. Differentiate with respect to \theta and set to zero (the score equation).
4. Solve, and check it's a maximum (second derivative negative).
5. Check the boundaries of the parameter space separately.

Step 5 is the one people skip, and it's where the Uniform example below bites.

Worked example 1: Bernoulli

x_1, \dots, x_n are 0/1 with \sum x_i = k successes. Estimate p.

L(p) = \prod p^{x_i}(1-p)^{1-x_i} = p^k(1-p)^{n-k}

\ell(p) = k\log p + (n-k)\log(1-p)

Differentiate and set to zero:

\ell'(p) = \frac{k}{p} - \frac{n-k}{1-p} = 0

k(1-p) = (n-k)p \implies k = np \implies \hat p = \frac{k}{n}

The sample proportion — reassuringly obvious, now derived rather than assumed.

Worked example 2: Normal

Estimate \mu and \sigma^2.

\ell(\mu, \sigma^2) = -\frac{n}{2}\log(2\pi) - \frac{n}{2}\log\sigma^2 - \frac{1}{2\sigma^2}\sum(x_i - \mu)^2

For \mu, only the last term matters:

\frac{\partial\ell}{\partial\mu} = \frac{1}{\sigma^2}\sum(x_i - \mu) = 0 \implies \hat\mu = \bar x

Note this also shows the MLE of \mu minimises \sum(x_i - \mu)^2 — least squares falls out of maximum likelihood under Normal errors, which is the link to §6.

For \sigma^2:

\frac{\partial\ell}{\partial\sigma^2} = -\frac{n}{2\sigma^2} + \frac{1}{2\sigma^4}\sum(x_i - \bar x)^2 = 0

\hat\sigma^2 = \frac{1}{n}\sum(x_i - \bar x)^2

Divide by n again — the MLE of the variance is biased. Maximum likelihood does not promise unbiasedness.

Worked example 3: Uniform — where calculus fails

X_i \sim \text{Uniform}(0, \theta).

L(\theta) = \prod_{i=1}^n \frac{1}{\theta}\mathbb{1}\{0 \le x_i \le \theta\} = \frac{1}{\theta^n}\mathbb{1}\{\theta \ge x_{(n)}\}

Differentiating gives -n/\theta^{n+1}, which is never zero. The score equation has no solution.

But look at the function. It is zero for \theta < x_{(n)} (impossible — you observed a value that large), then jumps up and decreases as 1/\theta^n. So it is maximised at the smallest permissible \theta:

\hat\theta_{\text{MLE}} = x_{(n)}

The maximum. Exactly the estimator that beat the method of moments two lessons ago on MSE. Maximum likelihood found it automatically, because it uses the whole likelihood rather than a couple of moments.

The lesson: always check the boundary. Setting the derivative to zero finds interior maxima only.

Why MLE is the default

Consistency. \hat\theta_{\text{MLE}} \xrightarrow{p} \theta under mild regularity conditions.

Asymptotic normality.

\sqrt n(\hat\theta - \theta) \xrightarrow{d} N\!\left(0, \frac{1}{I(\theta)}\right)

where I(\theta) is the Fisher information (next lesson but one). This gives standard errors and confidence intervals essentially for free.

Asymptotic efficiency. No consistent estimator has smaller asymptotic variance — the MLE attains the Cramér–Rao lower bound. It extracts all the information the data contains.

Invariance. If \hat\theta is the MLE of \theta, then g(\hat\theta) is the MLE of g(\theta), for any function g. Want the MLE of \sigma? Take \sqrt{\hat\sigma^2}. No re-derivation.

This last property is genuinely special — unbiasedness has nothing like it, as we saw when s failed to be unbiased for \sigma despite s^2 being unbiased for \sigma^2.

The costs: often biased in small samples, sometimes requires numerical optimisation, and it can be sensitive to model misspecification — the estimate is only as good as the assumed f(x \mid \theta).

Doing it in Python

The Bernoulli and Normal cases against their closed forms:

import numpy as np
from scipy.optimize import minimize_scalar

rng = np.random.default_rng(0)
data = (rng.random(1_000) < 0.3).astype(float)     # Bernoulli(0.3)

def neg_loglik(p):
    if not 0 < p < 1:
        return np.inf
    return -(data.sum() * np.log(p) + (len(data) - data.sum()) * np.log(1 - p))

opt = minimize_scalar(neg_loglik, bounds=(1e-9, 1 - 1e-9), method="bounded")
print("numerical MLE :", round(opt.x, 6))
print("closed form k/n:", round(data.mean(), 6))

The Uniform case, showing why the derivative approach fails:

import numpy as np

rng = np.random.default_rng(1)
theta_true = 5.0
x = rng.uniform(0, theta_true, size=50)

def loglik(theta):
    if theta < x.max():
        return -np.inf              # the data rules this out
    return -len(x) * np.log(theta)

grid = np.linspace(x.max() - 1, x.max() + 3, 9)
print(f"{'theta':>10} {'log-likelihood':>18}")
for t in grid:
    ll = loglik(t)
    print(f"{t:>10.4f} {ll:>18.4f}" + ("   <- observed max" if abs(t - x.max()) < 1e-9 else ""))

print(f"\nsample max (the MLE): {x.max():.6f}")
print("The likelihood jumps from -inf and then falls -- the peak is at the boundary,")
print("so setting the derivative to zero finds nothing.")

Maximising a two-parameter likelihood numerically, the general pattern:

import numpy as np
from scipy.optimize import minimize

rng = np.random.default_rng(2)
true_mu, true_sigma = 3.0, 2.0
x = rng.normal(true_mu, true_sigma, size=2_000)

def neg_loglik(params):
    mu, log_sigma = params                  # optimise log sigma to keep it positive
    sigma = np.exp(log_sigma)
    return -np.sum(-np.log(sigma) - 0.5 * ((x - mu) / sigma) ** 2)

res = minimize(neg_loglik, x0=[0.0, 0.0])
mu_hat, sigma_hat = res.x[0], np.exp(res.x[1])

print(f"numerical: mu {mu_hat:.4f}  sigma {sigma_hat:.4f}")
print(f"closed   : mu {x.mean():.4f}  sigma {x.std(ddof=0):.4f}")
print(f"true     : mu {true_mu}       sigma {true_sigma}")

And the invariance property, which has no analogue for unbiased estimators:

import numpy as np

rng = np.random.default_rng(3)
x = rng.normal(0, 3.0, size=10_000)

var_mle = x.var(ddof=0)
print("MLE of sigma^2 :", round(var_mle, 4))
print("MLE of sigma   :", round(np.sqrt(var_mle), 4), " (just the square root)")
print("MLE of 1/sigma :", round(1 / np.sqrt(var_mle), 4), " (just the reciprocal)")
print("\nAny function of the MLE is the MLE of that function -- no re-derivation.")

Your turn

1. X_i \sim \text{Exponential}(\lambda). Find the MLE of \lambda.

2. X_i \sim \text{Poisson}(\lambda). Find the MLE of \lambda.

3. Given the MLE of \sigma^2 is \frac{1}{n}\sum(x_i - \bar x)^2, what is the MLE of \sigma?

Solutions

1. f(x \mid \lambda) = \lambda e^{-\lambda x}, so

\ell(\lambda) = n\log\lambda - \lambda\sum x_i

\ell'(\lambda) = \frac{n}{\lambda} - \sum x_i = 0 \implies \hat\lambda = \frac{n}{\sum x_i} = \frac{1}{\bar x}

Check the second derivative: \ell''(\lambda) = -n/\lambda^2 < 0, so it's a maximum. ✓

Same as the method of moments estimator here — the two agree for some distributions and differ for others.

2. P(X = x) = \frac{\lambda^x e^{-\lambda}}{x!}, so

\ell(\lambda) = \left(\sum x_i\right)\log\lambda - n\lambda - \sum\log(x_i!)

The last term has no \lambda in it and can be dropped — a useful habit, since constants never affect the maximiser.

\ell'(\lambda) = \frac{\sum x_i}{\lambda} - n = 0 \implies \hat\lambda = \bar x

Sensible, since \lambda is the mean of a Poisson.

3. By invariance, with g(t) = \sqrt t:

\hat\sigma = \sqrt{\hat\sigma^2} = \sqrt{\frac{1}{n}\sum(x_i - \bar x)^2}

No new derivation needed.

Note this is biased for \sigma — but so is s = \sqrt{s^2}, the square root of the unbiased variance estimator. Neither is unbiased for \sigma, because square roots don't commute with expectation. Invariance is a property maximum likelihood has and unbiasedness simply does not.

Check yourself in code

Find the MLE of a Uniform(0,\theta) upper endpoint by direct search over a grid, and confirm it lands on the sample maximum rather than anywhere the derivative would suggest.

Print exactly this:

sample max 4.9037
grid MLE 4.9037
matches sample max: True

Use default_rng(1), \theta = 5, 50 samples. Search a grid of 100000 points from the sample max to the sample max plus 3. Round both values to 4 decimal places.

import numpy as np

rng = np.random.default_rng(1)
x = rng.uniform(0, 5.0, size=50)

print("sample max", round(x.max(), 4))

# Evaluate the log-likelihood -n*log(theta) on a grid from x.max() to
# x.max() + 3, take the argmax, and compare it with the sample maximum.

import numpy as np

rng = np.random.default_rng(1)
x = rng.uniform(0, 5.0, size=50)

print("sample max", round(x.max(), 4))

grid = np.linspace(x.max(), x.max() + 3, 100_000)
loglik = -len(x) * np.log(grid)
best = grid[np.argmax(loglik)]

print("grid MLE", round(best, 4))
print("matches sample max:", round(best, 4) == round(x.max(), 4))

---

Maximum likelihood picks the parameter under which the observed data was most probable. Work with the log, differentiate, solve — and check the boundary, because the peak isn't always interior. It's consistent, asymptotically Normal, asymptotically efficient, and invariant under transformation, which is why it's the default. It is not, however, unbiased.

Next: why some statistics carry all the information in a sample, so the rest of the data can be thrown away.