Using early rejection Markov chain Monte Carlo and Gaussian processes to accelerate ABC methods

Let $\boldsymbol{\theta}$ be the parameter of interest, which belongs to a parameter space $\Theta\subset\mathbb{R}^{p}$, and let $\pi(\boldsymbol{\theta})$ be a prior density for $\boldsymbol{\theta}$. The likelihood function, denoted as $p(\cdot\mid\boldsymbol{\theta})$, is costly to compute, while it is possible to simulate data from the model with a given parameter $\boldsymbol{\theta}$. Consider a given observed data set $\boldsymbol{y}=(\boldsymbol{y}_{(1)},\boldsymbol{y}_{(2)},\ldots,\boldsymbol{y}_{(T)})^{\top}\in\mathcal{Y}^{T}\subset\mathbb{R}^{d\times T}$, where $d$ denotes the dimension of data and $T$ denotes the number of repeated measurements. The likelihood-free inference is achieved by augmenting the target posterior from $\pi(\boldsymbol{\theta}\mid\boldsymbol{y})\propto p(\boldsymbol{y}\mid\boldsymbol{\theta})\pi(\boldsymbol{\theta})$ to

where the auxiliary parameter $\boldsymbol{x}$ is a (simulated) data set from $p(\cdot\mid\boldsymbol{\theta})$, $\Delta:\mathcal{Y}^{T}\times\mathcal{Y}^{T}\rightarrow\mathbb{R}^{+}$ is a discrepancy function that measures the difference between the synthetic data set $\boldsymbol{x}$ and the observed data set $\boldsymbol{y}$. The function $K_{\varepsilon}(\cdot)$ is defined as $K_{\varepsilon}(\cdot)=\frac{K(\cdot/\varepsilon)}{K(0)}$, where $K(\cdot)$ is a centered smoothing kernel density. In this article, we mainly focus on a kernel function $K(z)$ that takes its maximum value at zero and is zero for $z>1$. For example, the uniform, Epanechnikov, tricube, triangle, quartic (biweight), triweight and quadratic kernels (Epanechnikov, 1969; Cleveland and Devlin, 1988; Altman, 1992). The parameter $\varepsilon>0$ serves as a tolerance level, and it weights the posterior $\pi_{\varepsilon}(\boldsymbol{\theta},\boldsymbol{x}\mid\boldsymbol{y})$ higher in regions where $\boldsymbol{x}$ and $\boldsymbol{y}$ are similar.

The goal of approximate Bayesian computation is to estimate the marginal posterior distribution of $\boldsymbol{\theta}$, denoted as $\pi_{\varepsilon}(\boldsymbol{\theta}\mid\boldsymbol{y})$, by integrating out the auxiliary data set $\boldsymbol{x}$ from the modified posterior distribution. This is given by the following equation:

The ABC posterior provides an approximation of the true posterior distribution when the tolerance level $\varepsilon$ approaches to zero.

In practice, we set $\Delta(\boldsymbol{x},\boldsymbol{y})=\rho(\boldsymbol{S}(\boldsymbol{x}),\boldsymbol{S}(\boldsymbol{y}))$ for some distance metric $\rho(\cdot,\cdot)$ and a chosen vector of summary statistics $\boldsymbol{S}(\cdot)$: $\mathcal{Y}^{T}\rightarrow{\mathbb{R}^{d_{s}}}$ with some small $d_{s}$. For example, given an observed data set with $T$ independent and identically distributed data points, we can take $\boldsymbol{S}(\boldsymbol{y})=\sum_{i=1}^{T}\boldsymbol{y}_{i}/T$ and $\rho(\cdot,\cdot)=\left\|\cdot\right\|_{2}$. If the vector of summary statistic is sufficient for the parameter $\boldsymbol{\theta}$, then comparing the summary statistics of two data sets will be equivalent to comparing the data sets themselves (Brooks et al., 2011). In this case, the ABC posterior converges to the conditional distribution of $\boldsymbol{\theta}$ given $\boldsymbol{y}$ as $\varepsilon\rightarrow 0$; otherwise, the ABC posterior converges to the conditional distribution of $\boldsymbol{\theta}$ given $\boldsymbol{S}(\boldsymbol{y})$, that is, $\pi(\boldsymbol{\theta}\mid\boldsymbol{S}(\boldsymbol{y}))$. For convenience of description, we still denote $\pi_{\varepsilon}(\boldsymbol{\theta}\mid\boldsymbol{S}(\boldsymbol{y}))$ as $\pi_{\varepsilon}(\boldsymbol{\theta}\mid\boldsymbol{y})$.

Marjoram et al. (2003) proposed an ABC Markov chain Monte Carlo (MCMC) method based on the Metropolis–Hastings MCMC algorithm. The ABC-MCMC constructs a Markov chain that admits the target probability distribution as stationarity under mild regularity conditions. At $(n+1)$-th MCMC iteration, we repeat the following procedure.

1. (i)

   Generate a candidate value $\boldsymbol{\theta}^{*}\sim q(\cdot\mid\boldsymbol{\theta}_{n})$.

2. (ii)

   Generate a synthetic data set $\boldsymbol{x}^{*}=(\boldsymbol{x}_{(1)}^{*},\cdots,\boldsymbol{x}_{(T)}^{*})^{T}$ from the likelihood, $p(\cdot\mid\boldsymbol{\boldsymbol{\theta}^{*}})$.

3. (iii)

   Accept the proposed parameter value $(\boldsymbol{\theta}_{n+1},\boldsymbol{x}_{n+1})=(\boldsymbol{\theta}^{*},\boldsymbol{x}^{*})$ with probability

   $$\alpha[(\boldsymbol{\theta}_{n},\boldsymbol{x}_{n}),(\boldsymbol{\theta}^{*},\boldsymbol{x}^{*})]=\operatorname{min}\left\{1,\frac{\pi(\boldsymbol{\theta}^{*})q(\boldsymbol{\theta}_{n}\mid\boldsymbol{\theta}^{*})}{\pi(\boldsymbol{\theta}_{n})q(\boldsymbol{\theta}^{*}\mid\boldsymbol{\theta}_{n})}\frac{K_{\varepsilon}(\Delta(\boldsymbol{x}^{*},\boldsymbol{y}))}{K_{\varepsilon}(\Delta(\boldsymbol{x}_{n},\boldsymbol{y}))}\right\};$$

   (2)

   Otherwise, set $(\boldsymbol{\theta}_{n+1},\boldsymbol{x}_{n+1})=(\boldsymbol{\theta}_{n},\boldsymbol{x}_{n})$.

The MCMC algorithm targets the joint posterior distribution $\pi_{\varepsilon}(\boldsymbol{\theta},\boldsymbol{x}\mid\boldsymbol{y})$, with the marginal distribution of interest being the posterior distribution $\pi_{\varepsilon}(\boldsymbol{\theta}\mid\boldsymbol{y})$. The sampler generates a Markov chain sequence $\{(\boldsymbol{\theta}_{n},\boldsymbol{x}_{n})\}$ for $n\geq 0$, where $\boldsymbol{x}_{n}$ represents the synthetic data set, and $\boldsymbol{\theta}_{n}$ represents the parameter value.

For complex statistical models that are expensive to simulate, we aim to speed up inference by reducing the number of simulations. Based on a uniform kernel $K_{\varepsilon}(\Delta(\boldsymbol{x},\boldsymbol{y}))=\mathbb{I}(\Delta(\boldsymbol{x},\boldsymbol{y})\leq\varepsilon)$, Picchini (2014) proposed an early-rejection MCMC algorithm denoted by OejMCMC. Here, we generalize the uniform kernel to a more general kernel function. By the definition of the kernel function in Eq.(1), the acceptance probability $\alpha[(\boldsymbol{\theta}_{n},\boldsymbol{x}_{n}),(\boldsymbol{\theta}^{*},\boldsymbol{x}^{*})]$ (Eq. (2)) is always not greater than

Then we can reject some candidate parameters before simulations by generating $w\sim\mathcal{U}(0,1)$, and comparing $w$ with $\breve{\alpha}[(\boldsymbol{\theta}_{n},\boldsymbol{x}_{n}),\boldsymbol{\theta}^{*}]$. Specifically, If $w>\breve{\alpha}[(\boldsymbol{\theta}_{n},\boldsymbol{x}_{n}),\boldsymbol{\theta}^{*}]$, we immediately reject $\boldsymbol{\theta}^{*}$ before simulating synthetic data with $\boldsymbol{\theta}^{*}$, otherwise, simulate synthetic data and determine whether to receive candidate parameters according to Metropolis-Hastings acceptance probability.

In many cases, we may take a weak prior distribution or take a uniform distribution as prior, and use a symmetric proposal distribution, that is, for all $\boldsymbol{\theta},\boldsymbol{\theta}^{*}\in\Theta$, $\pi(\boldsymbol{\theta})\approx\pi(\boldsymbol{\theta}^{*})$ and $q(\boldsymbol{\theta}^{*}\mid\boldsymbol{\theta})=q(\boldsymbol{\theta}\mid\boldsymbol{\theta}^{*})$, $P(w>\breve{\alpha}[(\boldsymbol{\theta},\boldsymbol{x}),\boldsymbol{\theta}^{*}])\approx 0$. Hence, parameters are rarely rejected early, the method basically degenerates to the original ABC-MCMC. In Section 3, we will propose a new early rejection method based on discrepancy models.

Sisson et al. (2007); Del Moral et al. (2012) proposed ABC-SMC methods within the SMC framework of Del Moral et al. (2006). SMC is a class of importance sampling based methods for generating weighted particles from a target distribution $\pi$ by approximately simulating from a sequence of related probability distribution $\{\pi_{r}\}_{r=0:R}$, where the selected initial distribution $\pi_{0}$ is easy to approximate (e.g. the prior distribution), and the final distribution $\pi_{R}=\pi$ is the target distribution. The particles are moved from $r$-th target to $(r+1)$-th target, by using a Markov kernel.

In ABC-SMC, we select a sequence of intermediate target distributions $\{\pi_{r}=\pi_{\varepsilon_{r}}(\boldsymbol{\theta}\mid\boldsymbol{y})\}_{r=0:R}$ with a descending sequence of tolerance parameters $\{\varepsilon_{r}\}_{r=0:R}$. The initial distribution has a large tolerance and is easy to approximate, the final distribution has a small tolerance and is hence accurate. Let $\boldsymbol{\theta}_{r}^{(i)}$ denote the $i$-th particle after iteration $r$, let $\Delta_{r}^{(i)}$ denote the corresponding discrepancy and let $W_{r}^{(i)}$ denote the normalized weight. The set of weighted particles $\{\boldsymbol{\theta}_{r}^{(i)},W_{r}^{(i)}\}_{i=1:N}$ represents an empirical approximation of $\pi_{\varepsilon_{r}}(\boldsymbol{\theta}\mid\boldsymbol{y})$. At each iteration $r+1$, we conduct propagation, weighting and resampling to approximate $\pi_{\varepsilon_{r+1}}(\boldsymbol{\theta}\mid\boldsymbol{y})$. We refer reader to Appendix Section B.1 for more details of these three steps.

The main computational cost of ABC-MCMC comes from evaluating the acceptance ratio (Eq. 2), which requires simulating synthetic data sets. In this article, we propose to construct a pseudo MH acceptance probability

based on a discrepancy function $h$. Note that the pseudo MH acceptance probability $\hat{\alpha}[(\boldsymbol{\theta},\boldsymbol{x}),(\boldsymbol{\theta}^{*},\boldsymbol{x}^{*})]$ is not greater than

which does not depend on $\Delta(\boldsymbol{x}^{*},\boldsymbol{y})$. This allows us to reject some candidate parameters based on $\tilde{\alpha}[(\boldsymbol{\theta},\boldsymbol{x}),\boldsymbol{\theta}^{*}]$ before simulating the synthetic data. In the first stage of ejMCMC, we use $\tilde{\alpha}[(\boldsymbol{\theta},\boldsymbol{x}),\boldsymbol{\theta}^{*}]$ to early reject samples. Since $\tilde{\alpha}[(\boldsymbol{\theta},\boldsymbol{x}),\boldsymbol{\theta}^{*}]\leq\breve{\alpha}[(\boldsymbol{\theta}_{n},\boldsymbol{x}_{n}),\boldsymbol{\theta}^{*}]$, ejMCMC’s early rejection rate is always not lower than OejMCMC. Our ejMCMC algorithm degenerates to OejMCMC in case that $h(\boldsymbol{\theta})\equiv 0$. The proposed early rejection MCMC algorithm is shown in Algorithm 1.

If $h(\boldsymbol{\theta})$ satisfies that $h(\boldsymbol{\theta})\leq\Delta(\boldsymbol{x},\boldsymbol{y})$ in the ABC posterior region with positive support, the pseudo MH acceptance probability (Eq.(4)) is equal to the Metropolis-Hastings acceptance probability ${\alpha}[(\boldsymbol{\theta},\boldsymbol{x}),(\boldsymbol{\theta}^{*},\boldsymbol{x}^{*})]$ (Eq.(2)). However, the function $h(\boldsymbol{\theta})$ cannot strictly satisfy $h(\boldsymbol{\theta})\leq\Delta(\boldsymbol{x},\boldsymbol{y})$ in practice. Consequently, we will false reject samples in region $\{\boldsymbol{\theta}:h(\boldsymbol{\theta})>\Delta(\boldsymbol{x},\boldsymbol{y})\}$. Hence, the ejMCMC algorithm targets a modified posterior distribution based on a discrepancy function $h(\boldsymbol{\theta})$, which may not be equal to the original posterior distribution $\pi_{\varepsilon}(\boldsymbol{\theta},\boldsymbol{x}\mid\boldsymbol{y})$.

As we discussed above, the function $h(\boldsymbol{\theta})$ is a trade-off between speed and accuracy. In this article, we use a Gaussian process to model the functional relationship between the discrepancy and the parameters, and $h(\boldsymbol{\theta})$ is $a$ (for example, $a=0.05$) quantile prediction function of the deviation. More details of the selection and training of $h(\boldsymbol{\theta})$ will be described in Section 3.4.

Let $\Theta^{\prime}:=\{\boldsymbol{\theta}:h(\boldsymbol{\theta})<\varepsilon\}$ and $\Theta_{1}:=\{\boldsymbol{\theta}:\pi_{\varepsilon}(\boldsymbol{\theta}\mid\boldsymbol{y})>0\}$. Here $\Theta_{1}$ denotes the ABC posterior region with positive support. The following proposition shows that under mild assumptions, the ejMCMC algorithm satisfies the detailed balance condition.

This is a basic assumption for ejMCMC algorithm, without which the algorithm will collapse. All conclusions about the ejMCMC algorithm are based on this assumption.

Proposition 1 demonstrates that using the ejMCMC move does not violate the detailed balance condition of the Markov process. Moreover, Proposition 2 shows that the relationship between the posterior of the ejMCMC and that of the original ABC-MCMC.

Proposition 2 shows that the resulting Markov chain still converges to the true approximate Bayesian posterior under some conditions of $h(\cdot)$. This implies that the ejMCMC algorithm provides a promising framework to Bayesian posterior approximation, without compromising on accuracy or convergence guarantees under some mild conditions. In some other cases, there may exist some mismatch between the posterior of ejMCMC and ABC-MCMC. In this article, we utilize the $L_{1}$-distance shown in Eq. (7) to quantify the dissimilarity between two density functions.

Next we will prove the posterior concentration property of ejMCMC algorithm. Frazier et al. (2018) proved the posterior concentration property of the ABC posterior with the uniform kernel, that is, for any $\delta>0$, and for some $\varepsilon>0$,

where $d(\cdot,\cdot)$ is a metric on $\Theta$. Theorem 1 shows that the posterior concentration of the algorithm still holds, with some technical assumptions of function $h(\cdot)$. This property is important, because for any $A\subset\Theta$, ABC posterior probability $\Pi_{\varepsilon}\{A\mid\boldsymbol{y}\}$ will be different from the exact posterior probability. Without the guarantees of exact posterior inference, knowing that $\Pi_{\varepsilon}\{\cdot\mid\boldsymbol{y}\}$ will concentrate on the true parameter $\boldsymbol{\theta}_{0}$ can be used as an effective way to express our uncertainty about $\boldsymbol{\theta}$. The posterior concentration is related to the rate at which $\varepsilon$ goes to 0 and the rate at which the observed summaries $S(\boldsymbol{y})$ and the simulated summaries $S(\boldsymbol{x})$ converge to well defined limit counterparts $b(\boldsymbol{\theta}_{0})$ and $b(\boldsymbol{\theta})$. Here we consider $S(\boldsymbol{x})\rightarrow b(\boldsymbol{\theta})$ as $T\rightarrow\infty$, where $\boldsymbol{x}=(\boldsymbol{x}_{(1)},\ldots,\boldsymbol{x}_{(T)})^{\top}$ and $\boldsymbol{x}_{(t)}\sim p(\cdot\mid\boldsymbol{\theta})$ for $t=1,\ldots,T$, and consider $\varepsilon$ as a $T$-dependent sequence satisfying $\varepsilon_{T}\rightarrow 0$ as $T\rightarrow\infty$. Here, $\Pi(\cdot)$ is a probability measure with prior density $\pi(\cdot)$.

Remark. Assumption 2 controls the proportion of $\boldsymbol{\theta}$ satisfying $h(\boldsymbol{\theta})\geq\varepsilon$, that will be rejected before simulating, in a neighbourhood of $\boldsymbol{\theta}_{0}$. Assumption 3 is required for ABC with a more general kernel than a uniform kernel. This assumption controls the proportion of acceptable $(\boldsymbol{\theta},\boldsymbol{x})$ that may be rejected early. Intuitively, Assumptions 2 and 3 about $h(\boldsymbol{\theta})$ are to ensure that a small proportion of $\boldsymbol{\theta}$ near $\boldsymbol{\theta}_{0}$ will be rejected before simulating. As detailed in Section 3.4.1, our proposed discrepancy model $h(\boldsymbol{\theta})$ can always ensure that only a small proportion of $\boldsymbol{\theta}$ near the true value are early rejected by adjusting the quantile prediction function. Without these assumptions about the early rejection stage, posterior concentration concentration about $\boldsymbol{\theta}_{0}$ cannot occur.

We have previously discussed the quality of the ABC estimates. Now, we turn to investigate the efficiency of the algorithm.

Proposition 3 shows that when the kernel function $K(\cdot)$ is uniform and the posterior of ejMCMC and OejMCMC are consistent, the efficiency of ejMCMC is expected to be higher since the term $R_{ej}-R_{Oej}$ in Equation (8) is always positive. This is because ejMCMC targets the posterior region more effectively by rejecting more candidates in the early rejection step, potentially reducing the computational cost of generating synthetic data sets. Furthermore, when the function $h(\boldsymbol{\theta})$ satisfies $h(\boldsymbol{\theta})>\varepsilon$ in all regions with negligible posterior probabilities, the efficiency of ejMCMC is maximized. This means that the discrepancy model $h(\boldsymbol{\theta})$ is able to identify high probability ABC posterior regions, leading to higher early rejection rates and improved efficiency of the algorithm. The ejMCMC algorithm can play a role in correcting prior information before simulations in some sense, when the prior information is poor. In this article, we use

to quantify the early rejection efficiency of algorithms. This metric can be seen as a regularization of the early rejection rate, ensuring that the efficiency value is bounded between 0 and 1, with 1 being the best possible efficiency where all rejected parameters are early rejected. Proposition 4 shows the effect of the function $h(\cdot)$ on the early rejection efficiency of the ejMCMC algorithm.

Combining Propositions 2 and 4, the probability of $\Delta(\boldsymbol{x},\boldsymbol{y})\leq h(\boldsymbol{\theta})$ keeps a balance between efficiency and accuracy.

ABC-MCMC algorithm has several drawbacks. Firstly, the convergence speed of MCMC algorithm is often slow, especially in high-dimensional spaces. This may require a large number of iterations to obtain accurate approximations. Secondly, the MCMC algorithm is sensitive to the choice of initial values. Improper initial values can cause the Markov chain to get stuck in local optima, preventing effective exploration of the entire state space. In comparison with ABC-MCMC, ABC-SMC employs a sequential updating of particle weights to approximate the target distribution, which effectively mitigates the risk of getting stuck in local optima. Furthermore, the parallelization in SMC algorithms enhance their efficiency in complex applications. In this section, we introduce an early rejection ASMC (ejASMC) based on ABC adaptive SMC (ABC-ASMC) (Del Moral et al., 2012) and the ejMCMC method. The sequence of thresholds and MCMC kernels can be pre-specified before running the algorithm, but the algorithm may collapse with improperly selected thresholds and kernels. In this article, we adaptively tune the sequence of thresholds and ejMCMC proposal distributions.

Adaptive selection of thresholds: We adapt the scheme proposed in Del Moral et al. (2012) for selecting the sequence of thresholds, based on the key remark that the evaluation of weights $W_{r+1}^{(i)}={W}_{r}^{(i)}K_{\varepsilon_{r+1}}(\Delta(\boldsymbol{y},\boldsymbol{x}_{\boldsymbol{\theta}_{r}^{(i)}})/K_{\varepsilon_{r}}(\Delta(\boldsymbol{y},\boldsymbol{x}_{\boldsymbol{\theta}_{r}^{(i)}}))$ does not depend on particles of iteration $r+1$, $\{\boldsymbol{\theta}_{r+1}^{(i)},\boldsymbol{x}_{r+1}^{(i)}\}_{i=1:N}$. We select the threshold by controlling the proportion of unique alive particles, which is defined as

where $W_{r+1}^{(i)}$ is a function of the threshold $\varepsilon_{r+1}$, and $U$ is the largest subset of the set $\{1:N\}$ such that for any two distinct elements $i$ and $j$ in $U$, $\boldsymbol{\theta}_{r+1}^{(i)}\neq\boldsymbol{\theta}_{r+1}^{(j)}$. The proportion of unique alive particles is also intuitively a sensible measure of ‘quality’ of our SMC approximation. The threshold $\varepsilon_{r+1}$ is selected to make sure that the proportion of alive particles is equal to a given percentage of the particles number

for $\gamma\in(0,1)$. In practice, bisection is used to compute the root of Equation (10). The parameter $\gamma$ is a “quality” index for the resulting SMC approximation of the target. If $\gamma\rightarrow 1$ then we will move slowly towards the target but the SMC approximation is more accurate. However, if $\gamma\rightarrow 0$ then we can move very quickly towards the target but the resulting SMC approximation will be unreliable.

Adaptive MCMC proposal distribution: For each intermediate target, the SMC algorithm applies an MCMC move with invariant density $\pi_{\varepsilon_{r+1}}(\boldsymbol{\theta}\mid\boldsymbol{y})$. MCMC moves are implemented by accepting candidate parameters with a certain probability given by the MH ratio of ejMCMC algorithm. For $\boldsymbol{\theta}_{r}^{(i)}$ with $w_{r+1}^{(i)}>0$, we generate $\boldsymbol{\theta}^{*},\boldsymbol{x}^{*}$ according to a proposal $q_{r+1}(\boldsymbol{\theta}^{*}\mid\boldsymbol{\theta}_{r}^{(i)})p(\boldsymbol{x}^{*}\mid\boldsymbol{\theta}^{*})$, then accept the candidate with probability $\hat{\alpha}[(\boldsymbol{\theta}_{r}^{(i)},\boldsymbol{x}_{r}^{(i)}),(\boldsymbol{\theta}^{*},\boldsymbol{x}^{*})]$. We can adaptively determine the parameters of the proposal $q_{r+1}(\cdot\mid\cdot)$ based on the previous approximation of the target $\pi_{\varepsilon_{r}}$. In this article, we use Gaussian distribution as the proposal distribution, and the Gaussian covariance matrix is determined adaptively by computing the covariance matrix of the weighted particles $\{\boldsymbol{\theta}_{r}^{(i)},W_{r}^{(i)}\}_{i=1:N}$.

The output of ejASMC is the weighted particles $\{\boldsymbol{\theta}_{R}^{(i)},W_{R}^{(i)}\}_{i=1:N}$, the empirical distribution of which is the ABC posterior estimate. We set $\varepsilon_{0}=\infty$, then the initial target distribution is the prior distribution. Thus the initial weighted particles $\boldsymbol{\theta}_{0}^{(i)}\sim\pi(\boldsymbol{\theta})$ and $W_{0}^{(i)}=1/N$ for $i=1,\cdots,N$. The detailed pseudo-code of ejASMC is shown in Appendix Section B.2.

We first illustrate the advantage of our proposed method via a toy example. Suppose that the model of interest is a mixture of two normal distributions

where $\phi(y;\mu,\sigma^{2})$ denotes the probability density function of $\mathcal{N}(\mu,\sigma^{2})$ evaluated at $y$. Suppose that the prior distribution is $\pi(\theta)=\mathcal{U}(-6,6)$, and the observed data associated with this model is $y_{0}=1$. The discrepancy we use is $\Delta(y,y_{0})=|y-1|$. We choose a uniform kernel with the threshold $\varepsilon=0.6$ for ABC algorithms. Note that we only have one observed data point in this example.

Firstly, we run rejection sampling ABC $10^{6}$ iterations and treat the resulting posterior distribution as ground truth. We then compare our ejMCMC with OejMCMC and GP-ABC (Järvenpää et al., 2018). We collect $2,000$ training data points from prior, which are then used to train the GP discrepancy model. The resulting GP model is used to approximate the posterior of GP-ABC, and we name this as GP-ABC₁. We then run GP-ABC, OejMCMC and ejMCMC with same computational budgets ($N=100,000$ iterations). This GP-ABC with $N=100,000$ is named as GP-ABC₂. A random walk proposal $q(\theta\mid\theta_{n})\sim\mathcal{N}(\theta_{n},0.3^{2})$ is implemented for MCMC algorithms. In the implementation of ejMCMC, the discrepancy model $h_{a}(\theta)$ is a GP model trained with $a=0.05$. Figure 1 shows the Gaussian process modeling for identifying the relationship between the discrepancy $\Delta$ and parameter $\theta$. The x-axis denotes the parameter and the y-axis denotes the value of the discrepancy. Black dots denote the discrepancies, the red line is the mean function of GP and the area within two grey lines is the 95% predictive interval. In this example, the discrepancy for parameters is multi-modal, GP discrepancy model doesn’t fit the ‘shape’ of the discrepancy correctly. Figure 1 displays the estimated posterior distributions provided by GP-ABC₁, GP-ABC₂, ejMCMC, OejMCMC and the ground truth (rejection sampling). It indicates that both GP-ABC algorithms fail to capture a bimodal shape of the posterior, due to that they fail to capture the bimodal shape of discrepancy. This is also demonstrated in Järvenpää et al. (2018). Our ejMCMC is able to capture both peaks of the posterior. The $L_{1}$ distance ($D_{L_{1}}$) between the estimated posterior by ejMCMC and the corresponding true posterior is 0.0559, lower than the $D_{L_{1}}$ between OejMCMC and ground truth (0.0748). We also implement ejMCMC with online-GP (Stanton et al., 2021) discrepancy model for this toy example. It is available at https://github.com/caofff/ejMCMC.

In this subsection, we investigate the performance of ABC algorithms via an ordinary differential equation (ODE) example. We use the R package deSolve (Soetaert et al., 2010) to simulate differential equations. We generate ODE trajectories according to the following system,

where $\theta_{1}=2$ and $\theta_{2}=1$, and the initial conditions $x_{1}(0)=7$ and $x_{2}(0)=-10$. The observations $\boldsymbol{y}_{i}$ are simulated from a normal distribution with the mean $x_{i}(t|\boldsymbol{\theta})$ and the variance $\sigma_{i}^{2}$, where $\sigma_{1}=1$ and $\sigma_{2}=3$. We generate $121$ observations for each ODE function, equally spaced within $[0,60]$. We use root mean squared error (RMSE) to model the discrepancy between the observed data and the pseudo-data. We take the following prior distributions $\pi(\boldsymbol{\theta})$: $\theta_{1}\sim\mathcal{U}(1.8,2.2)$, $\theta_{2}\sim\mathcal{U}(0.8,1.2)$.

We conduct a performance comparison of two early rejection MCMC algorithms under different threshold levels. A smaller threshold value implies that the ABC posterior is closer to the actual posterior, and consequently, more simulation costs are required. In the ABC-MCMC algorithm, the convergence speed slows down with smaller thresholds. The 1% and 5% quantiles of the discrepancy correspond to the parameters that conform to the prior distribution are 4.03 and 4.68, respectively. Here we study the performance of ABC algorithms with $\varepsilon=4.8$, $4.5$, $4.1$ and $4.05$. We observe the Markov chains mix poor when the threshold value is smaller. For each level of tolerance, we first run MCMC algorithms $1,000,000$ iterations with 3 different initializations. The Gelman-Rubin diagnostic (Gelman et al., 1995) is used to assess the convergence of the Markov chain. The resulting posterior can be regarded as ground truth. For each case, we also run both ejMCMC and OejMCMC $100,000$ iterations. We collect $3,000$ training data from prior for ejMCMC algorithm. Appendix Figure C.1 illustrates the efficacy of the Gaussian Process (GP) model’s fit and the positioning of the training data points. The discrepancy can be fitted well by the Gaussian process. The choice of $a$ in the GP discrepancy model is a trade off between efficiency and accuracy. For each threshold $\varepsilon$, we run the ejMCMC algorithm with $a=0.5$, $0.2$, $0.1$, $0.05$ and $0.01$. Table 1 shows the results for $\varepsilon=4.5$ by varying $a$ (i.e. $h(\boldsymbol{\theta})$). This demonstrates that the selection of $h(\boldsymbol{\theta})$ keeps a balance between accuracy and computing speed. Results for the remaining three thresholds are shown in Appendix Section C.1. To balance the efficiency and accuracy, we choose $a=0.05$.

Table 2 shows the number of synthetic data generated ($N_{\text{sim}}$), number of collected posterior samples ($N_{\text{acc}}$), and the $L_{1}$-distance between true value and posterior estimates varying threshold parameter $\varepsilon$. Under the settings of uniform prior and Gaussian proposal distribution, the efficiency of OejMCM is 0. While our proposed algorithm can achieve a considerable computational acceleration by saving $20\%-80\%$ of data simulation. The accuracy of the posterior distribution estimates has no drop compared to OejMCMC. The efficiency of algorithms with threshold values $4.8$ and $4.5$ is higher than cases with that of $4.1$ and $4.05$. When the threshold is set to 4.1 or 4.05, the proportion of candidate parameters lying outside of the posterior region is low, resulting in a decrease of efficiency.

Figure 2 shows the posterior distributions of the ODE model provided by the three ABC algorithms varying the threshold parameter $\varepsilon$. The distributions provided by three methods are close, though ejMCMC slightly misses the tail of posteriors. This is because all candidate parameters satisfying $h(\boldsymbol{\theta})>\varepsilon$ in the posterior region are early rejected. However, the false rejection avoids the MCMC sampler “sticking” in some local region for long periods of time (marked by the red box in the Figure 2). We show an example of ABC-MCMC traceplots in Appendix Figure C.2.

In this study, we also employ the ejASMC algorithm to infer the posterior distribution of the parameters, with $N=2048$ samples generated for each iteration and the threshold parameters are updated by retaining 1024 (the scale parameter $\gamma=0.5$) active particles. The termination condition for the algorithm is set such that the number of simulated data exceeds 50,000. As shown in Figure 3, the posterior distribution becomes increasingly concentrated near the true value as the number of iterations increases. The threshold value $\varepsilon$ of the final iteration is $3.84$. ABC-MCMC algorithms achieve very low acceptance probability using this threshold. Here, the early rejection efficiency of ejASMC algorithm is 0.3820. $L_{1}$ distances between the marginal distributions of parameters and exact Bayesian posterior are 1.0626 and 1.0247, which are both smaller than the results provided by ejMCMC and OejMCMC (Appendix Table C.5), due to a smaller threshold. The posterior density plot provided by ejASMC are shown in Appendix Figure C.3.

Consider the parameter inference of a multidimensional stochastic differential equation (SDE) model for biochemical networks (Golightly and Wilkinson, 2010). In the dynamic system, there are eight reactions (${R_{1},R_{2},\cdots,R_{8}}$) that represent a simplified model of the auto-regulation mechanism in prokaryotes, based on protein dimerization that inhibits its own transcription. The reactions are defined as follows