################################# # R code for LAB4 - 30/10/2026 # # # # ! This is the pre-LAB version # # You must replace TODO properly# ################################# # ____________________________________________________________________________ # Numerical optimisation #### ## ............................................................................ ## A Binomial example #### # Let's generate some data (n = 100, p = 0.6) set.seed(13) n <- 100 p <- 0.6 y <- rbinom(n, 1, p) # ML estimate (closed form): used to check the numerical optimisers mean(y) # logit(theta) psi <- function(theta){ log(theta/(1-theta)) } # theta theta <- function(psi){ exp(psi)/(1 + exp(psi)) } ## Optimisation ---- # Since many optimisers allow only minimisation, # we implement the negative log-likelihood (also considering the reparametrisation) nllik_bin <- function(theta, data){ -sum(data) * log(theta) - (length(data) - sum(data)) * log(1-theta) } nllik_bin_rep <- function(param, data) nllik_bin(theta(param), data) # Score function score <- function(theta, data){ sum(data)/theta - (length(data) - sum(data))/(1 - theta) } # Observed information obs_info <- function(theta, data){ sum(data)/(theta^2) + (length(data) - sum(data))/((1 - theta)^2) } # Let's do a check of the score and the obs info by using numerical derivatives library(numDeriv) grad(nllik_bin, 0.5, data = y) -score(0.5, data = y) hessian(nllik_bin, 0.5, data = y) obs_info(0.5, data = y) # Minimise the negative log-likelihood bin.nlm_start1 <- nlm(f = nllik_bin, p = 0.5, data = y, hessian = TRUE) bin.nlm_start1 nllik_bin(bin.nlm_start1$estimate, data = y) obs_info(bin.nlm_start1$estimate, data = y) # Minimise the negative log-likelihood using optim bin.optim_start1 <- optim(par = 0.5, fn= nllik_bin, hessian = TRUE, data = y, method = "L-BFGS-B", lower = 1e-7, upper = 1-1e-7) bin.optim_start1 bin.optim_start2 <- optim(par = 0.001, fn = nllik_bin, data = y, method = "L-BFGS-B", lower = 1e-7, upper =1-1e-7) bin.optim_start2 optimHess(bin.optim_start2$par, nllik_bin, data = y) # Optimise under reparametrisation bin.nlm_start1_rep <- nlm(f = nllik_bin_rep , p = 0, data = y, hessian = TRUE) bin.nlm_start1_rep nllik_bin_rep (bin.nlm_start1_rep$estimate, data = y) c(bin.nlm_start1_rep$estimate, psi(bin.nlm_start1$estimate), psi(mean(y))) hessian(nllik_bin_rep, bin.nlm_start1_rep$estimate, data = y) # ____________________________________________________________________________ # Likelihood theory: group activity #### ## ===================================================================== ## Group activity - Likelihood theory (30 min) ## Model: waiting times between consecutive orders ~ Exponential, ## parametrized by the MEAN theta (minutes): f(y) = exp(-y/theta)/theta ## (R uses the rate: rexp(n, rate = 1/theta)) ## ===================================================================== ## HOW TO WORK ## - One Pilot types (rotate at every PART), the others read and check. ## - Lines with TODO are yours to complete. ## - Write the answers to the questions in comments: you will report them. ## ## TIMING ## 0-3 INDIVIDUAL predictions (paper, no R) ## 3-11 PART 1 The likelihood, the score, the information ## 11-19 PART 2 Properties of the MLE (simulation) ## 19-27 PART 3 Tests and confidence intervals ## 27-30 Discussion: were your predictions right? ## ===================================================================== ## --------------------------------------------------------------------- ## INDIVIDUAL PREDICTIONS (3 min) - write your answer before touching R. ## Data: n = 8 waiting times (see below), model Exponential with mean theta. ## P1. Is the log-likelihood of theta symmetric around its maximum? ## P2. Will the Wald interval (est +/- 1.96 SE) and the likelihood-ratio ## interval be the same? If not, which one reaches higher values? ## P3. For the RATE lambda = 1/theta, is the MLE 1/ybar unbiased? ## --------------------------------------------------------------------- y <- c(2.1, 9.4, 0.8, 5.6, 14.2, 3.3, 7.9, 1.7) # minutes between orders n <- length(y); S <- sum(y) ## ===================================================================== ## PART 1 - Likelihood, score, information (8 min) ## ===================================================================== ## Log-likelihood, up to an additive constant (as in the slides): ll_exp <- function(theta, n, S) -n * log(theta) - S / theta ## 1a. Plot the log-likelihood minus its maximum, over theta in (2, 20). ## Mark the MLE. th <- seq(2, 20, length.out = 400) lln <- sapply(th, ll_exp, n = n, S = S) plot(th, lln - max(lln), type = "l", lwd = 2, xlab = expression(theta), ylab = expression(l(theta) - max(l))) ## add a vertical line at the MLE (find it from the plot, then check in 1b) TODO ## 1b. MLE: numerical maximisation vs the closed form opt <- optimize(ll_exp, interval = c(0.1, 100), n = n, S = S, maximum = TRUE) opt$maximum TODO # closed form (the sample mean) ## 1c. Likelihood ratio: how much more supported is theta = 5 than theta = 10? TODO # exp(ll_exp(5, n, S) - ll_exp(10, n, S)) ## 1d. Another sample with the SAME n and the SAME sum, different values: y2 <- c(5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 10.0) c(length(y2), sum(y2)) ## Compute ll_exp for the two samples at theta = 3, 5, 8. What do you notice? ## 1e. Score and observed information. Derive them by hand from ll_exp ## (differentiate once and twice with respect to theta) and code them: U <- function(theta, n, S) TODO J <- function(theta, n, S) TODO # J = - second derivative ## 1f. At the MLE: the score and the standard error th_hat <- mean(y) U(th_hat, n, S) se_obs <- TODO # 1 / sqrt(J(th_hat, n, S)) c(se_obs, th_hat / sqrt(n)) ## QUESTIONS ## Q1. How many times is theta = 5 more supported than theta = 10? Is it a lot? ## Q2. What does 1d tell you about the data needed to build the likelihood? ## Link it to sufficient statistics (factorization in the slides). ## Q3. Why is the score 0 at the MLE? Is the log-likelihood symmetric (P1)? ## ===================================================================== ## PART 2 - Properties of the MLE by simulation (8 min) ## True theta = 5, samples of size n = 8, R = 10000 replications. ## ===================================================================== set.seed(2026) R <- 10000; theta_t <- 5 Y <- matrix(rexp(n * R, rate = 1 / theta_t), nrow = R) # each row = one sample Ssim <- rowSums(Y) ## 2a. Score and information EVALUATED AT THE TRUE VALUE, in each sample U_t <- U(theta_t, n, Ssim) J_t <- J(theta_t, n, Ssim) c(mean(U_t), var(U_t)) # expected score; variance of the score c(mean(J_t), n / theta_t^2) # mean observed information vs n/theta^2 ## 2b. Cramer-Rao. Two UNBIASED estimators of theta: ## the MLE (sample mean) and n * min(y_i). est_mle <- rowMeans(Y) est_min <- n * apply(Y, 1, min) rbind(mean = c(mean(est_mle), mean(est_min)), var = c(var(est_mle), var(est_min))) TODO # compute the Cramer-Rao bound theta^2 / n ## 2c. Invariance of the MLE: estimate the RATE lambda = 1/theta ## by plugging in: 1/ybar. Is it unbiased? (true lambda = 0.2) mean(1 / est_mle) mean((n - 1) / Ssim) # another estimator of lambda: what is it doing? ## 2d. (if time) Sampling distribution of the MLE for increasing n par(mfrow = c(1, 3)) for (nn in c(5, 50, 500)) { m <- replicate(5000, mean(rexp(nn, rate = 1 / theta_t))) hist(m, prob = TRUE, breaks = 40, main = paste("n =", nn), xlab = "MLE of theta") curve(dnorm(x, theta_t, theta_t / sqrt(nn)), add = TRUE, col = 2, lwd = 2) } par(mfrow = c(1, 1)) ## QUESTIONS ## Q4. What is E(U) at the true value? What is var(U)? Which identity of the ## slides is this? Same question for mean(J_t). ## Q5. Which of the two unbiased estimators has smaller variance? Does the MLE ## reach the Cramer-Rao bound? Does the other estimator violate it? ## Q6. Is 1/ybar unbiased for lambda? Why can the MLE of a function (plug-in ## rule) be biased even if the MLE of theta is not? ## Q7. In 2d, what changes when n grows? Which property of the MLE is this? ## ===================================================================== ## PART 3 - Tests and confidence intervals (8 min) ## Back to the observed data: n = 8, th_hat = ybar. ## ===================================================================== ## 3a. Four 95% confidence intervals for theta q95 <- qchisq(0.95, 1) W <- function(theta) 2 * (ll_exp(th_hat, n, S) - ll_exp(theta, n, S)) # LRT statistic ci_wald <- th_hat + c(-1, 1) * 1.96 * th_hat / sqrt(n) ci_wlog <- TODO # Wald on log(theta), then exp(): th_hat * exp(c(-1, 1) * 1.96 / sqrt(n)) ci_lrt <- c(uniroot(function(t) W(t) - q95, c(0.5, th_hat))$root, uniroot(function(t) W(t) - q95, c(th_hat, 100))$root) ci_exact <- 2 * S / qchisq(c(0.975, 0.025), df = 2 * n) # exact result rbind(wald = ci_wald, wald_log = ci_wlog, lrt = ci_lrt, exact = ci_exact) ## Draw the LRT interval on the plot of 1a: horizontal line at -q95/2 plot(th, lln - max(lln), type = "l", lwd = 2, xlab = expression(theta), ylab = expression(l(theta) - max(l))) abline(h = -q95 / 2, col = 2) abline(v = ci_lrt, col = 2, lty = 2) abline(v = ci_wald, col = 4, lty = 2) ## 3b. Three tests of H0: theta = 10 (minutes) vs H1: theta != 10 theta0 <- 10 W_lrt <- W(theta0) W_wald <- n * (th_hat - theta0)^2 / th_hat^2 # (th_hat - theta0)^2 * J(th_hat) W_score <- U(theta0, n, S)^2 / (n / theta0^2) # U(theta0)^2 / I(theta0), with I = n / theta^2 stat <- c(LRT = W_lrt, Wald = W_wald, Score = W_score) rbind(stat = stat, p_value = pchisq(stat, df = 1, lower.tail = FALSE)) ## 3c. Coverage in repeated sampling (n = 8, true theta = 5, 1e5 samples). ## A CI covers theta_t exactly when the corresponding test does NOT reject it. set.seed(1) Rr <- 1e5 Sx <- rgamma(Rr, shape = n, scale = theta_t) # sum of n exponentials hx <- Sx / n # MLEs cov_wald <- abs(hx - theta_t) <= 1.96 * hx / sqrt(n) cov_wlog <- TODO # abs(log(hx) - log(theta_t)) <= 1.96 / sqrt(n) cov_lrt <- 2 * n * (log(theta_t / hx) + hx / theta_t - 1) <= q95 # LRT statistic at theta_t <= q95 cov_exact <- (2 * Sx / theta_t >= qchisq(0.025, 2 * n)) & (2 * Sx / theta_t <= qchisq(0.975, 2 * n)) c(wald = mean(cov_wald), wald_log = mean(cov_wlog), lrt = mean(cov_lrt), exact = mean(cov_exact)) ## QUESTIONS ## Q8. Which interval is symmetric? Which are not? Does it match P1-P2? ## Q9. In 3b, do the three tests lead to the same decision at level 5%? ## Which interval agrees with the Wald test, which with the LRT? ## (Is theta0 = 10 inside each interval?) ## Q10. In 3c, which interval is closest to 95%? Why does the Wald interval ## undercover? Repeat with n = 30 and n = 100: what happens? ## Q11. The LRT interval for the rate lambda = 1/theta is just 1/(its interval ## for theta); this is NOT true for the Wald interval. What property of the ## LRT is this? Why do the LRT and Wald-on-log nearly coincide here? ## ===================================================================== ## DISCUSSION (3 min) - Reporter of each team: ## - Which of your three predictions was wrong, and what did it teach you? ## - If you had only one interval to report for n = 8, which one and why? ## =====================================================================