Roland Lamberti

Yohan Petetin

Senior Member, IEEE François Desbouvries

François Septier

Semi-independent resampling for particle filtering

Among Sequential Monte Carlo (SMC) methods, Sampling Importance Resampling (SIR) algorithms are based on Importance Sampling (IS) and on some (resampling-based) rejuvenation algorithm which aims at fighting against weight degeneracy. However this mechanism tends to be insufficient when applied to informative or high-dimensional models. In this paper we revisit the rejuvenation mechanism and propose a class of parameterized SIR-based solutions which enable to adjust the tradeoff between computational cost and statistical performances.

I. INTRODUCTION AND BACKGROUND

Bayesian filtering consists in estimating some variable x t from noisy measurements y 0:t = {y 0 , • • • , y t }. We assume that {(x t , y t )} t≥0 is a Hidden Markov Chain, i.e. that the joint density of (x 0:t , y 0:t ) reads p(x 0:t , y 0:t ) = p(x 0 ) t s=1 f s (x s |x s-1 ) t s=0 g s (y s |x s ). The problem can be traced back to Kalman [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] in the context of linear and Gaussian state space models. Approximate solutions for non linear and/or non Gaussian state space models include the extended Kalman filter [START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF]- [START_REF] Ristic | Beyond the Kalman Filter: Particle Filters for Tracking Applications[END_REF], the unscented Kalman filter [START_REF] Julier | Unscented filtering and nonlinear estimation[END_REF]- [START_REF] Menegaz | A systematization of the unscented Kalman filter theory[END_REF], or SMC methods (also called particle filters (PF)) [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear / non-Gaussian Bayesian tracking[END_REF], [START_REF] Gordon | Novel approach to nonlinear/ non-Gaussian Bayesian state estimation[END_REF], [START_REF]Sequential Monte Carlo Methods in Practice[END_REF], which propagate in time a discrete approximation p(x t |y 0:t ) = N i=1 w i t δ x i t of the posterior pdf p(x t |y 0:t ).

A. The classical SIR algorithm

Let Θ t = ϕ(x t )p(x t |y 0:t )dx t be a moment of interest of p(x t |y 0:t ). One iteration of an SMC algorithm can be decomposed in three steps.

Starting at time t -1 from {w i t-1 , x i t-1 } N i=1 , the first two steps consist in sampling (S.) N particles xi t from importance densities q i and weighting (W.) them so as to take into account the discrepancy between the target and importance densities; then Θ t is estimated as Θ SIS,N t = N i=1 wi t ϕ(x i t ) (superscript SIS will be justified below). Finally a third (optional) step consists in re-sampling (R.) the weighted particles, i.e. in redrawing each particle with a probability equal to its weight and assigning to the resampled particles the same weight 1 N . This yields the class of SIR algorithms [START_REF] Smith | Bayesian statistics without tears : a sampling-resampling perspective[END_REF] [START_REF] Gordon | Novel approach to nonlinear/ non-Gaussian Bayesian state estimation[END_REF] [10] [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear / non-Gaussian Bayesian tracking[END_REF] described by Algorithm 1.

Let us comment this algorithm. If resampling is totally absent, each time iteration reduces to the first two steps, i.e. is based on IS only. However such a sequential IS (SIS) algorithm is well known to fail in practice since after a few iterations most weights get close to zero. The third step (which can be performed whatever t or depending on some criterion such as the number of effective particles [START_REF] Kong | Sequential imputations and bayesian missing data problems[END_REF] [13] [14] [START_REF] Cornebise | Adaptive methods for sequential importance sampling with application to state-space models[END_REF])

Algorithm 1 The classical SIR algorithm Data: q(x t |x t-1 ), y t , {w i t-1 , x i t-1 } N i=1 for 1 ≤ i ≤ N do S. xi t ∼ q(x t |x i t-1 ); W. wi t ∝ w i t-1 ft(x i t |x i t-1 )gt(yt|x i t ) q(x i t |x i t-1 )
,

N i=1 wi t = 1; end for Θ SIS,N t = N i=1 wi t ϕ(x i t ); if R. then for 1 ≤ i ≤ N do l i ∼ Pr(L = l) = wl t , 1 ≤ l ≤ N ; end for Set {w i t , x i t } N i=1 = { 1 N , xli t } N i=1 . else Set {w i t , x i t } N i=1 = { wi t , xi t } N i=1 . end if
discards particles with low weights (such particles are likely to be never resampled) and is considered as a traditional rescue against weight degeneracy. On the other hand, this (R.) step introduces local extra variance [13, section 4.2.1], [16, p. 213], which in turn affects the variance of Θ SIS,N t at subsequent iterations. It has thus been proposed to control this extra variance term via alternative resampling schemes (see e.g. [17] [18] [19] • • • ). Yet despite many proposed refinements this generic SIR mechanism remains inefficient in informative models featuring very sharp likelihood functions (i.e., when g t (y t |x t ) is very small for most values of x t ), and in particular in high-dimensional state-space models [START_REF] Snyder | Obstacles to high-dimensional particle filtering[END_REF], [START_REF] Snyder | Particle filters, the optimal proposal and high-dimensional systems[END_REF].

B. The independent SIR algorithm

Recently it has thus been proposed to revisit the SIR algorithm [START_REF] Lamberti | Particle filters with independent resampling[END_REF] [23] [START_REF] Lamberti | Independent resampling sequential Monte Carlo algorithms[END_REF] and more precisely to come back to the rejuvenation mechanism (R.). The counterpart of this (R.) step is that it duplicates particles with high weights, which results in support degeneracy. Moreover given {w i t-1 , x i t-1 } N i=1 the samples {x j t } produced by Algorithm 1 are marginally distributed from some compound pdf qN t which takes into account the effects of the three elementary (S.), (W.) and (R.) steps, but are obviousy dependent [START_REF] Lamberti | Particle filters with independent resampling[END_REF] (a single particle can be resampled more than once); by contrast, given

{w i t-1 , x i t-1 } N i=1
the independent SIR Algorithm [START_REF] Lamberti | Particle filters with independent resampling[END_REF] [24] produces N i.i.d. draws from qN t . Note that Algorithm 2 below only decribes the rejuvenation step of the independent SIR algorithm, and replaces the "if R. then" part of Algorithm 1.

C. Scope of the paper

Algorithm 2 has displayed good results in severe situations [START_REF] Lamberti | Particle filters with independent resampling[END_REF] and can be combined with a post-resampling, secondstage reweigthing scheme due to its auxiliary particle filtering Algorithm 2 Indep. SIR algorithm (resampling step only)

Data: q(x t |x t-1 ), y t , {w i t-1 , x i t-1 } N i=1 ; for 1 ≤ j ≤ N do x1,j t ← xj t , w1,j t ← wj t . end for for 1 ≤ i ≤ N do R. l i ∼ Pr(L = l) = wi,l t , 1 ≤ l ≤ N ; Rejuvenation of the support for iteration i + 1 if (i < N ) then for 1 ≤ j ≤ N do xt i+1,j ∼ q(x t |x j t-1 ); w i+1,j t = w j t-1 ft(x i+1,j t |x j t-1 )gt(yt|x i+1,j t ) q(x i+1,j t |x j t-1 )
, end for wi+1,:

t ∝ w i+1,: t , N j=1 wi+1,j t = 1; end if end for Set {w i t , x i t } N i=1 = { 1 N , xi,l i t } N i=1 .
interpretation [23] [24]. However its rejuvenation mechanism involves the sampling of N 2 particles { xt i,j } N i,j=1 (plus N resampling steps), by contrast with the classical SIR algorithm which only samples N intermediate particles {x j t } N j=1 (and is also followed by N resampling steps). One can wonder whether this extra cost is indeed necessary, so the aim of this paper is to design an algorithm which is both efficient (in terms of computational cost) and effective (in terms of statistical results). The rest of this paper is organized as follows. Our algorithm is described in section II. Simulations are displayed in section III, and the paper ends with a conclusion.

II. SEMI-INDEPENDENT RESAMPLING

A. An intermediate resampling scheme

The classical and independent SIR resampling mechanisms can be reconciled in a common framework. In both schemes, one progressively builds N weighted sets x1,:

t , • • • , xN,: t (the N supports) and redraws one sample x i t out of each of them (see figure 1). The difference lies in the way xi,:

t is built from xi-1,:

t : in the classical SIR mechanism, xi,:
t is a copy of xi-1,:

t (so the resampling step amounts to redrawing N samples from the common support x1,: t , see Algorithm 1); in the independent SIR mechanism, a whole new support xi,:

t is drawn at each iteration i. In other words, from a computational point of view both schemes resample N particles from some intermediate set {x i,j t } N i,j=1 , but building that set requires N preliminary independent sampling steps in the classical case, while it requires N 2 independent sampling steps in the independent case.

In this paper we propose a resampling scheme which creates an intermediate set {x i,j t } N i,j=1 with more diversity than in the classical case, but at a reduced sampling cost as compared to the independent case. Starting from xi-1,j t , xi,j t can now either be a copy (to save cost) or a new sample (to enhance diversity). The algorithm is as follows. Fix the number k (with 0 ≤ k ≤ N ) of samples which will be redrawn at each iteration. At step i, uniformly draw a subset m i,

1:k = (m i,1 , • • • , m i,k ) of size    x1,1 t . . . x1,N t    x1,: t →    x2,1 t . . . x2,N t    x2,: t → • • • →    xN,1 t . . . xN,N t    xN,: t ↓ ↓ ↓ x 1 t x 2 t x N t Fig
k out of (1, • • • , N ) (m i,l
are the indices of the particles which will be redrawn). Next xi,j t ∼ q(x t |x j t-1 ) if j ∈ m i,1:k , and xi,j t = xi-1,j t if j / ∈ m i,1:k . Finally observe that the classical (resp. independent) SIR algorithm corresponds to the particular case k = 0 (resp. k = N ). The algorithm is summarized in Algorithm 3 below. Algorithm 3 Semi-ind. SIR algorithm (resampling step only)

Data: q(x t |x t-1 ), y t , {w i t-1 , x i t-1 } N i=1 for 1 ≤ j ≤ N do x1,j t ← xj t , w1,j t ← wj t . end for for 1 ≤ i ≤ N do R. l i ∼ Pr(L = l) = wi,l t , 1 ≤ l ≤ N ; Partial rejuvenation of the support for iteration i + 1 if (i < N ) then xi+1,: t ← xi,: t , w i+1,: t ← w i,: t ; for 1 ≤ j ≤ k do m j ∼ Pr(M = n|n ∈ 1:N \{m 1:j-1 }) = 1 N -j+1 ; xi+1,m j t ∼ q(x t |x m j t-1 ); w i+1,m j t = w m j t-1 ft(x i+1,m j t |x m j t-1 )gt(yt|x i+1,m j t ) q(x i+1,m j t |x m j t-1 )
; end for wi+1,:

t ∝ w i+1,: t , N j=1 wi+1,j t = 1; end if end for Set {w i t , x i t } N i=1 = { 1 N , xi,l i t } N i=1 .

B. Performances vs. computational cost

We now evaluate the performance of this procedure by comparing the variances of the estimates computed after the resampling step because they affect the variances of the estimates at subsequent iterations [START_REF] Douc | Comparison of resampling schemes for particle filtering[END_REF]. So let Θ .,N t = 1 N N i=1 ϕ(x .,i t ), where the generic notation x .,i t represents the points produced either by Algorithm 1, 2 or 3 (so we consider Θ SIR,N t , Θ I-SIR,N t and Θ SR,N,k t , where SR stands for semiresampling). We have the following proposition (the proof is given in the Appendix).

Proposition 1: Given the previous set of particles {x i 0:t-1 } N i=1 , for all k, 0 ≤ k ≤ N , we have:

E( Θ SR,N,k t ) = E( Θ I-SIR,N t ) = E( Θ SIR,N t ), (1a)
var( Θ I-SIR,N t ) ≤ var( Θ SR,N,k t ) ≤ var( Θ SIR,N t ), (1b)
var( Θ SR,N,k t ) ≤ var( Θ SR,N,k-1 t ). (1c)
So as the number k of intermediate redrawings increases from 0 to N , the conditional variance of the semi-independent resampling estimator Θ SR,N,k t decreases from the upper bound of inequality (1b) (if k = 0, Θ SR,N,0 t reduces to Θ SIR,N t ) to its lower bound (if k = N , Θ SR,N,N t reduces to Θ I-SIR,N t ). However remember from section II-A that N + (N -1) × k samples are needed for building Θ SR,N,k t ; so parameter k of the SR scheme enables to fix a compromise between variance reduction and computational budget.

C. A parallelized version

Finally Algorithm 3 can be transformed into a parallelized version, the non-sequential SR (NSSR) algorithm. At iteration i, instead of duplicating the N -k surviving particles from the previous support xi-1,: t (see Fig. 1), we propose to duplicate the Nk surviving particles directly from the initial set x1,: t of particles. The N -1 new supports can thus be produced in parallel, contrary to Algorithm 3 which by nature is sequential. Of course, this procedure alters the diversity of the final set of particles, as is illustrated by the following proposition.

Proposition 2: Let Θ N SSR,k t be the estimate built from the non-sequential semi-independent resampling procedure. Then given the previous set of particles {x i 0:t-1 } N i=1 , for all k, 0 ≤ k ≤ N , we have:

E( Θ NSSR,N,k t ) = E( Θ I-SIR,N t ) = E( Θ SIR,N t ), (2a)
var( Θ I-SIR,N t ) ≤ var( Θ NSSR,N,k t ) ≤ var( Θ SIR,N t ), (2b) var( Θ NSSR,N,k t ) ≤ var( Θ NSSR,N,k-1 t ), (2c)
var( Θ SR,N,k t ) ≤ var( Θ NSSR,N,k t ). (2d)
So we see that var( Θ NSSR,N,k t ) still decreases with k, but is always larger than var( Θ SR,N,k t ). As with Proposition 1, the variance inequalities still rely on Jensen's inequality, and the proof is omitted.

III. SIMULATIONS

We consider a tracking problem based on range-bearing measurements. The hidden state-vector contains the position and velocity of the target in cartesian coordinates, x t = [c x,t , ċx,t , c y,t , ċy,t ] T . We set f t (x t |x t-1 ) = N (x t ; Fx t-1 ; Q),

g t (y t |x t ) = N (y t ; c 2 x,t + c 2 y,t ; arctan cy,t cx,t T ; R), with R = diag(σ 2 ρ , σ 2 θ ), F = I 2 ⊗ 1 1 0 1 , Q = 10 × I 2 ⊗ 1 3 1 2 1 2
1 where ⊗ is the Kronecker product. We set q(x t |x t-1 ) = f t (x t |x t-1 ) and we compare the RMSEs averaged over 1000 MC runs.

A. Variance of SR procedures

We first analyze the behaviour of our algorithms as a function of k. We set N = 100, σ ρ = 0.1 and σ θ = π 1800 ; all MC runs use the same measurements. Fig. 2 displays the RMSE of Θ SR,N,k t , Θ NSSR,N,k t , Θ I-SIR,N t deduced from our resampling schemes and Θ SIS,N t (a resampling step is computed at each time step but the estimate is taken before this step). Of course, the performances of estimates based on the SR procedure improve when k incrases. Even for small values of k, the improvement is significant. It is also interesting to note that Θ SR,N,k t (resp. Θ NSSR,N t ) has the same performance as Θ I-SIR,N t when k ≥ N/2 (resp. k ≥ 4N/5).

B. RMSE in the informative case at equal cost

We now compare our estimates with existing improvements of the PF in informative models. In particular, the PF with MCMC resample move is a popular solution to introduce sample variety after resampling [START_REF] Gilks | Following a moving target -Monte Carlo inference for dynamic Bayesian models[END_REF]. Roughly speaking, the N particles which follow the (R.) step of Algorithm 1 are moved via an MCMC algorithm with k iterations (here an independent Metropolis-Hasting algorithm). Thus, our SR procedure has the same computational cost in terms of sampling steps as the SIR PF with MCMC moves. We also compare our estimates with those based on the classical SIR and I-SIR algorithms but with a given budget of total sampling (sampling + resampling) operations. We thus set N = 100 particles and k = N/2 for the computation of Θ SR,N,k t and the estimate based on the resample move PF, N = 72 for that of Θ I-SIR,N t and N + (N -1)k 2 = 2575 particles for that of Θ SIS,N t . The global sampling cost for all these algorithms is approximately (2N + N k). We also compute Θ NSSR,N,k t with N = 100 and k = 4N/5; its computation does not have the same computational cost but can be parallelized. The results are displayed in Fig. 3.

When the observations are very informative (σ ρ and σ θ are small), the classical solution tends to degenerate (it starts working when (σ ρ , σ θ ) = (0.15, π 1200 )), while our solutions are robust and present better performances. As the variance of the measurement noise increases, the different estimates tend to behave similarly; the classical SIR algorithm performs slightly better, which is not surprising since in this case it no longer suffers from the degeneracy phenomenon and the number of final samples used is far superior to the other solutions. We also observe that the resample move which uses differently the k extra samples does not perform well when compared to the SR procedure in very informative models, and is outperformed by our solutions when the observations are not informative. Finally, our SR algorithm with k = N 2 outperforms the (totally) independent resampling one when the budget is fixed. IV. CONCLUSION In this paper we revisited the resampling step of PF algorithms, and proposed a resampling scheme where each new final particle is resampled from a support which is partially rejuvenated with k new particles. This yields a class of parameterized solutions which encompasses the classical multinomial resampling technique (k = 0) and the independent resampling one (k = N ), enabling to tune the balance between variance and computational cost. Simulations showed that choosing k = N/2 leads to similar performances to the fully independent resampling procedure. Moreover, in very informative models our algorithm is not affected by the degeneration phenomenon, contrary to the classical SIR algorithm.

APPENDIX PROOF OF PROPOSITION 1

Let us consider a PF with resampling at time t. First, (1a) holds because the SIR, I-SIR and SR procedures all produce resampled particles which, given {x i 0:t-1 } N i=1 , are (marginally) sampled from the same distribution qN ; and (1b) is straightforward from (1c) and the fact that SR reduces to SIR (resp. I-SIR) when k = 0 (resp. k = N ). Let us address (1c). Since Θ SR,N,k

t = 1 N N i=1 ϕ(x SR,i t ), given {x i 0:t-1 } N i=1 N 2 var k ( Θ SR,N,k t ) = N i=1 var(ϕ(x i t ))+2 N i1,i2=1 i1<i2 cov k (ϕ(x i1 t ),ϕ(x i2 t ));
here index k in a (co)variance emphasizes the fact that it depends on k. The first term of the r.h.s. is independent of k (and coincides with var( Θ I-SIR,N t

)), so the difference between different values of k stems from the covariance terms. Next cov k (ϕ(x i1 t ), ϕ(

x i2 t )) = E k [ϕ(x i1 t )ϕ(x i2 t )]- E[ϕ(x i1 t )]E[ϕ(x i2 t )],
and again, the second term of the r.h.s. is independent of k. Finally for i

1 < i 2 , E k [ϕ(x i1 t )ϕ(x i2 t )] = E[E[ϕ(x i1 t )ϕ(x i2 t )| xt i1:i2,: ]] = E[E[ϕ(x i1 t )|x i1,: t ]E[ϕ(x i2 t )|x i2,: t ]] = E[ Θ SIS t (x i1,: t ) Θ SIS t (x i2,: t )] = E[E[ Θ SIS t (x
k (ϕ(x i1 t )ϕ(x i2 t )) = 1 (A k N ) i2-i1 m i 2 i 1 +1 (1:k) h(m i2 i1+1 (1:k)). (3)
It remains to compare (3) with the same expression with k ← k -1. We observe that (3) can be rewritten as for all m i2 i1+1 (k), and so that E[h(m i2 i1+1 (1:k))|m i2 i1+1 (1:k-1)] ≤ h(m i2 i1+1 (1:k -1)), whence (1c).

E k (ϕ(x i1 t )ϕ(x i2 t )) = 1 (A k-1 N ) i2-i1 m i 2 i 1 +1 (1:k-1) 1 N -k + 1 m i 2 i 1 +1 (k)

kFig. 2 :

 2 Fig. 2: RMSE as a function of k, tracking model.

5 PF with SR, k = N 2 BootstrapFig. 3 :

 523 Fig. 3: Tracking model, σρ ∈ [0.01, 0.3] and σ θ ∈ [ π 18000 , π 600 ].

i 2 i 1 +1

 21 (1:k) t ]|m i2 i1+1 (1:k)] = E[E 2 [ Θ SIS t (x i1,: t )|x i1,1:N \m i 2 i 1 +1 (1:k) ]|m i2 i1+1 (1:k)] = h(m i2 i1+1 (1:k)) = h(m i2 i1+1 (1:k -1), m i2 i1+1 (k))because given the trajectories from the previous time steps, particles from different supports are all marginally drawn from the same densities. Finally

E

h(m i2 i1+1 ( 1 :

 1 k -1), m i2 i1+1 (k)),where the second line coincides with the conditionnal expectationE[h(m i2 i1+1 (1:k))|m i2 i1+1 (1:k-1)]. Given m i2 i1+1 (1:k -1), the set xi1,1:N\m i 2 i 1 +1 (1:k) is included in xi1,1:N\m i 2 i 1 +1 (1:k-1) ; consequently, the Rao-Blackwell decomposition (E(E 2 (X|Y )) ≤ E(E 2 (X|Y, Z))) ensures that h(m i2 i1+1 (1:k -1), m i2 i1+1 (k)) ≤ h(m i2 i1+1 (1:k -1))

  . 1: The classical, independent and semi-independent resampling mechanisms. Each scheme draws N supports xi,:

	t and redraws one t out of each support. The difference lies in the way xi,: sample x i t is built from xi-1,: t is a copy of xi-1,j in the classical case; is t a new particle in the independent case; and can be either copied of
	redrawn in the intermediate, semi-independent case.

t : xi,j

  N\m i 2 i 1 +1(1:k) are the particles shared by supports xi1,: and xi2,: . Under this conditioning, Θ SIS t (x i1,: t ) and Θ SIS t (x i2,: t ) are independent so the general term is

		i1,: t ) Θ SIS t (x i2,: t )|m i2 i1+1 (1:k)]]
	where m i2 i1+1 (1 : k) represents all the indices redrawn from iterations i 1 + 1 to i 2 (the third equality holds because x i t is resampled from support xi,: t (see Fig. 1), so E(ϕ(x i t )) = Θ SIS,N t (x i,: t ) where Θ SIS,N was defined in section I-A). t The outer expectation in this last expression corresponds
	to a uniformly weighted sum over all possible values of
	m i2 i1+1 (1:k), i.e. over (A k N ) i2-i1 terms where A k N is the number of arrangements of k among N . Given m i2 i1+1 (1:k), the general term of this sum reads
	E[ Θ SIS t (x i1,: t ) Θ SIS t (x i2,: t )|m i2 i1+1 (1:k)] =
	E[E[ Θ SIS t (x i1,: t ) Θ SIS t (x i2,: t )|x i1,1:N \m t	i 2 i 1 +1 (1:k)	]|m i2 i1+1 (1:k)]
	where xi1,1:E[E[ Θ SIS t (x i1,: t )|x i1,1:N \m t	i 2 i 1 +1 (1:k)	]
	× E[ Θ SIS t (x i2,: t )|x i1,1:N \m


Roland Lamberti

Franc ¸ois Desbouvries Yohan Petetin

Franc ¸ois Septier

Independent Resampling Sequential Monte Carlo Algorithms

Keywords: Sequential Monte Carlo algorithms, Particle Filters, Importance Sampling, Auxiliary Particle Filter, Resampling

Sequential Monte Carlo algorithms, or Particle Filters, are Bayesian filtering algorithms which propagate in time a discrete and random approximation of the a posteriori distribution of interest. Such algorithms are based on Importance Sampling with a bootstrap resampling step which aims at struggling against weights degeneracy. However, in some situations (informative measurements, high dimensional model), the resampling step can prove inefficient.

In this paper, we revisit the fundamental resampling mechanism which leads us back to Rubin's static resampling mechanism. We propose an alternative rejuvenation scheme in which the resampled particles share the same marginal distribution as in the classical setup, but are now independent. This set of independent particles provides a new alternative to compute a moment of the target distribution and the resulting estimate is analyzed through a CLT. We next adapt our results to the dynamic case and propose a particle filtering algorithm based on independent resampling. This algorithm can be seen as a particular auxiliary particle filter algorithm with a relevant choice of the first-stage weights and instrumental distributions. Finally we validate our results via simulations which carefully take into account the computational budget.

I. INTRODUCTION

Let {X k ∈ R m } k≥0 (resp. {Y k ∈ R n } k≥0 ) be a hidden (resp. observed) process. Let X 0:k , say, denote {X i , 0 ≤ i ≤ k}, x 0:k = {x i , 0 ≤ i ≤ k}, and let p(x) (resp. p(x|y)), say, denote the probability density function (pdf) of random variable (r.v.) X (resp. of X given Y = y); capital letters are used for r.v. and lower case ones for their realizations. We assume that {(X k , Y k )} k≥0 is a Hidden Markov chain, i.e. that p(x 0:k , y 0:k ) = p(x 0 )

k i=1 f i (x i |x i-1 ) k i=0 g i (y i |x i ). (1)
Roughly speaking, pdf f k (x k |x k-1 ) describes the dynamical evolution of the Markovian hidden process {X k } k≥0 between time k-1 and time k while the likelihood g k (y k |x k ) describes the relation at time k between an observation R. Lamberti, Y. Petetin y k and the associated hidden state x k . We address the problem of computing a moment of some function f (.) w.r.t.

the filtering pdf p(x k |y 0:k ), i.e. the pdf of the hidden state given the past observations:

Θ k = f (x k )p(x k |y 0:k )dx k . (2)
As is well known, Θ k can be exactly computed only in very specific models, and one needs to resort to approximations in the general case. In this paper, we focus on a popular class of approximations called sequential Monte Carlo (SMC) algorithms or Particle Filters (PF), see e.g. [START_REF] Gordon | Novel approach to nonlinear/ non-Gaussian Bayesian state estimation[END_REF]- [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear / non-Gaussian Bayesian tracking[END_REF] PF propagate over time a set of N Monte Carlo (MC) weighted samples {w i k , x i k } N i=1 which defines a discrete approximation N i=1 w i k δ x i k of p(x k |y 0:k ) and enables to compute an estimate Θ k of Θ k :

Θ k = N i=1 w i k f (x i k ). (3)
More precisely, the computation of the set {w i k , x i k } N i=1 is based on the sequential application of the Importance Sampling (IS) mechanism [START_REF] Gewecke | Bayesian inference in econometric models using Monte Carlo integration[END_REF]. This mechanism consists in sampling particles according to an importance distribution and next weighting these samples in order to correct the discrepancy between the target and the importance distribution. However the direct sequential application of the IS mechanism in model [START_REF] Gordon | Novel approach to nonlinear/ non-Gaussian Bayesian state estimation[END_REF] fails in practice since after a few time steps most weights get close to to zero, while only a few particles have non neglictible weights.

Consequently IS alone becomes more and more inefficient since a lot of computational effort is devoted to sampling particles which will hardly contribute to the estimate Θ k in [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear / non-Gaussian Bayesian tracking[END_REF].

As is well known, a traditional rescue against weights degeneracy consists in resampling the particles (-either at each time step or depending on some criterion such as the number of efficient particles [START_REF] Kong | Sequential imputations and bayesian missing data problems[END_REF] [6] [START_REF] Liu | Metropolized independent sampling with comparisons to rejection sampling and importance sampling[END_REF] [START_REF] Cornebise | Adaptive methods for sequential importance sampling with application to state-space models[END_REF]), i.e. of redrawing each particle with a probability equal to its weight. This yields the class of Sampling Importance Resampling (SIR) algorithms [START_REF] Smith | Bayesian statistics without tears : a sampling-resampling perspective[END_REF] [1] [10] [11]. This resampling (i.e., bootstrap) mechanism has proved to be beneficial in the long run, but its instantaneous effects are mitigated; though the resampling step indeed discards particles with low weights (such particles are likely never to be resampled), particles with significant weights are resampled several times, which results in dependency among the resampled points and support shrinkage. Consequently, particle filters based on the resampling mechanism can give poor results in some Markovian models [START_REF] Gordon | Novel approach to nonlinear/ non-Gaussian Bayesian state estimation[END_REF], such as informative models where the likelihood g k (y k |x k ) is sharp. Our aim in this paper is thus to revisit this key rejuvenation scheme in order to design new PF algorithms which would keep the benefits of the resampling mechanism, while avoiding the local impoverishment of the resulting MC approximation of the filtering distribution.

To that end we begin with revisiting the SIR mechanism at one single time step k → k + 1. This leads us back to an analysis of Rubin's static SIR mechanism [12, §2] [13] [START_REF] Smith | Bayesian statistics without tears : a sampling-resampling perspective[END_REF] [14, §9.2], in which, roughly speaking, one obtains samples x j approximately drawn from a target distribution p by drawing intermediate samples {x i } N i=1 from an instrumental distribution q, and next selecting x j among {x i } N i=1 with a probability proportional to p(x i ) q(x i ) . We first observe that the samples {x j } produced by this SIR mechanism are dependent and marginally distributed from some compound pdf qN = φ(p, q, N ) which takes into account the effects of both pdfs p and q. Here the dependency is detrimental, because samples that would be i.i.d from qN would produce, whichever the number of sampled and resampled particles, a moment estimate with reduced variance; this result is further illustrated by a central limit [START_REF] Del Moral | Genealogical and interacting particle systems with applications, ser. Probability and its applications[END_REF] DRAFT theorem (CLT) which is compared to the existing CLTs for the static IS estimate (based on the pre-resampling samples {x i } N i=1 ), on the one hand, and for the SIR estimate (based on the post-resampling ones {x j } MN j=1 ), on the other hand.

We next propose a procedure to obtain i.i.d. samples from qN , which leads to the computation of two point estimates of Θ = f (x)p(x)dx. The first one is based on unweighted i.i.d. samples and is an improved version of the classical (i.e., dependent) SIR estimate; the second one is based on post-resampling-weighted i.i.d. samples and can be seen as new IS estimate, based on the compound pdf qN . Finally we adapt these results to the sequential computation of Θ k in model [START_REF] Gordon | Novel approach to nonlinear/ non-Gaussian Bayesian state estimation[END_REF]. We thus propose two new PF algorithms. One of them has an interesting interpretation in terms of Auxiliary Particle Filter (APF); more precisely, that algorithm naturally produces a relevant importance mixture distribution from which it is easy to sample. We finally illustrate our results via simulations, and carefully compare our algorithms with existing ones in terms of Root Mean Square Error (RMSE) and computational cost. The rest of this paper is organized as follows. Section II is devoted to the static case. In section III we address the sequential case, and derive new PF based on the results of section II. In section IV we perform simulations and discuss implementation issues, and we end the paper with a conclusion.

II. IS WITH RESAMPLING VIEWED AS A COMPOUND IS SCHEME

As recalled in the introduction, resampling from time to time is a standard rescue when applying IS in the sequential case. In this section we thus focus on one such time step k → k + 1. This amounts to revisiting Rubin's static SIR mechanism (see section (II-A)), which consists in resampling points {x i } MN i=1 from the weighted distribution N i=1 w i δ xi where xi i.i.d.

∼ q and the pre-resampling weights w i ∝ p (xi) q(xi) with

N i=1 w i = 1.
As is well known, when N → ∞ the resampled points {x i } MN i=1 become asymptotically i.i.d. from the target distribution p. For finite N however, these samples are dependent and drawn from some pdf qN which differs from p and can indeed be seen as a compound IS density qN = φ(p, q, N ) produced by the succession of the sampling (S), weighting (W) and resampling (R) steps. We discuss on the benefits of drawing independent samples from qN (see section II-B), and next on reweighting these independent samples with post-resampling weights w ′ i ∝ p(xi) qN (xi) (see section II-C). In all this section we assume the scalar case for simplicity. We end the section with a summary (see section II-D).

A. The dependent SIR mechanism

Let us begin with a brief review of Rubin's classical SIR sampling mechanism and of the properties of the sampled and resampled particles.

1) Properties of the sampled particles {x

i } N i=1 :
In the context of this paper we first recall the principle of IS. Let p(x) be a probability density function and assume that we want to compute

Θ = f (x)p(x)dx = E p (f (X)). (4)
In the Bayesian framework p(x) is generally only known up to a constant, i.e. p(x) ∝ p u (x) (subscript u is for unnormalized) and it is not possible to obtain samples directly drawn from p(x). A solution is to introduce November 6, 2018 DRAFT an importance distribution q(x) which satisfies q(x) > 0 when p(x) > 0 and to rewrite Θ as the ratio of two expectations w.r.t. q, Θ = f (x) pu(x) q(x) q(x)dx pu(x) q(x) q(x)dx

= E q (f (X) p(X) q(X) ) E q ( p(X) q(X) ) . (5)
Next, each expectation is approximated by a Monte Carlo method based on N i.i.d. samples (x 1 , • • • , xN ) drawn from q(.); the IS estimate of Θ is given by

Θ IS N = N i=1 w i f (x i ) = E p(f (X)) (6)
where

p(x) = N i=1 w i δ xi (x) (7)
and where w i (the i-th normalized importance weight) reads

w i = pu(x i ) q(x i ) N j=1 pu(x j ) q(x j ) = p(x i ) q(x i ) N j=1 p(x j ) q(x j ) . ( 8
)
As is well known [START_REF] Gewecke | Bayesian inference in econometric models using Monte Carlo integration[END_REF], under mild assumptions

Θ IS N a.s. → Θ, (9)
and a CLT is available too (D→ denotes the convergence in distribution):

√ N ( Θ IS N -Θ) D → N 0, E q p 2 (X) q 2 (X) (f (X) -Θ) 2 . ( 10
)
2) Properties of the resampled particles {x i } MN i=1 : From ( 9) and [START_REF] Douc | Comparison of resampling schemes for particle filtering[END_REF], p can be seen as a discrete approximation of the target density p, and one expects that for large N , (re)sampling from p would produce samples approximately drawn from p. This is the rationale of Rubin's SIR mechanism [12, §2], [START_REF] Gelfand | Sampling based approaches to calculating marginal densities[END_REF], [START_REF] Smith | Bayesian statistics without tears : a sampling-resampling perspective[END_REF], [14, §9.2]. More precisely, let us as above draw N i.i.d. samples xi from q, and next M N i.i.d samples x i from p in [START_REF] Liu | Metropolized independent sampling with comparisons to rejection sampling and importance sampling[END_REF]. It is indeed well-known (see [START_REF] Smith | Bayesian statistics without tears : a sampling-resampling perspective[END_REF] [START_REF] Rubin | Using the SIR algorithm to simulate posterior distributions[END_REF]) that when N → ∞, each r.v. x i produced by this mechanism converges in distribution to p(.), so Rubin's technique can be seen as a two-step sampling mechanism which transforms samples drawn from q into samples (approximately) drawn from p.

This convergence result can be completed by a CLT which involves the estimate of Θ based on the unweighted set {( 1MN , x i )} MN i=1 :

Θ SIR MN = 1 M N MN i=1 f (x i ). (11)
Let N → ∞, let M N be a non decreasing sequence with M N → ∞, and let lim N MN = α > 0 (possibly ∞); then under mild conditions (see e.g. [14, §9])

M N ( Θ SIR MN -Θ) D → N (0, var p (f (X)) + α -1 E q ( p 2 (X) q 2 (X) (f (X) -Θ) 2 ) ). (12)
If α → ∞ then the asymptotic variance tends to var p (f (X)), which shows that the SIR estimate asymptotically has the same behavior as a crude Monte Carlo estimate directly deduced from M N samples according to the target distribution p(.), provided the number N of intermediate samples is large compared to M N .

November 6, 2018 DRAFT However, for computational reasons, the number of samples N and M N should not be too large in practice.

Consequently we now focus on the samples produced by the SIR procedure from a non asymptotical point of view and we have the following result (the proof is given in the Appendix).

Proposition 1: Let us consider the samples {x i } MN i=1 produced by the SIR mechanism described above. Then these samples are identically distributed according to a pdf qN , with

qN (x) = N h N (x)q(x), (13)
h N (x) = p(x) q(x) p(x) q(x) + N -1 l=1 p(x l ) q(x l ) N -1 l=1 q(x l )dx l . ( 14
)
So for fixed sample size N , the SIR mechanism produces dependent samples {x i } MN i=1 distributed from qN (these samples are independent given the intermediate set {x i } N i=1 , but become dependent when this conditioning is removed). In practice, this dependency results in support shrinkage since, by construction, an intermediate sample xi can be resampled several times, and {x i } MN i=1 is a subset of {x i } N i=1 . For instance let M N = N . If we assume that w j = 1 for some j and w i = 0 for i = j, then x i = xj for all i. By contrast, if w i = 1/N for all i, then the average number of different samples {x i } N i=1 is approximately N/3 [START_REF] Petetin | Optimal SIR algorithm vs. fully adapted auxiliary particle filter: a non asymptotical analysis[END_REF]. Nevertheless the resampling step remains useful in a dynamic setup (see section III): even though locally it leads to an impoverishment of the diversity, this step is critical for recreating diversity at the next time step.

B. The independent SIR mechanism

Observe that the two factors in [START_REF] Gelfand | Sampling based approaches to calculating marginal densities[END_REF] reflect the effects of the sampling and resampling step: pdf q is used in the S step, while h N (x), which can be interpreted as the conditional expectation of a normalized importance weight when its associated particle is x, results from the (W,R) steps. So particles drawn from qN are likely to be in regions where 1) q is large (since these particles have first been sampled); and 2) which have also been resampled because their associated weight was large enough. Now our objective is to propose an alternative mechanism which, in the sequential case, will produce the same positive effect as the classical SIR mechanism (i.e. fighting against weight degeneracy by eliminating the samples with weak importance weights), while ensuring the diversity of the final support. Such a support diversity is ensured if we draw samples independently from the continuous pdf qN (.). We first study the potential benefits of this sampling mechanism (see section II-B1) and next discuss its implementation (see section II-B2).

1) Statistical properties:

Let us now assume that we have at our disposal a set of M N i.i.d. samples {x i } MN i=1 drawn from qN (.) defined in [START_REF] Gelfand | Sampling based approaches to calculating marginal densities[END_REF] [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. Before addressing the practical computation of such a set (see section II-B2), let us study its properties by considering the crude estimate of Θ based on these M N i.i.d samples:

Θ I-SIR MN = 1 M N MN i=1 f (X i ). (15)
(I in notation I-SIR stands for independent). Our aim is to compare Θ I-SIR MN to Θ SIR MN , and more generally Θ IS N , Θ SIR MN and Θ I-SIR MN . We first have the following result (the proof is given in the Appendix).

November 6, 2018 DRAFT Proposition 2: Let us consider the three estimates Θ IS N , Θ SIR MN and Θ I-SIR MN defined in ( 6), ( 11) and ( 15) respectively.

Then

E( Θ IS N ) = E( Θ SIR MN ) = E( Θ I-SIR MN ), (16)
var( Θ SIR MN ) = var( Θ I-SIR MN ) + M N -1 M N var( Θ IS N ). (17)
Equation [START_REF] Del Moral | Genealogical and interacting particle systems with applications, ser. Probability and its applications[END_REF]

N i=1 w i ( X1 , • • • , XN )f ( Xi ))
where Xi i.i.d ∼ q(.). However, we have the following CLT (the proof is given in the Appendix).

Theorem 1: Let us consider the independent SIR estimate defined in [START_REF] Petetin | Optimal SIR algorithm vs. fully adapted auxiliary particle filter: a non asymptotical analysis[END_REF]. Let assume that N → ∞, M N is a non decreasing sequence with M N → ∞ and lim

N →∞ N M N = α > 0. Then Θ I-SIR MN satisfies M N ( Θ I-SIR MN -Θ) D → N (0, var p (f (X))) . ( 18
)
Let us comment this result. First Theorem 1 enables again to compare Θ I-SIR MN to Θ SIR MN . Comparing ( 12) and ( 18) confirms [START_REF] Del Moral | Genealogical and interacting particle systems with applications, ser. Probability and its applications[END_REF], since the asymptotic variance of Θ I-SIR MN is always lower than that of Θ SIR MN . Also note that in the independent case the asymptotic variance of Θ I-SIR MN no longer depends on α > 0.

Next Theorem 1 also gives some elements for comparing Θ I-SIR MN to Θ IS N . Let for simplicity M N = N . Then the comparison of both estimates relies on that of the asymptotic variances in [START_REF] Douc | Comparison of resampling schemes for particle filtering[END_REF] and (18):

σ 2,IS ∞ (q) = E q p 2 (X) q 2 (X) (f (X) -Θ) 2 , ( 19
)
σ 2,I-SIR ∞ = var p (f (X)). (20)
For a given target pdf p(.) and function f (.), σ 2,IS ∞ (q) depends on the importance pdf q(.) and is well known [16, §2.9] [4, Theorem 3] to be minimum for q

⋆ (x) ∝ p(x)|f (x)-Θ|; for that q ⋆ , σ 2,IS ∞ (q ⋆ ) = ( |f (x)-Θ|p(x)dx) 2 ≤ var p (f (X)), so Θ IS N (q ⋆ ) outperforms Θ I-SIR N
for large values of N . On the other hand for other importance distributions σ 2,IS ∞ (q) may become larger than σ 2,I-SIR

∞

. Also note that the variances in ( 19) and ( 20) depend on function f (.); on the other hand, for large N , Θ I-SIR N has the same behavior as a crude estimate built from samples drawn from p(.) and so is adapted for a large class of functions f (.).

2) Sampling procedure: It remains to describe a procedure to obtain i.i.d. samples from qN . Algorithm 1 ensures that the final samples {x i , • • • , x MN } are drawn independently from qN .

Compared to the classical SIR procedure, the independent SIR algorithm described in Algorithm 1 relies on a sampling step of N × M N intermediate samples x and M N independent resampling steps. Consequently, for a given budget of sampling and resampling steps, the independent procedure should be compared with a classical SIR one in which we sample N × M N points and resample M N of them. In this last case, we obtain M N dependent samples drawn from qN×MN . First, using [START_REF] Rubin | Using the SIR algorithm to simulate posterior distributions[END_REF] with α = lim N →∞

N ×MN MN

= ∞, we see that both estimates Θ I-SIR MN and November 6, 2018 DRAFT Input: an importance distribution q, N and M N Result:

{x i } MN i=1 i.i.d. ∼ qN for 1 ≤ i ≤ M N do for 1 ≤ j ≤ N do S. xi,j ∼ q(.); W. w i,j ∝ p u (x i,j )/q(x i,j ), N j=1 w i,j = 1; end R. x i ∼ N j=1 w i,j δ xi,j end
Algorithm 1: The independent SIR algorithm Θ SIR MN with N × M N intermediate samples have the same asymptotic behavior. However the independent procedure can be easily parallelized because the resampling steps are by nature independent contrary to the SIR procedure where the N × M N intermediate samples are directly resampled.

C. Reweighting the independent samples?

We finally discuss the final weights which are attributed to the resampled particles. In the SIR procedure, each final sample is weighted by 1/M N . From an IS point of view, this weighting traduces the fact that the final samples become drawn from the target distribution p(.) and independent when N → ∞ [START_REF] Rubin | Using the SIR algorithm to simulate posterior distributions[END_REF]. Moreover the convergence results of Θ I-SIR MN to Θ (see e.g. [START_REF] Del Moral | Genealogical and interacting particle systems with applications, ser. Probability and its applications[END_REF] [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]) confirm that these weights are valid from an asymptotical point of view. In the independent SIR procedure, the only difference is that the final samples are independent, even from a non-asymptotical point of view. Now, if N is finite, one can wonder if weights 1/M N are optimal. In Algorithm 1, samples {X i } MN i=1 are independent and sampled from qN . Consequently, for a given N , qN can be seen as a post-resampling compound importance distribution qN = φ(p, q, N ), and a final sample x i should be weighted by a post-resampling weight proportional to p u (x i )/q N (x i ). This yields a new estimate Θ I-SIR-w of (4) (superscript w stands for weighted)

Θ I-SIR-w MN = MN i=1 pu(x i ) qN (x i ) MN j=1 pu(x j ) qN (x j ) f (x i ), (21)
which coincides with the IS estimate ( 6) with importance distribution qN (.). It is difficult to compare Θ I-SIR MN and Θ I-SIR-w MN because the expression of the weights in this last case depends on N . However, it is interesting to note that contrary to Θ I-SIR MN , M N impacts the bias of the estimate Θ I-SIR-w

MN

. For example, if we set N = 1 (so q N = q)

and M 1 is arbitrary then Θ I-SIR-w

M1

coincides with the IS estimate with M 1 i.i.d. samples drawn from q while the unweighted estimate Θ I-SIR

M1

is a crude estimate of f (x)q(x)dx and is not adapted for the estimation of Θ. More generally, using the delta method to approximate E( Θ I-SIR MN ) and E( Θ I-SIR-w

MN

) [START_REF] Liu | Monte Carlo strategies in scientific computing[END_REF] we observe that

E( Θ I-SIR MN ) = E( Θ IS N ) ≈ Θ- 1 N E q p 2 (X) q 2 (X) (f (x)-Θ) , (22)
E( Θ I-SIR-w MN ) ≈ Θ- 1 M N E qN p 2 (X) q2 N (X) (f (x) -Θ) . ( 23
)
November 6, 2018 DRAFT So for a fixed number of sampled points N , we see that in the unweighted case the bias of Θ I-SIR MN is independent of M N . By contrast, whichever N the bias of E( Θ I-SIR-w

MN

) tends to 0 as M N → ∞.

Finally, it remains to compute p u (x i )/q N (x i ) in practice. In general, qN in ( 13) is not available in close form because it relies on the integral h N (x) in ( 14). However, the N × M N intermediate samples which have been used

in Algorithm 1 can be recycled to approximate the conditional expectation h N (x). For a given x and using the intermediate samples xi,j of Algorithm 1, a crude Monte Carlo estimate of h N (x) reads

ĥN (x) = MN i=1 pu(x) q(x) pu(x) q(x) + N -1 j=1 pu(x i,j ) q(x i,j ) . ( 24
)
Importance weights pu(x) qN (x) in ( 21) can be approximated by pu(x)

N ĥN (x)q(x) . Note that the computation of these approximated weights do not require extra computational cost since p u (x i,j )/q(x i,j ) has already been computed in Algorithm 1 to obtain i.i.d. samples.

D. Summary

In summary, we now have at our disposal four estimates to compute Θ in (4) from an importance distribution q(.). Θ IS N and Θ SIR MN are deduced from the IS and Rubin's SIR mechanisms, respectively. Θ SIR MN relies on unweighted dependent samples from qN . Using unweighted independent samples from qN produces the estimate Θ I-SIR MN which outperforms Θ SIR MN and possibly Θ IS N ; it also becomes asymptotically independent of the choice of the initial importance distribution q(.) according to theorem 1. This estimate does not suffer from the support impoverishment caused by the resampling step. On the other hand it requires a larger computational cost which, however, can be exploited in order to associate to the i.i.d. samples post-resampling importance weights based on the qN (x). We thus obtain a weighted estimate Θ I-SIR-w MN which can be seen as the estimate deduced from the IS mechanism based on the compound IS distribution qN (x). We will compare these estimates via simulations and will take into account their computational cost in Section IV-A.

III. INDEPENDENT RESAMPLING BASED PF

We now adapt the results of Section II to the Bayesian filtering problem. In section III-A we briefly recall the principle of classical SIR algorithms which are based on dependent resampling. Our SIR algorithm with independent resampling and unweighted samples is proposed in section III-B. However, computing the post-resampling weights is more challenging here than in the static case because the pdf qN of the static case becomes a sum of N terms which should be computed for each final sample. So in section III-C we revisit the algorithm of section III-B in terms of APF. We first observe that the independent SIR algorithm can be seen as the first step of an APF algorithm since it implicitly draws samples from a mixture pdf. Making full use of the APF methodology enables us to weight our final samples.

A. Classical SIR algorithms (based on dependent resampling)

We now assume that we are given some hidden Markov model (1) and we briefly recall how Θ k in (2) can be computed recursively via PF. PF relies on the sequential application of the normalized IS mechanism described in November 6, 2018 DRAFT Section II-A for the target distribution p(x 0:k |y 0:k ) which is known up to a constant according to (1). Let q(x 0:k ) be an importance distribution (q(x 0:k ) can depend on y 0:k but this dependency is not written here to avoid notational burden). Starting from N weighted trajectories x i 0:k-1 sampled from q(x 0:k-1 ), we first extend each trajectory x i 0:k-1 by a particle xi k sampled from q(x k |x i 0:k-1 ) and next update the old weights w i k-1 via

w i k ∝ w i k-1 f k (x i k |x i k-1 )g k (y k |x i k ) q(x i k |x i 0:k-1 ) , N i=1 w i k = 1. (25)
Unfortunately, it is well-known that this direct sequential application of IS leads to weight degeneracy: after a few iterations only few weights w i k have a non null value [START_REF] Doucet | On sequential Monte Carlo sampling methods for Bayesian filtering[END_REF]. A traditional rescue consists in resampling, either systematically or according to some criterion such as the Effective Sample Size [START_REF] Kong | Sequential imputations and bayesian missing data problems[END_REF] [6] which is approximated by

1/ N i=1 (w i k ) 2 .
The corresponding algorithm is given in Algorithm 2 and we shall assume that the size N of the MC approximation remains constant thoughout the iterations. Finally Algorithm 2 enables to compute two estimates of Θ k :

Θ SIS N,k = N i=1 w i k f ( Xi k ), (26)
Θ SIR N,k = 1 N N i=1 f (X i k ). (27)
As is well known, the pre-resampling estimator Θ SIS N,k is preferable to the post-resampling one Θ SIR N,k and should be used in practice; but Θ SIR N,k is recalled here because it will be compared below to the independent resampling estimator (32).

Input: q(x k |x 0:k-1 ), y k , {w i k-1 , x i 0:k-1 } N i=1 Result: {w i k , x i 0:k } N i=1 for 1 ≤ i ≤ N do S. xk i ∼ q(x k |x i 0:k-1 ); W. w i k ∝ w i k-1 f k (x i k |x i k-1 )g k (y k |x i k ) q(x i k |x i 0:k-1 ) , N i=1 w i k = 1; end if Resampling then for 1 ≤ i ≤ N do R. l i ∼ Pr(L = l|{x j 0:k-1 , xj k } N j=1 ) = w l k Set x i 0:k = (x l i 0:k-1 , xl i k ), w i k = 1 N end else Set {x i k } N i=1 = {x i k } N i=1 end Algorithm 2:
The classical SIR algorithm (based on dependent resampling)

In practice, it remains to choose the conditional importance distribution q(x k |x 0:k-1 ). A popular solution consists

in choosing q(x k |x 0:k-1 ) = f k (x k |x k-1 )
, since this pdf is part of model ( 1) and is generally easy to sample from; November 6, 2018 DRAFT another one is the so-called optimal conditional importance distribution q(x k |x 0:k-1 ) = p(x k |x k-1 , y k ) which takes into account the new observation y k and for which weights w i k no longer depend on the sampled particles {x i k } N i=1 . The optimal conditional importance distribution is generally not available in closed form but some approximation techniques have been proposed, see e.g. [START_REF] Doucet | On sequential Monte Carlo sampling methods for Bayesian filtering[END_REF] [20] [START_REF] Saha | Gaussian proposal density using moment matching in SMC methods[END_REF]. The choice of the importance distribution will be not discussed in this paper and does not impact the proposed methodology. Finally, let us mention that convergence results are also available for the PF presented in Algorithm 2, see e.g. [START_REF] Crisan | A survey of convergence results on particle filtering methods for practitioners[END_REF] [23] [START_REF] Del Moral | Genealogical and interacting particle systems with applications, ser. Probability and its applications[END_REF] [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. Some of them are based on the recursive application of the CLTs recalled in Section II.

B. An alternative SIR algorithm (based on independent resampling)

Let us first adapt Proposition 1 to the sequential context. So we address the conditional distribution given

{x i 0:k-1 } N i=1 of
p i,k (x) = w i k-1 f k (x|x i k-1 )g k (y k |x), (28)
q i,k (x) = q(x|x i 0:k-1 ). ( 29
)
Then given the initial trajectories {x i 0:k-1 } N i=1 , the new samples {X i k } N i=1 are identically distributed according to a pdf qN,k which reads

qN,k (x) = N i=1 h i,k (x)q i,k (x), (30)
where h i,k (x) coincides with the conditional expectation (given

(X i k = x)) of the i-th importance weight at time k, h i,k (x) = pi,k(x) q i,k (x) p i,k (x) q i,k (x) + l =i p l,k (x l ) q l,k (x l ) l =i q l,k (x l )dx l . ( 31
)
Note that in this proposition we focus on the distribution of

X i k given {x i 0:k-1 } N i=1 . Given {x i 0:k-1 , xi k } N i=1 , the new samples {X i k } N i=1
are independent; when we remove the dependency in {x i k } N i=1 , {X i k } N i=1 become identically distributed according to qN,k but are dependent (a same particle can be resampled several times).

Since qN,k is a pdf, a procedure which would produce samples conditionally i.i.d. from qN,k would enable us to keep the advantage of the resampling step, i.e. to recreate diversity for the next time iteration while avoiding local impoverishment of the support. Except in a particular case which will be described later, sampling directly from qN,k (x) is difficult for an arbitrary conditional importance distribution q(x k |x 0:k-1 ). We thus propose a procedure similar to Algorithm 1 but adapted to the dynamical context. The SIR algorithm with independent resampling is given by Algorithm 3. Note that a difference with Algorithm 2 is that the distribution of the discrete index L i now depends on i.

We now propose a new estimate Θ I-SIR N,k of Θ which is based on the set {X i k } N i=1 produced by Algorithm 3:

Θ I-SIR N,k = 1 N f (X i k ). (32)
November 6, 2018 DRAFT

Input: q(x k |x 0:k-1 ), y k , {w i k-1 , x i 0:k-1 } N i=1 Result: {w i k , x i 0:k } N i=1 for 1 ≤ i ≤ N do for 1 ≤ j ≤ N do S. xk i,j ∼ q(x k |x i 0:k-1 ); W. w i,j k ∝ f k (x i,j k |x i k-1 )g k (y k |x i,j k ) q(x i,j k |x i 0:k-1 )
,

N j=1 w i,j k = 1; end R. l i ∼ Pr(L = l|{x j 0:n-1 , xi,j n } N j=1 ) = w i,l k Set x i 0:k = (x l i 0:k-1 , xi,l i k ), w i k = 1 N end Algorithm 3: A SIR algorithm based on independent resampling
Comparing (32) with [START_REF] Whiteley | Recent developments in auxiliary particle filtering[END_REF], remember that the samples {X i k } N i=1 share the same pdf qN,k , but that in (32) they are now independent given {x i

0:k-1 } N i=1 . Starting from a dataset {x i 0:k-1 } N i=1 , it is ensured that Θ I-SIR N,k outperforms Θ SIR k since E( Θ I-SIR N,k |{x i 0:k-1 } N i=1 ) = E( Θ SIR n |{x i 0:k-1 } N i=1 ), (33)
var( Θ SIR N,k |{x i 0:k-1 } N i=1 ) = var( Θ I-SIR n |{x i 0:k-1 } N i=1 ) + N -1 N var( Θ SIS N,k |{x i 0:k-1 } N i=1 ). (34)
Of course, computing Θ I-SIR N,k via the samples produced by Algorithm 3 requires an extra computational cost. This point will be discussed in detail in our Simulations section, but for the moment let us make two comments: first, this algorithm can be seen as an alternative resampling scheme which ensures the diversity of the resampled support without changing the conditional distribution of the final samples; if resampling needs to be performed rarely, then the independent resampling procedure may be used only when necessary. On the other hand, we will see that

Θ I-SIR N,k
can also provide an interesting alternative to Θ SIS N,k but requires an extra computational cost; so if we want to perform the independent resampling procedure at each time step we will decrease the number N of particles associated with Θ I-SIR N,k in order to reach the same computational cost associated with Θ SIS N,k . Remark 1: Note that the idea of using extra MC samples has already been proposed in the context of Island PFs [START_REF] Vergé | On parallel implementation of sequential monte carlo methods: the island particle model[END_REF]. The idea behind this class of techniques is to exploit parallel architectures, and the rationale is as follows.

Instead of considering a unique set of N particles, the method consists in dividing the population of N samples into N 1 sets of N 2 samples such as N 1 N 2 = N . It is well known that such a configuration does not improve the classical PF with N samples, but it has the advantage to split the associated computational cost when parallel architectures are available. In other words, the objective of the PFs is not to struggle against the support impoverishment.

C. Interpretation of the independent sampling scheme in terms of APF

At this point, we have seen that it was possible to obtain an estimate of Θ k based on i.i.d. samples from the conditional pdf qN,k . As in the static case, we now wonder whether the final weights 1/N used to compute Θ I-SIR N,k (see eq. ( 32)) are optimal when N is finite. To this end we would like to make use of the expression of qN,k to [START_REF] Del Moral | Genealogical and interacting particle systems with applications, ser. Probability and its applications[END_REF] DRAFT propose an alternative weighting mechanism. At first glance, the computation of a weight which would rely on (30)-(31) seems compromised because qN,k involves a sum of N terms which should be computed for each N final sample x i k . As we will see, the interpretation of the independent SIR algorithm as a particular first step of an APF algorithm will help circumvent this limitation. Let us first begin with a brief presentation of APF filters.

1) A brief presentation of APF:

In model ( 1), the filtering density at time k can be written in terms of that at time k -1,

p(x k |y 0:k ) ∝ g k (y k |x k ) f k (x k |x k-1 )p(x k-1 |y 0:k-1 )dx k-1 . ( 35
)
Plugging an MC approximation

{w i k-1 , x i k-1 } N i=1 of p(x k-1 |y 0:k-1 ) into (35) yields p(x k |y 0:k ) ∝ g k (y k |x k ) N i=1 w i k-1 f k (x k |x i k-1 ), ∝ N i=1 w i k-1 p(y k |x i k-1 )p(x k |x i k-1 , y k ), (36)
where p(y 36) leads to a particular SMC algorithm refered to as the FA-APF [START_REF] Pitt | Filtering via simulation : Auxiliary particle filter[END_REF]. However sampling directly from p(x k |y 0:k ) is not necessarily possible because p(y k |x i k-1 ) or p(x k |x k-1 , y k ) are often unavailable. To that end it has been proposed [START_REF] Pitt | Filtering via simulation : Auxiliary particle filter[END_REF] to obtain samples from an instrumental mixture pdf

k |x k-1 ) = f k (x k |x k-1 )g k (y k |x k )dx k-1 and p(x k |x k-1 , y k ) ∝ f k (x k |x k-1 )g k (y k |x k ). Sampling from p(x k |y 0:k ) in (
q(x k ) = N i=1 µ(x i 0:k-1 )τ (x k |x i 0:k-1 ) (37)
and to use IS in augmented dimension; finally APF aims at targeting the mixture pdf p(x k |y 0:k ) in (36) which, itself, targets the filtering distribution p(x k |y 0:k ). The resulting algorithm is displayed below.

Input: µ(x 0:k-1 ), τ (x k |x 0:k-1 ), y k , {w i k-1 , x i 0:k-1 } N i=1 Result: {w i k , x i 0:k } N i=1 for 1 ≤ i ≤ N do R. l i ∼ Pr(L = l|{x i 0:k-1 } N i=1 ) = µ(x l 0:k-1 ) S. x i k ∼ τ (x k |x l i 0:k-1 ); W. w i k ∝ w l i k-1 f k (x i k |x l i k-1 )g k (y k |x i k ) µ(x l i 0:k-1 )τ (x i k |x l i 0:k-1 ) , N i=1 w i k = 1; Set x i 0:k = (x l i 0:k-1 , x i k ) end Algorithm 4: The APF algorithm
Let us comment the choice of the instrumental distribution q(x k ) in (37). Compared to the SIS algorithm of paragraph III-A we see that there is an additional degree of freedom, µ(x 0:k-1 ), which is called the first stage weight; τ (x k |x i 0:k-1 ) refers to a given conditional importance distribution. Generally, the objective of the first stage weights is to avoid the computational waste induced by the resampling step of the SIR algorithm by pre-selecting trajectories at time k -1 which are in accordance with the new observation y k . Designing this pdf q(x k ) is critical and classical approximations of the predictive likelihood such as the likelihood taken at the mode of the transition November 6, 2018 DRAFT pdf (i.e. µ(

x i 0:k-1 ) ∝ w i k-1 g k (y k |φ(x i k-1 )) where φ(x i k-1 ) is the mode of f k|k-1 (x k |x i k-1
)) can actually damage the performance of the estimate. This is why it is often suggested in practice to build a first-stage weight as close as possible to w k-1 p(y k |x k-1 ), although this problem is generally difficult [START_REF] Johansen | A note on the auxiliary particle filter[END_REF] [27] due to the computation of the predictive likelihood p(y k |x k-1 ). It remains to choose the importance distribution τ (x k |x 0:k-1 ); as in the SIR algorithm, one generally tries to approximate the optimal importance distribution p(x k |x k-1 , y k ). Finally note that similarly to classical IS, the FA-APF setting is not necessarily optimal from an asymptotic point of view even if it performs very well in practice [START_REF] Douc | Optimality of the auxiliary particle filter[END_REF].

2) Independent resampling as the first step of a canonical APF algorithm: Let us now turn to the interpretation of our independent resampling procedure in terms of APF. Let us observe that qN,k in (30) can be rewritten as

qN,k (x) = N i=1 h i,k (x)q i,k (x)dx × h i,k (x)q i,k (x) h i,k (x)q i,k (x)dx (38)
and so can be seen as one particular mixture pdf q(x k ) in (37), in which the weights µ ind (x i 0:k-1 ) are given by h i,k (x)q i,k (x)dx and the components τ ind (x k |x i 0:k-1 ) by

h i,k (x)q i,k (x)
h i,k (x)q i,k (x)dx . We now verify that the couple of samples (l i , x i k ) produced by the independent resampling algorithm (Algorithm 3) can indeed be seen as an augmented sample according to qN,k (x) in (38):

• given {x j 0:k-1 } N j=1 and {x i,j k } N j=1 , Pr(L i = l) = w i,l k . Since xi,j k ∼ q i (x), the distribution of l i given {x j 0:k-1 } N j=1 becomes Pr(L i = l) = E(w l k |{x j 0:k-1 } N j=1 ) = h l,k (x)q l,k (x)dx; • given {x j 0:k-1 } N j=1 , {x i,j k } N j=1 and l i , x i k = xi,l i k .
Removing the dependency in {x i,j k } N j=1 , the distribution of x i k given {x j 0:k-1 } N j=1 and l i becomes

h l i ,k (x)q l i ,k (x)
h l i ,k (x)q l i ,k (x)dx . In summary, our independent resampling procedure is nothing but the first step of one particular APF algorithm, because the pdf qN,k (x) from which we draw i.i.d. samples (given {w i k-1 , x i 0:k-1 } N i=1 ) coincides with the mixture pdf (38), which itself constitutes a class of instrumental distributions q(x k ) in (37) parametrized by q(x k |x 0:k ).

In order to appreciate the relevance of that particular solution let us comment on the choice of the first-stage weights µ ind (x i 0:k-1 ) and distributions τ ind (x k |x i 0:k-1 ):

• at time k -1, trajectories {x i 0:k-1 } N i=1 are first resampled according to a first stage weight which coincides with the expectation of the importance weights w i k of the SIR algorithm defined in [START_REF] Pitt | Filtering via simulation : Auxiliary particle filter[END_REF]. In other words, these trajectories are preselected in such a way that the new importance weight w i k which would be affected in the weighting step of the SIR algorithm will tend to be large;

• once a trajectory x i 0:k-1 has been selected, it is not ensured that its associated weight w i k will indeed be large. By sampling according to a pdf proportional to h i,k (x)q i,k (x), the objective is to produce a sample in the region where h i,k (x) (the conditional expectation of the importance weight w i k , given that (X i k = x)) and the distribution q i,k (x) are large.

Consequently, the mixture pdf qN,k (x) appears as a natural instrumental candidate for the APF when the objective is to pre-select the trajectories and to extend them in accordance with the given conditional importance distributions q i,k (x) = q(x|x i 0:k-1 ) used in the SIR algorithm. If the SIR algorithm IS densities q i,k (x) coincide with the optimal importance distribution p(x|x i k-1 , y k ), then one can see easily that our canonical APF instrumental pdf (38) November 6, 2018 DRAFT reduces to the target mixture (36) (since h k,i in ( 31) is reduced to a term proportional to w i k-1 p(y k |x i k-1 )) and the independent SIR procedure to the FA-APF algorithm. In that case one can sample from qN,k very efficiently (since (38) is a known mixture) and the resulting estimate outperforms the SIR estimate Θ SIR N,k with optimal conditional importance distribution [START_REF] Cappé | Inference in Hidden Markov Models[END_REF] [START_REF] Petetin | Optimal SIR algorithm vs. fully adapted auxiliary particle filter: a non asymptotical analysis[END_REF]. In the case where the FA-APF algorithm is not available, it remains possible to sample from the mixture pdf qN,k (x) in (38) as soon as we can sample from the root pdf q i,k (x), even when µ ind (x i 0:k-1 ) cannot be computed, or one cannot sample from τ ind (x k |x i 0:k-1 ). 3) Reweighting the independent samples?: We can finally use this APF interpretation in order to reweight our conditional independent samples {x i k } N i=1 . Since qN,k can be seen as a mixture (37) with parameters µ ind (x i 0:k-1 ) and τ ind (x k |x i 0:k-1 ), µ ind (x i 0:k-1 ) × τ ind (x|x i 0:k-1 ) reduces to h i,k (x)q i,k (x). Finally when we target mixture (36), the second-stage weights associated with the independent samples x i k produced by Algorithm 3 read

w i k ∝ w l i k-1 f k|k-1 (x i k |x l i k-1 )g k (y k |x i k ) h l i ,k (x i k )q l i ,k (x i k ) , N i=1 w i k = 1. ( 39
)
We thus obtain a new estimate of Θ k ,

Θ I-SIR-w N,k = N i=1 w i k f (x i k ) ( 40
)
where w i k are defined in (39). The practical computation of these final weights relies on that of h i,k (x) in (31), which can be approximated by recycling the extra samples xi,j generated in Algorithm 3,

ĥl (x) = N i=1 p l,k (x) q l,k (x) p l,k (x)
q l,k (x) + j =l p j,k (x i,j ) q j,k (x i,j ) .

(41)

D. Summary

Let us summarize the discussions of section III. When the objective is to compute Θ k in (2) we have several options:

1) using the classical SIR algorithm (see Algorithm 2) in which we compute Θ SIS N,k defined in [START_REF] Johansen | A note on the auxiliary particle filter[END_REF]. The resampling step which follows the computation of this estimate produces a conditionally dependent unweighted set of particles sampled from qN,k ;

2) an alternative to avoid the local impoverishment induced by the traditional resampling step is to perform Algorithm 3 and to compute estimate Θ I-SIR N,k . This estimate is still based on an unweighted set of particles marginally sampled from qN,k but these samples have become conditionally independent;

3) finally, the samples produced by Algorithm 3 can also be seen as the result of a sampling procedure according to a partial APF instrumental mixture pdf (38). Using further the APF methodology with mixture qN,k it is possible to target mixture (36) which itself is an approximation of p(x k |y 0:k ). This leads to estimate Θ I-SIR-w N,k in (40), in which the weights (39) are estimated by recycling the extra samples produced by Algorithm 3.

These three estimates are now going to be compared (in terms of performances and computational cost) in the next section.

November 6, 2018 DRAFT

IV. SIMULATIONS

We now validate our discussions through computer-generated experiments. In section IV-A we first illustrate the results of Section II and we compare the classical resampling mechanism to the independent one with both unweighted and weighted samples. We also discuss the computational cost associated with our independent resampling mechanism.

In section IV-B we next perform simulations in the ARCH model. On the one hand, the FA-APF algorithm can be computed in this model [START_REF] Pitt | Filtering via simulation : Auxiliary particle filter[END_REF]. On the other hand, remember that our weighted estimate (40) can be interpreted as the estimate deduced from a particular APF which uses the instrumental mixture pdf qN,k in (38), from which it is always possible to sample from (with an extra computational cost). Thus the estimate deduced from the FA-APF algorithm is used as a benchmark and enables us to analyze the relevance of the instrumental pdf qN,k in the APF algorithm.

Next in section IV-C we compute our independent estimates for a target tracking problem with range-bearing measurements. Our estimates are compared to those obtained from the classical SIR algorithm, for a given computational budget measured via the number of sampling operations; this means that we compare

Θ I-SIR M,k and Θ I-SIR-w M,k
(M is the number of particles after the independent resampling step) to

Θ SIS N,k in which N = M 2 +M 2
. Thus all estimates are based on M 2 + M sampling operations (we do not distinguish if we sample according to a continuous or a discrete distribution). The relative performances of the estimates are analyzed in function of the parameters of the state-space model.

Finally in section IV-D we compute our estimates in models where the dimension m of the hidden state is large and we analyze their performances w.r.t. classical PF estimates in function of the dimension m and with a fixed number of sampling operations. Finally throughout this section our simulations are averaged over P = 1000 MC runs, we set f (x) = x in (2) and we use an averaged Root Mean Square Error (RMSE) criterion, defined as

RMSE( Θ) = 1 T T k=1 1 P P p=1 || Θ k,p -x k,p || 2 1/2 (42)
where x k,p is the true state at time k for the p-th realization, Θ k,p is an estimate of x k,p and T is the time length of the scenario.

A. Comparison of static sampling procedures

Let us first consider the (static) Bayesian estimation problem in which we look for computing

Θ = E(X|y) = xp(x|y)dx (43)
via the techniques described in Section II. We assume that p(x|y) is known up to a constant, p(x|y) ∝ p(x)p(y|x)

where p(x) = N (x; 0; σ 2 x ) and p(y|x) = N (y; x, σ 2 y ) with σ 2 x = 10 and σ 2 y = 3. We chose the IS distribution q(x) = p(x). For a given number of final samples N , we compute six estimates: the estimate Θ SIS N deduced from the IS mechanism with importance distribution q(.); the estimate Θ SIR which relies on extra samples to approximate the weight proportional to p(x, y)/q N (x).

In Fig. 1 we display the distance of each estimate w.r.t. the true expectation E(X|Y = y) in function of the number of samples N . As expected, the estimate Θ I-SIR N based on N independent samples drawn from qN outperforms the estimate Θ SIR N which is computed from N dependent samples drawn from qN . However, an interesting result is that Θ I-SIR N also outperforms Θ SIS N . It means that the distribution qN produced by the SIR mechanism is more adapted than the prior q(x) = p(x), which is not surprising since qN uses implicitly the observation y through the resampling mechanism of intermediate samples. Of course, the computation of Θ I-SIR N requires an extra computational cost but it is interesting to note that the size of the final support is the same in the three cases. We finally compare the estimates based on the same computational cost. When N increases, these estimates have the same asymptotical behavior. It can be seen that the estimate Θ SIR-2 N based on N samples drawn from qN 2 outperforms Θ I-SIR N . However, when our i.i.d. samples are weighted by a term proportional to p(x, y)/q N (x) in an IS perspective, our estimate Θ I-SIR-w N has the best performance whatever N . We finally note that contrary to the independent procedure, weighting the samples when they are dependent does not improve the performance when compared to the estimate based on dependent and unweighted samples; indeed, Θ SIR-w N is not any better than Θ SIR N . The performances of these algorithms are also presented in terms of RMSE (w.r.t. to the true value of X) in Table I.

B. Comparison with APF algorithms

We now focus on the interpretation of our independent resampling algorithm in terms of APF. We study the ARCH model which is a particular hidden Markov model [START_REF] Gordon | Novel approach to nonlinear/ non-Gaussian Bayesian state estimation[END_REF]

in which f k (x k |x k-1 ) = N (x k ; 0; β 0 + β 1 x 2 k-1
) and g k (y k |x k ) = N (y k ; x k ; R). We set R = 1, β 0 = 3 and β 1 = 0.75. In this model one can compute p(y

k |x k-1 ) = N (y k ; 0; R+β 0 +β 1 x 2 k-1 ) and p(x k |x k-1 , y k ) = N (x k ; β0+β1x 2 k-1 R+β0+β1x 2 k-1 y k ; R(β0+β1x 2 k-1 ) R+β0+β1x 2 k-1
); consequently, it is possible to obtain i.i.d. samples from the target mixture (36) and thus to compute the estimate Θ FA N,k based on the FA-APF algorithm. Remember that the FA-APF can also be seen as a particular case of our independent resampling Algorithm 3 in which the importance distribution q(x k |x 0:k-1 ) coincides with p(x k |x k-1 , y k ) (see section III-C2). However this setting can be implemented in specific models only, while Algorithm 3 can be used with any importance distribution q(x k |x 0:k-1 ), while keeping the same interpretation as the FA-APF (see our discussion in section III-C2). So we also compute our estimates Θ I-SIR N,k and Θ I-SIR-w N,k which can be seen as an estimate deduced from the APF in which the importance mixture (37) coincides with qN,k . We finally compute the estimate Θ APF N,k which is deduced from the APF with µ(x 0:k-1 ) ∝ w k-1 p(y k |x k-1 ) and τ (x k |x 0:k-1 ) = f k (x k |x k-1 ); with this configuration, the particles are pre-selected with the so-called optimal first stage weight and sampled from the transition pdf.

The RMSE of each estimate is displayed in Fig. 2(a) as a function of the number of samples N . Interestingly enough, our weighted independent resampling algorithm which produces Θ I-SIR-w N,k has the same performances as the FA-APF algorithm when N ≥ 15, without using the predictive likelihood p(y k |x k-1 ) nor the optimal importance distribution p(x k |x k-1 , y k ). It means that the mixture pdf qN which has been interpreted in section III-C2 is indeed as relevant as the target mixture (36); so in general models where the FA-APF is no longer computable, one can expect that our estimate Θ I-SIR-w N,k would give a performance close to that deduced from FA-APF. Indeed, one advantage of the mixture pdf qN deduced from the resampling mechanism is that its interpretation does not depend on the importance distribution q i,k which has been chosen and that it is possible to sample from it in general hidden Markov models [START_REF] Gordon | Novel approach to nonlinear/ non-Gaussian Bayesian state estimation[END_REF]. We also observe that re-weighting the final samples is beneficial w.r.t. attributing uniform

C. Tracking from range-bearing measurements

We now study the performance of our algorithms in a tracking scenario with range-bearing measurements. We look for estimating the state vector X k = [p x,k , ṗx,k , p y,k , ṗy,k ] T (position and velocity in Cartesian coordinates) of a target from noisy range-bearing measurements y k . The pdfs in model (1) associated with this tracking problem are November 6, 2018 DRAFT

f k (x k |x k-1 ) = N (x k ; Fx k-1 ; Q) and g k (y k |x k ) = N (y k ;   p 2 x,k + p 2 y,k arctan p y,k p x,k   ; R) where τ = 1, R =   σ 2 ρ 0 0 σ 2 θ   , F =         1 τ 0 0 0 1 0 0 0 0 1 τ 0 0 0 1         , Q = σ 2 Q         τ 3 3 τ 2 2 0 0 τ 2 2 τ 0 0 0 0 τ 3 3 τ 2 2 0 0 τ 2 2 τ         .
The conditional importance distribution used to sample particles is the transition pdf q(x k |x 0:k-1 ) = f k (x k |x k-1 ); so the importance weights w i k at time k are proportional to

w i k-1 g(y k |x i k ). We compute Θ SIS N,k (see (26)), Θ I-SIR M,k (see (32)), Θ I-SIR-w M,k (see (40)) with N = M 2 +M 2
to set the number of sampling operations. We also compare these estimates with Θ IPF N,k deduced from the Island PF with 5 islands and N/5 particles per island. The results are displayed for two set of parameters. Fig. 3 gives the best performances. Contrary to Θ SIS N,k and Θ IPF N,k our estimate does not suffer from the degeneration of the importance weights. Indeed when the measurements are informative (and so the likelihood is sharp), few importance weights have a non null value. However, the independent resampling procedure ensures the diversity of the final samples when we use uniform weights. Concerning Θ I-SIR-w M,k , remember that it relies on the MC approximation (41). A close analysis of (41) when the likelihood is sharp shows that the final weights tend to be null except that of the particle with the larger likelihood; consequently, in this case the estimate Θ I-SIR-w M,k is affected by the lack of diversity.

D. High dimensional problems

We finally study the impact of the dimension of the hidden state X k . We consider a state vector of dimension on weighted samples from q but using unweighted samples when the dimension is large ensures the diversity and gives better performances. Note that the dependent and independent SIR algorithms give approximately the same performance when m is low but the gap between the dependent and the independent SIR estimates increases with the dimension. V. CONCLUSION SMC algorithms in Hidden Markov models are based on the sequential application of the IS principle. However the direct sequential application of the IS principle leads to the degeneration of the weights, against which multinomial resampling has been proposed. This rejunevation scheme, which is now routinely used in SIR algorithms, enables to discard particles (or trajectories) with low weights, but particles with large weights will be resampled several times, which leads to dependency and support degeneracy. In this paper we thus revisited the resampling step used in the classical SIR algorithms. We first addressed the static case, showed that the particles sampled by Rubin's SIR mechanism are dependent samples drawn from some pdf qN , and proposed an alternative sampling mechanism which produces independent particles drawn from that same marginal pdf qN . This set of independent samples enables us to build a moment estimator which outperforms the classical SIR-based one, both from a non-asymptotical and an asymptotical points of view. Finally the succession of the sampling, weighting and resampling steps indeed transforms an elementary instrumental pdf q into a compound importance distribution qN = φ(p, q, N ), which leads us to reweight the (originally unweighted) resampled particles x i by post-resampling weights proportional to p(x i ) q(x i ) . Such post-resampling weights cannot be computed exactly, but can easily be estimated by recycling the extra MC samples which were needed for producing the independently resampled particles.

m = 4 × l, x k = [p 1 x,k , ṗ1 x,k , p 1 y,k , ṗ1 y,k , • • • , p l x,k , ṗl x,k , p l y,k , ṗl y,k ] T . Each component x l k = [p l x,k , ṗl x,k , p l y,k , ṗl y,k ] T evolves independently from all the other components, according to f k|k-1 (x k |x k-1 ) = N (x k ; Fx k-1 ; Q) where τ = 1, F =         1 τ 0 0 0 1 0 0 0 0 1 τ 0 0 0 1         , Q = σ 2 Q         τ 3 3 τ 2 2 0 0 τ 2 2 τ 0 0 0 0 τ 3 3 τ 2 2 0 0 τ 2 2 τ         . Each component is observed independently via g k (y k |x k ) = N (y k ; Hx k ; R) where H =   1 0 0 0 0 0 1 0   , R =   σ 2 x 0 0 σ 2 y   .
November 6, 2018 DRAFT

We next adapted this methodology to the dynamic case, in order to estimate a moment of interest in an hidden Markov model. The computation of the post-resampling weights is more challenging than in the static case, but reinterpreting our independent resampling scheme as the first step of a particular APF algorithm enables us to make full use of the APF methodology and so to reweight the final samples via the second-stage APF weights. Finally we validated our discussions by computer-generated experiments and carefully took into account the computational budget. Simulations in model where the FA-APF algorithm is computable show that the independent resampling gives a performance close to the FA-APF algorithm. Consequently, it confirms the relevance of the instrumental mixture pdf used implicitly by the independent resampling PF which can be used in any hidden Markov model since it not require to compute the predictive likelihood nor the optimal importance distribution. Finally independent PF gives very satisfying results when applied in highly informative models which are challenging for classical PF and limits the degeneration phenomenon in high dimensional models.

APPENDIX PROOF OF PROPOSITION 1

Let A be any Borel set. Let 1 A (x) = 1 if x ∈ A and 0 otherwise. Then for any l, 1 ≤ l ≤ M N , Pr(X l ∈ A) q(x j )dx 1:i-1,i+1:N ]q(x i )dx i = N i=1 A h N (x i )q(x i )dx i = A N h N (x)q(x)dx, so X l has pdf qN w.r.t. Lebesque measure.

= IR N [ N i=1 w i (x 1 , • • • , xN )1 A (x i )]

PROOF OF PROPOSITION 2

Let X i (for any i, 1 ≤ i ≤ M N be produced by the classical SIR mechanism. Then

N 2 N

 2 deduced from the SIR mechanism with N intermediate samples and M N = N final samples; our estimate Θ I-SIR N based on N unweighted independent November 6, 2018 DRAFT samples drawn from qN (see (15)); our estimate Θ I-SIR-w N based on N weighted independent samples from qN (see (21)). Remember that the computation of the independent resampling mechanism is based on the sampling of N 2 intermediate particles and N resampling steps and thus requires an extra computational cost w.r.t. the dependent one. Consequently, we also compute Θ SIR-based on the classical SIR procedure with N 2 intermediate samples and N (dependent) resampling steps; in other words this estimate relies on N dependent samples obtained from qN 2 . Finally, we would like to observe the effects of weighting the final samples in the dependent resampling case; so we compute Θ SIR-w N

N)Fig. 1 .

 1 Fig. 1. Static linear and Gaussian model -σ 2 x = 10, σ 2 y = 3 -Bayesian estimates of E(X|y) based on the independent resampling mechanism outperform the estimates based on the traditional IS and SIR mechanisms altough they require an extra computational cost. When the computational cost is fixed, the estimate based on weighted i.i.d. samples from qN outperforms the estimates based on identically distributed samples from qN 2 .

  Fig. 2. ARCH model -R = 1, β 0 = 3 and β 1 = 0.75 -(a) The estimate based on the independent resampling mechanism with a final reweighting has the same performances as the estimate deduced from the FA-APF. The final reweighting mechanism is beneficial when compared to the use of uniform weighs -(b) When N is large, final weights associated to the estimate Θ I-SIR-w N,k

  (a) corresponds to the case where σ Q = √ 10, σ ρ = 0.25 and σ θ = π 720 while Fig 3(b) corresponds to a very informative case where σ Q = √ 10, σ ρ = 0.05 and σ θ = π 3600 . For the first configuration, we observe that Θ I-SIR-w M,k outperforms the other estimates and improves Θ I-SIR M,k which does not rely on weighted samples. Compared to the classical SIS estimate, Θ I-SIR M,k gives better performance as long as the number of samples M is weak (M < 30, so N < 465) but is next outperformed when the number of samples is large. As shown in Fig. 3(b), when the observations become informative, Θ I-SIR M,k

Fig. 3 .

 3 Fig. 3. Target tracking model from range-bearing measurements -(a) the independent resampling procedure with final weighting outperforms the other estimates and is particularly interesting when the number of final samples is weak -(b) in the informative case, all estimates suffer from the degeneration of the importance weights except that based on the unweighted independent resampling algorithm. To achieve the same performances as Θ I-SIR-w M,k

Fig. 4 . 2 Q = 25 , σ 2 x = 4 and σ 2 y = 4 .

 4225224 Fig. 4. Multi-dimensional linear Gaussian model -σ 2 Q = 25, σ 2 x = 4 and σ 2 y = 4. The estimates of interest are compared as a function of the dimension m of the hidden state x k for a fixed number of sampling operations. The independent resampling mechanism limits the impact of the large dimension m and estimate Θ I-SIR 1000,k when m = 46 has the same performance as Θ I-SIR 500500,k when m = 32.

w i (x 1 ,IR N - 1 w i (x 1 ,

 111 • • • , xN )1 A (x i ) • • • , xN )

E

  (f (X i ))|x 1:N ) = Θ IS N . (44) So E( Θ SIR MN )|x 1:N ) = Θ IS N , and E( Θ SIR MN ) = E( Θ IS N ).On the other hand E( Θ I-SIR MN ) = E( Θ SIR MN ), whence[START_REF] Hesterberg | Advances in importance sampling[END_REF]. Nextvar( Θ SIR MN ) k =l Cov(f (X k ), f (X l )).

  the resampled particles X i k and we have the following result (the proof is omitted). Proposition 3: Let us consider the samples {X i k } N i=1 produced by the SIR mechanism of Algorithm 2. Let

in which X i ∼ qN for all i. The first term is equal to var( Θ I-SIR MN ). Let us compute the second term. For all k, l,

). Using (44) again, we conclude that Cov(f (X k ), f (X l )) = var( Θ IS N ), whence [START_REF] Del Moral | Genealogical and interacting particle systems with applications, ser. Probability and its applications[END_REF].

PROOF OF THEOREM 1

We first introduce the following notations:

and we will assume that

Our objective is to show that A N converges to a centered Gaussian distribution with variance var p (f (X)) and that B N converges to 0.

Convergence of B N

We have recalled (see [START_REF] Smith | Bayesian statistics without tears : a sampling-resampling perspective[END_REF]) that under mild assumptions

According to Theorem 9.1.10 in [START_REF] Cappé | Inference in Hidden Markov Models[END_REF],

) is bounded and so its upper bound is finite. According to the corollary of Theorem 25.12 in [START_REF] Billingsley | Probability and Measure, ser. Wiley Series in Probability and Statistics[END_REF], it is ensured that

To prove the convergence when N → ∞, we need a CLT for triangular arrays and we use the version presented in Theorem 9.5.13 of [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. The required assumptions are:

November 6, 2018 DRAFT

3) for any positive C,

Using again Theorem 9.1.10 of [START_REF] Cappé | Inference in Hidden Markov Models[END_REF] and Theorem 25.12 of [START_REF] Billingsley | Probability and Measure, ser. Wiley Series in Probability and Statistics[END_REF],

) which converges to Θ(f 2 1 |f |≥C ) and assumption 3) is satisfied. Consequently,

Combining ( 52), ( 54) and (49) we obtain [START_REF] Liu | Monte Carlo strategies in scientific computing[END_REF].