Group Importance Sampling for Particle Filtering and
MCMC
arXiv:1704.02771v3 [stat.CO] 28 Apr 2017

Luca Martino? , Vı́ctor Elvira† , Gustau Camps-Valls?
?
†

Image Processing Laboratory, Universitat de València (Spain).

IMT Lille Douai CRISTAL (UMR 9189), Villeneuve d’Ascq (France).

Abstract
Importance Sampling (IS) is a well-known Monte Carlo technique that approximates
integrals involving a posterior distribution by means of weighted samples. In this work,
we study the assignation of a single weighted sample which compresses the information
contained in a population of weighted samples. Part of the theory that we present as Group
Importance Sampling (GIS) has been employed implicitly in different works in the literature.
The provided analysis yields several theoretical and practical consequences. For instance, we
discuss the application of GIS into the Sequential Importance Resampling framework and
show that Independent Multiple Try Metropolis schemes can be interpreted as a standard
Metropolis-Hastings algorithm, following the GIS approach. We also introduce two novel
Markov Chain Monte Carlo (MCMC) techniques based on GIS. The first one, named Group
Metropolis Sampling method, produces a Markov chain of sets of weighted samples. All
these sets are then employed for obtaining a unique global estimator. The second one is the
Distributed Particle Metropolis-Hastings technique, where different parallel particle filters are
jointly used to drive an MCMC algorithm. Different resampled trajectories are compared and
then tested with a proper acceptance probability. The novel schemes are tested in different
numerical experiments such as learning the hyperparameters of Gaussian Processes, the
localization problem in a wireless sensor network and the tracking of vegetation parameters
given satellite observations, where they are compared with several benchmark Monte Carlo
techniques. Three illustrative Matlab demos are also provided.
Keywords: Importance Sampling, Markov Chain Monte Carlo (MCMC), Particle Filtering,
Particle Metropolis-Hastings, Multiple Try Metropolis, Bayesian Inference

1

Introduction

Monte Carlo methods are state-of-the-art tools for approximating complicated integrals involving
multidimensional distributions, which is often required in science and technology [17, 23, 24, 39].
The most popular classes of MC methods are the Importance Sampling (IS) techniques and the
Markov chain Monte Carlo (MCMC) algorithms [23, 39]. IS schemes produce a random discrete

approximation of the posterior distribution by a population of weighted samples [6, 28, 24, 39].
MCMC techniques generate a Markov chain (i.e., a sequence of correlated samples) with a preestablished target probability density function (pdf) as invariant density [23, 24].
In this work, we introduce theory and practice of a novel approach, called Group Importance
Sampling (GIS), where the information contained in different sets of weighted samples is
compressed by using only one, yet properly selected, particle, and one suitable weight. This
general idea supports the validity of different Monte Carlo algorithms in the literature: interacting
parallel particle filters [5, 34, 38], particle island schemes and related techniques [43, 44, 45],
particle filters for model selection [15, 33, 42], nested Sequential Monte Carlo (SMC) methods
[35, 36, 41] are some examples. We point out some consequences of the application of GIS in
Sequential Importance Resampling (SIR) schemes, allowing partial resampling procedures and
the use of different marginal likelihood estimators. Then, we show that the Independent Multiple
Try Metropolis (I-MTM) techniques and the Particle Metropolis-Hastings (PMH) algorithm can
be interpreted as a classical Independent Metropolis-Hastings method by the application of GIS.
Furthermore, we present two novel techniques based on GIS. The first one is the Group
Metropolis Sampling (GMS) algorithm that generates a Markov chain of sets of weighted samples.
All these resulting sets of samples are jointly exploited to obtain a unique particle approximation
of the target distribution. On the one hand, GMS can be considered an MCMC method since it
produces a Markov chain of sets of samples. On the other hand, the GMS can be also considered as
an iterated importance sampler where different estimators are finally combined in order to build a
unique IS estimator. This combination is obtained dynamically through random repetitions given
by MCMC-type acceptance tests. GMS is closely related to Multiple Try Metropolis (MTM)
techniques and Particle Metropolis-Hastings (PMH) algorithms [2, 3, 7, 11, 32, 30], as we discuss
below. The GMS algorithm can be also seen as an extension of the method in [8], for recycling
auxiliary samples in a MCMC method.
The second novel algorithm based on GIS is the Distributed PMH (DPMH) technique where
the outputs of several parallel particle filters are compared by an MH-type acceptance function.
The proper design of DPMH is a direct application of GIS. The benefit of DPMH is twofold:
different type of particle filters (for instance, with different proposal densities) can be jointly
employed, and the computational effort can be distributed in several machines speeding up the
resulting algorithm. As the standard PMH method, DPMH is useful for filtering and smoothing
the estimation of the trajectory of a variable of interest in a state-space model. Furthermore, the
marginal version of DPMH can be used for the joint estimation of dynamic and static parameters.
When the approximation of only one specific moment of the posterior is required, like GMS, the
DPMH output can be expressed as a chain of IS estimators. The novel schemes are tested in three
different numerical experiments: hyperparameter tuning for Gaussian Processes, localization in
a sensor network and filtering of Leaf Area Index (LAI), a parameter widely used to monitor
vegetation from satellite observations. The comparisons with other benchmark Monte Carlo
methods show the benefits of the proposed algorithms.1
The remainder of the paper has the following structure. Section 2 recalls some background
material. The basis of the GIS theory is introduced in Section 3. The applications of GIS in
1

Three illustrative Matlab demos are also provided at https://github.com/lukafree/GIS.git.

2

particle filtering and Multiple Try Metropolis algorithms are discussed in Section 4. In Section 5,
we introduce the novel techniques based on GIS. Section 6.1 provides the numerical results and
in Section 7 we discuss some conclusions.

2

Problem statement and background

In many applications, the goal is to infer a variable of interest, x = x1:D = [x1 , x2 , . . . , xD ] ∈ X ⊆
RD×ξ , where xd ∈ Rξ for all d = 1, . . . , D, given a set of related observations or measurements,
y ∈ RdY . In the Bayesian framework all the statistical information is summarized by the posterior
probability density function (pdf), i.e.,
π̄(x) = p(x|y) =

`(y|x)g(x)
,
Z(y)

(1)

where `(y|x) is the likelihood function, g(x) is the prior pdf and Z(y) is the marginal likelihood
(a.k.a., Bayesian evidence). In general, Z ≡ Z(y) is unknown and difficult to estimate in general,
so we assume to be able to evaluate the unnormalized target function,
π(x) = `(y|x)g(x).

(2)

The computation of integrals involving π̄(x) = Z1 π(x) is often intractable. We consider the Monte
Carlo approximation of these complicated integrals involving the target π̄(x) and an integrable
function h(x) with respect to π̄, i.e.,
Z
h(x)π̄(x)dx,
(3)
I = Eπ̄ [h(X)] =
X

where we denote X ∼ π̄(x). The basic Monte Carlo (MC) procedure consists in drawing N
P
independent samples from the target pdf, i.e., x1 , . . . , xN ∼ π̄(x), so that IbN = N1 N
n=1 h(xn ) is
an estimator of I [24, 39]. However, in general, direct methods for drawing samples from π̄(x) do
not exist so that alternative procedures are required. Below, we describe perhaps one of the most
popular.

2.1

Importance Sampling

Let us consider the use of a simpler proposal pdf, q(x), and rewrite the integral I in Eq. (3) as
Z
1
π(x)
I = Eπ̄ [h(X)] = Eq [h(X)w(X)] =
h(x)
q(x)dx,
(4)
Z X
q(x)
where w(x) = π(x)
: X → R. This suggests an alternative procedure. Indeed, we can draw N
q(x)
samples x1 , . . . , xN from q(x),2 and then assign to each sample the unnormalized weights
wn = w(xn ) =
2

π(xn )
,
q(xn )

n = 1, . . . , N.

We assume that q(x) > 0 for all x where π̄(x) 6= 0, and q(x) has heavier tails than π̄(x).

3

(5)

PN
1
If Z is known, an unbiased IS estimator [24, 39] is defined as IbN = ZN
n=1 wn h(xn ), where
xn ∼ q(x). If Z is unknown, defining the normalized weights, w̄n = PNwn w with n = 1, . . . , N , an
i=1 i
alternative asymptotically unbiased IS estimator is given by
IN =

N
X

w̄n h(xn ).

(6)

n=1

Both IbN and I N are consistent estimators of I in Eq. (3) [24, 39]. Moreover, an unbiased
R
P
estimator of marginal likelihood, Z = X π(x)dx, is given by Zb = N1 N
i=1 wi . More generally, the
N
pairs {xi , wi }i=1 can be used to build a particle approximation of the posterior distribution,
N
N
X
1 X
π
b(x|x1:N ) =
wn δ(x − xn ) =
w̄n δ(x − xn ),
N Zb
n=1

(7)

n=1

where δ(x) denotes the Dirac delta function. Table 1 summarizes the main notation of the work.
Note that the words sample and particle are used as synonyms along this work. Moreover, Table
14 shows the main used acronyms.

2.2

Concept of proper weighting

The standard IS weights in Eq. (5) are broadly used in the literature. However the definition of
properly weighted sample can be extended as suggested in [39, Section 14.2], [24, Section 2.5.4]
and in [16]. More specifically, given a set of samples, they are properly weighted with respect to
the target π̄ if, for any integrable function h,
Eq [w(xn )h(xn )] = cEπ̄ [h(xn )],

∀n ∈ {1, . . . , N },

(8)

where c > 0 is a constant value, also independent from the index n, and the expectation of the
left hand side is performed, in general, w.r.t. to the joint pdf of w(x) and x, i.e., q(w, x). Namely,
the weight w(x), conditioned to a given value of x, could even be considered a random variable.
Thus, in order to obtain consistent estimators, one can design any joint q(w, x) as long as the
restriction of Eq. (8) is fulfilled.

3

Group Importance Sampling: weighting a set of samples

In this section, we use the general definition in Eq. (8) for designing proper weights and summary
samples assigned to different sets of samples. Let us consider M sets of weighted samples,
NM
N2
1
S1 = {x1,n , w1,n }N
n=1 , S2 = {x2,n , w2,n }n=1 , ...., SM = {xM,n , wM,n }n=1 , where xm,n ∼ qm (x),
i.e., a different proposal pdf for each set Sm and in general Ni 6= Nj , for all i 6= j, i, j ∈ {1, ..., M }.
In some applications and different Monte Carlo schemes, it is convenient (and often required) to
em ,
compress the statistical information contained in each set using a pair of summary sample, x
and summary weight, Wm , m = 1, . . . , M , in such a way that the following expression
1
IeM = PM

j=1

Wj
4

M
X

m=1

Wm h(e
xm ),

(9)

Table 1: Main notation of the work.

x = [x1 , . . . , xD ]
π̄(x)
π(x)
π
b(x|x1:N )
e
x
wn = w(xn )
w̄n = w̄(xn )
w
em = w(e
e xm )
Wm
IN
IeN
Z
b Z
Z,

Variable of interest, x ∈ X ⊆ RD×ξ , with xd ∈ Rξ for all d
Normalized posterior pdf, π̄(x) = p(x|y)
Unnormalized posterior function, π(x) ∝ π̄(x)
Particle approximation of π̄(x) using the set of samples x1:N = {xn }N
n=1
e∼π
e ∈ {x1 , . . . , xN })
Resampled particle, x
b(x|x1:N ) (note that x
Unnormalized standard IS weight of the particle xn
Normalized weight associated to xn
em
Unnormalized proper weight associated to the resampled particle x
Summary weight of m-th set Sm
Standard self-normalized IS estimator using N samples
Self-normalized estimator using N samples and based on GIS theory
Marginal likelihood; normalizing constant of π(x)
Estimators of the marginal likelihood Z

is still a consistent estimator of I, for a generic integrable function h(x). Thus, although the
compression is lossy, we still have a suitable particle approximation π
b of the target π̄ as shown
below. In the following, we denote the importance weight of the n-th sample in the m-th group
m,n )
, the m-th marginal likelihood estimator as
as wm,n = w(xm,n ) = qπ(x
m (xm,n )
Nm
1 X
Zbm =
wm,n ,
Nm n=1

and the normalized weights within a set, w̄m,n =
m = 1, . . . , M .

wm,n
P Nm
j=1 wm,j

(10)
=

wm,n
bm ,
Nm Z

for n = 1, . . . , Nm and

Definition 1. A resampled particle, i.e.,
em ∼ π
x
bm (x) = π
b(x|xm,1:Nm ) =

Nm
X
n=1

w̄m,n δ(x − xm,n ),

(11)

em for the m-group. Note that x
em is selected within {xm,1 , . . . , xm,Nm }
is a summary particle x
according to the probability mass function (pmf ) defined by w̄m,n , n = 1, . . . , Nm .
It is possible to use the Liu’s definition in order to assign a proper importance weight to a
resampled particle [26], as stated in the following theorem.
em ∼ π
Theorem 1. Let us consider a resampled particle x
bm (x) = π
b(x|xm,1:Nm ). A proper
unnormalized weight following the Liu’s definition in Eq. (8) for this resampled particle is
w
em = Zbm , defined in Eq. (10).

e0m , x
e00m , resampled with
The proof is given in Appendix A. Note that two (or more) particles, x
5

e0m , x
e00m ∼ π
replacement from the same set and hence from the same approximation, x
bm (x), have the
00
0
b
e xm ) = Zm , as depicted in Figure 1. Note that the classical importance
same weight w(e
e xm ) = w(e
weight cannot be computed for a resampled particle, as explained in Appendix A and pointed out
in [22, 26, 35], [28, App. C1].
wm,4
wm,2
wm,1

Resample 3 times

wm,3
xm,1 xm,2

xm,3

x

e000
x
m

xm,4

bm
Z

bm
Z
= xm,2

e0m
x

=

e00m
x

x
= xm,4

Figure 1: Example of generation (one run) and proper weighting of 3 resampled
particles (with
1 P4
0
00
000
b
em , x
em and x
em , from the m-th group, where Nm = 4 and Zm = 4 n=1 wm,n .
replacement), x

Definition 2. The summary weight for the m-th group of samples is Wm = Nm w
em = Nm Zbm .

Particle approximation. Figure 2 represents graphically an example of GIS with M = 2 and
N1 = 4, N2 = 3. Given the M summary pairs {e
xm , w
em }M
m=1 in a common computational node,
we can obtain the following particle approximation of π̄(x), i.e.,
π
b(x|e
x1:M ) = PM

1

b
j=1 Nj Zj

M
X

m=1

em ),
Nm Zbm δ(x − x

involving M weighted samples in this case (see App.
corresponding specific GIS estimator in Eq. (9) is

(12)

B). For a given function h(x), the

M
X
1
e
Nm Zbm h(e
xm ).
IM = PM
b
N
Z
j=1 j j m=1

(13)

It is a consistent estimator of I, as we show in Appendix B. The expression in Eq. (13) can
be interpreted as a standard IS estimator where w(e
e xm ) = Zbm is a proper weight of a resampled
particle [26], and we give more importance to the resampled particles belonging to a set with
higher cardinality. See DEMO-2 at https://github.com/lukafree/GIS.git.
Combination of estimators. If we are only interested in computing the integral I for a specific
(m)
function h(x), we can summarize the statistical information by the pairs {I Nm , w
em } where
(m)
I Nm

=

Nm
X

w̄m,n h(xm,n ),

n=1

6

(14)

N1 Zb1
N1 = 4

w1,4

w1,2

N2 Zb2

e1
x

N2 = 3
w2,1

w1,1

e2
x

w2,3

w2,2

w1,3

x1,1 x1,2 x1,3

x1,4 x

x2,1

x2,2 x2,3

x

Figure 2: Graphical representation of GIS. In this case, M = 2 groups of N1 = 4 and N2 = 3 weighted
samples are summarized with a resampled particle and one summary weight w
em = Nm Zbm , m = 1, 2.

P
is the m-th partial IS estimator obtained by using Nm samples in Sm . Given all the S = M
j=1 Nj
weighted samples in the M sets, the complete estimator I S in Eq. (6) can be written as a convex
(m)
combination of the M partial IS estimators, I Nm , i.e.,
I S = PM

1

j=1

= PM

1

j=1

= PM

Nj Zbj

1

j=1

Nj Zbj
Wm

M X
Nm
X

wm,n h(xm,n ),

(15)

m=1 n=1
M
X

m=1

M
X

Nm Zbm

Nm
X

w̄m,n h(xm,n ),

(16)

n=1

(m)

Wm I Nm .

(17)

m=1

The equation above shows that the summary weight Wm measures the importance of the m-th
(m)
estimator I Nm . This confirms that Wm is a proper weight the group of samples Sm , and also
suggests another valid compression scheme.
Remark 1. In order to approximate only one specific moment I of π̄(x), we can summarize the
(m)
m-group with the pair {I Nm , Wm }M
m=1 , thus all the M partial estimators can be combined following
Eq. (17).
In this case, there is no loss of information w.r.t. storing all the weighted samples. However,
the approximation of other moments of π̄(x) is not possible. Figures 3-4 depict the graphical
representations of the two possible approaches for GIS.

4
4.1

GIS in other Monte Carlo schemes
Application in particle filtering

In Section 2.1, we have described the IS procedure in a batch0 way, i.e., generating directly a D)
dimensional vector x0 ∼ q(x) and then compute the weight π(x
. This procedure can be performed
q(x0 )
7

(a)

(b)

Central Node

Central Node
{e
x1 , W1 }

Wm

W1

S1

e1
x

...

{e
xm , Wm }

em
x

Sm

...

(1)

{I N1 , W1 }

{e
xM , WM }

Wm

W1

WM

(m)
I Nm

(1)
I N1

eM
x

S1

SM

(M )

(m)
{I Nm , Wm }

...

Sm

{I NM , WM }

WM
(M )

I NM
...

SM

Figure 3: Graphical overview of GIS in a parallel/distributed framework. (a) The central node obtains
em ∈ RD , and
all the pairs {e
xm , Wm }M
b(x|e
x1:M ) or I M . Note that only M particles, x
m=1 , and provides π
P
M scalar weights, Wm ∈ R, are transmitted, instead of S samples and S weights, with S = M
m=1 Nm .
(b) Alternatively, if we are interested only in a specific moment of the target, we can transmit the pairs
(m)
{I Nm , Wm }M
m=1 and then combine them as in Eq. (17).

Figure 4: Possible outputs of different GIS compression schemes. On the left, {Sm }M
m=1 , no compression
(m)

is applied. In the center, {I Nm , Wm }M
m=1 , we can perfectly reconstruct the estimator I S in Eq. (15) where
PM
(m)
S = m=1 Nm , but we cannot approximate other moments. Using {e
xNm , Wm }M
m=1 , we always obtain a
lossy compression, but any moments of π̄(x) can be approximated, as shown in Eqs. (12)-(13).

8

sequentially if the target density is factorized. In this case, the method is known as Sequential
Importance Sampling (SIS). It is the basis of particle filtering, along with the use of the resampling
procedure. Below, we describe the SIS method.
4.1.1

Sequential Importance Sampling (SIS)

Let us that recall x = x1:D = [x1 , x2 , . . . , xD ] ∈ X ⊆ RD×ξ where xd ∈ Rξ for all d = 1, . . . , D, and
let us consider a target pdf π̄(x) factorized as
D
Y
1
1
γ1 (x1 )
π̄(x) = π(x) =
γd (xd |x1:d−1 ),
Z
Z
d=2

(18)

where γ1 (x1 ) is a marginal pdf and γd (xd |x1:d−1 ) are conditional
pdfs. We can also consider a
Q
proposal pdf decomposed in the same way, q(x) = q1 (x1 ) D
q
(x
d=2 d d |xd−1 ). In a batch IS scheme,
(n)
given the n-th sample xn = x1:D ∼ q(x), we assign the importance weight
(n)

w(xn ) =
(n)

where β1

(n)

=

π(x1 )

(n)

(n)

(n)

(n)

D

γ1 (x1 )γ2 (x2 |x1 ) · · · γD (xD |x1:D−1 ) Y
π(xn )
βd ,
=
=
(n)
(n) (n)
(n) (n)
q(xn )
q1 (x1 )q2 (x2 |x1 ) · · · qD (xD |x1:D−1 )
d=1
(n)

(n)

=

(n)

γd (xd |x1:d−1 )

, with d = 1, . . . , D. Let us also denote the
Q
joint probability of [x1 , . . . , xd ] as π̄d (x1:d ) = Z1d πd (x1:d ) = Z1d γ1 (x1 ) dj=2 γj (xj |x1:j−1 ), so that
π̄D (x1:D ) ≡ π̄(x) and ZD ≡ Z. Thus, we can draw samples generating sequentially each component
Q
(n)
(n)
(n)
xd ∼ qd (xd |x1:d−1 ), d = 1, . . . , D, so that xn = x1:D ∼ q(x) = q1 (x1 ) D
d=2 qd (xd |xd−1 ), and
compute recursively the corresponding IS weight as in Eq. (19). The SIS technique is also given
in Table 2 by setting η = 0.
(n)

q(x1 )

and βd

(19)

(n)

(n)

qd (xd |x1:d−1 )

Remark
R 2. In SIS, there are two possible formulations of the estimator of the marginal likelihoods
Zd = Rdη πd (x1:d )dx1:d ,
N
N
1 X (n)
1 X (n) (n)
wd =
w β ,
N n=1
N n=1 d−1 d
" N
#
d
Y
X (n) (n)
=
w̄j−1 βj
.

Zbd =

Zd

j=1

(20)
(21)

n=1

In SIS, both estimators are equivalent Z d ≡ Zbd . See Appendix C for further details.
4.1.2

Sequential Importance Resampling (SIR)

The expression in Eq. (19) suggests a recursive procedure for generating the samples and
computing the importance weights, as shown in Steps 2a and 2b of Table 2. In Sequential
Importance Resampling (SIR), a.k.a., standard particle filtering, resampling steps are incorporated
9

during the recursion as in step 2(c)ii of Table 2 [12, 13]. In general, the resampling steps are
applied only in certain iterations in order to avoid the path degeneration, taking into account
[ of the Effective Sampling Size (ESS) [20, 27]. If ESS
[ is smaller than a
an approximation ESS
pre-established threshold, the particles are resampled. Two examples of ESS approximation are
(n)
w
(n)
1
[=
[ = PN 1 (n) and ESS
= PN d (i) . Note that, in both cases,
ESS
(n) where w̄d
n=1 (w̄d

)2

max w̄d

i=1

[ ≤ N.
0 ≤ ESS

wd

(22)

[ < ηN where η ∈ [0, 1].
Hence, the condition for the adaptive resampling can be expressed as ESS
SIS is given when η = 0 and SIR for η ∈ (0, 1]. When η = 1, the resampling is applied at
each iteration and in this case SIR is often called bootstrap particle filter [12, 13, 14]. If η = 0, no
resampling is applied, we only apply Steps 2a and 2b, and we have the SIS method described above,
(n)
that after D iterations is completely equivalent to the bath IS approach, since wn = w(xn ) ≡ wD
where xn = x1:D .
Partial resampling. In Table 2, we have considered that only a subset of R ≤ N particles are
resampled. In this case, step 2(c)iii including the GIS weighting is strictly required in order to
provide final proper weighted samples and hence consistent estimators. The partial resampling
procedure is an alternative approach to prevent the loss of particle diversity [26]. In the classical
description of SIR [40], we have R = N (i.e., all the particles are resampled) and the weight
recursion follows setting the unnormalized weights of the resampled particles to any equal value.
Since all the N particles have been resampled, the selection of this value has no impact in the
weight recursion and in the estimation of I.
Marginal likelihood estimators. Even in the case R = N , i.e., all the particle are resampled
as in the standard SIR method, without using the GIS weighting only the formulation Z d in Eq.
(n)
(21) provides a consistent estimator of Zd , since it involves the normalized weights w̄d−1 , instead
(n)
of the unnormalized ones, wd−1 .
Remark 3. If the GIS weighting is applied in SIR, both formulations Zbd and Z d in Eqs. (20)-(21)
provide consistent estimator of Zd and they are equivalent, Zbd ≡ Z d (as in SIS). See an exhaustive
discussion in Appendix C.
See DEMO-1 at https://github.com/lukafree/GIS.git.
GIS in Sequential Monte Carlo (SMC). The idea of summary sample and summary weight
have been implicitly used in different SMC schemes proposed in literature, for instance, for the
communication among parallel particle filters [5, 34, 38], and in the particle island methods
[43, 44, 45]. GIS also appears indirectly in particle filtering for model selection [15, 33, 42]
and in the so-called Nested Sequential Monte Carlo techniques [35, 36, 41].

10

Table 2: SIR with partial resampling

1. Choose N the number of particles, R ≤ N the number of particles to be resampled,
(n)
[ [27] and a constant
the initial particles x0 , n = 1, . . . , N , an ESS approximation ESS
value η ∈ [0, 1].
2. For d = 1, . . . , D:
(n)

(n)

(a) Propagation: Draw xd ∼ qd (xd |x1:d−1 ), for n = 1, . . . , N .

(b) Weighting: Compute the weights
(n)
wd

=

(n) (n)
wd−1 βd

d
Y

=

(n)

βj ,

n = 1, . . . , N,

(23)

j=1

(n)

(n)

where βd =

(n)

γd (xd |x1:d−1 )
(n) (n)
qd (xd |x1:d−1 )

.

[ < ηN then:
(c) if ESS
(j )

i. Select randomly a set of particles S = {xd r }R
r=1 where R ≤ N , jr ∈
{1, . . . , N } for all r, and jr 6= jk for r 6= k.
ii. Resampling: Resample R times within the set S according to the
(j )

probabilities w̄d r =

for r = 1, . . . , R.

(j )

wd r
PR
(jk ) ,
k=1 wd

iii. GIS weighting: Compute ZbS =
r = 1, . . . , R.
(n)

(j )

(j )

(j )

r
obtaining {x̄d r }R
= x̄d r ,
r=1 . Then, set xd

1
R

R
P

r=1

(jr )
(j )
wd−1
and set wd r = ZbS for all

(n)

3. Return {xn = x1:D , wn = wD }N
n=1 .

4.2

Multiple Try Metropolis schemes as a Standard MetropolisHastings method

The Metropolis-Hastings (MH) method is one of the most popular MCMC algorithm [23, 24, 39].
It generates a Markov chain {xt }∞
t=1 where π̄(x) is the invariant density. Considering a proposal
pdf q(x) independent from the previous state xt−1 , the corresponding Independent MH (IMH)
scheme is formed by the steps in Tableh 3 [39]. i
w(v0 )
Observe that α(xt−1 , v0 ) = min 1, w(x
in Eq. (24) involves the ratio between the
t−1 )
0
importance weight of the proposed sample v at the t-th iteration, and the importance weight
of the previous state xt−1 . Furthermore, note that at each iteration only one new sample v0 is

11

Table 3: The Independent Metropolis-Hastings (IMH) algorithm

1. Choose an initial state x0 .
2. For t = 1, . . . , T :
(a) Draw a sample v0 ∼ q(x).

(b) Accept the new state, xt = v0 , with probability




π(v0 )q(xt−1 )
w(v0 )
0
α(xt−1 , v ) = min 1,
= min 1,
,
π(xt−1 )q(v0 )
w(xt−1 )
where w(x) =

π(x)
q(x)

(24)

(standard importance weight). Otherwise, set xt = xt−1 .

3. Return {xt }Tt=1 .
generated and compared with the previous state xt−1 by the acceptance probability α(xt−1 , v0 )
(in order to obtain the next state xt ). The Particle Metropolis-Hastings (PMH) method3 [2]
and the alternative version of the Independent Multiply Try Metropolis technique [29] (denoted
as I-MTM2) are jointly described in Table 4. They are two MCMC algorithms where at each
iteration several candidates {v1 , . . . , vN } are generated. After computing the IS weights w(vn ), one
candidate is selected vj within the N possible values, i.e., j ∈ {1, . . . , N }, applying a resampling
n)
n)
= w(v
step according to the probability mass w̄n = PNw(vw(v
b0 , n = 1, . . . , N . Then the selected
NZ
i)
i=1
sample vj is tested with a proper probability α(xt−1 , vj ) in Eq. (25).
The difference between PMH and I-MTM2 is the procedure employed for the generation of
the N candidates and for the construction of the weights. PMH employs a sequential approach,
whereas I-MTM2 uses a standard batch approach [29]. Namely, PMH generates sequentially the
components vj,k of the candidates, vj = [vj,1 , . . . , vj,D ]> , and compute recursively the weights
as shown in Section 4.1. Since resampling steps are often used, then the resulting candidates
v1 , . . . , vN are correlated, whereas in I-MTM2 they are independent. I-MTM2 coincides with
PMH if the candidates are generated sequentially but without applying resampling steps, so that
I-MTM2 can be considered a special case of PMH.
Note that w(v
e j ) = Zb0 and w(x
e t−1 ) = Zbt−1 are the GIS weights of the resampled particles
vj and xt−1 respectively, as stated in Definition 1 and Theorem 1.4 Hence, considering the GIS
3
PMH is used for filtering and smoothing a variable of interest in state-space models (see, for instance, Figure
11). The Particle Marginal MH (PMMH) algorithm [2] is an extension of PMH employed in order to infer both
dynamic and static variables. PMMH is described in Appendix D.
N w(v
e j)
w(v
e j)
4
t
Note that the number of candidates per iteration is constant (N ), so that WWt−1
= N w(x
e t−1 ) = w(x
e t−1 ) .

12

Table 4: PMH and I-MTM2 techniques

1. Choose an initial state x0 and Zb0 .
2. For t = 1, . . . , T :

(a) Draw N particles v1 , . . . , vN from q(x) and weight them with the proper
importance weight w(vn ), n = 1, . . . , N , using a sequential approach (PMH),
P
or a batch approach (I-MTM2). Thus, denoting Zb0 = N1 N
n=1 w(vn ), we obtain
P
N
1
the particle approximation π
b(x|v1:N ) = N Zb0 n=1 w(vn )δ(x − vn ).

(b) Draw vj ∼ π
b(x|v1:N ).

(c) Set xt = vj and Zbt = Zb0 , with probability

"

#
Zb0
α(xt−1 , vj ) = min 1,
.
Zbt−1

(25)

Otherwise, set xt = xt−1 and Zbt = Zbt−1 .

3. Return {xt }Tt=1 .
theory, we can write

"

Zb0
α(xt−1 , vj ) = min 1,
Zbt−1

#




w(v
e j)
= min 1,
,
w(x
e t−1 )

(26)

which has the form of the acceptance function of the classical IMH method in Table 3. Therefore,
PMH and I-MTM2 algorithms can be also summarized as in Table 5.
Remark 4. The PMH and I-MTM2 algorithms take the form of the classical IMH method
employing the equivalent proposal pdf qe(x) in Eq. (27) (depicted in Figure 5; see also Appendix A),
e0 ∼ qe(x), within the acceptance function
and using the GIS weight w(e
e x0 ) of a resampled particle x
0
e ).
α(xt , x

5

Novel MCMC techniques based on GIS

In this section, we provide two examples of novel MCMC algorithms based on GIS. First of all,
we introduce a Metropolis-type method producing a chain of sets of weighted samples. Secondly,
we present a PMH technique driven by M parallel particle filters. In the first scheme, we exploit
the concept of summary weight and all the weighted samples are stored. In the second one, both
concepts of summary weight and summary particle are used. The consistency of the resulting
estimators and the ergodicity of both schemes is ensured.
13

Table 5: Alternative description of PMH and I-MTM2

1. Choose an initial state x0 .
2. For t = 1, . . . , T :
e0 ∼ qe(x), where
(a) Draw x

qe(x) =

Z

XN

"

N
Y
i=1

#

q(vi ) π
b(x|v1:N )dv1:N ,

(27)

is the equivalent proposal pdf associated to a resampled particle [22, 26].
e0 , with probability
(b) Set xt = x



w(e
e x0 )
e ) = min 1,
α(xt−1 , x
.
w(x
e t−1 )
0

(28)

Otherwise, set xt = xt−1 .
3. Return {xt }Tt=1 .

5.1

Group Metropolis Sampling

Here, we describe an MCMC procedure that yields a sequence of sets of weighted samples. All
the samples are then employed for a joint particle approximation of the target distribution. The
Group Metropolis Sampling (GMS) is outlined in Table 6. Figures 6(a)-(b) give two graphical
representations of GMS outputs (with N = 4 in both cases). Note that the GMS algorithm uses the
idea of summary weight for comparing sets. Let us denote as ρn,t the importance weights assigned
to the samples xn,t contained in the current set St . Given the generated sets St = {xn,t , ρn,t }N
n=1 ,
for t = 1, . . . , T , GMS provides the global particle approximation
T
N
1 X X ρn,t
π
b(x|x1:N,1:T ) =
δ(x − xn,t ),
P
T t=1 n=1 N
i=1 ρi,t
N
T
1 XX
ρ̄n,t δ(x − xn,t ).
=
T t=1 n=1

(29)
(30)

Thus, the estimator of a specific moment of the target is

T
N
T
1 XX
1 X e(t)
e
IN T =
ρ̄n,t h(xn,t ) =
I .
T t=1 n=1
T t=1 N

(31)

If the N candidates, v1 , . . . , vN , and the associated weights, w1 , . . . , wN , are built sequentially
by a particle filtering method, we have a Particle GMS (PGMS) algorithm (see Section 6.3) and
14

Figure 5: (Left) Graphical representation of the generation of one sample x0 from the equivalent
proposal pdf qe(x) in Eq. (27).(Right) Example of the equivalent density qe(x) (solid line) with
N = 2. The target, π̄(x), and proposal, q(x), pdfs are shown with dashed lines. See DEMO-3 at
https://github.com/lukafree/GIS.git.

marginal versions can be also considered (see Appendix D).
Table 6: Group Metropolis Sampling

b
1. Build an initial set S0 = {xn,0 , ρn,0 }N
n=1 and Z0 =

2. For t = 1, . . . , T :

1
N

PN

n=1

ρn,0 .

(a) Draw N samples, v1 , . . . , vN ∼ q(x) following a sequential or a batch procedure.
(b) Compute the weights, wn =
P
compute Zb0 = N1 N
n=1 wn .

π(vn )
,
q(vn )

n = 1, . . . , N ; define S 0 = {vn , wn }N
n=1 ; and

0
b
b0
(c) Set St = {xn,t = vn , ρn,t = wn }N
n=1 , i.e., St = S , and Zt = Z , with probability
"
#
b0
Z
α(St−1 , S 0 ) = min 1,
.
(32)
Zbt−1

Otherwise, set St = St−1 and Zbt = Zbt−1 .

N
P
(t)
(t)
3. Return {St }Tt=1 , or {IeN }Tt=1 where IeN =

n=1

ρn,t
PN
h(xn,t ).
i=1 ρi,t

Relation with IMH. The acceptance probability α in Eq. (32) is the extension of the
acceptance probability of IMH in Eq. (24), by considering the proper GIS weighting of a set
of weighted samples. Note that, in this version of GMS, all the sets contain the same number of
15

samples.
Relation with MTM methods. GMS is strictly related to Multiple Try Metropolis (MTM)
schemes [7, 31, 32, 29] and Particle Metropolis Hastings (PMH) techniques [2, 29]. The main
difference is that GMS uses no resampling steps at each iteration for generating summary samples,
indeed GMS uses the entire set. However, considering a sequential of a batch procedure for
generating the N tries at each iteration, we can recover an MTM (or the PMH) chain by the GMS
output applying one resampling step when St 6= St−1 ,

N
X

v
et ∼
ρ̄n,t δ(x − xn,t ), if St 6= St−1 ,
et =
x
(33)
n=1


et−1 ,
x
if St = St−1 ,
for t = 1, . . . , T . Namely, {e
xt }Tt=1 is the chain obtained by one run of the MTM (or PMH)
technique. Figure 6(b) provides a graphical representation of a MTM chain recovered by GMS
outputs.

Ergodicity. As also discussed above, (a) the sample generation, (b) the acceptance probability
function and hence (c) the dynamics of GMS exactly coincide with the corresponding steps of PMH
or MTM (with a sequential or batch particle generation, respectively). Hence, the ergodicity of
the chain is ensured [7, 32, 2, 29]. Indeed, we can recover the MTM (or PMH) chain as shown in
Eq. (33).
Recycling samples. The GMS algorithm can be seen as a method of recycling auxiliary
weighted samples in MTM schemes (or PMH schemes, if the candidates are generated by SIR). In
[8], the authors show how recycling and including the samples rejected in one run of a standard
MH method into a unique consistent estimator. GMS can be considered an extension of this
technique where N ≥ 1 candidates are drawn at each iteration.
Iterated IS. GMS can be also interpreted as an iterative importance sampling scheme where
an IS approximation of N samples is built at each iteration and compared with the previous
(t)
IS approximation. This procedure is iterated T times and all the accepted IS estimators IeN
are finally combined to provide a unique global approximation of N T samples. Note that the
temporal combination of the IS estimators is obtained dynamically by the random repetitions
due to the rejections in the MH test. Hence, the complete procedure for weighting the samples
generated by GMS can be interpreted as the composition of two weighting schemes: (a) by an
importance sampling approach building {ρn,t }N
n=1 and (b) by the possible random repetitions due
to the rejections in the MH test.
Consistency of the GMS estimator. Recovering the MTM chain {e
xt }Tt=1 as in Eq. (33),
P
T
1
xt ) is consistent, i.e., IeT converges almost-surely to I = Eπ̄ [h(x)] as
the estimator IeT = T t=1 h(e
P
e(t)
T → ∞, since {e
xt }Tt=1 is an ergodic chain [39]. Note that Eπb [h(e
xt )|St ] = N
n=1 ρ̄n,t h(xn,t ) = IN
16

PN
for St 6= St−1 , where π
b(x|x1:N,t ) =
xt )|St ] =
b [h(e
n=1 ρ̄n,t δ(x − xn,t ). If St = St−1 , then Eπ
(t−1)
(t)
(t−1)
(t)
e
e
e
e
Eπb [h(e
xt−1 )|St−1 ] = IN
and, since IN = IN , we have again Eπb [h(e
xt )|St ] = IN . Therefore, we
have
E[IeT |S1:T ] =

T
1X
Eπb [h(e
xt )|St ]
T t=1

(34)

T
1 X e(t) e
=
I = IN T .
T t=1 N

(35)

Thus, the GSM estimator IeN T in Eq. (31) can be expressed as IeN T = E[IeT S1:T ], where S1:T are
all the weighted samples obtained by GMS, and it is consistent for T → ∞ since IeT is consistent.
Furthermore, fixing T , IeN T is consistent for N → ∞ for standard IS arguments [24].
(b)

(a)

Figure 6: (a) Chain of sets St = {xn,t , ρn,t }N
n=1 generated by the GMS method (graphical representation
with N = 4). (b) Graphical examples of GMS outputs, St , St+1 , St+2 and St+3 , where St+2 = St+1 . The
weights of the samples are denoted by the size of the circles. A possible recovered MTM chain is also
eτ with τ = t, t + 1, t + 2, t + 3 and x
et+2 = x
et+1 .
depicted with solid line, where the states are x

5.2

Distributed Particle Metropolis-Hastings algorithm

The PMH algorithm is an MCMC technique particularly designed for filtering and smoothing a
dynamic variable in a state-space model [2, 29] (see for instance Figure 11). In PMH, different
trajectories obtained by different runs of a particle filter (see Section 4.1) are compared according
to suitable MH-type acceptance probabilities, as shown in Table 4. In this section, we show how
several parallel particle filters (for instance, each one consider a different proposal pdf) can drive
a PMH-type technique.
Q
The classical PMH method uses a single proposal pdf q(x) = q1 (x1 ) D
d=2 qd (xd |x1:d−1 ), employed
in single SIR method in order to generate new candidates before of the MH-type test (see Table
4). Let us consider the problem
of tracking a variable of interest x = [x1 , . . . , xD ]> ∈ RD×ξ with
QD
target pdf π(x) = π1 (x1 ) d=2 πd (xd |x1:d−1 ). We assume that M independent processing units are
available jointly with a central node as shown Fig. 7. We use M parallel particle filters, each one
17

Q
with a different proposal pdf, qm (x) = qm,1 (x1 ) D
d=2 qm,d (xd |x1:d−1 ), one per each processor. Then,
after one run of the parallel particle filters, we obtain M particle approximations π
bm (x). Since,
we aim to reduce the communication cost to the central node (see Figs. 3 and 7), we consider
em }, where x
em ∼ π
that each machine only transmits the pair {Zbm , x
bm (x) (we set N1 = . . . = NM ,
for simplicity). Applying the GIS theory, then it is straightforward to outline the method, called
Distributed Particle Metropolis-Hastings (DPMH) technique, shown in Table 7.
Table 7: Distributed Particle Metropolis-Hastings algorithm

1. Choose an initial state x0 and Zbm,0 for m = 1, . . . , M .
2. For t = 1, . . . , T :

(a) (Parallel Processors) Draw N particles vm,1 , . . . , vm,N from qm (x) and
weight them with IS weights w(vm,n ), n = 1, . . . , N , using a particle filter
(orPa batch approach), for each m = 1, . . . , M . Thus, denoting Zbm =
N
1
w(vm,n ), we obtain the M particle approximations π
bm (x) = π
b(x|vm,1:N ) =
n=1
N
PN
1
b
n=1 w(vm,n )δ(x − vm,n ).
NZ
m

em ∼ π
(b) (Parallel Processors) Draw x
b(x|vm,1:N ), for m = 1, . . . , M .

e ∈ {e
eM } according to the pmf
(c) (Central Node) Resample x
x1 , . . . , x

b
Z
PM m b ,
j=1 Zj

e∼π
eM ).
m = 1, . . . , M , i.e., x
b(x|e
x1 , . . . , x
e and Zbm,t = Zbm , for m = 1, . . . , M , with probability
(d) (Central Node) Set xt = x
"
PM b #
Zm
e) = min 1, PMm=1
.
(36)
α(xt−1 , x
bm,t−1
Z
m=1
Otherwise, set xt = xt−1 and Zbm,t = Zbm,t−1 , for m = 1, . . . , M .

The method in Table 7 has the structure of a Multiple Try Metropolis (MTM) algorithm using
different proposal pdfs [7, 32]. More generally, in step 2a, the scheme described above can even
employ different kinds of particle filtering algorithms. In step 2b, M total resampling steps are
performed,one per processor. Then, one resampling step is performed in the central node (step
2c). Finally, the resampled particle is accepted as new state with probability α in Eq. (36).

18

Ergodicity. The ergodicity of DPMH is ensured since it can be interpreted as a standard PMH
method considering a single particle approximation5
π
b(x|v1:M,1:N ) =

M
X

Zbm

PM b
m=1
j=1 Zj

π
b(x|vm,1:N ) =

M
X

m=1

W mπ
b(x|vm,1:N ),

(37)

e∼π
and then we resample once, i.e., draw x
b(x|v1:M,1:N ). Then, the proper weight of this resampled
P
M
b
particle is Zb = m=1
function
of the equivalent classical PMH method
i
h Zm , so
i that the
h acceptance
1 PM
P
bm
Z
b
e) = min 1, Z = min 1, 1 MPMm=1
, where Zbt−1 = M Zbm,t−1 (see Table 4).
is α(xt−1 , x
bt−1
Z

M

m=1

bm,t−1
Z

m=1

Using partial IS estimators. if we are interested in approximating only one moment of the
target pdf, as shown in Figures 3-4, at each iteration we can transmit the M partial estimators
PM b (m)
(m)
I N and combine them in the central node as in Eq. (17), obtaining IeN0 M = PM1 Zb
m=1 Zm I N .
j=1

j

(t)
Then, a sequence of estimators, IeN M , is created according to the acceptance probability α in Eq.
(36). Finally, we obtain the global estimator
T
1 X e(t)
e
I .
IN M T =
T t=1 N M

(38)

This scheme is depicted in Figure 7(b).

Benefits. The main advantage of the DPMH scheme is that the generation of samples can be
parallelized (i.e., fixing the computational cost, DPMH allows the use of M processors in parallel)
e0m , and
and the communication to the central node requires the transfer of only M particles, x
0
M weights, Zbm , instead of N M particles and N M weights. Figure 7 provides a general sketch
of DPMH. Its marginal version is described in Appendix D. Another benefit of DPMH is that
different types of particle filters can be jointly employed, for instance, different proposal pdfs can
be used.
Special cases and extensions. The classical PMH method is as a special case of the proposed
algorithm of Table 7 when M = 1. If the partial estimators are transmitted to the central node,
as shown in Figure 7(b), DPMH coincides with PGMS when M = 1. Adaptive versions of DPMH
can be designed in order select automatically the best proposal pdf among the M densities, based
b
of the weights W m = PMZm Zb , m = 1, . . . , M . For instance, Figure 12(b) shows that DPMH is
j=1

j

able to detect the best scale parameters within the M used values.
5

This particle approximation can be interpreted as being obtained by a single particle filter splitting the particles
in M disjoint sets and then applying the partial resampling described in Section 4.1, i.e., performing resampling
steps within the sets. See also Eq. (56).

19

(a)

(b)

Figure 7: Graphical representation of Distributed Particle Metropolis-Hastings (DPMH) method, (a)
for estimating a generic moment, or (b) for estimating of a specific moment of the target.

6

Numerical Experiments

In this section, we test the novel techniques considering several experimental scenarios and three
different applications: hyperparameters estimation for Gaussian Processes (D = 2, ξ = 1), a
localization problem jointly with the tuning of the sensor network (D = 8, ξ = 1), and the
online filtering of a remote sensing variable called Leaf Area Index (LAI; D = 365, ξ = 1). We
compare the novel algorithms with different benchmark methods such adaptive MH algorithm,
MTM and PMH techniques, parallel MH chains with random walk proposal pdfs, IS schemes, and
the Adaptive Multiple Importance Sampling (AMIS) method.

6.1

Hyperparameter tuning for Gaussian Process (GP) regression
models

We test the proposed GMS approach for the estimation of hyperparameters of a Gaussian process
(GP) regression model [4], [37]. Let us assume observed data pairs {yj , zj }Pj=1 , with yj ∈ R and
zj ∈ RL . We also denote the corresponding P ×1 output vector as y = [y1 , . . . , yP ]> and the L×P
input matrix as Z = [z1 , . . . , zP ]. We address the regression problem of inferring the unknown
function f which links the variable y and z. Thus, the assumed model is y = f (z) + e, where
e ∼ N (e; 0, σ 2 ), and that f (z) is a realization of a GP [37]. Hence f (z) ∼ GP(µ(z), κ(z, r)) where
µ(z) = 0, z, r ∈ RL , and we consider the kernel function
!
L
X
(z` − r` )2
κ(z, r) = exp −
.
(39)
2δ 2
`=1

20

Given these assumptions, the vector f = [f (z1 ), . . . , f (zP )]> is distributed as p(f |Z, δ, κ) =
N (f ; 0, K), where 0 is a P × 1 null vector, and Kij := κ(zi , zj ), for all i, j = 1, . . . , P , is a P × P
matrix. Therefore, the vector containing all the hyperparameters of the model is x = [δ, σ], i.e.,
all the parameters of the kernel function in Eq. (39) and standard deviation σ of the observation
noise. In this experiment, we focus on the marginal posterior density of the hyperparameters,
π̄(x|y, Z, κ) ∝ π(x|y, Z, κ) = p(y|x, Z, κ)p(x), which can be evaluated analytically, but we cannot
compute integrals involving it [37]. Considering a uniform prior within [0, 20]2 , p(x) and since
p(y|x, Z, κ) = N (y; 0, K + σ 2 I), we have


1
1
log [π(x|y, Z, κ)] = − y> (K + σ 2 I)−1 y − log det K + σ 2 I ,
2
2

where clearly K depends on δ [37]. The moments of this marginal posterior cannot be computed
analytically. Then, in order to compute the Minimum Mean Square Error (MMSE) estimator
bσ
b = [δ,
x
b], i.e., the expected value E[X] with X ∼ π̄(x|y, Z, κ), we approximate E[X] via Monte
Carlo quadrature. More specifically, we apply a the novel GMS technique and compare with
an MTM sampler, a MH scheme with a longer chain and a static IS method. For all these
methodologies, we consider the same number of target evaluations, denoted as E, in order to
provide a fair comparison.
We generated P = 200 pairs of data, {yj , zj }Pj=1 , according to the GP model above setting
δ ∗ = 3, σ ∗ = 10, L = 1, and drawing zj ∼ U([0, 10]). We keep fixed these data over the different
runs, and the corresponding posterior pdf is given in Figure 9(b). We computed the groundb = [δb = 3.5200, σ
truth x
b = 9.2811] using an exhaustive and costly grid approximation, in order
to compare the different techniques. For both GMS, MTM and MH schemes, we consider the
same adaptive Gaussian proposal pdf qt (x|µt , λ2 I) = N (x|µt , λ2 I), with λ = 5 and µt is adapted
considering the arithmetic mean of the outputs after a training period, t ≥ 0.2T , in the same
fashion of [25, 19] (µ0 = [1, 1]> ). First, we test both techniques fixing T = 20 and varying the
number of tries N . Then, we set N = 100 and vary the number of iterations T . Figure 8 (log-log
b averaged over 103 independent
plot) shows the Mean Square Error (MSE) in the approximation of x
runs. Observe that GMS always outperforms the corresponding MTM scheme. These results
confirm the advantage of recycling the auxiliary samples drawn at each iteration during an MTM
run. In Figure 9(a), we show the MSE obtained by GMS keeping invariant the number of target
evaluations E = N T = 103 and varying N ∈ {1, 2, 10, 20, 50, 100, 250, 103 }. As a consequence, we
have T ∈ {103 , 500, 100, 50, 20, 10, 4, 1}. Note that the case N = 1, T = 103 , corresponds to an
adaptive MH (A-MH) method with a longer chain, whereas the case N = 103 , T = 1, corresponds
to a static IS scheme (both with the same posterior evaluations E = N T = 103 ). We observe that
the GMS always provides smaller MSE than the static IS approach. Moreover, GMS outperforms
A-MH with the exception of two cases where T ∈ {1, 4}.

6.2

Localization of a target and tuning of the sensor network

We consider the problem of positioning a target in R2 using a range measurements in a wireless
sensor network [1, 21]. We also assume that the measurements are contaminated by noise with
different unknown power, one per each sensor. This situation is common in several practical
21

(a)
MSE

(b)
MSE
MTM
GMS

MTM
GMS

0

10

0

10

1

2

10

1

3

10
N

10

10

2

10

T

Figure 8: MSE (loglog-scale; averaged over 103 independent runs) obtained with the MTM and GMS
algorithms (using the same proposal pdf and the same values of N and T ) (a) as function of N with
T = 20 and (b) as function of T with N = 100.

(a)

(b)

Static IS

x2
A-MH (T = 1000)

0

10

0

10

1

10

2

N

10

3

10

x1

Figure 9: (a) MSE (loglog-scale; averaged over 103 independent runs) of GMS (circles) versus the
number of candidates N ∈ {1, 2, 10, 20, 50, 100, 250, 103 }, but keeping fixed the total number of posterior
evaluations E = N T = 1000, so that T ∈ {1000, 500, 100, 50, 20, 10, 4, 1}. The MSE values the extreme
cases N = 1, T = 1000, and N = 1000, T = 1, are depicted with dashed lines. In first case, GMS
coincides with an adaptive MH scheme (due the adaptation of the proposal, in this example) with a
longer chain. The second one represents a static IS scheme (clearly, using the sample proposal than
GMS). We can observe the benefit of the dynamic combination of IS estimators obtained by GMS. (b)
Posterior density π(x|y, Z, κ).

22

scenarios. Actually, even if the sensors have the same construction features, the noise perturbation
of each the sensor can vary with the time and depends on the location of the sensor. This occurs
owing to different causes: manufacturing defects, obstacles in the reception, different physical
environmental conditions (such as humidity and temperature) etc. Moreover, in general, these
conditions change along time, hence it is necessary that the central node of the network is able
to re-estimate the noise powers jointly with position of the target (and other parameters of the
models if required) whenever a new block of observations is processed. More specifically, let us
denote the target position with the random vector Z = [Z1 , Z2 ]> . The position of the target is then
a specific realization Z = z. The range measurements are obtained from NS = 6 sensors located
at h1 = [3, −8]> , h2 = [8, 10]> , h3 = [−4, −6]> , h4 = [−8, 1]> , h5 = [10, 0]> and h6 = [0, 10]> as
shown in Figure 10(a). The observation models are given by
Yj = 20 log (||z − hj ||) + Bj ,

j = 1, . . . , NS ,

(40)

where Bj are independent Gaussian random variables with pdfs, N (bj ; 0, λ2j ), j = 1, . . . , NS . We
denote λ = [λ1 , . . . , λNS ] the vector of standard deviations. Given the position of the target
z∗ = [z1∗ = 2.5, z2∗ = 2.5]> and setting λ∗ = [λ∗1 = 1, λ∗2 = 2, λ∗3 = 1, λ∗4 = 0.5, λ∗5 = 3, λ∗6 = 0.2]
(since NS = 6), we generate NO = 20 observations from each sensor according to the model in
Eq. (40). Then, we finally obtain a measurement matrix Y = [yk,1 , . . . , yk,NS ] ∈ RdY , where
dY = NO NS = 120, k = 1, . . . , NO . We consider uniform prior U(Rz ) over the position [z1 , z2 ]>
with Rz = [−30×30]2 , and a uniform prior over λj , so that λ has prior U(Rλ ) with Rλ = [0, 20]NS .
Thus, the posterior pdf is
π̄(x|Y) = π̄(z, λ|Y) = `(y|z1 , z2 , λ1 , . . . , λNS )

2
Y
i=1

p(zi )

NS
Y

p(λj ),

j=1




NS
NO Y
Y
1
1
q
exp − 2 (yk,j + 10 log (||z − hj ||)2  Iz (Rz )Iλ (Rλ ),
= 
2λj
2πλ2j
k=1 j=1


(41)

(42)

where x = [z, λ]> is a vector of parameters of dimension D = NS + 2 = 8 that we desire to infer,
and Ic (R) is an indicator variable that is 1 if c ∈ R, otherwise is 0.
Our goal is to compute the Minimum Mean Square Error (MMSE) estimator, i.e., the expected
value of the posterior π̄(x|Y) = π̄(z, λ|Y). Since the MMSE estimator cannot be computed
analytically, we apply Monte Carlo methods for approximating it. We compare GMS, the
corresponding MTM scheme, the Adaptive Multiple Importance Sampling (AMIS) technique [10],
and N parallel MH chains with a random walk proposal pdf. For all of them we consider Gaussian
proposal densities. For GMS and MTM, we set qt (x|µn,t , σ 2 I) = N (x|µt , σ 2 I) where is adapted
considering the empirical mean of the generated samples after a training period, t ≥ 0.2T [25, 19],
µ0 ∼ U([1, 5]D ) and σ = 1. For AMIS, we have qt (x|µt , Ct ) = N (x|µt , Ct ), where µt is as
previously described (with µ0 ∼ U([1, 5]D )) and Ct is also adapted using the empirical covariance
matrix, starting C0 = 4I. We also test the use of N parallel Metropolis-Hastings (MH) chains
(we also consider the case of N = 1, i.e., a single chain), with a Gaussian random-walk proposal
pdf, qn (µn,t |µn,t−1 , σ 2 I) = N (µn,t |µn,t−1 , σ 2 I) with µn,0 ∼ U([1, 5]D ) for all n and σ = 1.
23

We fix the total number of evaluations of the posterior density as E = N T = 104 . Note that,
generally, the evaluation of the posterior is the most costly step in MC algorithms (however,
AMIS has the additional cost of re-weighting all the samples at each iteration according to
the deterministic mixture procedure [6, 10, 16]). We recall that T denotes the total number
of iterations and N the number of samples drawn from each proposal at each iteration. We
consider x∗ = [z∗ , λ∗ ]> as the ground-truth and compute the Mean Square Error (MSE) in the
estimation obtained with the different algorithms. The results are averaged over 500 independent
runs and they are provided in Tables 8, 9, and 10 and Figure 10(b). Note that GMS outperforms
AMIS for each a pair {N, T } (keeping fixed E = N T = 104 ), and GMS also provides smaller MSE
values than N parallel MH chains (the case N = 1 corresponds to a unique longer chain). Figure
10(b) shows the MSE versus N maintaining E = N T = 104 for GMS and the corresponding MTM
method. This figure again confirms the advantage of recycling the samples in a MTM scheme.
Table 8: Results GMS.

MSE
N
T
E
MSE range

1.30
10
1000

1.24
20
500

1.22
50
200

1.21 1.22 1.19 1.31 1.44
100 200 500 1000 2000
100
50
20
10
5
N T = 104
Min MSE= 1.19 ——— Max MSE= 1.44

Table 9: Results AMIS [10].

MSE
N
T
E
MSE range

1.58
10
1000

1.57
20
500

1.53
50
200

1.48 1.42 1.29 1.48 1.71
100 200 500 1000 2000
100
50
20
10
5
N T = 104
Min MSE= 1.29 ——— Max MSE= 1.71

Table 10: Results N parallel MH chains with random-walk proposal pdf.

MSE
N
T
E
MSE range

1.42
1
104

1.31
5
2000

1.44
10
1000

2.32 2.73 3.21 3.18 3.15
50
100 500 1000 2000
200 100
20
10
5
4
N T = 10
Min MSE= 1.31 ——— Max MSE=3.21

24

(a)

(b)
0.2

10

10

MTM
GMS

x

2

5

0

−5

−10
−10

0.1

10

−5

0

x1

5

10

500

1000
N

1500

2000

Figure 10: (a) Sensor network: the location of the sensors (antennas) and the target (circle) in the
numerical example. The solid line represents the different unknown variances of the sensors. (b) MSE
(log-scale) versus the number of candidates N ∈ {50, 200, 500, 1000, 2000} obtained by GMS and the
corresponding MTM algorithm, keeping fixed the total number of evaluations E = N T = 104 of the
posterior pdf, so that T ∈ {200, 50, 20, 10, 5}.

6.3

Tracking of biophysical parameters

We consider the challenging problem of estimating biophysical parameters from remote sensing
(satellite) observations. In particular, we focus on the estimation of the Leaf Area Index (LAI).
It is important to track evolution of LAI through time in every spatial position on Earth because
LAI plays an important role in vegetation processes such as photosynthesis and transpiration,
and is connected to meteorological/climate and ecological land processes [9]. Let us denote LAI
as xd ∈ R+ (where d ∈ N+ also represents a temporal index) in a specific region at a latitude
of 42◦ N [18]. Since xt > 0, we consider Gamma prior pdfs over the evolutions of LAI and
Gaussian perturbations for the “in-situ” received measurements, yt . More specifically, we assume
the state-space model (formed by propagation and measurement equations),


(

(x
−b)/b
,
b
= c1d xd d−1
exp − xbd ,
gd (xd |xd−1 ) = G xd xd−1
b
(43)

1
`d (yd |xd )
= N (yd |xd , λ2 ) = √2πλ
exp − 2λ1 2 (yd − xd )2 ,
2
for d = 2, . . . , D, with initial probability g1 (x1 ) = G(x1 |1, 1), where b, λ > 0 and cd > 0 is a
normalizing constant. Note that the expected value of the Gamma pdf above is xd−1 and the
variance is b.
First Experiment. Considering that all the parameters of the model are known, the posterior
pdf is
"D
#" D
!
#
Y
Y
π̄(x|y) ∝ `(y|x)g(x) =
`d (yd |xd )
gd (xd |xd−1 ) g1 (x1 ) ,
(44)
d=2

d=2

25

with x = x1:D ∈ RD . For generating the ground-truth (i.e., the trajectory x∗ = x∗1:D =
[x∗1 , . . . , x∗D ]), we simulate the temporal evolution of LAI in one year (i.e., 1 ≤ d ≤ D = 365)
by using a double logistic function (as suggested in the literature [18]), i.e.,


1
1
+
+1 ,
(45)
xd = a1 + a2
1 + exp(a3 (d − a4 )) 1 + exp(a5 (d − a6 ))
with a1 = 0.1, a2 = 5, a3 = −0.29, a4 = 120, a5 = 0.1 and a6 = 240 as employed in
[18]. In Figure 11, the true trajectory x1:D is depicted with dashed lines. The observations

1
exp − 2λ1 2 (yd − xd )2 .
y = y2:D are then generated (each run) according to yd ∼ `d (yd |xd ) = √2πλ
2
First of all, we test the standard PMH, the particle version of GMS (PGMS), and DPMH
(fixing λ = 0.1). For DPMH, we use M = 4 parallel filters with different scale parameters
>
b = [b1 = 0.01, bP
2 = 0.05, b3 = 0.1, b4 = 1] . Figure 11 shows the estimated trajectories
t
1
eτ (averaged over 2000 runs) obtained by DPMH with N = 5 at
bt = x
x
b1:D,t = t τ =1 x
t ∈ {2, 10, 100}, in one specific run. Figure 12(a) depicts the evolution of the MSE obtained by
DPMH as a function of T and considering different values of N ∈ {5, 7, 10, 20}. The performance
of DPMH improves as T and N grow, as expected. DPMH detects the best parameters among
the four values in b, following the weights W m (see Figure 12(b)) and DPMH takes advantage of
this ability. Indeed, we compare DPMH with N = 10, T = 200, and M = 4 using b, with M = 4
different standard PMH and PGMS algorithms with N = 40 and T = 200 (clearly, each one driven
by a unique filter, M = 1) in order to keep the total number of evaluation of the posterior fixed,
E = N M T = 8 · 103 , each one using a parameter bm , m = 1, . . . , M . The results, averaged over
2000 runs, are shown in Table 11. In terms of MSE, DMPMH always outperforms the 4 possible
standard PMH methods. PGMS using two parameters, b2 and b3 , provides better performance,
but DPMH outperforms PGMS averaging the 4 different MSEs obtained by PGMS. Moreover,
due to the parallelization, in this case DPMH can save ≈ 15% of the spent computational time.
Second Experiment. Now we consider also the parameter λ unknown, so that the complete
variable of interest [x, λ] ∈ RD+1 . Then the posterior is π̄(x, λ|y) ∝ `(y|x, λ)g(x, λ) according
to the model Eq. (43), where g(x, λ) = g(x)gλ (λ) and gλ (λ) is a uniform pdf in [0.01, 5].
Then we test the marginal versions of the PMH and DPMH with qλ (λ) = gλ (λ) (see App.
D), for estimating [x∗ , λ∗ ] where x∗ = x∗1:D is given by Eq. (45) and λ∗ = 0.7. Figure
13 shows the MSE in estimation of λ∗ (averaged over 1000 runs) obtained by DPMMH as
function of T and different number of candidates, N ∈ {5, 10, 20} (with again M = 4 and
b = [b1 = 0.01, b2 = 0.05, b3 = 0.1, b4 = 1]> ). Table 12 compares the standard PMMH and
DPMMH for estimating λ∗ (we set E = N M T = 4 · 103 and T = 100). Averaging the results of
PMMH, we can observe that DPMMH outperforms the standard PMMH in terms of smaller MSE
and smaller computational time.

26

(a) t = 2

(b) t = 10

17

Est
True

16

(c) t = 100

17

Est
True

16

17

15

15

15

14

14

14

13

13

13

12

12

12

11

11

11

10

10

10

9
0

100

200
k

9
0

300

100

200
k

300

Est
True

16

9
0

100

200
k

300

Figure 11: Output of DPMH (with N = 5, λ = 0.1 and b = [0.01, 0.05, 0.1, 1]> ) at different iterations
(a) t = 2, (b) t = 10, and (c) t = 100, in one specific run. The true values, x∗ = x∗1:D , are shown dashed
bt = x
lines whereas the estimated trajectories by DPMH, x
b1:D,t , with solid lines.

(a)

(b)

1

1

10

N=5
N=7
N=10
N=20

0

N=5
0.5

10

W2
W3

W1

W4

0
0

0.2

0.6

0.8

1

1

−1

10

W2

0.5
−2

10

0.4

N=20

W3
W4

W1

0

10

1

10

0
0

2

T

10

0.2

0.4

bσm0.6

0.8

1

m

Figure 12: (a) MSE in estimation of the trajectory (averaged over 2000 runs) obtained by DPMH as
function T and different values of N ∈ {5, 7, 10, 20}. As expected, we can see that the performance of
bm
Z
DPMH improves as T and N grow. (b) Averaged values of the normalized weights W m = PM
b (with
j=1

Zj

N = 5 and N = 10) associated to each filter. DPMH is able to detect the best variances (b2 and b3 ) of
the proposal pdfs among the values b1 = 0.01, b2 = 0.05, b3 = 0.1 and b4 = 1 (as confirmed by Table 11).

27

N=5
N=10
N=20

0

10

−1

10

1

10

2

T

10

Figure 13: MSE in estimation of λ∗ = 0.7 (averaged over 1000 runs) obtained by DPMMH as function
T and different values of N ∈ {5, 10, 20}.
Table 11: Comparison among PMH, PGMS and DPMH with E = N M T = 8 · 103 and T = 200
(λ = 0.1), estimating the trajectory x∗ = x∗1:D .

Proposal Var
b1 = 0.01
b2 = 0.05
b3 = 0.1
b4 = 1
Average
Norm. Time

7

Standard PMH
N = 40
(M = 1)
MSE
0.0422
0.0130
0.0133
0.0178
0.0216
1

PGMS
N = 40
(M = 1)
MSE
0.0380
0.0100
0.0102
0.0140
0.0181
1

DPMH
N = 10
M =4
MSE
0.0108
0.0108
0.83

Conclusions

In this work, we have described the Group Importance Sampling (GIS) theory and its application
in other Monte Carlo schemes. We have considered the use of GIS in SIR (a.k.a., particle filtering),
showing that GIS is strictly required if the resampling procedure is applied only in a subset of the
current population of particles. Moreover we have highlighted that, in the standard SIR method,
if GIS is applied there exists two equivalent estimators of the marginal likelihood (one of them is
an estimator of the marginal likelihood only if the GIS weighting is used), exactly as in Sequential
Importance Sampling (SIS). We have also shown that the Independent Multiple Try Metropolis
(I-MTM) schemes and the Particle Metropolis-Hastings (PMH) algorithm can be interpreted as a
classical Metropolis-Hastings (MH) method taking into account the GIS approach.
Furthermore, two novel methodologies based on GIS have been introduced. One of them (GMS)
yields a Markov chain of weighted samples and can be also considered an iterative importance

28

Table 12: Comparison among PMMH and DMPMH with E = N M T = 4 · 103 and T = 100, for
estimating λ∗ = 0.7.

Proposal Var
b1 = 0.01
b2 = 0.05
b3 = 0.1
b4 = 1
Average
Norm. Time

PMMH
N = 40
(M = 1)
MSE
0.0929
0.0186
0.0401
0.0223
0.0435
1

DPMMH
N = 10
M =4
MSE
0.0234
0.0234
0.85

sampler. The second one (DPMH) is a distributed version of the PMH where different parallel
particle filters can be jointly employed. These filter cooperate for driving the PMH scheme. Both
techniques have been applied successfully in three different numerical experiments (tuning of the
hyperparameters for GPs, a localization problem in a wireless sensor network and the tracking
of the Leaf Area Index), comparing them with several benchmark methods. Marginal versions
of GMS and DPMH have been also discussed and tested in the numerical applications. Three
Matlab demos have been also given in order to facilitate the comprehension of the reader. As a
future line, we plan to design an adaptive DPMH scheme in order to select online the best particle
filters among the M run in parallel, and parsimoniously distribute the computational effort.

Acknowledgements
This work has been supported by the European Research Council (ERC) through the ERC
Consolidator Grant SEDAL ERC-2014-CoG 647423.

References
[1] A. M. Ali, K. Yao, T. C. Collier, E. Taylor, D. Blumstein, and L. Girod. An empirical study of
collaborative acoustic source localization. Proc. Information Processing in Sensor Networks
(IPSN07), Boston, April 2007.
[2] C. Andrieu, A. Doucet, and R. Holenstein. Particle Markov chain Monte Carlo methods. J.
R. Statist. Soc. B, 72(3):269–342, 2010.
[3] M. Bédard, R. Douc, and E. Mouline. Scaling analysis of multiple-try MCMC methods.
Stochastic Processes and their Applications, 122:758–786, 2012.
[4] C. Bishop. Pattern Recognition and Machine Learning. Springer, 2006.
29

[5] M. Bolić, P. M. Djurić, and S. Hong. Resampling algorithms and architectures for distributed
particle filters. IEEE Transactions Signal Processing, 53(7):2442–2450, 2005.
[6] M. F. Bugallo, L. Martino, and J. Corander. Adaptive importance sampling in signal
processing. Digital Signal Processing, 47:36–49, 2015.
[7] R. Casarin, R. V. Craiu, and F. Leisen. Interacting multiple try algorithms with different
proposal distributions. Statistics and Computing, 23(2):185–200, 2013.
[8] G. Casella and C. P. Robert. Rao-Blackwellisation of sampling schemes. Biometrika, 83(1):81–
94, 1996.
[9] J. M. Chen and T. A Black. Defining leaf area index for non-flat leaves. Plant, Cell and
Environment, 15:421–429, 1992.
[10] J. M. Cornuet, J. M. Marin, A. Mira, and C. P. Robert. Adaptive multiple importance
sampling. Scandinavian Journal of Statistics, 39(4):798–812, December 2012.
[11] R. V. Craiu and C. Lemieux. Acceleration of the Multiple Try Metropolis algorithm using
antithetic and stratified sampling. Statistics and Computing, 17(2):109–120, 2007.
[12] P. M. Djurić, J. H. Kotecha, J. Zhang, Y. Huang, T. Ghirmai, M. F. Bugallo, and J. Mı́guez.
Particle filtering. IEEE Signal Processing Magazine, 20(5):19–38, September 2003.
[13] A. Doucet, N. de Freitas, and N. Gordon, editors. Sequential Monte Carlo Methods in Practice.
Springer, New York (USA), 2001.
[14] A. Doucet and A. M. Johansen. A tutorial on particle filtering and smoothing: fifteen years
later. technical report, 2008.
[15] C. C. Drovandi, J. McGree, and A. N. Pettitt. A sequential Monte Carlo algorithm to
incorporate model uncertainty in Bayesian sequential design. Journal of Computational and
Graphical Statistics, 23(1):3–24, 2014.
[16] V. Elvira, L. Martino, D. Luengo, and M. F. Bugallo. Generalized multiple importance
sampling. arXiv:1511.03095, 2015.
[17] W. J. Fitzgerald. Markov chain Monte Carlo methods with applications to signal processing.
Signal Processing, 81(1):3–18, January 2001.
[18] J. L. Gomez-Dans, P. E. Lewis, and M. Disney. Efficient emulation of radiative transfer
codes using Gaussian Processes and application to land surface parameter inferences. Remote
Sensing, 8(2), 2016.
[19] H. Haario, E. Saksman, and J. Tamminen. An adaptive Metropolis algorithm. Bernoulli,
7(2):223–242, April 2001.

30

[20] J. H Huggins and D. M Roy. Convergence of sequential Monte Carlo based sampling methods.
arXiv:1503.00966, 2015.
[21] A. T. Ihler, J. W. Fisher, R. L. Moses, and A. S. Willsky. Nonparametric belief propagation for
self-localization of sensor networks. IEEE Transactions on Selected Areas in Communications,
23(4):809–819, April 2005.
[22] R. Lamberti, Y. Petetin, F. Septier, and F. Desbouvries. An improved sir-based sequential
monte carlo algorithm. In IEEE Statistical Signal Processing Workshop (SSP), pages 1–5,
2016.
[23] F. Liang, C. Liu, and R. Caroll. Advanced Markov Chain Monte Carlo Methods: Learning
from Past Samples. Wiley Series in Computational Statistics, England, 2010.
[24] J. S. Liu. Monte Carlo Strategies in Scientific Computing. Springer, 2004.
[25] D. Luengo and L. Martino. Fully adaptive Gaussian mixture Metropolis-Hastings algorithm.
Proceedings of IEEE International Conference on Acoustics, Speech, and Signal Processing
(ICASSP), 2013.
[26] L. Martino, V. Elvira, and F. Louzada. Weighting a resampled particle in Sequential Monte
Carlo. IEEE Statistical Signal Processing Workshop, (SSP), 122:1–5, 2016.
[27] L. Martino, V. Elvira, and M. F. Louzada. Effective Sample Size for importance sampling
based on the discrepancy measures. Signal Processing, 131:386–401, 2017.
[28] L. Martino, V. Elvira, D. Luengo, and J. Corander. Layered adaptive importance sampling.
Statistics and Computing, 27(3):599–623, 2017.
[29] L. Martino, F. Leisen, and J. Corander. On multiple try schemes and the Particle MetropolisHastings algorithm. viXra:1409.0051, 2014.
[30] L. Martino and F. Louzada. Issues in the Multiple Try Metropolis mixing. Computational
Statistics, 32(1):239–252, 2017.
[31] L. Martino, V. P. Del Olmo, and J. Read. A multi-point Metropolis scheme with generic
weight functions. Statistics & Probability Letters, 82(7):1445–1453, 2012.
[32] L. Martino and J. Read. On the flexibility of the design of multiple try Metropolis schemes.
Computational Statistics, 28(6):2797–2823, 2013.
[33] L. Martino, J. Read, V. Elvira, and F. Louzada. Cooperative parallel particle filters for on-line
model selection and applications to urban mobility. Digital Signal Processing, 60:172–185,
2017.
[34] J. Mı́guez and M. A. Vázquez. A proof of uniform convergence over time for a distributed
particle filter. Signal Processing, 122:152–163, 2016.
31

[35] C. A. Naesseth, F. Lindsten, and T. B. Schon. Nested Sequential Monte Carlo methods.
Proceedings of theInternational Conference on Machine Learning, 37:1–10, 2015.
[36] C. A. Naesseth, F. Lindsten, and T. B. Schon. High-dimensional filtering using nested
sequential monte carlo. arXiv:1612.09162, pages 1–48, 2016.
[37] C. E. Rasmussen and C. K. I. Williams. Gaussian processes for machine learning. MIT Press,
2006.
[38] J. Read, K. Achutegui, and J. Mı́guez. A distributed particle filter for nonlinear tracking in
wireless sensor networks. Signal Processing, 98:121 – 134, 2014.
[39] C. P. Robert and G. Casella. Monte Carlo Statistical Methods. Springer, 2004.
[40] D. B. Rubin. Using the SIR algorithm to simulate posterior distributions. in Bayesian
Statistics 3, ads Bernardo, Degroot, Lindley, and Smith. Oxford University Press, Oxford,
1988., 1988.
[41] R. Bassi Stern. A statistical contribution to historical linguistics. Phd Thesis, 2015.
[42] I. Urteaga, M. F. Bugallo, and P. M. Djurić. Sequential Monte Carlo methods under model
uncertainty. In 2016 IEEE Statistical Signal Processing Workshop (SSP), pages 1–5, 2016.
[43] C. Verg, C. Dubarry, P. Del Moral, and E. Moulines. On parallel implementation of sequential
Monte Carlo methods: the island particle model. Statistics and Computing, 25(2):243–260,
2015.
[44] C. Verg, P. Del Moral, E. Moulines, and J. Olsson. Convergence properties of weighted
particle islands with application to the double bootstrap algorithm. arXiv:1410.4231, pages
1–39, 2014.
[45] N. Whiteley, A. Lee, and K. Heine. On the role of interaction in sequential Monte Carlo
algorithms. Bernoulli, 22(1):494–529, 2016.

A

Proper weighting of a resampled particle

Let us consider the particle approximation of π̄ obtained by the IS approach drawing N particles
xn ∼ q(x),
N
N
X
1 X
w(xn )δ(x − xn ).
(46)
π
b(x|x1:N ) =
w̄(xn )δ(x − xn ) =
b
N
Z
n=1
n=1
QN
e0 ∼ π
Given the cloud of particle x1:N h∼ n=1 q(xin ), we have that x
b(x|x1:N ). Let us denote the
Q
N
e
e0 ,
joint pdf Q(x, x1:N ) = π
b(x|x1:N )
e(x) of a resampled particle x
i=1 q(xi ) . The marginal pdf q
32

e0 ∼ qe(x)), is
integrating out x1:N (i.e., x
Z
e x1:N )dx1:N
qe(x) =
Q(x,
XN
"N
#
Z
Y
=
π
b(x|x1:N )
q(xi ) dx1:N ,
XN

N
X

=

j=1

= π(x)

N
X
j=1

= π(x)

Z









X N −1

i=1

Z

X N −1

Z

X N −1



w(x) 
q(x)
N Zb

N
Y
i=1
i6=j







q(xi ) dx¬j  ,




N
1 Y


 q(xi ) dx¬j  ,
N Zb
i=1
i6=j



N
1 Y

 q(xi ) dx¬j ,
b
Z i=1

j ∈ {1, . . . , N }.

(47)

i6=j

e0 ∼ qe(x), is
Therefore, the standard IS weight of a resampled particle, x
w(e
x0 ) =

π(e
x0 )
.
qe(e
x0 )

(48)

However, generally qe(x) cannot be evaluated, hence the standard IS weight cannot be computed
[22, 26, 35], [28, App. C1]. An alternative is to use the Liu’s definition of proper weighting in Eq.
(8) and look for a weight function ρ(e
x) = ρ(e
x|x1:N ) such that
EQ(x,x
[ρ(x|x1:N )h(x)] = cEπ̄ [h(x)],
e
1:N )

(49)

i
hQ
N
e x1:N ) = π
where Q(x,
b(x|x1:N )
q(x
)
i . Below, we show that a suitable choice is
i=1

since it holds in Eq. (49).

b 1:N ) =
ρ(e
x|x1:N ) = Z(x

Proof. Note that
EQ(x,x
[ρ(x|x1:N )h(x)] =
e
1:N )
=

Z Z
X

XN

X

XN

Z Z

N
1 X
w(xi ),
N i=1

e x1:N )dxdx1:N ,
ρ(x|x1:N )h(x)Q(x,
"N
#
Y
h(x)ρ(x|x1:N )b
π (x|x1:N )
q(xi ) dxdx1:N .
i=1

33

(50)

(51)

Recalling that π
b(x|x1:N ) =

w(xn ) =

π(xn )
,
q(xn )

1
b
NZ

PN

j=1

b 1:N ) =
w(xj )δ(x − xj ), where Zb = Z(x

we can rearrange the expectation above as

EQ(x,x
[ρ(x|x1:N )h(x)] =
e
1:N )

=

Z

X

Z

X



Z
N
X 
h(x) 


ρ(x|x1:N )

X N −1

j=1



w(x) 
q(x)
N Zb

N
Y
i=1
i6=j

1
N

PN

n=1



w(xn ) and




q(xi ) dx¬j  dx,




N
1 Y




h(x)π(x) 
ρ(x|x1:N )
 q(xi ) dx¬j  dx,

b
N
−1
N Z i=1
X
j=1


N
X



Z

i6=j

where x¬j = [x1 , . . . , xj−1 , xj+1 , . . . , xN ]. If we choose ρ(x|x1:N ) = Zb and replace it in the
expression above, we obtain





Z
Z
N
N
Y
X 


b 1 
q(xi ) dx¬j  dx,
h(x)π(x)
Z
[ρ(x|x
)h(x)]
=
EQ(x,x



e
1:N
1:N )
N Zb i=1
X
X N −1
j=1
i6=j
Z
1
h(x)π(x)N dx,
=
N
ZX
=
h(x)π(x)dx
X

= cEπ̄ [h(x)],

(52)

where c = Z, that is the normalizing constant of π(x). Note that Eq. (52) coincides with (49).
2

B

Particle approximation by GIS

P
π(xm,n )
Let us consider S samples xm,n ∼ qm (x), where S = M
m=1 Nm , and weight them wm,n = qm (xm,n )
with m = 1, . . . , M and n = 1, . . . , Nm . Moreover, let us define two types of normalized weights,
one within the m-th group
wm,n
wm,n
w̄m,n = PN
=
,
Nm Zbm
k=1 wm,k

(53)

and the other one considering all the S samples,
wm,n
r̄m,n = PM PNj
j=1

k=1

wj,k

34

wm,n
= PM
.
b
j=1 Nj Zj

(54)

The complete particle approximation of the target distribution is
π
b(x|x1:M,1:N ) = PM

j=1

1
PNj

M X
Nm
X

=

m=1 n=1

k=1

wj,k

M X
Nm
X

m=1 n=1

wm,n δ(x − xm,n ),

r̄m,n δ(x − xm,n ).

(55)

Note that it can be also rewritten as
π
b(x|x1:M,1:N ) = PM

1

j=1

= PM
=

Nj Zbj

1

b
j=1 Nj Zj

M
X

m=1

M
X

m=1
M
X

m=1

Nm Zbm

N
X
n=1

w̄m,n δ(x − xm,n ),

Nm Zbm π
b(x|xm,1:N ),

W mπ
b(x|xm,1:N ),

where π
b(x|xm,1:N ) are the m-th particle approximation and W m =

b
N Z
PMm m b
N
j=1 j Zj

(56)

is the normalized

em ∼ π
weight of the m-th group. If we resample M times x
b(x|xm,1:N ) exactly one sample per
group, we obtain the particle approximation of Eq. (12), i.e.,
π
b(x|e
x1:M ) =

M
X

m=1

em ).
W m δ(x − x

(57)

Since π
b(x|x1:M,1:N ) is a particle approximation of the target distribution π̄ (converging to the
distribution for N → ∞), then π
b(x|e
x1:M ) is also a particle approximation of π̄ (converging for
N → ∞ and M → ∞).Therefore, any estimator of the moments of π̄ obtained using the summary
weighted particles as in Eq. (13) is consistent.

C

Estimators of the marginal likelihood in SIS and SIR

The classical IS estimator of the normalizing constant Zd =
is

R

Rd×η

πd (x1:d )dx1:d at the d-th iteration

N
N
1 X (n)
1 X (n) (n)
w =
w β ,
N n=1 d
N n=1 d−1 d
" d
#
N
1 X Y (n)
β
.
=
N n=1 j=1 j

Zbd =

35

(58)
(59)

An alternative formulation, denoted as Z d , is often used
#
" N
d
Y
X (n) (n)
Zd =
w̄j−1 βj
j=1

=

d
Y
j=1

(n)

where we have employed w̄j−1 =

n=1

" PN

(n)

n=1 wj

PN

(n)
n=1 wj−1

(n)

wj−1
PN
(i)
i=1 wj−1

#

(n)

and wj

= Zb1

d
Y
j=2

(n)

(60)
"

Zbj
Zbj−1

(n)

= wj−1 βj

#

= Zbd .

(61)

[13, 14]. Therefore, given Eq. (61),

in SIS these two estimators Zbd in Eq. (58) and Z d in Eq. (60) are equivalent approximations
of the d-th marginal likelthood Zd [33]. Furthermore, note that Z d can be written in a recursive
form as
" N
#
X (n) (n)
Z d = Z d−1
w̄d−1 βd .
(62)
n=1

C.1

Estimators of the marginal likelihood in particle filtering

Sequential Importance Resampling (SIR) (a.k.a., particle filtering) combines the SIS approach
with the application of the resampling procedure corresponding to step 2(c)ii of Table 2. If the
GIS weighting is not applied, in SIR only
" N
#
d
Y
X (n) (n)
Zd =
w̄j−1 βj
.
j=1

n=1

P
(n)
is a consistent estimator of Zd . In this case, Zbd = N1 N
n=1 wd is not a possible alternative without
using GIS. However, considering the proper GIS weighting of the resampled particles (the step
2(c)iii of Table 2), then Zbd is also a consistent estimator of Zd and it is equivalent to Z d . Below,
we analyze three cases considering a resampling applied to the entire set of particles:
• No Resampling (η = 0): this scenario corresponds to SIS where Zbd , Z d are equivalent as
shown in Eq. (61).

(n)
• Resampling at each iteration (η = 1): using the GIS weighting, wd−1 = Zbd−1 for all n
and for all d, and replacing in Eq. (58) we have
"
#
N
X
1
(n)
Zbd = Zbd−1
β
,
(63)
N n=1 d
" N
#
d
1 Y X (n)
=
β
.
(64)
N j=1 n=1 j

36

(n)

Since the resampling is applied to the entire set of particles, we have w̄d−1 =
Replacing it in the expression of Z d in (62), we obtain
#
" N
d
1 Y X (n)
,
Zd =
β
N j=1 n=1 j

1
N

for all n.

(65)

that coincides with Zbd in Eq. (64).

• Adaptive resampling (0 < η < 1): for the sake of simplicity, let us start considering a
unique resampling step applied at the k-th iteation with k < d. We check if both estimators
are equal at d-th iteration of the recursion. Due to Eq. (61), we have Z k ≡ Zbk ,6 since before
(n)
the k-th iteration no resampling has been applied. With the proper weighting wk = Zbk for
all n, at the next iteration we have
"
#
N
N
X
X
1
1
(n) (n)
(n)
,
w β
= Zbk
β
Zbk+1 =
N n=1 k k+1
N n=1 k+1
and using Eq. (62), we obtain

Z k+1 = Z k

"

#
"
#
N
N
X
X
1
1 (n)
(n)
,
βk+1 = Zbk
β
N
N n=1 k+1
n=1

so that the estimators are equivalent also at the (k + 1)-th iteration, Z k+1 ≡ Zbk+1 . Since
we are assuming no resampling steps after the k-th iteration and until the d-th iteration, we
have that Z i ≡ Zbi for i = k + 2, . . . , d due to we are in a SIS scenario for i > k (see Eq.
(61)). This reasoning can be easily extended for different number of resampling steps.

Figure 14 summarizes the expressions of the estimators in the extreme cases of η = 0 and
η = 1. Note that the operations of sum and product are inverted. See DEMO-1 at https:
//github.com/lukafree/GIS.git.

D

Particle Marginal
algorithms

Metropolis-Hastings

(PMMH)

Let us consider x = x1:D = [x1 , x2 , . . . , xD ] ∈ X ⊆ RD×ξ where xd ∈ Rξ for all d = 1, . . . , D and
an additional model parameter θ ∈ Rdθ to be inferred as well. Assuming a prior pdf gθ (θ) over θ,
and a factorized complete posterior pdf π̄(x, θ)
QD

π̄c (x, θ) ∝ πc (x, θ) = gθ (θ)π(x|θ),

where π(x|θ) = γ1 (x1 |θ) d=2 γd (xd |x1:d−1 , θ). Moreover, let us the denote as π
b(x|v1:N , θ) =
PN
1
0
w(v
|θ
)δ(x−v
)
a
particle
approximation
of
π(x|θ)
obtained
by
one
run
of a particle
n
n
b 0)
n=1
N Z(θ
6

We consider to compute the estimators before the resampling.

37

bd in two extreme scenarios:
Figure 14: Expressions of the marginal likelihood estimators Z d and Z
without resampling and applying resampling at each iterations. Note that in the formulations above the
operations of sum and product are inverted.

filter approach. The Marginal PMH (PMMH) technique is then summarized in Table 13. PMMH
is often used for both smoothing and parameter estimation in state-space models. Note that if
qθ (θ|θt−1 ) = gθ (θ) then the acceptance function becomes
"
#
b 0)
Z(θ
α = min 1,
.
(66)
b t−1 )
Z(θ
Table 13: Particle Marginal MH (PMMH) algorithm

1. Choose an initial state x0 and Zbm,0 for m = 1, . . . , M .
2. For t = 1, . . . , T :

b 0 ).
(a) Choose an initial state x0 , θ0 , Z(θ

(b) For t = 1, . . . , T :

PN
1
0
i. Draw θ 0 ∼ qθ (θ|θt−1 ) and vj ∼ π
b(x|v1:N , θ 0 ) = N Z(θ
b 0)
n=1 w(vn |θ )δ(x − vn )
(where π
b is obtained with one run of a particle filter).
ii. Set θt = θ 0 , xt = vj , with probability
"
#
b 0 )gθ (θ 0 )qθ (θt−1 |θ 0 )
Z(θ
α = min 1,
.
(67)
b t−1 )gθ (θt−1 )qθ (θ|θt−1 )
Z(θ
Otherwise, set θt = θ 0 and xt = xt−1 .

3. Return {xt }Tt=1 and {θt }Tt=1 .

38

Distributed Particle Marginal Metropolis-Hastings (DPMMH). We can easily design
a marginal version of DPMH in Section 5.2, drawing θ 0 ∼ qθ (θ|θt−1 ) and run M particle filters
addressing the target pdf π̄(x|θ 0 ). The algorithm follows the steps in Table 13 with the difference
that M parallel particle filters are used, and in this case the acceptance probability is
hP
i


M
bm (θ) gθ (θ 0 )qθ (θt−1 |θ 0 )
Z
m=1
.
i
(68)
α = min 1, hP
M
b
m=1 Zm (θt−1 ) gθ (θt−1 )qθ (θ|θt−1 )
Table 14: Main acronyms in the work.

GP
MSE
IS
SIS
SIR
PF
SMC
AMIS
MCMC
MH
IMH
MTM
I-MTM
I-MTM2
PMH
PMMH
GIS
GSM
PGSM
DPMH
DPMMH

Gaussian Process
Mean Square Error
Importance Sampling
Sequential Importance Sampling
Sequential Importance Resampling
Particle Filter
Sequential Monte Carlo
Adaptive Multiple Importance Sampling
Markov Chain Monte Carlo
Metropolis-Hastings
Independent Metropolis-Hastings
Multiple Try Metropolis
Independent Multiple Try Metropolis
Independent Multiple Try Metropolis (version 2)
Particle Metropolis-Hastings
Particle Marginal Metropolis-Hastings
Group Importance Sampling
Group Metropolis Sampling
Particle Group Metropolis Sampling
Distributed Particle Metropolis-Hastings
Distributed Particle Marginal Metropolis-Hastings

39

