# GPS-ABC: Gaussian Process Surrogate Approximate Bayesian Computation

###### Abstract

Scientists often express their understanding of the world through a computationally demanding simulation program. Analyzing the posterior distribution of the parameters given observations (the inverse problem) can be extremely challenging. The Approximate Bayesian Computation (ABC) framework is the standard statistical tool to handle these likelihood free problems, but they require a very large number of simulations. In this work we develop two new ABC sampling algorithms that significantly reduce the number of simulations necessary for posterior inference. Both algorithms use confidence estimates for the accept probability in the Metropolis Hastings step to adaptively choose the number of necessary simulations. Our GPS-ABC algorithm stores the information obtained from every simulation in a Gaussian process which acts as a surrogate function for the simulated statistics. Experiments on a challenging realistic biological problem illustrate the potential of these algorithms.

GPS-ABC: Gaussian Process Surrogate Approximate Bayesian Computation

Edward Meeds Informatics Institute University of Amsterdam tmeeds@gmail.com Max Welling Informatics Institute University of Amsterdam welling.max@gmail.com

## 1 Introduction

The morphogenesis of complex biological systems, the birth of neutrons stars, and weather forecasting are all natural phenomena whose understanding relies deeply upon the interaction between simulations of their underlying processes and their naturally observed data. Hypotheses positing the generation of observations evolve after critical evaluation of the match between simulation and observation.

This hypothesis–simulation–evaluation cycle is the foundation of simulation-based modeling. For all but the most trivial phenomena, this cycle is grossly inefficient. A typical simulator is a complex computer program with a potentially large number of interacting parameters that not only drive the simulation program, but are also often the variables of scientific interest because they play a role in explaining the natural phenomena. Choosing an optimal value setting for these parameters can be immensely expensive.

While a single, optimal parameter setting may be useful to scientists, often they are more interested in the distribution of parameter settings that provide accurate simulation results [20, 3, 18]. The interaction between parameters can provide insight regarding not only the properties of the simulation program, but more importantly, the underlying phenomena of interest. The main challenges that we address in this paper are 1) simulation-based modeling in a likelihood-free setting (we do not have a model in the typical machine learning sense, and therefore we do not have a standard likelihood function), and 2) running simulations is very expensive.

The first challenge is partly addressed by the approximate Bayesian computation (ABC) approach to sampling in likelihood-free scenarios [25, 4]. The ABC approach will be described in a later section, but in short it uses the distance between simulated and observed data as a proxy for the likelihood term in the parameter posterior. ABC provides the necessary framework to make progress in simulation-based modeling, but it is a very inefficient approach, even on simple problems.

The second challenge is approached with a surrogate model in mind. This means that every simulation (parameters and result) is stored and used to maintain a surrogate of the mapping from parameters to result. By carefully constructing an approximate Markov chain Monte Carlo (MCMC) sampler, we are able to sample from the parameter distribution in such a way that our sampling error is controlled and decreases over time. The main advantage of this approach is that by accepting some bias, we are able to very efficiently sample from the approximate posterior because parameters can sometimes be accepted within MCMC with high confidence by relying on the surrogate and thus avoiding expensive simulations.

In this paper we present a procedure for approximate Bayesian inference using a Gaussian process surrogate for expensive simulation-based models. In our approach, simulations that are run over the course of the inference procedure are incorporated into a GP model, gradually improving the surrogate over time. Using MCMC with a Metropolis-Hastings (MH) accept/reject rule, our method uses the surrogate and its uncertainty to make all MH decisions. The key insight is that the uncertainty in the acceptance probability is used to determine if simulations are required to confidently proceed to the MH step. If the uncertainty is high, and we are likely to make an error, then more simulations are run.

## 2 Approximate Bayesian Computation

The current state-of-the-art for simulation-based model inference is likelihood-free or approximate Bayesian computation methods [23, 15]. In this section we briefly introduce likelihood-free inference with a focus on MCMC inference procedures as our modeling approach will extend naturally from this work.

In likelihood-free sampling, we do not have a model in the typical sense. Instead, we have access to a simulator that generates pseudo-data that, given an accurate model of the world, look like the observations. The goal is to infer the parameters of the simulator which produce accurate pseudo-data. Importantly, we do not have access to a tractable likelihood function. We now describe in detail the likelihood-free set-up and in particular MCMC sampling techniques.

One of the primary goals of Bayesian inference is to infer the posterior distribution of latent variables using its prior and its data likelihood :

(1) |

where is a vector of latent parameters and is the observed data set. In simulation-based modeling, we do not have access to the likelihood . Instead our model of the world consists of a simulator that generates samples (where we indicate that the simulator was run with parameters and returns pseudo-data ). Simulation results are then compared with the observations through a distribution , which measures how similar is to . The distribution is parameterized by that controls the acceptable discrepancy between and . We can thus approximately infer the posterior distribution as

(2) |

This approximate posterior is only equal to the true posterior for where is the delta-function. For the exact posterior, one could apply a rejection sampling procedure that would repeatedly sample , run a simulation , then accept only if it equals . For continuous data, acts as a slack variable because equality cannot be achieved. However, we prefer small because this will improve our approximation to the true posterior. Unfortunately, there remains an unavoidable trade-off between approximation bias (large for large ) and rejection rate (large for small ).

### 2.1 Marginal and Pseudo-Marginal ABC

Instead of the rejection sampler described in the previous section (which is hopeless in high dimensions), we now describe two MCMC procedures, the marginal sampler and the pseudo-marginal sampler [2]. At every iteration of MCMC we propose a new . Next we sample samples from a simulator. From these samples we approximate the marginal likelihood as follows,

(3) |

We accept the proposed parameter with probability equal to,

(4) |

where the estimate of the marginal likelihood in the denominator (based on ) is carried over from the previous iteration. It can be shown that this algorithm is an instance of the more general pseudo-marginal procedure [2] because the estimate of the marginal likelihood is unbiased. From this we can immediately infer that this Markov chain converges to the posterior . Interestingly, there is an alternative view of this sampler [22] that interprets it as a Markov chain over the extended state which also leads to the conclusion that the samples will asymptotically follow the distribution .

Unfortunately, it is well known that pseudo-marginal samplers can suffer from slow mixing. In particular, when through a “lucky” draw our marginal likelihood estimate in Eqn. 3 attains a large value, then it is very difficult for the sampler to mix away from that state. To avoid this behavior it is sometimes beneficial to re-estimate the denominator (as well as the numerator) in every iteration. This procedure is more expensive and does not guarantee convergence to (unless ), but can result in much better mixing. We will call this the marginal LF MCMC method [2].

While for the pseudo-marginal approach we can interpret the fluctuations induced by estimating the from a finite sample set as part of randomly proposing a new state, this is no longer true for the approximate marginal MCMC. For the latter it is instructive to study the uncertainty in the acceptance probability due to these fluctuations: repeatedly estimating with samples will produce a distribution over . Clearly for small the distribution will be wider while for very large it will approach a delta peak. Our approach uses this distribution directly to determine the confidence in the MH accept step, allowing it to implicitly set based on the local uncertainty. This will be further discussed in next sections.

## 3 The Synthetic Likelihood

We next discuss the approximation introduced by Wood [28], who models the simulated pseudo-data at parameter value using a normal distribution, i.e. . We will later replace this with a Gaussian process. The procedure is very simple: we draw samples from our simulator and compute first and second order statistics,

(5) | |||||

(6) |

Using estimators and we set . If we moreover also use a Gaussian kernel,

(7) |

where is the dimension of , we can then analytically integrate over in Eqn 2 giving the synthetic-ABC likelihood:

(8) | |||||

(9) |

which has the satisfying result that likelihood of is the density under a Gaussian model at each simulation parameter setting .

This approximation has two advantages. First, we can take the limit . This implies that the bias introduced by the need to use a distribution to measure the similarity between simulations and observations is now removed. But this is traded off with the bias introduced by modeling the simulations from with a normal distribution. The second advantage was the main motivation in [28], namely that this procedure is more robust for extremely irregular probability distributions as encountered in chaotic or near chaotic simulation dynamics.

A marginal sampler based on a Gaussian approximation (Algorithm 1) has the following acceptance probability:

(10) |

As with the marginal sampler of Section 2, the fact that we estimate first and second order statistics from a finite sample set introduces uncertainty in the accept probability : another run of simulations would have resulted in different values for these statistics and hence of the accept probability. See Figure 1 for an illustration. In the following section we will analyze the distribution over and develop a method to decide how many simulations we need in order to be sufficiently confident that we are making the correct accept/reject decision. Random acceptance probability distributions have been studied in general [16] and for the specific case of Gaussian log-energy estimates [9].

### 3.1 MCMC with a Random Acceptance Probability

We now make explicit the role of randomness in the MCMC sampler with synthetic (normal) likelihoods. At each iteration of the MCMC sampler, we compute estimators as before using Eqns 5 and 6. To estimate the distribution over accept probabilities we would need sets of simulations, which would be too expensive. Instead, we use our Gaussian assumption to derive that the variance of the mean is times the variance in the sample ,

(11) |

and similarly for . This shortcut is important because it allows us to avoid a significant number of expensive simulations and replace them with samples from a normal distribution.

Given our samples of , we can compute samples of by inserting them into the expression for the randomized MH accept probability:

(12) |

We now derive a procedure to estimate the probability of making an error in an accept/reject decision and a threshold for actually making the decision. The error of making an incorrect decision can either be measured conditioned on (the uniformly distributed draw used in the MH decision), or unconditionally, by integrating over . First we start with the conditional error which trivially extends to the unconditional error by averaging.

If , then we accept the MH proposal and move to the proposed state. The probability of making an error in this case is (i.e. the probability we should actually reject):

(13) |

Similarly, if then we reject, and the error is (i.e. the probability we should actually accept):

(14) |

The total conditional error is therefore:

(15) |

and the total unconditional error is:

(16) |

which can again be estimated by Monte Carlo or grid values of . The analytic calculation of is the area under the cumulative distribution function of folded at (see Figure 1). This integral is also known as the mean absolute deviation [29] which is minimized at the median of (the value of where the CDF equals ), justifying our decision threshold (also determined by samples ).

With this in hand, we now have the necessary tools to construct an adaptive synthetic-likelihood MCMC algorithm that uses as a guide for running simulations (Algorithm 2). At the start of each MH step, simulations are run for both and ; estimators are computed; then are sampled. Based on these samples, the median and is computed. Note that this phase of the algorithm is very cheap; here we are sampling from bivariate Gaussian distributions to compute Monte Carlo estimates of and , so can be set high without a computational hit, though in practice should be fine. While , new simulations are run and the estimators updated, along with new draws of , etc. The user-defined error threshold is a knob which controls both the accuracy and computational cost of the MCMC. New simulations can be run at either or ; we run simulations at both locations, though selecting one over the other based on the higher individual mean uncertainty could result in fewer simulations. As increases, the uncertainty around decreases such that ; once this occurs, the MH is now confident and it proceeds using the usual acceptance test, with as the acceptance threshold.

In many cases, the actual number of simulations required at each step can be quite small, for example when one parameter setting is clearly better than an another (where the median is at or close to or ). Nevertheless, there remains a serious drawback to this algorithm for expensive simulations: all simulations are discarded after each MH step; a great waste considering the result is a single binary decision. Using a Gaussian process surrogate, described in the next section, we will remember all simulations and use them to gradually improve the surrogate and as a consequence, eventually eliminate the need to run simulations.

## 4 Gaussian Process Surrogate ABC

As mentioned in the introduction, in many scientific disciplines simulations can be extremely expensive. The algorithms up till now all have the downside that at each MCMC update a minimum number of simulations need to be conducted. This seems unavoidable unless we store the information of previous simulations and use them to make accept/reject decisions in the future. In particular, we can learn the mean and covariance of the synthetic likelihood as a function of and as such avoid having to perform simulations to compute them. There is a very natural tool that provides exactly this functionality, namely the Gaussian process (GP). For our purposes, the GP will “store” the simulation runs for all simulations conducted during the MCMC run. We will use the GP as a “surrogate” function for the simulated statistics from which will be able to estimate the marginal likelihood values away from regions where actual simulations were run. Importantly, the GP provides us with uncertainty estimates of the marginal likelihood which will inform us of the need to conduct additional experiments in order to make confident accept/reject decisions. Going from the synthetic likelihood model to the GP represents a change from frequentist statistics in favor of (nonparametric) Bayesian statistics. Gaussian processes have recently also become a popular tool in the machine learning literature as surrogates of expensive regression surfaces, such as log-likelihoods [17]; optimization surfaces [8, 24]; and simulations of physical systems [6].

Our surrogate model and algorithm follow directly from the synthetic-ABC approximation and randomized acceptance algorithm. The main difference between the two is that in this paper, we model the statistics as independent Gaussian processes (recall is the dimensionality of ). We note that it would be better to model the statistics using a single joint Gaussian process. This can be done using “co-Kriging” or variants thereof [12, 7] (see also Section 6). Although the independence assumption may lead to overconfidence (because it is assuming–falsely–independent evidence), it is also more robust in high-dimensions where the estimator of the full output covariance has high variance (it overfits). It may be that the mentioned multi-output GPs can provide an appropriate solution by tying the covariance structure across parameter space using a small number of kernel hyperparameters. For the current experiments we found that independent GPs worked well enough to illustrate the algorithm’s potential.

For each , we have the following bivariate conditional distribution:

(17) |

where the mean and covariance are determined by the set of training inputs and training outputs (using only statistic ). They are given by the usual expressions for the GP mean and covariance,

(18) |

and

where is a by vector of kernel evaluations for the ’th Gaussian process between and the input training set , is the th kernel matrix evaluated on the training data; is the data noise term for the ’th statistic (used below in the acceptance ratio), is the by training output data set and means selecting column from the training data, and is a single kernel evaluation at , for Gaussian process , and similarly for , , .

The GPS-ABC algorithm is now run as follows (see Algorithm 3). At each MH step, using Eqn 17, for each , independent draws of and are sampled from their bivariate Gaussian process posterior. Note that this is in fact different from the synthetic-ABC model in two ways 1) there are no default simulations at each MH step; instead, the current surrogate model is used to predict both the expectation and uncertainty in simulation output and 2) the Gaussian process takes into account the correlations between and . As before, Monte Carlo statistics are computed from acceptance probabilities as follows,

(19) |

The error and acceptance threshold are computed from the samples; if , then a procedure for acquiring training points (i.e. simulations) is run, with the objective of reducing uncertainty for this specific MH step. Again, as with the adaptive synthetic likelihood algorithm, computing the samples is very cheap. The procedure is then free to select any input location to run a simulation (whereas before we were forced to run at either and ), though the new simulation should be impactful for the current MH step. This means that we can choose locations other than and , perhaps trying to limit the number of future simulation runs required in the vicinity. Analogous to acquisition functions for Bayesian optimization [8], actively acquiring points has the implicit goals of speeding up MCMC, reducing MCMC error, and limiting simulations. Once a new training point is selected and run, the training input-output pair is added to all Gaussian processes and the model hyperparameters may or may not be modified (with a small number of optimization steps).

The key advantage of GPS-ABC is that with increasing frequency, we will not have to do any expensive simulations whatsoever during a MH step because the GP surrogate is sufficiently confident about the marginal likelihood surface in that region of parameter space.

### 4.1 Theoretical Aspects of GPS-ABC

Two of the main contributions of GPS-ABC are MCMC under uncertainty and the introduction of memory into the Markov chain; we consider these steps as the only way to reduce the number of expensive simulations and as such a necessary ingredient to GPS-ABC. Nevertheless, they present major differences from typical Bayesian inference procedures.

We now address two major theoretical aspects of the GPS-ABC algorithm: the approximate and adaptive nature of GPS-ABC. Although we have postponed formal proofs for future work we have outlined their main arguments below.

GPS-ABC is approximate because at each MH-step there is some probability that the chain will make an error, and that this corresponds to an error in the stationary distribution of the Markov chain (i.e. it is an approximation to the stationary distribution). In [14], another approximate MCMC algorithm is presented and it provides a proof for an upper bound on the error in the stationary distribution. The main argument is that if the MH-step error is small and bounded (along with a few other mild conditions), then the error in stationary distribution is bounded as well. We feel GPS-ABC fits nicely into this same proof framework.

GPS-ABC is also adaptive since the approximation to the stationary distribution changes as more training points are added to the Gaussian process (we are learning the surrogate as we run the MCMC). Two of the major requirements for a valid adaptive MCMC algorithm are diminishing adaptation and ergodicity [19]. GPS-ABC satisfies the former as the number of training points acquired over an MCMC rapidly decreases over time. When the adaptation slows and becomes insignificant, the Markov chain resembles the algorithm in [14], which, as we stated above, provides a proof of bounded convergence to the stationary distribution (and hence ergodicity); therefore GPS-ABC satisfies the latter requirement.

## 5 Experiments

We present three main experiments: 1) a toy Bayesian inference problem, 2) inference of parameters in a chaotic ecological system, and 3) inference of parameters in a stochastic biological system. The goal of the first experiment is to show the correctness of the adaptive synthetic likelihood algorithm (henceforth ASL-ABC) and GPS-ABC, and to demonstrate the computational improvements achieved by GPS-ABC over both ASL-ABC and Kernel-ABC. The goal of the second experiment is again to demonstrate computational savings but also illustrate the differences in posterior distributions found with and without the output independence assumption. Finally, the goal of the third experiment is to demonstrate the potential use of the posterior predictive distribution, along with the parameter posterior distribution, for model selection and comparison in same the spirit as [10].

For ASL-ABC, at each MH-step, there are initial simulations run at both and (akin to the marginal ABC sampler). Then, while the unconditional MH error is greater than , additional simulations are run. Some experiments will vary or or both. By default the full output covariance will be estimated by ASL-ABC, except for one part of experiment 2. For GPS-ABC, the meaning of is different; it represents the initial number of samples from the prior used to train the GPs. From that point onward, the GP runs one simulation at a time (while the error is greater than ). Although both algorithms permit , we have set for all experiments involving ASL-ABC and GPS-ABC.

### 5.1 Toy Problem

As an illustrative example we apply both ASL-ABC and GPS-ABC to estimate the posterior of the rate of an exponential distribution based on the mean statistic from observations from the true generating model with a Gamma prior over the rate parameter; this is the exponential example in [27]. Let be a vector of draws from an exponential distribution with rate , ; then let be the observation statistic. Our simulator is a two-step process: generate exponential draws with rate ; then return their average. I.e. the same generating process as the observations . With a Gamma prior , the goal of this experiment is to infer using ABC methods. The prior distribution has parameters ; a very broad prior over . In this controlled experiment, the same random seed was used for all MCMC runs, which fixes and -MAP (not quite and since they are random draws).

In Figure 2 we show the posterior distributions produced from different algorithms: Kernel-ABC, ASL-ABC, and GPS-ABC. Kernel-ABC uses the marginal-likelihood MCMC sampler (see Section 2.1, with ), using a Gaussian kernel, , and . At each MH step, values for are simulated at both and . The true posterior distribution is shown as a dashed line. All the runs produced K samples, the first K of which were discarded. All algorithms were initialized at (very far from the posterior mode). For ASL-ABC we ran two samplers using , , and ; and for GPS-ABC we initialized the GP with points from the prior, then ran the sampler with . All algorithms used the same proposal for , a normal with .

Model | ||||
---|---|---|---|---|

ASL-ABC | 572 | 135 | 100 | |

ASL-ABC | 2197 | 2060 | 2000 | |

GPS-ABC | 1.3 | 0.184 | 0.029 | |

Kernel-ABC | N/A | 20 | 20 | 20 |

The number of simulations used by the algorithms is shown in Table 1. The results for Kernel-ABC can be deceptive. Although Kernel-ABC performs quite well on this toy example (its over-dispersed results for small can be corrected using the pseudo-marginal sampler at the cost of a higher rejection rate), its performance should be taken with a grain of salt: Kernel-ABC has an extremely high rejection rate ( for ) and therefore mixes very slowly (as noted in Figure 2, the sampler was run longer to provide a smoother histogram). For this reason, though not pursued in this paper, a more useful measure of performance may be the rate of independent samples per simulation call. For the ASL-ABC the results are good even for large . Even though the distributions are more accurate when (not shown in figure), for , using an adaptive number of simulations produces many fewer simulations than for . For GPS-ABC the results show a significant advantage over ASL-ABC and Kernel-ABC. For the most accurate run, where , only simulations were required (compared to over for ASL-ABC). With , nearly the same accuracy is achieved with only simulations (compared to for ASL-ABC). The final GP models are shown in Figure 3, zoomed in to the region near the posterior mode (shown as dashed vertical). The mean and total uncertainty (2 stds) are shown with interior (yellow) and exterior (blue) shading, respectively; training points are shown as (red) circles. Notice how the GP model uncertainty contracts as decreases; this an effect produced by adding training points at uncertain locations.

### 5.2 Chaotic Ecological Systems

Adult blowfly populations exhibit dynamic behavior for which several competing population models exist. In this experiment, we use observational data and a simulation model from Wood [28], based on their improvement upon previous population dynamics’ theory. Population dynamics are modeled using (discretized) differential equations that can produce chaotic behavior for some parameter settings. We use model 4 from [28] as our simulator; the generated time-series are then processed into statistics explained below. An example blowfly replicate series is shown in Figure 4, along with several times-series generated from .

The population dynamics’ equation generates using the following update rule:

where and are sources of noise, and is an integer (not to be confused with the used as the MH acceptance threshold in our algorithms). In total, there are 6 parameters . See [28] for further details about the significance of the parameters.

We model the likelihood of statistics: 1) the mean of , 2) the mean minus the median, 3) the maximal peaks in the time-series with a smooth detection threshold, and 4) the log of the maximum of . Our first two statistics are used by [28], but they have 21 other statistics we do not use. The reasons we do not choose all are 1) we did not wish to inaccurately replicate the statistics using the code from [28], 2) our goal is sampling given a set of statistics (not finding the optimal statistics), and 3) time-series simulated with parameters based on our statistics are a reasonable representation of the observations (see Figure 4).

Like [28] we model the log distribution of , with the following broad (in the non-log space) : , , , , , . Time-series generated with parameters from this prior distribution produce extremely varied results, some are chaotic, some are degenerate, etc. Modeling this simulator is very challenging.

We compared the results from using ASL-ABC with , , and against GPS-ABC with , , and . (Recall that generating samples for GPS-ABC occurs only once at initialization, while for ASL-ABC twice at every MH-step.) We set for GPS-ABC to control the number of simulations at each time step. For both algorithms the proposal distribution was a zero-mean Gaussian with standard deviation equal to the prior’s value, with only a single dimension proposed at each MH step.

Results are shown in the form of scatter plots for in Figures 5 and 6. Note we are displaying the parameters in -space, not log-. Results for the ASL-ABC show interesting non-linear relationships between parameters. On close inspection of the results it can be seen that the posterior modes for the two models are very similar. However, the variance of the ASL-ABC is larger than for GPS-ABC. This is most likely due to the independence assumption made by GPS-ABC over the output statistics. We re-ran ASL-ABC with covariance estimators forced to be diagonal, and the results were very similar to GPS-ABC (see Figure 7). On the one hand the difference between full and diagonal covariance can be interpreted as overconfidence in the posterior, and on the other it can can be viewed as a robust approximation to the full posterior, because it is very difficult to determine, for the ASL-ABC, the degree to which the covariance estimator is influenced by outliers (and overfits, in a sense). By enforcing diagonal likelihoods for the ASL-ABC, similar variance in the posterior distributions are found. Despite the potential lack of a fully informative posterior, GPS-ABC results were achieved with only a few thousand simulation calls, whereas ASL-ABC used well over a million calls. If the cost per simulation is very expensive, the trade-off of accuracy versus computation leans heavily towards GPS-ABC.

### 5.3 Feedback Wright-Fisher model

The Wright-Fisher model is a standard stochastic model for reproduction in population genetics. In two experiments, we study the potential of our methods for model selection and comparison (rather than computational savings) using the ASL-ABC (similar results can be obtained using GPS-ABC). We apply our methods to the model in [30], where there are three types of cell populations: stem cells (SC), transit-amplifying cells (TAC), and terminal differentiated cells (TDC). A simulation begins with an initial cell population that is updated using stochastic equations having parameters . The true parameter setting is . A simulation consists of running the equations for time-steps, repeated times. Statistics are the mean and log-variance of these population sizes. Due to space considerations, only the results for the mean size are shown in the figures.

For this experiment we have slightly changed notation for the statistics. Here we call the true statistics (whereas previously we simply referred to these as the observed statistics ) and as the posterior predictive distribution defined as

(20) |

where is the target distribution of our MCMC sampler (, previously). Therefore the true statistics is the result of a single (arbitrarily chosen) simulation using . We set broad priors over .

In Figure 8, we show results of experiment 1. In this experiment, the number of parameters varies, allowing us to compare as the number of parameters that are inferred changes. From top to bottom, we infer , , and , respectively, and set any missing parameters to their corresponding value in . The solid curve is (the distribution of statistics generated at ). The true is a draw from this distribution (shown as vertical line). When the algorithm infers all three parameters, the posterior predictive distribution is centered at , whereas when it infers only a single parameter, it is shifted towards . If we inspect , we can see that the value for is the same for all three rows. For the first row, this can only be achieved by having a larger variance in the posterior predictive distribution; for the last row, this is achieved by decreasing the variance.

In experiment 2, we compare the result for the middle row in Figure 8 where is set to the true or “good” setting (), with set to a “bad” (). In Figure 9 top row, we show posterior predictive distribution for the latter case; notice that these posterior predictive distributions show evidence of large bias and large variance. The bottom row compares the posterior distributions for both cases. Here the bias is less pronounced, but the variance remains large.

The results of these two experiments demonstrate the potential of our approach for model checking, selection, and improvement described by Gelman et al. [10] (chapter 6). This type of analysis comes naturally for our algorithms for two reasons: the ABC statistics are natural choices for Gelman et al.’s model checking statistics and these statistics are pre-computed as part of the MH-step of our algorithms. We can now answer questions like: what is the likelihood of a set of observations using simulator versus simulator ?

## 6 Discussion and Future Work

We have presented a promising framework for performing inference in expensive, simulation-based models. Our algorithms improve current state-of-the-art ABC methods in that they require many fewer calls to the simulator.

Using GPs for surrogate modeling has an appealing elegance; as nonparametric Bayesian models, they naturally incorporate both model and pseudo-data uncertainty into the MCMC algorithms. However, there are several technical issues and modeling limitations with GPs used for surrogate modeling. Heteroskedatic noise is more likely the norm than the exception for complicated simulators. The blowfly simulator is a prime example of this. Improvements to our GPs may be achieved using an input-dependent noise model [11, 13], where the noise is an additional independent Gaussian process. Another limitation of our GP model is the output independence assumption. A more realistic assumption is a full covariance Gaussian process such as the convolution processes of [12, 7, 1]. The other major difficulty is setting the GP hyper-parameters and updating them during the ABC run.

Our GPS-ABC uses a random-walk proposal distribution which is inefficient for exploring the target distribution. Using GPs offers the opportunity to use other techniques to improve the mixing (and in turn computational cost). For example, in [17] a Hamiltonian Monte Carlo run on the GP surface is used to generate independent proposals. If applied to ABC, their algorithm would require the equivalent of a full simulation at the proposed location, whereas if we incorporated a similar technique, we would then test the GP uncertainty to determine if a simulation was required.

There has been a recent surge in interest in Bayesian optimization using Gaussian process (GP-BO) surrogates of the objective surface [8, 24]. GP-BO is often applied to problems where simulation or sometimes user feedback guides the surrogate’s construction. What is most interesting about GP-BO is its use of model uncertainty to actively determine the next simulation location implemented through acquisition functions. These ideas can probably be generalized to our GPS-ABC algorithm to further reduce simulation costs, while at the same time maintaining control of the MCMC error. The authors of [24] also develop distributed algorithms for optimization that may serve as inspiration for parallel ABC sampling algorithms.

## References

- Álvarez & Lawrence [2011] Álvarez, M. and Lawrence, L.D. Computationally efficient convolved multiple output gaussian processes. Journal of Machine Learning Research, 12:1425–1466, 2011.
- Andrieu & Roberts [2009] Andrieu, C. and Roberts, G. The pseudo-marginal approach for efficient monte carlo computations. The Annals of Statistics, 37(2):697–725, 2009.
- Bazin et al. [2010] Bazin, Eric, Dawson, Kevin J, and Beaumont, Mark A. Likelihood-free inference of population structure and local adaptation in a bayesian hierarchical model. Genetics, 185(2):587–602, 2010.
- Beaumont et al. [2002] Beaumont, Mark A, Zhang, Wenyang, and Balding, David J. Approximate bayesian computation in population genetics. Genetics, 162(4):2025–2035, 2002.
- Beaumont et al. [2009] Beaumont, Mark A, Cornuet, Jean-Marie, Marin, Jean-Michel, and Robert, Christian P. Adaptive approximate bayesian computation. Biometrika, 96(4):983–990, 2009.
- Bilionis et al. [2013] Bilionis, I., Zabaras, N., Konomi, B.A., and Lin, G. Multi-output separable gaussian process: Towards an efficient fully bayesian paradigm for uncertainty quantification. Journal of Computational Physics, 241:212–239, 2013.
- Boyle & Frean [2005] Boyle, P. and Frean, M. Dependent gaussian processes. Advances in Neural Information Processing Systems 17, 2005.
- Brochu et al. [2010] Brochu, E., Cora, Vlad M., , and de Freitas, Nando. A tutorial on bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning. Technical report, UBC Technical Report, 2010.
- Ceperley & Dewing [1999] Ceperley, D.M. and Dewing, M. The penalty method for random walks with uncertain energies. Journal of Chemical Physics, 110(20):9812–9820, 1999.
- Gelman et al. [2004] Gelman, A., Carlin, J.B., Stern, H.S., and Rubin, D.B. Bayesian Data Analysis. Chapman and Hall/CRC, New York, second edition, 2004.
- Goldberg et al. [1998] Goldberg, P.W., Williams, C.K.I., and Bishop, C.M. Regression with input-dependent noise: A gaussian process treatment. Advances in Neural Information Processing Systems 10, 1998.
- Higdon [2002] Higdon, D. Space and space-time modeling using process convolutions. In Anderson, C.W., Barnett, V., Chatwin, P.C, and Abdel, E.-S.H. (eds.), Quantitative Methods for Current Environmental Issues, pp. 37–56. Springer London, 2002. ISBN 978-1-4471-1171-9.
- Kersting et al. [2007] Kersting, Kristian, Plagemann, Christian, Pfaff, Patrick, and Burgard, Wolfram. Most likely heteroscedastic gaussian process regression. In Proceedings of the 24th international conference on Machine learning, pp. 393–400. ACM, 2007.
- Korattikara et al. [2014] Korattikara, A., Chen, Y., and Welling, M. Austerity in mcmc land: Cutting the metropolis-hastings budget. ICML, 2014. To appear.
- Marin et al. [2012] Marin, J.-M., Pudlo, P., Robert, C.P., and Ryder, R.J. Approximate bayesian computational methods. Statistics and Computing, 22:1167–1180, 2012.
- Nicholls et al. [2012] Nicholls, G.K., Fox, C., and Watt, A.M. Coupled mcmc with a randomized acceptance probability. Technical report, arXiv:1205.6857v1, 2012.
- Rasmussen [2003] Rasmussen, C.E. Gaussian processes to speed up hybrid monte carlo for expensive bayesian integrals. Bayesian Statistics, 7:651–659, 2003.
- Ratmann et al. [2007] Ratmann, O., Jorgensen, O., Hinkley, T., Stumpf, M., Richardson, S., and Wiuf, C. Using likelihood-free inference to compare evolutionary dynamics of the protein networks of h. pylori and p. falciparum. PLoS Computational Biology, 3(11):e230, 2007.
- Roberts & Rosenthal [2007] Roberts, G.O. and Rosenthal, J.S. Coupling and ergodicity of adaptive mcmc. J. Appl. Prob., 44:458–475, 2007.
- Schafer & Freeman [2012] Schafer, C.M. and Freeman, P.E. Likelihood-free inference in cosmology: Potential for the estimation of luminosity functions. In Statistical Challenges in Modern Astronomy V, pp. 3–19. Springer, 2012.
- Sisson et al. [2007] Sisson, SA, Fan, Y, and Tanaka, Mark M. Sequential monte carlo without likelihoods. Proceedings of the National Academy of Sciences, 104(6):1760, 2007.
- Sisson et al. [2010] Sisson, S.A., Peters, G.W., Briers, M., and Fan, Y. A note on target distribution ambiguity of likelihood-free samplers. Arxiv preprint arXiv:1005.5201v1 [stat.CO], 2010.
- Sisson & Fan [2010] Sisson, Scott A and Fan, Yanan. Likelihood-free markov chain monte carlo. Arxiv preprint arXiv:1001.2058, 2010.
- Snoek et al. [2012] Snoek, Jasper, Larochelle, Hugo, and Adams, Ryan Prescott. Practical bayesian optimization of machine learning algorithms. Advances in Neural Information Processing Systems 25, 2012.
- Tavare et al. [1997] Tavare, S., Balding, D.J., Griffiths, R.C., and Donnelly, P. Inferring coalescence times from dna sequence data. Genetics, 145(2):505–518, 1997.
- Toni et al. [2009] Toni, Tina, Welch, David, Strelkowa, Natalja, Ipsen, Andreas, and Stumpf, Michael PH. Approximate bayesian computation scheme for parameter inference and model selection in dynamical systems. Journal of the Royal Society Interface, 6(31):187–202, 2009.
- Turner & Van Zandt [2012] Turner, Brandon M and Van Zandt, Trisha. A tutorial on approximate bayesian computation. Journal of Mathematical Psychology, 56(2):69–85, 2012.
- Wood [2010] Wood, Simon N. Statistical inference for noisy nonlinear ecological dynamic systems. Nature, 466(7310):1102–1104, 2010.
- Xue & Titterington [2011] Xue, Jing-Hao and Titterington, D. Michael. The p-folded cumulative distribution function and the mean absolute deviation from the p-quantile. Statistics and Probability Letters, 81:1179–1182, 2011.
- Zhu et al. [2012] Zhu, X., Welling, M., Jin, F., and Lowengrub, J.S. Predicting simulation parameters of biological systems using a gaussian process model. Statistical Analysis and Data Mining, 5(6):509–522, 2012.