Maximum likelihood
Which version of the world makes your data least surprising?
Someone tosses a coin 500 times and gets 97 heads. Is the coin fair?
If it were fair (p = 0.5), the chance of seeing exactly 97 heads is about 9.3 \times 10^{-46}. Strictly mathematically speaking: not a whole lot. So the coin is almost certainly not fair. But then what is its probability of heads?
Maximum likelihood answers that by turning the question around. Instead of asking “given p, how likely is the data?”, we hold the data fixed and ask: which value of p would have made what we actually saw most likely? That value is the maximum likelihood estimate (MLE).
same formula, read backwards
The rest of this lesson shows that this one idea gives you logistic regression, standard errors, and even OLS as special cases.
The likelihood curve
For n tosses with k heads, the probability of the data is
L(p) = \binom{n}{k}\, p^{k}\,(1-p)^{n-k}.
Read as a function of p with the data fixed, this is the likelihood. Slide p in the app and watch how likely the observed 97 heads becomes. The red line is the slope at your current p; the slope tells you which way is uphill.
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 600
library(shiny)
library(bslib)
ink <- "#1B2559"; red <- "#D6334A"; blue <- "#2463EB"; grey <- "#9AA3B5"
base_par <- function() {
par(mar = c(4, 4.6, 1.6, 1), las = 1, bty = "l", col.axis = ink,
col.lab = ink, fg = grey, cex.axis = 0.85, col.main = ink,
font.main = 1, cex.main = 0.95)
}
sci <- function(v) {
if (v == 0) return("0")
e <- floor(log10(abs(v))); m <- v / 10^e
if (abs(e) < 4) return(formatC(v, digits = 4, format = "g"))
sprintf("%.2f × 10^%d", m, e)
}
ui <- page_sidebar(
sidebar = sidebar(
width = 260,
sliderInput("n", "Number of tosses n", min = 10, max = 1000, value = 500, step = 10),
sliderInput("k", "Heads observed k", min = 1, max = 499, value = 97, step = 1),
sliderInput("p", "Your guess for p", min = 0.001, max = 0.999, value = 0.5, step = 0.001),
radioButtons("scale", "Show", inline = FALSE,
choices = c("Log-likelihood ℓ(p)" = "log", "Likelihood L(p)" = "raw")),
actionButton("jump", "Jump to the top", class = "btn-primary")
),
layout_columns(
col_widths = c(7, 5),
card(card_header("How likely is the data, for each p?"),
plotOutput("curve", height = "300px"), uiOutput("readout")),
card(card_header("Zoom on the peak: sharper means more precise"),
plotOutput("zoom", height = "300px"), uiOutput("precision"))
)
)
server <- function(input, output, session) {
observeEvent(input$n, {
updateSliderInput(session, "k", max = input$n - 1,
value = min(input$k, input$n - 1))
})
observeEvent(input$jump, updateSliderInput(session, "p", value = round(input$k / input$n, 3)))
ll <- function(p, n, k) dbinom(k, n, p, log = TRUE)
score <- function(p, n, k) k / p - (n - k) / (1 - p)
output$curve <- renderPlot({
n <- input$n; k <- min(input$k, n - 1); p <- input$p
ps <- seq(0.001, 0.999, length.out = 600)
phat <- k / n
base_par()
if (input$scale == "log") {
y <- ll(ps, n, k); yp <- ll(p, n, k); slope <- score(p, n, k)
lo <- max(min(y), ll(phat, n, k) - 400)
plot(ps, y, type = "l", lwd = 2.5, col = ink, ylim = c(lo, ll(phat, n, k) + 5),
xlab = "p (probability of heads)", ylab = "log-likelihood")
} else {
y <- dbinom(k, n, ps); yp <- dbinom(k, n, p); slope <- yp * score(p, n, k)
plot(ps, y, type = "l", lwd = 2.5, col = ink,
xlab = "p (probability of heads)", ylab = "likelihood")
}
abline(v = phat, col = blue, lty = 2)
mtext(sprintf("top: p̂ = %.3f", phat), side = 3, at = phat, col = blue, cex = 0.85, line = 0.2)
h <- 0.07
segments(p - h, yp - slope * h, p + h, yp + slope * h, col = red, lwd = 2.5)
points(p, yp, pch = 19, col = red, cex = 1.5)
})
output$readout <- renderUI({
n <- input$n; k <- min(input$k, n - 1); p <- input$p
s <- score(p, n, k)
direction <- if (abs(s) < 2) "The slope is about zero: you are at the top. This p is the MLE."
else if (s > 0) "The slope is positive: uphill is to the right. Increase p."
else "The slope is negative: uphill is to the left. Decrease p."
tags$div(style = "font-size:.92rem;padding:0 .3rem;",
tags$p(style = "margin:.2rem 0;",
sprintf("At p = %.3f: L(p) = %s, ℓ(p) = log L(p) = %.1f", p, sci(dbinom(k, n, p)), ll(p, n, k))),
tags$p(style = "margin:.2rem 0;",
"Slope of ℓ (the score): ", tags$strong(style = paste0("color:", red), sprintf("%.1f", s))),
tags$p(style = "margin:.2rem 0;color:#4A5578;", direction))
})
output$zoom <- renderPlot({
n <- input$n; k <- min(input$k, n - 1); phat <- k / n
se <- sqrt(phat * (1 - phat) / n)
ps <- seq(max(0.001, phat - 5 * se), min(0.999, phat + 5 * se), length.out = 400)
rel <- ll(ps, n, k) - ll(phat, n, k)
rel4 <- ll(ps, 4 * n, 4 * k) - ll(phat, 4 * n, 4 * k)
base_par()
plot(ps, rel, type = "l", lwd = 2.5, col = ink, ylim = c(-17, 0.5),
xlab = "p", ylab = "ℓ(p) − ℓ(p̂)")
lines(ps, rel4, lwd = 2, lty = 2, col = blue)
abline(h = -1.92, col = grey, lty = 3)
ci <- range(ps[rel >= -1.92])
segments(ci[1], -1.92, ci[2], -1.92, col = red, lwd = 4)
legend("bottom", bty = "n", cex = 0.8, lwd = c(2.5, 2), lty = c(1, 2),
col = c(ink, blue), text.col = ink,
legend = c(sprintf("your data (n = %d)", n), sprintf("4× more data (n = %d)", 4 * n)))
})
output$precision <- renderUI({
n <- input$n; k <- min(input$k, n - 1); phat <- k / n
H <- n / (phat * (1 - phat))
se <- 1 / sqrt(H)
tags$div(style = "font-size:.92rem;padding:0 .3rem;",
tags$p(style = "margin:.2rem 0;", sprintf("Curvature at the top: H = %.0f", H)),
tags$p(style = "margin:.2rem 0;", "Standard error = 1/√H = ",
tags$strong(sprintf("%.4f", se))),
tags$p(style = "margin:.2rem 0;color:#D6334A;",
sprintf("95%% interval: %.3f to %.3f", phat - 1.96 * se, phat + 1.96 * se)))
})
}
shinyApp(ui, server)Try this
- Start at p = 0.5 with the log-likelihood view. The slope is strongly negative, so move left. Follow the slope until it reaches about zero. Where do you land? (The slides get 0.194, the true value being 0.2.)
- Switch to Likelihood L(p). The peak is in the same place, but the curve is a thin spike that is practically zero everywhere else: at p = 0.5 it is 9.3 \times 10^{-46}. The log view stretches those tiny numbers out so you can see the whole landscape. Same top, readable picture.
- In the right-hand panel, compare the solid curve with the dashed one (four times more tosses, same share of heads). The dashed peak is much sharper. Sharper peak = more information = smaller standard error.
- Set n = 20, k = 4. Same share of heads as before, but look how flat the peak is and how wide the interval gets.
The maths, decoded
Products are awkward, so we take logs. The log turns the product into a sum and doesn’t move the peak:
\ell(p) = \log L(p) = \text{const} + k\log p + (n-k)\log(1-p)
At the top the slope is zero. The slope of the log-likelihood has a name, the score:
S(p) = \frac{d\ell}{dp} = \frac{k}{p} - \frac{n-k}{1-p} = 0 \quad\Longrightarrow\quad \hat p = \frac{k}{n}.
So the MLE for a coin is just the share of heads. Reassuring: the fancy method agrees with common sense.
Formula decoder
| Symbol | In words |
|---|---|
| L(p) | likelihood: the probability of our data if the true value were p |
| \ell(p) | log-likelihood; easier to work with, same peak |
| S(p) | the score: the slope of \ell; positive means “go right” |
| \hat p | the MLE: the value of p where the score is zero |
| H | curvature: how sharply \ell bends at the top; its inverse gives the variance of \hat p |
Why logs? Because \log(a \times b) = \log a + \log b. A likelihood is a product of many small probabilities (one per observation). With a few thousand observations that product becomes smaller than the smallest number a computer can store, and it gets rounded to zero. Adding a few thousand logs is no problem at all. And since \log always increases, whatever p maximises L also maximises \log L.
Why derivatives? A derivative is just the slope of a curve at a point. At the peak of a hill the ground is flat, so “find where the derivative is zero” is maths-speak for “find the top”. When you can’t solve it by hand, a computer walks uphill numerically; that’s what R’s optim() does in the lab.
Curvature is information
The right-hand panel shows the key to uncertainty. If the log-likelihood is sharply peaked, only a narrow range of p values fits the data well: we know p precisely. If it is flat, many values fit almost equally well.
The curvature at the top is measured by the negative second derivative (the Hessian), and the standard error is its inverse square root:
\text{se}(\hat\theta) \approx \frac{1}{\sqrt{H_n(\hat\theta)}}, \qquad H_n(\theta) = -\frac{d^2 \ell}{d\theta^2}.
flat top → big variance
sharp top → small variance
For the coin, H = n / \big(\hat p(1-\hat p)\big), which grows with n. That’s the formal version of “more data, sharper peak”. The expected curvature per observation is called the Fisher information, J; under a correctly specified model the two ways of measuring it agree (the information matrix equality, J_\theta = H_\theta).
A real application: the Challenger disaster
On 28 January 1986 the space shuttle Challenger broke apart shortly after launch. The cause was a failure of rubber O-rings sealing joints in the booster rockets. Each flight has 6 of them. Before the launch, engineers had data from 23 earlier flights on how many O-rings showed damage, and the temperature at launch. The forecast for launch day was −0.6 °C, far colder than any previous flight.
Could the data have warned them? We want the probability p_i that a single O-ring fails at temperature \text{temp}_i.
A straight line p_i = \beta_0 + \beta_1\,\text{temp}_i doesn’t work: extend it to cold temperatures and it happily predicts probabilities above 1. Instead we pass the line through the logistic function, which squeezes any number into the range 0 to 1:
\log\frac{p_i}{1-p_i} = \beta_0 + \beta_1\,\text{temp}_i, \qquad y_i \sim \text{Binomial}(6,\, p_i).
Now the model is a complete probability story: for any (\beta_0, \beta_1) we can compute the probability of each flight’s damage count, multiply them, and get a likelihood. Maximum likelihood picks the pair that makes the 23 observed flights most likely.
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 560
library(shiny)
library(bslib)
ink <- "#1B2559"; red <- "#D6334A"; blue <- "#2463EB"; grey <- "#9AA3B5"
base_par <- function() {
par(mar = c(4, 4.4, 1.6, 1), las = 1, bty = "l", col.axis = ink,
col.lab = ink, fg = grey, cex.axis = 0.85, col.main = ink,
font.main = 1, cex.main = 0.95)
}
# Pre-Challenger flights: launch temperature (°C) and damaged O-rings out of 6
tempF <- c(53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70, 70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81)
damage <- c(5, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0)
temp <- round((tempF - 32) * 5 / 9, 1)
launch <- -0.6
loglik <- function(b0, b1) sum(dbinom(damage, 6, plogis(b0 + b1 * temp), log = TRUE))
mle <- coef(glm(cbind(damage, 6 - damage) ~ temp, family = binomial))
g0 <- seq(-2, 12, length.out = 90)
g1 <- seq(-0.8, 0.2, length.out = 90)
surface <- outer(g0, g1, Vectorize(loglik))
best <- loglik(mle[[1]], mle[[2]])
ui <- page_sidebar(
sidebar = sidebar(
width = 250,
sliderInput("b0", "β₀ (intercept)", min = -2, max = 12, value = 0, step = 0.05),
sliderInput("b1", "β₁ (effect of temperature)", min = -0.8, max = 0.2, value = 0, step = 0.005),
actionButton("climb", "Let R climb to the top", class = "btn-primary"),
uiOutput("numbers")
),
layout_columns(
col_widths = c(7, 5),
card(card_header("Share of O-rings damaged, by launch temperature"),
plotOutput("fit", height = "420px")),
card(card_header("The likelihood mountain, seen from above"),
plotOutput("map", height = "420px"))
)
)
server <- function(input, output, session) {
observeEvent(input$climb, {
updateSliderInput(session, "b0", value = round(mle[[1]], 2))
updateSliderInput(session, "b1", value = round(mle[[2]], 3))
})
output$fit <- renderPlot({
base_par()
plot(NA, xlim = c(-2, 28), ylim = c(0, 1), xlab = "temperature at launch (°C)",
ylab = "probability an O-ring fails")
rect(-2, 0, min(temp) - 0.3, 1, col = adjustcolor(grey, 0.15), border = NA)
text(4.8, 0.55, "no flights\nthis cold", col = "#6B7591", cex = 0.9)
abline(v = launch, lty = 2, col = ink)
mtext("Challenger: −0.6 °C", side = 3, at = launch, adj = 0, col = ink, cex = 0.8, line = 0.2)
tt <- seq(-2, 28, length.out = 300)
lines(tt, plogis(input$b0 + input$b1 * tt), col = blue, lwd = 3)
points(jitter(temp, 0.3), damage / 6, pch = 21, bg = ink, col = "white", cex = 1.5)
pl <- plogis(input$b0 + input$b1 * launch)
points(launch, pl, pch = 19, col = red, cex = 1.6)
})
output$map <- renderPlot({
base_par()
rel <- pmax(surface - best, -60)
# darker = higher likelihood; the peak is the darkest spot
image(g0, g1, rel, col = hcl.colors(40, "Blues 3", rev = TRUE),
xlab = "β₀", ylab = "β₁")
contour(g0, g1, surface - best, levels = c(-1, -3, -10, -30), add = TRUE,
col = "white", drawlabels = FALSE, lwd = 0.8)
points(mle[[1]], mle[[2]], pch = 3, col = "white", cex = 2.2, lwd = 3)
points(input$b0, input$b1, pch = 19, col = red, cex = 1.6)
})
output$numbers <- renderUI({
l <- loglik(input$b0, input$b1)
pl <- plogis(input$b0 + input$b1 * launch)
tags$div(style = "font-size:.9rem;margin-top:.8rem;",
tags$p(style = "margin:.2rem 0;", "Your log-likelihood: ",
tags$strong(style = paste0("color:", red), sprintf("%.1f", l))),
tags$p(style = "margin:.2rem 0;", sprintf("Best possible: %.2f", best)),
tags$hr(style = "margin:.5rem 0;"),
tags$p(style = "margin:.2rem 0;", "At −0.6 °C your model says each O-ring fails with probability ",
tags$strong(sprintf("%.2f", pl)), sprintf(", so about %.1f of 6 fail.", 6 * pl)))
})
}
shinyApp(ui, server)Try this
- Start at β₀ = 0, β₁ = 0: a flat 50% failure rate at every temperature. The red dot on the map is far from the peak (the darkest region, marked with a white cross). Move the sliders to climb towards it and watch the log-likelihood rise.
- Notice the long diagonal ridge on the map. Many combinations of β₀ and β₁ fit almost equally well, as long as they move together. That ridge is what a strong correlation between two estimates looks like.
- Press Let R climb to the top. The fitted curve predicts that at −0.6 °C nearly every O-ring fails. Look at the grey band, though: the model is extrapolating far beyond the data. The data say “cold is dangerous”, loudly; they can’t say exactly how dangerous at a temperature nobody had tried.
- Set β₁ to 0 and try to find the best β₀. That is the best model without temperature. Compare its log-likelihood with the full model’s: the big gap is the evidence that temperature matters.
The maths, decoded
Every flight contributes one binomial probability, and the log-likelihood adds them up:
\ell(\beta_0,\beta_1) = \sum_{i=1}^{n}\Big[\log\binom{n_i}{y_i} + y_i\log p_i + (n_i-y_i)\log(1-p_i)\Big], \qquad p_i = \frac{1}{1+e^{-(\beta_0+\beta_1\text{temp}_i)}}
Formula decoder
| Symbol | In words |
|---|---|
| n_i = 6 | O-rings on flight i, each assumed to fail independently |
| y_i | how many of them showed damage |
| p_i | failure probability for one O-ring at that flight’s temperature |
| \log\frac{p}{1-p} | the log-odds, which can be any number, so a straight line is fine there |
| \hat\beta_0, \hat\beta_1 | the peak of the mountain; R’s glm(..., family = binomial) finds it for you |
There’s no neat formula for the peak this time. R finds it numerically by walking uphill, the same thing you did with the sliders but faster.
The MLE is a random variable too
The estimate is a function of random data, so, exactly like \hat\beta in the regression lesson, it has a sampling distribution. The slides use waiting times between insurance claims, modelled as exponential with mean \lambda:
f(x\mid\lambda) = \frac{1}{\lambda}e^{-x/\lambda}, \qquad \hat\lambda = \bar X_n, \qquad \text{se}(\hat\lambda) \approx \frac{\lambda}{\sqrt n}.
careful: R’s rexp() wants the rate 1/λ, not the mean λ
The app simulates thousands of datasets from a known \lambda, computes the MLE for each, and compares the spread with what the curvature formula predicts.
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 520
library(shiny)
library(bslib)
ink <- "#1B2559"; red <- "#D6334A"; blue <- "#2463EB"; grey <- "#9AA3B5"
base_par <- function() {
par(mar = c(4, 4.4, 1.6, 1), las = 1, bty = "l", col.axis = ink,
col.lab = ink, fg = grey, cex.axis = 0.85, col.main = ink,
font.main = 1, cex.main = 0.95)
}
ui <- page_sidebar(
sidebar = sidebar(
width = 250,
sliderInput("lam", "True mean waiting time λ (days)", min = 1, max = 10, value = 5, step = 0.5),
sliderInput("n", "Claims per dataset n", min = 5, max = 500, value = 20, step = 5),
actionButton("go", "Simulate again", class = "btn-primary"),
uiOutput("compare")
),
layout_columns(
col_widths = c(6, 6),
card(card_header("15 datasets, 15 log-likelihood curves"),
plotOutput("curves", height = "360px")),
card(card_header("2,000 MLEs vs the curvature formula"),
plotOutput("hist", height = "360px"))
)
)
server <- function(input, output, session) {
sims <- reactive({
input$go
X <- matrix(rexp(2000 * input$n, rate = 1 / input$lam), nrow = 2000)
list(X = X, lhat = rowMeans(X))
})
output$curves <- renderPlot({
s <- sims(); lam <- input$lam; n <- input$n
grid <- seq(lam / 3, lam * 2.5, length.out = 300)
base_par()
plot(NA, xlim = range(grid), ylim = c(-8, 0.3),
xlab = "λ", ylab = "ℓ(λ) − ℓ(λ̂)")
for (i in 1:15) {
xb <- mean(s$X[i, ])
rel <- (-n * log(grid) - n * xb / grid) - (-n * log(xb) - n)
lines(grid, rel, col = adjustcolor(ink, 0.45), lwd = 1.4)
points(xb, 0, pch = 19, col = blue, cex = 0.8)
}
abline(v = lam, col = red, lwd = 2)
mtext("true λ", side = 3, at = lam, col = red, cex = 0.85, line = 0.2)
})
output$hist <- renderPlot({
s <- sims(); lam <- input$lam; n <- input$n
base_par()
h <- hist(s$lhat, breaks = 40, plot = FALSE)
se <- lam / sqrt(n)
xx <- seq(min(s$lhat), max(s$lhat), length.out = 300)
plot(h, freq = FALSE, col = adjustcolor(grey, 0.5), border = "white",
main = "", xlab = "λ̂ (the MLE, i.e. the sample mean)",
ylim = c(0, max(h$density, dnorm(lam, lam, se)) * 1.1))
lines(xx, dnorm(xx, lam, se), col = blue, lwd = 2.5)
abline(v = lam, col = red, lwd = 2)
})
output$compare <- renderUI({
s <- sims(); lam <- input$lam; n <- input$n
tags$div(style = "font-size:.9rem;margin-top:.8rem;",
tags$p(style = "margin:.2rem 0;", "Spread of the 2,000 MLEs: ",
tags$strong(sprintf("%.3f", sd(s$lhat)))),
tags$p(style = "margin:.2rem 0;", "Curvature formula λ/√n: ",
tags$strong(style = paste0("color:", blue), sprintf("%.3f", lam / sqrt(n)))),
tags$p(style = "margin:.2rem 0;color:#4A5578;",
"Blue curve: the normal distribution the theory predicts."))
})
}
shinyApp(ui, server)Try this
- With n = 20, each grey curve peaks in a different place: every dataset has its own MLE. The peaks scatter around the true λ.
- Look at the histogram with n = 5. It’s skewed to the right, and the blue normal curve fits poorly. The normal approximation is an asymptotic result: it needs enough data.
- Increase n to 200. The curves get sharper, the histogram gets narrower and more symmetric, and the simulated spread matches \lambda/\sqrt n closely. You just watched two theorems happen: consistency and asymptotic normality.
What makes the MLE special
In plain words, under standard conditions the MLE is:
- Consistent: with more and more data it homes in on the true value.
- Asymptotically normal: in large samples it is approximately N\big(\theta_0,\, (nJ_\theta)^{-1}\big), which is what the blue curve showed.
- Asymptotically efficient: no unbiased estimator has a smaller variance in large samples (it reaches the Cramér–Rao lower bound).
- Invariant: if \hat\theta is the MLE of \theta, then h(\hat\theta) is the MLE of h(\theta). The MLE of the claim rate is simply 1/\hat\lambda.
These rely on the model being correctly specified, the parameter being identified (a unique peak), and the true value not sitting on the edge of the allowed range. Keep the last one in mind; it comes back in the bootstrap lesson.
OLS is maximum likelihood in disguise
Take the regression model and assume normal errors, y_i = \beta_0 + \beta_1 x_i + \varepsilon_i with \varepsilon_i \sim N(0, \sigma^2). The log-likelihood is
\log L = \underbrace{-\tfrac{n}{2}\log 2\pi}_{\text{constant}} \;\underbrace{-\tfrac{n}{2}\log\sigma^2}_{\text{no } \beta\text{'s here}} \;-\;\frac{1}{2\sigma^2}\sum_{i=1}^n \big(y_i - \beta_0 - \beta_1 x_i\big)^2 .
The only part that involves the \beta’s is minus the sum of squared residuals. Maximising the likelihood is the same as minimising the RSS. The bowl from lesson 1 is this mountain turned upside down. The one small difference is the error variance: the MLE divides by n, while the usual OLS estimate divides by n - p. You will check both in the lab.
Check yourself
Next: the likelihood lab, where you write log-likelihoods yourself, maximise them with optim(), and get standard errors from the curvature.