Likelihood lab

Write a likelihood, let the computer climb it

The recipe for every exercise on this page is the same three steps, so it’s worth memorising:

  1. Write down the log-likelihood as an R function of the parameters.
  2. Hand it to an optimiser: optimize() for one parameter, optim() for several.
  3. Read off the peak (the estimate) and the curvature there (the standard error).

One R quirk to know. optim() minimises by default. So we give it the negative log-likelihood: the lowest point of -\ell is the highest point of \ell.

Exercise 1: the coin, numerically

97 heads in 500 tosses. Write the log-likelihood as a function of p and let optimize() find the top. dbinom(k, n, p, log = TRUE) gives the log of the binomial probability.

You want the maximum, so the argument should be TRUE.

loglik <- function(p) dbinom(97, 500, p, log = TRUE)
fit <- optimize(loglik, interval = c(0.001, 0.999), maximum = TRUE)
fit$maximum

Exercise 2: waiting times, and a trap

We have 100 waiting times between insurance claims. The lesson’s model is exponential with mean \lambda:

\ell(\lambda) = -n\log\lambda - \frac{\sum_i x_i}{\lambda}.

Write the negative of this and minimise it. The theory says the answer should equal mean(claims).

R’s rexp() and dexp() use the rate. rate = 1 / mean.

Flip the sign of each term: n * log(lambda) + sum(claims) / lambda. (Equivalently, -sum(dexp(claims, rate = 1 / lambda, log = TRUE)).)

n <- length(claims)
negloglik <- function(lambda) n * log(lambda) + sum(claims) / lambda
fit <- optimize(negloglik, interval = c(0.01, 50))
c(mle = fit$minimum, sample_mean = mean(claims))

Exercise 3: a standard error from the curvature

Ask optim() for the Hessian (the curvature at the minimum of -\ell) and turn it into a standard error with \text{se} = 1/\sqrt{H}. Theory says it should be about \hat\lambda/\sqrt{n}.

1 / sqrt(H)

fit <- optim(par = 1, fn = negloglik, method = "Brent",
             lower = 0.01, upper = 50, hessian = TRUE)
H  <- fit$hessian[1, 1]
se <- 1 / sqrt(H)
c(se_from_curvature = se, theory = fit$par / sqrt(n))

With several parameters the Hessian is a matrix and the standard errors are sqrt(diag(solve(H))). You’ll need that next.

Exercise 4: the Challenger model by hand

The 23 flights are loaded as temp (°C) and damage (O-rings damaged, out of 6). Write the negative log-likelihood of the logistic model, maximise it with optim(), and compare with glm(). plogis() is R’s logistic function, 1/(1 + e^{-z}).

You already computed the per-flight probabilities and called them p.

negloglik <- function(beta) {
  p <- plogis(beta[1] + beta[2] * temp)
  -sum(dbinom(damage, size = 6, prob = p, log = TRUE))
}
fit <- optim(c(0, 0), negloglik, method = "BFGS", hessian = TRUE)
by_hand <- fit$par
by_glm  <- coef(glm(cbind(damage, 6 - damage) ~ temp, family = binomial))
rbind(by_hand, by_glm)

Now the standard errors, from the curvature, next to glm()’s. Run it:

−0.98: that’s the long diagonal ridge on the likelihood map

Exercise 5: OLS is maximum likelihood

Fit a regression by maximum likelihood with normal errors and compare with lm(). We estimate \log\sigma rather than \sigma so the optimiser can’t wander into negative variances.

Under the model, y_i \sim N(\beta_0 + \beta_1 x_i, \sigma^2). So the mean is b0 + b1 * x.

negloglik <- function(theta) {
  b0 <- theta[1]; b1 <- theta[2]; sigma <- exp(theta[3])
  -sum(dnorm(y, mean = b0 + b1 * x, sd = sigma, log = TRUE))
}
fit <- optim(c(0, 0, 0), negloglik, method = "BFGS")
rbind(mle = fit$par[1:2], ols = coef(lm(y ~ x)))

The coefficients match, but the variance estimates differ slightly. Run this to see why:

The MLE divides by n and is slightly biased downward in small samples; lm() divides by n - p. With 80 observations the difference is small; with 10 it wouldn’t be.

Free play

Some ideas: refit the Challenger model with a probit link (pnorm instead of plogis) and compare predictions at −0.6 °C; or simulate exponential data with n = 5 and see how far \hat\lambda/\sqrt n is from the true spread.

Next lesson: The bootstrap.