An introduction to Rejection ABC

The Uniform Kernel: Rejection ABC

If we choose the ABC kernel to be a uniform kernel then we obtain Rejection-ABC. A uniform kernel is essentially an indicator function determining whether the true observation $y^{\ast }$ is within an $ϵ$-ball of the pseudo-observation $y\sim p(y\mid \theta )$ or not. $$p_{ϵ}(y^{\ast }\mid y)=\frac{1}{\text{Vol}(ϵ)}{\mathbb{I}}_{B_{ϵ}(y)}(y^{\ast })\propto \{\begin{array}{ll}1& y^{\ast }\in B_{ϵ}(y)\\ 0& y^{\ast }\notin B_{ϵ}(y)\end{array}$$

In the expression above $B_{ϵ}(y)$ is a ball of radius $ϵ$ centered around the auxiliary data $y$: it contains all the data $y^{\prime }\in \mathcal{Y}$ such that it is at a distance less than or equal to $ϵ$, our tolerance. The distance is computed using a distance metric $\rho (\cdot ,\cdot )$ which is usually the Euclidean one $$B_{ϵ}(y)=\{y^{\prime }\in \mathcal{Y}:\rho (y,y^{\prime })\le ϵ\}.$$ The indicator function ${\mathbb{I}}_{B_{ϵ}(y)}(\cdot )$ then is equal to $1$ whenever $y^{\ast }$ is within an $ϵ$-ball of $y$ and $0$ otherwise. Notice how we divide the indicator function by the $\text{Vol}(ϵ)$, i.e. the volume of $B_{ϵ}(y)$ so that the kernel is normalized and is thus a valid probability density function.

To sample $N$ parameters from $p(\theta \mid y^{\ast })$ then we can follow the Rejection-ABC algorithm.

Essentially the choice of an indicator function as our kernel implicitly defines an algorithm that either rejects or accepts samples. To see this, notice that this is one of those cases where in practice we can only evaluate an unnormalized kernel $${\tilde{p}}_{ϵ}(y^{\ast }\mid y)={\mathbb{I}}_{B_{ϵ}(y)}(y^{\ast })$$ so that each ${\tilde{w}}_n$ is either $1$ or $0$. Suppose out of $w_1,\dots ,w_N$ there are only $M\le N$ non-zero weights. Then their sum is clearly $M$ and so this algorithm can be interpreted similarly to Algorithm $1$ where $w_n=\frac{1}{M}$ when $\rho (y_n,y^{\ast })\le ϵ$.

Of course, this is extremely inefficient in high dimensions as the probability of “hitting the data” $y^{\ast }$ with our $ϵ$-ball $B_{ϵ}(y)$ becomes negligible. If you compare this with Algorithm 1 you can see that accepted samples are those with weight of one $w_n=1$ while rejected samples have zero weight $w_n=0$.

As a side note, the posterior distribution targeted by Rejection Sampling is sometimes denoted $p(\theta \mid \rho (y,y^{\ast })\le ϵ)$ rather than the more general $p_{ϵ}(\theta \mid y^{\ast })$.

Two Moon Example

Again, we come back to the Two Moon example. Here is an example of how we could code a rejection sampler. Notice how I have set a maximum number of iterations to avoid the loop running indefinitely.

def rejection_abc(prior, simulator, distance, epsilon, N, y_star, maxiter=50):
    """General Rejection ABC Algorithm."""
    samples = np.zeros((N, len(y_star)))  # Accepted samples
    
    n = 0   # Number of accepted samples so far
    i = 0   # Number of iterations so far

    while (n < N) and (i < maxiter):

      # Sample from the prior
      theta_n = prior()

      # Run the simulator with new parameter
      y_n = simulator(theta_n)

      # Accept only if data is within epsilon distance from ystar
      if distance(y_n, y_star) <= epsilon:
          samples[n] = theta_n
          n += 1
          
      # Increase iteration count
      i += 1
    return samples

We can see below how Rejection ABC, imposing a “hard constraint”, gets fewer and fewer samples accepted as our tolerance $ϵ$ becomes smaller.

Once again, feel free to play around with this example with the following notebook.