Worked solutions — Maximum Likelihood & Bayesian I

Exercise 1 — Bernoulli MLE

$y_i \sim \text{Bernoulli}(p)$ independently, $s = \sum y_i$, $n$ trials.

(a) Likelihood. Since the $y_i$ are independent,

(b) MLE. Log-likelihood: $\ell(p) = s\log p + (n-s)\log(1-p)$. Then

(c) Fisher information. $\ell''(p) = -s/p^2 - (n-s)/(1-p)^2$. At $\hat p = s/n$,

So Fisher information is $I(\hat p) = n/[\hat p(1-\hat p)]$, and Wald SE is $1/\sqrt{I(\hat p)} = \sqrt{\hat p(1-\hat p)/n}$.

(d) Red-light data. $s=30, n=50$. $\hat p = 0.6$; Wald SE $= \sqrt{0.6 \cdot 0.4/50} \approx 0.0693$; Wald 95% CI $\approx (0.464, 0.736)$; profile 95% CI $\approx (0.462, 0.728)$ from confint(mle(...)).

library(stats4)
nll <- function(p) -(30 * log(p) + 20 * log(1 - p))
fit <- mle(nll, start = list(p = 0.5), method = "L-BFGS-B",
 lower = 0.001, upper = 0.999)
summary(fit)
confint(fit) # profile-likelihood CI: (0.462, 0.728)

Exercise 2 — Profile likelihood for κ

The profile-likelihood for $\kappa$ is built by holding $\kappa$ fixed at each grid point and optimizing the other two parameters; the resulting curve has its minimum at the joint MLE.

library(stats4)
TheData <- read.csv("data/vonbert.csv")
Age <- TheData$age; Len <- TheData$length

profile_kappa <- function(k) {
 nll <- function(Linf, sigma) {
 pred <- Linf * (1 - exp(-k * Age))
 -sum(dnorm(Len, pred, sigma, log = TRUE))
 }
 fit <- mle(nll, start = list(Linf = 100, sigma = 10),
 method = "L-BFGS-B", lower = c(10, 0.01))
 -logLik(fit)[1]
}
grid <- seq(0.05, 0.45, by = 0.005)
prof <- sapply(grid, profile_kappa)
best <- min(prof)

plot(grid, prof - best, type = "l", lwd = 2,
 xlab = expression(kappa), ylab = expression(-Delta * log~L))
abline(h = 1.92, lty = 2, col = "red")
abline(v = grid[which.min(prof)], col = "grey")

ci <- range(grid[(prof - best) <= 1.92])
cat("Profile 95% CI for kappa:", round(ci, 3), "\n")

On the shipped data: $\hat\kappa \approx 0.24$ (pooled across sexes), profile 95% CI roughly $(0.205, 0.285)$. The curve is clearly asymmetric — Wald would have given a symmetric (and slightly wrong) CI.

Exercise 3 — Sex-specific VBGF revisited

(a) From the lab walk-through code (Model A = shared, Model B = full sex-specific):

LRT <- 2 * (-logLik(fitA)[1] - (-logLik(fitB)[1]))
pval <- pchisq(LRT, df = 2, lower.tail = FALSE)
c(LRT = LRT, p_value = pval)
AIC(fitA, fitB)

On the shipped data: $\Lambda \approx 75.4$ on 2 df, $p \approx 4\times10^{-17}$ — overwhelming. Both LRT and AIC decisively favor Model B (AIC $\approx 798$ for Model A vs $\approx 727$ for Model B, a drop of ~71).

(b) Intermediate model (shared $L_\infty$, sex-specific $\kappa$, 4 parameters): AIC $\approx 767$ — better than the shared Model A ($\approx 798$) but still well above the full Model B ($\approx 727$). Each sex-specific parameter helps, and the fully sex-specific model wins. Fitting Model B shows the difference is concentrated in asymptotic length $L_\infty$ ($\hat L_{\infty,1}\approx 102.5$ vs $\hat L_{\infty,2}\approx 82.6$) more than in growth rate $\kappa$ ($0.284$ vs $0.259$).

(c) Conclusion: the sexes differ mainly in asymptotic size $L_\infty$ (one sex reaches a substantially larger maximum length) and only modestly in $\kappa$; the full sex-specific model is preferred by both criteria. LRT requires nesting and delivers a p-value; AIC ranks (no nesting required) and here also points to Model B.

Exercise 4 — Zero-inflated Poisson

library(stats4)
fit_zip_and_pois <- function(lambda_true) {
 set.seed(99)
 n <- 300; pi_true <- 0.35
 z <- rbinom(n, 1, pi_true)
 y <- ifelse(z == 1, 0, rpois(n, lambda_true))

 nll_pois <- function(lam) -sum(dpois(y, lam, log = TRUE))
 nll_zip <- function(p, lam) {
 if (p <= 0 || p >= 1 || lam <= 0) return(1e10)
 pz <- p + (1 - p) * exp(-lam)
 ll <- ifelse(y == 0, log(pz), log(1 - p) + dpois(y, lam, log = TRUE))
 -sum(ll)
 }
 fp <- mle(nll_pois, start = list(lam = mean(y)), method = "L-BFGS-B", lower = 0.001)
 fz <- mle(nll_zip, start = list(p = 0.2, lam = 1), method = "L-BFGS-B",
 lower = c(0.001, 0.001), upper = c(0.999, 50))
 data.frame(lambda_true = lambda_true,
 obs_p0 = mean(y == 0),
 pois_p0 = exp(-coef(fp)),
 zip_p0 = coef(fz)["p"] + (1 - coef(fz)["p"]) * exp(-coef(fz)["lam"]),
 AIC_pois = AIC(fp),
 AIC_zip = AIC(fz))
}
do.call(rbind, lapply(c(1, 2.5, 8), fit_zip_and_pois))

Expected pattern:

(c) Why the small-$\lambda$ case is hardest: at $\lambda_{\text{true}}=1$ the extra zeros drag the plain Poisson's estimate down to $\hat\lambda \approx 0.62$, and that fitted Poisson already produces P(Y=0) $= e^{-0.62} \approx 0.54$ — not far from the inflated observed proportion of ~0.63. The plain Poisson absorbs most of the zero inflation just by shrinking $\lambda$. The likelihood signal that distinguishes Poisson from ZIP is weak. At $\lambda = 8$, plain Poisson predicts essentially zero zeros; any zero inflation is unmistakable.

Exercise 5 — Beta–Binomial posterior arithmetic

(a) Prior Beta(3, 7), data $y=22$, $n=40$. Posterior is Beta(3+22, 7+(40−22)) = Beta(25, 25).

library(HDInterval)
post <- rbeta(2e5, 25, 25)
c(mean = mean(post),
 median = median(post),
 BCI = quantile(post, c(.025, .975)),
 HPDI = hdi(post, .95))

(b) Posterior mean = 0.5, median = 0.5, BCI ≈ (0.36, 0.64), HPDI ≈ (0.36, 0.64). Symmetric posterior, so BCI ≈ HPDI.

(c) Beta(1, 1) prior + same data → Beta(23, 19), mean $\approx 0.548$, BCI ≈ (0.40, 0.70). The flat prior makes the posterior mean equal the sample mean shifted up by one in numerator and denominator. The skeptical prior pulled the posterior down toward 0.30; the flat prior didn't pull at all.

Exercise 6 — Metropolis for Gamma–Poisson

(a) Working in $\log\lambda$ ensures (i) the proposal is symmetric in the unconstrained space, (ii) you can't propose a negative $\lambda$, and (iii) tail behavior is well-handled. The Jacobian of the transformation contributes a $+\log\lambda$ term to the log-density (because $d\lambda/d(\log\lambda) = \lambda$).

(b–c) Code:

y <- c(3,1,4,2,0,5,2,3,1,4)
log_post <- function(log_lam) {
 lam <- exp(log_lam)
 dgamma(lam, 2, 1, log = TRUE) + sum(dpois(y, lam, log = TRUE)) + log_lam
}
set.seed(3)
ll <- numeric(3e4); ll[1] <- log(2); acc <- 0; s <- 0.30
for (i in 2:length(ll)) {
 prop <- rnorm(1, ll[i - 1], s)
 if (log(runif(1)) < log_post(prop) - log_post(ll[i - 1])) {
 ll[i] <- prop; acc <- acc + 1
 } else ll[i] <- ll[i - 1]
}
cat("acceptance =", round(acc / length(ll), 3), "\n")

lam <- exp(ll[-(1:3000)])
hist(lam, breaks = 60, freq = FALSE, main = "", xlab = expression(lambda))
curve(dgamma(x, 2 + sum(y), 1 + length(y)), add = TRUE, lwd = 2, col = "red")

Acceptance with $s = 0.30$ is about 58% here — above the ~30% target, so to tune down you would raise $s$ (bigger jumps get rejected more often). The posterior approximation matches the closed-form Gamma(27, 11) (mean 2.45) to within MCMC noise.

Exercise 7 — Logistic growth in Nimble

The full code is in Logistic_BayesianGrowth_Nimble.R. Typical posterior summaries (rounded; will vary by seed):

(c) $\Pr(r > 0.22 \mid \text{data}) \approx 0.025$ — the $r$ posterior is tight (sd $\approx 0.01$), so 0.22 sits about two SDs above the mean.

(d) $\Pr(\text{MSY} > 50 \mid \text{data}) \approx 1.0$ — essentially certain, since the entire 95% MSY interval (~57 to 67) lies above 50.

(e) The clearest BCI–HPDI separation is on $\sigma$, which is bounded below by zero and mildly right-skewed: the HPDI ($\approx 0.023$–$0.035$) sits a hair lower than the equal-tailed BCI and is marginally shorter. Why: equal-tail intervals split the 5% miscoverage symmetrically by quantile; HPDI captures the densest 95% of mass, which for a right-skewed distribution shifts the interval toward (but not past) the lower bound. The effect is modest here because $\sigma$ is well-determined — it would be dramatic for a strongly skewed posterior.

Exercise 8 — Bayesian bluebird survival

Prior pseudo-sample sizes: Beta(1,1) = 2, Beta(9,11) = 20, Beta(50,50) = 100. Data: $y=18, n=25$.

library(nimble); library(MCMCvis)
fit_one <- function(a, b) {
 code <- nimbleCode({ phi ~ dbeta(a, b); y ~ dbinom(phi, n) })
 out <- nimbleMCMC(code = code,
 constants = list(n = 25, a = a, b = b),
 data = list(y = 18),
 inits = list(phi = 0.5),
 monitors = "phi",
 nchains = 1, nburnin = 2000, niter = 12000,
 samplesAsCodaMCMC = TRUE)
 list(mean = mean(MCMCchains(out, "phi")),
 BCI = quantile(MCMCchains(out, "phi"), c(.025, .975)),
 p05 = mean(MCMCchains(out, "phi") > 0.5))
}
lapply(list(c(1,1), c(9,11), c(50,50)), function(p) fit_one(p[1], p[2]))

(c) The data have $n = 25$ trials. A Beta($\alpha,\beta$) prior is worth $\alpha+\beta$ trials' worth of information. So:

• Beta(1,1) is essentially weightless (2 vs 25) — posterior mirrors the MLE 18/25 = 0.72.
• Beta(9,11) is worth 20 trials, comparable to the data; the posterior splits the difference between 0.45 (prior mean) and 0.72.
• Beta(50,50) is worth 100 trials, 4× the data; the posterior is yanked strongly toward 0.50.