{'fraction_non_alphanumeric': 0.04070679735292111, 'fraction_numerical': 0.014031337411593478, 'mean_word_length': 4.239366138448707, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the use of backward simulation in particle Markov chain Monte Carlo methods

Recently, Andrieu, Doucet and Holenstein (2010) introduced a general framework for using particle filters (PFs) to construct proposal kernels for Markov chain Monte Carlo (MCMC) methods. This framework, termed Particle Markov chain Monte Carlo (PMCMC), was shown to provide powerful methods for joint Bayesian state and parameter inference in nonlinear/non-Gaussian state-space models. However, the mixing of the resulting MCMC kernels can be quite sensitive, both to the number of particles used in the underlying PF and to the number of observations in the data. In the discussion following (Andrieu et al., 2010), Whiteley suggested a modified version of one of the PMCMC samplers, namely the particle Gibbs (PG) sampler, and argued that this should improve its mixing. In this paper we explore the consequences of this modification and show that it leads to a method which is much more robust to a low number of particles as well as a large number of observations. Furthermore, we discuss how the modified PG sampler can be used as a basis for alternatives to all three PMCMC samplers derived in (Andrieu et al., 2010). We evaluate these methods on several challenging inference problems in a simulation study. One of these is the identification of an epidemiological model for predicting influenza epidemics, based on search engine query data.


Introduction
Particle Markov chain Monte Carlo (PMCMC) is an umbrella for a collection of methods introduced in [1]. The fundamental idea underlying these methods, is to use a sequential Monte Carlo (SMC) sampler, i.e. a particle filter (PF), to construct a proposal kernel for an MCMC sampler. The resulting methods were shown to be powerful tools for joint Bayesian parameter and state inference in nonlinear, non-Gaussian state-space models. However, to obtain reasonable mixing of the resulting MCMC kernels, it was reported that a fairly high number of particles was required in the underlying SMC samplers. This might seem like a very natural observation; in order to obtain a good proposal kernel based on a PF, we should intuitively use sufficiently many particles to obtain high accuracy in the PF. However, this is also one of the major criticisms against PMCMC. If we need to run a PF with a high number of particles at each iteration of the PMCMC sampler, then a lot of computational resources will be wasted and the method will be very computationally intense (or "computationally brutal" as Flury and Shephard put it in the discussion following [1]). It is the purpose of this work to show that this need not be the case. We will discuss alternatives to each of the three PMCMC methods derived in [1], that will function properly even when using very few particles. The basic idea underlying these modified PMCMC samplers was originally proposed by Whiteley [2].
To formulate the problem that we are concerned with, consider a general, discrete-time state-space model with state-space X, parameterised by θ ∈ Θ, The initial state has density π θ (x 1 ). We take a Bayesian viewpoint and model the parameter as a random variable with prior density p(θ). Given a sequence of observations y 1:T {y 1 , . . . , y T }, we wish to estimate the parameter θ as well as the system state x 1:T . That is, we seek the posterior density p(θ, x 1:T | y 1:T ).
During the last two decades or so, Monte Carlo methods for state and parameter inference in this type of nonlinear, non-Gaussian state-space models have appeared at an increasing rate and with increasingly better performance. State inference via SMC is thoroughly treated by for instance [4][5][6]. For the problem of parameter inference, the two overview papers [7,8] provide a good introduction to both frequentistic and Bayesian methods. Some recent results on Monte Carlo approaches to maximum likelihood parameter inference in statespace models can be found in [9][10][11][12]. Existing methods based on PMCMC will be discussed in the next section, where we also provide a preview of the material presented in the present work.

An overview of PMCMC methods
In [1], three PMCMC methods were introduced to address the inference problem mentioned above. These methods are referred to as particle Gibbs (PG), particle marginal Metropolis-Hastings (PMMH) and particle independent Metropolis-Hastings (PIMH).
Let us start by considering the PG sampler, which is an MCMC method targeting the joint density p(θ, x 1:T | y 1:T ). In an "idealised" Gibbs sampler (see e.g. [13] for an introduction to Gibbs sampling), we would target this density by the following two-step sweep, • Draw θ | x 1:T ∼ p(θ | x 1:T , y 1:T ).
The first step of this procedure can, for some problems be carried out exactly (basically if we use conjugate priors for the model under study). In the PG sampler, this is assumed to be the case. However, the second step, i.e. to sample from the joint smoothing density p θ (x 1:T | y 1:T ), is in most cases very difficult. In the PG sampler, this is addressed by instead sampling a particle trajectory x 1:T using a PF. More precisely, we run a PF targeting the joint smoothing density. We then sample one of the particles at the final time T , according to their importance weights, and trace the ancestral lineage of this particle to obtain the trajectory x 1:T . However, as we shall see in Section 4.2 this can lead to very poor mixing when there is degeneracy in the PF (see e.g. [5,14] for a discussion on PF degeneracy). This will inevitably be the case when T is large and/or N is small.
One way to remedy this problem is to append a backward simulator to the PG sampler, leading to a method that we denote particle Gibbs with backward simulation (PG-BSi). The idea of using backward simulation in the PG context was mentioned by Whiteley [2] in the discussion following [1]. To make this paper self contained, we will present a derivation of the PG-BSi sampler in Section 4. The reason for why PG-BSi can avoid the poor mixing of the original PG sampler will be discussed in Section 4.3. Furthermore, in Section 4.2 we shall see that the PG-BSi sampler can operate properly even with very few particles and vastly outperform PG in such cases. Now, as mentioned above, to apply the PG and the PG-BSi samplers we need to sample from the conditional parameter density p(θ | x 1:T , y 1:T ). This is not always possible, since conjugate priors can be both undesirable and unavailable (e.g. if there is a nonlinear parameter dependence in the model). For such cases, the PMMH of [1] can be used. Similarly to PG, this method targets the joint posterior density p(θ, x 1:T | y 1:T ), but parameter values are drawn from some proposal kernel q(θ | θ ) in a Metropolis-Hastings (MH) fashion. Hence, the method does not require the density p(θ | x 1:T , y 1:T ) to be available for sampling.
In this work, we propose an alternative method called Metropolis within particle Gibbs (MwPG), derived in Section 5. The MwPG sampler is based on the PG-BSi sampler, but instead of sampling exactly from p(θ | x 1:T , y 1:T ) we use an MH step to sample the parameters. This might at first seem superfluous and appear as a poor imitation of the PMMH sampler. However, as we shall see in Section 5.2, the MwPG sampler can vastly outperform PMMH, even when the two methods use the exact same proposal density to sample the parameters. The reason is that PMMH requires accurate estimates of the likelihood p θ (y 1:T ) to get a high acceptance probability [1]. This implies that also this method requires a high number of particles in the PF. On the contrary, as mentioned above, the PG-BSi sampler has the ability to function properly even for a very low number of particles. Since they share a common base, this desirable property of the PG-BSi sampler will "rub off" on the MwPG sampler.
Finally, the PIMH introduced in [1] is as a method targeting the joint smoothing density p θ (x 1:T | y 1:T ) when the parameter θ is assumed to be known. Hence, it is a competitor to existing particle smoothing algorithms, see e.g. [14][15][16][17]. The PIMH can be seen as a collapsed version of the PMMH, and will suffer from the same problem when we decrease the number of particles. To mitigate this, we propose to instead use a collapsed version of the PG-BSi sampler. The details are given in Section 6 and the proposed method is compared with the PIMH and a state of the art particle smoother [16] in Section 6.3.
Since the publication of the seminal paper [1], several contributions have been made which relate to the present work. We have published a preliminary version of the PG-BSi sampler in [18]. We have also applied this to the problem of Wiener system identification in [19]. In [20], a version of PG with backward simulation is used for inference in switching linear dynamical systems. In this work, it was observed that the backward sampling step could indeed drastically improve the performance of the PG sampler. In [3] a similar approach is used for multiple change-point models.
Backward simulation has also been used in the context of PIMH and PMMH [21]. However, it should be noted that their approach is fundamentally different from the MwPG sampler suggested in this work. In [21], the aim is to "improve the quality" of the sampled state trajectories. The acceptance probabilities in the MH samplers suggested in [21] are identical to those used in the original PMMH and PIMH samplers. Hence, they will suffer from the same problems when using few particles.

Particle filter and backward simulator
Before we go in to the details of the PG-BSi method in Section 4, we review the PF and the backward simulator. Throughout this section we assume that the parameter θ is fixed and discuss how to approximately sample from the joint smoothing density p θ (x 1:T | y 1:T ).

Auxiliary particle filter
The auxiliary particle filter (APF) [22] can be used to approximate (and approximately sample from) the joint smoothing distribution. The reason for why we choose an APF, rather than a "standard" PF, is that it is more general and that the notion of auxiliary variables to determine ancestor indices is well suited for the derivation of the PG-BSi sampler in Section 4.
We shall assume that the reader is familiar with the APF and make the presentation quite brief, mainly just to introduce the required notation. For readers not familiar with the PF or the APF, we refer to [4,5,22]. Let {x m 1:t−1 , w m t−1 } N m=1 be a weighted particle system targeting p θ (x 1:t−1 | y 1:t−1 ). That is, the particle system defines an empirical distribution, which approximates the target distribution. In the APF, we propagate this sample to time t by sampling from a proposal kernel, Here, the variable i t is the index to an "ancestor particle" x it t−1 and R θ t is a proposal kernel which proposes a new particle at time t given this ancestor. The factors ν it t−1 = ν θ t−1 (x it t−1 , y t ), known as adjustment multiplier weights, are used in the APF to increase the probability of sampling ancestors that better can describe the current observation. If these adjustment weights are identically equal to 1, we recover the standard PF.
Once we have generated N new particles (and ancestor indices) from the kernel M θ t , the particles are assigned importance weights according to w m t = where the weight function is given by, This results in a new weighted particle system {x m 1:t , w m t } N m=1 , targeting the joint smoothing density at time t.
The method is initialised by sampling from a proposal density x m 1 ∼ R θ 1 (x 1 ) and weighting the samples according to w m 1 = W θ 1 (x m 1 ) where the weight function is given by, . After a complete run of the APF we can sample a trajectory from the empirical joint smoothing density by choosing particle x m 1:T with probability proportional to w m T . Note that this particle trajectory consists of the ancestral lineage of x m T . Hence, if this lineage is defined by the indices b 1:T (with b T = m) we obtain a sample For later reference, we note that the proposal kernel used in the APF is not allowed to be completely arbitrary. In principle, we need to make sure that the support of the proposal covers the support of the target. Hence, we make the following assumption. Assumption 1. Let S = {θ ∈ Θ : p(θ) > 0}. Then, for any θ ∈ S and any t ∈ {1, . . . , T }, S θ t ⊂ Q θ t where we have defined,

Backward simulation
In the original PG sampler by [1], sample trajectories are generated as in (3). However, due to the degeneracy of the APF, the resulting Gibbs sampler can suffer from poor mixing. We will illustrate this in Section 4.2. As previously mentioned, we will focus on an alternative way to generate an approximate sample trajectory from the joint smoothing density, known as backward simulation. It was introduced as a particle smoother in [14,23]. In [2], it was suggested as a possible way to increase the mixing rate of a PG kernel. Consider a factorisation of the joint smoothing density as, Here, the backward kernel density can be written as, We note that the backward kernel depends on the filtering density p θ (x t | y 1:t ).
The key enabler of the backward simulator is that this density (in many cases) can be readily approximated by an (A)PF, without suffering from degeneracy.
(b) Set j t = m with probability w m t|T .
By using the APF to approximate the filtering density, we obtain an approximation of the backward kernel, Based on this approximation, we may sample a particle trajectory backward in time, approximately distributed according to the joint smoothing density. The procedure is given in Algorithm 1. The algorithm generates a set of indices j 1:T defining a backward trajectory according to Note that, since we only need to generate a single backward trajectory, the computational cost of the backward simulator is O(N ), i.e. of the same order as the APF. Hence, the computational complexity of the PG sampler is more or less unaffected by the introduction of the backward simulation step. Still, it is possible to reduce the complexity of the backward simulator even further by making use of rejection sampling [16].

Particle Gibbs with backward simulation
In this section we will derive the PG-BSi sampler and discuss its properties. This forms the basis for all the methods proposed in this work.

The PG-BSi sampler
Recall the idealised Gibbs sampler, briefly outlined in Section 2. The two steps of the method which are performed at each iteration are, As previously mentioned, the idea behind the PG sampler by [1] is, instead of sampling according to (5b), to run an APF to generate a sample trajectory x 1:T as in (3). In the PG-BSi sampler, we instead generate a sample trajectory by backward simulation, as in (4). However, in either case, to simply replace step (5b) of the idealised Gibbs sampler with such an approximate sampling procedure will not result in a valid approach. In [1] they solve this by replacing the APF with a so called conditional APF (CAPF), and show that the resulting PG sampler is a valid MCMC method. In this section, we make a similar derivation for the PG-BSi sampler.
The key in deriving the PG-BSi sampler is to consider a Markov chain in which the state consists of all the random variables generated by the particle filter and the backward simulator. Hence, the state of the Markov chain will not only be the particle trajectory x j1 1 , . . . , x j T T , but the complete set of particles, particle indices and backward trajectory indices. For this cause, let us start by determining the density of all the random variables generated by the APF. Let . . , i N t } be the particles and ancestor indices, respectively, generated by the APF at time t. Note that i m t is (the index of) the ancestor of particle x m t . For a fixed θ the APF then samples from, Note that the APF is simply a method which generates a single sample from ψ θ , on the extended space . It can be beneficial to adopt this view on the APF in order to understand the inner workings of the PMCMC samplers. By completing this with a backward simulation according to Algorithm 1, we generate the indices j 1 , . . . , j T defining a backward trajectory. The joint density of {x 1:T , i 2:T , j 1:T } is then given by, . (6) Thus, the state trajectory x j 1:T 1:T {x j1 1 , . . . , x j T T } will be distributed (jointly with the random indices j 1:T ) according to a marginal of (6), ψ ,θ (x j 1:T 1:T , j 1:T ), which in general is not equal to the joint smoothing density. In other words, the extracted trajectory x j 1:T 1:T is not distributed according to p θ (x 1:T | y 1:T ) and to simply use this trajectory in step (5b) of the idealised Gibbs sampler would not be a valid approach.
The key to circumventing this problem is to introduce a new, artificial target density, with two properties: 1. The artificial target density should admit the joint smoothing density as one of its marginals.
2. The artificial target density should be "similar in nature" to ψ θ , to enable the application of an APF in the sampling procedure.
Here, we define this artificial target density according to, Algorithm 2 CAPF -conditional auxiliary particle filter (conditioned on 2. For t = 2, . . . , T do: This construction is due to [21], who apply backward simulation in the context of PMMH and PIMH. The artificial target density (7) differs from the one used in [1] in the last factor. One of the important properties of the density φ is formalised in the following proposition. This proposition has the following important implication; if we can construct an MCMC method with stationary distribution φ, then the sub-chain given by the variables {θ, x j 1:T 1:T } will have stationary distribution p(θ, x 1:T | y 1:T ). The construction of such an MCMC sampler is enabled by the second important property of our artificial density, namely that it is constructed in such a way that we are able to sample from the conditional φ(x −j 1:T 1:T , i 2:T | θ, x j 1:T 1:T , j 1:T ). Here we have introduced the notation x −j 1:T This can be done by running a so called conditional APF (CAPF), introduced in [1] 1 . This method is given in Algorithm 2 and its relationship to φ is formalised in Proposition 2. .
Consequently, the CAPF given in Algorithm 2, conditioned on {θ, x j 1:T 1:T , j 1:T }, generates a sample from this conditional distribution.
Proof. See Appendix A.
Note that the conditional density considered in Proposition 2 only depends on known quantities, as opposed to the density φ defined in (7) which depends explicitly on the joint posterior p(θ, x 1:T | y 1:T ) (which in general is not known).
As a final component of the PG-BSi sampler we need a way to generate a sample state trajectory. As previously pointed out, we aim to do this by running a backward simulator. The fact that the backward simulator indeed generates a sample from one of the conditionals of φ is assessed in the following proposition. .
Consequently, the backward simulator given in Algorithm 1, conditioned on {θ, x 1:T , i 2:T }, generates a sample from this conditional distribution.
Proof. See Appendix A.
We are now ready to present the PG-BSi sampler, given in Algorithm 3. Based on the propositions and discussion above, the three steps of the method [2(a)-2(c)] can be interpreted as: Hence, each step of Algorithm 3 is a standard Gibbs update, targeting the distribution φ, and will thus leave φ invariant 2 . It follows that φ is a stationary distribution of the PG-BSi sampler.
It is worth to note that the sampling scheme outlined above is not a "complete" multistage Gibbs sweep, since it does not loop over all the variables of the model. More precisely, we do not sample new locations for the particles x j 1:T 1:T in any of the steps above (the same is true also for the original PG sampler). This might raise some doubts about ergodicity of the resulting Markov chain. However, as established in Theorem 1 below, the PG-BSi sampler is indeed ergodic (c) Run a backward simulator to generate j 1:T (r). Set x 1:T (r) to the corresponding particle trajectory.
(as long as the idealised Gibbs sampler is), meaning that it is a valid MCMC method. Intuitively, the reason for this is that the collection of variables that is "left out" is chosen randomly at each iteration. Proof. See Appendix A.
It should be noted that [2,20] motivate the use of backward simulation in the PG sampler in a slightly different way. They simply note that the backward simulator is a Markov kernel with invariant distribution given by the original artificial target proposed in [1]. Hence, they view the backward simulator as a way to resample the ancestor indices for the extracted trajectory x j 1:T 1:T .

Numerical illustration
We have argued that the PG-BSi sampler in Algorithm 3 should be preferable over the standard PG sampler in the sense that it produces a Gibbs sampler with better mixing properties. Since the PG-BSi sampler is the basis for all the methods proposed in this paper, we pause for a while and consider a numerical illustration of this improved mixing. Hence, consider the following nonlinear state-space model,  We generate a set of observations y 1:500 according to (8) with σ 2 v = 10 and σ 2 e = 1. The parameter priors are modelled as inverse gamma distributed with hyper-parameters a = b = 0.01. We then employ the PG sampler of [1] and the PG-BSi sampler of Algorithm 3 to estimate the posterior parameter distribution. Since the theoretical validity of the PG-BSi and the PG samplers can be assessed for any number of particles N ≥ 2 (see Theorem 1), it is interesting to see the practical implications of using very few particles. Hence, we run the methods four times on the same data with N = 5, 20, 1000 and 5000 (N = 5000 as was used in [1]). The parameters are initialised at θ(0) = 10 10 T in all experiments.
In Figure 1 we show the estimated posterior means of the standard deviations σ v and σ e (i.e. the square roots of the means of θ(1 : r)) vs. the number of MCMC iterations r. The left column shows the results for the PG sampler. When we use N = 5000 particles, there is a rapid convergence towards the true values, indicating that the method mixes well and quickly finds a region of high posterior probability. However, as we decrease the number of particles, the convergence is much slower. Even for N = 1000, the method struggles and for N = 20 and N = 5 it does not seem to converge at all (in a reasonable amount of time). The key observation is that the PG-BSi sampler (right column), on the other hand, seems to be more or less unaffected by the large decrease in the number of particles. In fact, the PG-BSi using just N = 5 performs equally well as the PG sampler using N = 1000 particles.

Why the big difference?
To see why the PG-BSi sampler is so much more insensitive to a decreased number of particles, consider a toy example where we have generated T = 50 observations from a 1st order linear system. Figure 2 (left) shows the particle tree generated by the CAPF at iteration r of a PG sampler, using N = 30 particles. Due to the degeneracy of the PF, there is only one distinct particle trajectory for t ≤ 32. In the PG sampler, we draw x 1:T (r) by sampling among the particles at time T and tracing the ancestral path of this particle. This trajectory is illustrated as a thick black line in the figure. At the next iteration of the PG sampler we run a CAPF, conditioned on this particle trajectory. This results in the tree shown in Figure 2 (right). Due to degeneracy, we once again obtain only a single distinct particle trajectory for t ≤ 29. Now, since we condition on the trajectory x 1:T (r), the particle tree generated at iteration r + 1 must contain x 1:T (r). Hence, all particle trajectories available at iteration r + 1 are identical to x 1:29 (r) up to time 29. Consequently, when we sample x 1:T (r + 1) at iteration r + 1, this trajectory will to a large extent be identical to x 1:T (r). This results in a poor exploration of the state space, which in turn means that the Gibbs kernel will mix slowly.
Based on this argument we also note that the mixing will be particularly poor when the degeneracy of the CAPF is severe. This will be the case if the length of the data record T is large and/or if the number of particles N is small. This is consistent with the results reported in the previous section.
The reason for why PG-BSi can circumvent the poor mixing of the PG kernel is that the backward simulator explores all possible particle combinations when generating a particle trajectory, and is thus not constrained to one of the ancestral paths. Consider the particles generated by the CAPF at iteration r of a PG-BSi sampler, shown in Figure 3 (left). A backward trajectory x 1:T (r), shown as a thick black line, is generated by a backward simulator. At iteration r + 1 of the PG-BSi sampler we will run a CAPF conditioned on x 1:T (r), generating the particles shown in Figure 3 (right). When we again apply a backward simulator   The "true" values are shown as thick solid lines. Note that the estimate of σ v does not converge to this "true" value, but this is not surprising since we consider just one data realisation. Hence, the posterior mean may very well differ from the "true" value.
to this collection of particles we will, with high probability, sample a trajectory x 1:T (r + 1) which is entirely different from x 1:T (r). This leads to a much better exploration of the state space and thus a faster mixing Gibbs kernel.

The MwPG sampler
Quite often, it is not possible to sample exactly from the conditional parameter density p(θ | x 1:T , y 1:T ). If this is the case, we can replace step 2(a) in Algorithm 3 with an MH update. This results in what we refer to as a Metropolis within particle Gibbs (MwPG) sampler. The name stems from the commonly used term Metropolis-within-Gibbs, which refers to a hybrid MH/Gibbs sampler (see [24,25]). Hence, let us assume that x 1:T is a fixed trajectory and consider Algorithm 4 MwPG -Metropolis within particle Gibbs 1. Initialise: Set θ(0), x 1:T (0) and j 1:T (0) arbitrarily. 2. For r ≥ 1, iterate: (a) MH step for sampling a parameter: • Sample θ ∼ q(θ | θ(r − 1), x 1:T (r − 1)).
(c) Run a backward simulator to generate j 1:T (r). Set x 1:T (r) to the corresponding particle trajectory.
the problem of sampling from p(θ | x 1:T , y 1:T ) using MH. We thus choose some proposal density q(θ | θ , x 1:T ), which may depend on the previous parameter θ and the state trajectory x 1:T . The proposal can also depend on the fixed measurement sequence y 1:T , but we shall not make that dependence explicit. We then generate a sample θ ∼ q(θ | θ , x 1:T ) and accept this with probability which is a standard MH update (see e.g. [13]). The acceptance probability can be computed directly from the quantities defining the model (1) since, By plugging this MH step into the PG-BSi sampler we obtain the MwPG method, presented in Algorithm 4. Due to the fact that the MH step will leave the target density invariant, the stationary distribution of the MwPG sampler is the same as for the PG-BSi sampler. MwPG can be seen as an alternative to the PMMH sampler by [1]. Note that, as opposed to case of PMMH, the acceptance probability (9) does not depend explicitly on the likelihood p(y 1:T | θ).

Numerical illustration -estimating variances
Let us return to the simulation example studied in Section 4.2. We use the same batch of data y 1:500 , generated from model (8) with σ 2 v = 10 and σ 2 e = 1. As before, we wish to estimate the noise variances and set θ = {σ 2 v , σ 2 e }, but now we apply the PMMH sampler of [1] and the MwPG sampler in Algorithm 4. As proposal kernel for the parameter, we use a Gaussian random walk with standard deviation 0.15 for σ 2 v and 0.08 for σ 2 e (both methods use the same proposal).
Running means for the estimated standard deviations, for both methods, are shown in Figure 4. The difference between MwPG and PMMH is similar to that between PG-BSi and PG, especially when we consider the effect of using few particles. Compared to the PG-BSi sampler operating on the same data (see Figure 1 and note the different scaling of the axes) the MwPG sampler is slower to converge. However, this is expected, since the latter does not make use of the conjugacy of the priors.
In Figure 5 we show scatter plots for the Markov chains after 50000 iterations, using a burnin of 10000 iterations. From these plots, it is clear that the estimated posterior parameter distribution for the MwPG sampler is rather insensitive to the number of particles used, and similar to what we get from PMMH using N = 5000. The same effect can be seen by analysing the autocorrelation functions, shown in Figure 6. Finally, the average acceptance probabilities for the two methods and for different number of particles are given in Table 1.
It is worth to note that, in the limit N → ∞, the PMMH will "converge to" an idealised, marginal MH sampler. The MwPG, on the other hand, "converges to" an idealised Metropolis within Gibbs sampler. Hence, for a very large number of particles, we expect that PMMH should outperform MwPG. Furthermore, and more importantly, whether or not MwPG is preferable over PMMH should be very much problem dependent. If there is a strong dependence between the states and the parameters, we expect that the MwPG sampler should suffer from poor mixing, since it samples the parameters conditioned on the states, and vice versa. In such cases, the PMMH might perform better, since it more closely resembles a marginal MH sampler.

Numerical illustration -estimating all parameters
So far, we have focused on estimating the noise variances in the model (8). However, Johannes, Polson and Yae [26] suggested to identify all the parameters θ = {β 1 , β 2 , β 3 , α, σ 2 v , σ 2 e } in model (8), as a more challenging test of PMCMC. They use the same parameter values as in Section 4.2, but with σ 2 v = 1 and  σ 2 e = 10, and apply an alternative, SMC based method with 300000 particles to estimate the parameters. Here, we address the same problem with MwPG using N = 5 particles and 60000 MCMC iterations.
We use vague normal distributed priors for β 1 , β 2 and β 3 , a uniform prior over [1,3] for α and inverse gamma priors with a = b = 0.01 for σ 2 v and σ 2 e . To put even more pressure on the method, Andrieu, Doucet and Holenstein [1, p. 338] suggested to use a larger data set. Therefore, we use T = 2000 samples, generated from the model (8) and apply the MwPG sampler. To propose values for α, we use a Gaussian random walk with standard deviation 0.02, constrained to the interval [1,3]. For the remaining parameters, we make use of conjugacy of the priors and sample exactly from the conditional posterior densities. Figure 7 show the estimated, marginal posterior parameter densities for the six parameters, computed after a burnin of 10000 iterations. The method appears to do a good job at finding the posterior density, even in this challenging scenario and using only 5 particles.

.1 Collapsed PG-BSi
Assume now that the parameter θ is known and that our objective is to solve the smoothing problem for the model (1). That is, we seek the joint smoothing density p θ (x 1:T | y 1:T ) or some marginal thereof. Several particle based methods have been presented in the literature, devoted to this problem, see e.g. [14][15][16][17]. Many of these have quadratic complexity in the number of particles. Though, quite recently some advancements have been made towards linear complexity, see [16,17]. In [1], the PIMH sampler is suggested as a PMCMC method targeting the joint smoothing density. The benefit of this is that the stationary distribution of the PIMH is the exact smoothing distribution, removing any bias that normally is present in more "traditional" particle smoothers. However, if we run the PIMH for R iterations using N particles, the computational complexity is O(N R). In order to obtain a reasonable acceptance probability N needs to be fairly high, rendering the method computationally demanding. It is desirable to find a similar method that works even with a low number of particles.
The alternative that we propose in this paper is very straightforward. We simply run a PG-BSi according to Algorithm 3, but since θ is fixed we skip step 2(a). This collapsed PG-BSi sampler will, similarly to PIMH, generate a Markov chain whose stationary distribution is the exact smoothing distribution. Though, a major difference between the two methods is that the collapsed PG-BSi sampler will accept all the generated particle trajectories, whereas the PIMH only will accept a portion of them in an MH fashion. The acceptance probability of the PIMH will typically deteriorate as we decrease the number of particles N or increase the length of the data record T (see [1, Figure 3]).
The complexity of the collapsed PG-BSi sampler is still O(N R), but since N is allowed to be fairly small in this case, we believe that this method could be a serious competitor to existing particle smoothers. Note that the method is exactly linear in the number of generated trajectories R. Furthermore, the PMCMC based smoothers are trivially parallelisable, simply by running multiple parallel chains.

Multiple trajectories
One possible extension, to further improve the performance of the collapsed PG-BSi, is to use multiple trajectories as suggested for the PIMH sampler in [21]. Assume that we wish to estimate E[h(x 1:T ) | y 1:T ] for some function h and that we run R iterations of the PG-BSi sampler (possibly with some burnin). The most straightforward estimator of h is then h(x 1:T (r)).
An alternative is to use the Rao-Blackwellised estimator To compute the expectation in this expression, we would have sum over all possible backward trajectories. Since there are N T possible trajectories, this is in general not feasible, unless h is constrained to be a function of only some small subset of the x:s, e.g. a single x t or a pair {x t , x t+1 }. However, a compromise would be to generate multiple backward trajectories in the backward simulator, providing an estimate of the above expectation. Hence, at iteration r of the PG-BSi sampler, instead of generating just a single backward trajectory as in Algorithm 3, we generate M trajectories,x m 1:T (r) for m = 1, . . . , M and use the estimatorĥ h(x m 1:T (r)).
See [21] for a more in-depth discussion on using multiple trajectories in the context of PIMH. Since the backward simulator generates conditionally independent and identically distributed (i.i.d.) samples (given {x 1:T , i 2:T }), we can set the state of the Markov chain to any of the generated backward trajectories, e.g. x 1:T (r) =x 1 1:T (r).

Numerical illustration
To illustrate the applicability of the collapsed PG-BSi for joint smoothing, we again consider the same batch of T = 500 samples from the model (8) as used in Section 4.2 and Section 5.2. We apply the PIMH, the collapsed PG-BSi and a state of the art particle smoother, known as fast forward filter/backward simulator (fast FFBSi) derived in [16]. The PMCMC samplers are run using N = 20 particles for M = 1000 iterations. The fast FFBSi is run using N = 1000 particles and backward trajectories. To obtain something that we can visualise, let us consider the marginal smoothing density at time t = 151, i.e. p θ (x 151 | y 1:T ). Figure 8 shows the corresponding cumulative distribution function 3 (CDF), as well as the empirical CDFs from the three methods.
One question of interest is if it is possible to avoid burnin, by using a good initialisation of the MCMC sampler. In the example above we did not use any burnin and the chain was initialised by a run of an (unconditional) PF using N = 20 particles, followed by a backward simulation. Another heuristic approach is to estimate the joint smoothing distribution by running a "standard" particle smoother, e.g. the fast FFBSi, with a moderate number of particles. We can then initiate the collapsed PG-BSi sampler by sampling from the estimated joint smoothing distribution, which is hopefully close to the stationary distribution. This heuristic can be of particular interest it we are able to run several parallel chains. We can then apply a fast FFBSi to generate M backward trajectories, which are conditionally i.i.d. samples from the empirical smoothing distribution, to initiate M independent chains.

Epidemiological model
As a final numerical illustration, we consider the identification of an epidemiological model using the MwPG sampler, based on search engine query observations. Seasonal influenza epidemics each year cause millions of severe illnesses and hundreds of thousands of deaths world-wide [27]. 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) epidemiological model with environmental noise and seasonal fluctuations [28,29], which is discretised according to the Euler-Maruyama method, Here, S t , I t and R t represent the number of susceptible, infected and recovered individuals at time t (months), respectively. The sampling time dt is chosen as 1/30 in our experiment, i.e. we sample the model once a day. The state of the system is comprised of x t = S t I t T . The total population size P = 10 6 and the host birth/death rate µ = 0.0012 are both assumed known. The seasonally varying transmission rate is given by β(t) =β(1 + α sin(2πt/12)) and the rate of recovery is given by γ. The process noise is assumed to be Gaussian according . In [29], the PMMH sampler is used to identify a similar SIR model, though with a different observation model than what we consider here (see below). A different Monte Carlo strategy, based on a particle filter with an augmented state space, for identification of an SIR model is proposed in [30]. Here, we use an observation model which is inspired by the Google Flu Trends project [27]. The idea is to use the frequency of influenza related search engine queries to infer knowledge about the dynamics of the epidemic. As proposed by [27], we use a linear relationship between the observations and the log-odds of infected individuals, i.e. We place a normal-inverse-gamma prior on the pair {ρ, τ 2 }, i.e. p(ρ, τ 2 ) = N (ρ; µ ρ , c ρ τ 2 )IG(τ 2 ; a τ , b τ ). The hyperparameters are chosen as µ ρ = 1, c ρ = 0.5 and a τ = b τ = 0.01. We place an inverse gamma prior on F 2 with hyperparameters a F = b F = 0.01 and a normal prior on γ with µ γ = 3 and σ 2 γ = 1. The above mentioned priors are conjugate to the model defined by (10), which is exploited in the MwPG sampler by drawing values from their conditional posteriors. For the parametersβ and α we use flat improper priors over the positive real line and sample new values according to MH moves, as described in Section 5. Samples are proposed independently forβ and α, according to Gaussian random walks with variances 5 · 10 −3 and 5 · 10 −4 , respectively.
We generate 8 years of data with daily observations. The number of infected individuals I t over this time period is shown in Figure 10. The first half of the data batch is used for estimation of the model parameters. For this cause, we employ the MwPG sampler using N = 20 particles for 100000 MCMC iterations. Histograms representing the estimated posterior parameter distributions are shown in Figure 9. We notice some raggedness in some of the estimated marginals. This implies that they have not fully converged, suggesting a fairly slow mixing of the Markov chain. Still, the estimated distributions seem to provide accurate information about the general shapes and locations of the posterior parameters distributions, and the true parameters fall well within the credible regions of the posteriors.
Finally, the estimated model is used to make 1 month ahead predictions of the disease activity for the subsequent 4 years, as shown in Figure 10. The predictions are made by Monte Carlo sampling from the posterior parameter distribution, followed by a run of 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.

Conclusions and future work
In the original PMCMC paper, Andrieu, Doucet and Holenstein writes "PM-CMC sampling is much more robust and less likely to suffer from this depletion problem [referring to the PF]. This stems from the fact that PMCMC methods do not require SMC algorithms to provide a reliable approximation of p θ (x 1:T | y 1:T ), but only to return a single sample approximately distributed according to p θ (x 1:T | y 1:T )". This observation is the very essence of the PMCMC methods and the reason for why they have shown to be such powerful tools for joint Bayesian state and parameter inference. In this paper we have argued that the use of backward simulation in PMCMC make this claim even stronger. Three methods have been presented, each similar to one of the methods derived in [1]. These are; PG-BSi (similar to PG), MwPG (similar to PMMH) and collapsed PG-BSi (similar to PIMH).
We have shown experimentally that all of these methods work properly, even when we use very few particles in the underlying PFs. This robustness to a decreased number of particles stems from a simple modification of the PG sampler, namely to add a backward simulation sweep to reduce the correlation between the particle trajectories sampled at any two consecutive iterations. This modification does not only enable the application of PMCMC methods with very few particles, but will also allow us to treat larger data sets.
It would be interesting to analyse how one should choose the number of particles. The mixing of the MCMC kernels will increase as we increase N , but the improvement seems to saturate rather quickly (see e.g. Figure 6). Based on this, we believe that it in general is better to use fewer particles and more MCMC iterations. However, deducing to what extent this statement is problem dependent is a topic for future work.
Another interesting avenue for future work is regarding theoretical aspects of the mixing of the PG-BSi sampler. That is, whether or not it is possible to show that PG-BSi always will mix faster than PG for the same number of particles. We believe that this is the case, and some insight into the matter can be gained by considering the maximal correlations for the two methods (see e.g. [31] for a general discussion).
When it comes to a comparison between the MwPG and the PMMH samplers, it is not as clear to deduce which method that is preferable over the other. In fact, we believe that this is very much problem dependent, and that the two methods are complementary. If there is a strong dependence between the states and the parameters, the MwPG sampler should suffer from poor mixing, since it samples these groups of variables conditioned on each other. Basically, this is the case if there is little noise in the system, or more severely, if the model is degenerate (e.g. if noise only enters on some of the states). In the latter case, the method might not even converge, since it is not clear that Assumption 2 is fulfilled. In such cases, the PMMH sampler might be a better choice. On the other hand, from the numerical examples that we have considered, it is clear that the MwPG sampler can substantially outperform the PMMH sampler in other cases. It would be of great interest to conduct a more in-depth comparison between these two methods, to provide some general guidelines for when to use which method.