Both the particle and Kalman ﬁlters attempt to approximate the minimum mean-square error (MMSE) estimate of the time-varying parameter. In this scenario, a prior model of the time evolution of the parameter of interest is assumed before the MMSE estimation takes place. The Kalman ﬁlter is the (optimal) MMSE estimator for a linear dynamical system with Gaussian noise. For a nonlinear system with nonGaussian noise, the particle ﬁlter approximates the mean of posterior distribution1at each discrete time step with a ﬁnite number of samples or particles. For these general nonlinear systems, the particle ﬁlter approaches the MMSE estimator as the number of particles approaches inﬁnity. 1 Beam Tracking in Wireless Communications One of the deﬁning characteristics of 5G and the future generation of wireless communications is the use of high frequencies (Millimeter wave, Terahertz wave, free-space optics) in the electromagnetic spectrum. The use of microwave frequencies served us well for the 2G, 3G and 4G systems. However, the microwave spectrum was not deemed suﬃcient to cater to high data-rate demands of the present and future. This drove the push for the deployment of higher frequency bands in the current and future generation of wireless communications. The high frequency signals will suﬀer a greater loss in energy per unit propagation distance in the atmospheric channel. Hence, there is a need to focus much of the energy of the high frequency signals in the direction of the receiver so that the receiver can receive enough energy to decode the transmitted symbol. These “directional” signals require the use of pointing and tracking subsystems in the receiver architecture so that the narrow, pencil-thin beams can be tracked eﬃciently in order to maintain the received signal-to-noise ratio above a certain threshold. There is a substantial number of studies in literature that have addressed beam tracking for millimeter wave systems [1–14]. Additionally, a number of articles deal exhaustively with ∗Correspondence: The author may be reached at muhammad.bashir@fulbrightmail.org orsalman.bashir@ alumni.purdue.edu for any comments regarding this article. 1The mean of posterior distribution is also know as the minimum mean-square error estimate. beam tracking in free-space optical communications as well [15–26]. As the reader may ﬁnd out, a signiﬁcant number of these beam tracking papers deal with Bayesian ﬁltering techniques such as diﬀerent versions of Kalman and particle ﬁlters. Even though there is no shortage of tutorials on particle ﬁlters, we believe that the majority of those tutorials do not provide a clear and easy introduction to the theory of particle ﬁlters. In this paper, we have tried to address this challenge by starting from the basics of Monte Carlo integration, and working our way up to the particle ﬁltering algorithm. We mention here that this study is not a detailed exposition on any specialized version of a particle ﬁlter, but rather an attempt to furnish an easy-to-follow introduction to the general framework of a particle ﬁltering algorithm. 1.1. Organization of This Study This article on an introduction to particle ﬁlters is organized as follows. In Section 2, we describe the state space model for a nonlinear dynamical system. In Section 3, we present the concepts behind probabilistic or Monte Carlo integration. The Monte Carlo integration is needed to approximate the posterior mean with the help of a ﬁnite number of randomly chosen samples at each time step of the ﬁltering process. Section 4 introduces the idea of importance sampling in which we sample from a proposal distribution instead of the target distribution to estimate the posterior mean. Since the sampling and storage complexity grows with time for an importance sampling ﬁlter, we introduce the idea of sequential importance sampling (SIS) in Section 5 that maintains a constant sampling complexity of O(N) at each stage of time n, whereNis the number of particles. A common version of the SIS ﬁlter adopts the prior density as a proposal for particle sampling. For this choice of the proposal density, the weight at each stage of the SIS ﬁlter is proportional to the likelihood function. However, in case the support of the prior and the likelihood densities do not overlap signiﬁcantly, most of the generated particles will possess insigniﬁcant weights. However, this problem can be overcome by resampling. In case the system suﬀers from sample impoverishment or loss of diversity of particles after resampling, the resampled particles can be perturbed about their original positions to generate diversity. All these ideas are discussed in Section 5. 2 State Space Model For the sake of brevity, lets denote observations from time 0 to time nasz0:n:=/bracketleftBig z0···zn/bracketrightBigT and, and the (random) parameter process X0:n:=/bracketleftBig X0···Xn/bracketrightBigT. A realization of this 2 random process is denoted by x0:n:=/bracketleftBig x0···xn/bracketrightBigT. The evolution of the state sequence in general is xn=fn(xn−1,un−1), n≥1 (1) where fn:Rnx×Rnu→Rnxis a possibly nonlinear function of state parameter xn−1andun−1. The quantity un−1is i.i.d noise process (possibly not gaussian) statistics. The integers nxand nurepresent dimensions of signal and noise processes, respectively. The measurement model of the dynamical system is given by zn=hn(xn,vn), n≥0, (2) where hn:Rnx×Rnv→Rnxis a possibly nonlinear function, vnis the measurement noise (possibly not gaussian) sequence and nvis its dimension. The prediction (orprior ) stage at time kis obtained via Chapman-Kolmogorov equation as p(xn|z0:n−1) =/integraldisplay∞ −∞p(xn|xn−1)p(xn−1|z0:n−1)dxn−1. (3) The observation znat timenis used to update the prior stage via Bayes’ rule p(xn|z0:n) =p(zn|xn)p(xn|z0:n−1) p(zn|z0:n−1), (4) where the normalizing constant p(zn|z0:n−1) =/integraldisplay∞ −∞p(zn|xn)p(xn|z0:n−1)dxn (5) depends on likelihood function p(zn|xn). Likelihood function is deﬁned by measurement model given by (2). From the point-of-view of ﬁltering , we are interested in computing posterior means (or expectations) E[xn|z0:n] =/integraldisplay∞ −∞xnp(xn|z0:n)dxn (6) at every instant of discrete where p(xn|z0:n) is termed the posterior density in (6). The posterior mean is also known as the minimum mean-square error (MMSE) estimate of the parameter in question. In case, the end goal is smoothing of data as opposed to ﬁltering, the posterior is deﬁned asp(x0:n|z0:n). The integral in (6) will become computationally intractable to compute since the dimen- sion of the integral is growing linearly with time-index n. In this case, one can use sub-optimal 3 Figure 1 . Estimation of area of circle with radius rcontained in a square using Monte Carlo technique. sample based methods—such as Monte Carlo techniques—in order to estimate posterior expec- tations using independent samples from the posterior density functions. For these systems, a sequential version of the sample based ﬁlter leads to a constant computational complexity at any stagenof the ﬁltering process. In the next section, we consider Monte Carlo integration techniques. 3 Monte Carlo approach 3.1. Monte Carlo methods As discussed earlier, computation of the posterior expectation can be a daunting task because of the large dimensionality of the parameter space. However, we may solve probabilistic expectations approximately by replacing the integral operation with a ﬁnite sum. In this regard, one has to obtain a ﬁnite number of independent samples from the posterior distribution and then compute a weighted sum of those samples where the weight of each sample is the value of the posterior density at that sample. This method relies on the weak law of large numbers which states that as the number of samples grow large, the “sample based” Monte Carlo expectation converges to the “true” expectation computed using integrals. We next illustrate the basic idea of Monte Carlo experiment by a very basic example. 3.2. A simple Monte Carlo technique We want to estimate the area of circle of a given radius rcontained inside a square of area Asas shown in Fig. 1. If we were to use Monte Carlo approach to compute the area of the circle, we would ﬁrst obtain NTindependent samples uniformly from the square. Thereafter, 4 we count the number of samples whose Euclidean distance from the center of the square is less than or equal to r( lets call that number Nr). We then have that Ac,Nr NT×As (7) is the area of the circle computed using the Monte Carlo technique. Then, by weak law of large numbers, we have Ac→πr2asNT→∞. (8) 3.3. Monte Carlo Integration Versus Deterministic Integration Consider E[X] =/integraltext∞ −∞xp(x)dx. Then, from the Bayesian perspective, assume p(x) to be the posterior density. The Bayesian approach integrates over high dimensional probability distributions to make inferences. Deterministic methods on the other hand involve numerical integration if integral is analytically intractable. The complexity of numerical integration increases exponentially with dimensionality of parameter space. On the other hand Monte Carlo technique positions its samples over regions of high probability density in order to approximate averages. Therefore computation resource is not wasted computing integrals at points that do not matter as much. Hence, the Monte Carlo approach can be more eﬃcient in terms of computational complexity than numerical integration. 3.4. Monte Carlo Integrals: Expectations As discussed earlier, Riemann integrals can be computed by Monte Carlo methods. This means that we can compute probabilistic expectations also with the help of Monte Carlo methods. SupposeXis a random variable with density p(x).Then the expectation of any function fofXcan be computed as I(f):= E[f(X)] =/integraldisplay∞ −∞f(x)p(x)dx. (9) However if we can obtain Ni.i.dXi’s fromp(x), then by strong law of large numbers, IN(f):=1 NN/summationdisplay i=1f(Xi) =1 NN/summationdisplay i=1/integraldisplay∞ −∞f(x)δ(x−Xi)dxa.s.−−→I (f) asN→∞ whereI(f) =/integraltext∞ −∞f(x)p(x)dx.In other words, if we deﬁne pN(x):=1 NN/summationdisplay i=1δ(x−Xi) (10) 5 thenpN(x)a.e.−−→p(x) asNgets large. It can be shown that IN(f) is unbiased estimate of I(f): E[IN(f)] =1 NN/summationdisplay i=1E[f(Xi)] =1 NN/summationdisplay i=1E[f(X)] = E[f(X)]. (11) Ifσ2 frepresents variance of f(X), then Var[IN(f)] =σ2 f N, (12) which goes to zero as N→∞ . 3.5. Monte Carlo Integration Example Lets suppose we want to compute the integral I=/integraldisplay Cg(x)dx. (13) Monte Carlo approach is to factorize the integrand as g(x) =f(x)p(x) wherep(x)≥0,and,/integraldisplay∞ −∞p(x)dx= 1. Then I=/integraldisplay Cg(x) =/integraldisplay Cf(x)p(x)dx= E[f(X)]≈1 NN/summationdisplay i=1f(Xi), (14) whereXiare independently drawn from p(x).Thus, we have (approximately) computed the integralIwith the help of a ﬁnite sum. 4 Importance Sampling It is not always possible to get independent samples from a given p(x) to computeIin (14). One possible solution is to use Importance Sampling . In importance sampling, instead of obtaining samples from p(x), we obtain samples from animportance sampling distribution ,q(x). In this case, the expectation is calculated as Ep[f(X)] =/integraldisplay∞ −∞f(x)p(x)dx=/integraldisplay∞ −∞f(x)p(x) q(x)q(x)dx= Eq/bracketleftBigg f(X)p(X) q(X)/bracketrightBigg , (15) 6 where Efdenotes the expectation with respect to a certain probability density function f(x). The density function p(x) is known as the target distribution. In literature, q(x) is also known as the proposal distribution. The following conditions are necessary for importance sampling to hold: •/integraltext∞ −∞q(x)dxshould be 1. •p(x)>0⇒q(x)>0∀x∈R. Then the expectation can be written as Ep[f(X)] = Eq/bracketleftBigg f(X)p(X) q(X)/bracketrightBigg ≈1 NN/summationdisplay i=1f(Xi)p(Xi) q(Xi). (16) We then have that, E[f(X0:n)|z0:n] =/integraldisplay∞ −∞f(x0:n)/bracketleftBiggp(x0:n|z0:n) q(x0:n|z0:n)/bracketrightBigg q(x0:n|z0:n)dx0:n. (17) Let ˜Wn:=p(x0:n|z0:n) q(x0:n|z0:n)=p(z0:n|x0:n)p(x0:n) p(z0:n)q(x0:n|z0:n). (18) However we need to know p(z0:n) =/integraltext∞ −∞p(z0:n|x0:n)p(x0:n)dx0:n, which usually is not available. Alternatively, we can set Wn:=p(z0:n|x0:n)p(x0:n) q(x0:n|z0:n). (19) Then E[f(X0:n)|z0:n] =1 p(z0:n)/integraldisplay∞ −∞f(x0:n)/bracketleftBiggp(z0:n|x0:n)p(x0:n) q(x0:n|z0:n)/bracketrightBigg q(x0:n|z0:n)dx0:n (20) =1 p(z0:n)/integraldisplay∞ −∞Wnf(x0:n)q(x0:n|z0:n)dx0:n=Eq[Wnf(X0:n)|z0:n] p(z0:n). SinceWnq(x0:n|z0:n) =p(z0:n|x0:n)p(x0:n) by (19), we have that E[f(X0:n)|z0:n] =Eq[Wnf(X0:n)|z0:n] p(z0:n)=Eq[Wnf(X0:n)|z0:n]/integraltext∞ −∞Wnq(x0:n|z0:n)dx0:n=Eq[Wnf(X0:n)|z0:n] Eq[Wn|z0:n]. (21) 7 Since we can sample from q(x0:n|z0:n), we have that q(x0:n|z0:n)≈1 N/summationtextN i=1δ(x0:n−x(i) 0:n), where x(i) 0:nare independent samples that are obtained from q(x0:n|z0:n). Substituting this estimate into (20) and using the sifting property of dirac delta functions along with the fact thatp(z0:n) =/integraltext∞ −∞Wnq(x0:n|z0:n)dxn, we obtain E[f(X0:n)|z0:n]≈1 1 N/summationtextN i=1W(i) n1 NN/summationdisplay i=1f(x(i) 0:n)p(z0:n|x(i) 0:n)p(x(i) 0:n) q(x(i) 0:n|z0:n)/bracehtipupleft /bracehtipdownright/bracehtipdownleft /bracehtipupright W(i) n=1 /summationtextN i=1W(i) nN/summationdisplay i=1f(x(i) 0:n)W(i) n. (22) Let us now deﬁne W(i) nasW(i) n:=W(i) n/summationtextN i=1W(i) n.Then, the ﬁnal estimate is E[f(X0:n)|z0:n]≈N/summationdisplay i=1W(i) nf(x(i) 0:n). (23) In this case, the posterior is approximated as p(x0:n|z0:n)≈N/summationdisplay i=1W(i) nδ(x0:n−x(i) 0:n). (24) We note from (22) that we need to sample Npoints from the n+ 1 dimensional space in order to compute the posterior expectation at each stage of time n2. Therefore, the dimensionality of the Monte Carlo expectation increases linearly with the discrete-time index n. In order to avoid this linear increase in computational complexity, we resort to a more eﬃcient version of importance sampling which is known as sequential importance sampling or recursive importance sampling . 5 Sequential Importance Sampling We have to satisfy three assumptions that will enable the sequential importance sampling algorithm to work. 5.1. Assumption on Importance Sampling Distribution The importance distribution q(x0:n|z0:n) can be modiﬁed to enable a “sequential” estimation of posterior p(x0:n|z0:n) in order to save computational resources. The sequential form can be realized ifq(x0:n|z0:n) admits the following form: q(x0:n|z0:n) =q(xn|x0:n−1,z0:n)q(x0:n−1|z0:n−1), (25) 2Here, we have assumed that the dimensionality of xnis 1. 8 which satisﬁes the probability chain rule decomposition q(x0:n|z0:n) =n/productdisplay k=0q(xn−k|xn−k−1,z0:n−k), (26) where the last term is q(x0|z0) and the second last term is q(x1|x0,z0:1).The weights are then calculated as Wn=p(z0:n|x0:n)p(x0:n) q(x0:n−1|z0:n−1)q(xn|x0:n−1,z0:n). (27) Multiplying and dividing (27) by p(z0:n−1|x0:n−1)p(x0:n−1) we have that Wn=p(z0:n−1|x0:n−1)p(x0:n−1) q(x0:n−1|z0:n−1)/bracehtipupleft /bracehtipdownright/bracehtipdownleft /bracehtipupright previous weight×p(z0:n|x0:n)p(x0:n) [p(z0:n−1|x0:n−1)p(x0:n−1)]q(xn|x0:n−1,z0:n). (28) Thus, we can write (28) as Wn=Wn−1×p(z0:n|x0:n)p(x0:n) [p(z0:n−1|x0:n−1)p(x0:n−1)]q(xn|x0:n−1,z0:n). (29) 5.2. Assumptions on Measurement Model We ﬁrst assume that the measurement model satisﬁes the i) conditional independence and the ii)Markov property. The conditional independence property implies that given that we know xnandxmforn/negationslash=m, the information about zmdoes not tell us anything about znand vice versa. The Markov property implies that the observation znonly depends on xn(please see the measurement model in (2)). This means that p(zn|x0:n) =p(zn|xn). We thus have that p(z0:n|x0:n) =n/productdisplay k=0p(zk|x0:n) =n/productdisplay k=0p(zk|xk) =p(zn|xn)n−1/productdisplay k=0p(zk|xk). (30) 5.3. Assumption on Evolution Model By chain rule of probability density functions, we know that p(x0:n) =n/productdisplay k=0p(xk|x0:k−1). (31) wherep(x0|x0:−1):=p(x0). We now assume that the evolution process is Markovian which means that xkis only dependent on xk−1. This implies that (31) can be rewritten as p(x0:n) =n/productdisplay k=0p(xk|xk−1) =p(xn|xn−1)n−1/productdisplay k=0p(xk)|xk−1). (32) 9 Now, we substitute (30) and (32) into (29) and we obtain Wn=Wn−1×p(zn|xn)/producttextn−1 k=0p(zk|xk) /producttextn−1 k=0p(zk|xk)×p(xn|xn−1)/producttextn−1 k=0p(xk)|xk−1) /producttextn−1 k=0p(xk)|xk−1)×1 q(xn|x0:n−1,z0:n) (33) =⇒Wn=Wn−1×p(zn|xn)p(xn|xn−1) q(xn|x0:n−1,z0:n). (34) Here, we note that we sample Npoints from the one-dimensional distribution q(xn|x(i) 0:n−1,z0:n) (assuming that the dimension of xnis one) at any stage nof the algorithm. This is in stark contrast to (22) where the Npoints were sampled from the n+ 1 dimensional distribution at thenth stage. Having said that, we still need to store Nnsamples at the nth stage (n samples corresponding to each x(i) 0:n−1). Hence, in this case, the sampling complexity is N, but the storage complexity is still Nn. 5.4. Summary of Sequential Importance Sampling algorithm 1. DrawNsamples from the proposed importance distribution: x(i) n∼q(xn|x(i) 0:n−1,z0:n). 2. Determine conditional distributions: p(x(i) n|x(i) n−1), p(zn|x(i) n) andq(x(i) n|x(i) 0:n−1,z0:n) 3. Calculate the unnormalized weights W(i) nas W(i) n=W(i) n−1p(zn|x(i) n)p(x(i) n|x(i) n−1) q(x(i) n|x(i) 0:n−1,z0:n). 4. Obtain normalized weights W(i) n. 5. Estimate posterior distribution as p(x0:n|z0:n)≈/summationtextN i=1W(i) nδ(x0:n−x(i) 0:n). 5.4.1. Initialization The initialization of this algorithm is done as follows: 1.x(i) 0∼q(x0),x(i) 1∼q(x1|x(i) 0,z0:1);i= 1,...,N. 2. W(i) 1=p(z1|x(i) 1)p(x(i) 1|x(i) 0) q(x(i) 1|x(i) 0,z0:1). 10 3. W(i) 1=W(i) 1 /summationtextN i=1W(i) 1. 5.5. Summary of Sequential Importance Sampling (SIS) Filter algorithm If only an estimate of xn(ﬁltering) is required instead of x0:n(smoothing) at time n, the se- quential importance sampling algorithm can be simpliﬁed further via another assumption on the importance sampling distribution q(·). If we assume that q(xn|x0:n−1,z0:n) =q(xn|xn−1,zn), then we only need to store x(i) n−1and generate the particle at next time step from q(xn|x(i) n−1,zn). In this case, the storage complexity is Nand the sampling complexity is also N. For this par- ticular case, the conditional mean of p(xn|z0:n) can be obtained as follows: 1. DrawNsamples from the proposed importance distribution: x(i) n∼q(xn|x(i) n−1,zn). 2. Determine conditional distributions: p(x(i) n|x(i) n−1), p(zn|x(i) n) andq(x(i) n|x(i) n−1,zn). 3. Calculate the unnormalized weights W(i) nas W(i) n=W(i) n−1p(zn|x(i) n)p(x(i) n|x(i) n−1) q(x(i) n|x(i) n−1,zn). 4. Obtain normalized weights W(i) n. 5. Estimate posterior distribution as p(xn|z0:n)≈/summationtextN i=1W(i) nδ(xn−x(i) n). 5.6. Degeneracy Problem In practice, iteration of steps 1 and 3 in the SIS ﬁlter algorithm (deﬁned in Section 5.5) may lead to problems of degeneracy where only a few particles will have signiﬁcant weights. Degeneracy is typically measured by eﬀective sample size Nedeﬁned as Ne:=1 /summationtextN i=1(W(i) n)2. (35) A large degeneracy implies a small value of Ne. A small Neimplies that we are spending too much computational eﬀort to update particles that contribute very small to the posterior approximation. 11 An example of degeneracy If we setq(xn|x(i) n−1,zn) =p(xn|x(i) n−1).Then W(i) n=W(i) n−1·p(zn|x(i) n). In this case, if the prior p(xn|x(i) n−1) is a much broader distribution than the likelihood p(zn|x(i) n), or if both the densities are not “aligned,” then only a few particles will have signiﬁcant weights, and a large number of particles will have weights close to zero. Thus, in this situation, we are eﬀectively using a very small number of particles to represent the posterior distribution. 5.7. Resampling IfNegoes below some preset threshold, we then need to resample. The motivation behind resampling is to retain particles with large weights and discard those with smaller weights. Resampling process involves re-selecting a new set of Nparticles from the current set of N particles, where a new particle is selected with probability equal to its normalized weight. This process is carried out “with replacement,” which means that a particle that has been selected before can be reselected in future. Thus, particles with higher normalized weights are likely to be selected more than once, whereas particles with lower weights might not be reselected at all after resampling is complete. After resampling, each particle is assigned a uniform weight of 1 N. We next look at a common algorithm used for resampling the particles. 5.7.1. Systematic Resampling Algorithm 1. Given a set of particles and weights {x(i) n,W(i) n, i= 1,...,N.} 2. Select a new sample from the set/braceleftBig x(i) n/bracerightBigN i=1with probabilityW(i) nand call it ˆx(i) n. The new sample is assigned the weight ˆW(i) n=1 N. 3. Repeat the resampling process until Nnumber of new particles have been sampled. 5.7.2. A Common Version of Sequential Importance Resampling (SIR) Filter One common version of SIR ﬁlter chooses the prior distribution as a proposal: q(xn|x(i) n−1,zn) = p(xn|x(i) n−1).Then W(i) n=p(zn|x(i) n), (36) and these updates are followed by resampling at each step. Thus the likelihood p(zn|x(i) n) solely decides the weight of the particle iat thenth stage. It is worth noting that there is no W(i) n−1 12 −10−5 0 5 1000.10.20.30.40.5 xnp(xn|x(i) n−1) p(zn|xn) −10−5 0 5 1000.10.20.30.40.5 xnp(xn|x(i) n−1) p(zn|xn) Two diﬀerent scenarios for weight propagation in particle ﬁlters. In the top ﬁgure most particles will have negligible weight: in the other scenario most particles will have signiﬁcant weight as the two densities are “aligned.”Figure 2 . For particles to have any meaningful weights, there should be a signiﬁcant overlap between the prior and likelihood density functions. 13 term in weight update equation (36) since after resampling at time n−1, the weights are set equal to1 N. There is a potential problem with this version of SIR ﬁlter. This has to do with the independence of the proposal or prior distribution p(xn|x(i) n−1) from observations zn. Since the particles are sampled from p(xn|x(i) n−1) and then weighted by the likelihood p(zn|x(i) n), only those particles x(i) nwill have considerable weight that are sampled from the region of signiﬁcant overlap of the two densities (please see Fig. 2 in this regard). Thus, it is more beneﬁcial to update states by taking observations or the likelihood function into account at each stage of timen. 5.7.3. Problems with Resampling: Sample Impoverishment or loss of diversity Particles with large weights are likely to be redrawn multiple times. In extreme case, all particles might “collapse” into a single particle. In this scenario, we are approximating the posterior with just one distinct sample which poorly characterizes the true posterior density. This phenomenon is known as sample impoverishment ordiversity loss of samples or particles. The solution to this problem is to use Bootstrap or SIR ﬁlter with Markov Chain Monte Carlo (MCMC) steps, which is the subject of discussion in the next topic. 5.8. Bootstrap Particle Filter with Markov Chain Monte Carlo (MCMC) steps We can use the SIR ﬁlter in conjunction with MCMC steps to alleviate sample impoverishment issues. After the particles have been resampled, they are perturbed about their locations using MCMC algorithm in order to distribute them to regions of higher posterior probabilities. This is one of the most common techniques to introduce “diversity” in resampled particles [27]. After obtaining resampled particles ˆx(i) n, we obtain a new set of particles ˜x(i) nas follows ˜x(i) n=ˆx(i) n+/epsilon1(i), (37) where /epsilon1(i)∼p(/epsilon1) =N(0,C/epsilon1/epsilon1), whereC/epsilon1/epsilon1is a diagonal matrix. The corresponding acceptance probability for the Metropolis-Hastings approach of MCMC is given by [27] A(˜x(i) n,ˆx(i) n) = min  p(zn|˜x(i) n)p(˜x(i) n|x(i) n−1) p(zn|ˆx(i) n)p(ˆx(i) n|x(i) n−1),1  . (38) This acceptance probability A(˜x(i) n,ˆx(i) n) is then compared to a value that is sampled from a uniform distribution between 0 and 1. If the acceptance probability is greater than the sampled value, we replace ˆx(i) nwith ˜x(i) n. Otherwise, we retain ˆx(i) nas our ﬁnal particle. 14 points sampled from p(xn|xn−1)p(zn|xn)Weighted particlesResamplingpoints sampled from p(xn+1|xn)p(zn+1|xn+1)p(zn+2|xn+2) particles sampled from p(xn+2|xn+1) xn Weight propagation of particles at diﬀerent instants of time.Figure 3 . Weight propagation in SIR ﬁlter 15 Since the posterior is proportional to p(zn|x(i) n)p(x(i) n|x(i) n−1), we note that the acceptance probability compares the fraction of two posterior densities (one corresponding to perturbed particle ˜x(i) n, and the other corresponding to resampled particle ˆx(i) n) with 1. This means that the larger the value of posterior corresponding to the perturbed particle, the greater is the probability of it replacing the resampled particle in the ﬁnal set. In this sense, we only choose a given perturbed particle to represent the posterior if it maximizes the posterior density. We can summarize the basic algorithm in the next section. 5.9. Bootstrap Particle Filter with MCMC steps algorithm 1Initialize: Draw x(i) 0fromp(x0) fori= 1,2,...,N . SetW(i) 0=1 N. 2Importance Sampling: Draw x(i) nfromp(xn|x(i) n−1). 3Weight update: W(i) n=p(zn|x(i) n). 4Weight normalization: W(i) n=W(i) n/summationtextN i=1W(i) n. 5Resampling decision: ComputeNeﬀ=1/summationtextN i=1W(i) n2. IfNeﬀ≤Nthreshold , then resample. Otherwise compute the approximate distribution as follows, p(xn|z1:n)≈N/summationdisplay i=1W(i) nδ(xn−x(i) n), and go to step 2. 6Resampling: Resample ˆx(i) nfrom x(i) nusing Systematic Resampling algorithm. 7Diversify: –Draw /epsilon1(i)fromN(0,C/epsilon1/epsilon1). –Compute ˜x(i) n=ˆx(i) n+/epsilon1(i). –Compute Diversiﬁcation Acceptance Probability as follows: A(˜x(i) n,ˆx(i) n) = min  p(zn|˜x(i) n)p(˜x(i) n|x(i) n−1) p(zn|ˆx(i) n)p(ˆx(i) n|x(i) n−1),1   –DrawukfromU(0,1). 16 –Setx(i) n=˜x(i) nifuk<A(˜x(i) n,ˆx(i) n).Else set x(i) n=ˆx(i) n. 8Distribution: Compute distribution as follows: p(xn|z0:n)≈N/summationdisplay i=1W(i) nδ(xn−x(i) n).