Ten problems on maximum likelihood and Bayesian inference

Post-release survival from a tagging study

You released 80 fish after barotrauma treatment and tracked them with acoustic tags; 47 were detected alive after 7 days. Treat each fish's outcome as an independent Bernoulli($\phi$), where $\phi$ is the probability of post-release survival.

1. Write the joint likelihood $L(\phi)$ for the 80 fish. Take the log, differentiate, and solve algebraically for the MLE $\hat\phi$.
2. Two confidence intervals, compute both and compare. Neither the term "Fisher information" nor the "Wald interval" was covered in lecture, so here is what they mean.
   • The observed Fisher information is the curvature of the log-likelihood at the MLE: the negative second derivative $I(\hat\phi) = -\ell''(\hat\phi)$ (for a model with several parameters it is the negative Hessian matrix). A sharper peak means more information and a smaller standard error, $\text{SE} = 1/\sqrt{I(\hat\phi)}$. The Wald 95% CI is then the symmetric interval $\hat\phi \pm 1.96\,\text{SE}$. Compute $I(\hat\phi)$ and this Wald interval. (For a fitted mle()/optim() object, R hands you the same SE via sqrt(diag(vcov(fit))), which inverts the Hessian for you.)
   • The profile-likelihood CI, the kind this course teaches, is the set of $\phi$ whose log-likelihood is within 1.92 units of the maximum, where $1.92 = \tfrac12\chi^2_{1,\,0.95} =$ 0.5 * qchisq(0.95, 1). Obtain it with stats4::mle() and confint().
   Report both intervals side by side and explain any asymmetry, that is, which one can be lopsided about $\hat\phi$, and why.
3. The manager asks: "What is the probability that true survival is at least 0.5?" Explain why maximum likelihood alone cannot answer that question directly, and name the framework (taught this afternoon) that can.

Von Bertalanffy growth from length-at-age data (vonbert.csv)

The file data/vonbert.csv holds length-at-age records with three columns, sex (1 or 2), age (years), and length. The Von Bertalanffy growth function (VBGF) is

where $L_\infty$ is asymptotic length and $\kappa$ is the growth coefficient. Load the data exactly as in the lab:

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

1. Define your own R function, named nll_vbgf(Linf, kappa, sigma), that takes the three parameters as arguments and returns the negative log-likelihood of the pooled data (both sexes). (nll_vbgf is not built in; you write it, exactly as in the lab walk-through.) Holding $L_\infty = 100$ and $\sigma = 10$ fixed, run a grid search for $\kappa$ over seq(0.05, 0.50, by = 0.005). Plot the NLL against $\kappa$ and report the grid-optimal $\hat\kappa$.
2. Now free all three parameters. Fit the full model with mle() (method "L-BFGS-B" with sensible lower bounds). Report $\hat L_\infty$, $\hat\kappa$, $\hat\sigma$, and the profile-likelihood 95% CIs from confint().
3. Test whether the sexes share a growth curve with a likelihood-ratio test. Fit a null model (one $L_\infty$, one $\kappa$, shared $\sigma$) and an alternative with sex-specific $L_\infty$ and $\kappa$ (shared $\sigma$). Compute $\Lambda = 2(\ell_{\text{alt}} - \ell_{\text{null}})$, compare to $\chi^2_2$, and report the p-value. Cross-check with AIC().
4. Explain why grid search does not scale to the full three-parameter model, and why it remains a useful sanity check on the optimizer.

Two-component Poisson mixture for whale blow counts (whale.csv)

The file data/whale.csv holds blow counts for individual whales (column sex codes each whale 1 or 2; column blows is the count). The population is a mix of males and females that blow at different average rates. Treat each whale's sex as unobserved and model the counts as a two-component mixture, where males blow at rate $\lambda_1$, females at rate $\lambda_2$, and a whale is male with probability $\pi$ (the proportion of males) and female with probability $1-\pi$. Your goal is to estimate the male rate $\lambda_1$, the female rate $\lambda_2$, and the proportions $\pi$ and $1-\pi$:

TheData <- read.csv("data/whale.csv")
blows <- TheData$blows   # the mixture ignores the sex column; fit is unsupervised

1. Explain, in one or two sentences, why the mixture MLE is not just the pair of single-component Poisson MLEs. Identify precisely which step in the usual derivation ($\log\prod = \sum\log$) breaks down.
2. Write the joint negative log-likelihood and fit it with mle(), starting from $\lambda_1 < \lambda_2$ (e.g. lam1 = 10, lam2 = 30, p = 0.4). Report the estimated male rate $\hat\lambda_1$, female rate $\hat\lambda_2$, proportion of males $\hat\pi$, and their standard errors.
3. Refit from at least 8 different starting values. How many distinct local optima appear? Report the global best (lowest NLL) and the fraction of starts that reached it.
4. The two-component mixture is invariant under relabeling: $(\lambda_1,\lambda_2,\pi)$ and $(\lambda_2,\lambda_1,1-\pi)$ are the same model. Explain how this shows up if you don't constrain $\lambda_1 < \lambda_2$, and why it would corrupt your standard errors.
5. Validate against the labels. Your mixture never used the sex column. It estimated $\hat\pi$, the proportion of whales in the low-rate component, from the counts alone. In this dataset sex is coded 1 = male, 2 = female, and males are the lower-rate group. Compute the actual proportion of males, mean(TheData$sex == 1), and compare it to $\hat\pi$. Also compare the observed mean blow count within each sex to $\hat\lambda_1$ and $\hat\lambda_2$. Does the unsupervised mixture recover the real sex structure?

Logistic (r/K) growth by maximum likelihood (Logistic.csv)

The file data/Logistic.csv has two columns: PopulationSize ($N$) and GrowthRate (the per-capita rate $\tfrac{1}{N}\tfrac{dN}{dt}$). The per-capita form of logistic growth is linear in $N$:

Logistic <- read.csv("data/Logistic.csv", header = TRUE)
Logistic <- Logistic[order(Logistic$PopulationSize), ]

1. The model is linear, with intercept $= r$ and slope $= -r/K$. Fit it with lm(GrowthRate ~ PopulationSize). Because a Normal-error linear model fit by least squares is its maximum-likelihood fit, the fitted intercept and slope are already the MLEs of $r$ and $-r/K$. Now use the invariance property of the MLE: if $\hat\theta$ is the MLE of $\theta$, then $g(\hat\theta)$ is the MLE of any transformation $g(\theta)$. That is exactly why you may push the fitted coefficients through a formula to get MLEs of the biological parameters, e.g. $\hat K = -\hat r / \widehat{\text{slope}}$ is the MLE of $K$. Recover $\hat r$ and $\hat K$ this way. (You will invoke the same invariance property again in part (c) to turn $\hat r,\hat K$ into the MLE of MSY.)
2. Now fit the same model directly by maximum likelihood with mle(), parameterized in $(r, K, \sigma)$ so that $\mu_i = r - (r/K)N_i$. Confirm you recover the same $\hat r$ and $\hat K$ as in (a), and report profile-likelihood 95% CIs. Start from this skeleton:
   R · direct MLE in (r, K, sigma)

   library(stats4)
   x <- Logistic$PopulationSize; y <- Logistic$GrowthRate
   nll_rk <- function(r, K, sigma) {
    mu <- r - (r / K) * x
    -sum(dnorm(y, mu, sigma, log = TRUE))
   }
   fit <- mle(nll_rk, start = list(r = 0.2, K = 1000, sigma = 0.03),
    method = "L-BFGS-B", lower = c(0.001, 100, 1e-4),
    upper = c(2, 4000, 2),
    control = list(parscale = c(0.2, 1000, 0.03)))
   summary(fit); confint(fit)

   Keep the control = list(parscale = ...) line. optim(), which mle() calls, takes one common step size in all three coordinates, and here the three parameters differ in scale by four orders of magnitude ($r \approx 0.2$, $K \approx 10^3$, $\sigma \approx 0.03$). A step large enough to move $K$ is catastrophic for $r$ and $\sigma$, so the optimizer converges on the first two and leaves $K$ sitting at its starting value. Drop the control argument and you get $\hat K = 1000.00$, exactly where you started, with no warning of any kind. parscale divides each parameter by its own typical magnitude so that one step size fits all three. If your $\hat K$ ever comes back equal to your starting value, suspect scaling first.
3. Using the invariance of the MLE, compute the plug-in MLE of the maximum sustainable yield $\text{MSY} = \hat r\hat K/4$ and of the population size at maximum growth $K/2$.
4. Maximum likelihood gives you point estimates and CIs for $r$ and $K$, but propagating that uncertainty into a derived quantity like MSY is awkward (delta method, bootstrap). Note in one sentence why. You will revisit this exact model in Problem 7, where the Bayesian posterior propagates the uncertainty for free.

Bayesian survival with an informative prior

You marked 50 adult animals in spring; 38 were re-encountered alive the following year. Literature suggests adult survival is around 0.55 with moderate uncertainty, which you encode as a Beta(11, 9) prior on $\phi$ (annual survival).

1. Justify the Beta(11, 9) prior. What are its prior mean and its effective (pseudo) sample size, and what do they represent?
2. Write the model in Nimble, with $y \sim \text{Binomial}(n = 50, \phi)$ and the Beta(11, 9) prior, and fit it with nimbleMCMC() using 3 chains.
3. Report the posterior mean, median, 95% BCI, and 95% HPDI for $\phi$, plus $\hat R$ and effective sample size.
4. Compute $\Pr(\phi > 0.7 \mid \text{data})$ from the MCMC draws. Then, because Beta–Binomial is conjugate, compare to the closed-form Beta$(11+38,\, 9+12)$ = Beta(49, 21) posterior to verify Nimble is doing what you think.

Bayesian detection rate for a species of concern

You ran 25 ten-minute point counts and recorded these per-count detection numbers:

y = c(0,1,0,2,0,1,0,0,3,1,0,0,1,2,0,1,0,1,0,0,2,0,1,1,0)

Model the counts as iid Poisson($\lambda$), where $\lambda$ is the mean detections per count.

1. Place a weakly informative Gamma($\alpha = 1$, $\beta = 1$) prior on $\lambda$. Write the model in Nimble and fit it.
2. Report the posterior mean and 95% BCI for $\lambda$.
3. Compute $\Pr(\lambda < 1 \mid y)$, a direct statement about how rare detections are, from the posterior draws.
4. Gamma–Poisson is conjugate, so verify your Nimble fit against the closed-form Gamma$(1 + \sum y_i,\; 1 + n)$ posterior (mean, 95% interval).

Bayesian logistic (r/K) growth with derived quantities (Logistic.csv)

Return to data/Logistic.csv from Problem 4, but now fit the per-capita logistic model in the Bayesian framework so you get full posterior distributions for $r$, $K$, and any derived quantity. The template Logistic_BayesianGrowth_Nimble.R in this folder is your scaffold.

Logistic <- read.csv("data/Logistic.csv", header = TRUE)
Logistic <- Logistic[order(Logistic$PopulationSize), ]

1. Write the model in Nimble with $y_i \sim \mathcal{N}\!\big(r - (r/K)x_i,\; \sigma^2\big)$ and priors $K \sim \text{Uniform}(0, 4000)$, $r \sim \text{Uniform}(0, 2)$, $\sigma \sim \text{Uniform}(0, 2)$ (use BUGS-style precision tau <- 1/sigma^2). Fit with 3 chains.
2. Report posterior means and 95% BCIs for $r$, $K$, and $\sigma$. Confirm the posterior of $r$ is consistent with the morning MLE from Problem 4.
3. Add the derived quantities $\text{MSY} = rK/4$ and $N_{\max} = K/2$ inside the model block. Report their posterior medians and 95% HPDIs. This is the "for free" uncertainty propagation promised in Problem 4(d).
4. Compute $\Pr(r > 0.22 \mid y)$ from the posterior draws, and plot the posterior of $\sigma$ with both its 95% BCI and 95% HPDI overlaid. Comment on why the two intervals differ for $\sigma$.
5. From a posterior to a decision. A point estimate hands a manager one number; a posterior tells them how much risk sits behind a proposed quota. Compute $\Pr(\text{MSY} < Q \mid y)$, the probability that an annual quota $Q$ is larger than what the population can actually sustain, for $Q = 55,\ 60,\ 62,\ 65$ and $70$. Report the five probabilities. Your point estimate of MSY is close to 62. Explain why handing the manager “62” and handing them these five numbers can lead to different decisions, and say which quota you would defend and on what grounds.
6. Is $K \sim \text{Uniform}(0, 4000)$ actually uninformative? A vague-looking prior is only vague if the data pin the parameter down, which is a claim you can test rather than assume. Do it in two steps, changing nothing but the bound (and keeping your starting values inside it, or NIMBLE will refuse to initialise):
   1. Refit the full dataset three times, with $K \sim U(0, 2000)$, $U(0, 4000)$ and $U(0, 20000)$. Report the posterior median and 95% BCI of $K$ and of MSY each time. How far does the bound move the answer?
   2. Now refit using only the first eight rows after the sort, head(Logistic, 8), the early low-density part of the curve where the population is nowhere near its ceiling, again under all three bounds. Report the same quantities.
   Something very different happens in the two cases. In two or three sentences, say what it is, what has to be true of your data before a bounded uniform prior deserves to be called “uninformative”, and what all this implies about the MSY you would have reported if the short time series were the only data you had.

Prior sensitivity in a data-poor problem

You observed 3 detections in 20 surveys. Compare four priors on the detection probability $p$: Beta(1, 1) (flat), Beta(1, 9) (skeptical, prior mean 0.10), Beta(20, 30) (informative around 0.40), and Beta(100, 100) (very informative around 0.50).

1. For each prior, write down the posterior in closed form (Beta–Binomial conjugacy) and its posterior mean.
2. Plot all four posteriors on the same axes. Which posterior is closest to the data MLE $\hat p = 3/20 = 0.15$, and which is most shifted away from it?
3. A published paper fits this same dataset and reports a posterior mean of 0.49. Which prior most likely produced that result, and what specific concern would you raise as a reviewer?

BCI vs HPDI for a skewed posterior

You fit a Poisson model to clutch sizes at $n = 12$ nests. Write $y_i$ for the clutch size (number of eggs) at nest $i$; the data enter the Gamma–Poisson posterior only through their total, $\sum_{i=1}^{12} y_i$, the total number of eggs summed across all 12 nests, here $\sum y_i = 38$. With a Gamma(2, 1) prior, the posterior on the mean clutch size $\lambda$ is Gamma$\bigl(2 + \textstyle\sum y_i,\; 1 + n\bigr)$ = Gamma$(2 + 38,\; 1 + 12)$ = Gamma(40, 13).

1. Compute the posterior mean, mode, median, 95% equal-tailed BCI, and 95% HPDI for $\lambda$ (use HDInterval::hdi()).
2. Plot the posterior density and overlay both intervals.
3. Explain why the BCI and HPDI are nearly identical here, but would separate substantially for a strongly skewed posterior, e.g. a Gamma(1, 2) obtained from a Gamma(1, 1) prior and a single observation $y = 0$. What geometric feature of the density anchors the HPDI in that case?

Identifiability and label switching in a Bayesian mixture

You refit the two-component whale-blow mixture from Problem 3 (males vs females) in Nimble instead of with mle(). With 3 chains, no constraint on the ordering of $\lambda_1$ and $\lambda_2$, and chains started from different labelings, each chain locks onto its own labeling and stays there. One chain reports "$\lambda_1$" $\approx 10$ while another reports "$\lambda_1$" $\approx 30$ for the whole run. Every individual trace plot looks perfectly well mixed, but the chains disagree about what "$\lambda_1$" means, and $\hat R$ is enormous.

The culprit is identifiability. A parameter, or a combination of parameters, is non-identifiable when two different parameter values produce exactly the same likelihood, so no amount of data can separate them. Here the mixture is invariant under relabeling: $(\lambda_1,\lambda_2,\pi)$ and $(\lambda_2,\lambda_1,\,1-\pi)$ give an identical likelihood. Swapping the "male" and "female" labels changes nothing the data can see. That exchange symmetry is what breaks the sampler.

TheData <- read.csv("data/whale.csv") # columns: sex (1/2), blows
blows <- TheData$blows                    # the mixture ignores the sex column

1. Explain geometrically what is happening. Why is $\hat R$ huge even though each chain has, in a sense, "converged" and the model fit is fine?
2. Reproduce the pathology. First fit the mixture with no ordering constraint, deliberately starting chain 2 from the mirrored labeling (rates swapped, labels flipped) so the chains disagree from the outset. Run this and record $\hat R$ for lam[1]. This is your "before" number:
   R · unconstrained mixture, the "before" run

   library(nimble); library(MCMCvis); library(coda)
   set.seed(10)                                # so everyone gets the same Rhat
   unc_code <- nimbleCode({
    lam[1] ~ dunif(0, 50)
    lam[2] ~ dunif(0, 50)     # nothing forces lam[1] < lam[2]
    pi    ~ dbeta(1, 1)
    pz[1] <- pi
    pz[2] <- 1 - pi
    for (i in 1:n) {
    z[i] ~ dcat(pz[1:2])
    y[i] ~ dpois(lam[z[i]])
    }
   })
   zinit <- ifelse(blows < 20, 1, 2)             # rough starting labels
   unc <- nimbleMCMC(unc_code,
    constants = list(n = length(blows)),
    data      = list(y = blows),
    inits     = list(list(lam = c(8, 28),  pi = 0.4,  z = zinit),
                     list(lam = c(28, 8),  pi = 0.6,  z = 3 - zinit),  # mirrored start
                     list(lam = c(10, 32), pi = 0.35, z = zinit)),
    monitors  = c("lam", "pi"),
    nchains = 3, nburnin = 2000, niter = 8000,
    samplesAsCodaMCMC = TRUE)
   MCMCsummary(unc)                            # Rhat for lam[1] is enormous

   Also plot the three traces for lam[1] and note that each one is individually flat and well mixed. The disagreement is between chains, not within any one of them.
3. Fix it with an ordering constraint, then confirm the fix. The clean cure is to make the two rates un-swappable by construction, forcing $\lambda_1 < \lambda_2$. One simple way is to give the lower (male) rate its own prior, build the higher (female) rate on top of it with a positive gap $\delta$, and give each whale a latent sex label. Complete and run this skeleton (3 chains, the same blows vector as Problem 3):
   R · Nimble mixture with λ₁ < λ₂ enforced

   library(nimble); library(MCMCvis); library(coda)
   set.seed(11)                                # so everyone gets the same Rhat
   mix_code <- nimbleCode({
    lam1  ~ dunif(0, 50)        # lower (male) rate
    delta ~ dunif(0, 50)        # positive gap, so lam[2] > lam[1]
    lam[1] <- lam1
    lam[2] <- lam1 + delta       # higher (female) rate
    pi    ~ dbeta(1, 1)        # P(a whale is component 1 = male)
    pz[1] <- pi
    pz[2] <- 1 - pi
    for (i in 1:n) {
    z[i] ~ dcat(pz[1:2])       # latent sex label (1 or 2)
    y[i] ~ dpois(lam[z[i]])
    }
   })
   zinit <- ifelse(blows < 20, 1, 2)             # rough starting labels
   out <- nimbleMCMC(mix_code,
    constants = list(n = length(blows)),
    data      = list(y = blows),
    inits     = list(list(lam1 = 8,  delta = 20, pi = 0.4,  z = zinit),
                     list(lam1 = 12, delta = 18, pi = 0.5,  z = zinit),
                     list(lam1 = 10, delta = 22, pi = 0.35, z = zinit)),
    monitors  = c("lam1", "delta", "pi"),
    nchains = 3, nburnin = 2000, niter = 8000,
    samplesAsCodaMCMC = TRUE)
   MCMCsummary(out)                             # check Rhat ~ 1

   What to look for: after the constraint every chain agrees on which rate is which. $\lambda_1$ sits near 10 and $\lambda_2 = $ lam1 + delta near 30 in all three chains, rather than one chain having them the other way round. $\hat R$ drops to $\approx 1$, and the posterior means land on the Problem 3 MLEs ($\lambda_1\approx 10$, $\lambda_2\approx 30$, $\pi\approx 0.42$). Report $\hat R$ for $\lambda_1$ before and after adding the constraint. The two numbers should be dramatically different.
4. If you relabel the posterior draws post hoc so that $\lambda_1 < \lambda_2$ always, what bias might you introduce, and when is such post hoc relabeling justified?