Most of the models on this blog share a property that is easy to miss because it is so convenient: you can write down the probability of the data. Once you have that, everything follows. optim finds the maximum, the Hessian gives standard errors, MCMC explores the posterior.
Plenty of ecological models do not have that property. An individual-based forest model, a stochastic population model observed through counting error, a coalescent genealogy: you can run these forwards all day, but the probability of your particular dataset is an integral over every unobserved path the simulator could have taken, and nobody can do that integral.
Approximate Bayesian computation is the response. It replaces “evaluate the likelihood at this parameter” with “simulate at this parameter and see whether the result looks like the data”. That sounds like a hack. It is not: done one particular way it is exactly correct, and the whole subject is about what you give up when you leave that way behind. This post builds rejection ABC by hand, checks it against a posterior we can compute in closed form, and then puts it on a model where no closed form exists.
Start where the answer is known. Twenty seeds are sown, 13 germinate, and we want the posterior for the germination probability under a flat Beta(1, 1) prior. Conjugacy gives Beta(14, 8) and we are done in one line. Pretend we are not.
The ABC algorithm in its oldest form is three lines: draw a parameter from the prior, simulate a dataset at that parameter, keep the parameter if the simulated dataset equals the observed one.
set.seed(288)
n_trial <- 20L; y_obs <- 13L; M1 <- 400000L
p_prior <- runif(M1) # draw from the Beta(1, 1) prior
y_sim <- rbinom(M1, n_trial, p_prior) # run the model forwards
p_acc <- p_prior[y_sim == y_obs] # keep the exact matches
length(p_acc)[1] 19054acc_rate <- length(p_acc) / M1
acc_theor <- 1 / (n_trial + 1)
post_m <- mean(p_acc); exact_m <- 14 / 22
post_s <- sd(p_acc); exact_s <- sqrt(14 * 8 / (22^2 * 23))
ks <- suppressWarnings(ks.test(p_acc, "pbeta", 14, 8))
round(c(acc_rate = acc_rate, acc_theory = acc_theor,
abc_mean = post_m, exact_mean = exact_m,
abc_sd = post_s, exact_sd = exact_s, ks_p = ks$p.value), 4) acc_rate acc_theory abc_mean exact_mean abc_sd exact_sd ks_p
0.0476 0.0476 0.6370 0.6364 0.1004 0.1003 0.3758 Those 19,054 accepted values are not approximate draws from the posterior. They are draws from the posterior, in the same sense that any rejection sampler produces exact draws. The one-line proof: the probability of proposing \(p\) and then matching is \(\pi(p) \times P(y^\ast = y \mid p)\), which is the prior times the likelihood, which is the posterior up to the constant that rejection sampling throws away for free.
Two numbers confirm it. The Kolmogorov-Smirnov test against Beta(14, 8) gives p = 0.376, and the acceptance rate is 0.04763 against a theoretical 0.04762. That second one is worth a moment. Under a uniform prior the acceptance rate is \(\int \binom{n}{y} p^y (1-p)^{n-y}\,dp = 1/(n+1)\) exactly, whatever y happens to be: 0.04762 here. It does not depend on the data at all, only on how much data there is. Twenty trials cost you a factor of 21. Two hundred trials would cost a factor of 201, and two hundred continuous measurements would cost you everything, which is the next section.
Now measure something continuous. Twenty-five plants, each height a draw from a normal with unknown mean and a standard deviation of 1 that we happen to know. The prior on the mean is normal with a standard deviation of 10, wide enough to be nearly flat over anything plausible.
The exact-match rule is now useless: the probability that a simulated dataset reproduces 25 real numbers to the last decimal is zero. Two changes rescue it, and they are the only two approximations in all of ABC.
First, compare summaries instead of raw data. Here the sample mean is sufficient, so nothing is lost by that step (the next post is about what happens when nothing is not what you lose). Second, accept when the summaries are close rather than identical: \(|S^\ast - S_{\text{obs}}| < \varepsilon\).
set.seed(2288)
n <- 25L; sigma_known <- 1; prior_sd <- 10
se <- sigma_known / sqrt(n) # sampling sd of the summary
y_data <- rnorm(n, mean = 2.0, sd = sigma_known)
s_obs <- mean(y_data)
# the conjugate answer, for grading
post_prec <- 1 / prior_sd^2 + n / sigma_known^2
post_mean <- (n * s_obs / sigma_known^2) / post_prec
post_sd <- 1 / sqrt(post_prec)
round(c(s_obs = s_obs, exact_mean = post_mean, exact_sd = post_sd), 4) s_obs exact_mean exact_sd
2.0841 2.0833 0.2000 M2 <- 2000000L
mu_prior <- rnorm(M2, 0, prior_sd)
s_sim <- rnorm(M2, mu_prior, se) # the summary of a simulated dataset
d <- abs(s_sim - s_obs)
eps_grid <- c(0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0)
sweep_tab <- do.call(rbind, lapply(eps_grid, function(e) {
keep <- mu_prior[d < e]
data.frame(eps = e, n_acc = length(keep), rate = length(keep) / M2,
abc_mean = mean(keep), abc_sd = sd(keep),
predicted_sd = sqrt(1 / (1 / prior_sd^2 + 1 / (se^2 + e^2 / 3))))
}))
print(sweep_tab, digits = 4, row.names = FALSE) eps n_acc rate abc_mean abc_sd predicted_sd
0.02 3146 0.001573 2.082 0.2000 0.2003
0.05 7783 0.003891 2.081 0.2029 0.2020
0.10 15617 0.007809 2.082 0.2081 0.2081
0.20 31468 0.015734 2.081 0.2316 0.2309
0.50 78379 0.039190 2.082 0.3501 0.3510
1.00 155687 0.077843 2.077 0.6093 0.6099
2.00 310533 0.155267 2.056 1.1681 1.1639The predicted_sd column is not a fudge factor, and it is the reason this example is here rather than a prettier one. Accepting whenever \(|S^\ast - S_{\text{obs}}| < \varepsilon\) is the same as saying that the observed summary was generated by the model and then shifted by a uniform error on \((-\varepsilon, \varepsilon)\). A uniform on that interval has variance \(\varepsilon^2/3\). So ABC with a uniform kernel does not fit your model: it fits your model plus one extra layer of made-up noise, and the posterior inherits the variance of that noise.
tol_cost <- within(sweep_tab, {
sd_ratio <- abc_sd / post_sd
rel_error <- abs(abc_sd - predicted_sd) / predicted_sd
})
print(tol_cost[, c("eps", "rate", "abc_sd", "predicted_sd", "rel_error", "sd_ratio")],
digits = 3, row.names = FALSE) eps rate abc_sd predicted_sd rel_error sd_ratio
0.02 0.00157 0.200 0.200 0.001416 1.00
0.05 0.00389 0.203 0.202 0.004390 1.01
0.10 0.00781 0.208 0.208 0.000241 1.04
0.20 0.01573 0.232 0.231 0.003051 1.16
0.50 0.03919 0.350 0.351 0.002562 1.75
1.00 0.07784 0.609 0.610 0.000959 3.05
2.00 0.15527 1.168 1.164 0.003602 5.84The prediction lands to within 0.4% across a hundredfold range of tolerance, which means the cost is not a vague warning but an equation you can solve before you start. At 0.05 the posterior is 1% wider than the truth and 0.39% of proposals survive. At 2 the posterior is 5.8 times too wide and 16% survive. There is no setting that is generous with both.
Notice also what the acceptance rate does in that table: it roughly doubles when the tolerance doubles. In one summary dimension the acceptance rate scales like \(\varepsilon\), so halving the extra variance is only twice as expensive. That linear exchange rate is the friendliest ABC ever gets, and it holds only for one summary. With k summaries the acceptance region is a ball of radius \(\varepsilon\) in k dimensions and its volume scales like \(\varepsilon^k\).
Time to earn the machinery. The Ricker map with process noise and Poisson counts (the noisy version is a standard test case for likelihood-free methods: Wood 2010) is:
\[N_t = N_{t-1}\,\exp\!\big(r - N_{t-1} + e_t\big), \qquad e_t \sim N(0, \sigma^2), \qquad y_t \sim \text{Poisson}(\phi N_t)\]
The counts y are all you see. The population N is latent, the noise e is latent, and the likelihood of a fifty-year count series is a fifty-dimensional integral over both. There is no closed form. A particle filter could estimate it, but nothing writes it down.
Simulating it, on the other hand, is four lines. Vectorise across proposals rather than looping over them one at a time and a hundred thousand simulated ecosystems take under a second.
sim_ricker <- function(log_r, sigma, phi, n, N0 = 1) {
M <- length(log_r)
N <- rep(N0, M)
Y <- matrix(0L, M, n)
for (t in 1:n) { # loop over TIME, not over proposals
N <- N * exp(log_r - N + rnorm(M, 0, sigma))
Y[, t] <- rpois(M, phi * N)
}
Y
}
summ <- function(Y) { # three summaries per simulated series
m <- rowMeans(Y)
sdv <- sqrt(rowMeans(Y^2) - m^2) * sqrt(ncol(Y) / (ncol(Y) - 1))
Yc <- Y - m
a1 <- rowSums(Yc[, -1, drop = FALSE] * Yc[, -ncol(Y), drop = FALSE]) / rowSums(Yc^2)
cbind(mean = m, cv = sdv / pmax(m, 1e-8), acf1 = a1)
}set.seed(3288)
true_log_r <- 1.5; true_sigma <- 0.3; phi_known <- 20; nyr <- 50L
y_ric <- as.vector(sim_ricker(true_log_r, true_sigma, phi_known, nyr))
S_obs <- summ(matrix(y_ric, nrow = 1))
round(S_obs, 4) mean cv acf1
[1,] 31.14 0.4398 -0.3743Those three summaries come from what the model is made of, not from a menu. The mean count tracks the equilibrium, which for this map is \(N^\ast = r\): the observed series gives 1.56 against a true 1.5. The coefficient of variation carries the process noise sitting on top of the Poisson noise. And the lag-1 autocorrelation reads the return rate of the map, which the linearisation puts at \(1 - N^\ast\), or -0.5 here, and which goes negative once density dependence overcompensates. That last summary lands at -0.37 rather than -0.5, for two reasons worth knowing: the linearisation only holds for small noise, and the Poisson layer piles uncorrelated variance on top, pulling the correlation towards zero exactly as measurement error attenuates a regression slope.
The distance needs one more decision. The three summaries live on wildly different scales (counts around 31, a ratio near 0.44, a correlation near -0.37), so a raw Euclidean distance would spend almost its whole budget matching the mean count and let the other two drift. Divide each summary by its spread across the prior draws first, which is what the standard implementations do (Beaumont 2002); the next post measures what happens if you do not.
set.seed(4288)
M3 <- 200000L
th_lr <- runif(M3, 0, 3) # prior on r
th_sg <- runif(M3, 0, 1) # prior on sigma
S <- summ(sim_ricker(th_lr, th_sg, phi_known, nyr))
scale_S <- apply(S, 2, mad) # spread under the prior predictive
Z <- sweep(sweep(S, 2, as.vector(S_obs), "-"), 2, scale_S, "/")
dist <- sqrt(rowSums(Z^2))
keep_frac <- 0.0025 # keep the closest 0.25 per cent
tol <- unname(quantile(dist, keep_frac))
sel <- dist <= tol
post <- data.frame(log_r = th_lr[sel], sigma = th_sg[sel])
c(kept = nrow(post), tolerance = round(tol, 4)) kept tolerance
500.0000 0.1517 ci <- function(x) round(c(mean = mean(x), quantile(x, c(0.025, 0.975))), 3)
rbind(log_r = ci(post$log_r), sigma = ci(post$sigma)) mean 2.5% 97.5%
log_r 1.554 1.383 1.725
sigma 0.347 0.244 0.474The posterior mean for the growth rate is 1.554 against a truth of 1.5, with a 95 per cent credible interval of [1.38, 1.72]. For the process noise it is 0.347 against 0.3, interval [0.24, 0.47]. Both cover, and the two are not trading off against each other: their correlation across the accepted draws is -0.052 (p = 0.24), so the mean count is doing the work for one parameter and the variability for the other. We never wrote a likelihood, and the only R functions involved were runif, rnorm, rpois and quantile.
Look again at what the fit cost: 200,000 simulations to buy 500 posterior draws. The other 99.75% of the work went in the bin, because rejection ABC proposes from the prior and the prior is mostly wrong. That waste is the entire reason the rest of this cluster exists, and it splits into two separate complaints that get confused constantly.
The first is statistical. Everything you accept is filtered through the summaries and the tolerance, so the thing you sampled is not the posterior given the data: it is the posterior given 3 numbers, blurred by a uniform kernel. More simulations will not fix either. Where those 3 numbers came from, and what it costs to pick them badly, is the subject of the next post.
The second is computational. The eps-posterior you settled for is a fixed target, and rejection sampling is simply a bad way to hit it. Proposing from the prior means proposing mostly into regions the data have already ruled out, and ABC-MCMC and sequential ABC both fix that by proposing from something better.
Keep the two apart, because they get confused constantly. Better sampling makes the answer cheaper. Only better summaries and smaller tolerances make it more correct.
Newsletter
Get new tutorials by email
New R and QGIS tutorials for ecologists, straight to your inbox. No spam; unsubscribe anytime.
By subscribing you agree to receive these emails and confirm your address once. See the privacy policy.