Likelihood free inference for Markov processes: a comparison

It has been well documented that the efficiency of random walk Metropolis algorithms is highly dependent on the choice of proposal kernel. A distribution which yields small deviations from the current state will ensure that a large number of moves are accepted but samples will be highly correlated. Large moves around the space on the other hand will often be rejected leading to the chain spending large amounts of time stuck at the same value. Under various assumptions about the target it has been shown that the optimal scaling for a Gaussian proposal kernel is

where $\Sigma$ is the covariance of the posterior distribution and $d$ is the number of parameters being estimated, (Roberts et al., 1997; Roberts & Rosenthal, 2001). The starting parameter vector, $\theta_{0}$, of the chain also has an effect on the efficiency of the sampler. A choice of $\theta_{0}$ which is far from a region of non–negligible posterior density will lead to a chain which takes a long time to move toward the target distribution, whereas a chain initialised close to stationarity will yield useful samples sooner. This burn–in period can sometimes consume a sizeable fraction of the computational budget. Particle MCMC as described in section 3.1 also relies on a sequential Monte Carlo algorithm for approximation of the likelihood, $\hat{\pi}(\mathcal{D}|\theta)$. The bootstrap filter requires multiple model realisations, via a set of “particles” in order to achieve this approximation. In addition the approximation has to be calculated at every iteration of the MCMC algorithm, hence clearly the number of particles in the filter will greatly affect the runtime of the resultant algorithm. A small number of particles will result in a shorter computation time for the likelihood approximation but leads to larger variability in the estimated likelihood. This increased variability leads to decreased efficiency of the inference scheme, as noted by Andrieu & Roberts (2009). A large number of particles, useful for consistent estimates of $\hat{\pi}(\mathcal{D}|\theta)$ will lead to slower posterior sampling in the chain. An optimal choice for the number of particles has been explored by Pitt et al. (2012) who suggest that the number of particles should be chosen such that the variance of the log–likelihood estimates is around 1. However Doucet et al. (2012) show that the efficiency of the scheme is good for variances between 0.25 and 2.25.

In practice, for the purpose of this comparison we choose an initial parameter vector, $\theta_{0}$, as a random sample from the posterior distribution. The number of particles used in the particle filter, $N$ is then chosen by repeated runs of a particle filter with increasing $N$ until $1.5<\operatorname{Var}(\hat{l}(\theta_{0}|\mathcal{D}))<1.8$. We then use the covariance matrix of the posterior, $\Sigma_{p}$ to inform our choice for $\Sigma_{q}$, the Gaussian random walk proposal variance,

During our first experiments with the pMCMC algorithm for these models we approached initialisation and tuning of the algorithm under the assumption of no knowledge of the posterior distribution of interest. However, this proved to be problematic as finding a sensible choice of $\theta_{0}$, number of particles, $N$, and proposal variance, $\Sigma_{q}$, often used a large proportion of the allocated computational budget. Under the computational restrictions imposed by our budget choice this made pMCMC look completely uncompetitve. This problem itself is interesting as it highlights a potential drawback of using pMCMC in practice. We discuss this issue further in section 5.1.5.

Initialisation of a sequential ABC algorithm as described in section 3.2 is somewhat less involved. This is due to the fact that optimal Gaussian proposal kernels for advancement to subsequent targets can be calculated during execution. In addition the sequence of tolerances is chosen adaptively throughout the algorithm. It remains that there is need to specify an initial tolerance value, $\epsilon_{0}$. One could argue that tuning the choice of metric and summary statistics to be used is also of interest. Discussion of how one might do this is beyond the scope of this article however, since we are limiting ourselves to the choice in equation 6 so as to ensure that as $\epsilon\rightarrow 0$ the resultant posterior approximation tends toward the true posterior distribution of interest, $\pi(\theta|\rho(\mathcal{D},\mathcal{D}^{\ast})<\epsilon)\rightarrow\pi(\theta|\mathcal{D})$. For an in depth discussion of methods in which to choose summary statistics see Blum et al. (2013). In order to choose a suitable $\epsilon_{0}$ for the scheme we simply calculate $\rho(\mathcal{D},\mathcal{D}^{\ast}|\theta)$ using a number of samples from $\pi(\theta)$. From this we take $\epsilon_{0}$ to be the value equivalent to the $1\%$-ile of the resultant distribution of distances.

The Lotka–Volterra predator prey system, Lotka (1925); Volterra (1926) is an example of a stochastic kinetic model that provides an ample starting point for investigation of parameter inference in models of this type. Although it is relatively simple, characterised by a set of 3 reactions on two species, it encompasses many of the difficulties associated with larger, more complex systems. Denoting the two species, prey, $X_{1}$ and predators, $X_{2}$, evolution of the system is governed by the following three reactions:

These reactions can be thought of as prey birth, a predator prey interaction resulting in the death of a prey and a predator birth, and predator death respectively. Reaction events are dependent on the current state of the system as well as the reaction rate parameters. Hence the trajectory of the evolution of species counts presents a Markov process on a discrete state space. This reaction network is summarised by its stoichiometry matrix , $S$ and hazard function $h(\mathbf{X},\theta)$:

Realisations from the model conditional on a vector of rate parameters, $\theta$, can be obtained exactly via algorithm 1, or approximated via a number of fast simulation algorithms.

For the purpose of making comparison we use a number of data sets over different observation regimes simulated using reaction rate vector $\theta=(1.0,0.005,0.6)$. In each case we corrupt the $X_{t}$ with a Gaussian error with mean 0 and variance $\sigma^{2}$, $\pi(d_{t}|X_{t},\sigma)\sim\mathcal{N}(X_{t},\sigma^{2})$. $X_{0}=(50,100)$ is used throughout. Plots of each of the data sets considered are in figure 1. Given this set of parameter values the model exhibits relatively stable oscillatory behavior for both species and provides an interesting starting point for our investigation. We shall use this model to explore posterior sampling efficiency given data sets of a range of sizes, under full and partial observation regimes, whilst also giving consideration to the effect of assuming known measurement or including this parameter in the set to be inferred. Data sets shown in figures 1b–d are denoted $\mathcal{D}^{1}$, $\mathcal{D}^{2}$ and $\mathcal{D}^{3}$ respectively. We introduce extra subscript notation such that $\mathcal{D}^{1}_{p}$ implies the data set $\mathcal{D}^{1}$ where predator observations have been discarded and $\mathcal{D}^{1}_{u}$ symbolises treatment of $\mathcal{D}^{1}$ under the assumption of unknown measurement error. In addition $\mathcal{D}^{1}_{\ast}$ will be used as a reference to the collection of data sets $\mathcal{D}^{1},\mathcal{D}^{1}_{u},\mathcal{D}^{1}_{p}$ and $\mathcal{D}^{1}_{u,p}$.

We now create a scenario in which prior parameter information is poor. We place uniform prior information on $\log(\theta)$,

and we place a Poisson prior distribution on the initial state

Where $\sigma^{2}$ is not assumed to be known we use

For each repeat we allow a computational budget of $10^{8}$ model realisations from the Direct method, algorithm 1. We choose this budget based on the fact that given the $\theta=(1.0,0.005,0.6)$ our simulator achieves $10^{4}$ simulations of length equivalent to $\mathcal{D}^{2}$ every 45-50 secs on our relatively fast Intel core i7-2600 clocked at 3.4 GHz. This yields an approximate total time spent simulating from the model of 14 hours plus some other comparably negligible computation costs for each individual inference run. Clearly improvement on simulation time can be made by parallelising the simulation of independent realisations from the direct method as well as other computational savings being made by clever optimisations in each algorithm. We have tried to disclude the effect of such algorithmic optimisation in our comparison as discussed in section 4. We include the information on approximate time here as a rough guide to practical implementation of inference for these types of models as well as the reasoning behind our particular budget choice.

Results for data sets $\mathcal{D}^{1}_{\ast}$ and $\mathcal{D}^{3}_{\ast}$ are contained as supplementary material. We report results for data set $\mathcal{D}^{2}_{\ast}$ , figure 1(c), shown in figure 2. Results in the supplementary material support those reported here, discrepancies between the two algorithms are reduced for $\mathcal{D}^{1}_{\ast}$ and exacerbated given the longer time series of $\mathcal{D}^{3}_{\ast}$. The plots in figure 2 show that for $\mathcal{D}^{2}$ there is a clear difference in the tails of the distribution. This feature is mirrored for $\mathcal{D}^{2}_{u}$. Additionally we see that, under the assumption of unknown measurement error ABC fails to identify the noise parameter, $\sigma$. Treatment of $\mathcal{D}^{2}_{p}$ yields a much smaller difference in the resultant posteriors between the two schemes however $\mathcal{D}^{2}_{u,p}$ reinforces the inability to infer $\sigma$ using the ABC SMC scheme. Each density plotted is the result of $10^{4}$ samples being collected. In the case of the ABC SMC this was ensured by retaining $10^{4}$ samples at each population, giving rise to a sequence of between 12 and 14 populations in the fully observed runs and 12 populations in each partially observed run, the decreasing choice of $\epsilon$ chosen as described in section 3.2. For pMCMC we ran the sampler for the full budget and then thinned the resultant collection of vectors such that the final sample contained 10000. In the case of $\mathcal{D}^{2}$ the average number of particles required was 132 which led to a chain that ran for approximately $7.5\times 10^{5}$ unthinned iterations which was then thinned by a factor of 75 to give the final sample.

From these results it would appear that given the use of the full budget pMCMC provides the better choice. Figure 3(a) shows the posterior learning experienced by the sequential ABC algorithm. The matching box plots for pMCMC is the posterior distribution using the pMCMC algorithm where we use only information gained under the same budget use as with ABC. i.e The budget used to obtain the first population of samples under ABC was $2.1\times 10^{6}$. The corresponding box-plot for population 1 under pMCMC is a snapshot of the chain having used up the same budget of $2.1\times 10^{6}$. It is interesting to note that the shape of the posterior distribution under pMCMC changes little over the sequence of populations leading one to believe that even with little computational expense the posterior distribution inferred by pMCMC is good. However this plot does not tell the whole story. Figure 3(b) shows estimates of the effective sample size of the posterior distributions for pMCMC in this sequence. Effective sample size is small for the populations with the lower computational expense showing that to obtain a good posterior sample from pMCMC, despite the fact that the overall shape of the distribution doesn’t change very much, we must still run for a long time, which is unsurprising. However, if we are only interested in obtaining posterior summaries, it would appear that the budget is less important. Figure 3(c-d) show posterior estimates of the mean and variance factored into the same computational expense groups as with figure 3(a) and figure 3(b). It is clear to see that for pMCMC the estimates are stable and remain reasonably constant. This is in stark contrast to estimates using the ABC approximations which appear less certain before tending toward those estimates gained using pMCMC. The overall trend here then seems to suggest that in this instance pMCMC is the favorable choice. The posterior estimates of mean and variance are stable even given a relatively short run time, and the shape of the distribution is also maintained. To obtain a large uncorrelated sample the chain must be run for a long time, although it is noted that running the ABC sampler for the same length of time does not yield better results. Posterior summaries appear consistent irrespective of the additional runtime. ABC here however is not so good. It could be argued that the posterior distributions are close given the full budget in some circumstances, notably $\mathcal{D}^{2}_{p}$, however with shorter runtime the approximation is much greater and hence inference is poorer as a result.

The results presented here have ignored the issue of tuning the pMCMC algorithm. The posterior inference suggest that pMCMC is the better choice for learning about model parameters but we have tuned the pMCMC at the start, relying on knowledge of the posterior. This knowledge comes with its own expense and is in some sense self-defeating, something that is not so much of an issue for the ABC SMC. The true cost of the pMCMC inference then is somewhat higher than shown here. The ABC SMC proposal variance and tolerance sequence are chosen adaptively and so the initialisation and tuning cost is small. We made the choice to tune pMCMC using information from the posterior earlier due to the fact that in our initial experiments for this model we found that in a number of cases the computational budget used in initialising the pMCMC was large, often larger than our allocated budget. Choosing a sensible $\theta_{0}$, $N$ and $\Sigma_{q}$ with little prior knowledge is difficult. Our initial attempts to find a $\theta_{0}$ involved sampling from the prior distribution, estimating the likelihoods and choosing the parameter vector which maximises the likelihood estimates. Due to the variability in the particle filter estimates away from the true parameter values this step typically involves a large number of particles and hence large expense. Conditional on this hopefully informative choice of $\theta_{0}$ we can then tune $N$, the number of particles in the bootstrap filter, such that each estimate has less expense, by running a number of filters starting with a small number of particles and steadily increasing $N$ until the $\operatorname{Var}(\hat{l})$ is suitably small. Again this has non-negligible expense. We want to try to find as small a value of $N$ as we can get away with for the main inference run and so this iterative procedure can be time consuming. On top of these two tuning steps we typically want to ensure that we choose a good $\Sigma_{q}$. This often involves a pilot MCMC run using a very small proposal variance and then using the variance of the resulting distribution to inform the choice of $\Sigma_{q}$. We found that, in practice, employing these steps to tune pMCMC was itself very expensive. Using 1000 particles in a particle filter to estimate likelihoods for 2000 parameter vectors drawn from the prior for $\mathcal{D}^{2}$ and repeating each 10 times and maximising over the average to choose $\theta_{0}$ was not enough to guarantee that the resulting draw had posterior support. This alone uses 20% of the allocated budget without then tuning the number of particles. The number of particles needed to satisfy the log–likelihood estimate variance criteria is often much larger in the tails of the posterior than at the true values. Add to this estimating $\Sigma_{q}$ and then the appropriate burn–in period and it is easy to see that this operation becomes very expensive. Owen et al. (2014) showed that initialising pMCMC with a random draw from an uninformative prior does not guarantee convergence to the stationary distribution. This was due to the large variability in log-likelihood estimates given by the particle filter in regions of negligible posterior support, rather than any flaw in the theory. It is suggested that a good thing to do in situations in which prior knowledge is poor is to run an ABC SMC algorithm as an aid to tuning and initialising a pMCMC algorithm. This approach is amenable to parallelisation and exploits relative strengths of the two approaches.

The Schlögl model is a well known test model which exhibits bimodal stability at certain parameter values. The system is characterised by a set of 4 reactions involving 3 species:

We can summarise this reaction network via its stoichiometry matrix $S$ and hazard function $h(\mathbf{X},\theta)$:

This system provides an interesting case for which to investigate how influential the measurement error is on the efficiency of posterior sampling for each algorithm.

Figure 4(a) shows 100 realisations from the model given $\theta=(3\times 10^{-7},1\times 10^{-4},0.000773,3.276)$ and $X_{0}=(250,10^{5},2\times 10^{5})$ highlighting the bimodal characteristics given this set of parameter values. This poses an interesting challenge when it comes to parameter inference since it is not necessarily the case that all parameter vectors in a region of space closely surrounding this give the same behaviour. For our investigation of the inference methods being discussed we choose one of these data traces at random and corrupt with Gaussian error, $d_{t}\sim\mathcal{N}(X_{t},\sigma^{2})$, $\sigma^{2}=1$. A second copy of the same underlying trace is then corrupted with the same error distribution but with $\sigma^{2}=10$. We denote these data set $\mathcal{D}^{S}_{1}$ and $\mathcal{D}^{S}_{10}$ respectively. The chosen observed data set is shown in figure 4(b).

In contrast to the simulation study for the Lotka–Volterra model we make use of a set of more informative prior distributions, that is a Gaussian distribution on the log scale, centered at the true values with relatively small standard deviation, $0.5$. In addition we assume knowledge of the initial count $\mathbf{X}_{0}$. For inference in this case we now restrict our focus to assuming known measurement error, and availability of observations of all 3 species in the model, these factors having been explored well in the previous example. The focus of this example is to determine whether the size of the measurement error informs our choice as to which algorithm is better.

Consistent with the results in section 5.1.4 ABC SMC yields posterior distribution with a larger mass in the tails. The plots in figure 5 show that under the larger measurement error, data set $\mathcal{D}^{S}_{10}$, comparative performance is similar to that seen in the previous example. However under small measurement error pMCMC struggles with this computational budget. The chain spends a large amount of time stuck at given parameter values. It is known that for likelihood free pMCMC to be efficient in this context the measurement error must be substantial. Development of pMCMC algorithms for informative observations are subject to ongoing research and typically involve bridging of the latent state conditional on the endpoints. Golightly & Wilkinson (2014) propose an approach to the problem of inference for a Markov jump process with informative observations on a discrete state space by conditioning the hazard function on the end points of the data observations.