DP2016-36

A Multiple-Try Extension of the Particle
Marginal Metropolis -Hastings (PMMH)
Algorithm with a n Independent
Proposa l
Takashi KAMIHIGASHI
Hiroyuki WATANABE
December 7, 2016

A Multiple-Try Extension of the Particle
Marginal Metropolis-Hastings (PMMH)
Algorithm with an Independent Proposal∗
Takashi Kamihigashi†

Hiroyuki Watanabe‡

December 6, 2016
Abstract
In this paper we propose a multiple-try extension of the PMMH
algorithm with an independent proposal. In our algorithm, I ∈ N
parameter particles are sampled from the independent proposal. For
each of them, a particle filter with K ∈ N state particles is run. We
show that the algorithm has the following properties: (i) the distribution of the Markov chain generated by the algorithm converges to theposterior of interest in total variation; (ii) as I increases to ∞, the acceptance probability at each iteration converges to 1 with probability
1; and (iii) as I increases to ∞, the autocorrelation of any order for
any parameter with bounded support converges to 0. These results
indicate that the algorithm generates almost i.i.d. samples from theposterior for sufficiently large I. Our numerical experiments suggestthat one can visibly improve mixing by increasing I from 1 to a small
number. This does not significantly increase computation time if acomputer with at least the same number of threads is used.
Keywords: multiple-try method; particle marginal Metropolis-Hastings;
Markov chain Monte Carlo; mixing; state space models∗

Financial support from the Japan Society for the Promotion of Science (JSPS KAKENHI Grant Number 15H05729) is gratefully acknowledged.
†
Research Institute for Economics and Business Administration (RIEB), Kobe University, Rokkodai, Nada, Kobe 657-8501 JAPAN. Email: tkamihig@rieb.kobe-u.ac.jp.
‡
Research Institute for Economics and Business Administration (RIEB), Kobe University, Rokkodai, Nada, Kobe 657-8501 JAPAN. Email: hwatanabe@rieb.kobe-u.ac.jp.

1

Introduction

Over the past few decades, Markov chain Monte Carlo (MCMC) algorithms
have become an indispensable tool for Bayesian inference. A standard MCMC
algorithm such as Metropolis-Hastings (MH) generates a Markov chain thatconverges in distribution to the posterior of interest as the number of iterations increases to ∞ (e.g., Tierney, 1994; Roberts and Rosenthal, 2004).
In this case the posterior can be approximated by the empirical distribution
provided that the number of iterations is large enough. In theory this holds
quite generally, but in practice the Markov chain generated by an MCMC
algorithm often suffers from poor mixing and fails to converge within a reasonable time frame. The problem is particularly severe for state space models,in which a large number of state (or hidden) variables are often treated as
unknown parameters.
As far as state space models are concerned, the problem can be mitigatedby using the particle marginal Metropolis-Hastings (PMMH) algorithm (Andrieu et al., 2009, 2010), which embeds a particle filter, or a sequential Monte
Carlo method (e.g., Doucet et al., 2000), within an MH algorithm. PMMH
belongs to the class of MCMC algorithms known as particle MCMC algorithms (Andrieu et al., 2010), which build on earlier work by Beaumont(2003) and Andrieu and Roberts (2009).
In the PMMH algorithm, given a sample of the static parameter vector from a proposal, the marginal likelihood of observed data is sequentiallyapproximated by Monte Carlo integration using (what we call) “state particles.” The state of the chain is updated with a probability based on theratio of marginal likelihood estimates. While PMMH has been applied andextended in a wide variety of contexts (e.g., Yuan et al., 2012; Kokkala andSärkkä, 2015; Golightly et al., 2015; Fearnhead and Meligkotsidou, 2015), it
leaves room for improvement in terms of static parameters, which are only
updated using a standard MH algorithm.
Earlier efforts to improve the mixing of parameters for static models include the multiple-try Metropolis (MTM) method (Liu et al., 2000). In thismethod, multiple “parameter particles” are drawn from a proposal. Oneof them is resampled according to importance weights, and is accepted to
update the state of the chain with a probability computed using additionalreference samples. This method is closely related to PMCMC algorithms—
in particular, the particle MH (PMH) algorithm—as discussed by Martino
et al. (2016). Pandolfi et al. (2010) provide a generalization of MTM, and1

Bédard et al. (2012) present a scaling analysis of MTM methods. Martinoand Louzada (2016) discuss mixing properties of various MTM algorithms.
While MTM and its extensions have been applied to state space models (e.g., So, 2006; Wang, 2009), studies in the literature that discuss both
PMMH and MTM are surprisingly limited. Exceptions include Andrieu etal. (2009, p. 57), who differentiate PMMH from MTM, and Storvik (2011),
who considers acceptance probabilities for a class of MCMC algorithms including PMMH and MTM. Related studies that do not explicitly mention
PMMH include Martino et al. (2016), who examine the relations between
PMH and MTM schemes, and Tran et al. (2016), who discuss extensions of
their adaptive sampler to MTM and a class of MCMC algorithms including
PMCMC.
In this paper we propose a new algorithm for state space models that combines PMMH and MTM under the assumption that the proposal for staticparameters is independent. We call this algorithm MTiPMMH (multiple-try
independent PMMH) for the reasons specified below. In the MTiPMMH
algorithm, a prespecified number of parameter particles are sampled fromthe independent proposal. For each of them, a particle filter is run with aprespecified number of state particles. Then one of the parameter particles isresampled according to their marginal likelihood estimates. The state of thechain is updated with a probability based on the ratio between the averages
of marginal likelihood estimates.
We show that the MTiPMMH algorithm has the following three asymptotic properties: (i) the distribution of the Markov chain generated by thealgorithm converges to the posterior of interest in total variation; (ii) as thenumber of parameter particles increases to ∞, the acceptance probability ateach iteration converges to 1 with probability 1; and (iii) as the number ofparameter particles increases to ∞, the autocorrelation of any order for any
parameter with bounded support converges to 0. These results are valid forany number of state particles, indicating that MTiPMMH generates almost
i.i.d. samples from the posterior with a sufficiently large number of parameterparticles.
To show these results, we start with a general static model and consider
the MH algorithm with an independent proposal, which we call iMH (independent MH). We extend iMH by introducing parameter particles. Martino
et al. (2016) argue that this extended algorithm (which they call I-MTM2) isa variant of the MTM method with an independent proposal. Based on theirargument, we regard the extended algorithm as a multiple-try (MT) exten2

sion of iMH; thus we call it MTiMH (multiple-try independent MH). While
Martino et al. (2016) only use MTiMH to compare MTM and PMH schemes
without establishing its theoretical properties, we show that MTiMH has thethree asymptotic properties listed above.
With these results for MTiMH in hand, we turn to a state space modeland consider the PMMH algorithm with an independent proposal, which wecall iPMMH (independent PMMH). Following Andrieu et al. (2010) in spirit,
we reformulate iPMMH as an iMH algorithm. This allows us to reformulate
MTiPMMH as an MTiMH algorithm and to show that the three asymptotic
results established for MTiMH hold also for MTiPMMH.
There are at least three advantages to using an independent proposal
rather than a more general Markovian proposal. First, a combination of
an independent proposal and high acceptance probabilities gurantees good
mixing. Note that high acceptance probabilities alone do not imply good
mixing. A case in point is Gibbs sampling, which always has an acceptanceprobability of 1, but often suffers from poor mixing. Second, an independent proposal is easy to construct. For example, one can approximate theempirical distribution of samples from a pilot run by a normal distribution,and use it as an independent proposal. We use this approach in some of
our numerical experiments. Third, an independent proposal makes it easy toutilize parallel computation. We present a parallelized version of MTiPMMH
for which computational cost does not significantly increase as the numberof parameter particles increases, provided that a computer with at least thesame number of threads is used.
The rest of the paper is organized as follows. In Section 2 we introduce
basic notation for state space models and related algorithms, and briefly
discuss the PMMH algorithm. In Section 3 we turn to static models andpresent the MTiMH algorithm. We establish our core results in this static
framework. In Section 4 we present the iPMMH and MTiPMMH algorithms.
We reformulate iPMMH as an iMH algorithm, and apply the results shownfor MTiMH to MTiPMMH. In Section 5 we conduct numerical experiments
by utilizing the simple stochastic volatility model employed by Fearnhead andMeligkotsidou (2015). In Section 6 we conclude the paper with discussion of
our results.

3

2Preliminaries

Given any (possibly multidimensional) random variables a and b, we letp(a, b) denote the joint density of (a, b), p(a) the marginal density of a, andp(a|b) the conditional density of a given b. These densities can be prior or posterior densities. A proposal density is denoted by q instead of p. Throughout
this paper all functions are assumed to Borel-measurable.

2.1

State Space Models

Let nθ , nx , ny ∈ N. Let Θ and X be Borel subsets of Rnθ and Rnx , respectively. Consider a state space model consisting of a state (or hidden) processny∞{xt }∞
t=1 ⊂ X , an observation process {yt }t=1 ⊂ R , and a vector θ ∈ Θ ofunknown static parameters, which are taken as random variables. The conditional distribution of xt+1 given x1:t , y1:t , and θ is given by p(xt+1 |xt , θ), wherex1:t = (x1 , x2 , . . . , xt ), etc. The conditional distribution of yt given x1:t , y1:t−1 ,and θ is given by p(yt |xt , θ). The parameter vector θ is distributed according
to prior p(θ). The initial value x1 is distributed as p(x1 |θ), which can be partof the model or a prior.
Given observed data y1:T with T ∈ N, Bayesian inference is typically
performed based on the posterior density p(θ, x1:T |y1:T ), which satisfiesp(θ, x1:T |y1:T ) = p(θ, x1:T , y1:T )/p(y1:T ) ∝ p(θ, x1:T , y1:T )
"T −1
# T
Y
Y
= p(θ)p(x1 |θ)p(xt+1 |xt , θ)p(yt |xt , θ).
t=1

2.2

(1)
(2)

t=1

The PMMH Algorithm

We start with some notation to present various algorithms in a concise manner. For m, n ∈ N with m ≤ n, we let[m : n] = {m, . . . , n},“m : n” = “m, m + 1, . . . , n,”“n : m” = “n, n − 1, . . . , m.”

(3)
(4)
(5)

a ∧ b = min{a, b}.

(6)

For a, b ∈ R+ , we let

4In what follows, the division operator / precedes the minimum operator ∧(e.g., see line 18 of Algorithm 1). Given w1 , . . . , wN ≥ 0 with w1 +· · ·+wN >0, let D(w1:N ) denote the discrete distribution on [1 : N ] with the probabilityof each n ∈ [1 : N ] given by
wn
PN

m=1 w

m

.

(7)

Algorithm 1 shows pseudocode for the PMMH (particle marginal MH)
algorithm. In line 2, θ∗ is initially set arbitrarily, and is subsequently drawn
from conditional proposal q(θ|θj−1 ). In line 4, xk1 is drawn from q(x1 |y1 , θ∗ );ak
in line 9, xkt+1 is drawn from q(xt+1 |yt+1 , xt t , θ∗ ). In line 11, kT∗ is drawn from
D(wT1:K ). In line 19, “w.p.” means “with probability.”
Andrieu et al. (2010, Theorem 4) show that Algorithm 1 generates asequence {(θj , ξj )}∞j=1 such that the distribution of (θj , ξj ) converges to theposterior p(θ, x1:T |y1:T ) in total variation as j ↑ ∞. This result relies on thefact that given θ, the product in line 14 is an unbiased estimator of p(y1:T |θ):
Z
1:K
1:K
u(θ, x1:K(8)
1:T , y1:T )p(x1:T |θ, y1:T )dx1:T = p(y1:T |θ),whereu(θ∗ , x1:K
1:T , y1:T ) =

T
Y

K −1

K
X

!
wtk

,

(9)

t=1

k=1
1:K
k k∈[1:K]
x1:T = {xt }t∈[1:T ] .

(10)

See Pit et al. (2012, Theorem 1) and Del Moral (2004, Proposition 7.4.1) for
proofs of (8).

3A Multiple-Try Extension of iMH

In this section we consider general static models and study an extension of
the MH algorithm with an independent proposal. The results shown here areapplied to state space models in Section 4.
The basis of our analysis is the MH algorithm with an independent proposal q(θ), which we call iMH. Algorithm 2 shows this algorithm as a functionof (q, g), where q and g are functions satisfying the following assumption,
which is maintained throughout this section.
5

Algorithm 1: PMMH (Particle Marginal Metropolis-Hastings)
1 for j = 1 : J(∈ Θ arbitrary if j = 12θ∗∼ q(θ|θj−1 )if j ≥ 2
3
for k = 1 : K
4
xk1 ∼ q(x1 |y1 , θ∗ )5
w1k = p(y1 |xk1 , θ∗ )p(xk1 |θ∗ )/q(xk1 |y1 , θ∗ )6
7
8

for t = 1 : T − 1
for k = 1 : K
akt ∼ D(wt1:K )ak

9

xkt+1 ∼ q(xt+1 |yt+1 , xt t , θ∗ )10k
= p(yt+1 |xkt+1 , θ∗ )p(xt+1 |xt t , θ∗ )/q(xt+1 |yt+1 , xt t , θ∗ )
wt+1

11
12
13
14
15
16
17
18
19

ak

k∗

kT∗ ∼ D(wT1:K ), x∗T = xTTfor t = (T − 1) : 1k∗k∗kt∗ = at t+1 , x∗t = xt t

Q 
P
k
w
h∗ = p(θ∗ ) Tt=1 K −1 K
k=1 t
if j = 1 then(θ1 , ξ1 , h1 ) = (θ∗ , x∗1:T , h∗ )else∗∗
r = 1 ∧ [h∗ /q(θ( |θj−1 )]/[hj−1 /q(θj−1 |θ )](θj , ξj , hi ) =(θ∗ , x∗1:T , h∗ )w.p. r(θj−1 , ξj−1 , hj−1 ) w.p. 1 − routput: θ1:J , ξ1:J

6

ak

Algorithm 2: iMH(q, g) (Independent Metropolis-Hastings)
1 for j = 1 : J
2θ∗ ∼ q(θ), w∗ = g(θ∗ )/q(θ∗ )3
if j = 1 then4(θ1 , w1 ) = (θ∗ , w∗ )
5
6
7

else
r = 1 ∧ w∗(/wj−1(θj , wj ) =(θ∗ , w∗ )w.p. r(θj−1 , wj−1 ) w.p. 1 − routput: θ1:J
Assumption 3.1. The Lebesgue measure of Θ ⊂ Rnθ is strictly positive,q : Θ → R++ is a probability density function on Θ, and g : Θ → R++ is an
integrable Borel function.
In what follows, when we refer to an algorithm like Algorithm 2, we
explicitly mention arguments like q and g only if necessary. The following
result is well known, but we state it for reference purposes.
Proposition 3.1. The Markov chain {θj }∞j=1 generated by Algorithm 2 (withJ = ∞) has a unique stationary distribution, which has a density representation π(θ) satisfying
π(θ) ∝ g(θ).

(11)

Furthermore, the distribution of θj converges to π(θ) in total variation asj ↑ ∞.
Proof. This result follows from Theorem 1 and Corollary 2 in Tierney (1994)
under Assumption 3.1.
Let us now consider an extension of iMH (Alg. 2). Algorithm 3 (Martino et al., 2016) extends iMH with I “parameter particles,” and one of theparameter particles is resampled according to the weights wi in line 3. The
average of these weights is used to compute the acceptance probability.
Martino et al. (2016) argue that Algorithm 3 is a variant of the MTM
method with an independent proposal. Based on their argument, we regard
7

Algorithm 3: MTiMH(q, g) (Multiple-Try iMH)
1 for j = 1 : J
2for i = 1 : I
3
θi ∼ q(θ), wi = g(θi )/q(θi )∗i∗ ∼ D(w1:I ), θ∗ = θiP5w∗ = I −1 Ii=1 wi
6
run lines 3–7 of iMH(q, g) (Alg. 2)
output: θ1:J
4

Algorithm 3 as a multiple-try extension of iMH (Alg. 2), and call it MTiMH(multiple-try iMH). It turns out that MTiMH (Alg. 3) is a special case ofiMH (Alg. 2) with an extended parameter space. This observation is usefulin showing theoretical properties of MTiMH; a similar approach is used byAndrieu et al. (2010) to establish various properties of PMMH (Alg. 1) byexpressing it as an MH algorithm with an extended space.
To reformulate MTiMH (Alg 3) as an iMH (Alg. 2) algorithm, let w0 = 0.
Consider w1:I given in line 3; for i ∈ [1 : I], define#
Pi−1 s Pi
s
w
w
.
(12)
, Ps=0Ui = Ps=0
I
I
s
s
w
s=0
s=0 w
In line 4 of MTiMH (Alg. 3), replace “i∗ ∼ D(w1:I )” byu ∼ U (0, 1),

∗

i =

I
X

1{u ∈ Ui }i,

(13)

i=1

where U (0, 1) denotes the uniform distribution on (0, 1), and 1{·} is theindicator function. Note that (13) implies that i∗ ∼ D(w1:I ); thus the distributional properties of MTiMH (Alg. 3) remain unchanged. We also defineθ̃ = (θ1:I , u),

(14)

I

(15)

Θ̃ = Θ × (0, 1),IY

q(θi ),(16)g̃(θ̃) = q̃(θ̃)w∗ ,(17)q̃(θ̃) =i=1

8

where w∗ is given in line 5 of MTiMH (Alg. 3). Consider iMH(q̃, g̃), which is
Algorithm 2 with q̃ and g̃ replacing q and g. Noticing that g̃(θ̃)/q̃(θ̃) = w∗ , itis easy to see that MTiMH(q, g) is equivalent to iMH(q̃, g̃) with i∗ determinedby (13).
It follows from Proposition 3.1 that the stationary density π̃(θ̃) of theMarkov chain {θ̃j }∞j=1 generated by iMH(q̃, g̃) satisfiesπ̃(θ̃) ∝ g̃(θ̃),(18)

and that the distribution of θ̃j converges to π̃(θ̃) in total variation. The
following result is shown by verifying that the marginal distribution of θ∗ (inline 4 of Alg. 3) with respect to π̃(θ̃) is the distribution π(θ) in (11).
Theorem 3.1. Algorithm 3 with J = ∞ generates a sequence {θj }∞j=1 suchthat the distribution of θj converges to π(θ) in total variation as j ↑ ∞.
Proof. See Appendix A.
Our next result shows that as the number I of parameter particles increases to ∞, the acceptance probability for each j ∈ [2 : J] converges to 1
a.s. (almost surely). The result uses Algorithm 4, which is a distributionallyequivalent version of Algorithm 3 in which the random seed is reset in a
specific way and the dependence of all variables on I and j is made explicit.
Theorem 3.2. Consider Algorithm 4. For each j ∈ [2 : J], we have rI,j → 1
a.s. as I ↑ ∞, where rI,j is computed in line 13 of Algorithm 4.
Proof. See Appendix B.
Recall from Subsection 2.1 that θ ∈ Rnθ . Let θ(n) denote the nth component of θ for n ∈ [1 : nθ ]; i.e., θ = (θ(1), . . . , θ(nθ )). Let Cor(θI,j−s (n), θI,j (n0 ))denote the correlation between θI,j−s (n) and θI,j (n0 ) for I ∈ N, j ∈ [1 : J], s ∈[2 : j − 1], and n, n0 ∈ [1 : nθ ]. The following result shows that this correlation converges to 0 as I ↑ ∞ provided that both θI,j−s (n) and θI,j (n0 ) havebounded support.
Theorem 3.3. Consider Algorithm 4. Suppose either that Θ is a boundedset (in which case we set n1 = nθ ), or that there exist sets Θ1 ⊂ Rn1 andΘ2 ⊂ Rn2 with n1 , n2 ∈ N and Θ = Θ1 × Θ2 such that Θ1 is a bounded set.
Then for any j ∈ [2 : J], s ∈ [1 : j − 1], and n, n0 ∈ [1 : n1 ], we havelim Cor(θI,j−s (n), θI,j (n0 )) = 0.

I↑∞

9

(19)

Algorithm 4: MTiMH(I) (MTiMH for Varying I)
1 reset the random seed to any constant independent of I
2 for j = 1 : J
3uj ∼ U (0, 1)
for j = 1 : J
5reset the random seed to uj
6
for i = 1 : I
i
i
i
i
7
θI,j
∼ q(θ), wI,j
= g(θI,j)/q(θI,j
)
4

∗

8
9
10
11
12
13
14

i
∗
1:I
= θI,j), θI,ji∗ ∼ D(wI,j
P
i
∗
wI,j
= I −1 Ii=1 wI,j
if j = 1 then∗∗)
, wI,1(θI,1 , wI,1 ) = (θI,1else∗
rI,j = 1 ∧ wI,j
/w
( I,j−1∗∗(θI,j, wI,j)w.p. rI,j(θI,j , wI,j ) =(θI,j−1 , wI,j−1 ) w.p. 1 − rI,j
output: θI,1:J

10

Algorithm 5: iPMMH(q) (Independent PMMH)
1 for j = 1 : J
2θ∗ ∼ q(θ)
3
run lines 3–14 of PMMH (Alg. 1)
4w∗ = h∗ /q(θ∗ )5
if j = 1 then6(θ1 , ξ1 , w1 ) = (θ∗ , x∗1:T , w∗ )8else
r = 1 ∧ w∗ /w(
j−1

9(θj , ξj , wj ) =7(θ∗ , x∗1:T , w∗ )w.p. r(θj−1 , ξj−1 , wj−1 ) w.p. 1 − routput: θ1:J , ξ1:J
Proof. See Appendix C.
This result shows that as far as bounded parameter components are concerned, autocorrelation of any finite order converges to 0 as I ↑ ∞. As
for unbounded components (in Θ2 ), one can, for example, apply a monotone transformation to make them bounded and apply Theorem 3.3 to thetransformed components.

4

A Multiple-Try Extension of iPMMH

The results shown in Section 3 apply to an MT extension of any schemethat can be reformulated as an iMH algorithm (Alg. 2). In this section wereturn to the framework of state space models discussed in Subsection 2.1
and extend a special case of the PMMH algorithm (Alg. 1). In particular we
consider the PMMH algorithm with an independent proposal, which is given
in Algorithm 5 as a function of the proposal q. The idea here is to express
this algorithm as an iMH algorithm, and to apply Theorems 3.1–3.3 to its
MT extension.
Throughout this section we assume the following in addition to the setup
specified in Subsection 2.1.
Assumption 4.1. For all θ ∈ Θ, we have p(θ), q(θ) ∈ (0, ∞). For all θ ∈ Θ
and x1 ∈ X , we have p(x1 |θ), q(x1 |θ) ∈ (0, ∞). For all θ ∈ Θ, t ∈ [1 : T − 1],11

Algorithm 6: MTiPMMH(q) (Multiple-Try iPMMH)
1 for j = 1 : J
2for i = 1 : I
3
θi ∼ q(θ)
4
run lines 3–14 of PMMH (Alg. 1) with θ∗ = θi
5wi = h∗ /q(θi ), xi1:T = x∗1:T∗∗i∗ ∼ D(w1:I ), θ∗ = θi , x∗1:T = xi1:TP7w∗ = I −1 Ii=1 wi
8
run lines 5–9 of iPMMH (Alg. 5)
output: θ1:J , ξ1:J
6

and xt , xt+1 ∈ X , we have p(xt+1 |xt , θ) ∈ (0, ∞). For all θ ∈ Θ, t ∈ [1 : T ],and xt ∈ X , we have p(yt |xt , θ) ∈ (0, ∞). Finally p(y1:T ) ∈ (0, ∞).
Let us now reformulate iPMMH(q) as an iMH (Alg. 2) algorithm. First
we rearrange state particles xkt and relabeling them as ztk as in Figure 1. To
be more precise, for k, k ∗ ∈ [1 : K], define∗
if k = 1,k∗µ(k, k ) = k − 1 if 2 ≤ k ≤ k ∗ ,
(20)


∗
k
if k > k .
For t ∈ [1 : T ] and k ∈ [1 : K], defineµ(k,kt∗ )ztk = xt

.

(21)

k∗

Note that zt1 = xt t , that ztk = xtk−1 if k ∈ [2 : kt∗ ], and that ztk = xkt if k > kt∗ ;see Figure 1 again. We define
1:K),θ̃ = (θ, z1:T

(22)

Θ̃ = Θ × (X K )T ,

(23)

1:K
q̃(θ̃) = q(θ)p̃(z1:T|θ, y1:T ),
1:K
1:Kg̃(θ̃) = p(θ)u(θ, z1:T , y1:T )p̃(z1:T|θ, y1:T ),

(24)
(25)

1:K
1:K
where p̃(z1:T|θ, y1:T ) is the conditional density of z1:T
given θ and y1:T , andu(·, ·, ·) is defined by (9).

12

x31

x32

x33

x34

x35

x21

x22

x23

x24

x25

x11

x12

x13

x14

x15

∗
(a) State particles xkt with k1:5
= {2, 2, 3, 1, 3}

x31

x32

x23

x34

x25

x11

x12

x13

x24

x15

x21

x22

x33

x14

x35

(b) State particles xkt rearranged with x∗1:t at the bottom

z31

z32

z33

z34

z35

z21

z22

z23

z24

z25

z11

z12

z13

z14

z15
µ(k,kt∗ )(c) State particles relabeled as ztk = xt

1:K
Figure 1: Transformation from x1:K
1:T to z1:T through (21)13Consider iMH(q̃, g̃), which is Algorithm 2 with q = q̃ and g = g̃. It iseasy to see that q̃(θ̃) ∈ (0, ∞) for all θ̃ ∈ Θ̃ by Assumption 4.1. Note thatgiven θ ∈ Θ we have
Z
1:K
1:K
1:K
u(θ, z1:T, y1:T )p̃(z1:T|θ, y1:T )dz1:T(26)
1:K, y1:T )|θ, y1:T ]= E[u(θ, z1:T

(27)

= E[u(θ, x1:K
1:T , y1:T )|θ, y1:T ]

(28)

Z
=

1:K
1:K
u(θ, x1:K
1:T , y1:T )p(x1:T |θ, y1:T )dx1:T= p(y1:T |θ),

(29)
(30)

1:K
, y1:T ) =where E[a|b] is the expectation of a given b, (28) holds since u(θ, z1:T
1:Ku(θ, x1:T , y1:T ); and (30) uses (8). Note from (25)–(30) thatZZg̃(θ̃)dθ̃ = p(θ)p(y1:T |θ)dθ = p(y1:T ) < ∞,(31)

where the last inequality holds by Assumption 4.1. It follows that iMH(q̃, g̃)satisfies Assumption 3.1.
Note from (24) and (25) that
1:K
1:K
)
g̃(θ∗ , z1:Tp(θ∗ )u(θ∗ , z1:T
, y1:T )h∗=== w∗ ,
1:Kq(θ∗ )q(θ∗ )q̃(θ∗ , z1:T)

(32)

where h∗ is given by line 14 of PMMH (Alg. 1), and w∗ is given by line1(as
4 of iPMMH(q) (Alg. 5). It follows that the sequences of θ and z1:T
components of θ̃) generated by iMH(q̃, g̃) are identical to those of θ and x∗1:T
generated by iPMMH(q) (Alg. 5).
Note from Proposition 3.1 that iMH(q̃, g̃) generates a sequence {θ̃j }∞j=1such that the distribution of θ̃j converges in total variation to a unique stationary distribution π̃(θ̃), which satisfies π̃(θ̃) ∝ g̃(θ̃). Theorem 4 in Andrieu
1
et al. (2010) implies that the marginal distribution of (θ, z1:T
) = (θ, x∗1:T ) withrespect to π̃(θ̃) is the desired posterior p(θ, x1:T |y1:T ). By Theorem 3.1, thestationary distribution and convergence properties of iMH(q̃, g̃) are inheritedby MTiMH(q̃, g̃) (recall Alg. 3).
Now consider Algorithm 6, which is a multiple-try extension of iPMMH(q).
We call it MTiPMMH (multiple-try PMMH). In this algorithm, I parameterparticles are sampled, and for each of them, K state particles are sampled for14

Algorithm 7: MTiPMMH(I) (MTiPMMH for Varying I)
1 reset the random seed to any constant independent of I
2 for j = 1 : J
3uj ∼ U (0, 1)
for j = 1 : J
5reset the random seed to uj
6
for i = 1 : I
7
θi ∼ q(θ)
8
run lines 3–14 of PMMH (Alg. 1) with θ∗ = θii= h∗ /q(θi ), xi1:T = x∗1:T
9
wI,j
4

∗

10
11
12
13
14
15
16

∗

1:I
∗
∗
i∗ ∼ D(wI,j), θI,j= θi , ξI,j
= xi1:TPi∗= I −1 Ii=1 wI,j
wI,j
if j = 1 then∗∗)
, wI,j(θI,1 , ξI,1 , wI,1 ) = (θ∗ , ξI,jelse∗
rI,j = 1 ∧ wI,j
/wI,j−1
(
∗
∗
∗
(θI,j
, ξI,j, wI,j)w.p. rI,j(θI,j , ξI,j , wI,j ) =(θI,j−1 , ξI,j−1 , wI,j−1 ) w.p. 1 − rI,j
output: θI,1:J , ξI,1:J

each t ∈ [1 : T ]. In line 6, one parameter particle is resampled according to
the weights wi given in line 5. Note that the structure of MTiPMMH is essentially the same as that of MTiMH (Alg. 3). Indeed, since iPMMH(q) canbe expressed as iMH(q̃, g̃) as discussed above, it follows that MTiPMMH(q)can be expressed as MTiMH(q̃, g̃). Thus we obtain the following results byapplying Theorems 3.1–3.3 to MTiMH(q̃, g̃).
Theorem 4.1. Algorithm 6 (with J = ∞) generates a sequence {(θj , ξj )}∞j=1such that the distribution of (θj , ξj ) converges to p(θ, x1,T |y1:T ) in total variation.
Theorem 4.2. Consider Algorithm 7. For each j ∈ [2 : J], we have rI,j → 1
a.s. as I ↑ ∞, where rI,j is computed in line 15 of Algorithm 7.

15

Theorem 4.3. Consider Algorithm 7. Suppose either that Θ is a boundedset (in which case we set n1 = nθ ), or that there exist sets Θ1 ⊂ Rn1 andΘ2 ⊂ Rn2 with n1 , n2 ∈ N and Θ = Θ1 × Θ2 such that Θ1 is a bounded set.
Then for any j ∈ [2 : J], s ∈ [1 : j − 1], and n, n0 ∈ [1 : n1 ], we havelim Cor(θI,j−s (n), θI,j (n0 )) = 0.

I↑∞

(33)

Theorem 4.1 can be shown more directly by following Andrieu et al.’s(2010, Theorem 4) proof or Wilkinson’s (2012, pp. 305–307) argument. It
is easy to see that Algorithm 7 is distributionally equivalent to Algorithm
6; recall Algorithm 4. With additional notation, Theorem 4.3 can easily be
extended to Θ × X T by replacing Θ with Θ × X T . Overall Theorems 4.1–4.3
indicate that as the number I of parameter particles increases to ∞, one canobtain almost i.i.d. samples from the posterior p(θ, x1:T |y1:T ). This is true
regardless of the proposal and the number K of state particles.
Before we conclude this section, let us briefly discuss how MTiPMMH(Alg. 6) can be parallelized. Note that lines 3–5 can easily be processed inparallel for all i ∈ [1 : I], but such a parallel scheme is rather inefficientsince for each j ∈ [1 : J], before line 6 can be executed, lines 3–5 must becompleted for all i ∈ [1 : I]. A more efficient but memory-intensive parallel
scheme is offered in Algorithm 8. In this algorithm, lines 2–5 can be processedin parallel for all i ∈ [1 : I], and the resampling and updating steps in lines
7–9 are executed using the data generated in lines 1–5. Depending on thecomputing environment, one can parallelize lines 1–5 even more efficiently
since in effect it suffices to draw I × J i.i.d. samples from q(θ) and run lines
4–5 for each sample.

5

Numerical Experiments

5.1

A Simple Stochastic Volatility Model

To illustrate our theoretical results, we consider the simple stochastic volatility model used by Fearnhead and Meligkotsidou (2015):
xt = γxt−1 + σx ηt , ηt ∼ N (0, 1),
yt = σy exp(xt )t , t ∼ N (0, 1),

(34)
(35)

where N (µ, σ 2 ) denotes the normal distribution with mean µ and varianceσ 2 . The initial value x1 is drawn from N (0, 1).
16

Algorithm 8: MTiPMMH for Parallelism
1 for i = 1 : I
2for j = 1 : J
3θji ∼ q(θ)
4
run lines 3–14 of PMMH (Algorithm 1) with θ∗ = θji∗
5wji = h∗ /q(θji ), xi,j
1:T = x1:T
for j = 1 : J∗∗
7i∗ ∼ D(wj1:I ), θ∗ = θji , x∗1:T = xi1:T,jP
8w∗ = I −1 Ii=1 wji
9
run lines 5–9 of iPMMH (Alg. 5)
output: θ1:J , ξ1:J
6

Though it is more common in the literature to divide xt by 2 in (35) (e.g.,
Shephard and Pitt, 2007), we simply use Fearnhead and Meligkotsidou’s(2015) setting, where all parameter values including priors are fully specified(thus we do not need to choose them arbitrarily). Following Fearnhead andMeligkotsidou (2015, Section 4.1.2), we fix T = 1000, and assume that thepriors for γ, σx2 , and σy2 are independent and given byγ ∼ N (0.9, 0.1)|(−1,1) ,1/σx2 ∼ Γ(1, 1/100),1/σy2 ∼ Γ(1, 1),(36)

where the distribution for γ is truncated to (−1, 1), and Γ(a, b) denotes the
gamma distribution with shape parameter a and scale parameter b.
Given simulated data y1:T (to be specified below), we consider two different proposals on θ = (γ, σx2 , σy2 ). The first proposal is simply the prior p(θ)
given by (36). We construct the second proposal q̂(θ) using samples from apilot run of a standard MCMC algorithm; see Appendix D for details.
For each of the two proposals, we consider two values, 4400 and 110000,
for J, where J is the number of iterations. If J = 4400, then we discard
the first 400 iterations as burn-in. If J = 110000, we discard the first 10000
iterations as burn-in. Note however that our theoretical results require noburn-in. In total we consider the four pairs of q(θ) and J given byq(θ) ∈ {p(θ), q̂(θ)},

J ∈ {4400, 110000}.

(37)

In each case, we increase I from 1 to I one by one, where I = 500 if J = 4400
17

and I = 10 if J = 110000. We fix K = 500 throughout, where K is thenumber of state particles; recall PMMH (Alg. 1).
Since it is costly to literally implement Algorithm 8 for many different
values of I, we modify it to reuse as much data as possible. The modified
version is shown in Algorithm 9. In lines 1–2 of this algorithm, basic data
1:I
1:I
, w1:J, and x1:I,1:J
for the largest possible value of I are generated andθ1:J
1:T
stored. Using the stored data, the steps for resampling i∗ and updating θjsummarized in line 4 are executed for each I = 1, . . . , I. This way we generate basic data only once, saving a significant amount of time. In addition,
this algorithm is consistent with Theorems 4.2 and 4.3; in fact these results
directly apply to Algorithm 9 with I = ∞.
We defineµ = ln σy2 ,βx = 1/σx2 .

(38)

Following Fearnhead and Meligkotsidou (2015), we present our simulation
results in terms of µ, βx , and γ.
In implementing Algorithm 9 (or any of our MTiPMMH algorithms), it
is useful to introduce an additional parameter, say M , to line 14 of PMMH
(Alg. 1) as follows:
!
T
K
Y
X
M K −1wtk .
(39)h∗ = p(θ∗ )
t=1

k=1

Theoretically the new parameter M cancels out in line 8 of iPMMH(q) (Alg.
5) as long as M > 0. In practice, however, it often enables us to avoid “0/0”by choosing an appropriate value for M . Note that the right-hand side of(39) involves the product of T expressions. Since T = 1000 here, if M = 1,
the product can easily evaluate to 0 numerically even if it is strictly positive
theoretically. This problem can mostly be avoided if I = 1, in which casemost computations can be done in terms of logarithms. This is not the casehere since the sum of I terms in line 8 of Algorithm 8 must be evaluated.
Accordingly, we replace line 14 of PMMH (Alg. 1) with (39). We set M = 5,
which seems to work well in the current setting (including (40) below).

5.2

Simulation Results

We simulated data x1:T and y1:T with the following parameter values:
γ = 0.99,

σx2 = (1 − γ 2 ),18

σy2 = 1.

(40)

The above values were borrowed from Fearnhead and Meligkotsidou (2015,
Section 4.1.2). The simulated data are common to all the results presented
below.
Figures 2–5 illustrate how the performance of Algorithm 8 depends on
I by implementing Algorithm 9 for the four different cases given by (37).
Since Algorithm 9 heavily uses lines from other algorithms, we use Algorithm
7 instead to explain these figures. In each of Figures 2–5, the top panel
plots the average of acceptance probabilities rI,j in line 15 of Algorithm 7
over j ∈ [1 : J]. The second panel plots the realized overall acceptance
rate in each run, i.e., the number of acceptance (or updating) events inline 16 of Algorithm 7 divided by J. The remaining three panels plot the
autocorrelation times for µ, βx , and γ, respectively.
The data for Figure 2 were generated by setting q(θ) = p(θ) and J = 4400
and by varying I from 1 to 500 one by one. The top panel shows thatthe average acceptance probability tends to increase as I increases. Thisobservation is consistent with Theorem 4.2, which implies that the average
acceptance probability converges a.s. to 1 as I increases to ∞. The second
panel shows that the realized overall acceptance rate in each run is essentially
identical to the corresponding average acceptance probability. The remaining
three panels show that the autocorrelation times for µ, βx , and γ quickly
decline as I increases from 1, and have a tendency to decrease as I further
increases. This observation is consistent with Theorem 4.3, which shows thatthe autocorrelation of any order for each parameter with bounded support
converges to 0 as I increases to ∞. Although µ and βx do not have bounded
support, their autocorrelation times also tend to decrease as I increases.
The data for Figure 3 were generated by setting q(θ) = q̂(θ) and J =4400 and by varying I from 1 to 500. Since, as mentioned above, q̂(θ) wasconstructed using samples from a pilot run of a standard MCMC sampler,q̂(θ) is expected to be substantially closer to the posterior p(θ|y1:T ) than theprior p(θ) is. One can see this effect in the top panel in Figure 3, which showsthat the average acceptance probability quickly rises as I increases from 1,
and is considerably higher for any I than the corresponding value in Figure
2. As in Figure 2, the bottom three panels show that the autocorrelation
times for µ, βx , and γ quickly decline as I increases from 1.
In Figures 2 and 3, due to a memory constraint, we had to choose a rather
small value for J in order to allow I to increase up to 500. In Figures 4 and
5, we chose a larger value for J by restricting the largest value for I to 10.
These figures exhibit patterns similar to those observed in Figures 2 and 3
19

Figure 2: q(θ) = p(θ), J = 4400

20

Figure 3: q(θ) = q̂(θ), J = 4400

21

Figure 4: q(θ) = p(θ), J = 11000022

Figure 5: q(θ) = q̂(θ), J = 11000023

Algorithm 9: MTiPMMH for Different Values of I
1 for i = 1 : I
2run lines 2–5 of Algorithm 8
for I = 1 : I
4run lines 6–9 of Algorithm 8
5(θI,1:J , ξI,1:J ) = (θ1:J , ξ1:J )output: θ1:I,1:J , ξ1:I,1:J
3

even if J = 110000 and I is varied from 1 to only 10. The figures suggestthat one can significantly improve mixing if one increases I from 1 to 10,
regardless the proposal used. This is useful since the time required to run
Algorithm 8 on a multithreaded computer does not increase significantly as
I is increased from 1 up to the number of threads available. Thus the timerequired to run Algorithm 8 with I = 10 is not substantially longer than thetime required to run the algorithm with I = 1, which is simply iPMMH.
Figure 6 shows the ACF (autocorrelation function) plots, trace plots, andhistograms for I = 1 and I = 500 using the same data as used for Figure
2. Figure 6(a) indicates that iPMMH suffers from high autocorrelation andpoor mixing. This is mostly because the prior p(θ) here is not an efficient
proposal. Even with this poor proposal, if I is chosen to be sufficiently
large, MTiPMMH exhibits low autocorrelation and good mixing, as shownin Figure 6(b). The histograms are rather rough here since J is only 4400.
On a standard PC, it may be easier to use a large value for J and asmall value for I. It may also be possible to significantly improve mixing
by constructing a proposal using samples from a pilot run. This is indeed
the case here, as shown in Figure 7(a), which uses the same data as used
for Figure 5. With q(θ) = q̂(θ) and J = 110000, iPMMH exhibits fairly low
autocorrelation and good mixing. Figure 7(b) shows that one can visibly
improve the performance by using MTiPMMH with I = 10. As discussed
above, this does not significantly increase computation time if Algorithm 8
is implemented on a PC with at least 10 threads.

24

(a) I = 1(b) I = 500
Figure 6: q(θ) = p(θ), J = 4400: ACF plots, trace plots, and histograms(dotted lines indicate true values)
25

(a) I = 1(b) I = 10
Figure 7: q(θ) = q̂(θ), J = 110000: ACF plots, trace plots, and histograms(dotted lines indicate true values)
26

6

Discussion

In this paper we have proposed a multiple-try extension of the PMMH algorithm with an independent proposal, which we call MTiPMMH. In thisalgorithm, I ∈ N parameter particles are sampled from the independent proposal. For each of them, a particle filter with K ∈ N state particles is run.
Then one of the parameter particles is resampled according to their marginal
likelihood estimates. The state of the chain is updated with a probability
based on the ratio between the averages of marginal likelihood estimates.
We have shown that MTiPMMH has the following three asymptotic properties: (i) the distribution of the Markov chain generated by the algorithm
converges to the posterior of interest in total variation; (ii) as I increases to∞, the acceptance probability at each iteration converges to 1 a.s.; and (iii)
as I increases to ∞, the autocorrelation of any order for any parameter withbounded support converges to 0. These results indicate that MTiPMMH
generates almost i.i.d. samples from the posterior for sufficiently large I.
It is well known that as the number of iterations increases to ∞, theMarkov chain generated by a standard MCMC algorithm converges in distribution to the posterior of interest. Our results offer another dimension inwhich one can approach the correct posterior. As they are valid for any proposal (satisfying regularity conditions) and any number K of state particles,
we believe that our asymptotic results become increasingly relevant as parallel computing with supercomputers, GPUs, and clouds becomes increasingly
accessible.
To show our theoretical results, we reformulated iPMMH as an iMH algorithm. The same procedure can in fact be used to reformulate PMMH as an
MH algorithm. Thus one can apply the standard MTM method to this MH
algorithm to obtain a multiple-try extension of PMMH without assuming an
independent proposal. However, the practical value of this extension is notnecessarily clear if one wishes to obtain almost i.i.d. samples asymptotically.
Note that if one is to obtain i.i.d. samples, then the acceptance probability
must always be equal to 1. This implies that the proposal must be effectively
independent; for otherwise it would inevitably cause autocorrelation.
Although our theoretical results require I to be sufficiently large, ournumerical experiments suggest that one can visibly improve mixing by increasing I from 1 to a small number like 10. This does not significantly
increase computation time if a computer with at least the same number of
threads is used. One can further improve the overall performance of the27

algorithm by constructing a proposal using samples from a pilot run of astandard MCMC algorithm. Our algorithm, which assumes an independent
proposal, makes it easy to utilize such a proposal.

Appendix A

Proof of Theorem 3.1

Consider iMH(q̃, g̃) (Alg. 2). We start by verifying Assumption 3.1 for q̃ andg̃. Note from (16) that q̃ : Θ̃ → R++ is a probability density. We haveZ
Z Zg̃(θ̃)dθ̃ =q̃(θ1:I , u)w∗ du dθ1:I
(41)
#
Z Z "
I
X
=
q̃(θ1:I , u)w∗1{u ∈ Ui } du dθ1:I
(42)
i=1

=

I Z
X
i=1
I Z
X

Z

q̃(θ1:I , u)w∗ 1{u ∈ Ui }du dθ1:I

" I
Y

#

wi 1:I
dθ
Iw∗
s=1
i=1
" I
#
I Z
X
Y
g(θi ) 1:I= I −1dθq(θs )q(θi )
s=1
i=1


I Z
I
X
Y

= I −1q(θs ) g(θi )dθ1:I
=

i=1

=I

−1

(43)

q(θ ) w∗
s

(44)

(45)

(46)

s∈[1:I],s6=i

I ZXg(θi )dθi < ∞,

(47)

i=1

P
where Ii=1 1{u ∈ Ui } = 1 in (42), (44) uses (16) and (12) and (17), (45)uses line 3 of Algorithm 3, and the last inequality holds by Assumption 3.1.
It follows that iMH(q̃, g̃) satisfies Assumption 3.1. Thus by Proposition 3.1,
it generates a sequence {θ̃j }∞j=0 such that the distribution of θ̃j converges tothe stationary density π̃(θ̃) satisfying (18) as j ↑ ∞.
Let p(θ∗ ) be the marginal density of θ∗ with respect to π̃(θ̃). We claimthatp(θ∗ ) ∝ g(θ∗ ).
28

(48)

To show this, let B ⊂ Θ be a Borel set. Recall from line 4 of MTiMH(q, g)∗(Alg. 3) that θ∗ = θi . Let P {θ∗ ∈ B} be the probability of θ∗ ∈ B withrespect to π̃(θ̃). Note from (18) that
Z Z∗∗
P {θ ∈ B} ∝g̃(θ1:I , u)1{θi ∈ B}du dθ1:I
(49)
Z Z
=

I
X

[1{i∗ = i}1{θi ∈ B}]du dθ1:I(50)

g̃(θ1:I , u)1{u ∈ Ui }1{θi ∈ B}du dθ1:I

(51)

1:I

g̃(θ , u)

i=1

=

I Z Z
X
i=1
I Z
X

" I
Y

#

wi
1{θi ∈ B}dθ1:I
∗
Iw
s=1
i=1
" I
#
Z
I
X
Y
g(θi )
−1
s
=I
q(θ )
1{θi ∈ B}dθ1:I
i
q(θ )
s=1
i=1
I Z
X
−1
=I
g(θi )1{θi ∈ B}dθi=q(θs ) w∗

(52)

(53)

(54)

i=1

Z
=

g(θ∗ )1{θ∗ ∈ B}dθ∗ ,(55)

where (52) uses (12) and (17), and (53) uses line 3 of Algorithm 3. Since(49)–(55) hold for any Borel set B ⊂ Θ, (48) follows.
To see that the marginal distribution of θj converges to π(θ) in total
variation, let π̃j (θ̃) be the distribution of θ̃j for j ∈ N. By Proposition 3.1,π̃j converges to π̃ in total variation as j ↑ ∞. Let j ∈ N, and let B ⊂ Θ bea Borel set. Let Pj {θ∗ ∈ B} be the probability of θ∗ ∈ B with respect toπ̃j (θ̃). We have
Pj {θ∗ ∈ B} − P {θ∗ ∈ B}
ZZ∗
= π̃j (θ̃)1{θ ∈ B}dθ̃ − π̃(θ̃)1{θ∗ ∈ B}dθ̃Z= (π̃j (θ̃) − π̃(θ̃))1{θ∗ ∈ B}dθ̃Z≤ |π̃j (θ̃) − π̃(θ̃)|dθ̃ → 0 as j ↑ ∞,
29

(56)
(57)
(58)
(59)

where the convergence holds by the total variation convergence of π̃j to π̃.
Since the right-hand side of the inequality in (59) does not depend on B, wehave
lim sup |Pj {θ∗ ∈ B} − P {θ∗ ∈ B}| = 0,

j↑∞ B

(60)

where the supremum is taken over all Borel sets B ⊂ Θ̃. It follows that themarginal distribution of θj converges to π(θ) in total variation.

Appendix B

Proof of Theorem 3.2

Note that the random seed is reset in lines 1 and 5 of Algorithm 3 in such a
way that for any I0 , I1 ∈ N with I0 ≤ I1 ,∀i ∈ [1 : I0 ],i∈[1:I ]wIi 1 ,j = wIi 0 ,j ,(61)i∈[1:I ]where {wIi 0 ,j }j∈[1:J]0 and {wIi 1 ,j }j∈[1:J]1 are given by line 7 of MTiMH(I) (Alg.
4) with I = I0 and I = I1 , respectively. In words, by increasing I from I0ito I1 , the values of wI,jchosen for I = I0 remain unchanged. We need thisproperty to apply the strong law of large numbers below.
Fix j ∈ [1 : J]. Note that for any I ∈ N,i∈[1:I]i{θI,j}j∈[1:J] is i.i.d. across i ∈ [1 : I] and j ∈ [1 : J].

Let I ∈ N. By (62), for any i ∈ [1 : I] we have
Z
Z
g(θ)
i
q(θ)dθ = g(θ)dθ < ∞,EwI,j =
q(θ)

(62)

(63)

where E is the expectation operator, and the inequality holds by integrability of g. Note from lines 11 and 14 of MTiMH(I) (Alg. 4) that wI,j−1 ≤j−1∗max{wI,j0 }j 0 =1 ; thus
∗
wI,j
j
∗
max{wI,j0 }j 0 =1≤ rI,j ≤ 1.

(64)

∗
i
For each j 0 ∈ N with j 0 ≤ j, we have wI,j
0 → EwI,j 0 a.s. as I ↑ ∞ by thestrong law of large numbers and (63). Hence letting I ↑ ∞ in (64), we haverI,j → 1 a.s.

30Appendix C

Proof of Theorem 3.3

If Θ is a bounded set, we set Θ1 = Θ. Since Θ1 is a bounded subset of Rn1 ,there exists θ ∈ (0, ∞) such that∀θ ∈ Θ,max |θ(n)| ≤ θ.

(65)n∈[1:n1 ]For I, j ∈ N and n ∈ [1 : n1 ], we define∗µI (n) = EθI,j
(n),

(66)

µI,j (n) = EθI,j (n),∗νI (n) = E(θI,j(n))2 ,

(67)
(68)

νI,j (n) = E(θI,j (n))2 ,
i
i
EwI,j
θI,j
(n)
µ(n) =
,
i
EwI,j

(69)

ν(n) =

i
i
(n))2
(θI,j
EwI,j
.
i
EwI,j

(70)
(71)

Note from (62) that µI and νI do not depend on j, and that µ and ν depend
neither i nor j.

C.1

Preliminary Lemmas

In this subsection we show some properties of the functions defined in (66)–(71).
Lemma C.1. For any n ∈ [1, n1 ] we have(a) lim µI (n) = µ(n),(b) lim νI (n) = ν(n).

I↑∞

I↑∞

(72)

Proof. Fix n ∈ [1, n1 ]. Throughout the proof we suppress dependence on n.
Line 8 of Algorithm 4 implies that
#
" PI
i
i
w
θ
i=1 I,j I,j.
(73)µI = EPIi
i=1 wI,j31

Note from (65) thatPIiii=1 wI,j θI,ji
i=1 wI,jPI
From (73) we have"
µI = E

I −1

PI

i
i
i=1 wI,j θI,jP
i
I −1 Ii=1 wI,j

#

≤ θ.

i
i
EwI,j
θI,j
→
= µ as I ↑ ∞,
i
EwI,j

(74)

(75)

where the convergence holds by the strong law of large numbers, (63), thedominated convergence theorem, and (74). We have verified (72)(a). Weii.
)2 replacing θI,jobtain (72)(b) by the same argument with (θI,jLemma C.2. For any j ∈ [1 : J] we havelim E(1 − rI,j ) = 0.

I↑∞

(76)

Proof. By Theorem 3.2, we have rI,j → 1 a.s. as I ↑ ∞. Hence 1 − rI,j → 0
a.s. as I ↑ ∞. Since 1 − rI,j ≤ 1 for all I ∈ N, we obtain (76) by thedominated convergence theorem.
Lemma C.3. For any n ∈ [1, n1 ] and j ∈ [1 : J] we have(a) lim µI,j (n) = µ(n),(b) lim νI,j (n) = ν(n).

I↑∞

I↑∞

(77)

Proof. Fix n ∈ [1, n1 ]. Throughout the proof we suppress dependence on∗. Hencen. Let I ∈ N. Note from line 11 of Algorithm 4 that θI,1 = θI,1
µI,1 = µI ; thus (77)(a) holds for j = 1 by Lemma C.1.
Let j ∈ [2 : J] and I ∈ N. From line 14 of Algorithm 4, we have∗
µI,j = E[(1 − rI,j )θI,j−1 + rI,j θI,j
].

(78)

By the triangle inequality, we have|µI,j − µ| ≤ |µI,j − µI | + |µI − µ|∗= |E(1 − rI,j )(θI,j−1 − θI,j)| + |µI − µ|

(79)
(80)

≤ E(1 − rI,j )2θ + |µI − µ|,(81)

where (80) uses (66). Hence (77)(a) holds by Lemmas C.2 and C.1.
∗ 2
We obtain (77)(b) by the same argument with (θI,j )2 and (θI,j) replacing∗θI,j and θI,j .
32

Lemma C.4. For any n ∈ [1, n1 ] and j ∈ [1 : J] we havelim V ar(θI,j (n)) = ν(n) − µ(n)2 > 0,

I↑∞

(82)

where V ar(θI,j (n)) is the variance of θI,j .
Proof. Fix n ∈ [1, n1 ]; we suppress dependence on n. Let j ∈ [1 : J] andI ∈ N. Recalling (69) and (67), we haveV ar(θI,j ) = E(θI,j )2 − (EθI,j )2 = νI,j − (µI,j )2 .

(83)

Hence the equality in (82) holds by Lemma C.3. To see the inequality, letii. Then in view of (70) and (71), it is straightforward to show
/EwI,jω i = wI,jthatiiEω i (θI,j− Eω i θI,j
)2 = ν − µ2 .

(84)i
6= µ a.s. by
The left-hand side is strictly positive since ω i > 0 and θI,j
Assumption 3.1. Thus the desired inequality follows.

C.2

Completing the Proof of Theorem 3.3

Fix j ∈ [2 : J], s ∈ [1 : j − 1], and n, n0 ∈ [1 : n1 ]. For I ∈ N definezI,j−s = θI,j−s (n) − µI,j−s (n),0zI,j−1= θI,j−1 (n0 ) − µI,j−1 (n0 ),

(85)
(86)

0
zI,j
= θI,j (n0 ) − µI,j (n0 ),

(87)

∗
∗
zI,j
= θI,j(n0 ) − µI,j (n0 ),0βI,j,s = EzI,j−s zI,j.

(88)

Cor(θI,j−s (n), θI,j (n0 )) = βI,j,s /[V ar(θI,j−s (n))V ar(θI,j (n0 ))]1/2 .

(90)

(89)

For any I ∈ N we have

By Lemma C.4, as I ↑ ∞, the denominator converges to{[ν(n) − µ(n)2 ][ν(n0 ) − µ(n0 )2 ]}1/2 > 0.

(91)

Thus to conclude (19), it suffices to show thatlim βI,j,s = 0.

I↑∞

33

(92)

To this end, let I ∈ N. Note from line 14 of Algorithm 4 that
0∗βI,j,s = E[(1 − rI,j )zI,j−s zI,j−1+ rI,j zI,j−s zI,j
].

(93)

∗
∗
Since EzI,j−s zI,j
= EzI,j−s EzI,j= 0 by independence and (85), we have∗0∗|βI,j,s | = |βI,j,s − EzI,j−s zI,j
| = |E(1 − rI,j )(zI,j−s zI,j−1− zI,j−s zI,j)|2

≤ E(1 − rI,j )2θ → 0 as I ↑ ∞,

(94)
(95)

where the convergence holds by Lemma C.2. Now (92) follows.

Appendix D

Construction of q̂(θ)

In this appendix we describe the pilot run and the associated proposal q̂(θ)
used in Section 5. Let φ(a|b, c) denote the normal density as a function ofa ∈ R parameterized by b ∈ R and c > 0 as follows:(a − b)2
1
exp −
φ(a|b, c) = √
.
(96)
2c
2πc
Consider an MH sampler for (34) and (35). Initially we draw γ, σx2 , and σy2
from the priors given in (36), x1 from N (0, 1), and x2:T from (34). In each
iteration afterwards, given x1:T , we defineν = σx−2

TX

x2t−1 ,

PT

δ = 10,t=2γ̂ = Pt=2
T

xt xt−12t=2 xt−1.

(97)

Then we apply the MH update to γ, σy2 , σx2 , x1 , . . . , xT one by one with thefollowing sampling schemes:γ̂ν + 0.9δ
1
γ∼φ γ
,
,
(98)
ν+δ
ν + δ (−1,1)
!
PT
2
−1
exp(−2α
x
)]
+
2
[y
T
+
2
t
1/σy2 ∼ Γ
, t=1 t
,
(99)
2
2
!
PT
2
2
2
(x
−
γx
)
+
x
(1
+
γ
)
+
0.02
T
+
2
t
t−1
1
1/σx2 ∼ Γ
, t=2
,
(100)
2
2
x1 ∼ φ(x1 |γx2 , σx2 ),

(101)
34



γ(xt+1 + xt−1 )σx2
xt ∼ φ xt,,(1 + γ 2 )(1 + γ 2 )xT ∼ φ(xT |γxT −1 , σx2 ).

∀t ∈ [2 : T − 1],

(102)
(103)

These distributions are components of the corresponding full conditionals, sothat the acceptance probability for each parameter and state variable takes
a rather simple form.
We run this algorithm for 110000 iterations, discard the first 10000 iterations as burn-in, and construct the proposal q̂(θ) using the remaining
samples as follows. First let mx and vx be the sample mean and variance ofthe 100000 samples of ln σx2 from this pilot run; define my and vy similarly.
We also define mγ and vγ as the sample mean and variance of the 100000
samples of γ. Let
m2γ (1 − mγ )α=
− mγ ,
vγ

β=

α
− α.
mγ

(104)

Now we specify the proposal q̂(θ) byγ ∼ Beta(α, β),ln(σx2 ) ∼ N (mx , vx ),ln(σy2 ) ∼ N (my , vy ),

(105)
(106)
(107)

where Beta(α, β) is the beta distribution with shape parameters α and β.

References
Andrieu, C., Doucet, A., Holenstein, R., 2009, Particle Markov chain Monte
Carlo for efficient numerical simulation, in Monte Carlo and Quasi-Monte
Carlo Methods, L’Ecuyer, P., Owen, A.B., eds., Springer, Berlin, pp. 45–
60.
Andrieu, C., Doucent, A., Holenstein, R., 2010, Particle Markov chain Monte
Carlo methods, Journal of the Royal Statistical Society: Series B, 72,
269–342.
Andrieu, C., Roberts, G.O., 2009, The pseudo-marginal approach for efficient
Monte Carlo computations, Annals of Statistics 2, 697–725.
Beaumont, M., 2003, Estimation of population growth or decline in genetically monitored populations, Genetics 164, 1139–1160.
35

Bédard, M., Douc, R., Moulines, E., 2012, Scaling analysis of mutiple-try MCMC
methods, Stochastic Processes and their Applications 122, 758–786.
Del Moral, P., 2004, Feynman-Kac Formulae: Genealogical and Interacting
Particle Systems with Applications, Springer, New York.
Doucet, A., Godsill, S., Andrieu, C., 2000, On sequential Monte Carlo sampling methods for Bayesian filtering, Statistics and Computing 10, 197–
208.
Fearnhead, P., 2011, MCMC for state-space models, in Handbook of Markov
Chain Monte Carlo, ed. by Brooks, S., Gelman, A., Jones, G.L., Meng,
X.-L., pp. 513–529.
Golightly, A., Henderson, D.A., Sherlock, C., 2015, Delayed acceptance particle MCMC for exact inference in stochastic kinetic models, StatisticalComputing 25, 1039–1055.
Kokkala, J., Särkkä, S., 2015, Combining particle MCMC with Rao-Blackwellized
Monte Carlo data association for parameter estimation in multiple target
tracking, Digital Signal Processing 47, 84–95.
Liu, J.S., Liang, F., Wong, W.H., 2000, The multiple-try method and local
optimization in Metropolis sampling, Journal of the American Statistical
Association 95, 121–134.
Martino, L., Louzada, F., 2016, Issues in the multiple-try Metropolis mixing, Computational Statistics, doi:10.1007/s00180-016-0643-9.
Martino, L., Leisen, F., Corander, J., 2016, On multiple try schemes and theparticle Metropolis-Hastings algorithm, memo.
Pandolfi, S., Bartolucci, F., Friel, N., 2010, A generalization of the MultipleTry Metropolis algorithm for Bayesian estimation and model selection,
Journal of Machine Learning Research: W&CP 9, 581–588.
Pitt, M.K., Silva, R.d.S, Giordani, P., Kohn, R., 2012, On some propertiesof Markov chain Monte Carlo simulation methds based on the particle
filter, Journal of Econometrics 171, 134–151.
Roberts, G.O., Rosenthal, J.S., 2004, General state space Markov chains and
MCMC algorithms, Probability Surveys 1, 20–71.
Shephard, N., Pitt, M.K., 1997, Likelihood analysis of non-Gaussian measurement time series, Biometrika 84, 653–667.
Storvik, G., 2011, On the flexibility of Metropolis-Hastings acceptance probabilities in auxiliary variable proposal generation, Scandinavian Journalof Statistics 38, 342–358.
36Tierney, L., 1994, Markov chains for exploring posterior distributions, Annals of Statistics 22, 1701–1728.
Tran, M.-N., Pitt, M.K., Kohn, R., 2016, Adaptive Metropolis-Hastings sampling using reversible dependent mixture propsals, Statistical Computing
26, 361–381.
Wang, G., 2009, Signal extraction from long-term ecological data using Bayesianand non-Bayesian state-space models, Ecological Informatics 4, 69–75.
Wilkinson, D.J., 2012, Stochastic Modelling for Systems Biology, 2nd Ed.,
CRC Press, Boca Raton.
Yuan, K., Girolami, M., Niranjan, M., 2012, Markov chain Monte Carlo Methods for state-space models with point process observations, Neural Computation 24, 1462–1486.

37