Tutorial on Particle Filters assembled and extended by

Tutorial on Particle Filters assembled and extended by Longin Jan Latecki Temple University, latecki@temple. edu using slides from Keith Copsey, Pattern and Information Processing Group, DERA Malvern; D. Fox, J. Hightower, L. Liao, D. Schulz, and G. Borriello, Univ. of Washington, Seattle Honggang Zhang, Univ. of Maryland, College Park Miodrag Bolic, University of Ottawa, Canada Michael Pfeiffer, TU Gratz, Austria

Outline z Introduction to particle filters – Recursive Bayesian estimation z Bayesian Importance sampling – Sequential Importance sampling (SIS) – Sampling Importance resampling (SIR) z Improvements to SIR – On-line Markov chain Monte Carlo z Basic Particle Filter algorithm z Example for robot localization z Conclusions

Particle Filters z Sequential Monte Carlo methods for on-line learning within a Bayesian framework. z Known as – Particle filters – Sequential sampling-importance resampling (SIR) – Bootstrap filters – Condensation trackers – Interacting particle approximations – Survival of the fittest

History z First attempts – simulations of growing polymers – z M. N. Rosenbluth and A. W. Rosenbluth, “Monte Carlo calculation of the average extension of molecular chains, ” Journal of Chemical Physics, vol. 23, no. 2, pp. 356– 359, 1956. First application in signal processing - 1993 – z N. J. Gordon, D. J. Salmond, and A. F. M. Smith, “Novel approach to nonlinear/non-Gaussian Bayesian state estimation, ” IEE Proceedings-F, vol. 140, no. 2, pp. 107– 113, 1993. Books – A. Doucet, N. de Freitas, and N. Gordon, Eds. , Sequential Monte Carlo Methods in Practice, Springer, 2001. – B. Ristic, S. Arulampalam, N. Gordon, Beyond the Kalman Filter: Particle Filters for Tracking Applications, Artech House Publishers, 2004. z Tutorials – M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-gaussian Bayesian tracking, ” IEEE Transactions on Signal Processing, vol. 50, no. 2, pp. 174– 188, 2002.

Problem Statement z Tracking the state of a system as it evolves over time z Sequentially arriving (noisy or ambiguous) observations z We want to know: Best possible estimate of the hidden variables

Solution: Sequential Update z Storing and processing all incoming measurements is inconvenient and may be impossible z Recursive filtering: –Predict next state pdf from current estimate –Update the prediction using sequentially arriving new measurements z Optimal Bayesian solution: recursively calculating exact posterior density

Particle filtering ideas z Particle filter is a technique for implementing recursive Bayesian filter by Monte Carlo sampling z The idea: represent the posterior density by a set of random particles with associated weights. z Compute estimates based on these samples and weights Posterior density Sample space

Tools needed Recall “law of total probability” and “Bayes’ rule”

Recursive Bayesian estimation (I) z Recursive filter: – System model: – Measurement model: – Information available:

Recursive Bayesian estimation (II) z Seek: – i = 0: filtering. – i > 0: prediction. – i<0: smoothing. z Prediction: – since:

Recursive Bayesian estimation (III) z Update: z where: – since:

Bayes Filters (second pass) Estimating system state from noisy observations System state dynamics Observation dynamics We are interested in: Belief or posterior density

From above, constructing two steps of Bayes Filters Predict: Update:

Assumptions: Markov Process Predict: Update:

Bayes Filter How to use it? What else to know? Motion Model Perceptual Model Start from:

Example 1 Step 0: initialization Step 1: updating

Example 1 (continue) Step 2: predicting Step 3: updating 1 Step 4: predicting

Classical approximations z Analytical methods: – Extended Kalman filter, – Gaussian sums… (Alspach et al. 1971) • Perform poorly in numerous cases of interest z Numerical methods: – point masses approximations, – splines. (Bucy 1971, de Figueiro 1974…) • Very complex to implement, not flexible.

Perfect Monte Carlo simulation z Recall that z Random samples are drawn from the posterior distribution. z Represent posterior distribution using a set of samples or particles. z Easy to approximate expectations of the form: – by:

Random samples and the pdf (I) z Take p(x)=Gamma(4, 1) z Generate some random samples z Plot histogram and basic approximation to pdf 200 samples

Random samples and the pdf (II) 500 samples 1000 samples

Random samples and the pdf (III) 5000 samples 200000 samples

Importance Sampling z Unfortunately it is often not possible to sample directly from the posterior distribution, but we can use importance sampling. z Let p(x) be a pdf from which it is difficult to draw samples. z Let xi ~ q(x), i=1, …, N, be samples that are easily generated from a proposal pdf q, which is called an importance density. z Then approximation to the density p is given by where

Bayesian Importance Sampling z By drawing samples distribution from a known easy to sample proposal we obtain: where are normalized weights.

Sequential Importance Sampling (I) z Factorizing the proposal distribution: z and remembering that the state evolution is modeled as a Markov process z we obtain a recursive estimate of the importance weights: z Factorizing is obtained by recursively applying

Sequential Importance Sampling (SIS) Particle Filter SIS Particle Filter Algorithm for i=1: N Draw a particle Assign a weight end (k is index over time and i is the particle index)

Derivation of SIS weights (I) z The main idea is Factorizing : and Our goal is to expand p and q in time t

Derivation of SIS weights (II)

Derivation of SIS weights (II) and under Markov assumptions

SIS Particle Filter Foundation z At each time step k z Random samples are drawn from the proposal distribution for i=1, …, N z They represent posterior distribution using a set of samples or particles z Since the weights are given by z and q factorizes as

Sequential Importance Sampling (II) z Choice of the proposal distribution: z Choose proposal function to minimize variance of z Although common choice is the prior distribution: We obtain then (Doucet et al. 1999):

Sequential Importance Sampling (III) z Illustration of SIS: z Degeneracy problems: – variance of importance ratios increases stochastically over time (Kong et al. 1994; Doucet et al. 1999). – In most cases then after a few iterations, all but one particle will have negligible weight

Sequential Importance Sampling (IV) z Illustration of degeneracy:

SIS - why variance increase z Suppose we want to sample from the posterior – choose a proposal density to be very close to the posterior density • Then • and z So we expect the variance to be close to 0 to obtain reasonable estimates – thus a variance increase has a harmful effect on accuracy

Sampling-Importance Resampling z SIS suffers from degeneracy problems so we don’t want to do that! z Introduce a selection (resampling) step to eliminate samples with low importance ratios and multiply samples with high importance ratios. z Resampling maps the weighted random measure on to the equally weighted random measure – by sampling uniformly with replacement from with probabilities z Scheme generates satisfies: children such that and

Basic SIR Particle Filter - Schematic Initialisation measurement Resampling step Importance sampling step Extract estimate,

Basic SIR Particle Filter algorithm (I) z Initialisation – – For sample – and set z Importance Sampling step – For sample and set – For compute the importance weights wik – Normalise the importance weights,

Basic SIR Particle Filter algorithm (II) z Resampling step – Resample with replacement particles: – from the set: – according to the normalised importance weights, z Set – proceed to the Importance Sampling step, as the next measurement arrives.

Resampling x

Generic SIR Particle Filter algorithm M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters …, ” IEEE Trans. on Signal Processing, 50( 2), 2002.

Improvements to SIR (I) z Variety of resampling schemes with varying performance in terms of the variance of the particles : – Residual sampling (Liu & Chen, 1998). – Systematic sampling (Carpenter et al. , 1999). – Mixture of SIS and SIR, only resample when necessary (Liu & Chen, 1995; Doucet et al. , 1999). z Degeneracy may still be a problem: – During resampling a sample with high importance weight may be duplicated many times. – Samples may eventually collapse to a single point.

Improvements to SIR (II) z To alleviate numerical degeneracy problems, sample smoothing methods may be adopted. – Roughening (Gordon et al. , 1993). • Adds an independent jitter to the resampled particles – Prior boosting (Gordon et al. , 1993). • Increase the number of samples from the proposal distribution to M>N, • but in the resampling stage only draw N particles.

Improvements to SIR (III) z Local Monte Carlo methods for alleviating degeneracy: – Local linearisation - using an EKF (Doucet, 1999; Pitt & Shephard, 1999) or UKF (Doucet et al, 2000) to estimate the importance distribution. – Rejection methods (Müller, 1991; Doucet, 1999; Pitt & Shephard, 1999). – Auxiliary particle filters (Pitt & Shephard, 1999) – Kernel smoothing (Gordon, 1994; Hürzeler & Künsch, 1998; Liu & West, 2000; Musso et al. , 2000). – MCMC methods (Müller, 1992; Gordon & Whitby, 1995; Berzuini et al. , 1997; Gilks & Berzuini, 1998; Andrieu et al. , 1999).

Improvements to SIR (IV) z Illustration of SIR with sample smoothing:

Ingredients for SMC z Importance sampling function – Gordon et al – Optimal – UKF pdf from UKF at – Gordon et al SIR – Liu & Chen Residual – Carpenter et al Systematic z Redistribution scheme – Liu & Chen, Doucet et al Resample when necessary z Careful initialisation procedure (for efficiency)

Particle filters z Also known as Sequential Monte Carlo Methods z Representing belief by sets of samples or particles are nonnegative weights called importance factors z Updating procedure is sequential importance sampling with re-sampling z

Example 2: Particle Filter Step 0: initialization Each particle has the same weight Step 1: updating weights. Weights are proportional to p(z|x)

Example 2: Particle Filter Step 2: predicting. Predict the new locations of particles. Step 3: updating weights. Weights are proportional to p(z|x) Step 4: predicting. Particles are more concentrated in the region where the person is more likely to be Predict the new locations of particles.

Compare Particle Filter with Bayes Filter with Known Distribution Updating Example 1 Example 2 Predicting Example 1 Example 2

Particle Filters

Sensor Information: Importance Sampling

Robot Motion

Sensor Information: Importance Sampling

Robot Motion

Global Localization of Robot with Sonar http: //www. cs. washington. edu/ai/Mobile_Robotics/mcl/animations/global-floor. gif

Tracking in 1 D: the blue trajectory is the target. The best of 10 particles is in red.

Matlab code: truex is a vector of 100 positions to be tracked.

Application Examples z Robot localization z Robot mapping z Visual Tracking –e. g. human motion (body parts) z Prediction of (financial) time series –e. g. mapping gold price to stock price z Target recognition from single or multiple images z Guidance of missiles z Contour grouping z Nice video demos: http: //www. cs. washington. edu/ai/Mobile_Robotics/mcl/

2 nd Book Advert z z z Statistical Pattern Recognition Andrew Webb, DERA ISBN 0340741643, Paperback: 1999: £ 29. 99 Butterworth Heinemann z Contents: – Introduction to SPR, Estimation, Density estimation, Linear discriminant analysis, Nonlinear discriminant analysis - neural networks, Nonlinear discriminant analysis - statistical methods, Classification trees, Feature selction and extraction, Clustering, Additional topics, Measures of dissimilarity, Parameter estimation, Linear algebra, Data, Probability theory.

Homework z Implement all three particle filter algorithms SIS Particle Filter Algorithm (p. 27) Basic SIR Particle Filter algorithm (p. 39, 40) Generic SIR Particle Filter algorithm (p. 42) z and evaluate their performance on a problem of your choice. z Groups of two are allowed. z Submit a report and a ready to run Matlab code (with a script and the data). z Present a report to the class.