Rapid Bayesian inference for expensive stochastic models

Bayesian analysis techniques are powerful tools for the quantification of uncertainty in parameters, models and predictions (Gelman et al., 2014). Unfortunately, for many stochastic models of practical interest, the likelihood function is intractable. ABC methods replace likelihood evaluation with an approximation based on stochastic simulations of the proposed model, this is captured directly in ABC rejection sampling (Pritchard et al., 1999; Tavaré et al., 1997) where $\mathcal{M}$ samples are generated from an approximate posterior, denoted by $p(\boldsymbol{\theta}\mid\rho(\mathcal{D},\mathcal{D}_{s})\leq\epsilon)$. Here $\mathcal{D}_{s}\sim s(\mathcal{D}\mid\boldsymbol{\theta})$ is a data generation process based on simulation of the model, $\rho(\mathcal{D},\mathcal{D}_{s})$ is a discrepancy metric, and $\epsilon$ is the discrepancy threshold. The resulting accepted parameter samples are distributed according to $p(\boldsymbol{\theta}\mid\rho(\mathcal{D},\mathcal{D}_{s})\leq\epsilon)\to p(\boldsymbol{\theta}\mid\mathcal{D})$ as $\epsilon\to 0$.

The average acceptance probability of a proposed parameter sample $\boldsymbol{\theta}^{*}$ is $\mathcal{O}(\epsilon^{d})$ (Fearnhead and Prangle, 2012), where $d$ is the dimensionality of the data space, $\mathbb{D}$. This renders rejection sampling computationally expensive or even completely prohibitive, especially for high-dimensional parameter spaces (Marjoram et al., 2003; Sisson et al., 2007). Summary statistics can reduce the data dimensionality, however, they will often incur information loss (Barnes et al., 2012; Blum et al., 2013; Fearnhead and Prangle, 2012). However, strategies including regression adjustment and marginal adjustment strategies can improve the accuracy of dimension reductions (Beaumont et al., 2002; Nott et al., 2014).

In the SMC-ABC method, importance resampling is applied to a sequence of $R$ ABC posteriors with discrepancy thresholds $\epsilon_{1}>\cdots>\epsilon_{R}$, with $\epsilon_{R}$ indicating the target ABC posterior. Given $\mathcal{M}$ weighted samples $\{(\boldsymbol{\theta}^{i},w^{i})\}_{i=1}^{\mathcal{M}}$, called particles, from the prior $p(\boldsymbol{\theta})$, particles are filtered through each ABC posterior using three main steps for each particle $\boldsymbol{\theta}^{i}$: (i) the particle is perturbed using a proposal kernel density $q_{r}(\boldsymbol{\theta}\mid\boldsymbol{\theta}^{i})$; (ii) an accept/reject step is performed; and (iii) importance weights are updated. Once all particles have been updated and reweighted, resampling of particles is performed to avoid particle degeneracy. For reference, the SMC-ABC algorithm as initially developed by Sisson et al. (2007) and Toni et al. (2009) is given in Algorithm 1. The number of particles, $\mathcal{M}$, and the number of intermediate distributions, $R$, influence the accuracy and performance, respectively, of the sampler. Setting $\mathcal{M}$ too small can lead to large estimator variability and particle degeneracy, and setting $R$ too small leads to large divergence between successive distributions that can result in high rejection rates.

Note that throughout all algorithms used in this manuscript, we assume that the initial set of weighted particles, $\{\boldsymbol{\theta}_{0}^{i},w_{0}^{i}\}_{i=1}^{\mathcal{M}}$, are independent, identically distributed samples from the prior, $p(\boldsymbol{\theta})$, and therefore have uniform weight, $w_{0}^{i}=1/\mathcal{M}$, for all $i=1,\ldots,\mathcal{M}$. However, the methods are general enough to deal with prior distributions that require importance sampling to draw truly weighted particles.

For a fixed choice of $R$, efficient use of SMC-ABC depends critically on the selection of appropriate proposal kernels and threshold sequences. The process of sampling at the target threshold, $\epsilon_{r}$, given the weights of the previous threshold, $\epsilon_{r-1}$, is described by Filippi et al. (2013)

where the data space, $\mathbb{D}$, has dimensionality $d$, $B(\mathcal{D},\epsilon_{r})$ is a $d$-dimensional ball centered on the data with radius $\epsilon_{r}$, and $\mathds{1}_{A}\left(x\right)$ denotes the indicator function with $\mathds{1}_{A}\left(x\right)=1$ if $x\in A$, otherwise $\mathds{1}_{A}\left(x\right)=0$. The normalization constant, $a_{r,r-1}$, can be interpreted as the average acceptance probability across all particles. We see this by noting that Equation (1) can be reduced to

Here the distribution $\eta_{r-1}(\boldsymbol{\theta}_{r})$ represent the proposal mechanism,

and $\mathbb{P}(\mathcal{D}_{s}\in B(\mathcal{D},\epsilon_{r})\mid\boldsymbol{\theta}_{r})$ is the probability that simulated data is within $\epsilon_{r}$ of the data $\mathcal{D}$ given a parameter value $\boldsymbol{\theta}_{r}$. Therefore, the normalizing constant is

that is, $a_{r,r-1}$ is the average acceptance probability.

From a computational perspective, the goal is to choose $\eta_{r-1}(\boldsymbol{\theta}_{r})$ to maximize $a_{r,r-1}$. However, this would not necessarily result in a $\xi_{r}\left(\boldsymbol{\theta}_{r}\mid\mathcal{D}\right)$ that is an accurate approximation to the true target ABC posterior $p(\boldsymbol{\theta}_{r}\mid\rho(\mathcal{D},\mathcal{D}_{s})\leq\epsilon_{r})$. To achieve this goal, we require $\eta_{r-1}(\boldsymbol{\theta}_{r})$ such that the Kullback-Leibler divergence (Kullback and Leibler, 1951), $D_{\text{KL}}\left(\xi_{r}\left(\cdot\mid\mathcal{D}\right);p(\cdot\mid\rho(\mathcal{D},\mathcal{D}_{s})\leq\epsilon_{r})\right)$, is minimized. Beaumont et al. (2009) and Filippi et al. (2013) demonstrate how the latter goal provides insight into how to optimally choose $\eta_{r-1}(\boldsymbol{\theta}_{r})$. The key is to note that $D_{\text{KL}}\left(\xi_{r}\left(\cdot\mid\mathcal{D}\right);p(\cdot\mid\rho(\mathcal{D},\mathcal{D}_{s})\leq\epsilon_{r})\right)$ can be decomposed as follows,

where $E(\boldsymbol{\theta}_{r})=\mathbb{E}\left[\log_{e}\left(\int_{\mathbb{D}}\mathds{1}_{B(\mathcal{D},\epsilon_{r})}\left(\mathcal{D}_{s}\right)s(\mathcal{D}_{s}\mid\boldsymbol{\theta}_{r})\,\text{d}\mathcal{D}_{s}\right)\right]$ is independent of $\eta_{r-1}(\boldsymbol{\theta}_{r})$. By rearranging Equation (5), we obtain

That is, minimizing $D_{\text{KL}}\left(\eta_{r-1}\left(\cdot\right);p(\cdot\mid\rho(\mathcal{D},\mathcal{D}_{s})\leq\epsilon_{r})\right)$ is equivalent to minimizing
$D_{\text{KL}}\left(\xi_{r}\left(\cdot\mid\mathcal{D}\right);p(\cdot\mid\rho(\mathcal{D},\mathcal{D}_{s})\leq\epsilon_{r})\right)$ and maximizing $a_{r,r-1}$ simultaneously. Therefore, any proposal mechanism that is closer, in the Kullback-Liebler sense, to $p(\cdot\mid\rho(\mathcal{D},\mathcal{D}_{s})\leq\epsilon_{r})$ is more efficient.

We apply the optimal adaptive scheme of Beaumont et al. (2009) and Filippi et al. (2013) for multivariate Gaussian proposals. That is, we set

where $n$ is the dimensionality of parameter space, $\boldsymbol{\Theta}$, and

Consider a fixed sequence of ABC posteriors for the stochastic model inference problem, $\{p(\boldsymbol{\theta}\mid\rho(\mathcal{D},\mathcal{D}_{s})\leq\epsilon_{r})\}_{r=1}^{R}$. We want to apply SMC-ABC (Supplementary Material) to efficiently sample from this sequence with adaptive proposal kernels, $\{q_{r}(\boldsymbol{\theta}^{*}\mid\boldsymbol{\theta})\}_{r=1}^{R}$ (Beaumont et al., 2009; Filippi et al., 2013). Our method exploits an approximate model to further improve the average acceptance probability.

The PC-SMC-ABC method is a promising modification to SMC-ABC that can accelerate inference for expensive stochastic models without introducing bias. However, other approaches can be used to obtain further computational improvements. Here, we consider an alternate approach to utilizing approximate models that aims to get the most out of a small sample of expensive stochastic simulations. Unlike PC-SMC-ABC, this method is generally biased, but it has the advantage of yielding a small and fixed computational budget. Specifically, we define a parameter $\alpha\in[0,1]$, such that $1/\alpha$ is the target computational speedup, for example, $\alpha=1/10$ should result in approximate $10$ times speedup. We apply the SMC-ABC method using $\tilde{\mathcal{M}}=\lfloor(1-\alpha)\mathcal{M}\rfloor$ particles based on the approximate model, and then use $\hat{\mathcal{M}}=\lceil\alpha\mathcal{M}\rceil$ particles based on the stochastic model to construct a hybrid population of $\mathcal{M}=\hat{\mathcal{M}}+\tilde{\mathcal{M}}$ particles that will represent the final inference on the stochastic model. The key idea is that we use the $\hat{\mathcal{M}}$ particles of the expensive stochastic model to inform a transformation on the $\tilde{\mathcal{M}}$ particles of the approximation such that they the emulate particles of expensive stochastic model. Here, $\lfloor\cdot\rfloor$ and $\lceil\cdot\rceil$ are, respectively, the floor and ceiling functions.

The Ornstein–Uhlenbeck process is a mean reverting stochastic process with many applications in finance, biology and physics (Uhlenbeck and Ornstein, 1930). The evolution of the continuous state $X_{t}$ is given by an Itō stochastic differential equation (SDE) of the form

where $\mu$ is the long-term mean, $\gamma$ is the rate of mean reversion, $\sigma$ is the process volatility, $W_{t}$ is a Wiener process and $x_{0}$ is a constant initial condition. Example realisations are shown in Figure 2(A).

We consider data consisting of $N$ independent realisations of the Ornstein–Uhlenbeck processs at time $T<\infty$, that is, $\mathcal{D}=[X_{T}^{1},X_{T}^{2},\ldots,X_{T}^{N}]$. This inference problem is analytically tractable since the Fokker–Planck equation can be solved to obtain a Gaussian distribution for the data,

For demonstration purposes we will assume a solution for the full Fokker–Planck equation is unavailable and perform ABC inference to estimate the volatility parameter, $\sigma$, using stochastic simulation with the Euler–Maruyama discretisation (Maruyama, 1955),

where $\xi_{t}$ is a standard normal variate and $\Delta t$ is a small time step. For the approximate model we take the stationary distribution of the Ornstein–Uhlenbeck process obtained by taking $T\to\infty$ and solving the steady state of the Fokker–Planck equation,

This kind of approximation will often be possible since the steady state Fokker–Planck equation is more likely to be tractable than the transient solution for most SDE models. As shown in Figure 2(B), the stationary solution is a better approximation for the true variance than for the true mean. Therefore, for small $T$ this approximation would be more appropriate as a preconditioner for inference of the volatility parameter, $\sigma$, in Equation (11) than for the long-time mean, $\mu$. In general, the performance expected from preconditioning will increase as $T$ increases since the stationary approximation will become more accurate (Supplementary Material).

Figure 3 demonstrates th results of applying PC-ABC-SMC (Algorithm 2) with $\mathcal{M}=1000$ particles, showing intermediate distributions, for inferring $D=\sigma^{2}/2$ given data at $T=1$, using stochastic simulation for the exact model (Equation (13)) with $\Delta t=0.01$, and the stationary distribution (Equation (14)) for the approximate model used in the preconditioning step.

In this example, all other parameters are treated as known with $x_{0}=10$, $\mu=1$ and $\gamma=2$. Data is generated using $D=10$ as in Figure 2. In each step, the preconditioner distribution for threshold $\epsilon_{r}$ (dashed line) is a better proposal distribution for the target (solid line) than that of the exact model at threshold $\epsilon_{r-1}$ (dotted line). The overall speedup factor is approximately $1.5$ for this example, and it continues to improve as $T$ increases since Equation (14) becomes an increasingly better approximation to Equation (12) (See Supplementary Material). For truly intractable problems, such as the lattice-based random walk models presented in Sections 3.4 and 3.5, we obtain superior performance gains of up to a factor of four.

The stochastic discrete model we consider is a lattice-based random walk model that is often used to describe populations of motile cells (Jin et al., 2016). The model involves initially placing a population of $N$ agents of size $\delta$ on a lattice, $L$ (Callaghan et al., 2006; Simpson et al., 2010), for example an $I\times J$ hexagonal lattice (Jin et al., 2016). This hexagonal lattice is defined by a set of indices
$L=\{(i,j):i\in[0,1,\ldots,I-1],j\in[0,1,\ldots,J-1]\}$, and a neighborhood function,

Lattice indices are mapped to Cartesian coordinates using

We define an occupancy function such that $C(\ell,t)=1$ if site $\ell$ is occupied by an agent at time $t\geq 0$, otherwise $C(\ell,t)=0$. This means that in our discrete model each lattice site can be occupied by, at most, one agent.

During each discrete time step of duration $\tau$, agents attempt to move with probability $P_{m}\in[0,1]$ and attempt to proliferate with probability $P_{p}\in[0,1]$. If an agent at site $\ell$ attempts a motility event, then a neighboring site will be selected uniformly at random. The motility event is aborted if the selected site is occupied, otherwise the agent will move to the selected site (Figure 4(A)–(C)).

For proliferation events, the local neighborhood average occupancy,

is calculated and a uniform random number $u\sim\mathcal{U}(0,1)$ is drawn. If $u>f(\hat{C}(\ell,t))$, where $f(\hat{C}(\ell,t))\in[0,1]$ is called the crowding function (Browning et al., 2017; Jin et al., 2016), then the proliferation event is aborted due to local crowding effects and contact inhibition. If $u\leq f(\hat{C}(\ell,t))$, then proliferation is successful and a daughter agent is placed at a randomly chosen unoccupied lattice site in $\mathcal{N}(\ell)$ (Figure 4(A),(D)–(E)). The evolution of the model is generated through repeating this process though $M$ time steps, $t_{1}=\tau,\,t_{2}=2\tau,\ldots,$ $t_{M}=M\tau$. This approach, based on the work by Jin et al. (2016), supports a generic proliferation mechanism since $f(\hat{C}(\ell,t))$ is an arbitrary smooth function satisfying $f(0)=1$ and $f(K)=0$, where $K>0$ is the carrying capacity density. However, in the literature there are also examples that include other mechanisms such as cell-cell adhesion (Johnston et al., 2013), directed motility (Binny et al., 2016), and Allee effects (Böttger et al., 2015).

Discrete models do not generally lend themselves to analytical methods, consequently, their application is intrinsically tied to computationally intensive stochastic simulations and Monte Carlo methods (Jin et al., 2016). As a result, it is common practice to approximate mean behavior using differential equations by invoking mean-field assumptions, that is, to treat the occupancy status of lattice sites as independent (Callaghan et al., 2006; Simpson et al., 2010). The resulting approximate continuum-limit descriptions (Supplementary Material) are partial differential equations (PDEs) of the form

where $\mathcal{C}(x,y,t)=\mathbb{E}\left[C(\ell,t)\right]$, $D=\lim_{\delta\to 0,\tau\to 0}P_{m}\delta^{2}/(4\tau)$ is the diffusivity, $\lambda=\lim_{\tau\to 0}P_{p}/\tau$ is the proliferation rate with $P_{p}=\mathcal{O}(\tau)$, and $f(\cdot)$ is the crowding function that is related to the proliferation mechanism implemented in the discrete model (Browning et al., 2017; Jin et al., 2016). For spatially uniform populations there will be no macroscopic spatial gradients on average, that is $\nabla\mathcal{C}(x,y,t)=\mathbf{0}$. Thus, $\mathcal{C}(x,y,t)$ is just a function of time, $\mathcal{C}(t)$, and the continuum limit reduces to an ordinary differential equation (ODE) describing the net population growth,

For many standard discrete models, the crowding function is implicitly $f(\mathcal{C})=1-\mathcal{C}$ (Callaghan et al., 2006). That is, the continuum limits in Equation (16) and Equation (17) yield the Fisher-KPP model (Edelstein-Keshet, 2005; Murray, 2002) and the logistic growth model (Tsoularis and Wallace, 2002; Warne et al., 2017), respectively. However, non-logistic growth models , for example, $f(\mathcal{C})=(1-\mathcal{C})^{n}$ for $n>1$, have also been considered (Jin et al., 2016; Simpson et al., 2010; Tsoularis and Wallace, 2002).

The Allee effect refers to the reduction in growth rate of a population at low densities. This is particularly well studied in ecology where there are many mechanisms that give rise to this phenomenon (Taylor and Hastings, 2005; Johnston et al., 2017). We introduce an Allee effect into our discrete model by choosing a crowding function of the form

where $\hat{C}(\ell,t)\in[0,1]$ is the local density at the lattice site $\ell\in L$, at time $t$, $K>0$ is the carrying capacity density, and $A$ is the Allee parameter which yields a weak Allee effect for $A\geq 0$ (Wang et al., 2019). Note that smaller values of $A$ entail a more pronounced Allee effect with $A<0$ leading to a strong Allee effect that can lead to species extinction (Wang et al., 2019). For simplicity, we only consider the weak Allee effect here, but our methods are general enough to consider any sufficiently smooth $f(\cdot)$.

Studies in ecology often involve population counts of a particular species over time (Taylor and Hastings, 2005). In the discrete model, the initial occupancy of each lattice site is independent, and hence there are no macroscopic spatial gradients on average. It is reasonable to summarize simulations of the discrete model at time $t$ by the average occupancy over the entire lattice, $\bar{C}(t)=(1/IJ)\sum_{\ell\in L}C(\ell,t)$. Therefore, the continuum limit for this case is given by (Wang et al., 2019)