Estimation and Maximum Likelihood Concept Card Template (Statistics)

Estimation and Maximum Likelihood Concept Card: [Your study question or course]

A concept card for point estimation: estimators and estimates, bias, variance, mean squared error, consistency and efficiency, the method of moments, and maximum likelihood (likelihood, log-likelihood, score, analytic MLEs for Bernoulli, Poisson and normal models, numerical MLE with R optim and Python scipy.optimize, standard errors from the observed information, the likelihood ratio test and AIC). It is for students in a mathematical statistics or inference course, master's and PhD students who fit models by maximum likelihood in a thesis, and analysts who want to know what their software is doing when it says "log-likelihood". It speeds up learning by pairing every formula with runnable code on the same tiny data set, a fully worked example, two faded ones, and a self-test with a separate answer key. The worked examples use fictional data; replace them with your own. Source basis: Fisher (1922, 1925), Pearson (1894), Wilks (1938), Akaike (1974) and the NIST/SEMATECH e-Handbook, with all text and examples written fresh.

Document control

The idea in one paragraph

An estimator is a rule that turns data into a guess about an unknown parameter; the number it gives for your sample is the estimate. Good rules are, on average, close to the truth (low bias), do not jump around much from sample to sample (low variance), and home in on the truth as the sample grows (consistency). Maximum likelihood is the most widely used rule: choose the parameter value under which the observed data would have been most probable. It gives sensible estimates for almost any model you can write down, comes with standard errors from the curvature of the log-likelihood, and leads directly to likelihood ratio tests and AIC. Logistic, Poisson and mixed models are all fitted this way.

Intuition

Imagine 14 of 20 Lakemoor University (fictional) students pass a statistics quiz, and you want the pass probability p. Try candidate values: if p were 0.4, getting 14 passes would be fairly unlikely; if p were 0.95, getting as many as 6 failures would be very unlikely. The value that makes the observed result most likely is p = 0.7, the peak of the log-likelihood curve below. A sharp peak means the data rule out nearby values quickly (precise estimate, small SE); a flat peak means many values fit almost equally well (large SE). That curvature is the "information" in the data.

A dartboard picture separates bias and variance. Bias is how far the centre of your cloud of darts is from the bullseye; variance is how spread out the cloud is. MSE combines the two, and sometimes a slightly biased thrower with a tight cloud beats an unbiased thrower with a wide one, as the simulation later shows.

Formal definition

Estimators and their properties

An estimator $\hat{\theta}$ is a function of the random sample, so it has its own sampling distribution. Its bias is the gap between its average and the truth; its variance is its spread; MSE combines them.

$\text{Bias}(\hat{\theta}) = E[\hat{\theta}] - \theta, \qquad \text{MSE}(\hat{\theta}) = E\left[(\hat{\theta} - \theta)^2\right] = \text{Var}(\hat{\theta}) + \text{Bias}(\hat{\theta})^2$

• Unbiased: bias is zero for every θ.
• Consistent: $\hat{\theta}$ gets arbitrarily close to θ with high probability as n grows.
• Efficient: among unbiased estimators, it has the smallest variance. The Cramér-Rao bound says no unbiased estimator can have variance below $1/I(\theta)$, where $I(\theta)$ is the Fisher information.

Method of moments

Set sample moments equal to their population formulas and solve for the parameters (Pearson 1894). It is quick and often a good starting value for numerical MLE. For a gamma distribution with shape k and rate r, the mean is k/r and the variance is k/r², so

$\hat{k}_{MM} = \frac{\bar{x}^2}{\hat{\sigma}^2}, \qquad \hat{r}_{MM} = \frac{\bar{x}}{\hat{\sigma}^2}, \qquad \hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2$

Maximum likelihood

The likelihood is the joint probability (or density) of the observed data, read as a function of the parameter. Logs turn the product into a sum, which is easier to maximise and numerically stable.

$L(\theta) = \prod_{i=1}^n f(x_i \mid \theta), \qquad \ell(\theta) = \sum_{i=1}^n \log f(x_i \mid \theta), \qquad \hat{\theta}_{ML} = \arg\max_{\theta} \ell(\theta)$

The score is the slope of the log-likelihood; the MLE usually solves score = 0. The observed information is the negative curvature at the peak; its inverse estimates the variance of the MLE.

$U(\theta) = \frac{\partial \ell}{\partial \theta}, \qquad J(\hat{\theta}) = -\frac{\partial^2 \ell}{\partial \theta^2}\Big|_{\hat{\theta}}, \qquad SE(\hat{\theta}) \approx \frac{1}{\sqrt{J(\hat{\theta})}}$

With several parameters, $J$ is the Hessian matrix of $-\ell$ and the SEs are the square roots of the diagonal of its inverse. Under regular conditions (parameter not on a boundary, model correct, n reasonably large) the MLE is consistent, approximately normal and asymptotically efficient, and it is invariant: the MLE of $g(\theta)$ is $g(\hat{\theta})$.

Analytic MLEs

Bernoulli (x successes in n trials): $\ell(p) = x\log p + (n - x)\log(1-p)$. Setting $U(p) = x/p - (n-x)/(1-p) = 0$ gives

$\hat{p} = \frac{x}{n}, \qquad J(\hat{p}) = \frac{n}{\hat{p}(1-\hat{p})}, \qquad SE(\hat{p}) = \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}$

Poisson (counts $x_1, \dots, x_n$): $\ell(\lambda) = \sum x_i \log\lambda - n\lambda - \sum \log x_i!$. Setting $U(\lambda) = \sum x_i/\lambda - n = 0$ gives

$\hat{\lambda} = \bar{x}, \qquad J(\hat{\lambda}) = \frac{\sum x_i}{\hat{\lambda}^2} = \frac{n}{\hat{\lambda}}, \qquad SE(\hat{\lambda}) = \sqrt{\frac{\hat{\lambda}}{n}}$

Normal (mean μ, variance σ²): solving both score equations gives

$\hat{\mu} = \bar{x}, \qquad \hat{\sigma}^2_{ML} = \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2$

The MLE of σ² divides by n, not n − 1, so it is biased downwards by the factor (n − 1)/n. It is still consistent, because that factor tends to

1. For the 12 Brightwell Bakery loaf weights used in other cards, R gives $\hat{\mu}$ = 497.55, $\hat{\sigma}_{ML}$ = 2.1654 and s = 2.2617.

Likelihood ratio test and AIC

To compare a simpler (nested) model with a fuller one, compare their maximised log-likelihoods (Wilks 1938). Under H0 (the simpler model is true) the statistic is approximately chi-square with df equal to the number of extra parameters.

$\Lambda = 2\left[\ell(\hat{\theta}_{full}) - \ell(\hat{\theta}_{0})\right] \sim \chi^2_{df}$

AIC (Akaike 1974) trades fit against complexity for any set of models fitted to the same data, nested or not; smaller is better.

$\text{AIC} = -2\,\ell(\hat{\theta}) + 2k$

Worked example A: fully solved

Question. Harbour City Library (fictional) logs lost-property reports per day. Over ten days the Central branch had the counts below. Estimate the daily rate, give an SE and 95% CI, and test whether the Harbour branch has the same rate.

Step 1: model. Counts of independent events per day at a constant rate: Poisson. Step 2: analytic MLE. From the derivation, $\hat{\lambda} = \bar{x}$ = 38 / 10 = 3.8 reports per day. Step 3: SE. $SE = \sqrt{3.8/10}$ = 0.616. Step 4: Wald CI. 3.8 ± 1.96 × 0.616 = [2.59, 5.01].

Step 5: likelihood ratio test. H0: both branches share one rate (one parameter, $\hat{\lambda}$ = 52 / 20 = 2.6). Full model: separate rates (two parameters, 3.8 and 1.4). One extra parameter, so df = 1. Outputs on this card are trimmed to the key lines (for example the header of R's analysis-of-deviance table is dropped).

central <- c(3, 5, 2, 4, 6, 3, 4, 5, 2, 4)
harbour <- c(1, 2, 0, 3, 1, 2, 2, 1, 0, 2)
lam <- mean(central); c(lambda_hat = lam, se = sqrt(lam / length(central)))
ll0 <- sum(dpois(c(central, harbour), mean(c(central, harbour)), log = TRUE))
ll1 <- sum(dpois(central, mean(central), log = TRUE)) +
       sum(dpois(harbour, mean(harbour), log = TRUE))
LR <- 2 * (ll1 - ll0)
round(c(ll0 = ll0, ll1 = ll1, LR = LR, p = pchisq(LR, 1, lower.tail = FALSE)), 5)
# The same test as a Poisson GLM
branch <- rep(c("central", "harbour"), each = 10)
anova(glm(c(central, harbour) ~ branch, family = poisson), test = "LRT")

lambda_hat         se
 3.8000000  0.6164414
      ll0       ll1        LR         p
-37.53596 -31.78191  11.50811   0.00069
       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)
NULL                      19     23.998
branch  1   11.508        18     12.490 0.0006929 ***

import numpy as np
import statsmodels.api as sm
from scipy import stats

central = np.array([3, 5, 2, 4, 6, 3, 4, 5, 2, 4])
harbour = np.array([1, 2, 0, 3, 1, 2, 2, 1, 0, 2])
print(central.mean(), round(np.sqrt(central.mean() / 10), 4))
y = np.concatenate([central, harbour]); grp = np.repeat([0, 1], 10)
m0 = sm.GLM(y, np.ones(20), family=sm.families.Poisson()).fit()
m1 = sm.GLM(y, sm.add_constant(grp), family=sm.families.Poisson()).fit()
LR = 2 * (m1.llf - m0.llf)
print(round(m0.llf, 5), round(m1.llf, 5), round(LR, 5), round(stats.chi2.sf(LR, 1), 5))

3.8 0.6164
-37.53596 -31.78191 11.50811 0.00069

Step 6: interpret. The Central branch logs about 3.8 reports a day (95% CI 2.59 to 5.01). Allowing separate rates raises the log-likelihood by 5.75, giving Λ = 11.51 on 1 df, p < .001: the branches' rates differ. The GLM deviance drop (23.998 − 12.490 = 11.508) is exactly the same LR statistic, because deviance is −2 times the log-likelihood up to a constant.

Report: A Poisson model estimated 3.8 lost-property reports per day at Central (95% CI [2.59, 5.01]) and 1.4 at Harbour; a likelihood ratio test favoured separate rates, χ²(1) = 11.51, p < .001.

Check before trusting it. The Poisson model assumes the variance equals the mean. If counts are more spread out (overdispersion), the SEs are too small; see stats-count-regression.

Worked example B: half solved

Question. Fifteen waiting times (minutes) at Lakemoor University's IT help desk (fictional): 2.1, 4.8, 1.3, 3.6, 7.2, 2.9, 5.5, 1.8, 3.1, 9.4, 2.4, 4.0, 6.1, 1.6, 3.3. Waiting times are positive and right-skewed, so fit a gamma distribution. The gamma shape has no closed-form MLE, so we maximise numerically, starting from the method of moments.

x <- c(2.1, 4.8, 1.3, 3.6, 7.2, 2.9, 5.5, 1.8, 3.1, 9.4, 2.4, 4.0, 6.1, 1.6, 3.3)
v <- mean((x - mean(x))^2)
mom <- c(shape = mean(x)^2 / v, rate = mean(x) / v); round(mom, 4)
# Work on the log scale so that optim cannot propose negative values
nll <- function(par) -sum(dgamma(x, shape = exp(par[1]), rate = exp(par[2]), log = TRUE))
fit <- optim(log(mom), nll, method = "BFGS", hessian = TRUE)
est <- exp(fit$par); round(c(est, logLik = -fit$value), 4)
se <- est * sqrt(diag(solve(fit$hessian)))   # delta method back to the natural scale
round(se, 4)
ll_exp <- sum(dexp(x, 1 / mean(x), log = TRUE))           # exponential = gamma with shape 1
round(c(AIC_exp = -2 * ll_exp + 2, AIC_gamma = 2 * fit$value + 4), 3)
LR <- 2 * (-fit$value - ll_exp); round(c(LR = LR, p = pchisq(LR, 1, lower.tail = FALSE)), 4)

 shape   rate
3.1734 0.8054
   shape     rate   logLik
  3.4524   0.8762 -31.0026
 shape   rate
1.2047 0.3291
  AIC_exp AIC_gamma
   73.135    66.005
    LR      p
9.1302 0.0025

import numpy as np
from scipy import stats, optimize
from statsmodels.tools.numdiff import approx_hess

x = np.array([2.1, 4.8, 1.3, 3.6, 7.2, 2.9, 5.5, 1.8, 3.1, 9.4, 2.4, 4.0, 6.1, 1.6, 3.3])
v = x.var()                                   # ddof=0, the moment estimator
mom = np.array([x.mean()**2 / v, x.mean() / v]); print(np.round(mom, 4))
def nll(par):
    shape, rate = par
    if shape <= 0 or rate <= 0:
        return np.inf
    return -stats.gamma.logpdf(x, a=shape, scale=1 / rate).sum()
fit = optimize.minimize(nll, mom, method="Nelder-Mead", options={"xatol": 1e-8, "fatol": 1e-10})
print(np.round(fit.x, 4), round(-fit.fun, 4))
H = approx_hess(fit.x, nll)                   # observed information matrix
print(np.round(np.sqrt(np.diag(np.linalg.inv(H))), 4))
a, loc, scale = stats.gamma.fit(x, floc=0); print(round(a, 4), round(1 / scale, 4))
ll_exp = stats.expon.logpdf(x, scale=x.mean()).sum()
print("AIC exp", round(-2 * ll_exp + 2, 3), "AIC gamma", round(2 * fit.fun + 4, 3))

[3.1734 0.8054]
[3.4524 0.8762] -31.0026
[1.2047 0.3291]
3.4524 0.8762
AIC exp 73.135 AIC gamma 66.005

Note the parameterisations: R's dgamma takes a rate, while scipy's gamma takes a scale = 1 / rate; scipy's fit also estimates a location unless you fix floc=0. Mixing them up is a classic error.

1. Method-of-moments estimates: shape = [ ], rate = [ ].
2. MLE: shape = [ ], rate = [ ]; maximised log-likelihood = [ ].
3. Implied mean waiting time = shape / rate = [ ] minutes. Compare it with $\bar{x}$ = 3.94: [same / different]? (Hint: the gamma MLE always matches the sample mean.)
4. 95% Wald CI for the shape: 3.4524 ± 1.96 × [ ] = [ ].
5. Exponential versus gamma: Λ = [ ] on [ ] df, p = [ ]; AIC prefers [ ] by [ ] points.
6. Report: "A gamma model (shape [ ], rate [ ] per minute) fitted the waiting times better than an exponential model, χ²(1) = [ ], p = [ ], AIC [ ] versus [ ]."

Exercise C: unsolved

1. Brightwell Bakery (fictional) records the minutes between customer arrivals: 2, 4, 3, 7, 4. Assume an exponential model with density $f(x) = \lambda e^{-\lambda x}$. Derive the MLE of λ, compute it, and give its SE from the observed information.
2. Compute the maximised log-likelihood.
3. For the quiz example (14 passes in 20), compute the 95% Wald CI for p and the likelihood ratio statistic for H0: p = 0.5 (use the log-likelihood values −12.217 at 0.7 and −13.863 at 0.5).
4. By the invariance property, what is the MLE of the odds p / (1 − p)?

Simulation you can run

Bias versus MSE. Draw 20,000 samples of n = 5 from a normal population with σ² = 4. Compare the MLE of the variance (divide by n) with the unbiased version (divide by n − 1).

set.seed(2026)
sims <- replicate(20000, {
  y <- rnorm(5, 0, 2)
  c(mle = mean((y - mean(y))^2), unbiased = var(y))
})
round(rowMeans(sims), 3)              # average estimate (truth = 4)
round(rowMeans((sims - 4)^2), 3)      # MSE

     mle unbiased
   3.208    4.010
     mle unbiased
   5.750    8.004

import numpy as np
rng = np.random.default_rng(2026)
ys = rng.normal(0, 2, size=(20000, 5))
mle, unb = ys.var(axis=1, ddof=0), ys.var(axis=1, ddof=1)
print(np.round([mle.mean(), unb.mean()], 3), np.round([((mle - 4)**2).mean(), ((unb - 4)**2).mean()], 3))

[3.188 3.986] [5.724 7.915]

Theory predicts an average of $\frac{n-1}{n}\sigma^2$ = 3.2 for the MLE and 4 for the unbiased estimator, with MSEs of 5.76 and 8.00. Both simulations land close (R and Python use different random generators, so their numbers differ slightly). The lesson: the MLE is biased, yet it has the smaller mean squared error, because its lower variance more than pays for its bias. "Unbiased" is not the same as "best".

Misconceptions

Connections to other concepts

• Distributions (stats-distributions-sheet): every likelihood starts from a distribution's pmf or pdf.
• Sampling distributions (stats-sampling-distributions): an estimator's bias and variance describe its sampling distribution.
• Confidence intervals (stats-confidence-intervals): Wald intervals are estimate ± 1.96 × SE from the information; profile-likelihood intervals invert the LR test.
• Hypothesis testing (stats-hypothesis-testing): Wald, score and likelihood ratio tests are three large-sample tests built from the same log-likelihood.
• Regression models (stats-logistic-regression, stats-count-regression, stats-mixed-models): all fitted by maximum likelihood (or REML), and their output reports log-likelihood, AIC and deviance.
• Bayesian inference (stats-bayesian-inference): posterior ∝ prior × likelihood; with a flat prior the posterior mode equals the MLE.

Self-test

1. Define bias, variance and MSE, and state the identity that links them.
2. Derive the MLE of p for x successes in n Bernoulli trials.
3. Why does a sharper log-likelihood peak imply a smaller SE?
4. What is the MLE of σ² for a normal sample, is it biased, and is it consistent?
5. Two nested models have log-likelihoods −120.4 (3 parameters) and −116.1 (5 parameters). Compute Λ, its df, and the AIC of each model.
6. Give one reason to start a numerical optimiser from method-of-moments estimates.

My notes and next steps

Derivation I can now do from memory: [list them, for example Poisson MLE and its SE].

Sources and further reading

Structure aligned with a standard first course in mathematical statistics and with the GAISE College Report (ASA 2016) emphasis on simulation; based on the original papers below. All wording and examples are original.

• Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716–723.
• Fisher, R. A. (1922). On the mathematical foundations of theoretical statistics. Philosophical Transactions of the Royal Society of London, Series A, 222, 309–368.
• Fisher, R. A. (1925). Theory of statistical estimation. Proceedings of the Cambridge Philosophical Society, 22(5), 700–725.
• Pearson, K. (1894). Contributions to the mathematical theory of evolution. Philosophical Transactions of the Royal Society of London, Series A, 185, 71–110.
• Wilks, S. S. (1938). The large-sample distribution of the likelihood ratio for testing composite hypotheses. Annals of Mathematical Statistics, 9(1), 60–62.
• Nelder, J. A., and Mead, R. (1965). A simplex method for function minimization. The Computer Journal, 7(4), 308–313.
• NIST/SEMATECH e-Handbook of Statistical Methods, section 1.3.6.5.2 Maximum Likelihood, https://www.itl.nist.gov/div898/handbook/ (accessed 2026-10-25). Public domain.
• GAISE College Report ASA Revision Committee (2016). Guidelines for Assessment and Instruction in Statistics Education (GAISE) College Report 2016. American Statistical Association.
• Next steps: Pawitan, Y. (2001). In All Likelihood: Statistical Modelling and Inference Using Likelihood. Oxford University Press; Casella, G., and Berger, R. L. (2002). Statistical Inference (2nd ed.). Duxbury.
• Further reading: de Smith, M. J. (2024). Statistical Analysis Handbook, 2024 edition. StatsRef, https://www.statsref.com/; OpenIntro Statistics (CC BY-SA).

Answer key

Example B.

1. Shape 3.1734, rate 0.8054.
2. Shape 3.4524, rate 0.8762; log-likelihood −31.0026.
3. 3.4524 / 0.8762 = 3.94 minutes, the same as $\bar{x}$.
4. 3.4524 ± 1.96 × 1.2047 = [1.09, 5.81].
5. Λ = 9.13 on 1 df, p = .0025; AIC prefers the gamma model by 7.13 points (66.005 versus 73.135).
6. "A gamma model (shape 3.45, rate 0.88 per minute) fitted the waiting times better than an exponential model, χ²(1) = 9.13, p = .003, AIC 66.0 versus 73.1."

Exercise C.

1. $\ell(\lambda) = n\log\lambda - \lambda\sum x_i$; $U = n/\lambda - \sum x_i = 0$ gives $\hat{\lambda} = 1/\bar{x}$ = 1 / 4 = 0.25 per minute. $J = n/\hat{\lambda}^2$, so $SE = \hat{\lambda}/\sqrt{n}$ = 0.25 / √5 = 0.112.
2. $5\log 0.25 - 0.25 \times 20$ = −11.93.
3. Wald CI 0.7 ± 1.96 × 0.1025 = [0.499, 0.901]. Λ = 2 × (−12.217 − (−13.863)) = 3.29 on 1 df, p = 0.070: not significant at α = 0.05.
4. 0.7 / 0.3 = 2.33.

Self-test.

1. Bias = E[estimator] − truth; variance = spread of the estimator across samples; MSE = average squared error = variance + bias².
2. $\ell(p) = x\log p + (n-x)\log(1-p)$; derivative $x/p - (n-x)/(1-p) = 0$ gives $\hat{p} = x/n$.
3. The SE is roughly one over the square root of the curvature; a sharp peak means high curvature (much information), so the SE is small.
4. $\frac{1}{n}\sum(x_i - \bar{x})^2$; biased (expectation $\frac{n-1}{n}\sigma^2$) but consistent.
5. Λ = 2 × (−116.1 + 120.4) = 8.6 on 2 df (p ≈ 0.014). AIC = 246.8 (3 parameters) and 242.2 (5 parameters): the larger model is preferred.
6. Good starting values speed convergence and reduce the risk of the optimiser wandering into an invalid region or a local optimum.

Revision history