Approximate Bayesian Computation

Last time

Last time, we used rejection sampling to sample from a posterior distribution.

We saw that sometimes we could use a summary of the data instead of the full dataset for rejection sampling, but only if it was “sufficient.”

We saw that the procedure ran into trouble when we had many observations to match and that it didn’t work with continuous data.

Today

Today we will:

• Introduce ABC to solve the problems from last time
• Investigate what summaries to retain and discuss diagnostics for whether there is more information left
• Use regression to improve ABC analysis (semi-automatic ABC + regression adjustment)

First example: normal location model

\[ \theta \sim \text{Uniform}[-3,3], \qquad X \mid \theta \sim N(\theta, 1) \]

Suppose we observe \(x_{\text{obs}} = 0.8\).

Normal example, step 1 in code

set.seed(681)
simulate_normal <- function() {
  theta <- runif(1, -3, 3)
  X <- rnorm(1, mean = theta,
             sd = 1)
  return(list(theta = theta,
              X = X))
}

import numpy as np
rng = np.random.default_rng(681)

def simulate_normal():
    theta = rng.uniform(-3, 3)
    X = rng.normal(theta, 1)
    return {"theta": theta, "X": X}

Try rejection sampling

rejection_sample_normal <-
  function(x_obs = 0.8) {
    while (TRUE) {
      run <- simulate_normal()
      if (run$X == x_obs)
        return(run$theta)
    }
  }

X <- replicate(100000,
  simulate_normal()$X)
sum(X == 0.8)   # 0

def rejection_sample_normal(
        x_obs=0.8):
    while True:
        run = simulate_normal()
        if run["X"] == x_obs:
            return run["theta"]

X = np.array(
    [simulate_normal()["X"]
     for _ in range(100_000)])
print(np.sum(X == 0.8))   # 0

Doesn’t work! \(P(X = 0.8) = 0\), so rejection_sample_normal() never returns.

Approximate Bayesian computation

The ABC algorithm, given a model \(\theta \sim \pi, X \sim P_\theta\), a summary statistic \(S\) and a distance function \(d\):

1. Draw \(\theta^* \sim \pi\)
2. Draw \(X^* \sim P_{\theta^*}\)
3. If \(d(S(X^*), S(x_{\text{obs}})) < \epsilon\), keep \(\theta^*\).
4. Repeat and study the retained parameters.

The two approximations in ABC

ABC as described above is rejection sampling on the condition \(d(S(X^*), S(x_{\text{obs}})) < \epsilon\).

Therefore, the retained \(\theta^*\)’s are samples from the posterior distribution of \(P(\theta \mid d(S(X), S(x_{\text{obs}})) < \epsilon)\).

Our target posterior is \(P(\theta \mid X = x_{\text{obs}})\).

• Approximation 1: Suppose we take \(S(x) = x\). We are conditioning on \(d(X, x_{\text{obs}}) < \epsilon\) rather than on \(X = x_{\text{obs}}\).
• Approximation 2: If \(S\) is not the identity function, we are comparing a summary of the data rather than the full data.

What does \(\epsilon\) do?

Go back to the simple model:

\[ \theta \sim \text{Uniform}[-3,3], \qquad X \mid \theta \sim N(\theta, 1) \]

with \(x_{\text{obs}} = 0.8\).

Take \(S(x) = x\) (the identity summary, so no information loss), \(d(x,y) = |x - y|\).

Compare

\[ \epsilon=1.5,\qquad 0.2,\qquad 0.05. \]

Tolerance changes the distribution we target

Compare the tolerances with Q–Q plots

Which comparison departs most from the equality line? What do we gain by reducing \(\epsilon\) from 0.2 to 0.05?

Here 0.2 and 0.05 are both small relative to the measurement noise standard deviation of 1. Their targets are similar, while the smaller tolerance retains fewer draws.

What if we have twenty measurements?

Now suppose we observe twenty independent calibration measurements:

\[ X_i\mid\theta\stackrel{\mathrm{iid}}{\sim}N(\theta,1), \qquad i=1,\ldots,20. \]

The prior and the unknown parameter stay the same.

A direct extension retains a dataset only when every simulated measurement is close to its observed counterpart:

\[ \max_i |X_i-x_{\mathrm{obs},i}|<\epsilon. \]

How often do you expect that to happen?

Matching twenty measurements in code: R

set.seed(681)
x_obs <- rnorm(20, 0.8, 1)   # observed data
X <- rnorm(20, 0.8, 1)       # simulated at theta = 0.8

plot_match <- function(x_obs, X, eps) {
  near <- abs(X - x_obs) < eps
  i <- seq_along(x_obs)
  plot(i, x_obs, pch = 4, ylim = range(x_obs - eps, x_obs + eps, X),
       xlab = "measurement i", ylab = "value",
       main = paste(sum(near), "of", length(X), "within", eps))
  segments(i, x_obs - eps, i, x_obs + eps)
  points(i, X, pch = 16, col = ifelse(near, "darkgreen", "gray"))
}
plot_match(x_obs, X, eps = 0.5)

Matching twenty measurements in code: Python

import matplotlib.pyplot as plt
rng = np.random.default_rng(681)
x_obs = rng.normal(0.8, 1, 20)   # observed data
X = rng.normal(0.8, 1, 20)       # simulated at theta = 0.8

def plot_match(x_obs, X, eps):
    near = np.abs(X - x_obs) < eps
    i = np.arange(1, len(x_obs) + 1)
    plt.errorbar(i, x_obs, yerr=eps, fmt="x", color="black")
    plt.scatter(i, X, color=np.where(near, "darkgreen", "gray"))
    plt.xlabel("measurement i"); plt.ylabel("value")
    plt.title(f"{near.sum()} of {len(X)} within {eps}")
    plt.show()
plot_match(x_obs, X, eps=0.5)

Even \(\theta = 0.8\) rarely matches every measurement

Each simulated measurement has to land inside its own bar. Here only a few do, so this run is rejected.

Recall from last time: sufficient statistics

Definition: A statistic \(S(X)\) is sufficient for \(\theta\) if the conditional distribution of the full data \(X\), given \(S(X)\), does not depend on \(\theta\).

If \(S\) is a sufficient statistic, then \(P(\theta \mid X = x_{\text{obs}})\) is the same as \(P(\theta \mid S(X) = S(x_{\text{obs}}))\).

If we choose \(S\) in ABC to be a sufficient statistic, the only approximation left is the tolerance (\(\epsilon\)).

If we choose \(S\) in ABC not a sufficient statistic, then the approximation could be good or bad depending on how much information the statistic retains.

Example of \(S\) in the normal model

For independent Normal measurements with known variance, the sample mean is sufficient for \(\theta\).

After we know the mean, the remaining configuration of the measurements provides no additional information about \(\theta\).

Consequently,

\[ \theta\mid X=x_{\mathrm{obs}} \quad\text{and}\quad \theta\mid\overline X=\overline x_{\mathrm{obs}} \]

have the same distribution under the same prior.

How might we check that the sample mean is sufficient?

Most of the time we won’t be able to derive whether a summary is sufficient. Can we check it by simulation instead?

Start from the definition: if \(S\) is sufficient, the distribution of \(X\) conditional on \(S(X)=S(x_{\text{obs}})\) does not depend on \(\theta\).

Under our joint model (\(\theta\) from the prior, \(X\) from the simulator), this means

\[ \theta\;\perp\!\!\!\perp\;X\mid S(X). \]

Once we know \(S(X)\), nothing else about \(X\) tells us anything more about \(\theta\).

What conditional independence predicts

For any feature \(h(X)\) and any sets \(A\) and \(B\),

\[ \begin{aligned} &P\bigl(\theta\in A,\ h(X)\in B \mid S(X)=S(x_{\text{obs}})\bigr)\\ &\quad= P\bigl(\theta\in A \mid S(X)=S(x_{\text{obs}})\bigr)\,P\bigl(h(X)\in B \mid S(X)=S(x_{\text{obs}})\bigr). \end{aligned} \]

For example, among runs with \(S(X)=S(x_{\text{obs}})\), cut \(\theta\) at its median and \(h(X)\) at its median.

Each of the four boxes should then contain about \(1/4\) of the runs.

If some boxes hold noticeably more than others, \(h(X)\) carries information about \(\theta\) that \(S(X)\) does not.

Which features should we look at?

We usually do not know in advance which feature of \(X\) might carry leftover information.

One option is to use random probes:

• Draw a random direction \(a\) and use the projection \(h(X)=a^\top X\).
• Cutting the projection at its median gives a random half-space \(\{x : a^\top x > c\}\).
• Repeat with many random directions.

The recipe

1. Simulate a large bank of \((\theta_i, X_i)\) from the prior and the simulator.
2. Compute \(S(X_i)\) and keep the runs with \(S(X_i)\) close to \(S(x_{\text{obs}})\).
3. Draw a random direction \(a\) and compute \(a^\top X_i\) for the kept runs.
4. Count the kept runs in the four \((\theta,\ a^\top X)\) boxes.
5. Repeat steps 3–4 for many directions.

We will compare two summaries of the twenty measurements:

• \(S(X)\) = mean of all twenty (sufficient)
• \(S(X)\) = mean of the first ten (not sufficient)

Checking a summary in code: R

set.seed(682)
n <- 20; N <- 200000
theta <- runif(N, -3, 3)
X <- matrix(rnorm(N * n, mean = theta), nrow = N)
check_summary <- function(S_X, S_x_obs, width = 0.05) {
  keep <- abs(S_X - S_x_obs) < width   # S(X) near S(x_obs)
  a <- rnorm(n)                      # random direction
  proj <- X[keep, ] %*% a            # projection a'X
  th <- theta[keep]
  table(th > median(th), proj > median(proj)) / sum(keep)
}

check_summary(rowMeans(X), mean(x_obs))
check_summary(rowMeans(X[, 1:10]), mean(x_obs[1:10]))

Row \(i\) of X is twenty measurements simulated with \(\theta_i\). Table rows: \(\theta\) above its median; columns: \(a^\top X\) above its median.

Checking a summary in code: Python

rng = np.random.default_rng(682)
n, N = 20, 200_000
theta = rng.uniform(-3, 3, N)
X = rng.normal(theta[:, None], 1, (N, n))
def check_summary(S_X, S_x_obs, width=0.05):
    keep = np.abs(S_X - S_x_obs) < width  # S(X) near S(x_obs)
    a = rng.normal(size=n)              # random direction
    proj = X[keep] @ a                  # projection a'X
    th = theta[keep]
    t_hi, p_hi = th > np.median(th), proj > np.median(proj)
    return np.histogram2d(t_hi, p_hi, bins=2)[0] / keep.sum()

check_summary(X.mean(axis=1), x_obs.mean())
check_summary(X[:, :10].mean(axis=1), x_obs[:10].mean())

Four random directions

Top: every box holds 1/4 of the runs, in every direction. Bottom: the boxes are uneven by different amounts in different directions.

Many random directions

For the sufficient summary, every direction stays within sampling noise. For the first-ten mean, many directions show dependence, though some do not.

Caveats

• Continuous summaries: \(S(X)=S(x_{\text{obs}})\) has probability zero, so we condition on a band. Within a band, \(\theta\) still affects where \(S(X)\) falls, and any projection related to \(S(X)\) picks this up. A narrower band reduces this but keeps fewer runs.
• Passing is not proof: We only examined linear projections and median splits, for finitely many directions. Other features could still depend on \(\theta\).
• Local check: We examined runs near \(S(x_{\text{obs}})\), under this prior.
• Monte Carlo noise: With finitely many runs, the boxes will not be exactly \(1/4\).

The distance \(d\)

If you have multiple summaries, you need a way to combine them.

• This is the same as in SMM
• Same ideas apply, although you wouldn’t be doing the two-step algorithm
• You want the distance to ensure that the summaries are roughly on the same scale
• You might want to take into account covariance between the summaries

The ABC approximations

First two are specific to ABC, last one is going to be a property of any sampling-based approach

As the tolerance decreases, ABC can approach the posterior given the summary.

Recovering the full-data posterior also requires that the summary preserve the relevant information.

Established approaches to choosing summaries

• Collect a large number of potenteial summaries, compare posteriors using different subsets. If posterior doesn’t change, when a summary is dropped, that is evidence that it doesn’t provide extra information. (Approximate sufficiency.)
• Simulate a large number of parameter-dataset pairs to learn a summary that predicts the parameter (semi-automatic ABC).

More in Prangle, Summary Statistics in ABC (appeared in “Handbook of Approximate Bayesian Computation,” but the arxiv version is open access.)

Semi-automatic ABC: Learn the summary

Choose candidate features \(z(X)\), which can include nonlinear transformations of the data.

1. Simulate training pairs \((\theta_i,X_i)\).
2. Fit a regression predicting \(\theta_i\) from \(z(X_i)\).
3. Freeze the fitted function \(\widehat m\) and define \(S(X)=\widehat m(z(X))\).
4. Run ABC with this summary on fresh simulations.

For a vector parameter, fit one prediction for each component.

Fearnhead and Prangle (2012) + arxiv version

Why predict the parameter?

Under squared-error loss, the best predictor from the full data is

\[ m(X)=E[\theta\mid X]. \]

If we could use this exact function as a summary, then

\[ E[\theta\mid m(X)]=m(X). \]

It preserves the posterior mean. It need not preserve posterior spread, tails, or multiple modes.

A learned summary for the rain model

Candidate features could include rainy-day count, weather changes, endpoint weather, and longest rainy run.

Use simulations to learn how these features jointly predict \(\theta\). Fitted by least squares on 150,000 simulated six-day records:

The learned summary \(S(X)\) is the fitted prediction. For DRRRRD, \(S(x_{\text{obs}})=0.50\).

The learned summary on held-out simulations

Even a good predictor has error because different parameter values can generate similar data.

Regression adjustment: Correct imperfect matches

The tolerance means that ABC mixes together \(\theta\)’s corresponding to summaries that don’t exactly equal the observed.

Write \(d_i=S(X_i)-S(x_{\mathrm{obs}})\), and imagine that

\[ \theta_i=\alpha+\beta\,d_i+e_i, \]

(this is not actually the case, \(\theta_i\) comes from a different process, but perhaps it holds approximately).

In that case, for every \(\theta_i\), we can estimate what \(\theta_i\) would have been if it had corresponded to \(X_i = x_{\text{obs}}\):

\[ \theta_i^{\mathrm{adj}}=\hat \alpha + \hat \beta \cdot 0 + \hat e_i=\widehat\alpha+\hat e_i. \] (because if \(S(X_i) = S(x_{\text{obs}})\), \(d_i = 0\)).

Beaumont, Zhang, and Balding (2002)

What the adjustment does

Each draw slides along the fitted line to \(S(x_{\mathrm{obs}})\), keeping its residual.

The low draw on the left moves up and the high draw on the right moves down, so they end up closer together: the adjustment removes spread from mixing over the band.

Example: Twenty measurements, wide tolerance

The posterior standard deviation is 0.22, so a tolerance of 0.5 is wide. Unadjusted, the kept draws have standard deviation 0.37.

Adjusted, they have standard deviation 0.22: the same as tolerance 0.05, which keeps ten times fewer draws.

When can adjustment help?

• When the tolerance is wide relative to the posterior spread, or the kept summaries are lopsided around \(S(x_{\mathrm{obs}})\).
• It needs a regression that describes how the posterior changes across the band. Curvature or changing spread calls for a nonlinear adjustment; bounded parameters might require a transformation.
• It cannot recover information absent from the summaries.

Two different uses of regression

The methods can be combined: learn a summary, run ABC, then adjust the retained parameters.

ABC can also compare empirical distributions

For IID observations, we could compare observed and simulated empirical distributions using MMD.

The kernel determines which differences the comparison emphasizes.

K2-ABC uses an MMD-based discrepancy to weight simulations, avoiding an explicit small vector of hand-chosen moments.

Kernel and bandwidth choices still affect inference.

Park, Jitkrittum, and Sejdinovic (2016)

Wrapping up

What we saw

• Two types of approximation to the desired posterior, one to do with the tolerance and one to do with the summaries chosen.
• Pretty much everything we talked about re summaries in SMM applies. If you don’t have enough you won’t be able to identify parameters, they need to be informative, if you have more than one you need to scale them appropriately.
• We saw a couple of ways to attempt to diagnose sufficiency/approximate sufficiency.
• Regression can be used to create summary statistics automatically and to “correct” for an overly large \(\epsilon\).
• Other SMM extensions, like MMD, can also be incorporated in ABC.

Further reading

• Sisson, Fan, and Beaumont, An Overview of Approximate Bayesian Computation. The rejection algorithm and its interpretation.

• Prangle, Summary Statistics in Approximate Bayesian Computation.

• Fearnhead and Prangle (2012), Constructing Summary Statistics for ABC: Semi-automatic ABC. Learning summaries by parameter prediction.

• Beaumont, Zhang, and Balding (2002), Approximate Bayesian Computation in Population Genetics. Local regression adjustment.

• Park, Jitkrittum, and Sejdinovic (2016), K2-ABC. MMD-based comparisons.