A Review of Multiple Try MCMC algorithms for Signal Processing
Luca Martino
stat.CO stat.ML






Hermite spectral collocation methods for fractional PDEs in unbounded domains
    [
    December 30, 2023
=============================================================================














Many applications in signal processing require the estimation of some parameters of interest given a set of observed data. More specifically, Bayesian inference needs the computation of a-posteriori estimators which are often expressed as complicated multi-dimensional integrals. Unfortunately, analytical expressions for these estimators cannot be found in most real-world applications, and Monte Carlo methods are the only feasible approach. A very powerful class of Monte Carlo techniques is formed by the Markov Chain Monte Carlo (MCMC) algorithms. They generate a Markov chain such that its stationary distribution coincides with the target posterior density. In this work, we perform a thorough review of MCMC methods using multiple candidates in order to select the next state of the chain, at each iteration. With respect to the classical Metropolis-Hastings method, the use of multiple try techniques foster the exploration of the sample space. We present different Multiple Try Metropolis schemes, Ensemble MCMC methods, Particle Metropolis-Hastings algorithms and the Delayed Rejection Metropolis technique.  We highlight limitations, benefits, connections and differences among the different methods, and compare them by numerical simulations.


Keywords: Markov Chain Monte Carlo, Multiple Try Metropolis, Particle Metropolis-Hastings, Particle Filtering, Monte Carlo methods, Bayesian inference

















§ INTRODUCTION









Bayesian methods have become very popular in signal processing over the last years <cit.>. 
They require the application of sophisticated Monte Carlo techniques, such as Markov chain Monte Carlo (MCMC) and particle filters, for the efficient computation of a-posteriori estimators <cit.>. More  specifically, the MCMC algorithms generate a Markov chain such that its stationary distribution coincides with the posterior probability density function (pdf) <cit.>. Typically, the only requirement is to be able to evaluate the target function, where the knowledge of the normalizing constant is usually not needed. 

The most popular MCMC method is undoubtedly the Metropolis-Hastings (MH) algorithm  <cit.>. The MH technique is a very simple method, easy to be applied: this is the reason of its success. In MH, at each iteration, one new candidate is generated  from a proposal pdf and then is properly compared with the previous state of the chain, in order to decide the next state. 
However, the performance of MH are often not satisfactory. For instance, when the posterior is multimodal, or when the dimension of the space increases, the correlation among the generated samples is usually high and, as a consequence, the variance of the resulting estimators grows. To speed up the convergence and reduce the “burn-in” period of the MH chain, several extensions have been proposed in literature. 


 In this work, we provide an exhaustive review of more sophisticated MCMC methods that, at each iteration, consider different candidates as possible new state of the chain.  More specifically, at each iteration different samples are compared by certain weights and then one of them is selected as possible future state.
 The main advantage of these algorithms is that they foster the exploration of a larger portion of the sample space, decreasing the correlation among the states of the generated chain. In this work, we describe different algorithms of this family, independently introduced in literature. 
The main contribution is to present them under the same the framework and notation, remarking differences, relationships, limitations and strengths. 
 
  All the discussed techniques yield an ergodic chain converging to the posterior density of interest (in the following, referred also as target pdf). 

 
The first scheme of this MCMC class, called Orientational Bias Monte Carlo (OBMC) <cit.>, was proposed in the context of molecular simulation. Later on a more general algorithm, called Multiple Try Metropolis (MTM), was introduced <cit.>.[MTM includes OBMC as a special case (see Section <ref>).] 
The MTM algorithm has been extensively studied and generalized in different ways <cit.>. Other techniques, alternative to the MTM schemes, are the so-called the Ensemble MCMC (EnMCMC) methods <cit.>. They follow a similar approach to MTM but employ a different acceptance function for selecting the next state of the chain. With respect to (w.r.t.) a generic MTM scheme, EnMCMC does not require any generation of auxiliary samples (as in a MTM scheme employing a generic proposal pdf) and hence, in this sense, EnMCMC are less costly. 

In all the previous techniques, the candidates are drawn in a batch way and compared jointly. In the Delayed Rejection Metropolis  (DRM) algorithm <cit.>, in case of rejection of the novel possible state, the authors suggest to perform an additional acceptance test considering a new candidate. If this candidate is again rejected, the procedure can be iterated until reaching a desired number of attempts. The main benefit of DRM is that the proposal pdf can be improved at each intermediate stage. However, the acceptance function progressively becomes more complex so that the implementation of DRM for a great number of attempts is not straightforward (compared to the implementation of a MTM  scheme with a generic number of tries).
  
 In the last years, other Monte Carlo methods which combine particle filtering and MCMC have become very popular in the signal processing community. For instance, this is the case of the Particle Metropolis Hastings (PMH) and the Particle Marginal Metropolis Hastings (PMMH) algorithms, which have been widely used in signal processing in order to make inference and smoothing about dynamical and static parameters in state space models <cit.>. PMH can be interpreted as a MTM scheme where the different candidates are generated and weighted by the use of a particle filter <cit.>. In this work, we present PMH and PMMH and discuss their connections and differences with the classical MTM approach.  Furthermore, we describe a suitable procedure for recycling some candidates in the final Monte Carlo estimators, called Group Metropolis Sampling (GMS) <cit.>. The GMS scheme can be also seen as a way of generating a chain of sets of weighted samples. Finally, note that other similar and related techniques can be found within the so-called data augmentation approach <cit.>.
 
The remaining of the paper is organized as follows. Section <ref> recalls the problem statement and some background material, introducing also the required notation. The basis of MCMC and the Metropolis-Hastings (MH) algorithm are presented in Section <ref>. Section <ref> is the core of the work, which describes the different MCMC using multiple candidates. Section <ref> provides some numerical results, applying different techniques in a hyperparameter tuning problem for a Gaussian Process regression model, and in a localization problem considering a wireless sensor network. Some conclusions are given in Section <ref>. 








§ PROBLEM STATEMENT AND PRELIMINARIES




In many signal processing applications, the goal consists in inferring a variable of interest, θ=[θ_1,…,θ_D]∈𝒟⊆ℝ^D, given a set of observations or measurements, y∈ℝ^P. In the Bayesian framework, the total knowledge about the parameters, after the data have been observed, is represented by the posterior probability density function (pdf) <cit.>, i.e.,

    π̅(θ)   =   p(θ|  y)= ℓ( y|θ) g(θ)/Z( y), 
       =   1/Zπ(θ),

where ℓ( y|θ) denotes the likelihood function (i.e., the observation model), g(θ) is the prior probability density function (pdf) and Z=Z( y) is the marginal likelihood (a.k.a., Bayesian evidence) <cit.> and π(θ)=ℓ( y|θ) g(θ).

In general Z is unknown, and it is possible to evaluate only the unnormalized target function,
π(θ) ∝π̅(θ).

The analytical study of the posterior density π̅(θ) is often unfeasible and integrals involving π̅(θ) are typically intractable <cit.>. For instance, one might be interested in the estimation of

    I=E_π[f(θ)]   =   ∫_𝒟 f(θ)π̅(θ) dθ, 
       =   1/Z∫_𝒟 f(θ)π(θ) dθ,

where f(θ) is a generic integrable function w.r.t. π̅. 

Dynamic and static parameters.

In some specific application, the variable of interest θ can be split in two disjoint parts, θ=[ x,λ], where one, x, is involved into a dynamical system (for instance, x is the hidden state in a state-space model) and the other, λ, is a static parameter (for instance, an unknown parameter of the model). The strategies for making inference about x and λ should take into account the different nature of the two parameters (e.g., see Section <ref>). The main notation and acronyms are summarized in Tables <ref>-<ref>.









 §.§ Monte Carlo integration


In many practical scenarios, the integral I cannot be computed in a closed form, and numerical approximations are typically required. Many deterministic quadrature methods are available in the literature <cit.>. However, as the dimension D of the inference problem grows (θ∈ℝ^D), the deterministic quadrature schemes become less efficient. In this case, a common approach consists of approximating the integral I in Eq. (<ref>) by using Monte Carlo (MC) quadrature <cit.>. Namely, considering T independent and identically distributed (i.i.d.) samples drawn from the posterior target pdf, i.e. θ_1,…,θ_T ∼π̅(θ),[In this work, for simplicity, we use the same notation for denoting a random variable or one realization of a random variable.] we can build the consistent estimator

    I_T=1/T∑_t=1^T f(θ_t).

I_T converges in probability to I due to the weak law of large numbers. The approximation above I_T is known as a direct (or ideal) Monte Carlo estimator if the samples θ_t are independent and identically distributed  (i.i.d.) from π̅. Unfortunately, in many practical applications, direct methods for drawing independent samples from π̅(θ) are not available. Therefore, different approaches are required, such as the Markov chain Monte Carlo (MCMC) techniques. 








§ MARKOV CHAIN MONTE CARLO (MCMC) METHODS



A MCMC algorithm generates an ergodic Markov chain with invariant (a.k.a., stationary) density given by the posterior pdf π̅(θ) <cit.>. Specifically, given a starting state θ_0, a sequence of correlated samples is generated, 
θ_0→θ_1→θ_2→ ....  →θ_T. Even if the samples are now correlated, the estimator

    I_T=1/T∑_t=1^T f(θ_t)

is consistent, regardless the starting vector θ^(0) <cit.>.[Recall we are assuming that the Markov chain is ergodic and hence the starting value is forgotten.] With respect to the direct Monte Carlo approach using i.i.d. samples, the application of an MCMC algorithm entails a loss of efficiency of the estimator I_T, since the samples are positively correlated, in general. In other words, to achieve a given variance obtained with the direct Monte Carlo estimator, it is necessary to generate more samples.

 Thus, in order to improve the performance of an MCMC technique we have to decrease the correlation among the states of the chain.[For the sake of simplicity, we use all the generated states in the final estimators, without removing any burn-in period <cit.>.]
 





 §.§ The Metropolis-Hastings (MH) algorithm


One of the most popular and widely applied MCMC algorithm is  the Metropolis-Hastings (MH) method <cit.>. Recall that  we are able to evaluate point-wise a function proportional to the target, i.e., π(θ)∝π̅(θ).  A proposal density (a pdf which is easy to draw from) is denoted as q(θ|θ_t-1)>0, with θ,θ_t-1∈ℝ^D.  In Table <ref>, we describe the standard MH algorithm in detail.





The algorithm returns the sequence of states {θ_1, θ_2, …, θ_t, …, θ_T} (or a subset of them removing the burn-in period if an estimation of its length is available). We can see that the next state θ_t can be the proposed sample θ' (with probability α) or the previous state θ_t-1 (with probability 1-α).  Under some mild regularity conditions, when t grows, the pdf of the current state θ_t converges to the target density π̅(θ)  <cit.>. The MH algorithm satisfies the so-called detailed balance condition which is sufficient to guarantee that the output chain is ergodic and has π̅ as stationary distribution <cit.>. Note that the acceptance probability α can be rewritten as 

    α(θ_t-1,θ') = min[1, π(θ')  q(θ_t-1|θ')/π(θ_t-1)q(θ'|θ_t-1)]= min[1,w(θ'|θ_t-1)/w(θ_t-1|θ')].

where we have denoted w(θ'|θ_t-1)=π(θ')/q(θ'|θ_t-1) and w(θ_t-1|θ')=π(θ_t-1)/q(θ_t-1|θ') in a similar fashion of the importance sampling weights of θ' and θ_t-1 <cit.>. If the  proposal pdf is independent from the previous state, i.e., q(θ|θ_t-1)=q(θ), the acceptance function depends on the ratio of the importance weights w(θ')=π(θ')/q(θ') and w(θ_t-1)=π(θ_t-1)/q(θ_t-1), as shown in Table <ref>. We refer to this special MH case as the Independent MH (I-MH) algorithm. It is strictly related to other techniques described in the following (e.g., see Section <ref>).












§ MCMC USING MULTIPLE CANDIDATES



In the standard MH technique described above, at each iteration one new sample θ' is generated to be tested with the previous state θ_t-1 by the acceptance probability α(θ_t-1,θ'). Other generalized MH schemes generate several candidates at each iteration to be tested as new possible state. In all these schemes, an extended acceptance probability α is properly designed in order to guarantee the ergodicity of the chain. Figure <ref> provides a graphical representation of the difference between MH and the techniques using several candidates.




Below, we describe the most important examples of this class of MCMC algorithms. In most of them, a single MH-type test is performed at each iteration whereas in other methods a sequence of tests is employed. 
Furthermore, most of these techniques use an Importance Sampling (IS) approximation of the target density <cit.> in order to improve the proposal procedure employed within a MH-type algorithm. Namely, they build an IS approximation, and then draw one sample from this approximation (resampling step). Finally, the selected sample is compared with the previous state of the chain, θ_t-1, according to a suitable generalized acceptance probability α. It can be proved that all the methodologies presented in this work yield a ergodic chain with the posterior π̅ as invariant density.   






 §.§ The Multiple Try Metropolis (MTM) algorithm













The  Multiple Try Metropolis (MTM) algorithms are examples of this class of methods, where N samples θ^(1),θ^(2),…, θ^(N) (called also “tries” or “candidates”) are drawn from the proposal pdf q(θ), at each iteration <cit.>. Then, one of them is selected according to some suitable weights. Finally, the selected candidate is accepted or rejected as new state according to a generalized probability function α. 

The MTM algorithm is given in Table <ref>. For the sake of simplicity, we have considered the use of the importance weights w(θ|θ_t-1)=π(θ)/q(θ|θ_t-1), but there is not a unique possibility, as also shown below  <cit.>. In its general form, when the proposal depends on the previous state of the chain q(θ|θ_t-1), the MTM requires the generation of N-1 auxiliary samples, v^(i), which are employed in the computation of the acceptance function α. They are needed in order to guarantee the ergodicity. Indeed,  the resulting MTM kernel satisfies the detailed balance condition, so that the chain is reversible <cit.>. Note that for N=1, we have θ^(j)=θ^(1), v^(1)=θ_t-1 and the acceptance probability of the MTM method becomes

    α(θ_t-1,θ^(1))   =   min(1,w(θ^(1)|θ_t-1)/w( v^(1)|θ^(1))), 
       =   min(1,w(θ^(1)|θ_t-1)/w(θ_t-1|θ^(1)))=min(1,π(θ^(1)) q(θ_t-1|θ^(1))/π(θ_t-1) q(θ^(1)|θ_t-1)),

that is the  acceptance probability of the classical MH technique. Several variants have been studied, for instance, with correlated tries and considering the use of different proposal pdfs <cit.>.





The MTM method in Table <ref> needs at step <ref> the generation of N-1 auxiliary samples and at step <ref> the computation of their weights (and, as a consequence, N-1 additional evaluation of the target pdf are required), that are only employed in the computation of the acceptance function α. 

 
 




  §.§.§ Generic form of the weights



The importance weights are not the unique possible choice.
It is possible to show that the MTM algorithm generates an ergodic chain with invariant density π̅, if the weight function w(θ|θ_t-1) is chosen with the form

    w(θ|θ_t-1)=π(θ)q(θ_t-1|θ) ξ(θ_t-1,θ),

where

    ξ(θ_t-1,θ)= ξ(θ,θ_t-1),     ∀θ,θ_t-1∈𝒟.

For instance, choosing ξ(θ_t-1,θ)=1/q(θ|θ_t-1) q(θ_t-1|θ), we obtain the importance weights w(θ|θ_t-1)=π(θ)/q(θ|θ_t-1) used above. If we set ξ(θ_t-1,θ)=1, we have w(θ|θ_t-1)=π(θ)q(θ_t-1|θ). Another interesting example can be employed if the proposal is symmetric, i.e., q(θ|θ_t-1)=q(θ_t-1|θ). In this case, we can choose ξ(θ_t-1,θ)=1/q(θ_t-1|θ) and then w(θ|θ_t-1)=w(θ)=π(θ), i.e., the weights only depend on the value of the target density at θ. Thus, 
MTM contains the Orientational Bias Monte Carlo (OBMC) scheme <cit.> as a special case, when a symmetric proposal pdf is employed, and then one candidate is chosen with weights proportional to the target density, i.e., w(θ|θ_t-1)=π(θ).




  §.§.§ Independent Multiple Try Metropolis (I-MTM) schemes


The MTM method described in Table <ref> requires to draw 2N-1 samples at each iteration (N candidates and N-1 auxiliary samples) and N-1 are only used in the acceptance probability function. The generation of the auxiliary points

    v^(1),…,  v^(j-1),  v^(j+1), …,  v^(N)∼ q(θ|θ^(j)),

can be avoided if the proposal pdf is independent from the previous state, i.e., q(θ|θ_t-1)=q(θ). Indeed, in this case, we should draw N-1 samples again from q(θ) at the step <ref> of Table <ref>. Since we have already drawn N samples from q(θ) at step <ref> of Table <ref>, we can set 

    v^(1)=θ^(1),…, v^(j-1)=θ^(j-1),  v^(j)=θ^(j+1)… v^(N-1)=θ^(N),

 without jeopardizing the ergodicity of the chain (recall that v^(j)=θ_t-1). Hence, we can avoid step <ref> and  the acceptance function can be rewritten as  
   
    α(θ_t-1,θ^(j))
                        	   =   min(1,w(θ^(j))+∑_n=1,n≠ j^N w(θ^(n))/w(θ_t-1)+∑_n=1,n≠ j^N w(θ^(n))).
   
The I-MTM algorithm is provided in Table <ref>.  

An I-MTM method requires only N new evaluations of the target pdf at each iteration, instead of 2N-1 new evaluations in the generic MTM scheme in Table <ref>.
Note that we can also write α(θ_t-1,θ^(j)) as 
   
    α(θ_t-1,θ^(j))   =   min(1,Z_1/Z_2),
 
 where we have denoted 
   
    Z_1=1/N∑_n=1^N w(θ^(n)),     Z_2=1/N(w(θ_t-1)+∑_n=1,n≠ j^N w(θ^(n))).
 
 

  
   
Alternative version (I-MTM2). From the IS theory, we know that Z_1=1/N∑_n=1^N w(θ^(n)) is an unbiased estimator of the normalizing constant Z of the target π (a.k.a, Bayesian evidence or marginal likelihood). It suggests to replace Z_2 with other unbiased estimators of Z (without jeopardizing the ergodicity of the chain). For instance, instead of recycling the samples generated in the same iteration as auxiliary points as in Eq. (<ref>),  we could reuse samples generated in the previous iteration t-1. This alternative version of I-MTM method (I-MTM2) is given in Table <ref>. Note that, in both cases I-MTM and I-MTM2, the selected candidate θ^(j) is drawn from the following particle approximation of the target π̅,
 
    π(θ|θ^(1:N))= ∑_i=1^N w̅(θ^(i))δ(θ-θ^(i)),     w̅_i=w̅(θ^(i))=w(θ^(i))/∑_n=1^N w(θ^(n)),

 i.e., θ^(j)∼π(θ|θ^(1:N)).  The acceptance probability α used in I-MTM2 can be also justified considering a proper IS weighting of a resampled particle <cit.> and using the expression (<ref>) related to the standard MH method, as discussed in <cit.>.  Figure <ref> provides a graphical representation of the I-MTM schemes.






  

























  §.§.§ Reusing candidates in parallel I-MTM chains




Let us consider to run C independent parallel chains yielded by an I-MTM scheme. In this case, we have NC evaluations of the target function π and C resampling steps performed at each iteration (so that we have NCT total target evaluations and CT total resampling steps). 

In literature, different authors have suggested to recycle the N candidates, θ^(1),θ^(2),…, θ^(N)∼ q(θ), in order to reduce the number of evaluations of the target pdf <cit.> . The idea is to performs C-times the resampling procedure considering the same set of candidates, θ^(1),θ^(2),…, θ^(N) (a similar approach was proposed in <cit.>). Each resampled candidate is then tested as possible future state of one chain.
 In this scenario, The number of target evaluations per iteration is only N (hence, the total number of evaluation of π is NT). However, the resulting C parallel chains are no longer independent, and there is a lose of performance w.r.t. the independent chains. There exists also the possibility of reducing the total number of resampling steps, as suggested in the Block Independent MTM scheme <cit.> (but  the dependence among the chains grows even more). 






 §.§ Particle Metropolis-Hastings (PMH) method







Assume that the variable of interest is formed by only a dynamical variable, i.e., θ= x=x_1:D=[x_1…,x_D]^⊤ (see Section <ref>). This is the case of inferring a hidden state in state-space model, for instance. More generally, let assume that we are able to factorize the target density as

    π̅( x)∝π( x)   =   γ_1(x_1) γ_2(x_2|x_1) ⋯γ_D(x_D|x_1:D-1), 
       =   γ_1(x_1) ∏_d=2^D  γ_d(x_d|x_d-1).

The Particle Metropolis Hastings (PMH) method <cit.> is an efficient MCMC technique, proposed independently from the MTM algorithm, specifically designed for being applied in this framework. Indeed, we can take advantage of the factorization of the target pdf and consider a proposal pdf decomposed in the same fashion

    q( x)=q_1(x_1) q_2(x_2|x_1) ⋯  q_D(x_D|x_1:D-1)=q_1(x_1) ∏_d=2^D  q_d(x_d|x_d-1).

Then, as in a batch IS scheme, given an n-th sample x^(n)=x_1:D^(n)∼ q( x) with x_d^(n)∼ q_d(x_d|x_d-1), we assign the importance weight  

    w( x^(n))=w_D^(n)=π( x^(n))/q( x^(n))=γ_1(x_1^(n)) γ_2(x_2^(n)|x_1^(n)) ⋯γ_D(x_D^(n)|x_1:D-1^(n))/q_1(x_1^(n)) q_2(x_2^(n)|x_1^(n)) ⋯  q_D(x_D^(n)|x_1:D-1^(n)).

 The previous expression suggests a recursive procedure for computing the importance weights: starting with w_1^(n)=π(x_1^(n))/q(x_1^(n)) and then
 
    w_d^(n)=w_d-1^(n)β_d^(n)=∏_j=1^dβ_j^(n),        d=1,…,D,

 where we have set
 
    β_1^(n)=w_1^(n)    β_d^(n)=γ_d(x_d^(n)|x_1:d-1^(n))/q_d(x_d^(n)|x_1:d-1^(n)),
 
for d=2,…,D. This method is usually referred as Sequential Importance Sampling (SIS).  If resampling steps are also employed at some iteration, the method is called Sequential Importance Resampling (SIR), a.k.a., particle filtering (PF) (see Appendix <ref>). PMH uses a SIR approach for providing the particle approximation π( x| x^(1:N))= ∑_i=1^N w̅_D^(i)δ( x- x^(i)) where w̅_D^(i)=w_D^(i)/∑_n=1^Nw_D^(n) and w_D^(i)=w( x^(i)), obtained using Eq. (<ref>) (with a   proper weighting of a resampled particle <cit.>). Then, one particle is drawn from this approximation, i.e., with a probability proportional to the corresponding normalized weight. 
Estimation of the marginal likelihood Z in particle filtering. SIR combines the SIS approach with the application of resampling procedures. In SIR, a consistent estimator of Z is given by 

    Z=∏_d=1^D[∑_n=1^Nw̅_d-1^(n)β_d^(n)],     w̅_d-1^(i)=w_d-1^(i)/∑_n=1^Nw_d-1^(n).

Due to the application of the resampling, in SIR the standard estimator 

    Z=1/N∑_n=1^N w_D^(n)=1/N∑_n=1^N w( x^(n)),

is a possible alternative only if a proper weighting of the resampled particles is applied <cit.> (otherwise, it is not an estimator of Z).  If a proper weighting of a resampled particle is employed, both Z and Z are equivalent estimators of Z <cit.>. Without the use of resampling steps (i.e., in SIS),  Z and Z are always equivalent estimators <cit.>. See also Appendix <ref>.
The complete description of PMH is provided in Table <ref> considering the use of Z. At each iteration, a particle filter is run in order to provide an approximation by N weighted samples of the measure of the target. Then, a sample among the N weighted particles is chosen by one resampling step. This selected sample is then accepted or rejected as next state of the chain according to an MH-type acceptance probability, which involves two estimators of marginal likelihood Z. PMH is also related to other popular method in molecular simulation called Configurational Bias Monte Carlo (CBMC) <cit.>. 
 
  


















  §.§.§ Relationship among I-MTM2, PMH and I-MH



A simple look at I-MTM2 and PMH shows that they are strictly related <cit.>. Indeed, the structure of the two algorithms coincides. The main difference lies that the candidates in PMH are generated sequentially, using a SIR scheme. If no resampling steps are applied, then I-MTM2 and PMH are exactly the same algorithm, where the candidates are drawn in a batch setting or sequential way. Hence, the application of resampling steps is the main difference between the generation procedures of PMH and I-MTM2. Owing to the use of resampling, the candidates { x^(1),…, x^(N)} proposed by PMH are not independent (differently from I-MTM2). As an example, Figure <ref> shows N=40 particles (with D=10) generated and weighted by SIS and SIR procedures (each path is a generated particle x^(i)=x_1:10^(i)). The generation of correlated samples can be also considered in MTM methods without jeopardizing  the ergodicity of the chain, as simply shown for instance in <cit.>, for instance. Another difference is the use of Z or Z. However, if a proper weighting of a resampled particle is employed, both estimators coincide  <cit.>. Furthermore, both I-MTM2 and PMH can be considered as I-MH schemes where a proper importance sampling weighting of a resampled particle is employed <cit.>. Namely, I-MTM2 and PMH are equivalent to an I-MH technique using the following complete proposal pdf, 

    q(θ)=∫_𝒟^Nπ(θ|θ^(1:N)) [∏_i=1^N q(θ^(i))]dθ^(1:N),

 where π is given in Eq. (<ref>), i.e., θ^(j)∼q(θ), and then considering the generalized (proper) IS weighting, w(θ^(j))=Z^*, w(θ_t-1)=Z_t-1  <cit.>. For further details see Appendix <ref>.   














  §.§.§ Particle Marginal Metropolis-Hastings (PMMH) method




Assume now that the variable of interest if formed by both dynamical and static variables, i.e., θ=[ x, λ]^⊤. For instance, this is the case of inferring both, an hidden state x in state-space model, and static parameters λ of the model. The Particle Marginal Metropolis-Hastings (PMMH) technique is a extension of PMH which addresses this problem.

Let us consider x=x_1:D=[x_1,x_2,…,x_D]∈ℝ^d_x, and an additional model parameter λ∈ℝ^d_λ to be inferred as well (θ=[ x, λ]^⊤∈ℝ^D, with D=d_x+d_λ). Assuming a prior pdf g_λ(λ) over λ, and a factorized complete posterior pdf π̅(θ)=π̅( x,λ), 

    π̅( x,λ)∝π( x,λ)   =   g_λ(λ) π( x|λ),

where π( x|λ)=γ_1(x_1|λ) ∏_d=2^Dγ_d(x_d|x_1:d-1,λ). For a specific value of λ, we can use a particle filter approach, obtaining the approximation π(|λ)=∑_n=1^N w̅_D^(n)δ(-^(n)) and the estimator Z(λ), as described above. The PMMH technique is then summarized in Table <ref>. The pdf q_λ(λ|λ_t-1) denotes the proposal density for generating possible values of λ.
 Observe that, with the specific choice q_λ(λ|λ_t-1)=g_λ(λ), then the acceptance function becomes

    α = min[1, Z(λ^*)/Z(λ_t-1)].

Note also that PMMH  w.r.t. to λ can be interpreted as MH method where an the posterior cannot be evaluated point-wise. Indeed, Z(λ) approximates the marginal likelihood p( y|λ) <cit.>.













 §.§ Group Metropolis Sampling



The  auxiliary weighted samples in the I-MTM schemes (i.e., the N-1 samples drawn at each iteration that are not selected to be compared with the previous state θ_t-1) can be recycled providing a consistent and more efficient estimators <cit.>.

The so-called  Group Metropolis Sampling (GMS) method is shown in Table <ref>.
 GMS yields a sequence of sets of weighted samples 𝒮_t={θ_n,t,ρ_n,t}_n=1^N, for t=1,…,T, where we have denoted with  ρ_n,t the importance weights assigned to the samples θ_n,t (see Figure <ref>). All the samples are then employed for a joint particle approximation of the target. Alternatively, GMS can directly provide an approximation of a specific moment of the target pdf (i.e., given a particular function f). The estimator of this specific moment provided by GMS is 

    I_NT=1/T∑_t=1^T ∑_n=1^N ρ_n,t/∑_i=1^N ρ_i,t f(θ_n,t)=1/T∑_t=1^T I_N^(t).













 
 Unlike in the I-MTM schemes, no resampling steps are performed in GMS. However, we can recover an  I-MTM chain from the GMS output applying one resampling step when 𝒮_t≠𝒮_t-1, i.e.,       

    θ_t= {   θ_t∼∑_n=1^N ρ_n,t/∑_i=1^N ρ_i,tδ(θ-θ_n,t),       𝒮_t≠𝒮_t-1, 
       θ_t-1,                         𝒮_t=  𝒮_t-1,
    .

for t=1,…,T. More specifically, {θ_t}_t=1^T is a Markov chain obtained by one run of an I-MTM2 technique. 

The consistency of the GMS estimators is discussed in Appendix <ref>.

GMS can be also interpreted as an iterative IS scheme where an IS approximation of N samples is built at each iteration and compared with the previous IS approximation. This procedure is iterated T times and all the accepted IS estimators I_N^(t) are finally combined to provide a unique global approximation of NT samples. Note that the temporal combination of the IS estimators is obtained dynamically by the random repetitions due to the rejections in the acceptance test. 

The complete weighting procedure in GMS can be interpreted as the composition of two weighting schemes: (a) by an IS approach building {ρ_n,t}_n=1^N and (b) by the possible random repetitions due to the rejections in the acceptance test.
Figure <ref> depicts a graphical representation of the GMS outputs as chain of sets 𝒮_t={θ_n,t,ρ_n,t}_n=1^N.







 §.§ Ensemble MCMC algorithms




Another alternative procedure, often referred as Ensemble MCMC (EnMCMC) methods (a.k.a., called Locally weighted MCMC), involving several tries at each iteration <cit.>. Related techniques has been proposed independently in different works <cit.>. 
First, let us define the joint proposal density

    q(θ^(1),…,θ^(N)|θ_t): 𝒟^N→ℝ,

and, considering N+1 possible elements, 𝒮={θ^(1), …, θ^(N),θ^(N+1)}, we define the D× N matrix

    Θ_ k=[θ^(1),…,θ^(k-1),θ^(k+1),…,θ^(N+1)],

with columns all the vectors in 𝒮 with the exception of θ^(k). For simplicity, in the followings we abuse of the notation writing q(θ^(1),…,θ^(N)|θ_t)=q(Θ_ N+1|θ_t), for instance. One simple example of joint proposal pdf is

    q(θ^(1),…,θ^(N)|θ_t)=∏_n=1^N q(θ^(n)|θ_t),

i.e., considering independence among θ^(n)'s (and having the same marginal proposal pdf q). More sophisticated joint proposal densities can be employed. A generic EnMCMC algorithm is outlined in Table <ref>.


































Note that with respect to the generic MTM method, EnMCMC does not require to draw auxiliary samples and weights. Therefore, in EnMCMC a smaller number of evaluation of target is required w.r.t.  a generic MTM scheme. 






  §.§.§ Independent Ensemble MCMC





In this section, we present an interesting special case, which employs a single proposal pdf q(θ) independent on the previous state of the chain, i.e.,

    q(θ^(1),…,θ^(N)|θ_t)=q(θ^(1),…,θ^(N))=
    ∏_n=1^N q(θ^(n)).

In this case, the technique can be simplified as shown below. At each iteration, the algorithm described in Table <ref> generates N new samples θ^(1),θ^(2),…, θ^(N)  and then resample the new state θ_t within a set of N+1 samples, {θ^(1),…,θ^(N),θ_t+1} (which includes the previous state), according to the probabilities
 
    α(θ_t-1,θ^(j))=w̅_j=w(θ^(j))/∑_i=1^N w(θ^(i))+ w(θ_t-1),      j=1,…,N+1,

where w(θ)=π(θ)/q(θ) denotes the importance sampling weight. Note that Eq. (<ref>) for N=1 becomes 

    α(θ_t-1,θ^(j))   =   w(θ^(j))/ w(θ^(i))+ w(θ_t-1), 
       =   π(θ^(j))/q(θ^(j))/π(θ^(j))/q(θ^(j))+ π(θ_t-1)/q(θ_t-1),
       =   π(θ^(j))q(θ_t-1)/π(θ^(j))q(θ_t-1)+ π(θ_t-1)q(θ^(j)),

that is the Barker's acceptance function (see <cit.>).





As discussed in <cit.>, <cit.>, <cit.>, the density of a resampled candidate becomes closer and closer to π̅ as N grows, i.e., N→∞.
Hence, the performance of  I-EnMCMC clearly improves with N→∞ (see Appendix <ref>).
 The I-EnMCMC algorithm produces an ergodic chain with invariant density π̅, by resampling N+1 samples  at each iteration (N new samples θ^(1),…, θ^(N) from q and setting θ^(N+1)=θ_t-1). Figure <ref> summarizes the steps of I-EnMCMC.













 §.§ Delayed Rejection Metropolis (DRM) Sampling




An alternative use of different candidates in one iteration of a Metropolis-type method is given in <cit.>. 
The idea behind the proposed algorithm, called Delayed Rejection Metropolis (DRM) algorithm, is the following. As in a standard MH method, at each iteration, one sample is proposed θ^(1)∼ q_1(θ|θ_t-1) and accepted with probability

    α_1(θ_t-1,θ^(1))= min(1,π(θ^(1))q_1(θ_t-1|θ^(1))/π(θ_t-1) q_1(θ^(1)|θ_t-1)).

If θ^(1) is accepted then θ_t=θ^(1) and the chain is moved forward. If  θ^(1) is rejected, the DRM method suggests of drawing another samples  θ^(2)∼ q_2(θ|θ^(1),θ_t-1) (considering a different proposal pdf q_2 taking into account possibly the previous candidate θ^(1)) and accepted with a suitable acceptance probability

    α_2(θ_t-1,θ^(2))= 
            =min(1,π(θ^(2)) q_1(θ^(1)|θ^(2))q_2(θ_t-1|θ^(1),θ^(2))(1-α_1(θ^(2),θ^(1)))  /π(θ_t-1) q_1(θ^(1)|θ_t-1) q_2(θ^(2)|θ^(1),θ_t-1)(1-α_1(θ_t-1,θ^(1))) ).
 
      
 The acceptance function α_2(θ_t-1,θ^(2)) is designed in order to ensure the ergodicity of the chain. If θ^(2) is rejected we can set θ_t=θ_t-1 and perform another iteration of the algorithm, or continue with this iterative strategy  drawing θ^(3)∼ q_3(θ|θ^(2),θ^(1),θ_t-1) and test it with a proper probability α_3(θ_t-1,θ^(3)).

  The DRM algorithm with only 2 acceptance stages is outlined in Table <ref> and summarized in Figure <ref>.
   
   Note that the proposal pdf can be improved at each intermediate stage (θ^(1)∼ q_1(θ|θ_t-1), θ^(2)∼ q_2(θ|θ^(1),θ_t-1) etc.), using the information provided by the previous generated samples and the corresponding target evaluations.

The idea behind DRM of creating a path of intermediate points, then improving the proposal pdf, and hence fostering larger jumps have been also considered in other works <cit.>.   












§ SUMMARY: COMPUTATIONAL COST, DIFFERENCES AND CONNECTIONS



The performance of the algorithms described above improves as N grows, in general: the correlation among samples vanishes to zero, and the acceptance rate of new state approaches one (see Section <ref>).[Generally, an acceptance rate close to 1 is not an evidence of good performance for an MCMC algorithm. However, for the techniques tackled in this work, the situation is different: as N grows, the procedure used for proposing a novel possible state (involving N tries, resampling steps etc.) becomes better and better, yielding a better approximation of the target pdf. See Appendix <ref> for further details.] As N increases, they become similar and similar to an  exact sampler drawing independent samples directly from the target density (for MTM, PMH and EnMCMC schemes the explanation is given in Appendix <ref>). However, this occurs at the expense of an additional computational cost. 
 
In Table <ref>, we summarize the total number of target evaluations, E, and the total number of samples 
used in the final estimators, Q (without considering to remove any burn-in period). The generic MTM algorithm has the greatest number of target evaluations. However, a random-walk proposal pdf can be used in a generic MTM algorithm  and, in general, it fosters the exploration of the state space.  In this sense, the generic EnMCMC seems to be preferable w.r.t. MTM, since E=NT and the random-walk proposal can be applied. A disadvantage of the EnMCMC schemes is that their acceptance function seems worse in terms of Peskun's ordering <cit.> (see numerical results in Section <ref>). Namely,  fixing the number of N tries, the target π the proposal q pdfs, the MTM schemes seem to provide greater acceptance rates than the corresponding EnMCMC techniques. This is theoretically proved for N=1 <cit.>, and the difference vanishes to zero as N grows.
The GMS technique, like other strategies <cit.>, has been proposed to recycle samples or re-use target evaluations, in order to increase Q (see also Section <ref>).




 In PMH, the components of the different tries are drawn sequentially and they are correlated due to the application of the resampling steps. In DMR, each candidate is drawn in a batch way (all the components jointly) but the different candidates are drawn in a sequential manner (see Figure <ref>), θ^(1) then θ^(2) etc.    
The benefit of this strategy is that the proposal pdf can be improved considering the previous generated tries. Hence, if the proposal takes into account the previous samples, DMR generates correlated candidates as well. The main disadvantage  of DRM is that the implementation for a generic N>2 is not straightforward.


 I-MTM and  I-MTM2 differs for the acceptance function employed.
Furthermore, The main difference between the I-MTM2 and PMH schemes is the use of resampling steps during the generation the different tries. For this reason, the candidates of PMH are correlated (unlike in I-MTM2).
I-MTM2 and PMH can be interpreted as I-MH methods using a sophisticated proposal density q(θ) in Eq. (<ref>), and an extended IS weighting procedure is employed. Note that, indeed, q cannot be evaluated pointwise, hence a standard IS weighting strategy cannot be employed. 










§ NUMERICAL EXPERIMENTS



We test different MCMC using multiple candidates in different numerical experiments. In the first example, an exhaustive  comparison among several techniques with an independent proposal is given. We have considered different number of tries, length of the chain, parameters of the proposal pdfs and also different dimension of the inference problem. In the second numerical simulation,  
we compare different particle methods. The third one regards the hyperparameter selection for a Gaussian Process (GP) regression model. The last two examples are  localization problems in a wireless sensor network (WSN): in the fourth one some parameters of the WSN are also tuned, whereas in last example a real data analysis is performed.  






 §.§ A first comparison of efficiency



In order to compare the performance of different techniques, in this section we consider a multi-modal, multidimensional Gaussian target density. More specifically, we have 
 
    π̅(θ)=1/3∑_i=1^3𝒩(θ|μ_i,Σ_i),     θ∈ℝ^D,

where μ_1=[μ_1,1,…,μ_1,D]^⊤, μ_2=[μ_2,1,…,μ_2,D]^⊤, μ_3=[μ_3,1,…,μ_3,D]^⊤, with μ_1,d=-3, μ_2,d=0, μ_3,d=2 for all d=1,…,D. Moreover, the covariance matrices are diagonal, Σ_i=δ_i  I_D (where I_D is the D× D identity matrix), with δ_i=0.5 for i=1,2,3. Hence, given a random variable Θ∼π̅(θ), we know analytically that E[Θ]=[θ̅_1,…,θ̅_D] ^⊤ with θ̅_d= -1/3 for all d, and {[Θ]}=[ ξ_1,…, ξ_D] ^⊤ with ξ_d= 85/18 for all d=1,…,D.
 
 
 We apply I-MTM, IMTM2 and I-EnMCMC in order to estimate  all the expected values and all the variances of the marginal target pdfs. Namely, for a given dimension D,  we have to estimate all {θ̅_d}_d=1^D and  
 {ξ_d}_d=1^D, hence 2D values. The results are averaged over 3000 independent runs. At each run, we compute an averaged square error obtained in the estimation of the 2D values and then calculate the Mean Square Error (MSE) averaged over the 3000 runs. For all the techniques, we consider a Gaussian proposal density q(θ)= 𝒩(θ|μ,σ^2 I_D) with μ=[μ_1=0,…,μ_D=0]^⊤ (independent from the previous state) and different values of σ are considered.  




  §.§.§ Description of the experiments


We perform several experiments varying the number of tries, N, the length of the generated chain, T,  the dimension of the inference problem, D, and the scale parameter of the proposal pdf, σ. In Figures <ref>(a)-(b)-(c)-(d)-(e), we show the MSE (obtained by different techniques)  as function of N, T, D and σ, respectively. In Figure <ref>(d), we only consider I-MTM with N∈{1,5,100,1000} in order to show the effect of using different tries N in different dimensions D. Note that I-MTM with N=1 coincides with I-MH.  

Let us denote as ϕ(τ) the auto-correlation function of the states of the generated chain. Figures <ref>(a)-(b)-(c) depicts the normalized auto-correlation function ϕ̅(τ)=ϕ(τ)/ϕ(0) (recall that ϕ(0)≥ϕ(τ) for τ≥ 0) at different lags τ=1,2,3, respectively. Furthermore, given the definition of the Effective Sample Size (ESS) <cit.>, 

    ESS=T/1+2∑_τ=1^∞ϕ̅(τ),

 in Figure <ref>(d), we show of the ratio ESS/T (approximated; cutting off the series in the denominator at lag τ=10), as function of N.[Since the MCMC algorithms yield positive correlated sequences of states, we have ESS< T in general.]
 Finally, in Figures <ref>(a)-(b), we provide the Acceptance Rate (AR) of a new state (i.e., the expected number of accepted jumps to a novel state), as function of N and D, respectively. 




  §.§.§ Comment on the results



Figures <ref>(a)-(b), <ref> and <ref>(a), clearly show that the performance improves as N grows, for all the algorithms.  
The MSE values and the correlation decrease, and  the ESS and the AR grow.
I-MTM seems to provide the best performance. Recall that for N=1,  I-EnMCMC becomes an I-MH with Baker's acceptance function and I-MTM becomes the I-MH in Table <ref> <cit.>. For N=1, the results confirm the Peskun's ordering about the acceptance function for a MH method <cit.>. Observing the results, the Peskun's ordering appears valid also for the multiple try case, N>1. 
I-MTM2 seems to have worse performance than I-MTM for all N. With respect to  I-EnMCMC, I-MTM2 performs better for smaller N. The difference among the MSE values obtained by the samplers becomes smaller as N grows, as shown in Figure <ref>(a)-(b) (note that in the first one D=1, in the other D=10, and the range of N is different).
 The comparison among I-MTM, I-MTM2 and I-EnMCMC  seems not to be affected by changing T and σ, as depicted in Figures <ref>(c)-(e). Namely, the MSE values change but the ordering of the methods (e.g., best and worst) seems to depend mainly on N. 
Obviously, for greater D, more tries are required in order to obtain good performance (see Figures <ref>(d) and<ref>(b), for instance). Note that for N→∞,   I-MTM, I-MTM2 and I-EnMCMC perform similarly to an exact sampler drawing T independent samples from π̅(θ): the correlation ϕ(τ) among the samples approaches zero (for all τ), ESS approaches T and AR approaches 1. 



  
 
 



 

 
 





 §.§ Numerical experiment comparing particle schemes



In this section,  in order to clarify the differences between batch and particle schemes, we consider again a multidimensional Gaussian target density, that can be express as

    π̅(θ)=π̅(θ_1,…, θ_D)=∏_d=1^D 𝒩(θ_d|μ_d,σ^2),

with θ=θ_1:D∈ℝ^D, D=10, with μ_1:3=2,  μ_4:7=4, μ_8:10=-1 (shown dashed line in Figure <ref>(d)), and σ=1/2.
The target pdf π̅(θ)=𝒩(θ|μ,σ^2 I_D) (with μ=μ_1:10) is formed by independent Gaussian pieces, and can factored as in Eq. (<ref>).  The factorization allows the use of a sequential proposal procedure and hence particle schemes can be employed, even if there is not correlation between θ_d and θ_d-1.
[For the particle schemes, we consider θ= x as a dynamic parameter (see Section <ref>).]


We apply I-MTM, I-MTM2, PMH, and a variant of PMH, denote as var-PMH which uses the corresponding acceptance probability of I-MTM in Eq. (<ref>) instead of the acceptance function of the classical PMH in Eq. (<ref>). The goal is to estimate the vector μ. We compute the MSE in estimating the vector μ=μ_1:10, averaging the results over 500 independent simulations. The components of μ are shown in Figure <ref>(d) with a dashed line.

For all the techniques, we employ a sequential construction of the N candidates  (using the chain rule, see below): in PMH and var-PMH the resampling is applied at each iteration whereas in I-MTM and I-MTM2 no resampling is applied. More specifically, the proposal density for all the methods is 

    q(θ)=q_1(θ_1)∏_d=2^D q_d(θ_d|θ_d-1),

where q_1(θ_1)=𝒩(θ_1|-2, 4) and q_d(θ_d|θ_d-1)=𝒩(θ_d|θ_d-1, σ_p^2), but PMH and var-PMH employ resampling steps so that the generated tries are correlated (whereas in  I-MTM and I-MTM2 the generated candidates are independent).





 We test all the techniques considering different value of number of tries N and number of iterations of the chain T. 
 Figures <ref>(a)-(b) show the MSE    
as function of number of iterations T, keeping fixed the number of tries N=3. Figure <ref>(a) reports the results of the MTM schemes whereas Figure <ref>(b) reports the results of the PMH schemes. Figure <ref>(c) depicts the MSE as function of N (with T=2000), for the PMH methods. Note that the use of only N=3 particles and the application of the resampling at each iteration is clearly a disadvantage for the PMH schemes. If the resampling is applied very often (as in this case), a greater number of N is advisable (such as N=100 or N=1000). Hence, the results confirm that applying a resampling step at each iteration is not optimal and that a smaller rate of resampling steps could improve the performance <cit.>. The results also confirm that the use of an acceptance probability of type in Eq. (<ref>) provides smaller MSE, i.e., I-MTM and var-PMH perform better than I-MTM2 and PMH, respectively. This is more evident for small number of candidates N. 
 When N grows,  the performance of  PMH  and var-PMH methods becomes similar, since the acceptance probability approaches 1, in both cases.  Figure <ref>(d) depicts 35 different states θ_t=θ_1:10,t at different iteration indices t, obtained with var-PMH (N=1000 and T=1000) and the values μ_1:10 are given in dashed line.


 





 §.§ Hyperparameter tuning for Gaussian Process (GP) regression models




 
Let us assume observed data pairs {y_j, z_j}_j=1^P, with y_j∈ℝ and z_j∈ℝ^L. We also denote the corresponding P× 1 output vector as y=[y_1,…,y_P]^⊤ and the L× P input matrix as Z=[ z_1,…, z_P]. 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 Gaussian Process (GP) <cit.>. Hence f( z) ∼𝒢𝒫(μ( z),κ( z, r)) where μ( z)=0, z, r∈ℝ^L, and we consider the kernel function 

    κ( z, r)=exp(-∑_ℓ=1^L(z_ℓ-r_ℓ)^2/2δ^2).

Given these assumptions, the vector f=[f( z_1),…, f( z_P)]^⊤ is distributed as 
p(  f|  Z,δ, κ)=𝒩(  f;  0,  K),
where 0 is a P× 1 null vector, and K_ij:=κ( z_i, z_j), for all i,j=1,…,P, is a P× P matrix.
The vector containing all the hyperparameters of the model is
θ=[δ, σ],
i.e., all the parameters of the kernel function in Eq. (<ref>) and standard deviation σ of the observation noise. In this experiment, we focus on the marginal posterior density of the hyperparameters,
π̅(θ|  y,   Z, κ)∝π(θ|  y,   Z, κ)=p(  y|θ,   Z, κ) p(θ),
which can be evaluated analytically, but we cannot compute integrals involving it <cit.>. Considering a uniform prior within [0,20]^2, p() and since p(  y|θ,   Z, κ)=𝒩(  y;  0, K+σ^2   I), we have 

    log[π(θ|  y,   Z, κ)] = -1/2  y^⊤ (  K+σ^2   I)^-1  y -1/2log[(  K +σ^2   I)]+C,

 where C>0, and clearly K depends on δ <cit.>. The moments of this marginal posterior cannot be computed analytically. Then, in order to compute the Minimum Mean Square Error (MMSE) estimator θ=[δ,σ], i.e., the expected value E[Θ] with 
Θ∼π̅(θ|  y,   Z, κ), we approximate E[Θ] via Monte Carlo quadrature. More specifically, we apply I-MTM2, GMS, 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, {y_j, z_j}_j=1^P, according to the GP model above setting δ^*=3, σ^*=10, L=1, and drawing z_j∼𝒰([0,10]). We keep fixed these data over the different runs. We computed the ground-truth θ= [δ= 3.5200,σ= 9.2811] using an exhaustive and costly grid approximation, in order to compare the different techniques. 

For I-MTM2, GMS, and MH schemes, we consider the same adaptive Gaussian proposal pdf q_t(θ|μ_t,λ^2  I)=𝒩(θ|μ_t,λ^2  I), with λ=5 and μ_t is adapted considering the arithmetic mean of the outputs after a training period, t≥ 0.2 T, in the same fashion of <cit.> (μ_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 <ref> (log-log plot) shows the Mean Square Error (MSE) in the approximation of θ averaged over 10^3 independent runs. Observe that GMS always outperforms the corresponding I-MTM2 scheme. These results confirm the advantage of recycling the auxiliary samples drawn at each iteration during an I-MTM2 run.  In Figure <ref>, we show the MSE obtained by GMS keeping invariant the number of target evaluations E=NT=10^3 and varying N∈{1,2,10,20,50,100,250,10^3}. As a consequence, we have T∈{10^3,500,100,50,20,10,4,1}. Note that the case  N=1, T=10^3, corresponds to an adaptive MH (A-MH) method with a longer chain, whereas the case N=10^3, T=1, corresponds to a static IS scheme (both with the same posterior evaluations E=NT=10^3). 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}.










 §.§ Localization of a target in a wireless sensor network


We consider the problem of positioning a target in ℝ^2 using a range measurements in a wireless sensor network (WSN) <cit.>.  
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 scenarios. Indeed, 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=[Z_1,Z_2]^⊤. The position of the target is then a specific realization Z= z. 
The range measurements are obtained from N_S=6 sensors located at h_1=[3, -8]^⊤, h_2=[8,10]^⊤, h_3=[-4,-6]^⊤, h_4=[-8,1]^⊤, h_5=[10,0]^⊤ and h_6=[0,10]^⊤ as shown in Figure <ref>(a). The observation models are given by

    Y_j=20log(|| z- h_j ||)+B_j,    j=1,…, N_S,
   
where B_j are independent Gaussian random variables with pdfs, 𝒩(b_j;0,ζ_j^2), j=1,…, N_S.  We denote 
ζ=[ζ_1,…,ζ_N_S] the vector of standard deviations.

Given the position of the target z^*=[z_1^*=2.5,z_2^*=2.5]^⊤ and setting ζ^*=[ζ_1^*=1,ζ_2^*=2,ζ_3^*=1,ζ_4^*=0.5,ζ_5^*=3,ζ_6^*=0.2] (since N_S=6 and D=N_S+2=8), we generate N_O=20 observations from each sensor according to the model in Eq. (<ref>).
 Then, we finally obtain a measurement matrix Y=[y_k,1,…, y_k,N_S] ∈ℝ^d_Y, where d_Y=N_ON_S=120, k=1,…,N_O.  

We consider uniform prior 𝒰(ℛ_z) over the position [z_1,z_2]^⊤  with ℛ_z=[-30× 30]^2, and a uniform prior over ζ_j, so that ζ has prior 𝒰(ℛ_ζ) with ℛ_ζ=[0,20]^N_S. Thus, the posterior pdf is 

    π̅(θ|Y)=π̅( z,ζ|Y)= ℓ(y|z_1,z_2,ζ_1,…,ζ_N_S)∏_i=1^2p(z_i)∏_j=1^N_S p(ζ_j), 
       =   [∏_k=1^N_O∏_j=1^N_S1/√(2πζ_j^2)exp(-1/2ζ_j^2(y_k,j+10log(|| z- h_j ||)^2) ]𝕀_z(ℛ_z) I_ζ(ℛ_ζ),

where θ=[ z,ζ]^⊤ is a vector of parameters of dimension D=N_S+2=8 that we desire to infer, and 𝕀_c(ℛ) is an indicator variable that is 1 if c∈ℛ, otherwise is 0. 

Our goal is to compute the Minimum Mean Square Error (MMSE) estimator, i.e., the expected value of the posterior π̅(θ|Y)=π̅( z,ζ|Y) (recall that D=8). 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 <cit.>, 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 q_t(θ|μ_n,t,σ^2 I)=𝒩(θ|μ_t,σ^2 I) which is adapted considering the empirical mean of the generated samples after a training period, t≥ 0.2 T <cit.>, μ_0∼𝒰([1,5]^D) and σ=1. For AMIS, we have q_t(θ|μ_t, C_t)=𝒩(θ|μ_t, C_t), where μ_t is as previously described (with μ_0∼𝒰([1,5]^D)) and  C_t is also adapted using the empirical covariance matrix, starting C_0= 4   I. 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, q_n(μ_n,t|μ_n,t-1,σ^2  I)=𝒩(μ_n,t|μ_n,t-1,σ^2  I) with μ_n,0∼𝒰([1,5]^D) for all n and σ=1. 
 


We fix the total number of evaluations of the posterior density as E=NT=10^4. 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 <cit.>).  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  θ^*=[ 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 <ref>, <ref>, and <ref> and Figure <ref>(b). Note that GMS outperforms AMIS for each a pair {N,T} (keeping fixed E=NT=10^4), and GMS also provides smaller MSE values than N parallel MH chains (the case N=1 corresponds to a unique longer chain). Figure <ref>(b) shows the MSE versus N maintaining E=NT=10^4 for GMS and the corresponding MTM method. This figure again confirms the advantage of recycling the samples in a MTM scheme.



















 §.§ Localization with real data


In this section, we describe a numerical experiment involving real data. More specifically, we consider a localization problem <cit.>. We have carried out an experiment with a network consisting of four nodes. Three of them are placed at fixed positions and play the role of sensors that measure the strength of the radio signals transmitted by the target. The other node plays the role of the target to be localized. All nodes are bluetooth devices (Conceptronic CBT200U2A) with a nominal maximum range of 200 m. 
We consider a square monitored area of 4 × 4 m and place the sensors at fixed positions h_1=[0.5,1], h_2=[3.5,1] and h_3=[2,3], with all coordinates in meters. The target is located at z=[z_1=2.5,z_2=2].  The measurement provided by the i-th sensor is denoted as a random variable Y_i, considering the following model 

    Y_i=κ-20log[||z-h_i||]+B_i,

where B_i are again independent Gaussian random variables with pdfs 𝒩(b_i;0,ζ^2), for all i=1,2,3. Differently from the previous section, we estimate in advance the following parameters of the model, κ̂≈ -26.58 and ζ̂≈ 4.73, using a least square fitting. We obtain N_O=5 measurements from each sensor (d_Y=3N_O=15),
and we consider a uniform prior on the 4 × 4 m  area. Given these measurements, we approximate the expected value E[ Z] of the corresponding posterior π̅( z) (here θ= z) using a thin deterministic bivariate grid, obtaining the ground truth ≈ [3.17,2.62]^⊤. 
We test an MH method and MTM scheme using with a random walk Gaussian proposal pdf, q( z|ẕ_t-1)= 𝒩( z| z_t-1,σ^2 I_2), with σ=1, T∈{1000,5000}, and N∈{10, 100, 1000} tries for MTM (clearly, N=1 for MH).  We also test a Metropolis-adjusted Langevin algorithm (MALA), where the proposal is Gaussian random walk density with mean z_t-1+β∇[logπ( z)] and  ∇[logπ( z)] denotes the gradient of logπ( z) <cit.>.  The covariance matrix of MALA Gaussian proposal is  σ^2 I_2 (as the other techniques) and the drift parameter β=σ^2/2. We compute the MSE in estimating E[ Z]≈ [3.17,2.62]^⊤ and averaged the results over 2000 independent runs (at each run, we take the mean of the  square error values of each component). The results are shown in Table <ref>. Recall that MALA uses the additional information of the gradient. Note that, a MALA-type proposal pdf can also be used in a MTM scheme. 
 The use of multiple tries improves the mixing of the Markov chain and speeds up the convergence.









§ CONCLUSIONS



We have provided a thorough review of MCMC methods using multiple candidates in order to select the next state of the chain. We have presented and compared different Multiple Try Metropolis, Ensemble MCMC and Delayed Rejection Metropolis schemes. We have also described the Group Metropolis Sampling technique which generates a chain of set of weighted samples, so that some candidates are properly reused in the final estimators. Furthermore, we have shown how the Particle Metropolis-Hastings algorithm can be interpreted as an MTM scheme using a particle filter for generating the different weighted candidates. Several connections and differences have been pointed out. 

Finally, we have tested several techniques in different numerical experiments: two toy examples in order to provide an exhaustive comparison among the methods, a numerical example regarding the hyperparameter selection for a Gaussian Process (GP) regression model, and two localization problems, one of them involving a real data analysis.  












§ ACKNOWLEDGEMENTS


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



















§ DISTRIBUTION AFTER RESAMPLING




Let us also denote as θ∈{θ^(1)…,θ^(N)}, a generic sample after applying one multinomial resampling step according to the normalized IS weights w̅_n, n=1,…,N.

The density of θ is given by






    q(θ)=∫_𝒟^Nπ(θ|θ^(1:N)) [∏_i=1^N q(θ^(i))]dθ^(1:N),

where 

    π(θ|θ^(1:N))=∑_j=1^N w̅_j δ(θ-θ^(j)).
 
We also define also the matrix
m_ n=[θ^(1),…, θ^(n),θ^(n+1),…, θ^(N)], 
containing all the samples except for the n-th.
After some straightforward rearrangements,  Eq. (<ref>) can be rewritten as

    q(θ)= ∑_j=1^N(∫_𝒟^N-1π(θ)/∑_n=1^Nπ(θ^(n))/q(θ^(n))[∏_n=1 
     n ≠ j^Nq(θ^(n))]
    		 d m_ j).

Finally, we can write

    q(θ) = π(θ) ∑_j=1^N∫_𝒟^N-11/NZ[∏_n=1 
     n ≠ j^N q(θ^(n))]  d m_ j,

where Z= 1/N∑_n=1^Nπ(θ^(n))/q(θ^(n)) that is that  IS estimator of Z. The equation above represents the density of a resampled particle θ∈{θ^(1)…,θ^(N)}. Note that if Z=Z then q(θ) =π(θ). Clearly, for a finite value of N, there exists a discrepancy between q(θ) and π̅(θ), but this discrepancy decreases as N grows.




§ PARTICLE FILTERING



In this appendix, we recall some basic concepts and more recent results about the Sequential Importance Sampling (SIS) and Sequential Importance Resampling (SIR) methods, important for the description of the techniques above. 
SIS procedure. We assume  that the variable of interest is formed by only a dynamical variable, i.e., θ= x=x_1:D=[x_1…,x_D]^⊤ (for simplicity, consider x_d∈ℝ) and  the target can be factorized as

    π̅( x)∝π( x)   =   γ_1(x_1) ∏_d=2^D  γ_d(x_d|x_d-1).

Given a proposal of type q( x)=q_1(x_1) ∏_d=2^D  q_d(x_d|x_d-1), and a sample x^(n)=x_1:D^(n)∼ q( x) with x_d^(n)∼ q_d(x_d|x_d-1), we assign the importance weight  

    w( x^(n))=w_D^(n)=π( x^(n))/q( x^(n))=γ_1(x_1^(n)) γ_2(x_2^(n)|x_1^(n)) ⋯γ_D(x_D^(n)|x_1:D-1^(n))/q_1(x_1^(n)) q_2(x_2^(n)|x_1^(n)) ⋯  q_D(x_D^(n)|x_1:D-1^(n)).

 The weight above can be compute with a recursive procedure for computing the importance weights: starting with w_1^(n)=π(x_1^(n))/q(x_1^(n)) and then
 
    w_d^(n)=w_d-1^(n)β_d^(n)=∏_j=1^dβ_j^(n),        d=1,…,D,

 where we have set
 
    β_1^(n)=w_1^(n)    β_d^(n)=γ_d(x_d^(n)|x_1:d-1^(n))/q_d(x_d^(n)|x_1:d-1^(n)),
 
for d=2,…,D.  Let also define the partial target pdfs

    π_d(x_1:d)=γ_1(x_1) ∏_i=2^d  γ_i(x_i|x_i-1)












SIR procedure. In SIR, a.k.a., standard particle filtering, resampling steps are incorporated during the recursion as shown of Table <ref> <cit.>. In general, the resampling steps are applied only in certain iterations in order to avoid the path degeneration, taking into account an approximation ESS of the Effective Sampling Size (ESS) <cit.>. If ESS is smaller than a pre-established threshold, the particles are resampled. Two examples of ESS approximation are ESS=1/∑_n=1^N (w̅_d^(n))^2 and ESS=1/maxw̅_d^(n) where w̅_d^(n)=w_d^(n)/∑_i=1^N w_d^(i) (note that 1 ≤ESS≤ N).
Hence, the condition for the adaptive resampling can be expressed as ESS < η N where η∈ [0,1]. 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 <cit.>. If η= 0, no resampling steps are applied, and we have the SIS method described above.


Note that in Table <ref>, we have employed a proper weighting for resampling particles <cit.>,

    w_d^(1)=w_d^(2)=…= w_d^(n)=Z_d.

Generally, it is remarked that w_d^(1)=w_d^(2)=…= w_d^(n) but a specific value is not given. If a different value c≠Z_d is employed, i.e., w_d^(1)=…= w_d^(n)=c, the algorithm is still valid but the weight recursion loses part of the statistical meaning. This is the reason why the marginal likelihood estimator Z=Z_D=1/N∑_n=1^Nw_D^(n) is consistent only, if a proper weighting after resampling is used <cit.>.  















 §.§ Marginal likelihood estimators in SIS


In SIS, there are two possible formulations of the estimator of the marginal likelihoods Z_d =∫_ℝ^dπ_d(x_1:d) d x_1:d, 

    Z_d   =   1/N∑_n=1^N w_d^(n)= 1/N∑_n=1^N w_d-1^(n)β_d^(n),  
    Z_d   =   ∏_j=1^d[∑_n=1^Nw̅_j-1^(n)β_j^(n)].

In SIS, both estimators are equivalent Z_d≡Z_d. 


Indeed, the classical IS estimator of the normalizing constant Z_d at the d-th iteration is

    Z_d   =   1/N∑_n=1^N w_d^(n)=1/N∑_n=1^N w_d-1^(n)β_d^(n),  
       =   1/N∑_n=1^N[∏_j=1^d β_j^(n)].
 

An alternative formulation, denoted as Z_d, is often used  

    Z_d   =   ∏_j=1^d[∑_n=1^Nw̅_j-1^(n)β_j^(n)]  
       =   ∏_j=1^d[∑_n=1^N w_j^(n)/∑_n=1^N w_j-1^(n)]=
    Z_1∏_j=2^d[  Z_j/Z_j-1]=Z_d.
 
where we have employed w̅_j-1^(n)=w_j-1^(n)/∑_i=1^N w_j-1^(i) and w_j^(n)=w_j-1^(n)β_j^(n) <cit.>. 

Furthermore, note that Z_d can be written in a recursive form as

    Z_d=Z_d-1[∑_n=1^Nw̅_d-1^(n)β_d^(n)].
 








 §.§ Marginal likelihood estimators in SIR




If  a proper weighting after resampling is applied in SIR, both formulations Z_d  and Z_d in Eqs. (<ref>)-(<ref>) provide consistent estimator of Z_d and they are equivalent, Z_d≡Z_d (as in SIS). 


If a proper weighting is not applied, only  

    Z_d=∏_j=1^d[∑_n=1^Nw̅_j-1^(n)β_j^(n)].
 
is a  consistent estimator of Z_d, in SIR. In this case, Z_d=1/N∑_n=1^N w_d^(n) is not a possible alternative (without using a proper weighting after resampling). However, considering the proper weighting of the resampled particles, then Z_d is also a consistent estimator of Z_d and it is equivalent to Z_d. Below, we analyze three cases:


  * No Resampling (η=0): this scenario corresponds to SIS where Z_d, Z_d are equivalent as shown in Eq. (<ref>). 

  * Resampling at each iteration (η=1):  using the proper weighting, w_d-1^(n)=Z_d-1 for all n and for all d, and replacing in Eq. (<ref>)  we have  

    Z_d   =   Z_d-1[1/N∑_n=1^Nβ_d^(n)], 
       =   1/N∏_j=1^d [ ∑_n=1^Nβ_j^(n)].

Since after resampling all particles have the same weight, we have w̅_d-1^(n)=1/N for all n. Replacing it in the expression of Z_d in (<ref>), we obtain

    Z_d=1/N∏_j=1^d [  ∑_n=1^Nβ_j^(n)],

that coincides with Z_d in Eq. (<ref>).

  * 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. (<ref>), we have Z_k≡Z_k,[We consider to compute the estimators before the resampling.] since before the k-th iteration no resampling has been applied. 
With the proper weighting w_k^(n)=Z_k for all n, at the next iteration we have

    Z_k+1   =   1/N∑_n=1^N w_k^(n)β_k+1^(n)=Z_k[1/N∑_n=1^N β_k+1^(n)],

and using Eq. (<ref>),  we obtain 

    Z_k+1   =   Z_k[∑_n=1^N1/Nβ_k+1^(n)]=Z_k[1/N∑_n=1^N β_k+1^(n)],

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

 





§ CONSISTENCY OF GMS ESTIMATORS




Dynamic of GMS.

We have already seen that we can recover an I-MTM chain from the GMS outputs applying one resampling step for each t when 𝒮_t≠𝒮_t-1, i.e.,       

    θ_t= {   θ_t∼∑_n=1^N ρ_n,t/∑_i=1^N ρ_i,tδ(θ-θ_n,t),       𝒮_t≠𝒮_t-1, 
       θ_t-1,                         𝒮_t=  𝒮_t-1,
    .

for t=1,…,T. The sequence {θ_t}_t=1^T is a chain obtained by one run of an I-MTM2 technique. Note that 
(a) the sample generation, (b) the acceptance probability function and hence (c) the dynamics of GMS exactly coincide with the corresponding steps of I-MTM2 (or PMH; depending on candidate generation procedure). Hence, the ergodicity of the recovered chain is ensured.

Parallel chains from GMS outputs. 

As described in Section <ref>, we can extend the consideration above for generation C parallel I-MTM2 chains. Indeed, we resample C times instead of only one, i.e.,  

    θ_t^(c)= {   θ_t∼∑_n=1^N ρ_n,t/∑_i=1^N ρ_i,tδ(θ-θ_n,t),       𝒮_t≠𝒮_t-1, 
       θ_t-1,                         𝒮_t=  𝒮_t-1, .

for c=1,…,C, where the super-index denotes the c-th chain (similar procedures have been suggested in <cit.>).
Clearly, the resulting C parallel chains are not independent, and there is  an evident loss of performance w.r.t. the case of independent chains. However, at each iteration, the number of target evaluations per iteration is only N instead of NC. 
Note that that each chain in ergodic, so that each estimator I_T^(c)=1/T∑_t=1^Tg(θ_t^(c)) is consistent (i.e., convergence to the true value for T→∞). As a consequence, the arithmetic mean of consistent estimators,

    I_C,T=1/C∑_c=1^CI_T^(c)=1/CT∑_t=1^T  ∑_c=1^C g(θ_t^(c)),

is also consistent, for all values of C≥ 1. 

GMS as limit case. 

Let us consider the case 𝒮_t≠𝒮_t-1 (the other one is trivial), at some iteration t. In this scenario, the samples of the C parallel I-MTM2 chains, θ_t^(1),θ_t^(2),...,θ_t^(C), are obtained by resampled independently C samples from the set {θ_1,t,…,θ_N,t} according to the normalized weights  ρ̅_n,t=ρ_n,t/∑_i=1^N ρ_i,t, for n=1,…,N. Recall that the samples θ_t^(1),θ_t^(2),...,θ_t^(C), will be used in the final estimator I_C,T in Eq. (<ref>). 


Let us denote as #j the number of times that a specific candidate θ_j,t (contained in the set {θ_n,t}_n=1^N) has been selected as state of one of C chains, at the t iteration. As C→∞, The fraction  #j/C approaches exactly the corresponding weights ρ̅_j,t. Then, for C→∞, we have that the estimator in Eq. (<ref>) approaches the GMS estimator, i.e.,

    lim_C→∞I_C,T= 1/T∑_t=1^T ∑_n=1^N ρ̅_n,t g(θ_n,t).

 Since I_C,T as T→∞ is consistent for all values of C, then the GMS estimator is also consistent (and it can be obtained as C→∞).
 





