Ensemble MCMC: Accelerating Pseudo-Marginal MCMC for State Space Models using the Ensemble Kalman Filter

TL;DR

This paper introduces ensemble MCMC, a method that uses the ensemble Kalman filter to accelerate Bayesian inference in state space models, reducing computational costs while maintaining accuracy.

Contribution

It replaces the particle filter with the ensemble Kalman filter within MCMC, providing a faster alternative for likelihood estimation in complex models.

Findings

eMCMC significantly reduces computational time

Maintains reasonable accuracy across various models

Extensions further improve efficiency

Abstract

Particle Markov chain Monte Carlo (pMCMC) is now a popular method for performing Bayesian statistical inference on challenging state space models (SSMs) with unknown static parameters. It uses a particle filter (PF) at each iteration of an MCMC algorithm to unbiasedly estimate the likelihood for a given static parameter value. However, pMCMC can be computationally intensive when a large number of particles in the PF is required, such as when the data is highly informative, the model is misspecified and/or the time series is long. In this paper we exploit the ensemble Kalman filter (EnKF) developed in the data assimilation literature to speed up pMCMC. We replace the unbiased PF likelihood with the biased EnKF likelihood estimate within MCMC to sample over the space of the static parameter. On a wide class of different non-linear SSM models, we demonstrate that our new ensemble MCMC (eMCMC) method can significantly reduce the computational cost whilst maintaining reasonable accuracy. We also propose several extensions of the vanilla eMCMC algorithm to further improve computational efficiency. Computer code to implement our methods on all the examples can be downloaded from https://github.com/cdrovandi/Ensemble-MCMC.

Tables

Table 1. Table 1: Efficiency comparisons for the four non-linear population ecology models.

Model	filter	$N$	ESS	Time (h)	ESS/Time
Ricker	BPF	50000	920	36.8	25
Ricker	EnKF	250	2400	0.14	17000
Ricker	EnKF + correlation	25	2100	0.07	30000
theta-logistic	BPF	50000	410	40.8	10
theta-logistic	EnKF	250	500	0.48	1040
theta-logistic	EnKF + correlation	25	540	0.06	9000
mate-limited	BPF	50000	770	37.6	20
mated-limited	EnKF	200	1460	0.35	4200
mated-limited	EnKF + correlation	25	1800	0.07	25700
flexible-allee	BPF	50000	750	37.0	20
flexible-allee	EnKF	250	1750	0.26	6700
flexible-allee	EnKF + correlation	25	1600	0.08	20000

Table 2. Table 2: Summaries for the Lotka Volterra application: number of particles N 𝑁 N , standard deviation of the noise in the log-posterior ( τ 𝜏 \tau ) at the posterior median, acceptance rate, multivariate effective sample size (ESS), wall clock time in seconds and ESS per second.

Filter	$N$	$\tau$	Acc. rate	ESS	Time (s)	ESS$/$Time
	${𝒟}_1$ (51 obs. every 1 time unit)
APF	55	1.4	0.11	1117	26298	0.042
EnKF	150	1.4	0.11	1762	12299	0.143
	${𝒟}_2$ (26 obs. every 2 time units)
APF	350	1.4	0.08	1156	165015	0.007
EnKF	65	1.4	0.11	2054	5282	0.389

Table 3. Table 3: Summaries for the autoregulatory example: number of particles N 𝑁 N , standard deviation of the noise in the log-posterior ( τ 𝜏 \tau ) at the posterior median, acceptance rate, multivariate effective sample size (ESS), wall clock time in seconds and ESS per second.

Filter	$N$	$\tau$	Acc. rate	ESS	Time (s)	ESS$/$Time
	${𝒟}_1$ (${\sigma }_1={\sigma }_2=1$)
APF	400	1.3	0.15	3348	72081	0.046
EnKF	200	1.5	0.13	2972	15862	0.187
	${𝒟}_2$ (${\sigma }_1={\sigma }_2=0$)
APF	2000	1.4	0.11	3314	403456	0.0082
EnKF	370	1.4	0.11	3176	34112	0.0931

Table 4. Table 4: Parameters of the Liley model.

Parameter	Prior	Value for simulation
${\Gamma }_{ee}$	$𝒰(0.1,2)$	0.10631
${\Gamma }_{ei}$	$𝒰(0.1,2)$	0.64105
${\Gamma }_{ie}$	$𝒰(0.1,2)$	0.46477
${\Gamma }_{ii}$	$𝒰(0.1,2)$	0.28663
${\gamma }_{ee}$	$𝒰(100,1000)$	291.5
${\gamma }_{ei}$	$𝒰(100,1000)$	697.76
${\gamma }_{ie}$	$𝒰(10,500)$	458.67
${\gamma }_{ii}$	$𝒰(10,500)$	82.33
${\overline{p}}_{ee}$	$𝒰(0,10000)$	6603.4
${\overline{p}}_{ei}$	$𝒰(0,10000)$	2625.7
$\sigma$	$𝒰(0,10)$	0.01

Equations

$$y_t|x_t\sim \mathcal{N}(P{x}_t,S)$$

$$p(x,y|\theta)=p(x_{0})\prod_{t=1}^{T}\big{[}p(x_{t}|x_{t-1},\theta)p(y_{t}|x_{t},\theta)\big{]}.$$

$$L(\theta)=\int p(x_{0})\prod_{t=1}^{T}\big{[}p(x_{t}|x_{t-1},\theta)p(y_{t}|x_{t},\theta)\big{]}dx.$$

$$r=\frac{\hat{L}({\theta }^{\ast })\pi ({\theta }^{\ast })q({\theta }_{i-1}|{\theta }^{\ast })}{\hat{L}({\theta }_{i-1})\pi ({\theta }_{i-1})q({\theta }^{\ast }|{\theta }_{i-1})}.$$

$$p_{\text{enkf}}(x_t|y_{1:t-1})=\mathcal{N}(x_t;{\hat{\mu }}_{t|t-1},{\hat{\Sigma }}_{t|t-1})$$

$$\left(\begin{array}{c}{X}_t\\ Y_t\end{array}\right)\sim N\left\{\left(\begin{array}{c}{\hat{\mu }}_{t|t-1}\\ P{\hat{\mu }}_{t|t-1}\end{array}\right),\left(\begin{array}{cc}{\hat{\Sigma }}_{t|t-1}& {\hat{\Sigma }}_{t|t-1}{P}^{\prime }\\ P{\hat{\Sigma }}_{t|t-1}& P{\hat{\Sigma }}_{t|t-1}{P}^{\prime }+S\end{array}\right)\right\}.$$

$$p_{\text{enkf}}(x_t|y_{1:t})=\mathcal{N}(x_t;{\hat{\mu }}_{t|t},{\hat{\Sigma }}_{t|t})$$

$${\hat{K}}_t={\hat{\Sigma }}_{t|t-1}{P}^{\prime }(P{\hat{\Sigma }}_{t|t-1}{P}^{\prime }+S{)}^{-1}.$$

$$L(\theta )=p(y_1|\theta )\prod _{t=2}^{T}p(y_t|y_{1:t-1},\theta ).$$

$$p_{\text{enkf}}^N(y_t|y_{1:t-1},\theta )=\mathcal{N}(y_t;P{\hat{\mu }}_{t|t-1},P{\hat{\Sigma }}_{t|t-1}{P}^{\prime }+S)$$

$$L_{\text{enkf}}^N(\theta )=\prod _{t=1}^T{p}_{\text{enkf}}^N(y_t|y_{1:t-1},\theta ).$$

$$p_{\text{enkf}}^N(\theta |y)\propto L_{\text{enkf}}^N(\theta )p(\theta ),$$

$$c(k,v)$$

$$\hat{\mathcal{N}}(y;\mu ,\Sigma )$$

$$\psi {\left(M_N-(y-{\mu }_N)(y-{\mu }_N{)}^{\top }/(1-1/N)\right)}^{(N-d-3)/2},$$

$$d{X}_t$$

$$\alpha (X_t,\theta )$$

$$\Sigma$$

$$x_{i+1}=x_i+\alpha (x_i,\theta )\Delta t+{\Sigma }^{1/2}\sqrt{\Delta t}{z}_{i+1},$$

$$\text{Pr}\left\{X_{1,t+dt}=x_{1,t}+1,X_{2,t+dt}=x_{2,t}|x_t\right\}$$

$$\text{Pr}\left\{X_{1,t+dt}=x_{1,t}-1,X_{2,t+dt}=x_{2,t}+1|x_t\right\}$$

$$\text{Pr}\left\{X_{1,t+dt}=x_{1,t},X_{2,t+dt}=x_{2,t}-1|x_t\right\}$$

$$h(X_t,c)=(c_1{x}_{1,t},c_2{x}_{1,t}{x}_{2,t},c_3{x}_{2,t}{)}^{\prime }.$$

$$Y_{t}|X_{t}\sim\mathcal{N}\left\{\left(\begin{array}[]{c}x_{1,t}\\ x_{2,t}\end{array}\right)\,,\,\left(\begin{array}[]{cc}\sigma_{1}^{2}&0\\ 0&\sigma_{2}^{2}\end{array}\right)\right\}.$$

$$\theta ={\left(\log c_1,\log c_2,\log c_3,\log {\sigma }_1,\log {\sigma }_2\right)}^{\prime }.$$

$${\mathcal{R}}_1:$$

$${\mathcal{R}}_3:$$

$${\mathcal{R}}_5:$$

$${\mathcal{R}}_7:$$

$$h(x_t,c)=(c_1{x}_{1,t}{x}_{5,t},c_2{x}_{2,t},c_3{x}_{1,t},c_4{x}_{3,t},c_5{x}_{4,t}(x_{4,t}-1)/2,c_6{x}_{5,t},c_7{x}_{3,t},c_8{x}_{4,t}{)}^{\prime }.$$

$$c=(0.1,0.7,0.35,0.2,0.1,0.9,0.3,0.1{)}^{\prime }.$$

$$\left(\frac{d}{dt}+{\gamma }_{ek}\right)\left(\frac{d}{dt}+{\bar{\gamma }}_{ek}\right)I_{ek}(t)$$

$$\left(\frac{d}{dt}+{\gamma }_{ik}\right)\left(\frac{d}{dt}+{\bar{\gamma }}_{ik}\right)I_{ik}(t)$$

$${\left(\frac{d}{dt}+v\Lambda \right)}^2{\Phi }_{ek}(t)$$

$${\Phi }_{ek}(t)$$

$$y_i=h_e(i\cdot \Delta t)+z_i$$

Full text

Keywords: data assimilation, ensemble Kalman filter, particle filter, particle MCMC, pseudo-marginal MCMC, state space models

1 Introduction

Particle Markov chain Monte Carlo (pMCMC, Andrieu et al.,, 2010) is now a popular method for performing Bayesian statistical inference on challenging state space models (SSMs) with unknown static parameters. The appeal of particle MCMC is that it is a pseudo-marginal method (Andrieu and Roberts,, 2009), which attempts to mimic the ideal sampler that proposes directly over the space of the static parameters and integrates out the hidden states. Furthermore it is an exact approximation, exactly targeting the true posterior distribution.

Each static parameter proposal in pMCMC is evaluated using a particle filter (Gordon et al.,, 1993). Particle filters were originally proposed to solve the state space filtering problem: inferring the state parameters at a given time under known static parameters. To do so they propagate a set of particles through the state space model, and use a weighting and resampling process to concentrate on the particles with significant posterior weights. Using particle filters in pMCMC is costly. Firstly, each particle filter involves processing the entire dataset. Secondly, a particle filter can require a large number of particles, especially when the data is highly informative and/or the model is misspecified. This is because there must be enough particles to randomly propagate forwards to produce good matches to unlikely data. Thus, despite the popularity of particle MCMC, it is generally a highly computationally intensive method.

Data assimilation (DA) is a field of research originating in the geosciences, initially based on the problem of numerical weather prediction. The task most commonly addressed in this field is the estimation of the state of a dynamical system, based on a dynamic model (usually a system of partial differential equations) and noisy and/or indirect measurements of this state. In this paper we take inspiration from the DA literature to propose a new approach to estimating the posterior distribution of static parameters in SSMs.

The field of DA has evolved in parallel to other fields in which SSMs play an important role, such as target tracking, economics and statistical ecology. The distinguishing feature of problems in DA is the large dimension of the state space. For example, in numerical weather prediction the state space consists of a representation of the state of the atmosphere across the globe, which for modern applications can have dimension $d_{x}$ of the order of $10^{9}$ (van Leeuwen,, 2015). The traditional approach to estimating the dynamic state in DA is to use approaches that solely estimate the mode of the state posterior (e.g. 4DVar) or Kalman filters that make use of approximations so as to avoid storing the full state covariance, whose size scales quadratically in the state dimension. Such methods have huge practical importance and are still deployed in DA applications, but more recent research has focussed on methods that improve the accuracy of state estimation when using nonlinear dynamics. As in other fields where this the case, particle filters are an important methodology.

Particle filters are not the usual method of choice in DA. The reason is their degeneracy when used on states of high dimension (Snyder et al.,, 2008). This degeneracy arises due to the limitations of importance sampling in high dimensions: the variance of importance sampling estimators depends on the distance between the target and proposal distributions, and this distance grows with dimension such that the variance is only controlled by using a number of importance points that is exponential in the dimension (see Agapiou et al.,, 2017 for a review). To combat this degeneracy, the approach usually taken in the particle filtering literature is to introduce diversity into the sample though using MCMC updates of the state (Beskos et al.,, 2014). However, in many problems in DA, MCMC updates are often not available due to the use of an intractable dynamic model, and where available may have a low acceptance rate. An alternative means of maintaining diversity is given by the ensemble Kalman filter (EnKF) (Evensen,, 1994; see Katzfuss et al.,, 2016 for a tutorial). This approach propagates a set of particles (often refered to as “ensemble members”) through the dynamic model in the same way as the bootstrap particle filter, but instead uses these particles to approximate a Gaussian representation of the state distribution. This method approximates the Kalman filter and provides a means to avoid storing and manipulating the state covariance matrix. Further, it has also been shown to perform well when applied to non-linear dynamic models. Its performance is often superior to the particle filter in cases where particle filter suffers from degeneracy, including but not limited to the case of a state space of high dimensions.

Methods from DA have also been applied to the situation of inferring an unknown static parameter simultaneously with the state. The standard approach is to augment the state vector with the static parameter, then to apply one of the previously mentioned filters to this augmented state (see, for example, Evensen, (2007). In this case, the EnKF assumes that both parameters and states follow a linear Gaussian state space model. When this assumption is unreasonable, another approach is to combine the EnKF likelihood with a particle representation of the static parameter (Stroud et al.,, 2018; Katzfuss et al.,, 2019). To mitigate degeneracy, the static parameter is allowed to dynamically vary, by adding Gaussian noise to each parameter particle. This step can be further refined via the kernel resampling strategy of Liu and West, (2001). However, additional tuning parameters must be specified (e.g. to control the smoothness of the kernel) and the particle approximation can be sensitive to these choices, and in in particular, the number of particles used (Vieira and Wilkinson,, 2016).

In this paper we propose a new method, called ensemble MCMC (eMCMC), which is a compromise between pMCMC and EnKF. It can be viewed as a pMCMC algorithm, where each use of a particle filter is replaced with using the EnKF (with fixed static parameters). The motiviation is to reduce the number of particles/members required and to cope better with informative or surprising data. Thus we reduce the computational cost relative to pMCMC while improving the accuracy relative to EnKF with time-varying static parameters.

Moreover, we propose several extensions of eMCMC to further improve computational efficiency and reduce the bias in the EnKF estimate of the likelihood. Some of the extensions may be of interest to the data assimilation community more generally. The basic idea underlying eMCMC is also suggested in Katzfuss et al., (2019). However, Katzfuss et al., (2019) consider only a simple example. Here we demonstrate that the method can be successful on a wider variety of more challenging applications and also develop several extensions just discussed.

The rest of this article is structured as follows. In Section 2 we provide the necessary background on state space models, pMCMC and EnKF to understand our method. Our ensemble MCMC approach together with extensions is described in Section 3. Section 4 shows the results of our approach on a wide class of different models. We discuss limitations, further extensions and possible future research in Section 5.

2 Background

This section describes relevant existing work. Section 2.1 defines state space models. Section 2.2 describes pseudo-marginal Metropolis Hastings and the bootstrap particle filter, which can be used to perform inference for these models. Section 2.3 introduces the EnKF.

2.1 State Space Models

A state space model is a model for sequential data. It introduces a Markov chain of latent states $x_{1},\ldots,x_{T}$. Independent noisy observations $y_{t}$ are available that depend on the state $x_{t}$. Let $x$ denote the collection of all latent states and $y$ the collection of all observations. The model can be defined using an evolution distribution for $x_{t+1}|x_{t},\theta$ and an observation distribution $y_{t}|x_{t},\theta$. Here $\theta$ is a vector of parameters controlling the model’s behaviour. We also specify a distribution for an initial state $x_{0}$. For more background on state space models see for example Särkkä, (2013).

Throughout the paper we will make some standard assumptions about state space models. We will assume that each $x_{t}$ and $y_{t}$ are random vectors with support $\mathbb{R}^{d_{x}}$ and $\mathbb{R}^{d_{y}}$ respectively. In this section we assume the distributions above – evolution, observation and initial state – have densities $p(x_{t+1}|x_{t},\theta),p(y_{t}|x_{t},\theta)$ and $p(x_{0})$. The material immediately generalises to the case where some or all of these distributions have probability mass functions instead (by interpreting these as densities with respect to counting measure). This is required in several of our examples. A point mass can be used for the initial state distribution if the initial state is known.

As we shall see, the EnKF is restricted to certain observation models. Hence in this paper we focus on one particular case,

[TABLE]

where $P$ is a $d_{y}\times d_{x}$ matrix and $S$ is a variance matrix, possibly a function of $\theta$. (We assume conditional independence of the $y_{t}$’s given $x$ and $\theta$.) The EnKF can also be used where $P$ is replaced by a time dependent matrix $P_{t}$. The case where $P=I$ gives a complete observation regime, in the sense that all components of $x_{t}$ have a corresponding noisy observation. In contrast, a partial observation regime only allows observation of a subset of the components e.g. by taking $P$ to be a projection matrix.

In practice we may wish to model states at a finer time discretisation than that at which the observation data is available. For example, consider the case where we only have observations $y_{t}$ at $t=k,2k,\ldots,kL$. This can easily be converted into the framework described above, by defining a state space model with $x^{*}_{\tau}=x_{\tau k}$ and $y^{*}_{\tau}=y_{\tau k}$ for $\tau=1,2,\ldots,L$.

The joint density of the latent states and observations in a state space model is:

[TABLE]

The likelihood can be found by marginalisation i.e. integrating out the latent states $x$,

[TABLE]

(If there is an observation $y_{0}$, a factor $p(y_{0}|x_{0},\theta)$ can easily be included in (2) and (3).)

Bayesian inference assigns a prior $p(\theta)$ to the parameters and targets the posterior $p(\theta|y)\propto p(\theta)L(\theta)$. The likelihood $L(\theta)$ typically cannot be evaluated as it is a high dimensional integral. One strategy to perform inference is to instead consider an augmented target density (often of interest in its own right), the joint posterior $p(\theta,x|y)\propto p(\theta)p(x,y|\theta)$. The posterior for $\theta$ can then be obtained by marginalisation.

2.2 Pseudo-marginal MCMC, particle filters, and particle MCMC

Monte Carlo algorithms are designed to sample from a target distribution, often a Bayesian posterior distribution. Markov chain Monte Carlo (MCMC) does so using a Markov chain which converges to the target distribution in the long run. Performing each update in MCMC typically requires likelihood calculations, which are not possible for models with intractable likelihoods. However it is often possible to produce unbiased likelihood estimates. Algorithm 1, pseudo-marginal Metropolis Hastings (PMMH) (Andrieu and Roberts,, 2009), makes use of these to perform parameter inference. Unbiased likelihood estimates for state space models can be produced by particle filter algorithms. Algorithm 2 presents the basic bootstrap particle filter (BPF) used in this paper, but there are many variations. For more details see for example Doucet and Johansen, (2011), Särkkä, (2013) and Fearnhead and Künsch, (2018). For proof that the particle filter likelihood estimate is indeed unbiased see Del Moral, (2004) and Pitt et al., (2010).

Combining the PMMH algorithm with a particle filter can target $p(\theta|y)$ for state space models. Andrieu et al., (2010) extend this approach to give particle MCMC (pMCMC), which targets the joint posterior $p(\theta,x|y)$; we refer the reader to this paper for a full description of pMCMC.

2.2.1 PMMH tuning

PMMH using BPF likelihood estimates has several tuning choices. This section sets out the approach we use to make these choices in this paper. Our choices are consistent with the theoretical analyses of Sherlock et al., (2015) and Doucet et al., (2015), who derive tuning recommendations under two different sets of simplifying assumptions.

We select the number of particles $N$ for the PF prior to running Algorithm 1 so that the estimated log-likelihood at a representative parameter value has a standard deviation of roughly 1.5. The parameter value used should have good support under the posterior; we typically use marginal posterior medians from exploratory analyses.

We use a normal random walk proposal distribution: $\theta^{*}\sim\mathcal{N}(\theta_{i-1},\Sigma)$. We take $\Sigma$ to be an estimate of the posterior variance, again taken from exploratory analyses. Sherlock et al., (2015) and Doucet et al., (2015) provide guidance for scaling the variance matrix by a scalar to improve performance, finding that this was helpful for high dimensional target distributions for instance. We did not find this necessary for our analyses of low dimensional targets, but make use of this approach in Section 4.5 where an 11-dimensional target is considered.

2.3 Ensemble Kalman Filter

Here, we give a brief overview of the ensemble Kalman filter (Evensen,, 1994) and refer the reader to Katzfuss et al., (2016) and the references therein for further details.

Consider the task of sampling the filtering density $p(x_{t}|y_{1:t})$ where $y_{1:t}=(y_{1},\ldots,y_{t})$. (We omit explicit conditioning on the parameter vector $\theta$ throughout this section.) The EnKF generates approximate draws from $p(x_{t}|y_{1:t})$ via a sequence of forecasting and updating steps. Suppose that a sample $\{x_{t-1}^{(1)},\ldots,x_{t-1}^{(N)}\}$ (known as the filtering ensemble) is available at time $t-1$ from $p(x_{t-1}|y_{1:t-1})$. The forecast ensemble $\{\tilde{x}_{t}^{(1)},\ldots,\tilde{x}_{t}^{(N)}\}$ is obtained by drawing $\tilde{x}_{t}^{(i)}\sim p(\cdot|x_{t-1}^{(i)}),i=1,\ldots,N$. The forecast density $p(x_{t}|y_{1:t-1})$ is then approximated by

[TABLE]

where $\mathcal{N}(\cdot;\mu,\Sigma)$ denotes the multivariate Gaussian density with mean $\mu$ and variance matrix $\Sigma$. The quantities $\hat{\mu}_{t|t-1}$ and $\hat{\Sigma}_{t|t-1}$ are typically taken to be the sample mean and variance computed from the forecast ensemble (some extensions of the EnKF use alternative estimates; see Katzfuss et al., (2016) for some common approaches). Now, given the linear Gaussian form of (1), the joint distribution of $X_{t}$ and $Y_{t}$ (given $y_{1:t-1}$) can be obtained approximately as

[TABLE]

Hence, conditioning on $Y_{t}=y_{t}$ gives

[TABLE]

where $\hat{\mu}_{t|t}=\hat{\mu}_{t|t-1}+\hat{K}_{t}(y_{t}-P\hat{\mu}_{t|t-1})$, $\hat{\Sigma}_{t|t}=(I_{d_{x}}-\hat{K}_{t}P)\hat{\Sigma}_{t|t-1}$ and $\hat{K}_{t}$ is an estimate of the Kalman gain, that is

[TABLE]

It is then straightforward to generate samples from (5) to be used as the filtering ensemble at the next time point. However, rather than explicitly calculate the filtering density in (5), the standard implementation of the EnKF (see e.g. Katzfuss et al.,, 2016) performs a shifting step, which is equivalent under the Gaussianity assumption (4) (and a Gaussian prior for $x_{0}$). For each particle (known in this context as an ensemble member), we compute $x_{t}^{(i)}=\tilde{x}_{t}^{(i)}+\hat{K}_{t}(y_{t}-\tilde{y}_{t}^{(i)})$, where $\tilde{y}_{t}^{(i)}\sim\mathcal{N}(P\tilde{x}_{t}^{(i)},S)$ is a pseudo-observation. Note that the shifting step only requires a draw from a $d_{y}$ variate Gaussian distribution per particle, rather than a draw from a $d_{x}$ variate Gaussian if (5) is sampled directly. Moreover, the shifting approach does not make the strong assumption that the forecast ensemble is Gaussian distributed. There are also other schemes for performing the shifting that we do not consider here. The randomness in the shifting step leads to the variant of the EnKF described being referred to as the “stochastic” EnKF; a commonly used alternative is to use a deterministic shift of the ensemble members (Tippett et al.,, 2003). In what follows, we typically use the stochastic shifting step.

Given a sample $\{x_{0}^{(i)},\ldots,x_{0}^{(N)}\}$ from the state prior, the EnKF recursively alternates between computing the forecast ensemble, and shifting each ensemble member, to give approximate draws from the filtering density $p(x_{t}|y_{1:t})$, $t=1,\ldots,T$. We state a version of the EnKF based on a shifting step below as Algorithm 3.

The EnKF is most easily understood in the context of a linear Gaussian state space model. In this special case, the filtering distribution $p_{\textrm{enkf}}(x_{t}|y_{1:t-1})$ converges to the true filtering distribution as the number of ensemble members $N\to\infty$. Essentially, the EnKF converges to the Kalman filter. For finite $N$ and a linear state space model, the EnKF approximates the Kalman filter by replacing the mean and variance of the forecast distribution with their sample equivalents. The resulting dimension reduction (that only requires storing and manipulating $d_{x}$-vectors) avoids the potentially expensive calculation and storage of the forecast variance matrix. Moreover, several studies (e.g. Lei et al.,, 2010; Houtekamer et al.,, 2014; Katzfuss et al.,, 2019) have found that the EnKF shifting step works well for non-Gaussian evolution densities. We therefore consider the use of the EnKF likelihood inside a Metropolis-Hastings scheme. We provide a motivation and give details of our proposed approach in the next section.

3 Ensemble MCMC

It is well known that as the variance of the likelihood estimator increases, the acceptance probability of the pseudo-marginal MH scheme rapidly decreases to 0 (Pitt et al.,, 2012), resulting in slow mixing behaviour of the parameter chains. As discussed in Section 2.2.1, a value of $N$ (the number of particles) can be chosen to balance mixing performance and computational cost. Nevertheless, in scenarios where the stochasticity inherent in the state process dominates the observation variance, the number of particles required to maintain a reasonable likelihood variance is likely to render BPF-driven PMMH computationally infeasible. Methods that aim to alleviate this problem include the use of an auxiliary particle filter (see e.g. Golightly and Wilkinson,, 2015), which requires careful exploitation of the model structure in order to propagate particles conditional on the observations. Our proposed approach is simple to implement and, for the simplest implementation, does not require the specification of any additional tuning parameters.

Here we outline our proposed ensemble MCMC (eMCMC) algorithm. In essence, this is PMMH using the EnKF as fast replacement for the BPF to estimate the likelihood $L(\theta)$. We describe the eMCMC algorithm in more detail now, and in the following subsections we discuss some extensions to improve its efficiency.

First we derive a likelihood estimate based on EnKF calculations. Recall that the (marginal) likelihood can be factorised as

[TABLE]

From (4) it follows that an EnKF approximation of $p(y_{t}|y_{1:t-1},\theta)$ is

[TABLE]

which can easily be computed for each $t=1,\ldots,T$, with the notational convention that $p(y_{1}|\theta)=p(y_{t}|y_{1:0},\theta)$. The need for the explicit dependence on $N$ for this likelihood estimate will become clearer later in this section. The overall approximation to the likelihood is given by

[TABLE]

The EnKF including likelihood estimation is given by Algorithm 3. The ensemble MCMC scheme is then implemented by running Algorithm 1 with $\hat{L}$ replaced by $\hat{L}_{\textrm{enkf}}^{N}$. One issue in implementing eMCMC is how to perform tuning. Due to the absence of specialised theory, we use the same tuning guidance as for PMMH with BPF likelihood estimates, described above in Section 2.2.1.

It is worth emphasising that, unlike pMCMC, the eMCMC posterior

[TABLE]

does not in general equal the posterior $\pi(\theta|y)$ exactly. The reason is that, unlike the BPF, the EnKF gives a biased estimator of $L(\theta)$, precluding its use for exact approximate inference. Nevertheless, as noted by Stroud et al., (2010), Stroud et al., (2018) and Katzfuss et al., (2019) among others, the variance of the likelihood estimator under the EnKF can be relatively small, suggesting that use of EnKF inside a Metropolis-Hastings scheme is likely to be of practical use, particularly in scenarios when the BPF is computationally prohibitive.

In fact, even when the forecast ensemble is exactly Gaussian distributed for all $t$, the EnKF posterior still does not target the exact posterior, since $\mathcal{N}(y_{t}\,;\,P\hat{\mu}_{t|t-1}\,,\,P\hat{\Sigma}_{t|t-1}P^{\prime}+S)$ is a biased estimate of the idealised normal density should we be able to take $N\rightarrow\infty$. Thus, for finite $N$, the eMCMC target is not the idealised eMCMC target, $p_{\textrm{enkf}}^{\infty}(\theta|y)$. However, we find empirically that our method appears to be weakly dependent on $N$. Given this, we suggest to choose $N$ to maximise the computational efficiency by borrowing similar advice from the pseudo-marginal literature (as described in Section 2.2.1). Interestingly, there is an exactly unbiased estimator of a normal density given a sample from it, and we exploit this in Section 3.3. We discuss the unbiased version and other extensions below.

3.1 Randomised Quasi-Monte Carlo

For this subsection, all quantities are conditioned on $\theta$ so we drop it for notational convenience. At iteration $t$ of the EnKF we are interested in estimating $\mu_{t|t-1}$ and $\Sigma_{t|t-1}$, so that we can approximate the conditional likelihood $p(y_{t}|y_{1:t-1})$ with a Gaussian density. These moments can be estimated via firstly performing the shifting step at $t-1$ and conditional on the result simulating from the forward evolution density. This is effectively an approximate sample from the joint distribution $p(x_{t}|x_{t-1})p(x_{t-1}|y_{1:{t-1}})$. Then $\mu_{t|t-1}$ and $\Sigma_{t|t-1}$ are estimated from the $N$ ensemble members.

Often it is possible to write the simulation from a standard statistical distribution as a function of a uniform random number. For example, to simulate from a $y\sim\mathcal{N}(\mu,\sigma^{2})$ distribution we can compute the following, $y=\mu+\sigma\cdot\Phi^{-1}(u)$ where $u\sim\mathcal{U}(0,1)$ and $\Phi^{-1}(u)$ is the quantile function of the standard normal density. Assume that we can write the evolution density as a function of $m$ uniform random variates. Then, we require $d_{y}+m$ uniform random numbers to approximately simulate from $x_{t},x_{t-1}|y_{1:{t-1}}$ ($d_{y}$ for the shifting step and $m$ for simulating the evolution density). Given $N$ particles, we use $N\times(d_{y}+m)$ uniform random numbers for estimating $\mu_{t|t-1}$ and $\Sigma_{t|t-1}$. The naive approach is to draw these via pseudo-random numbers. However, significant variance reduction could be achieved by simulating from the $(d_{y}+m)$-dimensional object $N$ times using randomised quasi-Monte Carlo (RQMC, Niederreiter,, 1992). QMC is well known to generate a sequence of numbers that have superior space filling properties in the unit hypercube compared to pseudo-random numbers. The randomised component ensures that expectations can be estimated unbiasedly. Random numbers from the joint distribution of interest, $x_{t},x_{t-1}|y_{1:{t-1}}$, can be achieved via transforming the RQMC numbers as recently discussed. We use this approach to bring down the variance of the estimators of $\mu_{t|t-1}$ and $\Sigma_{t|t-1}$, which hopefully reduces the variance of the estimator for $p(y_{t}|y_{1:t-1})$. For generating the RQMC numbers in this paper, we use the scrambled Sobol’s net, i.e. the scrambled $(t,m,s)$-net in base $b=2$.

RQMC has recently received increasing attention in the statistics community. Tran et al., (2017) document a faster convergence in their Variational Bayes updating procedure when the noisy gradient is computed using RQMC, Drovandi and Tran, (2018) use it to reduce the variance of expected utility estimation within Bayesian optimal design and Gerber and Chopin, (2015) show the efficiency of RQMC in particle filtering. However their application to particle filtering requires considerable ingenuity (using a Hilbert curve method to perform resampling). It is interesting to note the ease with which RQMC can be exploited in the EnKF in comparison.

3.2 Correlated eMCMC

As mentioned above, the EnKF requires generating random numbers for the shifting step and simulating the evolution density. The former can be generated by standard normal random variates and the Cholesky factorisation of the covariance matrix. We assume in this section that the evolution density can be simulated either directly or indirectly via a suitable transformation with standard normal random numbers. Denote the collection of these random numbers required in the EnKF as $u$.

Deligiannidis et al., (2018) and Dahlin et al., (2015) develop the correlated pseudo-marginal MCMC method where they consider the joint target density $p(\theta,u|y)$ where $u$ are random numbers required to estimate the likelihood unbiasedly, $p(y|\theta,u)$. It is easy to show that the $\theta$-marginal of the joint distribution is the posterior of interest, $p(\theta|y)$. Assume that $u$ are independent standard normal random variates. The idea of correlated pseudo-marginal is to induce correlation in successive likelihood estimates in MCMC by correlating the $u$ random numbers. This can have the effect of mitigating “sticky” behaviour often seen in pseudo-marginal chains since, in the correlated scheme, if the likelihood is overestimated at the current iteration, it is also likely to be overestimated at the next. The joint proposal distribution of the correlated pseudo-marginal method is given by $q(\theta^{*},u^{*}|\theta,u)=q(\theta^{*}|\theta)\mathcal{N}(u^{*};\sqrt{1-\sigma_{u}^{2}}u,\sigma_{u}^{2}\mathcal{I})$, where $\mathcal{I}$ is the identity matrix. The proposal for $u$ is the Crank-Nicolson proposal and it is invariant with respect to the marginal distribution of $u$. $\sigma_{u}^{2}$ is an additional tuning parameter that is typically set to be small so that $u^{*}$ is highly correlated with $u$.

Here we consider applying this correlated pseudo-marginal approach to our eMCMC method, with the motivation that a smaller ensemble size $N$ can be used, reducing computational cost. Note that BPF driven pMCMC requires additional modification to accommodate this approach, as it did for RQMC. Essentially, the resampling step has the effect of breaking down correlation between successive likelihood estimates. To alleviate this problem, the particles can be sorted before propagation e.g. using a Hilbert sorting procedure (Deligiannidis et al.,, 2018) or simple Euclidean sorting (Choppala et al.,, 2016). The random numbers used in the resampling step itself should also be updated using the Crank-Nicolson proposal. Since the eMCMC scheme does not use resampling, incorporating correlation is straightforward.

3.3 Unbiased Ensemble Kalman Filter Likelihood

As mentioned earlier, even if the sample from the forecast distribution was exactly Gaussian for some $t$, the corresponding EnKF likelihood estimate for $y_{t}$ would not be unbiased. In general, for some data $y$ and a sample of size $N$ from a Gaussian distribution, $x=x_{1},\ldots,x_{N}\sim\mathcal{N}(\mu,\sigma)$, the density estimator $\mathcal{N}(y;\mu_{N},\Sigma_{N})$ is not an unbiased estimator of $\mathcal{N}(y;\mu,\Sigma)$ where $\mu_{N}$ and $\Sigma_{N}$ are the sample mean covariance computed from the sample $x$. Given the bias present in the EnKF likelihood estimate, even when the Gaussian assumption is correct, the EnKF posterior, unlike standard pseudo-marginal, theoretically depends on $N$.

Even though we demonstrate empirically in Section 4 that the eMCMC posterior seems to be only weakly dependent on $N$, we present a new approach now that will likely be less sensitive to $N$. Interestingly, there does exist an unbiased estimator of a Gaussian density given only an iid sample from the same Gaussian density. Using the notation of Ghurye and Olkin, (1969), let

[TABLE]

and for a square matrix $A$ write $\psi(A)=|A|$ if $A>0$ and $\psi(A)=0$ otherwise, where $|A|$ is the determinant of $A$ and $A>0$ means that $A$ is positive definite. The result of Ghurye and Olkin, (1969) shows that an exactly unbiased estimator of ${\cal N}(y;\mu,\Sigma)$ is (in the case where $y$ is Gaussian and $N>d+3$ where $d$ is the dimension of $y$)

[TABLE]

where $M_{N}=(N-1)\Sigma_{N}$. We propose to replace the standard Gaussian density estimator in the EnKF likelihood estimator with this alternative estimator. Note that this estimator has also been used in Price et al., (2018) for approximating intractable likelihoods in simulation-based likelihood-free estimation problems.

We refer to the method when we use the unbiased Gaussian density estimator in the EnKF likelihood estimator as ueMCMC. We stress that this approach still does not target the true posterior, but at least will not depend on the number of ensemble members $N$ when the Gaussian assumption is correct, i.e. the target is exactly the idealised approximation, $p_{\textrm{enkf}}^{\infty}(\theta|y)$. Even though the forecast density is unlikely to be exactly Gaussian in practice, we do expect ueMCMC to be less sensitive to $N$ compared with eMCMC. We note that this method might be particularly useful when combined with the correlated approach in Section 3.2, since it might be sufficient to use a very small $N$ to achieve reasonable computational efficiency, but the small $N$ may produce bias in the eMCMC posterior compared to the idealised eMCMC posterior.

3.4 Early rejection

Prangle et al., (2018) apply pMCMC in the setting of approximate Bayesian computation (ABC). They outline a method for rejecting proposed values of $\theta$ that have a small estimated likelihood without running the whole particle filter. In Everitt and Sibly, (2019) it is shown that this approach can be extended to pMCMC when using the BPF. A similar approach may be used in eMCMC. Suppose that the likelihood estimate from the EnKF is implemented sequentially, as the EnKF is running. Recall from (8) that the EnKF likelihood estimate is a product, $L_{\textrm{enkf}}^{N}(\theta)=\prod_{t=1}^{T}\alpha_{t}$. Here $\alpha_{t}=\mathcal{N}\left(y_{t};P\hat{\mu}_{t\mid t-1},P\Sigma_{t\mid t-1}P^{\prime}+S\right)$ is calculated in iteration $t$ of the EnKF. We are guaranteed that an upper bound on $\alpha_{t}$ is given by $B(\theta^{*}):=\mathcal{N}\left(\mathbf{0};\mathbf{0},S(\theta^{*})\right)$. Thus we have that $\alpha_{t}/B\leq 1$. This fact ensures that $\hat{r}_{\text{enkf}}^{(\tau)}:=\prod_{t=1}^{\tau}\alpha_{t}/B$ is an upper bound on $L_{\textrm{enkf}}^{N}(\theta)/B^{T}$ which can be calculated at iteration $\tau$ of the EnKF.

We use this property to propose an “early rejection” algorithm. The idea is that during an EnKF run, as soon as $\hat{r}_{\text{enkf}}$ drops below a certain threshold, we are sure that the MCMC proposal $\theta^{*}$ will not be accepted. Hence we can save time by immediately terminating the EnKF run. Algorithm 4 describes a single iteration of the resultant MCMC algorithm, which involves reorganising the order of calculation of the acceptance probability and likelihood estimate from our standard eMCMC algorithm. This early rejection approach is employed in Section 4.5, where a computationally expensive model is studied.

4 Results

Here we demonstrate the potential of our method on several examples with different kinds of complexity. We select the number of particles $N$ and MCMC proposal variance as described in Section 2.2.1. Given that the different methods have different target distributions, we tune the random walk covariance matrix individually for each method.

In terms of accuracy we compare the approximate eMCMC and the ‘exact’ pMCMC approach visually. We note that in many applications it might not be critical to obtain samples from the exact posterior given the potential for model misspecification and/or high accuracy not being important for the analysis aims. When we deem the eMCMC approximation to be reasonable enough, we compare the statistical efficiency of the two methods using the multivariate effective sample size (ESS) of Vats et al., (2019). The overall efficiency considers the statistical efficiency and computing time simultaneously.

In Section 5 we provide suggestions on how ‘exact’ posterior sampling can be achieved whilst still using the EnKF. However, the statistical efficiency gains of these approaches will be reduced compared to eMCMC.

4.1 Population Ecology Example

4.1.1 Model and inference task

Peters et al., (2010) consider a set of competing non-linear state-space population models in ecology and apply them to several datasets. Denoting the observation at time $t$ as $y_{t}$ and the corresponding hidden state as $n_{t}$, the four models we consider are defined below:

Ricker model: $\log n_{t+1}=\log n_{t}+\beta_{0}+\beta_{1}n_{t}+\epsilon_{t}.$ 2. 2.

Theta-logistic model: $\log n_{t+1}=\log n_{t}+\beta_{0}+\beta_{2}n_{t}^{\beta_{3}}+\epsilon_{t}.$ 3. 3.

Mate-limited model: $\log n_{t+1}=2\log n_{t}+\beta_{0}+\beta_{1}n_{t}-\log(\beta_{4}+n_{t})+\epsilon_{t}.$ 4. 4.

Flexible-Allee model: $\log n_{t+1}=\log n_{t}+\beta_{0}+\beta_{1}n_{t}+\beta_{5}n_{t}^{2}+\epsilon_{t}.$

Here $\epsilon_{t}\sim\mathcal{N}(0,\sigma_{w}^{2})$. The observation process is assumed to be Gaussian, $y_{t}|n_{t}\sim\mathcal{N}(\log n_{t},\sigma_{e}^{2})$. See Peters et al., (2010) for a justification and some qualitative analyses of these models. The parameters are assumed independent a priori and have the following specifications: $\beta_{0},\beta_{1},\beta_{3},\beta_{5}\sim\mathcal{N}(0,1)$, $\beta_{4},\sigma_{w},\sigma_{e}\sim\mathcal{E}\!xp(1)$ and $\log n_{0}$ has an improper uniform prior over the real line.

Here we re-analyse the nutria dataset, a time series of female nutria abundance in East Anglia at monthly intervals, considered in Peters et al., (2010) and some references therein. The data is shown in Figure 1.

4.1.2 Inference

It is likely that all the considered models are misspecified but we would like a robust method for fitting them in order to compare the models and investigate possibilities for extending the models. We find that all models have particular difficulty in capturing the sudden drop in abundance between months 107 and 108. Further, there appears to be only small observation error. The consequence for the bootstrap filter is a very small ESS and high variance estimates of the likelihood unless a very large number of particles is used.

For eMCMC, we only require $N=250$ (Ricker, Flexible-Allee, theta-logistic) and $N=200$ (mate-limited) particles. In contrast, we use $N=50000$ for pMCMC. For some of the models, the standard deviation of the estimated log-likelihood is still larger than 1.5 even with this large number of particles. However, we find that when these occur the distribution of the log-likelihood estimator with the BPF has a skew-left distribution, which is less problematic for pMCMC getting stuck at overestimated log-likelihood values. We find that the pMCMC acceptance rates remain reasonable with $N=50000$ particles.

The MCMC acceptance rates for the four models are 15%, 4%, 11% and 10% (eMCMC), and 8%, 3%, 6% and 5% (pMCMC), respectively. The acceptance rates are lower for the theta-logistic model as the posterior distribution is far more irregular compared to the other three models (see Figure 3).

Based on Figures 2, 4 and 5 eMCMC obtains estimated univariate posterior distributions that are remarkably similar to pMCMC. There is more difference for the theta-logistic model (Figure 3) but they remain broadly similar. Further, the Monte Carlo error is greater for this model, potentially exaggerating the differences.

The two methods are compared in terms of computational efficiency on the four models in Table 1. It is evident that the eMCMC approach is producing a two order of magnitude improvement in terms of computational efficiency and still produces reasonable approximations of the posterior.

We test eMCMC for a range of $N$ values in 100-1000 and find the univariate posteriors to show little sensitivity to $N$ (results not shown). For $N=100$, the MCMC acceptance rate drops substantially and reducing $N$ further is likely to significantly reduce the statistical efficiency of MCMC due to the high-variance likelihood estimates. Therefore, it is difficult to test the sensitivity of the results to small $N$.

However, using the correlated extension (with $\sigma_{u}=0.1$) allows us to use small $N$ and maintain statistically efficient results. Similar MCMC acceptance rates as eMCMC with $250$ particles can be achieved using only $N=25$ particles. The efficiency results can be seen in Table 1. It is evident that the correlation further improves the computational efficiency in this example. The resulting approximate marginal posteriors compared to eMCMC with $N=1000$ are shown in Figures 15-18 for the four models in Appendix A. It is evident that similar approximate posteriors are obtained even with vastly different $N$ values. However, for all models there is a noticeable bias in the approximate posterior for $\sigma_{w}$. We also run the unbiased version of Section 3.3 with the correlated extension, again for $N=25$. The same figures in the appendix demonstrate that ueMCMC is able to reduce the bias in the approximate posterior for $\sigma_{w}$. The largest difference between the results for $N=1000$ and $N=25$ occurs for the theta-logistic model. For $N=25$, the unbiased version seems to offer some correction for $\theta_{w}$ and $\theta_{e}$ but produces similar results to the biased version for the other parameters. We find that with the unbiased version the ESS remains similar, but the overall efficiency is slightly reduced. The reduction in computational efficiency mainly comes here from the extra time to compute the unbiased multivariate normal density estimator (the ESS is roughly the same). We note that for applications where simulating the transition density consumes the majority of the computation, the additional time associated with computing the unbiased estimator will be significantly less noticeable.

Finally, we investigate improvements that can be obtained in this example when using the RQMC extension. Here we use $N=50$ ensemble members for each model. It is evident that the RQMC extension is producing similar marginal posteriors compared to eMCMC with $N=1000$ particles (see Appendix B). The ESS values for the four models are roughly 2500, 450, 1700, and 2100, which are competitive with standard eMCMC using a significantly larger, $N=200-250$, number of particles (see Table 1). However, the ESS/Time scores for the four models are only roughly 900, 130, 500 and 620 with the RQMC extension. Given that simulation of the transition density is trivial in this example, the cost associated with generating the RQMC samples is significant and consequently the ESS/time score with the RQMC extension is substantially reduced. However, in complex examples where simulating the transition density is expensive, the cost associated with RQMC will be far less noticeable.

4.2 Lorenz Example

4.2.1 Model and inference task

The Lorenz 63 dynamical system (Lorenz,, 1963) is a classic low dimensional example of chaotic behaviour. An Itô stochastic differential equation (SDE) version from Vrettas et al., (2015) is

[TABLE]

Here $X_{t}$ is a vector of the random variables $X_{1,t},X_{2,t},X_{3,t}$, and $W_{t}$ is a vector of three standard uncorrelated Brownian motion processes. Note that $\Sigma^{1/2}$ is interpreted as a matrix square root. We assume independent observations $Y_{i,t}\sim\mathcal{N}(X_{i,t},\sigma_{\text{obs}}^{2})$ are made at a grid of prespecified $t$ values for $i=1,2,3$.

Exact simulation of SDEs is extremely challenging, so it is common to work with an Euler-Maruyama discretisation (see e.g. Wilkinson,, 2018). For the Lorenz model above this gives,

[TABLE]

where each $z_{i+1}$ is an independent $\mathcal{N}(0,I_{3})$ realisation. Then $x_{i}$ is an approximation to $X_{t}$ for $t=i\Delta t$.

Following Vrettas et al., (2015) we simulate data from this discretised model under $\theta=(10,28,8/3)$, $\sigma_{i}^{2}=10$ for $i=1,2,3$, $\sigma_{\text{obs}}^{2}=2$ and $\Delta t=0.01$. The initial conditions are $x_{0}=(0,0,0)$. We make observations at $i=20,40,\ldots,600$, corresponding to $t=0.2,0.4,\ldots,6$. Figure 6 shows our data.

4.2.2 Log likelihoods

First we compare log-likelihood estimates produced by the EnKF and BPF. We run each method 5 times for $\theta_{1}=1,2,\ldots,20$. The other parameters are held constant at their true values. We use 100 particles for both filtering methods. The average run-times were roughly half as long for EnKF – $0.019$s – compared to BPF $0.046$s.

Figure 7 shows the results. For any $\theta_{1}$ value, the log likelihood estimates are more variable under BPF than EnKF. Variability becomes particularly large under BPF when $\theta_{1}$ is far from its true value. Such high variance is problematic in pMCMC, as it is likely to cause chains to become stuck.

Figure 7 suggests that the EnKF and BPF produce similar expected likelihood estimates when $\theta_{1}$ is close to its true value. It is hard to draw any conclusions for other $\theta_{1}$ values, as the BPF expected likelihood will be strongly driven by the upper tail of its log-likelihood estimates, and this would take a very large number of simulations to estimate well.

4.2.3 Inference

Here we assume $\sigma_{\text{obs}}$ is known, and attempt to infer $\theta_{i}$ and $\sigma_{i}$ for $i=1,2,3$. We assume these parameters have independent exponential prior distributions with rate 0.1. We ran the EnKF and BPF at the true parameter for 30 times for various choices of $N$ and calculated empirical variances. Based on these values we select $N=500$ for eMCMC and $N=2500$ for pMCMC.

We ran our algorithms targeting the log transformed parameters. Both pMCMC and eMCMC achieve acceptance rates in the range 10% to 20% indicating reasonable mixing. Trace plots also suggested good mixing, with no evidence of chains becoming stuck in the same state for a large number of iterations. The ESS values for the MCMC outputs were 390 (eMCMC) and 197 (pMCMC). Run times were 689s and 10,992s for eMCMC and pMCMC respectively. Interestingly, the pMCMC run time is roughly 15 times that of eMCMC despite using only 5 times as many particles.

Figure 8 shows the resulting marginal posterior estimates. The eMCMC posterior approximation is similar to the gold standard pMCMC results, but there are some noticable differences for some parameters e.g. the $\theta_{1}$ and $\theta_{3}$ posterior marginals are shifted downwards. Posterior correlations were small for both MCMC methods (all below $0.35$ in magnitude).

We also ran RQMC and correlated variants of eMCMC (using $\sigma_{u}=0.1$ for the latter). For RQMC, initial tuning based on variance of the log likelihood selected $N=500$, as for eMCMC. For correlated eMCMC we used a reduced number of particles, $N=100$. Posterior marginals are shown in Figure 8 and are extremely similar to eMCMC results. RQMC eMCMC produced a similar acceptance rate and ESS value (305) to eMCMC, but the cost of QMC sampling increased the run time to 3,073s (roughly a 5 times increase). Correlated eMCMC increased the acceptance rate (to 26%) and ESS (to 417) while also reducing the run time (to 195s).

4.3 Lotka Volterra Example

4.3.1 Model and inference task

The Lotka-Volterra predator-prey model (e.g. Boys et al.,, 2008) describes the continuous time evolution of the non-negative integer-values process $X_{t}=(X_{1,t},X_{2,t})^{\prime}$ where $X_{1,t}$ denotes prey and $X_{2,t}$ denotes predator. Starting from an initial value, $X_{t}$ evolves according to a Markov jump process (MJP) parameterised by stochastic rate constants $c=(c_{1},c_{2},c_{3})^{\prime}$ and characterised by the instantaneous rate or hazard function $h(x_{t},c)=(h_{1}(x_{t},c_{1}),h_{2}(x_{t},c_{2}),h_{3}(x_{t},c_{3}))^{\prime}$. Transitions over $(t,t+dt]$ take the form of one of three types (prey reproduction, prey death / predator reproduction, predator death) with associated probabilities given by

[TABLE]

The hazard function for this system is

[TABLE]

It is then relatively simple to generate realisations of this process via Gillespie’s direct method (Gillespie,, 1977), where at time $t$, the dwell time between transition events is drawn from an exponential distribution with rate $h_{0}(x_{t},c)=\sum_{i=1}^{3}h_{i}(x_{t},c_{i})$ and the transition is type $i$ with probability proportional to $h_{i}(x_{t},c_{i})$.

We assume that the MJP is observed with Gaussian error so that

[TABLE]

As all parameters of interest must be strictly positive, we consider inference for

[TABLE]

We consider two synthetic data sets ($\mathcal{D}_{1}$ and $\mathcal{D}_{2}$) simulated with rate parameters $c=(0.5,0.0025,0.3)^{\prime}$ and initial condition $x_{0}=(71,79)^{\prime}$. We further assume $\sigma_{1}=\sigma_{2}=1$ to be unknown. To allow the analysis of two data-poor scenarios, data set $\mathcal{D}_{1}$ has 51 equally spaced observations on $[0,50]$ and data set $\mathcal{D}_{2}$ is constructed by thinning $\mathcal{D}_{1}$ to give 26 equally spaced observations on $[0,25]$.

4.3.2 Inference

We compare the performance of EnKF to the gold standard auxiliary particle filter (APF) driven pMCMC scheme described in Golightly and Wilkinson, (2015). In brief, state particles are propagated using Gillespie’s direct method, with the hazard function replaced by an approximate conditioned hazard, derived from a linear Gaussian approximation to the MJP. Full details of this approach, including the calculation of the particle filter weights can be found in Golightly and Wilkinson, (2015).

We follow the practical advice given in Section 2.2.1 to choose the number of particles / ensemble members $N$ and the scaling of the innovation variance in the random walk proposal distribution. We assume independent uniform $U(-8,8)$ priors for the components of $\theta$ and ran both eMCMC and pMCMC for $10^{5}$ iterations. Since the EnKF treats the state as continuous, eMCMC used a reflecting barrier at 0 to avoid the state of the system going negative.

The results are summarised by Table 2 and Figures 9–11. We see that for both data sets, the output of eMCMC is consistent with the true values that produced the data and, more importantly, the ground truth posterior based on the output of pMCMC. For data set $\mathcal{D}_{1}$, eMCMC required more particles than pMCMC but gives better overall efficiency (as measured by the ESS per second) since sampling from the propagation construct in the auxiliary particle filter is relatively expensive. We see an increase of about a factor of 3. For data set $\mathcal{D}_{2}$, the number of particles required by pMCMC must be increased, since the propagation construct is based on a linear Gaussian approximation of the true (but unknown) hazard function of the conditioned MJP. The construct breaks down as observations are made sparsely in time (and the dynamics of the conditioned process are nonlinear between observations). Ensemble MCMC on the other hand seems to work well, requiring even fewer particles than for $\mathcal{D}_{1}$. We see an increase in overall efficiency (compared to pMCMC) of a factor of around 55.

4.4 Autoregulatory Network Example

4.4.1 Model and inference task

A commonly used mechanism for auto-regulation in prokaryotes which has been well-studied and modelled is a negative feedback mechanism whereby dimers of a protein repress its own transcription (e.g. Arkin et al.,, 1998). A simplified model for such a prokaryotic auto-regulation, based on this mechanism of dimers of a protein coded for by a gene repressing its own transcription into RNA, can be found in Golightly and Wilkinson, (2005) (see also Golightly and Wilkinson,, 2011).

Let $X_{t}=(X_{1,t},X_{2,t},X_{3,t},X_{4,t},X_{5,t})^{\prime}$ denote the number of copies of the unbound gene $X_{1,t}$, bound gene $X_{2,t}$, RNA $X_{3,t}$, protein $X_{4,t}$ and dimers of the protein $X_{5,t}$. We assume that $X_{t}$ evolves according to a Markov jump process. The possible transitions can be succinctly described by the pseudo-reaction list

[TABLE]

where, for example, occurence of $\mathcal{R}_{1}$ at time $t$ reduces $X_{1,t}$ and $X_{5,t}$ by 1, increases $X_{2,t}$ by 1, and leaves the remaining components unchanged. The associated hazard function is

[TABLE]

We consider here two challenging synthetic data sets, each consisting of 101 observations at integer times on $X_{3,t}$ (RNA) and total protein counts, $X_{4,5}+2X_{5,t}$ so that $X_{1,t}$, $X_{2,t}$, $X_{4,t}$ and $X_{5,t}$ are not observed exactly. Moreover, as in Section 4.3, we corrupt the observations by adding independent, Gaussian $\mathcal{N}\{0,\textrm{diag}(\sigma_{1}^{2},\sigma_{2}^{2})\}$ innovations to each count. We fix $\sigma_{1}=\sigma_{2}=1$ for data set $\mathcal{D}_{1}$ and $\sigma_{1}=\sigma_{2}=0$ for data set $\mathcal{D}_{2}$. Following Golightly and Wilkinson, (2005), we use the rate constants

[TABLE]

We assume that the initial condition $x_{0}=(5,5,8,8,8)^{\prime}$, the measurement error variances and the rate constants of the reversible dimerisation reactions ($c_{5}$ and $c_{6}$) are known leaving $\theta_{i}=\log c_{i}$, $i=1,2,3,4,7,8$ as the object of inference.

4.4.2 Inference

We again compare the performance of eMCMC to the gold standard auxiliary particle filter driven pMCMC scheme described in Golightly and Wilkinson, (2015). The number of particles / ensemble members $N$ was chosen as in Section 2.2.1. We assign independent Gamma $Ga(1,0.5)$ priors to each unknown rate constant and ran eMCMC and pMCMC for $2\times 10^{5}$ iterations. Note that when running eMCMC for data set $\mathcal{D}_{2}$, the values of $\sigma_{1}$ and $\sigma_{2}$ result in no shifting of the ensemble members, rendering this step ineffectual. We therefore modified eMCMC for this scenario by setting the measurement error variance to be “small” throughout the algorithm’s execution. Specifically, we found that setting $\sigma_{1}^{2}=\sigma_{2}^{2}=0.01$ gave reasonable mixing, at the expensive of introducing additional bias into the eMCMC posterior.

Table 3 and Figures 12–13 summarise the results. It is clear that eMCMC gives output that is consistent with the true values that produced the data and output from pMCMC, which exactly targets the posterior of interest. We therefore compare overall efficiency of eMCMC and pMCMC in terms of effective sample size per second, as reported in Table 3. For data set $\mathcal{D}_{1}$, eMCMC requires half the number of particles of pMCMC and gives a comparable ESS value. In terms of overall efficiency, eMCMC outperforms pMCMC by around a factor of 4. For data set $\mathcal{D}_{2}$, pMCMC requires around 2000 particles, due to the strict requirement of particle trajectories having to “hit” the observations to recieve a non zero weight. Ensemble MCMC on the other hand is able to give a comparable ESS value with just 370 particles. Consequently, for this example, eMCMC outperforms pMCMC by around a factor of 11.

4.5 Neuroscience Example

4.5.1 Model

We investigate the following realistic Neural Population Model (NPM) for brain activity. This model (see e.g. Bojak and Liley,, 2005) is known as the Liley model, and Bayesian inference for the parameters of this model has previously been described by Maybank et al., (2017). Here a high-level description of the model is presented; more detail can be found in this latter paper. The model consists of the following differential equations, where $k=e,i$ for excitatory and inhibitory contributions:

[TABLE]

where the Kronecker delta $\delta_{ek}$ admits only excitatory noise input $p$ (white noise with zero mean and fixed standard deviation) to this stochastic differential equation system, and where $S$ is a sigmoidal activation function.

We model an electroencephalogram (EEG) time-series as noisy observations of the $h_{e}$ variable of this NPM, assuming that the EEG observations are linearly proportional to $h_{e}$ with some added observational noise.

[TABLE]

where $\Delta t$ is some constant time-step and the $z_{i}$ are iid normal random variables $z_{i}\sim\mathcal{N}\left(0,\sigma^{2}\right)$ for $i=0,\ldots,n-1$.

In this paper an input noise of variance $10^{8}$ was used to simulate data, and the dynamics (consisting of 14 state variables) were simulated using the Euler-Maruyama method with step size $2.5\times 10^{-3}$. Table 4 gives the prior distributions for the parameters that were treated as unknown, giving the uniform priors that restrict the parameters to ranges found to be plausible in Bojak and Liley, (2005); other parameters were fixed to values chosen from the ranges given by Bojak and Liley, (2005).

4.5.2 Results

We compared the performance of eMCMC and pMCMC on data simulated from the Liley model for parameters that result in quasi-linear dynamics about a stable fixed point. The work in Maybank et al., (2017) suggests that the accuracy of the parameter posterior is likely to be improved by using a method that is suitable for non-linear systems (such as particle MCMC) compared to using a linearised approach such as the extended Kalman filter or the approach introduced in Maybank et al., (2017). In both MCMC approaches we use a Metropolis-Hastings approach, with a truncated multivariate normal proposal for the 11 parameters with covariance chosen to be $2.562^{2}/11$ times the estimated posterior covariance from pilot runs (this scaling being recommended by Sherlock et al.,, 2015). We considered a situation that is challenging for a particle filter, with a relatively small measurement noise of $\sigma=0.01$.

We study a simulated data set, generated from the model using Euler-Maruyama approximation. The data, shown in Figure 14 has a length of 4s and a sampling frequency of 50Hz, and was generated for the parameters $(\Gamma_{ee}=0.10631,\Gamma_{ei}=0.64105,\Gamma_{ie}=0.46477,\Gamma_{ii}=0.28663,\gamma_{ee}=291.5,\gamma_{ei}=697.76,\gamma_{ie}=458.67,\gamma_{ii}=82.33,\bar{p}_{ee}=6603.4,\bar{p}_{ei}=2625.7,\sigma=0.01)$.

We ran 40 chains of $1000$ iterations of pMCMC and eMCMC on this data, all initialised from the parameters at which the data was generated then run for an additional 500 iterations. Based on the scheme described in Section 2.2.1, we chose 1000 particles for the BPF in pMCMC, and 100 ensemble members for the EnKF in eMCMC. Both algorithms were implemented with early rejection schemes, as detailed in Section 3.4. In both cases the early rejection results in a reduction in computational cost of approximately a factor of two; with this scheme each iteration of pMCMC took an average of 1415s, compared to the average of 91s for eMCMC. The mean acceptance rate for pMCMC was $0.33\%$, compared to $0.93\%$ for eMCMC, indicating that eMCMC is (by this measure) is approximately three times as efficient whilst being more than 15 times faster. Pilot runs on longer simulated time series suggest that the efficiency of eMCMC (relative to pMCMC) improves as the length of the time series increases, but in these cases the computational cost of pMCMC was too large to permit a rigorous comparison. Kernel density estimates of the marginal posterior of each parameter are shown in Figure 14: we observe that the posteriors obtained by both methods are similar.

5 Discussion

In this paper we replace the BPF with the EnKF within a particle MCMC algorithm. We have demonstrated on a variety of examples that significant computational gains can be achieved without sacrificing much on posterior accuracy.

If exact posterior inferences are essential, there are likely to be ways to exploit our EnKF approach to improve computational performance. For example, it could be used as the cheap approximate likelihood within a delayed-acceptance MCMC algorithm (e.g. Sherlock et al.,, 2017 and Golightly et al.,, 2015) or importance sampling scheme (Franks and Vihola,, 2017). Alternatively, we might bridge our approximate posterior with the true posterior using sequential Monte Carlo. Further, our approach could be used in pilot MCMC runs to more quickly identify the regions of the parameter space with non-negligible posterior support and assist MCMC tuning generally. Particularly in the posterior tails we find that the EnKF likelihood estimator have significantly lower variance than the BPF likelihood estimator.

It is important to note that there will likely be many applications where the EnKF approximation may not be appropriate. The approach relies on being able to approximate the filtering distribution reasonably well with a Gaussian density. However, our paper illustrates that there are a wide class of models where our approach can provide reasonable accuracy. Further, Katzfuss et al., (2019) present a hierarchical approach for allowing non-Gaussian observation densities with EnKF methods, which would also be applicable to our approach.

We also did some comparisons of our approach with the particle EnKF (pEnKF) method of Katzfuss et al., (2019) (see their Algorithm 4). Briefly, their approach is a sequential algorithm that evolves a population of static parameters (particles) over time, where each particle has an associated ensemble for the hidden state. The approach uses the EnKF approximation of the likelihood to re-weight the particles. The ensemble of latent states is propagated by the transition density and the particles are propagated via resampling and jittering step. We found this algorithm to be fast. However, it requires the user to choose several aspects of the algorithm and we found the posterior approximations to be significantly less accurate then what we obtain here. However, we suggest that our approach could be incorporated into the SMC2 algorithm of Chopin et al., (2013), which uses an MCMC kernel for jittering particles and thus preserves the current target. We leave that for further research.

We did not consider posterior inference for the hidden states in this paper. It might be possible to combine our method with the ensemble Kalman smoother of van Leeuwen and Evensen, (1996), but this requires further investigation.

In this paper we compared the most commonly used particle filter (the BPF, except in the Markov jump process examples) and EnKF (the stochastic EnKF) within MCMC algorithms. In future work it would be interesting to compare extensions to both approaches. Extensions to the BPF are familiar to many in computational statistics (e.g. adaptive resampling, MCMC rejuvenation moves, the auxiliary PF) and the improvements they can bring to particle MCMC algorithms are relatively well understood. In the paper we have seen how ideas previously used in the particle MCMC context (i.e. using Quasi Monte Carlo, and the correlated approach) can also be exploited in the EnKF case. Other extensions and alternatives to the stochastic EnKF from the DA literature are also possible and have the potential to provide further improvements in the particle MCMC setting. Examples are: the deterministic EnKF (Tippett et al.,, 2003), which uses a deterministic rather than a stochastic transformation in the shift step, which may further reduce the variance of the likelihood estimates (but which may introduce further bias for nonlinear models); the equivalent weights particle filter (van Leeuwen,, 2010), which uses a deterministic transform in each step of a PF to avoid degeneracy (but which may, again, introduce bias); or for state spaces of high dimension, strategies such as variance inflation and localisation (Katzfuss et al.,, 2016).

Acknowledgments

CD and DP are grateful to RGE for providing support to present at a Bayesian workshop at the University of Reading, where this research project was instigated. RGE thanks Philip Maybank for the implementation of the solver for the Liley model, and Philip Maybank and Ingo Bojak for invaluable discussions about this model.

Appendix A - Correlated and unbiased results for the population ecology example

Appendix B - RQMC results for the population ecology example