arXiv:1503.00966v2 [math.ST] 10 Apr 2017

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLING:
BOUNDS, ADAPTIVE RESAMPLING VIA ∞-ESS,
AND AN APPLICATION TO PARTICLE GIBBS
JONATHAN H. HUGGINS AND DANIEL M. ROY
Abstract. Sequential Monte Carlo (SMC) algorithms were originally designed
for estimating intractable conditional expectations within state-space models,
but are now routinely used to generate approximate samples in the context of
general-purpose Bayesian inference. In particular, SMC algorithms are often
used as subroutines within larger Monte Carlo schemes, and in this context,
the demands placed on SMC are different: control of mean-squared error is
insufficient—one needs to control the divergence from the target distribution
directly. Towards this goal, we introduce the conditional adaptive resampling
particle filter, building on the work of Gordon, Salmond, and Smith (1993),
Andrieu, Doucet, and Holenstein (2010), and Whiteley, Lee, and Heine (2016).
By controlling a novel notion of effective sample size, the ∞-ESS, we establish
the efficiency of the resulting SMC sampling algorithm, providing an adaptive
resampling extension of the work of Andrieu, Lee, and Vihola (2013). We apply
our results to arrive at new divergence bounds for SMC samplers with adaptive resampling as well as an adaptive resampling version of the Particle Gibbs
algorithm with the same geometric-ergodicity guarantees as its nonadaptive
counterpart.

1. Introduction
2. Preliminaries
3. Conditional αSMC
4. Main Results
5. Bounding the Divergence of SMC Samplers
6. The α-Particle Gibbs Sampler
Appendix A. Additional Proofs
Acknowledgments
References

1
8
10
12
18
20
23
33
33

1. Introduction
Sequential Monte Carlo (SMC) methods are a popular class of algorithms for approximate inference [Doucet et al., 2000, 2001, Del Moral et al., 2006, Kantas et al.,
2009, Doucet and Johansen, 2010, Künsch, 2013]. In the context of Bayesian inference, SMC produces a particle approximation to the posterior distribution as well
as an unbiased estimate of the marginal likelihood. Traditionally, particle approximations were built to estimate conditional expectations, and the analysis of SMC
methods focused on this operator perspective, by bounding the mean squared error
of the resulting estimates.
Date: Tuesday 11th April, 2017.
1

2

J. H. HUGGINS AND D. M. ROY

Increasingly, SMC methods are being used to produce approximate samples,
usually in the inner loop of other approximate inference algorithm. A key example is the class of particle Markov chain Monte Carlo (PMCMC) methods, which
aim to combine the best features of SMC and MCMC approaches by using SMC
as a proposal mechanism for a Metropolis–Hastings (“particle MH”) or approximate Gibbs (“particle Gibbs”) sampler [Holenstein, 2009, Andrieu et al., 2010].
Characterizing the efficiency of PMCMC methods is an active area of investigation [Andrieu and Roberts, 2009, Andrieu and Vihola, 2014, Andrieu et al., 2013,
Chopin and Singh, 2013, Lindsten et al., 2014, Lee and Latuszynski, 2014].
When SMC methods are employed for sampling, convergence guarantees from
the operator perspective are insufficient. In this work, we take up the sampling perspective, and study the distribution of a sample drawn from the SMC particle filter
approximation. Building off the work of Gordon et al. [1993], Andrieu et al. [2010,
2013], and Whiteley et al. [2016], we use conditional filters to derive a minorization
condition, lower bounding the density of the approximate sample’s distribution in
terms of that of the target distribution. The analysis extends to conditional SMC
as well as to adaptive resampling versions. One of our key contributions is a novel
notion of effective sample size, the ∞-ESS, which we use to establish the efficiency
of the adaptive SMC sampling algorithm. Thus, our results are both a sampling
analogue to the operator work of Whiteley et al. [2016] and an adaptive resampling
extension to the sampling work of Andrieu et al. [2013]. We apply our results to
arrive at new divergence bounds for SMC samplers with adaptive resampling as
well as an adaptive resampling version of the Particle Gibbs algorithm with the
same geometric-ergodicity guarantees as its nonadaptive counterpart.
In the remainder of this section we provide an overview of our contributions
for the special case of the the conditional adaptive resampling particle filter: we
introduce the conditional adaptive resampling particle filter, present our main theoretical results characterizing its performance, and describe an application to a novel
adaptive resampling Particle Gibbs algorithm.
1.1. (Conditional) adaptive resampling particle filters. We follow a similar
setup and notation to Del Moral [2004]. Let (ξt )t≥1 be an inhomogeneous Markov
chain on the measurable space (E, E) with transition kernels (Mt )t≥2 and initial
distribution M1 . Denote expectations with respect to the Markov chain by E[·].
Let gt : E → R+ , for t ≥ 1, be a sequence of E-measurable potential functions
Qt
on E, let g0 ≡ 1, and write gs:t (xs:t ) , τ =s gτ (xτ ). For t = 1, 2, . . . , define the
measure π1:t on E t given by
π1:t (dx1:t ) , γ1:t (dx1:t )/Zt ,

(1.1)

where
γ1:t (dx1:t ) ,

t
Y

gs (xs )Ms (xs−1 , dxs ) and Zt , γ1:t (1).

s=1

(We have written M1 (x0 , dx1 ) for M1 (dx1 ).) Equivalently,
π1:t (φ) ,

E [φ(ξ1:t )g1:t (ξ1:t )]
,
Zt

φ : E t → R measurable,

(1.2)

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

3

where Zt , E [g1:t (ξ1:t )] is the normalization constant.1
Towards the goal of efficiently approximating π1:t , we introduce a novel sequential Monte Carlo algorithm: the i-times conditional adaptive resampling particle
filter (ci ARPF), which is a generalization of the adaptive resampling particle filter [Gordon et al., 1993, Künsch, 2013] and the conditional SMC algorithm used
in particle Gibbs [Andrieu et al., 2010]. (In Section 3, we will introduce a further
generalization.) The integer parameter i ≥ 0 determines the number of fixed tra1
i
jectories y1:t
, . . . , y1:t
∈ E t required by the algorithm, which operates by generating
1:N
n N
a collection X̃ 1:t , {X̃1:t
}n=1 of N > i particles with corresponding nonnega1:N
tive weights Wt
, {Wtn }N
n=1 . When i = 0, we recover the standard (unconditional) adaptive resampling particle filter; when i = 1, we recover a generalization
of the conditional SMC algorithm that includes adaptive resampling. For time
s = 1, . . . , t, the measure π1:s is approximated by
i,N
π1:s
,

N
X

Wsn gs (Xsn )
δ n .
PN
k
k X̃1:s
k=1 Ws gs (Xs )
n=1

1:N
The ci ARPF algorithm iteratively constructs X̃ 1:N
as follows: The
1:t and Wt
first i particles are deterministically set to match the fixed trajectories:

Xsn = ysn ,

n
n
X̃1:s
= y1:s
,

s = 1, . . . , t and n = 1, . . . , i.

At time s = 1, the remaining N − i particles X1n , for n = i + 1, . . . , N , are sampled
independently and identically from M1 . The corresponding (length 1) trajectories
are
n
X̃1:1
= X1n ,

n = i + 1, . . . , N.

Furthermore, for all n = 1, . . . , N , W1n = 1.
The remaining particle trajectories are generated as follows: First, we introduce a cutoff parameter η ∈ [0, 1] and an effective sample size (ESS) function
ESS : RN
+ → [1, N ]. The ESS function measures how uniform the current weights
1:N
Ws are. Typically ESS(Ws1:N ) = 1 indicates that all but one weight is zero and
ESS(Ws1:N ) = N indicates all the weights are equal.
For each time s = 2, . . . , t:
1:N
• If ESS(Ws−1
) ≤ ηN , a resampling step is introduced. For n = 1, . . . , N ,
the weights are set to a common value

Wsn = Ws ,

N
1 X k
k
Ws−1 gs−1 (Xs−1
)
N
k=1

and, for n = i + 1, . . . , N , particle n’s “ancestor” at time s, denoted Ans , is
sampled independently, such that Ans = k, for k = 1, . . . , N , with probability
k
k
Ws−1
gs−1 (Xs−1
)
.
N Ws
1In the state-space setting, the potential g would be the conditional density (i.e., likelihood)
t
of the observation vt at time t as a function of unobserved state xt : i.e., gt (xt ) = pt (vt | xt ).
Then π1:t would be the posterior distribution of the unobserved state sequence given the observed
sequence.

4

J. H. HUGGINS AND D. M. ROY
1:N
• If ESS(Ws−1
) > ηN , then the algorithm does not resample the particles.
For n = 1, . . . , N , the weights are copied, i.e.,
n
n
Wsn = Ws−1
gs−1 (Xs−1
),

and, for n = i + 1, . . . , N , a record is made that particle n was its own
ancestor by setting Ans = n.
• Having sampled the ancestors, the algorithm propagates the particles forAn
s
ward. For n = i + 1, . . . , N , Xsn is sampled from Ms (Xs−1
, ·), and the
corresponding trajectories are set to
An

n
s
X̃1:s
= hX1:s−1
, Xsn i .
∗
In the final step of the algorithm, a single particle X̃1:t
is sampled from the full
i,N
approximation π1:t , and the algorithm yields an estimate of the normalization
constant Zt ,

Ẑt ,

N
1 X n
W gt (Xtn ).
N n=1 t

Let Ei,N
[·] denote the expectation operator with respect to the ci ARPF, and
y 1:i
1:t
write
i,N
i,N
P i,N (y 1:i
1:t , dx1:t ) , Py 1:i (dx1:t ) , Ey 1:i [δX̃ ∗ (dx1:t )]
1:t

∗
X̃1:t

1:t

1:t

y 1:i
1:t

for the law of
when the i fixed trajectories are
∈ (E t )i . We can now
i
describe in more precise terms how the c ARPF kernel P i,N generalizes several
0,N
well-known SMC kernels. When i = 0, π1:t
is the standard adaptive SMC particle
∗
approximation of π1:t and X̃1:t is a single sample from the particle approximation. When i = 1 and resampling is done at every step by taking η = 1, P 1,N
is exactly the conditional SMC kernel used in particle Gibbs samplers [Holenstein,
2009, Andrieu et al., 2010, 2013]. For general η ∈ (0, 1), we obtain a novel adaptive
resampling variant that we study in the sequel. In particular, under mild regularity
conditions, P 1,N defines a Markov kernel with invariant distribution π1:t .
1.2. Controlling ci ARPF efficiency with ∞-ESS. We can analyze the quality
of the ci ARPF kernel P i,N by quantifying the extent to which high-probability
sets under the target distribution also have high probability under the kernel. The
following theorem establishes a minorization condition for the i-times conditional
filter in terms of the (i + 1)-times conditional filter:
1
i
Theorem 1.1. For all t ≥ 1, i ≥ 0, N > i, and y1:t
, . . . , y1:t
∈ Et,
Z
Zt
i+1
t
π1:t (dy1:t
), S ⊆ E t measurable. (1.3)
P i,N (y 1:i
1:t , S) ≥ (1 − i/N )
i+1,N
S E 1:i+1 [Ẑt ]
y
1:t

The integral appearing in Eq. (1.3) has no simple form in general, but in many
settings we will be able to obtain a lower bound on the integrand that does not depend on the fixed trajectories y 1:i+1
1:t . In those cases, the integral is simply replaced
by this uniform lower bound. For i = 0, we are then immediately able to control
∗
numerous measures of divergence between π1:t and P 0,N (i.e., the law of X̃1:t
). For
example, in the case of total variation distance, we have the following corollary to
Theorem 1.1:

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

5

t
Corollary 1.2. If E1,N
y1:t [Ẑt /Zt ] ≤ Bt,N for all y1:t ∈ E , then
−1
dT V (π1:t , P 0,N ) ≤ 1 − Bt,N
.

(1.4)

For i = 1, a uniform lower bound assumption implies a minorization condition
on the kernel P 1,N (y1:t , dx1:t ), which in turn implies fast mixing of the Markov
chain with kernel P 1,N :
1
2
Corollary 1.3. If E2,N
[Ẑt /Zt ] ≤ Bt,N for all y1:t
, y1:t
∈ E t , then the Markov
y 1:2
1:t

chain with transition kernel P 1,N (y1:t , dx1:t ) is uniformly ergodic in total variation
distance and has invariant distribution π1:t . In particular, for all y1:t ∈ E t and
k ≥ 1,
k

−1
(1.5)
(1 − 1/N )t ,
dT V (π̃yk1:t , π1:t ) ≤ 1 − Bt,N
where π̃yk1:t , δy1:t [P 1,N ]k is the law of the Markov chain, with initial state y1:t ,
after k transitions.

In order to apply the corollaries, it remains to bound Ei,N
[Ẑt /Zt ]. Such a bound
y 1:i
1:t
was obtained for the nonadaptive conditional SMC kernel in Andrieu et al. [2013].
However, in our general adaptive resampling setting, one must make a careful choice
of effective sample size function. To this end, we introduce a generalized notion
of effective sample size, which includes several existing definitions as special cases.
p
For p ∈ (1, ∞], let p∗ , p−1
be the conjugate exponent of p (so 1/p + 1/p∗ = 1).
The p-effective sample size (p-ESS) of the weight vector w1:N ∈ RN
+ is
kw1:N kp1∗
.
kw1:N kpp∗
The following proposition highlights some elementary properties of p-ESS.
ESSp (w1:N ) ,

Proposition 1.4. The p-ESS has the following properties:
(1) For all p ∈ (1, ∞], 1 ≤ ESSp (w1:N ) ≤ N . The lower bound is achieved if and
only if all but one of the weights is zero. The upper bound is achieved if and
only if all the weights are equal.
(2) For 1 < p < q ≤ ∞, ESSp (w1:N ) ≥ ESSq (w1:N ) ≥ N −(1−q∗ /p∗ ) ESSp (w1:N ),
with equality if and only if K weights are equal and the rest are zero.
Part (1) demonstrates that the p-ESS satisfies basic properties one would expect
of a measure of effective sample size. Part (2) places the family of p-ESS measures
in a linear order: the larger the value of p, the more stringent the notion of effective
sample size.
The standard definition of effective sample size is precisely the 2-ESS. Whiteley et al.
[2016] provided a rigorous justification for the use of 2-ESS from the operator perspective: if adaptive resampling is used to guarantee that the 2-ESS does not fall
below ζN , for some fixed parameter ζ ∈ (0, 1], then the error bounds on the operator approximation match those of the nonadaptive sampler with ζN particles.
More formally, let ℑt (y1:t ) = yt be the projection onto the t-th component. Under
appropriate regularity conditions, for every bounded measurable φ : E → R and
real r ≥ 1,
i1/r
h
a(r)b(φ)
0,N
(φ ◦ ℑt ) − π(φ ◦ ℑt )|r
≤ √
sup ESS2 (Wt1:N ) ≥ ζN =⇒ sup E0,N |π1:t
,
ζN
t≥1
t≥1

6

J. H. HUGGINS AND D. M. ROY

where a(r) and b(φ) are explicit functions.
To upper bound Ei,N
[Ẑt /Zt ], however, we will require a lower bound on the
y 1:i
1:t
∞-ESS, which by Proposition 1.4(2) is a more stringent notion of effective sample
size than 2-ESS:
Assumption 1.A. There exists ζ ∈ (0, 1] such that ESS∞ (Ws1:N ) ≥ ζN for all
1 ≤ s ≤ t.
For the choice ESS = ESS∞ (i.e., under Assumption 1.A), we can bound the
estimate of the normalization constant. Let Gs,t (xs ) , E[gs:t−1 (ξs:t−1 ) | ξs = xs ]
for s = 1, . . . , t and G0,t , E[g1:t−1 (ξ1:t−1 )]. We now arrive at our second main
result:
1
i
Theorem 1.5. If Assumption 1.A holds, then for all i, t ≥ 1, N > i, y1:t
, . . . , y1:t
∈
t
E ,
Pt Pi
i
Zt−1 s=1 j=1 G0,s Gs,t+1 (y1:t
) − ζi
i,N
+ Θ(N −2 ).
(1.6)
Ey 1:i [Ẑt /Zt ] ≤ 1 +
1:t
ζN

Two possible further assumptions both lead to uniform bounds on Ei,N
[Ẑt ].
y 1:i
1:t

Assumption 1.B. The potentials satisfy gs , supx∈E gs (x) < ∞ for all 1 ≤ s ≤ t.
Assumption 1.C. There exists a constant β > 0 such that for any t, s ∈ N,
sup
x∈E

G0,t Gt,t+s (x)
≤ β.
G0,t+s

Corollary 1.6. Under the same conditions as Theorem 1.5, if Assumption 1.B
holds then
#
"
t
t
Y
i
i,N
−1
Ey 1:i [Ẑt /Zt ] ≤ 1 + Zt
−1
(1.7)
1+
gs
1:t
ζN
s=1
while if Assumption 1.C holds then
Ei,N
[Ẑt /Zt ]
y 1:i
1:t


t
β
≤ 1+
.
ζN

(1.8)

Combining Corollaries 1.2 and 1.6 yields the following guarantees for the ARPF
sampler:
Theorem 1.7. If Assumptions 1.A and 1.B hold then
Qt
tZt−1 s=1 g s
+ Θ(N −2 )
dT V (π1:t , P 0,N ) ≤
Q
ζN + tZt−1 ts=1 gs

(1.9)

while if Assumptions 1.A and 1.C hold then
dT V (π1:t , P 0,N ) ≤

βt
+ Θ(N −2 ).
ζN + βt

(1.10)

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

7

1.3. Applications to Particle Gibbs. In the language of state-space models,
the setting described so far involves approximating the posterior distribution of a
Markov chain given indirect stochastic observations of the chain’s values. However,
it is often the case that the chain and the potentials are controlled by a global
parameter θ ∈ Θ for which there is a prior distribution ̟(dθ). Replace Ms by
Msθ and gs by gsθ , then parameterize the other quantities defined previously in
terms of Ms and gs by θ. Let (Y, Y) , (E t , B(E t )). We will suppress much of the
time dependence when possible to make the notation less cluttered. The target
distribution on the product space (Θ × Y, B(Θ × Y )) is
π(dθ × dy) , γ(dθ × dy)/Z,

(1.11)

where
γ(dθ × dy) , ̟(dθ)

t
Y

gsθ (ys )Msθ (ys−1 , dys ) and Z , γ(1).

(1.12)

s=1

Let πθ (dy) and πy (dθ) denote the disintegrations of π along Θ and along Y , respectively.
The particle Gibbs sampler approximates the two-stage Gibbs kernel
Π(θ, y, dϑ × dz) , πy (dϑ)πϑ (dz).

(1.13)

In many settings, such as non-linear or non-Gaussian state-space models, it is possible to sample from πy (dϑ), but difficult to sample from πϑ (dz). The idea is to
replace πϑ (dz) with an SMC-based approximation Πϑ (y, dz) that leaves πϑ (dz)
invariant, leading to a kernel of the form πy (dϑ)Πϑ (y, dz).
We introduce the adaptive resampling particle Gibbs (ARPG) sampler, which
1,N
employs the cARPF kernel Pϑ,y
(dz) to approximate the conditional distribution
πϑ (dz) that would be used in a standard Gibbs sampler. The ARPG kernel is thus
given by
1,N
ΠN (θ, y, dϑ × dz) , πy (dϑ)Pϑ,y
(dz).

(1.14)

Theorem 1.1 and Corollaries 1.3 and 1.6 together with the results of Andrieu et al.
[2013] yield guarantees on the ergodicity properties of the cARPF kernel and the
ARPG sampler. Once instances of N are replaced by ζN , the guarantees essentially
match those provided by Andrieu et al. [2013] for the standard PG sampler.
Theorem 1.8. If Assumption 1.A holds, then the cARPF kernel and ARPG sampler have the following properties:
(1) If Assumption 1.B holds then there exists ǫt,N = 1 − Ct /N such that for any
θ ∈ Θ, y ∈ Y , and k ≥ 1,
dT V (δy [Pθ1,N ]k , πθ ) ≤ (1 − ǫt,N )k .

(1.15)

(2) If Assumption 1.C holds and N ≥ t/C + 1 for any fixed C > 0, then for any
t ≥ 1, Eq. (1.15) holds with


C
(1.16)
ǫt,N ≥ exp − (2β + ζ) .
ζ
(3) If either Assumption 1.B or Assumption 1.C holds, then whenever the Gibbs
sampler is geometrically ergodic the ARPG sampler is geometrically ergodic as
well.

8

J. H. HUGGINS AND D. M. ROY

At a high level, the results we have obtained highlight the role of the expected
value of Ẑt in the mixing properties of conditional SMC Markov chains and particle Gibbs (PG) samplers: In order to show geometric ergodicity for adaptive
resampling particle Gibbs samplers, it suffices to establish bounds on the expected
value of Ẑt under the twice-conditional filter, and the growth of the expectation as
t increases determining how well the particle Gibbs algorithm scales. Similarly, a
bound on the expected value of Ẑt under the once-conditional filter implies a bound
on dT V (π1:t , P 0,N ). Hence, as a slogan, good performance of (adaptive resampling)
particle Gibbs is equivalent to good performance of (adaptive) SMC for sampling.
2. Preliminaries
In this section, we fix some additional notation, introduce a few key additional
definitions, and then present αSMC [Whiteley et al., 2016], a generalization of the
adaptive resampling particle filter described in the introduction.
For a positive integer K, let [K] , {1, 2, . . . , K}. If xi , . . . , xj are elements of a
P
sequence,
write xi:j , hxi , xi+1 , . . . , xj i. We use the following conventions: ∅ = 0,
Q
∅ = 1, and 0/0 = 0.
Let (S, S), (S ′ , S ′ ) be measurable spaces. Then K : S × S ′ → R is a kernel if
K(·, B) is a (S, S)-measurable function for all B ∈ S ′ and K(x, ·) is measure on
(S ′ , S ′ ) for Rall x ∈ S. For a measure µ on (S, S)Rand kernels K, K ′ : S × S → R, let
µK(dy) , µ(dx)K(x, dy) and KK ′ (x, dz) , K(x, dy)K ′ (y, dz). We will often
use measures and kernels
function φ : S → R, let
R as operators. For a measurable
R
µ(φ) , Eξ∼µ [φ(ξ)] = φ(x)µ(dx) and K(φ)(x) , φ(y)K(x, dy). For measures
µ, ν on (S, S), we will write µ ≪ ν to denote that µ is absolutely continuous with
respect to ν, in which case we will write
R dµ/dν for the ν-almost everywhere (ν-a.e.)
unique function f satisfying µ(A) = A f dν, for all A ∈ S. When the choice is
clear from context, we may write B(S) for the σ-algebra of the space S.
For probability measures µ and ν on (S, S), the total variation distance between
µ and ν is
dT V (µ, ν) , sup |µ(A) − ν(A)|.

(2.1)

A∈S

If µ ≪ ν, then the Kullback–Liebler (KL) divergence is
dKL (µ||ν) , µ(log dµ/dν)

(2.2)

dχ2 (µ||ν) , ν([dµ/dν − 1]2 ) = µ(dµ/dν) − 1.

(2.3)

and the χ2 divergence is

Finally, we note that, when there is little risk of confusion, we will ignore
measure-theoretic niceties such as the distinction between equality and a.e.-equality.
Recall that (ξt )t≥1 is an inhomogeneous Markov chain on (E, E) with transition
kernels (Mt )t≥2 and initial distribution M1 , and that E[·] denotes expectation with
respect to the Markov chain. We will write M1 (x0 , ·) for M1 (·) when convenient
and, for all t ≥ 1 and xt−1 ∈ E, we will assume that Mt (xt−1 , ·) has a density with
respect to some common σ-finite dominating measure (which we denote by dx).
We will abuse notation and write Mt (xt−1 , xt ) for the density of Mt (xt−1 , ·) as xt .
Recall that, for each t ≥ 1, gt : E → R+ denotes a E-measurable potential function,
Qt
with g0 ≡ 1. Finally, recall that gs:t (xs:t ) , τ =s gτ (xτ ).

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

9

2.1. Target distributions. We now introduce some additional target distributions. (We will also repeat the definition of γ1:t and π1:t for completeness.)
Let φ1:t : E t → R and φt : E → R denote generic measurable functions. For each
t ≥ 1, the unnormalized predictive and updated measures are defined, respectively,
by
′
γ1:t
(φ1:t ) , E [φ1:t (ξ1:t )g1:t−1 (ξ1:t−1 )]

and

γ1:t (φ1:t ) , E [φ1:t (ξ1:t )g1:t (ξ1:t )]

and

γt (φt ) , E [φt (ξt )g1:t (ξ1:t )] .

with corresponding marginal versions
γt′ (φt ) , E [φt (ξt )g1:t−1 (ξ1:t−1 )]

Our ultimate goal is to approximate the normalized predictive and updated measures along with their marginal versions:
γ1:t (φ1:t )
,
Zt
γt (φt )
,
πt (φt ) ,
Zt

′
γ1:t
(φ1:t )
,
Zt′
γ ′ (φt )
ηt (φt ) , t ′ ,
Zt

η1:t (φ1:t ) ,

π1:t (φ1:t ) ,

and

(2.4)
(2.5)

where Zt , γt (1) and Zt′ , γt′ (1) are normalization constants.
2.2. The αSMC algorithm. In the introduction, adaptation in the particle filter was implemented via a simple multinomial resampling step, triggered when the
effective sample size fell below a fixed threshold. For the remainder of the article,
we will consider a more general mechanism for adaptation captured by the αSMC
algorithm introduced by Whiteley et al. [2016]. The αSMC algorithm can produce sequential importance sampling (SIS), sampling importance resampling (SIR,
also known as the bootstrap filter), and numerous other SMC variants as special
cases. Not only does the αSMC formulation aid in analyzing adaptive resampling
strategies, it provides a useful framework for devising novel adaptive schemes with
attractive computational properties, such as admitting parallelization even on resampling steps. In the remainder of this section, we outline the (unconditional)
αSMC algorithm. In the following section, we introduce a novel i-times conditional
version of αSMC, which will include the ci ARPF as a special case.
The αSMC algorithm, which is given as Algorithm 1, provides a flexible resampling mechanism: at each time t, a stochastic matrix αt−1 is chosen from a set AN
of N × N matrices. We denote the value in the n-th row and k-th column of αt−1
by αnk
t−1 . The αSMC estimators are
0,N
π1:t
,

N
X

Wtn gt (Xtn )
δ n ,
PN
k
k X̃1:t
n=1
k=1 Wt gt (Xt )

0,N
η1:t
,

N
X

n=1

Wtn

PN

k=1

Wtk

δX̃ n ,
1:t

πt0,N ,

N
X

Wtn gt (Xtn )
δXtn ,
PN
k g (X k )
W
t
t
t
n=1
k=1

and ηt0,N ,

N
X

Wtn
δXtn ,
PN
k
n=1
k=1 Wt

and the estimators of the normalization constants Zt and Zt′ are
Ẑt ,

N
1 X n
W gt (Xtn )
N n=1 t

and Ẑt′ ,

N
1 X n
W .
N n=1 t

Expectations with respect the law of the αSMC algorithm are written as E0,N [·].

10

J. H. HUGGINS AND D. M. ROY

Algorithm 1 αSMC
for n = 1, . . . , N do
Sample X1n from M1
n
Set X̃1:1
← X1n
n
Set W1 ← 1
end for
for t = 2, 3, . . . do
1:N
Select αt−1 from AN according to some function of X 1:t−1
and A1:N
1:t−2
for n = 1, . . . , N do
P
N
k
k
Set Wtn ← k=1 αnk
t−1 Wt−1
Dgt−1 (Xt−1 )
EN 
k
k
αnk
t−1 Wt−1 gt−1 (Xt−1 )
n
Sample At−1 from Multi
Wn
t

k=1

n

At−1
Sample Xtn Dfrom Mt (Xt−1
E , ·)
n
At−1
n
n
Set X̃1:t ← X1:t−1 , Xt

end for
end for

SIS, SIR, and the standard adaptive algorithm are obtained as special cases of
αSMC as follows. SIS is recovered when αt−1 = I N , the N × N identity matrix,
while SIR is recovered when αt−1 = 11/N , the N × N matrix with all entries equal
to 1/N . The adaptive particle filter (APF) algorithm is obtained by setting αt−1
1:N
to 11/N if ESS2 (Wt−1
) ≤ ζN and to I N otherwise, where ζ ∈ (0, 1] is fixed.
3. Conditional αSMC
It is useful both algorithmically and analytically to generalize αSMC in such a
way that one or more particle trajectories is fixed ahead of time. The result, which
we will refer to as conditional αSMC, is a strict generalization of the conditional
adaptive particle filter given in the introduction. We will use conditional αSMC
to study the properties of (unconditional) αSMC, to design novel adaptive particle
Gibbs algorithms, and to analyze their mixing properties.
For this section, fix t ≥ 1, i ≥ 0, and N > i. The i-times conditional αSMC
(ci αSMC) process (or simply the cαSMC process when i = 1) is defined on the
space (E N × [N ]N × [N ]i )t−1 × E N × [N ] × [N ]i , and is essentially equivalent to
αSMC with the first i particle trajectories, but not their lineages,
Q fixed a priori.
′
If f 1:i ∈ [N ]i are indices of the first i particles, let D(f 1:i ) , j6=j ′ 1(f j 6= f j )
be the function that indicates whether the indices are distinct. As in αSMC, the
1:N
1:N
N t
matrix αt−1 ∈ AN is a function of X 1:N
1:t−1 and A1:t−2 . We have x1:t ∈ (E ) ,
1:i
i t
1:N
N t−1
∗
f 1:t ∈ ([N ] ) , a1:t−1 ∈ ([N ] ) , and at ∈ [N ], and use the notation
w1n , 1,

wtn ,

N
X

k
k
αnk
t−1 wt−1 gt−1 (xt−1 ),

(3.1)

k=1

and
1:N
rn (k|ws−1
, x1:N
1:s−1 ) ,

k
k
αnk
s−1 ws−1 gs−1 (xs−1 )
.
wsn

(3.2)

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

11

Algorithm 2 Conditional αSMC
Require: Fixed trajectory y1:t
Sample F11 uniformly from {1, . . . , N }
for n = 1, . . . , N do
if n = F11 then
Set X1n ← y1
else
Sample X1n from M1
end if
n
Set X̃1:1
← X1n
n
Set W1 ← 1
end for
for t = 2, 3, . . . do
1:N
1:N
Select αt−1 from AN according to some
! function of X 1:t−1 and A1:t−2


N
kF 1
Sample Ft1 from Multi
αt−1t−1
k=1

for n = 1, . . . , N do
PN
k
k
Set Wtn ← k=1 αnk
t−1 Wt−1 gt−1 (Xt−1 )
if n = Ft1 then
1
Set Ant−1 ← Ft−1
Set Xtn ← yt
else
D
EN 
k
k
αnk
t−1 Wt−1 gt−1 (Xt−1 )
Sample Ant−1 from Multi
Wn
t

k=1

An
t−1
Mt (Xt−1
, ·)

Xtn

Sample
from
end if
D An
E
t−1
n
Set X̃1:t
← X1:t−1
, Xtn
end for
end for
D
EN 
W k g (X k )
PN t tn t n
Sample A∗t from Multi
W g (X )
n=1

t

t

t

k=1

1
i
For fixed trajectories y1:t
, . . . , y1:t
∈ E t , the law of the ci αSMC process is given by
1:i
1:i
i
1:i
Pi,N
[X11:N ∈ dx1:N
1 , F1 = f1 ] , C1 D(f1 )
y1:i
1:t

N
i
Y
Y
1
fj
M1 (dxn1 ), (3.3)
δyj (dx11 )
1
N
1:i
j=1
n∈f
/ 1

for s = 2, . . . , t,
1:N
1:N
1:i
1:i
Pi,N
[Xs1:N ∈ dx1:N
s , As−1 = as−1 , Fs = fs |
y 1:i
1:t

1:N
1:N
1:N
1:i
1:i
X 1:N
1:s−1 = x1:s−1 , A1:s−2 = a1:s−2 , Fs−1 = fs−1 ]

, Csi D(fs1:i )
×

Y

n∈f
/ s1:i

i
Y

fjfj

fj

fj

j
s s−1
s
= fs−1
)
δysj (dxss )1(as−1
αs−1

j=1
an

s−1
1:N
n
rn (ans−1 |ws−1
, x1:N
1:s−1 )Ms (xs−1 , xs )

(3.4)

12

J. H. HUGGINS AND D. M. ROY

and
a∗

Pi,N
[A∗t
y1:i
1:t

=

a∗t

| X 1:N
1:t

=

1:N
x1:N
1:t , A1:t−1

=

a1:N
1:t−1 ]

a∗

w t gt (xt t )
.
, PN t
n
n
n=1 wt gt (xt )

(3.5)

The Csi terms are normalization constants that ensure the expressions are valid
A∗
∗
probabilities. Let X̃1:t
, X1:tt , let Ei,N
[·] denote the expectation operator with
y 1:i
respect to the ci αSMC, and write

1:t

i,N
i,N
P i,N (y 1:i
1:t , dx1:t ) , Py 1:i (dx1:t ) , Ey 1:i [δX̃ ∗ (dx1:t )]
1:t

1:t

1:t

∗
X̃1:t
.

for the law of
The normalization constants Csi arise because the lineages f 1:i
1:t of the fixed trai
jectories y 1:i
are
kept
distinct.
The
c
αSMC
kernel
enforces
distinct
lineages for
1:t
j
j′
′
the fixed trajectories since in general y1:t 6= y1:t for j 6= j and, from an algorithmic
standpoint, allowing overlapping lineages could lead to a substantial increase in
complexity, both in terms of implementation and computation. The distinct lineage requirement is enforced by the D(fs1:i ) terms. Since there is at most one fixed
trajectory when i = 0 or 1, Cs0 = Cs1 = 1 for all s ∈ [t].
Algorithm 2 provides pseudocode to sample from the law of the cαSMC process,
which is a necessary part of the particle Gibbs sampler described in Section 6.
Sampling from the law of the ci αSMC process for i > 1 is more delicate, but
unnecessary since these are only used for analytical purposes.
Remark 3.1. To recover the ci ARPF described in the introduction, let αt−1 = 11/N
1:N
if ESS(Wt−1
) ≤ ζN and let αt−1 = I N otherwise. Then note that by the symmetry
of 11/N and I N , instead of sampling them, we can set the lineage for the j-th fixed
trajectory to j: that is, set fsj = j for all j ∈ [i] and s ∈ [t].
▽
4. Main Results
We are now ready to undertake our study of the i-times conditional αSMC kernel
1:i
, dx1:t ). Specifically, our aim is to understand the conditions under which
P i,N (y 1:t
i
the c αSMC kernel is close to π1:t (dx1:t ). Formally, we will establish a minorization
condition for the i-times conditional filter in terms of the expected value of Ẑt under
the (i + 1)-times conditional filter. The remainder of the section presents conditions
under which the expected value of Ẑt can be bounded. One of the key assumptions
is that adaptation controls the ∞-ESS.
Of particular interest are the cases i = 0, which corresponds to the αSMC filter,
and i = 1, which corresponds to the conditional αSMC kernel. The former provides
approximate samples from π1:t . The latter can be used to define a Markov chain
with invariant distribution π1:t , producing an adaptive resampling particle Gibbs
sampler. We consider both these applications in, respectively, Sections 5 and 6.
4.1. A minorization condition for the ci αSMC kernel. For the remainder of
the article, we will work under the following assumption:
Assumption 4.D. For all N ≥ 1, all α ∈ AN are doubly stochastic.

Remark 4.1. Assumption 4.D is the same as Assumption (B++ ) in [Whiteley et al.,
2016], although there, the condition is stated as assuming each α admits the uniform
distribution as an invariant measure.
▽

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

13

Let
κN ,

max

n6=m,α∈AN

N
X

αkn αkm

k=1

and

κ′N , max{κN , 1/N }.

(4.1)

Our first main result provides control over how much the measure Pyi,N
1:i differs from
1:t
π1:t . The theorem gives a stronger result when i = 0 and gives a simpler result
when i = 1, by expressing the lower bound on P 1,N (y1:t , S) in terms the more
transparent quantity κ′N instead of the normalization terms Cs2 . For example, if
AN = {I N } then κN = 0, while if AN = {11/N }, then κN = 1/N , so in either case
κ′N = 1/N .
Theorem 4.2. If Assumption 4.D holds, then for all t ≥ 1, i ≥ 0, N > i, S ⊆ E t
1
i
measurable, and y1:t
, . . . , y1:t
, x1:t ∈ E t ,
Z
Zt
P i,N (y 1:i
,
S)
≥
π1:t (dx1:t ).
(4.2)
Qt
1:t
i+1,N
i+1
i
C
[
Ẑ
E
S y 1:i ,x
s /Cs ]
t
s=1
1:t
1:t

In particular, in the case of i = 0, we have
 
Zt
Zt
dP 0,N
≥ 1,N
(x1:t ) = E1,N
x1:t
dπ1:t
Ẑt
Ex1:t [Ẑt ]

(4.3)

and in the case of i = 1, we have
P

1,N

(y1:t , S) ≥

Z

S

Zt (1 − κ′N )t

E2,N
y1:t ,x1:t [Ẑt ]

π1:t (dx1:t ).

(4.4)

Remark 4.3. By identical arguments, Theorem 4.2 also holds in the marginal and
predictive cases. In the predictive cases, however, Ẑt′ and Zt′ replace, respectively, Ẑt
and Zt . In the predictive case, under Assumption 4.D, Ẑt′ = Ẑt−1 and Zt′ = Zt−1 , so
later results pertaining to Ei,N
[Ẑt ], such as Proposition 4.4 and Theorem 4.7, apply
y1:i
1:t
′
to Zt as well. The fact that Zt′ = Zt−1 follows immediately from the definitions.
To show that Ẑt′ = Ẑt−1 , apply Assumption 4.D:
X
XX
X
k
k
−1
k
k
Ẑt′ = N −1
Wtn = N −1
αnk
Wt−1
gt−1
= Ẑt−1 .
▽
t−1 Wt−1 gt−1 = N
n

n

k

k

We will prove Theorem 4.2 in two parts: first for the case of i = 0, then for
i ≥ 1. For the i = 0 case (corresponding to vanilla αSMC), we begin by writing
the joint density of the αSMC process as
1:N
ψ(x1:N
1:t , a1:t−1 )

,

N
Y

n=1

!

M1 (xn1 )

t Y
N
Y

s=2 n=1

!

an
s−1
1:N
n
, x1:N
rn (ans−1 |ws−1
1:s−1 )Ms (xs−1 , xs )

(4.5)

.

1:N
1:N
1
Under Assumption 4.D, the density of the cαSMC process with law P1,N
y1:t [X 1:t , A1:t−1 , F1:t ]
f1

t
can be written in the following “collapsed” form, by implicitly identifying x1:t
with

14

J. H. HUGGINS AND D. M. ROY

y1:t :
1:N
1
ψ̃(x1:N
1:t , a1:t−1 , f1:t )
1
Qt
fs1 fs−1
1:N
ψ(x1:N
1:t , a1:t−1 )
s=2 Is αs−1
=
(4.6)
1
1
f1 Q
1 |w1:N , x1:N )M (xfs−1 , xfs )
N M1 (x11 ) ts=2 rfs1 (fs−1
s
s
1:s−1
s−1
s−1


t
Y
Y
n
1 1
1 Y
a
f
f
s−1
s s−1
1:N
n 
Is αs−1
rn (ans−1 |ws−1
, x1:N
M1 (xn1 )
=
1:s−1 )Ms (xs−1 , xs ) ,
N
1
1
s=2

n6=fs

n6=f1

fs1

1
where Is , 1(as−1 = fs−1
).

Proof of Theorem 4.2, i = 0 case. Consider the density
f1

1:N
1
1:N
1:N
1
t
π̃1:t (x1:N
1:t , a1:t−1 , f1:t ) , π1:t (x1:t )ψ̃(x1:t , a1:t−1 , f1:t ).

Then
f1

f1

1:N
t
t
ψ(x1:N
1:t , a1:t−1 )gt (xt )wt
P
N
1:N
1
n
n
π̃1:t (x1:N
1:t , a1:t−1 , f1:t )
n=1 gt (xt )wt
1
fs−1
f 1 Qt
fs1
ft1
ft1
1
1:N
M1 (x11 ) s=2 rfs1 (fs−1
|ws−1
, x1:N
1:s−1 )Ms (xs−1 , xs )gt (xt )wt
=
1
PN
fs1 fs−1
ft1 Qt
π1:t (x1:t
) s=2 Is αs−1
N −1 n=1 gt (xnt )wtn
1
1
1
1
fs1 fs−1
fs−1
fs−1
fs−1
f 1 Qt
f1
f1
f1
M1 (x11 ) s=2 αs−1
ws−1
gs−1 (xs−1
)Ms (xs−1
, xss )gt (xt t )wt t
=
1
PN
fs1 fs−1
ft1 Qt
f 1 Qt
π1:t (x1:t
) s=2 ws s s=2 Is αs−1
N −1 n=1 gt (xnt )wtn
1
Qt
fs−1
fs1
fs1
s=1 gs (xs )Ms (xs−1 , xs )
=
Qt
P
ft1
n
n
)N −1 N
π1:t (x1:t
s=2 Is
n=1 gt (xt )wt
1
Zt
,
=
Qt
Ẑt s=2 Is

Using the convention that 0/0 = 0, it follows that

P 0,N (dx1:t )
(
)
a∗
a∗
X Z
gt (xt t )wt t
1:N
1:N
ψ(x1:t , a1:t−1 ) PN
=
δ at (dx1:t ) dx1:N
1:t
n )wn x1:t
g
(x
t
∗
t
t
n=1
a1:N
,a
1:t−1 t
(
ft1
ft1 Qt
1:N
X Z
ψ(x1:N
1:t , a1:t−1 )gt (xt )wt
s=2 Is
=
1:i PN
1:N
1:N
π̃1:t (x1:t , a1:t−1 , f 1:t ) n=1 gt (xnt )wtn
1
a1:N
1:t−1 ,f1:t
)
1:N
1
1:N
× π̃1:t (x1:N
1:t , a1:t−1 , f1:t )δxft (dx1:t ) dx1:t
1:t

X

=

1
a1:N
1:t−1 ,f1:t

=

(

X

Z (

Zt

Z

Zt

1
a1:N
1:t−1 ,f1:t

Ẑt

)

1:N
1
π̃1:t (x1:N
t (dx1:t )
1:t , a1:t−1 , f1:t )δxf1:t

Ẑt

dx1:N
1:t

−ft1
ft
1:N
1
ψ̃(x1:N
1:t , a1:t−1 , f1:t )δx1:t (dx1:t )dx1:t

)

π1:t (dx1:t ).

(4.7)

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

15



The result follows from Lemma A.1.
We defer the proof of Theorem 4.2 in the i ≥ 1 case to Appendix A.2.

4.2. Bounding Ẑt under the ci αSMC kernel. In order to apply Theorem 4.2,
we must be able to control the quantity Ei,N
[Ẑt ]. As an initial step toward this
y 1:i
1:t
goal, we consider the SIR case:
Assumption SIR. For all s ∈ [t − 1], αs = 11/N .
For SIR, we can rewrite Ei,N
[Ẑt ] in an equivalent but more explicit form (Propoy 1:i
1:t

sition 4.4). Our goal will then be to provide general conditions under which Ei,N
[Ẑt ]
y 1:i
1:t
can be rewritten in a similar manner (Theorem 4.7).
Recall that G0,t , E[g1:t−1 (ξ1:t−1 )] and Gs,t (xs ) , E[gs:t−1 (ξs:t−1 ) | ξs = xs ] for
s ∈ [t] and xs ∈ E. For t ≥ 1, 1 ≤ ℓ ≤ s ≤ t + 1, let
Tt,ℓ,s , {hτ1 , . . . , τℓ i : t − s + 1 < τ1 < · · · < τℓ = t + 1}

and, for τ ∈ Tt,ℓ,s , define

Cℓ (τ , y1:t ) ,

ℓ−1
Y

Gτi ,τi+1 (yτi ).

(4.8)

(4.9)

i=1

We will sometimes write Cℓy (τ ) or Cℓτ (y1:t ) instead of Cℓ (τ , y1:t ). The following is
a straightforward generalization of [Andrieu et al., 2013, Proposition 9].
Proposition 4.4. If Assumption SIR holds, then for all t ≥ 1, i ≥ 1, N ≥ i,
i
1
∈ Et,
, . . . , y1:t
y1:t
Ei,N
[Ẑt ] =
y 1:i
1:t

t+1
1 X
(N − i)t+1−ℓ
Nt
ℓ=1

X

G0,τ1

τ ∈Tt,ℓ,t+1

In particular, in the case of i = 1, we have
E1,N
y1:t [Ẑt ] =

t+1
1 X
(N − 1)t+1−ℓ
Nt
ℓ=1

ℓ−1
Y

i
X

Gτm ,τm+1 (yτjm ).

(4.10)

m=1 j=1

X

G0,τ1 Cℓ (τ , y1:t ).

(4.11)

τ ∈Tt,ℓ,t+1

In order to obtain a version of Proposition 4.4 for the general ci αSMC case,
we will require that the algorithm enforce a lower bound on a carefully chosen
notion of effective sample size called ∞-ESS. The ∞-ESS is a member of a family
of effective sample size measures we call p-ESS, which also includes two commonly
used definitions as special cases.
p
be the conjugate exponent
Definition 4.5. For parameter p ∈ [1, ∞], let p∗ , p−1
of p (so 1/p+1/p∗ = 1). The p-effective sample size (p-ESS) of the weight vector
w1:N ∈ RN
+ is
  1:N p∗
 kw1:N k1
p>1
kw
kp
1:N
(4.12)
ESSp (w ) ,
1:N
kw
k
1
 QN
p = 1.
n wn /kw1:N k1
n=1 (w

)

⊳

The following proposition highlights some elementary properties of p-ESS and
subsumes Proposition 1.4 (see Appendix A.3 for a proof).

16

J. H. HUGGINS AND D. M. ROY

Proposition 4.6. The p-ESS has the following properties:
(1) For all p ∈ [1, ∞], 1 ≤ ESSp (w1:N ) ≤ N . The lower bound is achieved if and
only if all but one of the weights is zero. The upper bound is achieved if and
only if all the weights are equal.
(2) For 1 < p < q ≤ ∞, ESSp (w1:N ) ≥ ESSq (w1:N ) ≥ N −(1−q∗ /p∗ ) ESSp (w1:N ),
with equality if and only if K weights are equal and the rest are zero.
(3) The 1-ESS satisfies
ESS1 (w1:N ) = lim ESSp (w1:N ) = ESSent (w1:N ) , eH(w

1:N

p↓1

where H(w1:N ) , −

wn
n kw 1:N k1

P

)

,

(4.13)

n

log kww
is the entropy.
1:N k
1

Parts (1) and (2) generalize their counterparts in Proposition 1.4 to all p ∈ [0, 1],
including the case p = 1. Part (3) shows that the 1-ESS corresponds to the entropic
ESS, which is a common choice of ESS in applications [Cornebise et al., 2008].
In order to obtain a bound on Ei,N
[Ẑt ], we will require a lower bound on the
y 1:i
1:t
∞-ESS of the weights, as formalized in Assumption 1.A. Our development follows
that of Whiteley et al. [2016], who used the 2-ESS lower bound guarantee to bound
the L2 norm of the weights in terms of their L1 norm. Similarly, we will use the
∞-ESS lower bound guarantee to bound the sup-norm of the weights in terms of
their L1 norm. Specifically, under Assumption 1.A, we have
ζN ≤ ESS∞ (Ws1:N ) =

kWs1:N k1
kWs1:N k1
=
,
kWs1:N k∞
supn Wsn

(4.14)

kW 1:N k

1
s
and so, for all n ∈ [N ] and s ∈ [t], Wsn ≤ ζN
. We can use this upper bound on
n
Ws to prove a result that is very similar to Proposition 4.4, but permits an arbitrary
adaptation scheme satisfying Assumptions 1.A and 4.D. (See Appendix A.4 for a
proof):

Theorem 4.7. If Assumptions 1.A and 4.D hold, then for all t ≥ 1, i ≥ 1, N ≥ i,
1
i
y1:t
, . . . , y1:t
∈ Et,
Ei,N
[Ẑt ]
y1:i

(4.15)

1:t

≤

t+1
X
1
N (ζN )t−1

X

(ζN )t+1−ℓ

ℓ=1 τ ∈Tt,ℓ,t+1



N −i
ζN

1(τ1 >1)

G0,τ1

ℓ−1
Y

i
X

Gτm ,τm+1 (yτjm ).

m=1 j=1

In particular, in the case of i = 1, we have
E1,N
y1:t [Ẑt ] ≤

t+1
X
1
N (ζN )t−1

X

ℓ=1 τ ∈Tt,ℓ,t+1

(ζN )t+1−ℓ



N −1
ζN

1(τ1 >1)

G0,τ1 Cℓy (τ ).
(4.16)

The gap between Proposition 4.4 and Theorem 4.7 is that most of the factors of
N − i in the former are replaced by factors of N in the latter. Luckily we are interested in the i = 1, 2 cases, so we expect the differences between the two quantities
to be fairly small. The following result, which is immediate upon expanding the
left-hand sides of Eqs. (4.10) and (4.15) and keeping only Ω(1/N ) terms, formalizes
this intuition:

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

17

1
i
Corollary 4.8. If Assumption SIR holds, then for all i, t ≥ 1, N > i, y1:t
, . . . , y1:t
∈
t
E ,
Pt Pi
i
Zt−1 s=1 j=1 G0,s Gs,t+1 (y1:t
) − ti
Ei,N
[
Ẑ
/Z
]
=
1
+
+ Θ(N −2 ).
(4.17)
t
t
y 1:i
1:t
N
1
i
If Assumptions 1.A and 4.D hold, then for all i, t ≥ 1, N > i, y1:t
, . . . , y1:t
∈ Et,
P
P
i
Zt−1 ts=1 ij=1 G0,s Gs,t+1 (y1:t
) − ζi
[
Ẑ
/Z
]
≤
1
+
+ Θ(N −2 ).
(4.18)
Ei,N
t
t
1:i
y 1:t
ζN

4.3. Quantitative bounds. Recall Assumptions 1.B and 1.C, either of which can
be used in conjunction with Theorem 4.7 to obtain uniform, quantitative bounds
on Ei,N
[Ẑt ] by following the approach of Andrieu et al. [2013]:
y 1:i
1:t

Assumption 1.B. The potentials satisfy gs , supx∈E gs (x) < ∞ for all 1 ≤ s ≤ t.
Assumption 1.C. There exists a constant β > 0 such that for any t, s ∈ N,
sup
x∈E

G0,t Gt,t+s (x)
≤ β.
G0,t+s

Assumption 1.C is implied by a standard “strong mixing” condition which is
often employed in SMC analyses (e.g., [Del Moral, 2004, Whiteley et al., 2016]).
See Andrieu et al. [2013] for details.
Proposition 4.9. If αs = 11/N for s ∈ [t − 1] and Assumption 1.B holds, then for
1
i
all t ≥ 1, i ≥ 1, N ≥ i, y1:t
, . . . , y1:t
∈ Et,
t #
"


t
Y
i
.
(4.19)
gs − 1 1 − 1 −
Ei,N
[Ẑt /Zt ] ≤ 1 + Zt−1
y 1:i
1:t
N
s=1

1
i
If Assumptions 1.A, 1.B and 4.D hold, then for all t ≥ 1, i ≥ 1, N ≥ i, y1:t
, . . . , y1:t
∈
t
E ,
"
#
t
t
Y
i
i,N
−1
Ey 1:i [Ẑt /Zt ] ≤ 1 + Zt
gs
1+
−1 .
(4.20)
1:t
ζN
s=1

Proof. The proof Eq. (4.19) is a straightforward generalization of that for [Andrieu et al.,
2013, Proposition 12] with some additional bookkeeping for i (instead of 2) fixed
trajectories. As for Eq. (4.20), we have
Ei,N
[Ẑt ] ≤
y 1:i
1:t

t+1
X

X

(ζN )−ℓ+1 G0,τ1

ℓ=1 τ ∈Tt,ℓ,t+1
t
Y

t+1 
X

ℓ−1
Y

i
X

Gτm ,τm+1 (yτjm )

m=1 j=1


t
gs
iℓ−1 (ζN )−ℓ+1
≤ Zt +
ℓ
−
1
s=1
ℓ=2
t  
t
Y X
t
gs
(ζN/i)−ℓ
= Zt +
ℓ
s=1
ℓ=1
#
"
t
t
Y
i
= Zt +
−1 .
gs
1+
ζN
s=1

(4.21)



18

J. H. HUGGINS AND D. M. ROY

Proposition 4.10. If αs = 11/N for s ∈ [t − 1] and Assumption 1.C holds, then
1
i
for all t ≥ 1, i ≥ 1, N ≥ i, y1:t
, . . . , y1:t
∈ Et,

t
i(β − 1)
.
(4.22)
Ei,N
[
Ẑ
/Z
]
≤
1
+
t
t
y 1:i
1:t
N

1
i
If Assumptions 1.A, 1.C and 4.D hold, then for all t ≥ 1, i ≥ 1, N ≥ i, y1:t
, . . . , y1:t
∈
t
E ,
t

iβ
Ei,N
[
Ẑ
/Z
]
≤
1
+
.
(4.23)
t
t
1:i
y 1:t
ζN

Proof. The proof of Eq. (4.22) is a simple generalization of that for [Andrieu et al.,
2013, Proposition 14]. As for Eq. (4.23), observe that for s ∈ [t + 1], G0,t+1 =
G
= G0,s ηs (Gs,t+1 ), so we can write for ℓ ∈ [t], τ ∈ Tt,ℓ,t+1 ,
G0,t+1 G0,t+1
0,s
Zt = G0,t+1 = G0,τk

ℓ−1
Y

ητi (Gτi ,τi+1 ).

(4.24)

i=1

Combined with Assumption 1.C and writing Ḡs,t , supx∈E Gs,t (x),
t+1
X

(ζN )−ℓ+1

X

G0,τ1

τ ∈Tt,ℓ,t+1

ℓ=1

≤ Zt + Zt

t+1
X
ℓ=2

t+1 
X

(ζN )−ℓ+1

ℓ−1
i
YX

Gτm ,τm+1 (yτjm )

i=1 j=1

X

τ ∈Tt,ℓ,t+1

ℓ−1 i
G0,τ1 Y X Ḡτi ,τi+1
G0,τ1 i=1 j=1 ητi (Gτi ,τi+1 )


t
= Zt
(ζN )−ℓ+1 (iβ)ℓ−1
ℓ−1
ℓ=1

t
iβ
= Zt 1 +
.
ζN

(4.25)



To compare Eqs. (4.19) and (4.20), consider the Θ(1/N ) terms, which are, respectively,
Qt
Qt
tiZt−1 s=1 gs
ti[Zt−1 s=1 g s − 1]
and
.
N
ζN
Thus, up to a −ti term and a factor of 1/ζ, the two bounds are of the same leading
order in 1/N . The −it is likely an artifact of the analysis while the 1/ζ term
accounts for there being only ζN “effective particles.” The differences between
Eqs. (4.22) and (4.23) are identical.
5. Bounding the Divergence of SMC Samplers
0,N
∗
Recall that P 0,N (dx1:t ) is the distribution of X̃1:t
∼ π1:t
, a single sample from
0,N
0,N
the αSMC estimator of π1:t . Equivalently, P
(dx1:t ) = E0,N [π1:t
](dx1:t ) ,
0,N
0,N 0,N
E [π1:t (dx1:t )] is the expected value of the random measure π1:t . As a first
application of our results from Section 4, we consider bounding the distance between the measures π1:t and P 0,N . That is, for some divergence d(µ||ν) between
measures, can we bound d(π1:t ||P 0,N )? To the best of our knowledge, there has

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

19

been minimal investigation of this question, with [Del Moral, 2004, Chapter 8] a
notable exception. For example, under Assumption SIR, the bound
c
(5.1)
dKL (P 0,N ||πt ) ≤ ,
N
can be extracted as a special case of a more general propagation-of-chaos result
[Del Moral, 2004, Theorem 8.3.2].
Let F1 be the set of functions f : R+ → R that are monotonically increasing or
decreasing and satisfy f (1) = 0. We consider the class of monotonic divergences of
the form

dµ1 
di,f (µ1 ||µ2 ) , µi f ◦
, i ∈ {1, 2}, f ∈ F1 .
(5.2)
dµ2

Table 1 lists some common divergences that can be written this way.
The following result characterizes the divergence between π1:t and P 0,N , only
assuming that f is concave.
−1 Pt
Proposition 5.1. Let Rt (y1:t ) , E1,N
y1:t [Ẑt /Zt ] and St , Zt
s=1 G0,s π1:t (Gs,t+1 )−
ζ. If Assumption 4.D holds, then for all concave f ∈ F1 ,
d1,f (π1:t ||P 0,N ) ≤ f (π1:t (Rt )).

In particular, if Assumption SIR holds, then


St − t
+ Θ(N −2 )
dKL (π1:t ||P 0,N ) ≤ log 1 +
N
S
−
t
t
dχ2 (π1:t ||P 0,N ) ≤
+ Θ(N −2 )
N
while if Assumptions 1.A and 4.D hold, then


St − ζ
−2
0,N
+ Θ(N )
dKL (π1:t ||P
) ≤ log 1 +
ζN
St − ζ
dχ2 (π1:t ||P 0,N ) ≤
+ Θ(N −2 ).
ζN
Proof. The general statement follows by applying Jensen’s inequality to Eq. (5.2),
then using Theorem 4.2. The special cases correspond to using the KL divergence
(version 1) and χ2 distance rows of Table 1 and applying Corollary 4.8.

Similar results for (sequential) importance sampling are included in Appendix A.5.
We can also consider the divergence between π1:t and P 0,N when E1,N
y1:t [Ẑt ] is
uniformly bounded:
Proposition 5.2. If Assumption 4.D holds and E1,N
y1:t [Ẑt /Zt ] ≤ Bt,N for all y1:t ∈
E t , then for all increasing f ∈ F1
di,f (π1:t ||P 0,N ) ≤ f (Bt,N )
and for all decreasing f ∈ F1 ,
−1
di,f (P 0,N ||π1:t ) ≤ f (Bt,N
).

20

J. H. HUGGINS AND D. M. ROY

Table 1. Divergences of the form Eq. (5.2). The operator (a)+
gives the positive part of a ∈ R.
Symbol
Name
KL divergence (version 1)
dKL
KL divergence (version 2)
dKL
dχ2
χ2 distance
total variation distance (version 1) dT V
total variation distance (version 2) dT V

i
1
2
1
2
2

f
a 7→ log a
a 7→ − log a
a 7→ a − 1
a 7→ (a − 1)+
a 7→ (1 − a)+

In particular,
dKL (π1:t ||P 0,N ) ≤ log Bt,N

dχ2 (π1:t ||P 0,N ) ≤ Bt,N − 1
Bt,N − 1
≤ Bt,N − 1.
dT V (π1:t , P 0,N ) ≤
Bt,N

Proof. The general statements follow immediately from Eq. (5.2) and Theorem 4.2.
The special cases correspond to using the KL divergence (version 1), χ2 distance,
and total variation distance (version 2) rows of Table 1. The second total variation
inequality holds since Bt,N ≥ 1.

The bounds in Proposition 5.2 for KL divergence, χ2 distance, and total variation distance are asymptotically equivalent if Bt,N → 1 as N → ∞. Combining
Proposition 5.2 with, for example, Proposition 4.10, yields quantitative bounds for
SIR and αSMC:
Corollary 5.3. If Assumptions 1.C and SIR hold, then
t(β − 1)
N
t(β
− 1)
+ O(N −2 )
dχ2 (π1:t ||P 0,N ) ≤
N
t(β − 1)
dT V (π1:t , P 0,N ) ≤
+ O(N −2 ).
N + t(β − 1)

dKL (π1:t ||P 0,N ) ≤

(5.3)
(5.4)
(5.5)

If Assumptions 1.A, 1.C and 4.D hold, then

tβ
ζN
tβ
+ O(N −2 )
dχ2 (π1:t ||P 0,N ) ≤
ζN
tβ
dT V (π1:t , P 0,N ) ≤
+ O(N −2 ).
ζN + tβ

dKL (π1:t ||P 0,N ) ≤

(5.6)
(5.7)
(5.8)

6. The α-Particle Gibbs Sampler
As a second application of our results from Section 4, we consider the mixing
properties of particle Gibbs with adaptive resampling. Recall from Section 1.3 that
we introduce a global parameter θ ∈ Θ with prior distribution ̟(dθ). Replace
Ms by Msθ and gs by gsθ , then parameterize the other quantities defined previously

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

21

in terms of Ms and gs by θ. The target distribution on the product space (Θ ×
Y, B(Θ × Y )), Y , E t , is
π(dθ × dy) , γ(dθ × dy)/Z,

(6.1)

where
γ(dθ × dy) ,

t
Y

gsθ (ys )Msθ (ys−1 , dys )̟(dθ)

and Z , γ(1).

(6.2)

s=1

Particle Gibbs samplers have kernels of the form πy (dθ)Πθ (y, dz), where Πθ (y, dz)
is an SMC-based kernel with invariant distribution πθ . The standard PG sampler
employs the iterated conditional SMC (i-cSMC) kernel [Andrieu et al., 2013]: that
1,N
is, Πθ = Pθ,y
and require Assumption SIR to hold.
We now introduce the novel α-particle Gibbs (αPG) sampler, which employs the
1,N
1,N
iterated conditional αSMC (i-cαSMC) kernel Pθ,y
, so Πθ = Pθ,y
. In Appendix A.6
we prove that the i-cαSMC kernel is reversible with respect to πθ and hence has
invariant distribution πθ .
The first step to proving mixing results for the i-cαSMC kernel and the αPG
sampler is to use Theorem 4.2 to obtain a sufficient condition for the i-cαSMC
transition kernel to satisfy a minorization condition.
Proposition 6.1. If Assumption 4.D holds and E2,N
θ,y 1:2 [Ẑt /Zt ] ≤ Bt,N for all θ ∈ Θ
and y 1 , y 2 ∈ Y , then
1,N
Pθ,y
(y, dx) ≥ εt,N πθ (dx),

where εt,N ,

(6.3)

(1−κ′N )t
Bt,N .

The constant ǫt,N , which determines mixing speed, can be found explicitly using
the quantitative bounds from Section 4.3. For example, using Assumption 1.C we
obtain the following:
Corollary 6.2. If Assumptions 1.A, 1.C and 4.D hold, then for all y ∈ Y ,
1,N
Pθ,y
(y, dx) ≥ εt,N πθ (dx),

(6.4)

where
εt,N ,
Furthermore, if ζN ≥

2βt
C(1−κ′N )−κ′N t



1 − κ′N
2β
1 + ζN

t

.

(6.5)

for some constant C > 0, then
ǫt,N ≥ e−C .

(6.6)

In particular, assuming κ′N ≤ B/N for some constant B ≥ 1, if N ≥ t/C + B, then
 C

ǫt,N ≥ exp − (2β + ζB) .
ζ

(6.7)

22

J. H. HUGGINS AND D. M. ROY

Proof. The first part follows from Propositions 6.1 and 4.10. For the second part,
we then have


−t

2β −t
1 + ζN
1
2β
′
= 1+
+ κN
(6.8)
ǫt,N ≥
1 − κ′N
1 − κ′N ζN
−t

C
≥ e−C .
(6.9)
≥ 1+
t

The final part follows after noting that if κ′N ≤ B/N , then






1
1
2β
2β
2β
1
′
+ κN ≤
+ B/N =
+B .
1 − κ′N ζN
1 − B/N ζN
N −B ζ

(6.10)


Remark 6.3. In the case of the i-cSMC kernel, Corollary 6.2 is almost as tight as
[Andrieu et al., 2013, Corollary 14]: the former result replaces β − 1 with β.
▽
The minorization condition Eq. (6.3) implies uniform ergodicity and a number
of other types of convergence guarantees for the i-cSMC process. The following
generalizes [Andrieu et al., 2013, Theorem 1], which applies only to the i-cSMC
kernel and the PG sampler.
Theorem 6.4. Assume that Assumptions 1.A and 4.D hold.
I. Let N ≥ 2, and consider the i-cαSMC process with kernel P = Pθ1,N .
1. P is reversible with respect to π and defines a positive operator,
2. If the potentials are bounded then there exists ǫt,N = 1 − Ct /N such that
(i) for all y ∈ Y , P (y, dz) ≥ ǫt,N πθ (dz),
(ii) for every measure ν ≪ πθ and k ≥ 1,
dχ2 (νP k ||πθ ) ≤ dχ2 (ν||πθ )(1 − ǫt,N )k ,

(6.11)

(iii) for every y ∈ Y and k ≥ 1,

dT V (δy P k , πθ ) ≤ (1 − ǫt,N )k ,

(6.12)

ǫt,N ≥ εB,C,ζ > 0.

(6.13)

3. If Assumption 1.C also holds and there is a constant B > 0 such that
κ′N ≤ B/N , then for every C > 0, there exists εB,C,ζ > 0 such that for
N ≥ t/C + B and all t > 1,
II. If there exists β ≥ 1 such that, for all t, s ∈ N,
π- ess sup
θ,x

Gθ0,t Gθt,t+s (x)
≤ β,
Gθ0,t+s

(6.14)

or if
π- ess sup
θ

Qt

g θs
< ∞,
γθ (1)
s=1

(6.15)

then the αPG chain is geometrically ergodic whenever the Gibbs sampler is
geometrically ergodic.
Proof. Part I.1 follows from Lemma A.4. Parts I.2-3 follow from Proposition 6.1,
Corollary 6.2, and [Andrieu et al., 2013, Proposition 31]. Part II follows from
[Andrieu et al., 2013, Section 7].


SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

23

Remark 6.5. Part I.3 means that if Assumption 1.C holds, then scaling N linearly
with t ensures a uniform convergence rate, as measured by χ2 -divergence or total
variation distance.
▽
Appendix A. Additional Proofs
A.1. Technical Lemma.
Lemma A.1. Let X and Y be random elements in Borel spaces (S, S) and (T, T ),
respectively, let ψ : S × T → R+ be a measurable, and let µ be the distribution of
X. If

then ν ≪ µ and

ν = E[ψ(X, Y )δX ],

(A.1)

dν
(X) = E[ψ(X, Y ) | X] a.s.
dµ

(A.2)

Proof. Because S is Borel, there exists an f satisfying f (X) = E[ψ(X, Y ) | X] a.s.
It follows from the chain rule of conditional expectation and then some elementary
manipulations that, for all A ∈ S,
Z
f (x)µ(dx),
ν(A) = E[f (X)δX (A)] = E[f (X)1A (X)] =
A

and so f is a version of the Radon–Nikodym derivative dν/dµ.



A.2. Proof of Theorem 4.2, i ≥ 1 case. First observe that we can write the
ci αSMC kernel as
h
i
i,N
δX̃ ∗ (dx1:t )
(A.3)
P i,N (y 1:i
1:t , dx1:t ) = Ey 1:i
1:t
1:t
hP
i
1:N
1:N
= Ei,N
(A.4)
k1:t ∈[N ]t Ik1:t (X 1:t , A1:t , dx1:t ) ,
y 1:i
1:t

where

∗
1:N
Ik1:t (x1:N
1:t , a1:t , dx1:t ) , δx̃kt (dx1:t )1(kt = at )
1:t

t−1
Y

1(ks = aks s+1 ).

(A.5)

s=1

Next note that
N
X

k1 =1

=

1:i
1:i
t
1(x̃k1:t
∈ S)Pi,N
[X11:N ∈ dx1:N
1 , F1 = f1 ]
y 1:i

N
X

k1 =1

≥

(A.6)

1:t

t
1(x̃k1:t
∈ S)C1i D(f11:i )

N Z
N C1i X
C1i+1 k =1 E
1

N
i
Y
Y
1
fj
M1 (dxn1 )
δyj (dx11 )
1
N
1:i
j=1

t
1(x̃k1:t
∈ S)C1i+1 D(f11:i , k1 )

×

N
Y

n∈f
/ 11:i ,k1

(A.7)

n∈f
/ 1

i
Y

j=1

M1 (dxn1 )M1 (dx1 )

1
fj 1
δyj (dx11 ) δx1 (dxk11 )
1
N
N
(A.8)

24

J. H. HUGGINS AND D. M. ROY

≥

N Z
N C1i X
i+1
1:i
1:i
t
1(x̃k1:t
[X11:N ∈ dx1:N
∈ S)Pi+1,N
= k1 ]M1 (dx1 )
1 , F1 = f1 , F1
y 1:i
1:t ,x1:t
C1i+1 k =1 E
1
(A.9)

For the remainder of the proof, to keep notation compact when writing laws, instead of writing, e.g., Xs1:N ∈ x1:N
or Fs1:i = fs1:i , whenever a random variable is
s
instantiated to be (the differential) of the lowercase version of itself, we will write
only the random variable: for example, Xs1:N or Fs1:i . Now, for s = 2, . . . , t,
N
X

ks =1

1:N
1:N
1:i
t
s
1(x̃k1:t
∈ S)1(ks−1 = aks−1
)Pi,N
[Xs1:N , As−1
, Fs1:i | X 1:N
]
1:s−1 , A1:s−2 , Fs−1(A.10)
y1:i
1:t

N
X

,

ks =1

t
s
1(x̃k1:t
∈ S)1(ks−1 = aks−1
)Csi D(fs1:i )

×

≥

Y

n∈f
/ s1:i

N
X

i
Y

fjfj

fj

fj

j
s s−1
s
= fs−1
)
αs−1
δysj (dxss )1(as−1

j=1

an

s−1
n
1:N
, x1:N
rn (ans−1 |ws−1
1:s−1 )Ms (xs−1 , xs )

Csi

k k

s s−1
i+1
ks =1 Cs αs−1

(A.11)
Z

E

t
∈ S)Csi+1 D(fs1:i , ks )
1(x̃k1:t

i
Y

fjfj

fj

fj

j
s s−1
s
= fs−1
)
αs−1
δysj (dxss )1(as−1

j=1
k k

s s−1
s
)
δxs (dxks s )1(ks−1 = aks−1
× αs−1

×

Y

n∈f
/ s1:i ,ks

(A.12)
an

s−1
1:N
1:N
rn (ans−1 |ws−1
, x1:s−1
)Ms (xs−1
, xns )

k

s−1
1:N
, x1:N
× rks (ks−1 |ws−1
1:s−1 )Ms (xs−1 , dxs )

=

N
X

Csi

i+1 ks ks−1
ks =1 Cs αs−1

Z

E

k

s−1
1:N
t
1(x̃k1:t
, x1:N
∈ S)rks (ks−1 |ws−1
1:s−1 )Ms (xs−1 , dxs ) (A.13)

1:N
i+1
1:i
i+1
1:i
[Xs1:N , A1:N
= ks | X 1:N
Pi+1,N
s−1 , Fs , Fs
1:s−1 , A1:s−2 , Fs−1 , Fs−1 = ks−1 ]
y1:i ,x1:t
1:t

Using Eqs. (A.9) and (A.14), we have (note that the terms such as those involving
a0 should be ignored)
X

k1:t ∈[N ]t

t
1(x̃k1:t
∈ S)1(kt = a∗t )

t
Y

s=2

1:N
1:i
s
1(ks−1 = aks−1
)Pi,N
[X 1:N
1:t , A1:t , F 1:t ]
y 1:i
1:t

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

X

=

k1:t ∈[N ]t

i,N
1:N
1:N
t
1(x̃k1:t
∈ S)1(kt = a∗t )Pi,N
[A∗t | X 1:N
, F11:i ]
1:t , A1:t−1 ]Py 1:i [X1
y1:i
1:t

t
Y

×

≥

X

k1:t ∈[N ]t

s=2

Z

×
×

Et

1:t

1:N
1:N
1:i
1:i
s
[Xs1:N , A1:N
1(ks−1 = aks−1
)Pi,N
s−1 , Fs | X 1:s−1 , A1:s−2 , Fs−1 ]
y1:i
1:t

(A.14)

t
N 1(x̃k1:t
∈ S)1(kt = a∗t )
Qt
Qt−1 ks ks−1
( s=1 Csi+1 /Csi )( s=1 αs−1
)

t
Y

s=2
t
Y

s=2

k

a∗
a∗
wt t gt (xt t )
PN
n
n
n=1 wt gt (xt )

s−1
1:N
, x1:N
rks (ks−1 |ws−1
1:s−1 )Ms (xs−1 , dxs )

1:t

1:N
1:i
i+1
[X 1:N
Pi+1,N
1:t , A1:t−1 , F 1:t , F1:t
y1:i
1:t ,x1:t

Z

X

Zt

S k ∈[N ]t
1:t

(A.15)

1:N
i+1
1:i
i+1
1:i
Pi+1,N
[Xs1:N , A1:N
= ks | X 1:N
s−1 , Fs , Fs
1:s−1 , A1:s−2 , Fs−1 , Fs−1 = ks−1 ]
y1:i ,x1:t

i+1
1:i
1:i
[X11:N ∈ dx1:N
× Pi+1,N
= k1 ]M1 (dx1 )
1 , F1 = f1 , F1
y1:i
1:t ,x1:t
X Z
N 1(x1,t ∈ S)
=
Qt
P
i+1
n
n
t
( s=1 Cs /Csi ) N
n=1 wt gt (xt )
k1:t ∈[N ]t E

=

25

Ẑt

Qt

i+1
i
s=1 Cs /Cs

(A.16)

= k1:t ]γ1:t (dx1:t )

1:N
1:i
i+1
Pi+1,N
[X 1:N
1:t , A1:t−1 , F 1:t , F1:t = k1:t ]π1:t (dx1:t ),
y1:i ,x1:t
1:t

(A.17)

from which Eq. (4.2) follows.
To prove Eq. (4.4), first note that under Assumption 4.D, the normalization
constants for the c2 αSMC process are given by
C12 ,
and, for s = 2, . . . , t,
Cs2
Thus, Cs2 ≤

1
1−κN

,

1−

N
X

N
N −1

(A.18)

1
2
kfs−1
kfs−1
αs−1
αs−1

k=1

for s = 2, . . . , t and hence Cs2 ≤

!−1

1
1−κ′N

.

(A.19)

for all s ∈ [t].

A.3. Proof of Proposition 4.6. For (1), the fact that ESS1 (w1:N ) = ESSent (w1:N )
is a straightforward algebraic manipulation. To prove the limit equality, observe
that, using the Taylor series for xp and log(1 + x), we have
!1/(1−p)
kw1:N kp1
(A.20)
lim PN
n p
p→1
n=1 (w )
!1/(1−p)
P∞
1:N
k1 (p − 1)k logk (kw1:N k1 )/k!
k=0 kw
= lim
(A.21)
PN P∞
k
n
k
n
p→1
n=1
k=0 w (p − 1) log (w )/k!

26

J. H. HUGGINS AND D. M. ROY

= lim

x→∞

P∞

1:N
k1 x−k logk (kw1:N k1 )/k!
k=0 kw
PN P∞
n −k logk (wn )/k!
n=1
k=0 w x

!x

(A.22)

!x
P
k
−k
1:N
exp(log(1 + ∞
x
log
(kw
k
)/k!))
1
k=1
(A.23)
= lim
P∞ PN
−k logk (wn )/k!))
x→∞
exp(log(1 + k=1 n=1 wn kw1:N k−1
1 x
P∞
P∞
exp(x m=1 (−1)m+1 [ k=1 x−k logk (kw1:N k1 )/k!]m )
= lim
PN
P
P∞
n
1:N k−1 x−k logk (wn )/k!]m )
m+1 [ ∞
x→∞ exp(x
1
n=1 w kw
k=1
m=1 (−1)
(A.24)
exp(log(kw1:N k1 ) + Θ(x−1 ))


Q
N
x→∞
n )w n /kw 1:N k1 + Θ(x−1 )
exp log
(w
n=1

= lim

= QN



kw1:N k1

n=1 (w

n )w n /kw 1:N k1

.

(A.25)

(A.26)

To prove the remaining parts, we make repeated use of the following:
1:N
Fact. For 1 ≤ r < s ≤ ∞, and any vector w1:N ∈ RN
ks ≤ kw1:N kr ≤
+ , kw
1/r−1/s
1:N
N
kw ks , with the lower (upper) bound achieved if and only if w1:N has
one non-zero entry (w1:N has all equal entries).

For (2), apply the Fact with r = 1, s = p > 1, and note that in this case
1/r − 1/s = 1 − 1/p = 1/p∗ . We then have 1 ≤ kw1:N k1 /kw1:N kp ≤ N 1/p∗ , proving
the result for p > 1. For p = 1, the result follows from part (1) and elementary
properties of the entropy.
For (3), in the case that p > 1, note that
kw1:N kq1∗ −p∗ ≥ N (q∗ −p∗ )/q∗ kw1:N kqq∗ −p∗ = N 1−p∗ /q∗ kw1:N kqq∗ −p∗
= N −p∗ (1/p−1/q) kw1:N kqq∗ −p∗ ,

(A.27)

where the final equality follows since

1 − p∗ /q∗ = 1 − p∗ (1 − 1/q) = 1 − p∗ + p∗ /q = −p∗ /p + p∗ /q.

(A.28)

We conclude that
 1:N p∗
kw1:N kp1∗
kw k1
≥
kw1:N kp
N p∗ (1/p−1/q) kw1:N kpq ∗

(A.29)

≥

(A.30)

kw1:N kp1∗
kw1:N kq1∗ −p∗
p
∗
N p∗ (1/p−1/q) kw1:N kq N −p∗ (1/p−1/q) kvkqq∗ −p∗
 1:N q∗
kw k1
=
kw1:N kq
 1:N q∗
kw k1
≥
kw1:N kp
 1:N p∗ −q∗  1:N p∗
kw k1
kw kp
=
kw1:N k1
kw1:N kp
 1:N p∗
kw k1
,
≥ N −(p∗ −q∗ )/p∗
kw1:N kp

(A.31)
(A.32)
(A.33)
(A.34)

where the first, third, and fourth inequalities follow from the Fact and the second
follows from Eq. (A.27).

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

27

The case of p = 1 follows from the p > 1 case and part (1).
A.4. Proof of Theorem 4.7. We prove the result for i = 1. The general case
follows from straightforward modifications.
For t ≥ 1, let Qt (xt−1 , dxt ) , gt−1 (xt−1 )Mt (xt−1 , dxt ), and for 0 ≤ s < t, let
Qs,t , Qs+1 Qs+2 · · · Qt ,

(A.35)

so Qt,t+1 = Qt . By convention Qt,t (xt , dyt ) = δxt (dyt ) and Q0,t (dxt ) is a measure,
not a probability kernel. Notice that for s ∈ [t], xs ∈ E, and φt : E → R,
Qs,t (xs )(φt ) = E[φt (ξt )gs:t−1 (ξs:t−1 ) | ξs = xs ]

(A.36)

and Q0,t (φt ) = M1 Q1,t (φt ). Generalizing these identities, we will abuse notation
and write, for s ∈ [t], xs ∈ E, and φs,t : E t−s → R,
Qs,t (xs )(φs,t ) , E[φs:t (ξs:t )gs:t−1 (ξs:t−1 ) | ξs = xs ]

(A.37)

and Q0,t (φ1:t ) , M1 Q1,t (φ1:t ). Note that Gs,t (y) = Qs,t (y)(1) for s ∈ [t − 1] and
G0,t = Q0,t (1).
We will use the abbreviated notation Qks,t (·) = Qs,t (·)(Xsk ) or Qs,t (·)(xks ), Gks,t =
Gs,t (Xsk ) or Gs,t (xks ), Gys,t = Gs,t (ys ), gsk = gs (Xsk ) or gs (xks ), and gsy = gs (ys ). The
variables are Xsk inside expectations and xks outside expectations. Throughout the
proof, when limits of a sum are not specified, the sum is from 1 to N .
1:N
1:N
1
, where by conven, and F1:s
, A1:s−1
Let Fs be the σ-algebra generated by X 1:s
tion we let F0 be the trivial σ-algebra. The proof relies on the following lemma.
Lemma A.2. If y1:t ∈ E t , then
(1) for s = 2, . . . , t and any functions φns : E → R, n ∈ [N ],
#
"
X
1,N
n
n
Ey1:t
φs (Xs ) | Fs−1
n

=

X

f 1f 1

f1

s s−1
αs−1
φss (ys ) +

X X X
fs1 n6=fs1

fs1

f 1f 1

s s−1
αs−1

k

(A.38)

k
αnk
s−1 ws−1 k
Qs−1,s (φns );
n
ws

(2) for τ ∈ [t − s],
"
#
X
X
1 X n n
1,N
n n
n
Ey1:t
Ws Gs,s+τ | Fs−1 ≤
ws−1 gs−1 Gys,s+τ +
ws−1
Gns−1,s+τ ;
ζN
n
n
n
(A.39)
and
(3) for s = 1, . . . , t − 1,
where,
At−s , (ζN )−s+1

X
n

Bt−s , (ζN )−s+1

X
n

i
h
N E1,N
y1:t Ẑt | Ft−s ≤ At−s + Bt−s ,


n
n 
wt−s
gt−s



n 
wt−s

s
X
X

ℓ=1 τ ∈Tt,ℓ,s

s
X
X

ℓ=1 τ ∈Tt,ℓ,s

(A.40)


(ζN )s−1−ℓ Gyt−s+1,τ1 Cℓy (τ ) (A.41)


(ζN )s−ℓ Gnt−s,τ1 Cℓy (τ ) .

(A.42)

28

J. H. HUGGINS AND D. M. ROY

Proof. For (1),
"
#
X
1,N
n
n
Ey1:t
φs (Xs ) | Fs−1
n

=

fs1

=

1
fs1 fs−1
αs−1

XX

k6=fs1

−f 1
as−1s

XX
fs1

1
fs1 fs−1

αs−1

−f 1
as−1s

=

X

=

X

Y

1,N
1:N
rk (aks−1 |ws−1
, x1:N
1:s−1 )Ey1:t

"
X
n

φns (Xsn ) | Fs−1 , A1:N
s−1

=

1
a1:N
s−1 , Fs



kak
ak
ak
s−1
s−1
s−1
i
h
Y αs−1
X
n
1
ws−1
gs−1 (xs−1
)
as−1

φfss (ys ) +
E φns (ξs ) | ξs−1 = xs−1
k
w
s
1
1

k6=fs

n6=fs

f 1f 1

f 1f 1

f1

X X X

s s−1
αs−1

f 1f 1

f1

X X X

s s−1
αs−1

s s−1
αs−1
φss (ys ) +

fs1 n6=fs1

fs1

s s−1
αs−1
φss (ys ) +

fs1 n6=fs1

fs1

φns (x)

k

f 1f 1

k

wsn Gs,s+τ (x),

For (2), choosing
=
"
#
X
E1,N
Wsn Gns,s+τ | Fs−1
y1:t

k
k

αnk
s−1 ws−1 gs−1 (xs−1 )  n
E φs (ξs ) | ξs−1 = xks−1
n
ws

k
αnk
s−1 ws−1 k
Qs−1,s (φns )
wsn

we have
(A.43)

n

=

X

f 1f 1

f1

s s−1
αs−1
ws s Gys,s+τ +

=

X

1
fs1 fs−1
f1
αs−1
ws s

+

≤ Gys,s+τ

X
fs1

f 1f 1

s s−1
αs−1

f 1f 1

s s−1
αs−1

k

X X X
fs1 n6=fs1

fs1

=

X X X
fs1 n6=fs1

fs1

Gys,s+τ

k
αnk
s−1 ws−1 k
Qs−1,s (wsn Gs,s+τ )
wsn

(A.44)
1
fs1 fs−1
k
k
αs−1
αnk
s−1 ws−1 Gs−1,s+τ

(A.45)

k

1
kws1:N k1 X X X fs1 fs−1
k
k
αs−1 αnk
+
s−1 ws−1 Gs−1,s+τ
ζN
1
n

fs

(A.46)

k

Gys,s+τ X n X X nk k
ws +
αs−1 ws−1 Gks−1,s+τ
ζN
n
n

(A.47)

k

X
Gys,s+τ X X nk k
k
k
αs−1 ws−1 gs−1
+
ws−1
Gks−1,s+τ
ζN
n
k
k
X
1 X k
y
k
k
ws−1 gs−1
Gs,s+τ +
ws−1
Gks−1,s+τ ,
=
ζN

=

k

=

fs1

(A.48)
(A.49)

k

where the inequality follows from Assumption 1.A, and we have repeatedly used
Assumption 4.D.
To show (3), we start by using (2) with s = t and τ = 1:
"
#
X
1 X k k y ζN X m m
1,N
n n
wt−1 gt−1 gt +
w G
(A.50)
Ey1:t
Wt gt | Ft−1 =
ζN
ζN m t−1 t−1,t+1
n
k

= At−1 + Bt−1 ,

(A.51)

Hence, (3) holds for s = 1. We now assume that the bound holds for some
s ∈ {1, . . . , t − 2} and establish that it also holds for s + 1. Using the inductive

#

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

29

hypothesis,
i
i
h
h
i
h
1,N
1,N
N E1,N
y1:t Ẑt | Ft−s−1 = Ey1:t N Ey1:t Ẑt | Ft−s | Ft−s−1

(A.52)

≤ E1,N
y1:t [At−s + Bt−s | Ft−s−1 ] .

(A.53)

Using (2), we have
(A.54)
A , E1,N
y1:t [At−s | Ft−s−1 ]


"
#
s
X
X
X
y
−s+1 
s−1−ℓ y
1,N
n
n
= (ζN )
(ζN )
Gt−s+1,τ1 Cℓ (τ ) Ey1:t
Wt−s gt−s
n

ℓ=1 τ ∈Tt,ℓ,s

(A.55)



≤ (ζN )−s 
×

X

s
X
X

ℓ=1 τ ∈Tt,ℓ,s



(ζN )s−1−ℓ Gyt−s+1,τ1 Cℓy (τ )

y
n
n
gt−s−1
gt−s
wt−s−1

+ ζN

n

X

n
wt−s−1
Gnt−s−1,t−s+1

n

(A.56)

!

and
(A.57)
B , E1,N
y1:t [Bt−s | Ft−s−1 ]


#
"
s
X
X
X
−s+1 
n
= (ζN )
(ζN )s−ℓ E1,N
Wt−s
Gnt−s,τ1 | Ft−s−1 Cℓy (τ )
y1:t
n

ℓ=1 τ ∈Tt,ℓ,s

(A.58)



≤ (ζN )−s 

s
X
X

(ζN )s−ℓ

X

n
n
wt−s−1
gt−s−1
Gyt−s,τ1

n

ℓ=1 τ ∈Tt,ℓ,s

+ ζN

X

n
wt−s−1
Gnt−s−1,τ1

n

(A.59)

Hence,
A + B ≤ (ζN )−s

X
n

+ (ζN )−s



s
X
X
n
n

wt−s−1
gt−s−1
(ζN )s−1−ℓ Gyt−s,t−s+1 Gyt−s+1,τ1 Cℓy (τ )
ℓ=1 τ ∈Tt,ℓ,s

X
n

+ (ζN )−s

X
n

+ (ζN )−s

X
n



s
X
X
n
n

(ζN )s−ℓ Gyt−s,τ1 Cℓy (τ )
wt−s−1
gt−s−1
ℓ=1 τ ∈Tt,ℓ,s



n

wt−s−1



n

wt−s−1

s
X
X

ℓ=1 τ ∈Tt,ℓ,s

s
X
X

ℓ=1 τ ∈Tt,ℓ,s



(ζN )s−ℓ Gnt−s−1,t−s+1 Gyt−s+1,τ1 Cℓy (τ )


(ζN )s−ℓ+1 Gnt−s−1,τ1 Cℓy (τ ) .

(A.60)

!



Cℓy (τ ) .

30

J. H. HUGGINS AND D. M. ROY

Summing the parenthesized double sums of the first two terms yields
s
s
X
X
X
X
(ζN )s−ℓ Gyt−s,τ1 Cℓy (τ )
(ζN )s−1−ℓ Gyt−s,t−s+1 Gyt−s+1,τ1 Cℓy (τ ) +
ℓ=1 τ ∈Tt,ℓ,s

ℓ=1 τ ∈Tt,ℓ,s

=

s+1
X

X

(ζN )s−ℓ Gyt−s,τ1 Cℓy (τ ) +

ℓ=1 τ ∈Tt,ℓ,s+1
τ1 =t−s+1

s+1
X

(ζN )s−ℓ Gyt−s,τ1 Cℓy (τ )

X

ℓ=1 τ ∈Tt,ℓ,s+1
τ1 >t−s+1

(A.61)
=

s+1
X

X

(ζN )s−ℓ Gyt−s,τ1 Cℓy (τ ),

(A.62)

ℓ=1 τ ∈Tt,ℓ,s+1

so the first two terms are equal to At−(s+1) . Summing the parenthesized double
sums of the last two terms yields
s
s
X
X
X
X
(ζN )s−ℓ+1 Gnt−s−1,τ1 Cℓy (τ )
(ζN )s−ℓ Gnt−s−1,t−s+1 Gyt−s+1,τ1 Cℓy (τ ) +
ℓ=1 τ ∈Tt,ℓ,s

ℓ=1 τ ∈Tt,ℓ,s

=

s+1
X

X

(ζN )s−ℓ+1 Gnt−s−1,τ1 Cℓy (τ ) +

ℓ=1 τ ∈Tt,ℓ,s+1
τ1 =t−s+1

s+1
X

(ζN )s−ℓ+1 Gnt−s−1,τ1 Cℓy (τ )

X

ℓ=1 τ ∈Tt,ℓ,s+1
τ1 >t−s+1

(A.63)
=

s+1
X

X

(ζN )s−ℓ+1 Gnt−s−1,τ1 Cℓy (τ ),

(A.64)

ℓ=1 τ ∈Tt,ℓ,s+1



so the last two terms are equal to Bt−(s+1) .
Using part (3) of Lemma A.2 with s = t − 1, we have
h i
1,N
N E1,N
y1:t Ẑt ≤ Ey1:t [A1 + B1 ] .

(A.65)

Therefore,

−t+2 1,N
Ey1:t
E1,N
y1:t [A1 ] = (ζN )

"
X
n

#  t−1
X
g1n 

X

ℓ=1 τ ∈Tt,ℓ,t−1


t−1
X
= (ζN )−t+2 (Gy1,2 + (N − 1)G0,2 ) 



(ζN )t−2−ℓ Gy2,τ1 Cℓy (τ )
X

′
ℓ=1 τ ∈Tℓ,t−1

(A.66)


(ζN )t−2−ℓ Gy2,τ1 Cℓy (τ )

(A.67)

and

t−1
X
−t+2 
[B
]
=
(ζN
)
E1,N
1
y
1:t

#



Gn1,τ1 Cℓy (τ ) (A.68)

(ζN )t−1−ℓ E1,N
y1:t

X

(ζN )t−1−ℓ (Gy1,τ1 + (N − 1)G0,τ1 )Cℓy (τ ) .

ℓ=1 τ ∈Tt,ℓ,t−1


t−1
X
= (ζN )−t+2 

"
X

X

ℓ=1 τ ∈Tt,ℓ,t−1

n



(A.69)

SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

31

Hence, using arguments analogous to those from the proof of Lemma A.2 and the
fact that G0,1 = 1 yields
E1,N
y1:t [A1 + B1 ]
= (ζN )−t+2

t
X
X

(ζN )t−1−ℓ Gy1,τ1 Cℓy (τ )

(A.70)

ℓ=1 τ ∈Tt,ℓ,t

+

t
X
X
N −1
(ζN )t−ℓ G0,τ1 Cℓy (τ )
(ζN )−t+2
ζN

(A.71)

ℓ=1 τ ∈Tt,ℓ,t

= (ζN )−t+2

t+1
X

X

(ζN )t−ℓ G0,τ1 Cℓy (τ )

ℓ=1 τ ∈Tt,ℓ,t+1
τ1 =1

(A.72)

t+1

X
N −1
(ζN )−t+2
+
ζN

X

(ζN )

G0,τ1 Cℓy (τ )



1(τ1 >1)

t−ℓ

ℓ=1 τ ∈Tt,ℓ,t+1
τ1 >1

−t+2

= (ζN )

t+1
X

X

t−ℓ

(ζN )

ℓ=1 τ ∈Tt,ℓ,t+1

N −1
ζN

G0,τ1 Cℓy (τ ).

(A.73)

A.5. Divergence of importance samplers. The key quantity in this section is
the variance of the potentials:
h


2 i
Vt , V Zt−1 g1:t (ξ1:t ) = E Zt−1 g1:t (ξ1:t ) − 1

(A.74)

Theorem A.3. If αs = I N for all s ∈ [t − 1], then


Vt
dKL (π1:t ||P
) ≤ log 1 +
N
Vt
dχ2 (π1:t ||P 0,N ) ≤ .
N
0,N

(A.75)
(A.76)

Proof. By Theorem 4.2 and Jensen’s inequality
 
dP 0,N
Zt
(x1:t ) = E1,N
x1:t
dπ1:t
Ẑt
≥

PN
E1,N
x1:t [
k=1

(A.77)
N

k )]
Zt−1 g1:t (X̃1:t
N
.
=
N − 1 + Zt−1 g1:t (x1:t )

(A.78)
(A.79)

32

J. H. HUGGINS AND D. M. ROY

By definition of the χ2 divergence,

dπ1:t
−1
dP 0,N


dπ1:t dπ1:t
= M1:t
−1
dP 0,N dM1:t


N − 1 + Zt−1 g1:t −1
Zt g1:t − 1
≤ M1:t
N

S,N
dχ2 (π1:t ||π̄1,t
) = π1:t



M1:t (Zt−1 g1:t )2 − 1
N
V[Zt−1 g1:t (ξ1:t )]
.
=
N
=

(A.80)
(A.81)
(A.82)
(A.83)
(A.84)

The bound of the KL divergence follows from the elementary inequality dKL (µ||ν) ≤
log(1 + dχ2 (µ||ν)).

A.6. Invariant distribution of the i-cαSMC kernel.
1,N
Lemma A.4. Pθ,y
(dz) is reversible with respect to πθ (dz).

Proof. We mostly suppress dependence on θ since θ is fixed. We will show that the
cαSMC kernel is Gibbs sampler for the artificial joint density given in Eq. (4.7),
which we recall is
f1

1:N
1:N
1
1
1:N
t
π̃1:t (x1:N
1:t , a1:t−1 , f1:t ) , π1:t (x1:t )ψ̃(x1:t , a1:t−1 , f1:t ).

(A.85)

1:N
1
In particular, letting ω1:t , (x1:N
1:t , a1:t−1 , f1:t ), by definition,
f1

t
Pi,N
[dω 1:t ] = π̃1:t (dω 1:t | x̃1:t
= y1:t ).
y 1:i

(A.86)

1:t

Furthermore,
f1

a∗

t
π̃1:t (x̃1:t
= x̃1:tt | ω 1:t ) = π̃1:t (ft1 = a∗t | ω1:t )
∗
a∗
a∗
a∗
a∗
a∗ Q
a∗
s−1
s as−1
, xs s−1 )gs (xs s )αs−1
M1 (x1 1 )g1 (x1 1 ) ts=2 Is Ms (xs−1
∝
∗
∗
a∗ Qt
1:N , x1:N )M (xas−1 , xas−1 )
M1 (x1 1 ) s=2 ra∗s (a∗s−1 |ws−1
s
s s−1
1:s−1

a∗

a∗

= wt t gt (xt t )

1:N
1
∝ Pi,N
[a∗ | x1:N
1:t , a1:t−1 , f1:t ].
y1:i t
1:t

Reversibility now follows easily:
Z
1,N
Py (dz)π(dy) = π̃1:t (dz | ω 1:t )π̃1:t (dω 1:t | y)π(dy)
Z
= π̃1:t (dz | ω 1:t )π̃1:t (dy | ω1:t )π̃1:t (dω 1:t )
Z
= π̃1:t (dy | ω1:t )π̃1:t (dω 1:t | z)π(dz)
= Pz1,N (dy)π(dz).



SEQUENTIAL MONTE CARLO AS APPROXIMATE SAMPLINGS

33

Acknowledgments
The authors would like to thank Arnaud Doucet for critical feedback and numerous helpful suggestions; Cameron Freer for feedback on early versions of this
work; Josh Tenenbaum for discussions that helped to inspire this work; Vikash
Mansinghka for suggesting we investigate the expected value of SMC estimators;
and an anonymous referee for suggestions on improving the presentation of the results. JHH was supported by the U.S. Government under FA9550-11-C-0028 and
awarded by the DoD, Air Force Office of Scientific Research, National Defense Science and Engineering Graduate (NDSEG) Fellowship, 32 CFR 168a. This research
was carried out in part while DMR held a Research Fellowship at Emmanuel College, Cambridge, with funding also from a Newton International Fellowship through
the Royal Society, an NSERC Discovery Grant, Connaught Award, and U.S. Air
Force Office of Scientific Research grant #FA9550-15-1-0074.
References
Christophe Andrieu and G O Roberts. The pseudo-marginal approach for efficient
Monte Carlo computations. The Annals of Statistics, 37(2):697–725, April 2009.
Christophe Andrieu and Matti Vihola. Convergence properties of pseudo-marginal
Markov chain Monte Carlo algorithms. The Annals of Applied Probability, 2014.
Christophe Andrieu, Arnaud Doucet, and Roman Holenstein. Particle Markov
chain Monte Carlo methods. Journal of the Royal Statistical Society: Series B
(Statistical Methodology), 72(3):269–342, 2010.
Christophe Andrieu, Anthony Lee, and Matti Vihola. Uniform Ergodicity of the
Iterated Conditional SMC and Geometric Ergodicity of Particle Gibbs samplers.
arXiv.org, December 2013.
Nicolas Chopin and Sumeetpal S Singh. On the particle Gibbs sampler. arXiv.org,
April 2013.
Julien Cornebise, E Moulines, and Jimmy Olsson. Adaptive methods for sequential importance sampling with application to state space models. Statistics and
Computing, 18(4):461–480, 2008.
Pierre Del Moral. Feynman-Kac Formulae: Genealogical and Interacting Particles
Systems with Applications. Springer, New York, 2004.
Pierre Del Moral, Arnaud Doucet, and Ajay Jasra. Sequential Monte Carlo samplers. Journal of the Royal Statistical Society: Series B (Statistical Methodology),
68(3):411–436, 2006.
Arnaud Doucet and Adam M Johansen. A tutorial on particle filtering and smoothing: fifteen years later. In D Crisan and B Rozovsky, editors, Handbook of Nonlinear Filtering. Cambridge University Press, Cambridge, 2010.
Arnaud Doucet, Simon J Godsill, and Christophe Andrieu. On sequential Monte
Carlo sampling methods for Bayesian filtering. Statistics and Computing, 10(3):
197–208, 2000.
Arnaud Doucet, Nando de Freitas, and Neil Gordon. Sequential Monte Carlo in
Practice. Springer, New York, 2001.
Neil J Gordon, David J Salmond, and Adrian FM Smith. Novel approach to
nonlinear/non-Gaussian Bayesian state estimation. IEE Proceedings F (Radar
and Signal Processing), 140(2):107–113, 1993.
Roman Holenstein. Particle Markov Chain Monte Carlo. PhD thesis, University
Of British Columbia, Vancouver, 2009.

34

J. H. HUGGINS AND D. M. ROY

Nicholas Kantas, Arnaud Doucet, Sumeetpal Sindhu Singh, and Jan Marian Maciejowski. An overview of sequential Monte Carlo methods for parameter estimation in general state-space models. In 15th IFAC Symposium on System
Identification, pages 774–785, 2009.
Hans R Künsch. Particle filters. Bernoulli, 19(4):1391–1403, September 2013.
Anthony Lee and Krzysztof Latuszynski. Variance bounding and geometric ergodicity of Markov chain Monte Carlo kernels for approximate Bayesian computation.
Biometrika, 101(3):655–671, 2014.
Fredrik Lindsten, Randal Douc, and E Moulines. Uniform ergodicity of the Particle
Gibbs sampler. arXiv.org, January 2014.
Nick Whiteley, Anthony Lee, and Kari Heine. On the role of interaction in sequential
Monte Carlo algorithms. Bernoulli, 22(1):494–529, 2016.
Massachusetts Institute of Technology
URL: http://jhhuggins.org/
E-mail address: jhuggins@mit.edu
University of Toronto
URL: http://danroy.org/
E-mail address: droy@utstat.toronto.edu

