Particle Gibbs with Ancestor Sampling

Introduction

Monte Carlo methods are one of the standard tools for inference in statistical models as they, among other things, provide a systematic approach to the problem of computing Bayesian posterior probabilities. Sequential Monte Carlo (SMC) and Markov chain Monte Carlo (MCMC) methods in particular have found application to a wide range of data analysis problems involving complex, high-dimensional models. These include state-space models (SSMs) which are used in the context of time series and dynamical systems modeling in a wide range of scientific fields. The strong assumptions of linearity and Gaussianity that were originally invoked for SSMs have indeed been weakened by decades of research on SMC and MCMC.

These methods have not, however, led to a substantial weakening of a further strong assumption, that of Markovianity. It remains a major challenge to develop efficient inference algorithms for models containing a latent stochastic process which, in contrast with the state process in an SSM, is non-Markovian. Such non-Markovian latent variable models arise in various settings, either from direct modeling or via a transformation or marginalization of an SSM. We discuss this further in Section 6; see also [5, Section 4].

In this paper we present a new tool in the family of Monte Carlo methods which is particularly useful for inference in SSMs and, importantly, in non-Markovian latent variable models. However, the proposed method is by no means limited to these model classes. We work within the framework of particle MCMC (PMCMC) which is a systematic way of combining SMC and MCMC, exploiting the strengths of both techniques. More specifically, PMCMC samplers make use of SMC to construct efficient, high-dimensional MCMC kernels with certain invariance properties. These kernels can then be used as off-the-shelf components in MCMC algorithms and other inference strategies relying on Markov kernels, such as Markovian stochastic approximation methods. PMCMC has, in a relatively short period of time, found many applications in areas such as hydrology , finance , systems biology , and epidemiology , to mention a few.

Our method builds on the particle Gibbs (PG) sampler proposed by . In PG, the aforementioned Markov kernel is constructed by running an SMC sampler in which one particle trajectory is set deterministically to a reference trajectory that is specified a priori. After a complete run of the SMC algorithm, a new trajectory is obtained by selecting one of the particle trajectories with probabilities given by their importance weights. The effect of the reference trajectory is that the resulting Markov kernel leaves its target distribution invariant, regardless of the number of particles used in the underlying SMC algorithm.

However, PG suffers from a serious drawback, which is that the mixing of the Markov kernel can be very poor when there is path degeneracy in the underlying SMC sampler . Unfortunately, path degeneracy is inevitable for high-dimensional problems, which significantly reduces the applicability of PG. This problem has been addressed in the generic setting of SSMs by adding a backward simulation step to the PG sampler, yielding a method denoted as PG with backward simulation (PGBS) . It has been found that this considerably improves mixing, making the method much more robust to a small number of particles as well as growth in the size of the data .

Unfortunately, however, the application of backward simulation is problematic for models with more intricate dependencies than in SSMs, such as non-Markovian latent variable models. The reason is that we need to consider complete trajectories of the latent process during the backward simulation pass (see Section 6 for details). The method proposed in this paper, which we refer to as particle Gibbs with ancestor sampling (PGAS), is geared toward this issue. PGAS alleviates the problem with path degeneracy by modifying the original PG kernel with a so called ancestor sampling step, thereby achieving the same effect as backward sampling, but without an explicit backward pass.

The PGAS Markov kernel is constructed in Section 3, extending the preliminary work that we have previously published in . It is also illustrated how ancestor sampling can be used to mitigate the problems with path degeneracy which deteriorates the performance of PG. In Section 4 we establish the theoretical validity of the PGAS approach, including a novel uniform ergodicity result. We then show specifically how PGAS can be used for inference and learning of SSMs and of non-Markovian latent variable models in Sections 5 and 6, respectively. The PGAS algorithm is then illustrated on several numerical examples in Section 7. As part of our development, we also propose a truncation strategy specifically for non-Markovian models. This is a generic method that is also applicable to PGBS, but, as we show in the simulation study in Section 7, the effect of the truncation error is much less severe for PGAS than for PGBS. Indeed, we obtain up to an order of magnitude increase in accuracy in using PGAS when compared to PGBS in this study. Finally, in Section 8 we conclude and point out possible directions for future work.

Sequential Monte Carlo

Let $\gamma_{\theta,t}(x_{1:t})$, for $t=1,\,\dots,\,T$, be a sequence of unnormalized densitiesWith respect to some dominating measure which we denote simply as $dx_{1:t}$. on the measurable space $(\mathsf{X}^{t},\mathcal{X}^{t})$, parameterized by $\theta\in\Theta$. Let $\bar{\gamma}_{\theta,t}(x_{1:t})$ be the corresponding normalized probability densities:

where $Z_{\theta,t}=\int\bar{\gamma}_{\theta,t}(x_{1:t})\,dx_{1:t}$ and where it is assumed that $Z_{\theta,t}>0,\forall\theta\in\Theta$. For instance, in the (important) special case of an SSM we have $\bar{\gamma}_{\theta,t}(x_{1:t})=p_{\theta}(x_{1:t}\mid y_{1:t})$, $\gamma_{\theta,t}(x_{1:t})=p_{\theta}(x_{1:t},y_{1:t})$, and $Z_{\theta,t}=p_{\theta}(y_{1:t})$. We discuss this special case in more detail in Section 5.

To draw inference about the latent variables $x_{1:T}$, as well as to enable learning of the model parameter $\theta$, a useful approach is to construct a Monte Carlo algorithm to draw samples from $\bar{\gamma}_{\theta,T}(x_{1:T})$. The sequential nature of the problem suggests the use of SMC methods; in particular, particle filters (PFs) .

We start by reviewing a standard SMC sampler, which will be used to construct the PGAS algorithm in the consecutive section. We will refer to the index variable $t$ as time, but in general it might not have any temporal meaning. Let $\{x_{1:t-1}^{i},w_{t-1}^{i}\}_{i=1}^{N}$ be a weighted particle system targeting $\bar{\gamma}_{\theta,t-1}(x_{1:t-1})$. That is, the weighted particles define an empirical point-mass approximation of the target distribution given by

This particle system is propagated to time $t$ by sampling $\{a_{t}^{i},x_{t}^{i}\}_{i=1}^{N}$ independently from a proposal kernel,

Note that $M_{\theta,t}$ depends on the complete particle system up to time $t-1$, $\{x_{1:t-1}^{i},w_{t-1}^{i}\}_{i=1}^{N}$, but for notational convenience we shall not make that dependence explicit. Here, $a_{t}^{i}$ is the index of the ancestor particle of $x_{t}^{i}$. In this formulation, the resampling step is implicit and corresponds to sampling these ancestor indices. When we write $x_{1:t}^{i}$ we refer to the ancestral path of particle $x_{t}^{i}$. That is, the particle trajectory is defined recursively as

Once we have generated $N$ ancestor indices and particles from the proposal kernel (3), the particles are weighted according to $w_{t}^{i}=W_{\theta,t}(x_{1:t}^{i})$ where the weight function is given by

for $t\geq 2$. The procedure is initialized by sampling from a proposal density $x_{1}^{i}\sim r_{\theta,1}(x_{1})$ and assigning importance weights $w_{1}^{i}=W_{\theta,1}(x_{1}^{i})$ with $W_{\theta,1}(x_{1})=\gamma_{\theta,1}(x_{1})/r_{\theta,1}(x_{1})$. The SMC sampler is summarized in Algorithm 1.

It is interesting to note that the joint law of all the random variables generated by Algorithm 1 can be written down explicitly. Let

refer to all the particles and ancestor indices, respectively, generated at time $t$ of the algorithm. It follows that the SMC sampler generates a collection of random variables $\{\mathbf{x}_{1:T},\mathbf{a}_{2:T}\}\in\mathsf{X}^{NT}\times\{1,\,\dots,\,N\}^{N(T-1)}$. Furthermore, $\{a_{t}^{i},x_{t}^{i}\}_{i=1}^{N}$ are drawn independently (conditionally on the particle system generated up to time $t-1$) from the proposal kernel $M_{\theta,t}$, and similarly at time $t=1$. Hence, the joint probability density function (with respect to a natural product of $dx$ and counting measure) of these variables is given by

The PGAS kernel

We now turn to the construction of PGAS, a family of Markov kernels on the space of trajectories $(\mathsf{X}^{T},\mathcal{X}^{T})$. We will provide an algorithm for generating samples from these Markov kernels, which are thus defined implicitly by the algorithm.

Before stating the PGAS algorithm, we review the main ideas of the PG algorithm of and we then turn to our proposed modification of this algorithm via the introduction of an ancestor sampling step.

In a standard PF, the samples $\{a_{t}^{i},x_{t}^{i}\}$ are drawn independently from the proposal kernel (3) for $i=1,\,\dots,\,N$. When sampling from the PG kernel, however, we condition on the event that the reference trajectory $x_{1:T}^{\prime}$ is retained throughout the sampling procedure. To accomplish this, we sample according to (3) only for $i=1,\,\dots,\,N-1$. The $N$th particle and its ancestor index are then set deterministically as $x_{t}^{N}=x_{t}^{\prime}$ and $a_{t}^{N}=N$. This implies that after a complete pass of the algorithm, the $N$th particle path coincides with the reference trajectory, i.e., $x_{1:T}^{N}=x_{1:T}^{\prime}$.

The fact that $x_{1:T}^{\prime}$ is used as a reference trajectory in the SMC sampler implies an invariance property of the PG kernel which is of key relevance. More precisely, as show by [6, Theorem 5], for any number of particles $N\geq 1$ and for any $\theta\in\Theta$, the PG kernel leaves the exact target distribution $\bar{\gamma}_{\theta,T}$ invariant. We return to this invariance property below, when it is shown to hold also for the proposed PGAS kernel.

2 Ancestor sampling

As noted above, the PG algorithm keeps the reference trajectory $x_{1:T}^{\prime}$ intact throughout the sampling procedure. While this results in a Markov kernel which leaves $\bar{\gamma}_{\theta,T}$ invariant, it has been recognized that the mixing properties of this kernel can be very poor due to path degeneracy .

To address this fundamental problem we now turn to our new procedure, PGAS. The idea is to sample a new value for the index variable $a_{t}^{N}$ in an ancestor sampling step. While this is a small modification of the algorithm, the improvement in mixing can be quite considerable; see Section 3.3 and the numerical evaluation in Section 7. The ancestor sampling step is implemented as follows.

At time $t\geq 2$, we consider the part of the reference trajectory $x_{t:T}^{\prime}$ ranging from the current time $t$ to the final time point $T$. The task is to artificially assign a history to this partial path. This is done by connecting $x_{t:T}^{\prime}$ to one of the particles $\{x_{1:t-1}^{i}\}_{i=1}^{N}$. Recall that the ancestry of a particle is encoded via the corresponding ancestor index. Hence, we can connect the partial reference path to one of the particles $\{x_{1:t-1}^{i}\}_{i=1}^{N}$ by assigning a value to the variable $a_{t}^{N}\in\{1,\,\dots,\,N\}$. To do this, first we compute the weights

The sampling procedure outlined above is summarized in Algorithm 2 and the family of PGAS kernels is formally defined below. Note that the only difference between PG and PGAS is on line 8 of Algorithm 2 (where, for PG, we would simply set $a_{t}^{N}=N$). However, as we shall see, the effect of this small modification on the mixing of the kernel is quite significant.

For any $N\geq 1$ and any $\theta\in\Theta$, Algorithm 2 maps $x_{1:T}^{\prime}$ stochastically into $x_{1:T}^{\star}$, thus implicitly defining a Markov kernel $P_{\theta}^{N}$ on $(\mathsf{X}^{T},\mathcal{X}^{T})$. The family of Markov kernels $\{P_{\theta}^{N}:\theta\in\Theta\}$, indexed by $N\geq 1$, is referred to as the PGAS family of kernels.

3 The effect of path degeneracy on PG and on PGAS

We have argued that ancestor sampling can considerably improve the mixing of PG. To illustrate this effect and to provide an explanation of its cause, we consider a simple numerical example. Further empirical evaluation of PGAS is provided in Section 7. Consider the stochastic volatility model,

where the state process $\{x_{t}\}_{t\geq 1}$ is latent and observations are made only via the measurement process $\{y_{t}\}_{t\geq 1}$. Similar models have been used to generalize the Black-Scholes option pricing equation to allow for the variance to change over time .

For simplicity, the parameter $\theta=(a,\sigma)=(0.9,0.5)$ is assumed to be known. A batch of $T=400$ observations are simulated from the system. Given these, we seek the joint smoothing density $p(x_{1:T}\mid y_{1:T})$. To generate samples from this density we employ both PG and PGAS with varying number of particles ranging from $N=5$ to $N=1\thinspace 000$. We simulate sample paths of length $1\thinspace 000$ for each algorithm. To compare the mixing, we look at the update rate of $x_{t}$ versus $t$, which is defined as the proportion of iterations where $x_{t}$ changes value. The results are reported in Figure 1, which reveals that ancestor sampling significantly increases the probability of updating $x_{t}$ for $t$ far from $T$.

The poor update rates for PG is a manifestation of the well known path degeneracy problem of SMC samplers (see, e.g., ). Consider the process of sampling from the PG kernel for a fixed reference trajectory $x_{1:T}^{\prime}$. A particle system generated by the PG algorithm (corresponding to Algorithm 2, but with line 8 replaced with $a_{t}^{N}=N$) is shown in Figure 2 (left). For clarity of illustration, we have used a small number of particles and time steps, $N=20$ and $T=50$, respectively. By construction the reference trajectory (shown by a thick blue line) is retained throughout the sampling procedure. As a consequence, the particle system degenerates toward this trajectory which implies that $x_{1:T}^{\star}$ (shown as a red line) to a large extent will be identical to $x_{1:T}^{\prime}$.

What is, perhaps, more surprising is that PGAS is so much more insensitive to the degeneracy issue. To understand why this is the case, we analyze the procedure for sampling from the PGAS kernel $P_{\theta}^{N}(x_{1:T}^{\prime},\cdot)$ for the same reference trajectory $x_{1:T}^{\prime}$ as above. The particle system generated by Algorithm 2 (with ancestor sampling) is shown in Figure 2 (right). The thick blue lines are again used to illustrate the reference particles, but now with updated ancestor indices. That is, the blue line segments are drawn between $x_{t-1}^{a_{t}^{N}}$ and $x_{t}^{\prime}$ for $t\geq 2$. It can be seen that the effect of ancestor sampling is that, informally, the reference trajectory is broken into pieces. It is worth pointing out that the particle system still collapses; ancestor sampling does not prevent path degeneracy. However, it causes the particle system to degenerate toward something different than the reference trajectory. As a consequence, $x_{1:T}^{\star}$ (shown as a red line in the figure) will with high probability be substantially different from $x_{1:T}^{\prime}$, enabling high update rates and thereby much faster mixing.

Theoretical justification

We begin by stating a theorem, whose proof is provided later in this section, which shows that the invariance property of PG is not violated by the ancestor sampling step.

For any $N\geq 1$ and $\theta\in\Theta$, the PGAS kernel $P_{\theta}^{N}$ leaves $\bar{\gamma}_{\theta,T}$ invariant:

An apparent difficulty in establishing this result is that it is not possible to write down a simple, closed-form expression for $P_{\theta}^{N}$. In fact, the PGAS kernel is given by

Recall that the particle trajectory $x_{1:T}^{k}$ is the ancestral path of the particle $x_{T}^{k}$. That is, we can write

and similarly for the ancestor indices. By construction, $\phi_{\theta}$ is nonnegative and integrates to one, i.e., $\phi_{\theta}$ is a probability density function on $\Omega$. We refer to this density as the extended target density.

The factorization into a marginal and a conditional density is intended to reveal some of the structure inherent in the extended target density. In particular, the marginal density of the variables $\{x_{1:T}^{b_{1:T}},b_{1:T}\}$ is defined to be equal to the original target density $\bar{\gamma}_{\theta,T}(x_{1:T}^{b_{1:T}})$, up to a factor $N^{-T}$ related to the index variables $b_{1:T}$. This has the important implication that if $\{\mathbf{x}_{1:T},\mathbf{a}_{2:T},k\}$ are distributed according to $\phi_{\theta}$, then, by construction, the marginal distribution of $x_{1:T}^{b_{1:T}}$ is $\bar{\gamma}_{\theta,T}$.

By constructing an MCMC kernel with invariant distribution $\phi_{\theta}$, we will thus obtain a kernel with invariant distribution $\bar{\gamma}_{\theta,T}$ (the PGAS kernel) as a byproduct. To prove Theorem 1 we will reinterpret all the steps of the PGAS algorithm as partially collapsed Gibbs steps for $\phi_{\theta}$. The meaning of partial collapsing will be made precise in the proof of Lemma 2 below, but basically it refers to the process of marginalizing out some of the variables of the model in the individual steps of the Gibbs sampler. This is done in such a way that it does not violate the invariance property of the Gibbs kernel, i.e., each such Gibbs step will leave the extended target distribution invariant. As a consequence, the invariance property of the PGAS kernel follows. First we show that the PGAS algorithm in fact implements the following sequence of partially collapsed Gibbs steps for $\phi_{\theta}$.

Given $x_{1:T}^{\prime,b_{1:T}^{\prime}}\in\mathsf{X}^{T}$ and $b_{1:T}^{\prime}\in\{1,\,\dots,\,N\}^{T}$:

Draw $\mathbf{x}_{1}^{-b_{1}^{\prime}}\sim\phi_{\theta}(\,\cdot\mid x_{1:T}^{\prime,b_{1:T}^{\prime}},b_{1:T}^{\prime})$ and, for $t=2$ to $T$, draw:

Draw $k\sim{}\phi_{\theta}(\,\cdot\mid\mathbf{x}_{1:T}^{-b_{1:T}},\mathbf{a}_{2:T},x_{1:T}^{\prime,b_{1:T}^{\prime}})$.

Algorithm 2 is equivalent to the partially collapsed Gibbs sampler of Procedure 1, conditionally on $x_{1:T}^{\prime,b_{1:T}^{\prime}}=x_{1:T}^{\prime}$ and $b_{1:T}^{\prime}=(N,\,\dots,\,N)$.

By marginalizing this expression over $\{\mathbf{x}_{t+1:T}^{-b_{t+1:T}},\mathbf{a}_{t+1:T}^{-b_{t+1:T}}\}$ we get

Hence, we can sample from (14a) and (14b) by drawing $x_{1}^{i}\sim r_{\theta,1}(\cdot)$ for $i\in\{1,\,\dots,\,N\}\setminus b_{1}$ and $\{a_{t}^{i},x_{t}^{i}\}\sim M_{\theta,t}(\cdot)$ for $i\in\{1,\,\dots,\,N\}\setminus b_{t}$, respectively. Consequently, with the choice $b_{t}=N$ for $t=1,\,\dots,\,T$, the initialization at line 1 and the particle propagation at line 5 of Algorithm 2 correspond to sampling from (14a) and (14b), respectively.

Next, we consider the ancestor sampling step. Recall that $a_{t}^{b_{t}}$ identifies to $b_{t-1}$. We can thus write

To simplify this expression, note first that we can write

By using the definition of the weight function (5), this expression can be expanded according to

Plugging the trajectory $x_{1:t-1}^{b_{1:t-1}}$ into the above expression, we get

Expanding the numerator in (15) according to (18) results in

Consequently, with $b_{t}=N$ and $x_{t:T}^{b_{t:T}}=x_{t:T}^{\prime}$, sampling from (19) corresponds to the ancestor sampling step of line 8 of Algorithm 2. Finally, analogously to (19), it follows that $\phi_{\theta}(k\mid\mathbf{x}_{1:T},\mathbf{a}_{2:T})\propto w_{T}^{k}$, which corresponds to line 12 of Algorithm 2. ∎

Next, we show that Procedure 1 leaves $\phi_{\theta}$ invariant. This is done by concluding that the procedure is a properly collapsed Gibbs sampler; see . Marginalization, or collapsing, is commonly used within Gibbs sampling to improve the mixing and/or to simplify the sampling procedure. However, it is crucial that the collapsing is carried out in the correct order to respect the dependencies between the variables of the model.

The Gibbs sampler of Procedure 1 is properly collapsed and thus leaves $\phi_{\theta}$ invariant.

Consider the following sequence of complete Gibbs steps:

Draw $\{\mathbf{x}_{1}^{-b_{1}^{\prime}},\underline{\mathbf{x}}_{2:T}^{-b_{2:T}^{\prime}},\underline{\mathbf{a}}_{2:T}^{-b_{2:T}^{\prime}}\}\sim\phi_{\theta}(\,\cdot\mid x_{1:T}^{\prime,b_{1:T}^{\prime}},b_{1:T}^{\prime})$ and, for $t=2$ to $T$, draw:

Draw $k\sim{}\phi_{\theta}(\,\cdot\mid\mathbf{x}_{1:T}^{-b_{1:T}^{\prime}},\mathbf{a}_{2:T},x_{1:T}^{\prime,b_{1:T}^{\prime}})$.

In the above, all the samples are drawn from conditionals under the full joint density $\phi_{\theta}(\mathbf{x}_{1:T},\mathbf{a}_{2:T},k)$. Hence, it is clear that the above procedure will leave $\phi_{\theta}$ invariant. Note that some of the variables above have been marked by an underline. It can be seen that these variables are in fact never conditioned upon in any subsequent step of the procedure. That is, the underlined variables are never used. Therefore, to obtain a valid sampler it is sufficient to sample all the non-underlined variables from their respective marginals. Furthermore, from (14b) it can be seen that $\{\mathbf{x}_{t}^{-b_{t}},\mathbf{a}_{t}^{-b_{t}}\}$ are conditionally independent of $a_{t}^{b_{t}}$, i.e., it follows that the complete Gibbs sweep above is equivalent to the partially collapsed Gibbs sweep of Procedure 1. Hence, the Gibbs sampler is properly collapsed and it will therefore leave $\phi_{\theta}$ invariant. ∎

denote the law of the random variables generated by Procedure 1, conditionally on $x_{1:T}^{\prime,b_{1:T}^{\prime}}=x_{1:T}^{\prime}$ and on $b_{1:T}^{\prime}$. Using Lemma 2 and recalling that $\phi_{\theta}(x_{1:T}^{b_{1:T}},b_{1:T})=N^{-T}\bar{\gamma}_{\theta,T}(x_{1:T}^{b_{1:T}})$ we have

By Lemma 1 we know that Algorithm 2, which implicitly defines $P_{\theta}^{N}$, is equivalent to Procedure 1 conditionally on $x_{1:T}^{\prime,b_{1:T}^{\prime}}=x_{1:T}^{\prime}$ and $b_{1:T}^{\prime}=(N,\,\dots,\,N)$. That is to say,

However, the law of $x_{1:T}^{\star}$ in Algorithm 2 is invariant to permutations of the particle indices. That is, it does not matter if we place the reference particles on the $N$th positions, or on some other positions, when enumerating the particlesA formal proof of this statement is given for the PG sampler in . The same argument can be used also for PGAS.. This implies that for any $b_{1:T}^{\prime}\in\{1,\,\dots,\,N\}^{T}$,

Plugging (23) into (21) gives the desired result,

2 Ergodicity

To show ergodicity of the PGAS kernel we need to characterize the support of the target and the proposal densities. Let,

with obvious modifications for $t=1$. The following is a minimal assumption.

For any $\theta\in\Theta$ and $t\in\{1,\,\dots,\,T\}$ we have $\mathcal{S}_{t}^{\theta}\subseteq\mathcal{Q}_{t}^{\theta}$.

Assumption (A4.2) basically states that the support of the proposal density should cover the support of the target density. Ergodicity of PG under Assumption (A4.2) has been established by Andrieu et al. . The same argument can be applied also to PGAS.

Assume (A4.2). Then, for any $N\geq 2$ and $\theta\in\Theta$, $P_{\theta}^{N}$ is $\bar{\gamma}_{\theta,T}$-irreducible and aperiodic. Consequently,

To strengthen the ergodicity results for the PGAS kernel, we use a boundedness condition for the importance weights, given in assumption (A4.2) below. Such a condition is typical also in classical importance sampling and is, basically, a slightly stronger version of assumption (A4.2).

For any $\theta\in\Theta$ and $t\in\{1,\,\dots,\,T\}$, there exists a constant $\kappa_{\theta}<\infty$ such that $\|W_{\theta,t}\|_{\infty}\leq\kappa_{\theta}$.

Assume (A4.2). Then, for any $N\geq 2$ and $\theta\in\Theta$, $P_{\theta}^{N}$ is uniformly ergodic. That is, there exist constants $R_{\theta}<\infty$ and $\rho_{\theta}\in[0,1)$ such that

We show that $P_{\theta}^{N}$ satisfies a Doeblin condition,

for some constant $\varepsilon_{\theta}>0$. Uniform ergodicity then follows from [19, Proposition 2]. To prove (25) we use the representation of the PGAS kernel in (9),

Here, the inequality follows from bounding the weights in the normalization by $\kappa_{\theta}$ and by simply discarding the $N$th term of the sum (which is clearly nonnegative). The last equality follows from the fact that the particle trajectories $\{x_{1:T}^{i}\}_{i=1}^{N-1}$ are equally distributed under Algorithm 2.

where the inequality follows analogously to (26). Now, let

Then, by iteratively making use of (27) and changing the order of integration, we can bound (26) according to

With $N\geq 2$ and since $Z_{\theta,T}>0$ the result follows. ∎

PGAS for state-space models

SSMs comprise an important special case of the model class treated above. In this section, we illustrate how PGAS can be used for inference and learning of these models. We consider the nonlinear/non-Gaussian SSM

and $x_{1}\sim\mu_{\theta}(x_{1})$, where $\theta\in\Theta$ is a static parameter, $x_{t}$ is the latent state and $y_{t}$ is the observation at time $t$, respectively. Given a batch of measurements $y_{1:T}$, we wish to make inferences about $\theta$ and/or about the latent states $x_{1:T}$.

Consider first the Bayesian setting where a prior distribution $\pi(\theta)$ is assigned to $\theta$. We seek the parameter posterior $p(\theta\mid x_{1:T})$ or, more generally, the joint state and parameter posterior ${p(\theta,x_{1:T}\mid y_{1:T})}$. Gibbs sampling can be used to simulate from this distribution by sampling the state variables $\{x_{t}\}$ one at a time and the parameters $\theta$ from their respective posteriors. However, it has been recognized that this can result in poor mixing, due to the often high autocorrelation of the state sequence. The PGAS kernel offers a different approach, namely to sample the complete state trajectory $x_{1:T}$ in one block. This can considerably improve the mixing of the sampler . Due to the invariance property of the kernel (Theorem 1), the validity of the Gibbs sampler is not violated. We summarize the procedure in Algorithm 3.

PGAS is also useful for maximum-likelihood-based learning of SSMs. A popular strategy for computing the maximum likelihood estimator

is to use the expectation maximization (EM) algorithm . EM is an iterative method, which maximizes $\log p_{\theta}(y_{1:T})$ by iteratively maximizing an auxiliary quantity: $\theta[n]=\operatorname*{arg\,max}_{\theta\in\Theta}Q(\theta,\theta[n-1])$, where

When the above integral is intractable to compute, one can use a Monte Carlo approximation or a stochastic approximation of the intermediate quantity, leading to the MCEM and the SAEM algorithms, respectively. When the underlying Monte Carlo simulation is computationally involved, SAEM is particularly useful since it makes efficient use of the simulated values. The SAEM approximation of the auxiliary quantity is given by

where $\alpha_{n}$ is the step size and, in the vanilla form of SAEM, $x_{1:T}[n]$ is drawn from the joint smoothing density $p_{\theta[n-1]}(x_{1:T}\mid y_{1:T})$. In practice, the stochastic approximation update (32) is typically made on some sufficient statistic for the complete data log-likelihood; see for details. While the joint smoothing density is intractable for a general nonlinear/non-Gaussian SSM, it has been recognized that it is sufficient to sample from a uniformly ergodic Markov kernel, leaving the joint smoothing distribution invariant . A practical approach is therefore to compute the auxiliary quantity according to the stochastic approximation (32), but where $x_{1:T}[n]$ is simulated from the PGAS kernel $P_{\theta{[n-1]}}^{N}(x_{1:T}[n-1],\cdot)$. This particle SAEM algorithm, previously presented in , is summarized in Algorithm 4.

2 Sampling from the PGAS kernel

Sampling from the PGAS kernel, i.e., running Algorithm 2, is similar to running a PF. The only non-standard (and nontrivial) operation is the ancestor sampling step. For the learning algorithms discussed above, the distribution of interest is the joint smoothing distribution. Consequently, the unnormalized target density is given by $\gamma_{\theta,T}(x_{1:T})=p_{\theta}(x_{1:T},y_{1:T})$. The ancestor sampling weights in (7) are thus given by

This expression can be understood as an application of Bayes’ theorem. The importance weight $w_{t-1}^{i}$ is the prior probability of the particle $x_{t-1}^{i}$ and the factor $f_{\theta}(x_{t}^{\prime}\mid x_{t-1}^{i})$ is the likelihood of moving from $x_{t-1}^{i}$ to $x_{t}^{\prime}$. The product of these two factors is thus proportional to the posterior probability that $x_{t}^{\prime}$ originated from $x_{t-1}^{i}$.

and take $x_{1:T}^{\star}=x_{1:T}^{j_{1:T}}$ as the output from the algorithm. In the above, the conditioning on the forward particle system $\{\mathbf{x}_{1:T},\mathbf{a}_{2:T}\}$ is implicit.

Let the Markov kernel on $(\mathsf{X}^{T},\mathcal{X}^{T})$ defined by this procedure be denoted as $P_{\mathsf{BS},\theta}^{N}$. An interesting question to ask is whether or not the PGAS kernel $P_{\theta}^{N}$ and the PGBS kernel $P_{\mathsf{BS},\theta}^{N}$ are probabilistically equivalent. It turns out that, in some specific scenarios, this is indeed the case.

Assume that PGAS and PGBS both target the joint smoothing distribution for an SSM and that both methods use the bootstrap proposal kernel in the internal particle filters, i.e., $r_{\theta,t}(x_{t}\mid x_{1:t-1})=f_{\theta}(x_{t}\mid x_{t-1})$. Then, for any $x_{1:T}^{\prime}\in\mathsf{X}^{T}$ and $B\in\mathcal{X}^{T}$, $P_{\theta}^{N}(x_{1:T}^{\prime},B)=P_{\mathsf{BS},\theta}^{N}(x_{1:T}^{\prime},B)$.

Proposition 1 builds upon [30, Proposition 5], where the equivalence between a (standard) bootstrap PF and a backward simulator is established. In Appendix A, we adapt their argument to handle the case with conditioning on a reference trajectory and ancestor sampling. The conditions of Proposition 1 imply that the weight functions (5) in the internal particle filters are independent of the ancestor indices. This is key in establishing the above result and we emphasize that the equivalence between the samplers does not hold in general for models outside the class of SSMs. In particular, for the class of non-Markovian latent variable models, discussed in the subsequent section, we have found that the samplers have quite different properties.

PGAS for non-Markovian models

A very useful generalization of SSMs is the class of non-Markovian latent variable models,

Similarly to the SSM (29), this model is characterized by a latent process $x_{t}\in\mathsf{X}$ and an observed process $y_{t}\in\mathsf{Y}$. However, it does not share the conditional independence properties that are central to SSMs. Instead, both the transition density $f_{\theta}$ and the measurement density $g_{\theta}$ may depend on the entire past history of the latent process. In Sections 6.2 and 6.3, we discuss the ancestor sampling step of the PGAS algorithm specifically for these non-Markovian models. We consider two approaches for efficient implementation of this step, first by using Metropolis-Hastings within PGAS and then by using a truncation strategy for the ancestor sampling weights. First, however, to motivate the present development we review some application areas in which this type of models arise.

In Bayesian nonparametrics the latent random variables of the classical Bayesian model are replaced by latent stochastic processes, which are typically non-Markovian. This includes popular models based on the Dirichlet process, e.g., , and Gaussian process regression and classification models . These processes are also commonly used as components in hierarchical Bayesian models, which then inherit their non-Markovianity. An example is the Gaussian process SSM , a flexible nonlinear dynamical systems model, for which PGAS has been successfully applied .

Another typical source of non-Markovianity is by marginalization over part of the state vector (i.e., Rao-Blackwellization ) or by a change of variables in an SSM. This type of operations typically results in a loss of the Markov property, but they can, however, be very useful. For instance, by expressing an SSM in terms of its “innovations” (i.e., the driving noise of the state process), it is possible to use backward and ancestor sampling in models for which the state transition density is not available to us. This includes many models for which the transition is implicitly given by a simulator or degenerate models where the transition density does not even exist . We illustrate these ideas in Section 7. See also [5, Section 4] for a more in-depth discussion on reformulations of SSMs as non-Markovian models.

Finally, it is worth to point out that many statistical models which are not sequential “by nature” can be conveniently viewed as non-Markovian latent variable models. This includes, among others, probabilistic graphical models such as Markov random fields; see [5, Section 4].

2 Forced move Metropolis-Hastings

To employ PGAS (or in fact any backward-simulation-based method; see ) we need to evaluate the ancestor sampling weights (7) which depend on the ratio,

Assuming that $g_{\theta}$ and $f_{\theta}$ can both be evaluated in constant time, the computational cost of computing the backward sampling weights (7) will thus be $O(NT)$. This step can easily become the computational bottleneck when applying the PGAS algorithm to a non-Markovian model.

A simple way to reduce the complexity is to employ Metropolis-Hastings (MH) within PGAS. Let

denote the law of the ancestor index $a_{t}^{N}$, sampled at line 8 of Algorithm 2. From Lemma 1, we know that this step of the algorithm in fact corresponds to a Gibbs step for the extended target distribution (11). To retain the correct limiting distribution of the PGAS kernel, it is therefore sufficient that $a_{t}^{N}$ is sampled from a Markov kernel leaving (37) invariant (resulting in a standard combination of MCMC kernels; see, e.g., ).

Let $q(k^{\prime}\mid k)$ be an MH proposal kernel on $\{1,\,\dots,\,N\}$. We can thus propose a move for the ancestor index $a_{t}^{N}$, from $N$ to $k^{\prime}$, by simulating $k^{\prime}\sim q(\cdot\mid N)$. With probability

the sample is accepted and we set $a_{t}^{N}=k^{\prime}$, otherwise we keep the ancestry $a_{t}^{N}=N$. Using this approach, we only need to evaluate the ancestor sampling weights for the proposed values, bringing the total computational cost down from $O(NT^{2})$ to $O(NT+T^{2})$. While still quadratic in $T$, this reduction can be very useful whenever $N$ is moderately large.

Since the variable $a_{t}^{N}$ is discrete-valued, it is recommended to use a forced move proposal in the spirit of . That is, $q$ is constructed so that $q(k\mid k)=0,\forall k$, ensuring that the current state of the chain is not proposed anew, which would be a wasteful operation. One simple choice is to let $q(k^{\prime}\mid k)$ be uniform over $\{1,\,\dots,\,N\}\setminus k$. In the subsequent section, we discuss a different strategy for reducing the complexity of the ancestor sampling step, which can also be used to design a better proposal for the forced move MH sampler.

3 Truncation of the ancestor sampling weights

Numerical evaluation

In this section we illustrate the properties of PGAS in a simulation study. First, in Section 7.1 we consider a simple linear Gaussian SSM and investigate the improvement in mixing offered by ancestor sampling when PGAS is compared with PG. We do not consider PGBS in this example since, by Proposition 1, PGAS and PGBS are probabilistically equivalent in this scenario.

When applied to non-Markovian models, however, Proposition 1 does not apply since the weight function will depend on the complete history of the particles. PGAS and PGBS will then have different properties as is illustrated empirically in Section 7.2 where we consider inference in degenerate SSMs reformulated as non-Markovian models. Finally, in Section 7.3 we use a similar reformulation and apply PGAS for identification of an epidemiological model for which the transition kernel is not available.

Consider a first-order linear Gaussian state-space (LGSS) model,

with initial state $p(x_{1})=\mathcal{N}(x_{1};0,q/(1-a^{2}))$ and unknown parameters $\theta=(a,q,r)$. This system is used as a proof of concept to illustrate the superior mixing of PGAS when compared to PG. For this system it is possible to implement an ideal Gibbs sampler, i.e., by iteratively sampling from the posterior parameter distribution $p(\theta\mid x_{1:T},y_{1:T})$ and from the full joint smoothing distribution $p_{\theta}(x_{1:T}\mid y_{1:T})$. This is useful for comparison, since the ideal Gibbs sampler is the baseline for both PG samplers.

To further investigate the robustness of PGAS we repeat the same experiment with a larger data batch consisting of $T=2\thinspace 000$ samples. The results are given in Figure 4 (bottom row). The effect can be seen even more clearly in this more challenging scenario. The big difference in mixing between the two samplers can be understood as a manifestation of how they are affected by path degeneracy. These results are in agreement with the discussion in Section 3.3.

2 Degenerate LGSS models

Many dynamical systems are most naturally modeled as degenerate in the sense that the transition kernel of the state-process does not admit any density with respect to a dominating measure. It is problematic to use (particle-filter-based) backward sampling methods for these models, owing to the fact that the backward kernel of the state process will also be degenerate. As a consequence, it is not possible to approximate the backward kernel using the forward filter particles.

To illustrate how this difficulty can be remedied by a change of variables, consider an LGSS model of the form

Since the Gaussian process noise enters only on the first part of the state vector (or, equivalently, the process noise covariance matrix is rank deficient) the state transition kernel is degenerate. However, for the same reason, the state component $z_{t}$ is $\sigma(x_{1:t})$-measurable and we can write $z_{t}=z_{t}(x_{1:t})$. Therefore, it is possible to rephrase (41) as a non-Markovian model with latent process given by $\{x_{t}\}_{t\geq 1}$.

As a first illustration, we simulate $T=200$ samples from a fourth-order, single output system with poles at $-0.65$, $-0.12$, and $0.22\pm 0.10i$. We let $\dim(x_{t})=1$ and $Q=R=0.1$. For simplicity, we assume that the system parameters are known and seek the joint smoothing distribution $p(x_{1:T}\mid y_{1:T})$. In the non-Markovian formulation it is possible to apply backward-simulation-based methods, such as PGAS and PGBS, as described in Section 6. The problem, however, is that the non-Markovianity gives rise to an $O(T^{2})$ computational complexity. To obtain more practical inference algorithms we employ the weight truncation strategy (39).

Again, we compute the posterior means of $x_{1:T}$ (discarding $1\thinspace 000$ samples) and RMSE values relative the true posterior mean. Box plots over the different systems are shown in Figure 6. Since the process noise only enters on one of the state components, the mixing tends to deteriorate as we increase the model order. Figure 3 shows how the probability distributions on $\{1,\,\dots,\,N\}$ change as we increase the truncation level, in two representative cases for a 5th and a 20th order system, respectively. By using an adaptive level, we can obtain accurate results for systems of different dimensions, without having to change any settings between the runs.

3 Epidemiological model

As a final numerical illustration, we consider identification of an epidemiological model using PGAS. Seasonal influenza epidemics each year cause millions of severe illnesses and hundreds of thousands of deaths world-wide . Furthermore, new strains of influenza viruses can possibly cause pandemics with very severe effects on the public health. The ability to accurately predict disease activity can enable early response to such epidemics, which in turn can reduce their impact.

We consider a susceptible/infected/recovered (SIR) model with environmental noise and seasonal fluctuations . The model, specified by a stochastic differential equation, is discretized according to the Euler-Maruyama method, yielding

where $v_{t}\sim\mathcal{N}(0,1/\sqrt{dt})$ and $dt$ is the sampling time. Here, $S_{t}$, $I_{t}$ and $R_{t}$ represent the number of susceptible, infected and recovered individuals at time $t$ (months), respectively. The total population size $\mathcal{P}=10^{6}$ and the host birth/death rate $\mu=0.0012$ are assumed known. The seasonally varying transmission rate is given by $\beta_{t}=R_{0}(\gamma+\mu)(1+\alpha\sin(2\pi t/12))$ where $R_{0}$ is the basic reproductive ratio, $\gamma$ is the rate of recovery and $\alpha$ is the strength of seasonality.

Furthermore, we consider an observation model which is inspired by the Google Flu Trends project . The idea is to use the frequency of influenza-related search engine queries to infer knowledge of the dynamics of the epidemic. Let $Q_{k}$ be the proportion of influenza-related queries counted during a time interval $(\Delta(k-1),\Delta k]$. Following , we use a linear relationship between the log-odds of the relative query counts and the log-odds of the proportion of infected individuals,

where $\bar{I}_{k}$ is the mean value of $I_{t}$ during the time interval $(\Delta(k-1),\Delta k]$ and $\operatorname{logit}(p)=\log(p/(1-p))$. As in we consider weekly query counts, i.e., $\Delta=7/30$ (assuming for simplicity that we have 30 days in each month). Using this value of $\Delta$ as sampling time will, however, result in an overly large discretization errors. Instead, we sample the model (42) $m=7$ times per week: $dt=\Delta/m$.

In , the particle marginal MH sampler is used to identify a similar SIR model, though with a different observation model. A different Monte Carlo strategy, based on a particle filter with an augmented state space, for identification of an SIR model is proposed in . We suggest to use the PGAS algorithm for joint state and parameter inference in the model (42)–(43). However, there are two difficulties in applying PGAS directly to this model. Firstly, the transition kernel of the state process, as defined between consecutive observation time points $\Delta(k-1)$ and $\Delta k$, is not available in closed form. Secondly, since the state is three-dimensional, whereas the driving noise $v_{t}$ is scalar, the transition kernel is degenerate. To cope with these difficulties we (again) suggest collapsing the model to the driving noise variables. Let $V_{k}=\begin{pmatrix}v_{\Delta(k-1)}&v_{\Delta(k-1)+dt}&\cdots&v_{\Delta k-dt}\end{pmatrix}^{T}$. It follows that the model (42)–(43) can be equivalently expressed as the non-Markovian latent variable model,

for some likelihood function $g_{\theta}$ (see (46)). A further motivation for using this reformulation is that the latent variables $V_{k}$ are a priori independent of the model parameters $\theta$. This can result in a significant improvement in mixing of the Gibbs sampler, in particular when there are strong dependencies between the system state the parameters .

We generate 8 years of data with weekly observations. The number of infected individuals $I_{t}$ over this time period is shown in Figure 8. The first half of the data batch is used for estimation of the model parameters. We run PGAS with $N=10$ for $50\thinspace 000$ iterations (discarding the first $10\thinspace 000$). For sampling the system parameters $(\gamma,R_{0},\alpha,F)$, we use Metropolis-Hastings steps with a Gaussian random walk proposal, tuned according to an initial trial run. For $(\rho,\sigma^{2})$, we exploit the conjugacy of the normal-inverse-gamma prior to the likelihood (43) and sample the variables from their true posterior. The innovation variables $V_{1:T}$ are sampled from the PGAS kernel by Algorithm 2 (no truncation is used for the ancestor sampling weights). Since the latter step is the computational bottleneck of the algorithm, we execute ten MH steps for $\theta$, for each draw from the PGAS kernel.

It is worth pointing out that, while the sampler effectively targets the collapsed model (44), it is most straightforwardly implemented using the original state variables from (42). With $x_{k}=(S_{\Delta k},I_{\Delta k},R_{\Delta k})^{\mathsf{T}}$ we can simulate $x_{k+1}$ given $x_{k}$ according to (42) which is used in the underlying particle filter. The innovation variables $V_{k}$ need only be taken into account for the ancestor sampling step. Let $V_{1:T}^{\prime}$ be the reference innovation trajectory. To compute the ancestor sampling weights (7) we need to evaluate the ratios,

Using (43), the observation likelihood can be written as

Histograms representing the estimated posterior parameter distributions are shown in Figure 7. As can be seen, the true system parameters fall well within the credible regions. Finally, the identified model is used to make one-month-ahead predictions of the disease activity for the subsequent four years, as shown in Figure 8. The predictions are made by sub-sampling the Markov chain and, for each sample, running a particle filter on the validation data using 100 particles. As can be seen, we obtain an accurate prediction of the disease activity, which falls within the estimated 95 % credibility intervals, one month in advance.

Discussion

PGAS is a novel approach to PMCMC that provides the statistician with an off-the-shelf class of Markov kernels which can be used to simulate, for instance, the typically high-dimensional and highly autocorrelated state trajectory in a state-space model. This opens up for using PGAS as a key component in different inference algorithms, enabling both Bayesian and frequentist parameter inference as well as state inference. However, PGAS by no means limited to inference in state-space models. Indeed, we believe that the method can be particularly useful for models with more complex dependencies, such as non-Markovian, nonparametric, and graphical models.

The PGAS Markov kernels are built upon two main ideas. First, by conditioning the underlying SMC sampler on a reference trajectory the correct stationary distribution of the kernel is enforced. Second, ancestor sampling enables movement around the reference trajectory which drastically improves the mixing of the sampler. In particular, we have shown empirically that ancestor sampling makes the mixing of the PGAS kernels robust to a small number of particles as well as to large data records.

Ancestor sampling is basically a way of exploiting backward simulation ideas without needing an explicit backward pass. Compared to PGBS, a conceptually similar method that does require an explicit backward pass, PGAS has several advantages, most notably for inference in non-Markovian models. When using the proposed truncation of the backward weights, we have found PGAS to be more robust to the approximation error than PGBS, yielding up to an order-of-magnitude improvement in accuracy. An interesting topic for future work is to further investigate the effect on these samplers by errors in the backward weights, whether these errors arise from a truncation or some other approximation of the transition density function. It is also worth pointing out that for non-Markovian model PGAS is simpler to implement than PGBS as it requires less bookkeeping. It can also be more memory efficient since it does not require storage of intermediate quantities that are needed for a separate backward simulation pass. See for a related discussion on path storage in the particle filter.

Other directions for future work include further analysis of the ergodicity of PGAS. While the established uniform ergodicity result is encouraging, it does not provide information about how fast the mixing rate improves with the number of particles. Finding informative rates with an explicit dependence on $N$ is an interesting, though challenging, topic for future work. It would also be interesting to further investigate empirically the convergence rate of PGAS for different settings, such as the number of particles, the amount of data, and the dimension of the latent process.

Appendix A Proofs

Let $B\in\mathcal{X}^{T}$ be a measurable rectangle: $B=\times_{t=1}^{T}B_{t}$ with $B_{t}\in\mathcal{X}$ for $t=1,\,\dots,\,T$. Then,

For $t=T-1$, the induction hypothesis holds since $b_{T}$ and $j_{T}$ are equally distributed (both are drawn from the discrete distribution induced by the weights $\{w_{T}^{i}\}_{i=1}^{N}$). Let

where we recall that $w_{t}^{i}=W_{\theta,t}(x_{t}^{i})$ and where the last equality follows from (34). Consider,

Using the Markov property of the generated particle system and the tower property of conditional expectation, we have

Hence, since the function $x_{t+1}\mapsto\Lambda_{t}(x_{t+1},h_{t})$ is bounded, we can use the induction hypothesis to write (49) as

A.2 Proof of Proposition 2

With $M=T-t+1$ and $w(k)=w_{t-1}^{k}$, the distributions of interest are given by

respectively. Let $\varepsilon_{s}\triangleq\max_{k,l}\left(h_{s}(k)/h_{s}(l)-1\right)\leq A\exp(-cs)$ and consider

It follows that the KL divergence is bounded according to,

References