Maximum likelihood in the morning, Bayes in the afternoon

Morning · Why we care about likelihood

Every model you fit in this course is, secretly, a statement of the form $Y \sim f(\theta)$: data $Y$ are random draws from a distribution whose shape is controlled by parameters $\theta$. Once you write down $f$, the data become a function of $\theta$:

Same formula, same numbers. Only the interpretation flips. As a function of $y$ it is a probability (or probability density). As a function of $\theta$ with $y$ held fixed at the value you actually observed, it is the likelihood. The likelihood is not a probability over $\theta$. It is a measure of how compatible different $\theta$ values are with the data you have.

Read the full plain-language summary for Day 4 →. The same ideas are worked through slowly, with examples and no math.

Core concepts you will use all morning

Walk-through 1 · Bernoulli warmup: proportion of successes

Cleanest possible example. We watch 50 cars approach an intersection on a yellow-then-red light; 30 run the red, 20 stop. Each car is an independent Bernoulli trial with unknown probability $p$. The data $y = (y_1, \dots, y_{50})$ are 0s and 1s, with $\sum y_i = 30$.

The likelihood is

Take the log:

Differentiate, set to zero:

The proportion of successes is the MLE under Bernoulli/Binomial sampling. Same answer if we used a single Binomial$(n=50, p)$ likelihood. Only the binomial coefficient differs, and it doesn't depend on $p$.

p_grid <- seq(0.01, 0.99, by = 0.001)
L <- p_grid^30 * (1 - p_grid)^20
nLL <- -30 * log(p_grid) - 20 * log(1 - p_grid)

par(mfrow = c(1, 2))
plot(p_grid, L, type = "l", lwd = 2, xlab = "p", ylab = "Likelihood")
abline(v = 0.6, col = "red")
plot(p_grid, nLL, type = "l", lwd = 2, xlab = "p", ylab = "-log L")
abline(v = 0.6, col = "red")

Walk-through 2 · Grid search: Von Bertalanffy growth

The Von Bertalanffy growth function (VBGF) describes asymptotic growth in fishes:

where $a$ is age, $L_a$ is length at age $a$, $L_\infty$ is asymptotic maximum length, and $\kappa$ is the growth coefficient. Given $n$ aged individuals, the likelihood is

The original course code holds $\sigma = 10$ and $L_\infty = 100$ fixed and grids over $\kappa$ to show how the likelihood changes with the parameter. Here is that grid search on a grid fine enough to locate the minimum ($\kappa$ from $0.05$ to $0.50$ in steps of $0.005$), complete and runnable.

# Load the Vonbert data (columns: sex, age, length) and keep the sex == 1 fish.
TheData <- read.csv(here::here("Day4_Likelihood_BayesI", "data", "vonbert.csv"))
Sex1 <- TheData[TheData$sex == 1, ]
age1 <- Sex1$age; len1 <- Sex1$length

# negative log-likelihood as a function of (Linf, kappa, sigma)
nll_vbgf <- function(Linf, kappa, sigma) {
 pred <- Linf * (1 - exp(-kappa * age1))
 -sum(dnorm(len1, mean = pred, sd = sigma, log = TRUE))
}

# Grid search over kappa with Linf and sigma fixed (as in the original lab).
kappa_grid <- seq(0.05, 0.50, by = 0.005)
nll_vec <- sapply(kappa_grid, function(k) nll_vbgf(Linf = 100, kappa = k, sigma = 10))
plot(kappa_grid, nll_vec, type = "l", lwd = 2,
 xlab = expression(kappa), ylab = "Negative log-likelihood")
abline(v = kappa_grid[which.min(nll_vec)], col = "red")

Walk-through 3 · Numerical optimization: mle() and optim()

The stats4::mle() wrapper around optim() takes a function of named parameters and a list of starting values. It returns the MLE, the Hessian-based standard errors, and lets you compute Wald or profile confidence intervals. We stay on the same sex == 1 subgroup as Walk-through 2, so the two results are directly comparable. The block below rebuilds age1/len1 from scratch so it also runs on its own.

library(stats4)

TheData <- read.csv(here::here("Day4_Likelihood_BayesI", "data", "vonbert.csv"))
Sex1 <- TheData[TheData$sex == 1, ]
age1 <- Sex1$age; len1 <- Sex1$length

nll_vbgf3 <- function(Linf, kappa, sigma) {
 pred <- Linf * (1 - exp(-kappa * age1))
 -sum(dnorm(len1, mean = pred, sd = sigma, log = TRUE))
}

fit <- mle(nll_vbgf3,
 start = list(Linf = 100, kappa = 0.2, sigma = 10),
 method = "L-BFGS-B",
 lower = c(10, 0.001, 0.01),
 upper = c(500, 2, 100))
summary(fit)
confint(fit) # profile-likelihood CIs (the right kind for nonlinear models)
-logLik(fit)[1] # minimum NLL at the MLE

Walk-through 4 · Two-component Poisson mixture (whales)

Bowhead whale blow counts come from a mix of males and females, and the two sexes blow at different average rates. Treat each whale's sex as an unobserved (latent) label and model the counts as a two-component mixture, a male component with rate $\lambda_1$ and a female component with rate $\lambda_2$. A whale is male with probability $\pi$ (the proportion of males in the population) and female with probability $1-\pi$. The goal is to estimate all three quantities, the male blow rate $\lambda_1$, the female blow rate $\lambda_2$, and the sex proportions $\pi$ (males) and $1-\pi$ (females). The likelihood for one count $y$ is a mixture:

Three parameters, no closed form. We minimize the joint negative log-likelihood with mle().

TheData <- read.csv(here::here("Day4_Likelihood_BayesI", "data", "whale.csv"))   # columns: sex (1/2, unused here), blows
blows <- TheData$blows

nll_mix <- function(lam1, lam2, p) {
 comp <- p * dpois(blows, lam1) + (1 - p) * dpois(blows, lam2)
 -sum(log(comp))
}

fit_mix <- mle(nll_mix,
 start = list(lam1 = 10, lam2 = 30, p = 0.4),
 method = "L-BFGS-B",
 lower = c(0.01, 0.01, 0.001),
 upper = c(Inf, Inf, 0.999))
summary(fit_mix)

Walk-through 5 · Likelihood ratio tests & profile likelihood

LRT: do the two sex groups share a growth curve?

Fit two nested Vonbert models to the full dataset (sex 1 and sex 2 combined):

• Model A (null, 3 parameters): one $L_\infty$, one $\kappa$, one $\sigma$, shared across sexes.
• Model B (alternative, 5 parameters): sex-specific $L_\infty$ and $\kappa$, shared $\sigma$.

Likelihood-ratio statistic: $\Lambda = 2(\ell_B - \ell_A)$. Under $H_0$: $\Lambda \sim \chi^2_{\Delta\,\text{df} = 2}$.

# Capitalised Sex/Age/Len = ALL 100 fish, both sexes. (Walk-throughs 2 and 3
# used lower-case age1/len1, the 50 fish with sex == 1. Do not mix them up.)
TheData <- read.csv(here::here("Day4_Likelihood_BayesI", "data", "vonbert.csv"))
Sex <- TheData$sex; Age <- TheData$age; Len <- TheData$length

nll_A <- function(Linf, kappa, sigma) {
 pred <- Linf * (1 - exp(-kappa * Age))
 -sum(dnorm(Len, pred, sigma, log = TRUE))
}
nll_B <- function(Linf1, kappa1, Linf2, kappa2, sigma) {
 pred <- ifelse(Sex == 1,
 Linf1 * (1 - exp(-kappa1 * Age)),
 Linf2 * (1 - exp(-kappa2 * Age)))
 -sum(dnorm(Len, pred, sigma, log = TRUE))
}

fitA <- mle(nll_A, start = list(Linf = 100, kappa = 0.2, sigma = 10),
 method = "L-BFGS-B", lower = c(10, 0.001, 0.01))
fitB <- mle(nll_B, start = list(Linf1 = 100, kappa1 = 0.2, Linf2 = 100, kappa2 = 0.2, sigma = 10),
 method = "L-BFGS-B", lower = c(10, 0.001, 10, 0.001, 0.01))

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

# AIC ranking on the same models
AIC(fitA, fitB)

Profile likelihood confidence intervals

Most likelihood surfaces in ecology are not symmetric around the MLE, so Wald CIs (MLE ± 1.96·SE) lie. The profile likelihood CI is the set of parameter values for which

i.e. all values whose deviance is within 1.92 NLL units of the minimum. confint() on an mle object gives exactly this. You can also draw it.

nll_p <- function(p) -(30 * log(p) + 20 * log(1 - p))
phat <- 0.6
best <- nll_p(phat)
pp <- seq(0.40, 0.80, length.out = 400)
plot(pp, nll_p(pp) - best, type = "l", lwd = 2,
 ylab = expression(-Delta * log~L), xlab = "p")
abline(h = 1.92, lty = 2, col = "red")
abline(v = phat, col = "grey")

Read the 95% CI off where the curve crosses the horizontal line at 1.92.

Bonus walk-through 6 · MLE for logistic regression, by hand vs glm()

This example shows that glm(y ~ x, family=binomial) is not magic. It is iteratively minimizing the same negative log-likelihood you can write down in five lines of R. Simulating data lets us check that glm() and your hand-rolled MLE agree to many decimal places.

Simulate $n=200$ presence/absence observations of a songbird detected as a function of forest cover $x$:

library(stats4)
set.seed(123)
n <- 200
x <- rnorm(n)
beta0_true <- -0.5; beta1_true <- 1.2
p <- 1 / (1 + exp(-(beta0_true + beta1_true * x)))
y <- rbinom(n, size = 1, prob = p)

# --- 1. Hand-rolled MLE ----------------------------------------------
nll_logit <- function(beta0, beta1) {
 eta <- beta0 + beta1 * x
 mu <- 1 / (1 + exp(-eta)) # inverse-logit
 -sum(dbinom(y, size = 1, prob = mu, log = TRUE))
}
fit_mle <- mle(nll_logit, start = list(beta0 = 0, beta1 = 0))

# --- 2. The R one-liner ----------------------------------------------
fit_glm <- glm(y ~ x, family = binomial(link = "logit"))

# --- 3. Side-by-side -------------------------------------------------
cbind(MLE = coef(fit_mle), GLM = coef(fit_glm), truth = c(beta0_true, beta1_true))
# Standard errors:
cbind(MLE = sqrt(diag(vcov(fit_mle))),
 GLM = sqrt(diag(vcov(fit_glm))))

The estimates agree precisely. The standard errors agree because both methods invert the same observed Fisher information.

Bonus walk-through 7 · Zero-inflated Poisson (a non-standard mixture)

You count fledglings per nest. Some nests fail entirely (true structural zeros, such as nest abandoned or depredated). The successful nests produce a Poisson count of fledglings, including occasional Poisson zeros. The mixture has two sources of zeros and is called a zero-inflated Poisson (ZIP):

Here $\pi$ is the probability of structural zero and $\lambda$ is the Poisson mean among non-structural nests. Base R's glm() has no zero-inflated family, so ZIP is not a one-line glm(...) call, but ready-made ZIP fitters do exist: pscl::zeroinfl() is the standard one, and the VGAM and glmmTMB packages fit ZIP models too. We write the negative log-likelihood by hand anyway, because it is a single line and doing so is the transferable skill for the non-standard models that have no canned fitter:

library(stats4)
set.seed(11)
n <- 300
pi_true <- 0.35; lam_true <- 2.5
z <- rbinom(n, 1, pi_true) # 1 = structural zero
y <- ifelse(z == 1, 0, rpois(n, lam_true))

# Plain Poisson MLE underestimates the mean and misses the zero inflation:
mean(y); var(y)
mean(y == 0); dpois(0, mean(y)) # observed zeros >> Poisson-predicted zeros

# ZIP negative log-likelihood:
nll_zip <- function(p, lambda) {
 if (p <= 0 || p >= 1 || lambda <= 0) return(1e10)
 prob_zero <- p + (1 - p) * exp(-lambda)
 ll_pos <- log(1 - p) + dpois(y, lambda, log = TRUE)
 ll <- ifelse(y == 0, log(prob_zero), ll_pos)
 -sum(ll)
}
fit_zip <- mle(nll_zip, start = list(p = 0.2, lambda = 1),
 method = "L-BFGS-B",
 lower = c(0.001, 0.001),
 upper = c(0.999, 50))
summary(fit_zip)
confint(fit_zip)

Morning exercises

Four exercises on the likelihood half of the day. Work them with a partner; open each reveal only after you have tried. You are not expected to finish these during the session.

  s  n  phat   Wald CI            profile CI
 30 50  0.60   (0.464, 0.736)     (0.462, 0.728)
  1 20  0.05   (-0.046, 0.146)    (0.003, 0.202)
  0 20  0.00   (0.000, 0.000)     (0.000, 0.092)

library(stats4)

# (d) red-light data: s = 30 successes in n = 50 trials
s <- 30; n <- 50
nll <- function(p) -(s * log(p) + (n - s) * log(1 - p))
fit <- mle(nll, start = list(p = 0.5), method = "L-BFGS-B", lower = 0.001, upper = 0.999)
phat <- coef(fit)
se   <- sqrt(diag(vcov(fit)))         # Wald SE from the curvature (inverse Hessian)
c(phat = phat, se = se)
c(lower = phat - 1.96 * se, upper = phat + 1.96 * se)   # Wald 95% CI
confint(fit)                                            # profile 95% CI

# (e) a rare event: compare both intervals as successes become scarce
wald_vs_profile <- function(s, n) {
  nll  <- function(p) -(s * log(p) + (n - s) * log(1 - p))
  phat <- s / n
  se   <- sqrt(phat * (1 - phat) / n)
  grid <- seq(0.0001, 0.9999, by = 0.0001)
  keep <- grid[nll(grid) - min(nll(grid)) <= 1.92]      # profile interval, read off a grid
  round(c(s = s, n = n, phat = phat,
          wald_lo = phat - 1.96 * se, wald_hi = phat + 1.96 * se,
          prof_lo = min(keep), prof_hi = max(keep)), 3)
}
rbind(wald_vs_profile(30, 50), wald_vs_profile(1, 20), wald_vs_profile(0, 20))

library(stats4)
TheData <- read.csv(here::here("Day4_Likelihood_BayesI", "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")

              proportion of zeros           fitted lambda        ZIP             AIC
              observed  Poisson   ZIP       Poisson   ZIP        pi      Poisson     ZIP   Poisson - ZIP
lambda = 1      0.563    0.512   0.563       0.670   0.927     0.278      681.0    672.4        8.6
lambda = 2.5    0.357    0.190   0.357       1.663   2.335     0.288     1104.6   1023.0       81.6
lambda = 8      0.310    0.004   0.310       5.500   7.968     0.310     2214.3   1363.1      851.2

library(stats4)

fit_both <- function(lam_true, pi_true = 0.35, n = 300) {
  set.seed(11)
  z <- rbinom(n, 1, pi_true)                     # 1 = structural zero
  y <- ifelse(z == 1, 0, rpois(n, lam_true))

  # (a) plain Poisson: one parameter
  nll_pois <- function(lambda) -sum(dpois(y, lambda, log = TRUE))
  fit_p <- mle(nll_pois, start = list(lambda = mean(y)),
               method = "L-BFGS-B", lower = 0.001, upper = 100)

  # (a) ZIP: two parameters, same NLL as the walk-through
  nll_zip <- function(p, lambda) {
    if (p <= 0 || p >= 1 || lambda <= 0) return(1e10)
    prob_zero <- p + (1 - p) * exp(-lambda)
    ll <- ifelse(y == 0, log(prob_zero), log(1 - p) + dpois(y, lambda, log = TRUE))
    -sum(ll)
  }
  fit_z <- mle(nll_zip, start = list(p = 0.2, lambda = 1),
               method = "L-BFGS-B", lower = c(0.001, 0.001), upper = c(0.999, 50))

  lam_p <- coef(fit_p)[["lambda"]]
  p_z <- coef(fit_z)[["p"]]; lam_z <- coef(fit_z)[["lambda"]]
  round(c(lambda_true = lam_true,
          obs_zeros   = mean(y == 0),
          pois_zeros  = dpois(0, lam_p),                   # (b) zeros predicted by Poisson
          zip_zeros   = p_z + (1 - p_z) * exp(-lam_z),     # (b) zeros predicted by ZIP
          pois_lambda = lam_p, zip_lambda = lam_z, zip_pi = p_z,
          AIC_pois = AIC(fit_p), AIC_zip = AIC(fit_z),
          dAIC = AIC(fit_p) - AIC(fit_z)), 3)
}
t(sapply(c(1, 2.5, 8), fit_both))

Afternoon · From likelihood to posterior

You spent the morning maximizing $L(\theta \mid y)$, finding the single $\theta$ that makes the data most plausible. Bayesian inference goes further. It gives you a distribution over $\theta$ that combines the data with what you knew (or assumed) before. Same likelihood, plus a prior, equals a posterior:

The marginal $p(y) = \int p(y\mid\theta) p(\theta)\,d\theta$ does not depend on $\theta$, so it acts only as a normalizing constant. For most realistic models it is impossible to compute in closed form, and that is exactly why MCMC exists. MCMC samples from the unnormalized product on the right, dodging the integral entirely.

Bayes' rule, worked out for a single coin

To anchor the algebra, fit a binomial model to a single observation: $y = 7$ successes in $n = 10$ trials. Use a flat prior $p \sim \text{Uniform}(0,1)$, which, note, is exactly a $\text{Beta}(1,1)$ distribution (a constant density on $[0,1]$):

Recognize this kernel? It is a Beta(8, 4) density (up to a constant). The normalizing constant exists in closed form via the Beta function:

This is the magic of conjugacy: the prior family is chosen so that the posterior stays in the same family. The flat $\text{Uniform}(0,1)$ prior is a $\text{Beta}(1,1)$, so this really is a Beta prior + Binomial likelihood → Beta posterior, and the posterior $\text{Beta}(8,4) = \text{Beta}(1+7,\, 1+3)$ simply adds the 7 successes and 3 failures to the prior's two parameters.

theta <- seq(0.001, 0.999, length.out = 500)
prior <- dbeta(theta, 1, 1)
likelihood <- dbinom(7, size = 10, prob = theta)
posterior <- dbeta(theta, 8, 4)
plot(theta, posterior, type = "l", lwd = 2, col = "darkred",
 ylab = "density (or scaled likelihood)", xlab = expression(theta))
lines(theta, prior, lty = 2)
lines(theta, likelihood / max(likelihood) * max(posterior), col = "darkblue")
legend("topleft", lty = c(2, 1, 1), col = c("black", "darkblue", "darkred"),
 legend = c("Prior Beta(1,1)", "Likelihood (scaled)", "Posterior Beta(8,4)"), bty = "n")

Conjugacy 1 · Beta–Binomial in full

For $y \mid p \sim \text{Bin}(n, p)$ and $p \sim \text{Beta}(\alpha, \beta)$,

The posterior just adds successes to $\alpha$ and failures to $\beta$.

The mean of a $\text{Beta}(\alpha, \beta)$ distribution shows why that is the right way to read it:

The prior mean has exactly the form of a success proportion, $\alpha$ successes out of $\alpha + \beta$ trials. The posterior mean is the same proportion after the real data are added: $y$ more successes on top, and $n$ more trials underneath. So the prior behaves as if you had already watched $\alpha + \beta$ trials and seen $\alpha$ successes. The posterior mean is also a weighted average of the prior mean and the observed proportion $y/n$,

so the more real trials you have compared with $\alpha + \beta$, the more the data dominate.

Worked example: spotted-owl nest survival. You monitored 12 nests; 9 fledged. Choose three priors and see how each posterior moves.

y <- 9; n <- 12
theta <- seq(0.001, 0.999, length.out = 500)
priors <- list(
 flat = c(1, 1), # Uniform
 skeptical = c(2, 8), # prior centered at 0.2
 informative = c(30, 10) # tight, centered at 0.75
)
par(mfrow = c(1, 3))
# names(priors) returns the list's labels, c("flat", "skeptical", "informative"),
# so the loop runs once per prior with nm holding that prior's name.
for (nm in names(priors)) {
 # Double brackets pull out the element itself: priors[["flat"]] is the vector c(1, 1).
 # Single brackets, priors["flat"], would instead return a one-item list still wrapping it.
 ab <- priors[[nm]]                              # ab[1] = alpha, ab[2] = beta
 post <- c(ab[1] + y, ab[2] + n - y)          # add successes to alpha, failures to beta
 plot(theta, dbeta(theta, post[1], post[2]), type = "l", lwd = 2,
 main = paste(nm, ": Beta(", post[1], ", ", post[2], ")", sep = ""),
 xlab = expression(theta), ylab = "posterior density")
 lines(theta, dbeta(theta, ab[1], ab[2]), lty = 2, col = "grey")
 abline(v = y / n, col = "red")
}
par(mfrow = c(1, 1))

Conjugacy 2 · Gamma–Poisson

For $y_1, \dots, y_n \mid \lambda \sim \text{Pois}(\lambda)$ and $\lambda \sim \text{Gamma}(\alpha, \beta)$ (shape-rate),

Again, the prior plays the role of pseudo-observations: $\alpha$ events accumulated over $\beta$ units of effort. After observing $n$ new units and $\sum y_i$ events, you just add them in. Note that the bookkeeping differs from the Beta–Binomial case. Here the prior's weight is carried by $\beta$ alone, the pseudo-effort that your $n$ new units are added to, and not by $\alpha+\beta$. That is a claim about weight, not about influence: $\alpha$ still moves the answer, because it sets the prior mean $\alpha/\beta$. Writing the posterior mean as a weighted average separates the two jobs,

So $\alpha$ decides what the prior claims and $\beta$ decides how loudly. Hold $\beta = 1$ and raise $\alpha$ from 2 to 200 in the tick example below. The prior's weight stays pinned at $1/11$, yet the posterior mean climbs from 2.45 to 20.45, an extravagant belief asserted with the authority of a single bird.

Worked example: a parasitologist counts $y = (3, 1, 4, 2, 0, 5, 2, 3, 1, 4)$ ticks on 10 birds (total 25 in 10 birds). Prior $\lambda \sim \text{Gamma}(2, 1)$ (weakly informative around mean 2 per bird).

y <- c(3, 1, 4, 2, 0, 5, 2, 3, 1, 4); n <- length(y)
alpha0 <- 2; beta0 <- 1
alpha1 <- alpha0 + sum(y)
beta1 <- beta0 + n
c(post_mean = alpha1 / beta1,
 post_sd = sqrt(alpha1 / beta1^2),
 MLE = mean(y)) # compare to MLE

lam <- seq(0, 6, length.out = 400)
plot(lam, dgamma(lam, alpha1, beta1), type = "l", lwd = 2,
 xlab = expression(lambda), ylab = "posterior")
lines(lam, dgamma(lam, alpha0, beta0), lty = 2, col = "grey")
abline(v = mean(y), col = "red")

Why conjugacy is not enough, and what MCMC buys you

Conjugacy gives elegant closed forms for one-parameter problems. Real models have:

• Nonlinear functions of parameters (logistic growth: $\mu = r - r/K \cdot x$ depends on the ratio $r/K$).
• Many parameters interacting (random effects, hierarchical means and variances).
• Awkward likelihoods (mixtures, occupancy, capture-recapture, state-space).

For these models the marginal $p(y) = \int p(y,\theta)\,d\theta$ is intractable. MCMC sidesteps the integral by generating dependent samples from the unnormalized posterior; the long-run frequency of those samples is the posterior itself. These are the samplers you will meet today and tomorrow:

So a single Nimble model can mix samplers: each parameter gets whichever one fits it. configureMCMC(model, print = TRUE) lists what Nimble chose.

The Metropolis algorithm in 30 lines of R

Same Beta–Binomial spotted-owl example, but pretending we don't know the closed form. The Metropolis recipe at each iteration:

1. Propose $\theta^\star \sim \text{Normal}(\theta_{\text{current}}, s)$ (symmetric proposal).
2. Compute the acceptance ratio $\alpha = \min\!\left(1, \dfrac{p(\theta^\star)\, L(\theta^\star \mid y)}{p(\theta_{\text{current}})\, L(\theta_{\text{current}} \mid y)}\right)$.
3. Draw $u \sim \text{Uniform}(0,1)$. If $u < \alpha$, accept; otherwise stay.

We will write it twice. First exactly as the recipe above reads, with densities multiplied, ratio taken, and compared with a uniform draw. Then again on the log scale, which is how every real sampler is written, once you have seen why.

y <- 9; n <- 12; alpha0 <- 2; beta0 <- 2

# The posterior is proportional to prior x likelihood. That product is all
# Metropolis needs: the normalising constant cancels in the ratio below.
post <- function(p) {
  if (p <= 0 || p >= 1) return(0)     # p is a probability; outside (0,1) it is impossible
  dbeta(p, alpha0, beta0) * dbinom(y, size = n, prob = p)
}

set.seed(2)
niter <- 20000; chain <- numeric(niter); chain[1] <- 0.5
accept <- 0; s <- 0.35                # proposal SD

for (i in 2:niter) {
  prop  <- rnorm(1, mean = chain[i - 1], sd = s)   # 1. propose
  alpha <- post(prop) / post(chain[i - 1])         # 2. acceptance ratio
  if (runif(1) < alpha) {                          # 3. accept w.p. min(1, alpha)
    chain[i] <- prop; accept <- accept + 1
  } else chain[i] <- chain[i - 1]
}
cat("acceptance rate =", round(accept / niter, 3), "\n")   # 0.365

That runs, and its posterior mean of 0.686 sits right on the closed-form Beta(11, 5) mean of 0.688. Note there is no min(1, ...) anywhere. If alpha exceeds 1 then runif(1) < alpha is certain, so the comparison already is the minimum.

n <- 200                                               # then try 700 and 1400
y_nests <- rep(c(1, 0), times = c(0.75 * n, 0.25 * n))   # 1 = fledged, 0 = failed
prod(dbinom(y_nests, size = 1, prob = 0.75))
# n =  200  ->  1.4e-49
# n =  700  ->  1.1e-171
# n = 1400  ->  0            underflowed

log_post <- function(p) {
  if (p <= 0 || p >= 1) return(-Inf)        # log(0) = -Inf: impossible
  dbeta(p, alpha0, beta0, log = TRUE) +     # log(prior x likelihood)
  dbinom(y, size = n, prob = p, log = TRUE) #   = log prior + log likelihood
}

niter <- 20000; chain <- numeric(niter); chain[1] <- 0.5
accept <- 0; s <- 0.35
set.seed(2)
for (i in 2:niter) {
  prop      <- rnorm(1, mean = chain[i - 1], sd = s)
  log_alpha <- log_post(prop) - log_post(chain[i - 1])  # ratio -> difference
  if (log(runif(1)) < log_alpha) {                      # u < alpha, both logged
    chain[i] <- prop; accept <- accept + 1
  } else chain[i] <- chain[i - 1]
}

burn <- 2000
cat("acceptance rate =", round(accept / niter, 3), "\n")
hist(chain[-(1:burn)], breaks = 60, freq = FALSE, main = "Posterior of p",
     xlab = "p", col = "grey90", border = "white")
# Overlay the closed-form Beta(alpha0+y, beta0+n-y) posterior
curve(dbeta(x, alpha0 + y, beta0 + n - y), col = "red", lwd = 2, add = TRUE)

Run both. With the same seed they accept and reject at exactly the same steps and produce a bit-for-bit identical chain. The logs change nothing about the algorithm. They only change the arithmetic it can survive.

The main event · logistic growth with Nimble

We now fit the ecological example the morning's lm() couldn't honestly answer, posterior distributions for the logistic-growth parameters $r$, $K$, and $\sigma$, plus derived quantities $K/2$ (population size at maximum growth) and $\text{MSY} = rK/4$ (theoretical maximum sustainable yield).

The model is the per-capita form of $dN/dt = rN(1 - N/K)$:

The complete Nimble workflow

library(nimble); library(MCMCvis); library(HDInterval); library(coda)
Logistic <- read.csv(here::here("Day4_Likelihood_BayesI", "data", "Logistic.csv"), header = TRUE)
Logistic <- Logistic[order(Logistic$PopulationSize), ]
N_grid <- seq(0, 1500, by = 10); N_grid[1] <- 1

logistic_code <- nimbleCode({
 # priors
 K ~ dunif(0, 4000)
 r ~ dunif(0, 2)
 sigma ~ dunif(0, 2)
 tau <- 1 / sigma^2 # precision = 1 / variance

 # likelihood
 for (i in 1:n) {
 mu[i] <- r - (r / K) * x[i]
 y[i] ~ dnorm(mu[i], tau)
 }

 # derived quantities (scalars -- monitoring these is free)
 N_at_maxdNdt <- K / 2
 MSY <- r * K / 4
})

# Nimble takes the inputs to a model as two lists:
# constants : sizes, indices, fixed covariates (things that are not random variables)
# data : observed stochastic nodes (y[i] here)
constants <- list(
 n = nrow(Logistic),
 x = as.double(Logistic$PopulationSize)
) # N_grid stays in R -- it is a plotting grid, not part of the model
data <- list(y = as.double(Logistic$GrowthRate))

inits_list <- list(
 list(K = 1500, r = 0.20, sigma = 1.0),
 list(K = 1000, r = 0.15, sigma = 0.1),
 list(K = 900, r = 0.30, sigma = 0.01)
)

set.seed(10)
samples <- nimbleMCMC(
 code = logistic_code,
 constants = constants,
 data = data,
 inits = inits_list,
 monitors = c("K", "r", "sigma", "N_at_maxdNdt", "MSY"),
 nchains = 3,
 nburnin = 5000,
 niter = 25000,
 thin = 1,
 samplesAsCodaMCMC = TRUE
)

MCMCtrace(samples, params = c("K", "r", "sigma"), pdf = FALSE, Rhat = TRUE, n.eff = TRUE)
MCMCsummary(samples, params = c("K", "r", "sigma", "N_at_maxdNdt", "MSY"), round = 3)

print(object.size(samples), units = "MB") # ~2.3 Mb, not 94.3

# P(r > 0.22 | data)
r_chain <- as.numeric(MCMCchains(samples, params = "r"))
K_chain <- as.numeric(MCMCchains(samples, params = "K"))
mean(r_chain > 0.22)

library(HDInterval)
N_grid <- seq(0, 1500, by = 10); N_grid[1] <- 1

# One grid point at a time, so the full 60,000 x 151 matrix never exists.
dNdt_summary <- sapply(N_grid, function(Nj) {
 d <- r_chain * Nj * (1 - Nj / K_chain)
 c(median = median(d), hdi(d, 0.95)) # median, lower, upper
})
med_dNdt <- dNdt_summary["median", ]
lo_dNdt <- dNdt_summary["lower", ]
hi_dNdt <- dNdt_summary["upper", ]

plot(N_grid, med_dNdt, type = "l", lwd = 2,
 ylim = range(med_dNdt, lo_dNdt, hi_dNdt),
 xlab = "Population size N", ylab = "Population growth rate dN/dt")
abline(h = 0, col = "grey")
lines(N_grid, lo_dNdt, lty = "dashed")
lines(N_grid, hi_dNdt, lty = "dashed")

Two kinds of credible interval: BCI vs HPDI

A 95% Bayesian credible interval (BCI, equal-tailed) cuts 2.5% off each tail of the posterior. A 95% highest posterior density interval (HPDI) is the shortest interval containing 95% of the posterior. By construction, every point inside has higher posterior density than every point outside.

• For symmetric posteriors (a typical sample-size-large $r$ posterior here), BCI and HPDI are nearly identical.
• For skewed posteriors (often $\sigma$, $K$ when the prior bound bites, ratios of parameters), HPDI is shorter and shifted toward the mode.
• HPDI is the right interval when you want to say "the 95% most-credible values of $\theta$."
• BCI is the right interval when you want a quantile-symmetric summary (e.g. for forecasts where the median matters).

# mu[i] = r - (r/K)*x[i] is deterministic, so rebuild its 60,000 posterior
# draws from the saved r and K columns rather than monitoring 50 extra nodes.
mu_draws <- sapply(constants$x, function(xi) r_chain - (r_chain / K_chain) * xi)

BCI <- apply(mu_draws, 2, quantile, probs = c(.025, .50, .975)) # 3 x n
HPDI <- apply(mu_draws, 2, function(v) hdi(v, .95)) # 2 x n

plot(constants$x, data$y, pch = 19,
 xlab = "Population size N", ylab = "Per-capita growth (1/yr)")
lines(constants$x, BCI[2, ], lwd = 2)
lines(constants$x, BCI[1, ], lty = 2, col = "red")
lines(constants$x, BCI[3, ], lty = 2, col = "red")
lines(constants$x, HPDI[1, ], lty = 2, col = "blue")
lines(constants$x, HPDI[2, ], lty = 2, col = "blue")

Bonus walk-through 8 · Bayesian binomial estimation of survival

This example shows beta-binomial conjugacy from end to end, with the prior justification spelled out and a comparison between (i) the closed-form posterior, (ii) a Nimble fit, and (iii) Metropolis. All three give the same answer. They have to, because the model is the same.

Setting. You banded 25 juvenile Western Bluebirds in spring; 18 were resighted alive the following spring. The biological literature suggests annual juvenile survival in this species is "around 0.45 with a fair amount of variability." We encode this prior as Beta(9, 11), which has mean $9/20 = 0.45$ and an effective sample size of 20 trials.

y <- 18; n <- 25; alpha0 <- 9; beta0 <- 11
alpha1 <- alpha0 + y # 27
beta1 <- beta0 + n - y # 18
c(mean = alpha1 / (alpha1 + beta1),
 median = qbeta(0.5, alpha1, beta1),
 BCI_lo = qbeta(0.025, alpha1, beta1),
 BCI_hi = qbeta(0.975, alpha1, beta1))

surv_code <- nimbleCode({
 phi ~ dbeta(9, 11) # prior
 y ~ dbinom(phi, n) # likelihood
})

set.seed(7)
surv_samps <- nimbleMCMC(
 code = surv_code,
 constants = list(n = 25),
 data = list(y = 18),
 inits = list(list(phi = 0.5), list(phi = 0.3), list(phi = 0.7)),
 monitors = "phi",
 nchains = 3, nburnin = 2000, niter = 10000,
 samplesAsCodaMCMC = TRUE)
MCMCsummary(surv_samps, round = 3)

# P(phi > 0.5 | data)?
mean(MCMCchains(surv_samps, "phi") > 0.5)
# closed-form check:
1 - pbeta(0.5, alpha1, beta1)

Bonus walk-through 9 · Bayesian Normal regression, with posterior predictive checks

Body mass as a function of body length in Pacific giant salamanders. We will (i) fit the regression in Nimble, (ii) use the posterior to predict new salamanders, and (iii) use those predictions to check whether the model can produce data that look like the real data. This is the workflow you should use on Bayesian models.

library(nimble); library(MCMCvis)

set.seed(21)
n <- 120
length_mm <- runif(n, 80, 250)
true_int   <- -15      # grams
true_slope <- 0.35     # grams per mm
true_sigma <- 4        # grams of scatter around the line
mass_g <- true_int + true_slope * length_mm + rnorm(n, 0, true_sigma)
plot(length_mm, mass_g, pch = 19, xlab = "length (mm)", ylab = "mass (g)")

reg_code <- nimbleCode({
  # priors, chosen on the scale of the data
  beta0 ~ dnorm(0, sd = 100)             # mean mass at length 0: very wide
  beta1 ~ dnorm(0, sd = 1)               # grams added per mm: a few tenths is plausible
  sigma ~ T(dnorm(0, sd = 20), 0, )      # scatter in grams, kept positive

  # likelihood
  for (i in 1:n) {
    mu[i] <- beta0 + beta1 * length_mm[i]
    mass[i] ~ dnorm(mu[i], sd = sigma)
  }
})

reg_samps <- nimbleMCMC(
  code = reg_code,
  constants = list(n = n, length_mm = length_mm),
  data = list(mass = mass_g),
  inits = list(list(beta0 = 0,   beta1 = 0,   sigma = 10),
               list(beta0 = -30, beta1 = 0.5, sigma = 5),
               list(beta0 = 10,  beta1 = 0.2, sigma = 15)),
  monitors = c("beta0", "beta1", "sigma"),
  nchains = 3, nburnin = 3000, niter = 15000,
  samplesAsCodaMCMC = TRUE)
MCMCsummary(reg_samps, round = 3)

Step 1 · The posterior predictive distribution of one new salamander

The posterior tells you which values of $\beta_0$, $\beta_1$ and $\sigma$ are plausible. It does not yet tell you what a new salamander will weigh. Two kinds of uncertainty go into that: uncertainty about where the line is, and the scatter of individual salamanders around the line. The posterior predictive distribution includes both, and the recipe is short. For each posterior draw, take that draw’s parameter values and simulate one new observation from the likelihood, $\text{Normal}(\beta_0 + \beta_1 x,\ \sigma)$.

# One row per saved posterior draw; each row is one plausible set of parameter values
post <- MCMCchains(reg_samps, params = c("beta0", "beta1", "sigma"))
head(post)
S <- nrow(post)                        # 36,000 draws

# A new salamander that is 150 mm long
mu_150  <- post[, "beta0"] + post[, "beta1"] * 150        # average mass of 150 mm salamanders, one value per draw
new_150 <- rnorm(S, mean = mu_150, sd = post[, "sigma"])   # one simulated 150 mm salamander per draw

quantile(mu_150,  c(0.025, 0.5, 0.975))   # uncertainty about the AVERAGE 150 mm salamander
quantile(new_150, c(0.025, 0.5, 0.975))   # where ONE new 150 mm salamander is likely to fall

hist(new_150, breaks = 60, col = "grey85", border = "white",
     main = "Posterior predictive distribution, 150 mm", xlab = "mass (g)")
abline(v = quantile(mu_150, c(0.025, 0.975)), col = "blue", lwd = 2)    # interval for the average
abline(v = quantile(new_150, c(0.025, 0.975)), col = "red", lwd = 2)    # interval for one individual

The two intervals answer different questions. The interval for mu_150, about 36.4 to 38.0 g, is where the average 150 mm salamander sits; it is narrow because 120 salamanders pin the line down well. The interval for new_150, about 29.3 to 45.1 g, is where one new 150 mm salamander is likely to fall; it is much wider because it also includes the individual scatter $\sigma$. A manager asking “what will the next salamander we catch weigh?” needs the second.

Step 2 · Fake datasets from the posterior

Do the same thing for all 120 real lengths at once and you get a whole fake dataset: the data the model says you could have collected. Each posterior draw gives a different fake dataset.

# A whole fake dataset: one posterior draw, one simulated mass for each of the 120 real lengths
s <- 1
fake <- rnorm(n, mean = post[s, "beta0"] + post[s, "beta1"] * length_mm, sd = post[s, "sigma"])

par(mfrow = c(1, 2))
plot(length_mm, mass_g, pch = 19, main = "real data", xlab = "length (mm)", ylab = "mass (g)")
plot(length_mm, fake,   pch = 19, main = "one fake dataset", xlab = "length (mm)", ylab = "mass (g)")
par(mfrow = c(1, 1))

# Now 1,000 fake datasets, each from a different posterior draw (one per row)
nfake <- 1000
draws <- sample(S, nfake)
fake_data <- matrix(NA, nrow = nfake, ncol = n)
for (k in 1:nfake) {
  s <- draws[k]
  fake_data[k, ] <- rnorm(n, mean = post[s, "beta0"] + post[s, "beta1"] * length_mm,
                          sd = post[s, "sigma"])
}

Step 3 · The posterior predictive check

If the model is adequate, the real data should look like just another one of the fake datasets. To check that, pick a summary you care about, compute it on the real data and on every fake dataset, and see where the real value falls. Here the summary is the standard deviation of mass: can the model reproduce how spread out the real masses are?

# Does the model reproduce the spread of the real masses?
real_sd <- sd(mass_g)                   # one number from the real data
fake_sd <- apply(fake_data, 1, sd)      # the same summary for each of the 1,000 fake datasets

hist(fake_sd, breaks = 40, col = "grey85", border = "white",
     main = "SD of mass in 1,000 fake datasets", xlab = "SD of mass (g)")
abline(v = real_sd, col = "red", lwd = 2)

mean(fake_sd >= real_sd)                # posterior predictive p-value

The real SD, 17.8 g, sits in the middle of the fake SDs, and about 53% of fake datasets are at least as spread out as the real one. Your exact numbers will differ slightly from run to run.

Afternoon exercises

Three exercises on the Bayesian half of the day: conjugacy and two Nimble fits. Work them with a partner; open each reveal only after you have tried. You are not expected to finish these during the session.

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

library(nimble); library(MCMCvis); library(coda)

surv_code <- nimbleCode({
 phi ~ dbeta(a, b) # a, b passed in as constants: one model, three priors
 y ~ dbinom(phi, n)
})

fit_prior <- function(a, b) {
 set.seed(4)
 out <- nimbleMCMC(surv_code,
 constants = list(n = 25, a = a, b = b),
 data = list(y = 18),
 inits = list(list(phi = 0.3), list(phi = 0.5), list(phi = 0.8)),
 monitors = "phi", nchains = 3, nburnin = 2000, niter = 22000,
 samplesAsCodaMCMC = TRUE)
 phi <- as.numeric(MCMCchains(out, "phi"))
 c(mean = mean(phi),
 quantile(phi, c(.025, .975)),
 p_gt_0.5 = mean(phi > 0.5),
 Rhat = MCMCsummary(out)$Rhat[1])
}

priors <- list(c(1,1), c(9,11), c(50,50), c(50,100))
round(t(sapply(priors, function(ab) fit_prior(ab[1], ab[2]))), 3)

# Exact check -- no MCMC needed, because Beta-Binomial is conjugate:
# posterior is Beta(a + 18, b + 7)
qbeta(c(.025, .975), 19, 8) # 0.5221 0.8567
pbeta(0.5, 19, 8, lower.tail = FALSE) # 0.9855