Erdogmus, Deniz

Lehn-Schiøler, Tue

Principe, Jose C.

Using a Parzen density estimator any distribution can be approximated arbitrarily close by a sum of kernels. In particle filtering this fact is utilized to estimate a probability density function with Dirac delta kernels; when the distribution is discretized it becomes possible to solve an otherwise intractable integral. In this work we propose to extend the idea and use any kernel to approximate the distribution. The extra work involved in propagating small kernels through the nonlinear function can be made up for by decreasing the number of kernels needed, especially for high dimensional problems. A further advantage of using kernels with nonzero width is that the density estimate becomes continuous

Online Research Database In Technology

PARZEN PARTICLE FILTERSTue Lehn-SchiølerISP, Technical University of Denmarkemail: tls@imm.dtu.dkDeniz Erdogmus & Jose C. PrincipeCNEL, University of Floridaemail: [deniz,principe]@cnel.ufl.eduABSTRACTUsing a Parzen density estimator any distribution can be ap-proximated arbitrarily close by a sum of kernels. In particlefiltering this fact is utilized to estimate a probability den-sity function with Dirac delta kernels; when the distributionis discretized it becomes possible to solve an otherwise in-tractable integral. In this work we propose to extend theidea and use any kernel to approximate the distribution. Theextra work involved in propagating small kernels throughthe nonlinear function can be made up for by decreasingthe number of kernels needed, especially for high dimen-sional problems. A further advantage of using kernels withnonzero width is that the density estimate becomes contin-uous.1. INTRODUCTIONThe filtering problem can be formulated asxk = f(xk−1) + vk−1 (1a)zk = h(xk) + wk (1b)where v and w are the process noise and the observationnoise. The state transition density is fully specified by f andthe process noise distribution and the observation likelihoodis fully specified by h and the observation noise distribution.p(xk|xk−1) = pv(xk − f(xk−1)) (2a)p(zk|xk) = pw(zk − h(xk)) (2b)The problem is to find an update formula from p(xk−1|z1:k−1)to p(xk|z1:k), where z1:k denotes all observation {z1, . . . , zk}up to time k. The Bayesian approach [1] gives the followingupdate:p(xk|z1:k) = (3)p(zk|xk)∫p(xk|xk−1)p(xk−1|z1:k−1)dxk−1∫p(zk|xk)p(xk|z1:k−1)dxkThe problem can be broken down to two subproblems. 1)Find the propagation of a probability density function (pdf)through a nonlinearity. 2) Modify the pdf according to therecorded measurements zk. Figure 1(a) illustrates stage oneof the problem.In equation (4) multiplying with p(zk|xk) can be seenas stage two, and performing the multi-dimensional integra-tionp(xk|z1:k−1) =∫p(xk|xk−1)p(xk−1|z1:k−1)dxk−1(4)as stage one. In general the integral can not be calculatedanalytically, hence, we need some way of estimating thedistribution p(xk|z1:k−1) .Algorithms fall into four categories: Extended KalmanFilters, Gaussian Sum Filters, Sigma-Point Kalman Filtersand Sequential Monte Carlo Methods (Particle Filters) [2].Another way to categorize the methods is Gaussian belief(Extended Kalman filters, sigma point filters, moment match-ing), mixture of Gaussians (Gaussian-sum filter, pseudo-Bayes) and non-parametric methods (Particle filters) [3]. Inthe Extended Kalman Filter, the distributions are assumedGaussian, but, the functions are not linear. The functionsf and h are linearized around the previous state xk−1 us-ing a second order Taylor expansion and then the standardKalman equations are used. The result is a Gaussian distri-bution for p(xk|z1:k) (see figure 1(b)). For nonlinear sys-tems the solution is better than a normal Kalman filter, andit is accurate to first order. The Unscented Kalman Filter(Sigma Point Filter) [4] propagates points one standard de-viation from the previous state xk−1 through the nonlin-earity, then uses the points weighted appropriately (Gaus-sian quadrature like) to estimate mean and co-variance of aGaussian. Finally this is used in the standard Kalman equa-tions. It is accurate to the second order.If the process noise distribution is approximated by amixture of Gaussians the family of Gaussian sum filtersarises [5]. In the mixture of Gaussians each mixture compo-nent is propagated through an extended Kalman filter. Thestate update f is linearized around the means of each mix-ture component and h is linearized around the predictedvalue for the mean of each mixture (f(xk−1)). The resultingdistribution is again a mixture of Gaussians. If the processnoise is also non-Gaussian, this too can be approximatedV - 7810-7803-8484-9/04/$20.00 ©2004 IEEE ICASSP 2004à á(a) A density propagated through anonlinear function.(b) In the extended Kalman filter,the distribution is assumed to be aGaussian modified by a lineariza-tion of the nonlinearity around theprevious state.(c) In the particle filter the densitiesare approximated by discrete sam-ples. These samples can be sentthrough the nonlinearity to give anestimate of the output density.Fig. 1. Propagation of a pdf through a nonlinearity 1(a)and different approximations to the propagated distribution1(b) and 1(c). This is the prediction step corresponding toequation (1a), the resulting pdf is then modified to matchthe measurements equation (1b).with a mixture of Gaussians. However, in this case the num-ber of mixing components increases quickly.Nonparametric methods are an entirely different approachto nonlinear filtering. In the Particle Filter it is assumedthat the distributions p(xk|z1:k) and p(xk−1|z1:k−1) fromequation (4) can be estimated by discrete distributions (fig-ure 1(c)). Samples are drawn from the posterior distributionusing importance sampling and a proposal distribution. Inthe generic particle filter the transition probability(p(xk|xk−1)) is used as proposal, but other proposals hasbeen proposed in e.g. the extended Kalman particle filterand the unscented particle filter [6], in these methods a filter(ekf or ukf) is calculated for each particle and the resultingmixture of Gaussians is used as proposal distribution for theparticle filter. In an attempt to combine the particle filter andthe Gaussian sum filter the Gaussian Sum Particle Filteringwas proposed [7]. In this approach both the density and theprocess noise is considered mixture of Gaussians, in eachtime step samples are drawn from the mixture approximat-ing p(xk−1|z1:k−1). These samples are propagated throughthe nonlinearity and used to offset the means in a mixturedescribing p(xk|xk−1), then samples are drawn from thisdistribution too. In this way a discrete approximation ofp(xk|z1:k−1) is obtained and the sample mean and covari-ance of the new mixtures can be estimated. Unfortunatelythe number of mixtures explode, to avoid this mixtures withsmall weight can be thrown away. In a similar manner, theGaussian mixture Sigma Point Particle Filter [2] uses a bankof sigma point filters to update p(xk|z1:k−1) then samplesare drawn from the mixture and the importance weights arecalculated before a Gaussian mixture is fitted to produce theposterior estimate.In this paper an algorithm based on the Parzen densityestimator is presented. The algorithm is best categorized asnon-parametric, since it can be seen as a direct extensionto the particle filter. The basic concept is to improve theperformance of the particle filter by using a better densityestimate.The algorithm is similar to the Gaussian Sum ParticleFilter and the Kernel Filter [8], however, it is derived in adifferent manner that allows use of any kernel type. Thederivation of the algorithm uses a sample mean estimate ofthe integral p(xk|z1:k−1) and a particle filter like update ofthe weights. In section 2 the algorithm is derived and insection 3 experimental results are provided.2. KERNEL METHODWith a Parzen density estimator [9, 10] a distribution can beapproximated arbitrarily close by a number of identical ker-nels centered on points chosen from the distribution. In theparticle filter the kernels are delta functions, but informationcan be gained by using a broader kernel.V - 782á áThe distribution at time k − 1 can be approximated by:p(xk−1|z1:k−1) ≈∑Ni wik−1K(Aik−1(xk−1−xik−1)), whereAi is a transformation matrix used to keep track of distor-tions of the kernel. Each kernel can be propagated throughthe mapping p(xk|xk−1) by using a local linearization, yield-ing a continuous output distribution p(xk|z1:k), this is againa sum of kernels but the kernels are no longer identical (inthe sense that they are from the same family of functions,yet they have different parameters).Using the kernel representation equation (4) can be writ-ten as:N∑iwik−1∫pv(xk−f(xk−1))K(Aik−1(xk−1−xik−1))dxk−1(5)Each part of the sum can be handled individually, and underthe assumption that the kernels are small compared to thedynamic in the nonlinearity, f can be locally linearized. Bylinearizing f around xik−1 the jacobian J|xik−1 =∂f∂x |xik−1is introduced and the following change of variables can beemployed: xˆk−1 = xk − f(xik−1)−J|xik−1(xk−1−xik−1).Inserting this in the integral from equation (5) yields:∣∣∣J|xik−1∣∣∣−1 ∫ [pv(xˆk−1) (6)K(Aik−1J|−1xik−1(xk − f(xik−1)− xˆk−1))] dxˆk−1This integral is an expectation over the process noiseEpv[K(Aik−1J|−1xik−1(xk − f(xik−1)− xˆk−1)]and can beapproximated by a sample mean. In the extreme case a sin-gle sample drawn from pv can be used, and the result is atranslation of the kernel by the noise sample:Epv[K(Aik−1J|−1xik−1(xk − f(xik−1)− xˆk−1))]≈ K(Aik−1J|−1xik−1(xk − f(xik−1)− vk−1)),vk−1 ∼ pvWriting p(xk|z1:k) =∑Ni wikK(Aik(xk − xik)) It is pos-sible to identify the centers of the kernels xik = f(xik−1) +vk−1 and the transformation matrix Aik = Aik−1J|−1xik−1.By considering equation (4), (5) and (6) the weight updatecan be found to be wik = wik−1p(zk|xik)∣∣∣J|xik−1∣∣∣−1. Thisderivation holds for any kernel, however, for simplicity, inthis paper the kernels are considered Gaussian.For a Gaussian a change of variables can be employedsuch that the update of Aik can be replaced with an updateof the covariance matrix as follows:Aik = Aik−1 (7a)Σik = J|xik−1Σik−1J|Txik−1(7b)The transformation matrix is absorbed in the covariance ma-trix.Fig. 2. Problem with nonlinear state transition, nonlinearobservation process and Gaussian noise. Note that the per-formance of a Parzen particle filter with≈ 10 kernels equalsthat of a normal particle filter with ≈ 20 kernels.The transformation matrix A (or Σ in the Gaussian case)is distorted in each iteration, to avoid to much distortiona resampling schema can be applied. With a suitable fre-quency the distribution can be re-approximated by a Parzenestimator by drawing samples from p(xk|z1:k), choosing Aor Σ to take their initial values and setting the weights equal.Earlier attempts use the kernels in the resampling phasewhere the shape and kernel size are selected based on theparticle statistics (e.g. covariance) [11] . However, the pro-posed method iterates these properties of the kernel throughthe system equations, thus there is no need for optimizationof kernel parameters at every step. In addition, the approx-imation of the integral stochastically using a sample drawnfrom pv includes an inherent resampling step at every it-eration, which allows the particle filter accuracy to survivelonger than the standard version.3. EXPERIMENTSIn this section the performance of the Parzen particle filterwill be compared to the performance of the standard particlefilter1. The method is tested on a one dimensional problem:xk =xk−12+ 25xk−1(1 + x2k−1)+ 8 cos(1.2k) + vk(8a)zk = 10 arctan(xk10) + wk (8b)Where vk and wk are drawn from Gaussian distributionsG(0, 1) (figure 2) and from gamma distributions Γ(3, 2) (fig-1Code to reproduce the results can be found at www.imm.dtu.dk/∼tlsV - 783á áFig. 3. Problem with nonlinear state transition, nonlinearobservation process and gamma distributed noise.ure 3).The Parzen particle filter and the generic particle filterhas been used on 100 time series generated using equa-tion (8).In figure 2 the mean square error is plotted as a func-tion of the number of kernels. It can be seen that with fewkernels the methods perform equally good (or bad), but asthe number of kernels increases the kernel method becomesbetter. It can be seen that for this one dimensional examplethe kernel methods perform equally well, but the number ofparticles can be reduced drastically by improving the den-sity estimate. It is expected that this effect will be even moreimpressive in higher dimensional problems.4. CONCLUSIONA novel algorithm for nonlinear filtering is presented, thealgorithm is based on Parzen density estimates and particlefilter like propagation of the kernel through local lineariza-tions of the nonlinear function.It is shown that the improved density estimate help per-formance both with Gaussian and non-Gaussian noise. Inthis work only the special case with a Gaussian kernel is ex-amined, however it is expected that a broader kernel wouldbe well suited for long tailed noise, since it will be morelikely to get the particles spread out.The basic formulation for the arbitrary kernel case hasbeen derived and performances of various kernel choiceswill be compared in a future publication.5. REFERENCES[1] N. Gordon, D. Salmond, and A. F. M. Smith, “Novelapproach to non-linear and non-gaussian bayesianstate estimation,” IEE Proceedings-F, vol. 140, pp.107–113, 1993.[2] R. van der Merwe and E. Wan, “Gaussian mixturesigma-point particle filters for sequential probabilisticinference in dynamic state-space models,” Speech andSignal Processing (ICASSP), April 2003.[3] T. Minka, “Bayesian inference in dynamic models –an overview,www.stat.cmu.edu/∼ minka/dynamic.html,” 2003.[4] S. J. Julier and J. Uhlmann, “A new extension of thekalman filter to nonlinear systems,” 1997.[5] D. L. Alspach and H. W. Sorensen, “Nonlinearbayesian estimation using gaussian sum approxima-tions,” IEEE Transactions on Automatic Control, vol.17, no. 4, pp. 439–448, August 1972.[6] R. van der Merwe, N. de Freitas, A. Doucet, andE. Wan, “The unscented particel filter,” Advances inNeural Information Processing Systems 13, Nov 2001.[7] J.H. Kotecha and P.M. Djuric, “Gaussian sum parti-cle filtering for dynamic state space models,” Acous-tics, Speech, and Signal Processing, 2001. Proceed-ings. 2001 IEEE International Conference on, vol. 6,pp. 3465 –3468 vol.6, 2001.[8] M. Hurzeler and H.R. Kunsch, “Monte carlo approx-imations for general state space models,” Journal ofComputational and Graphical Statistics, vol. 7, no. 2,pp. 175–193, 1998.[9] E. Parzen, “On estimation of a probability densityfunction and mode,” Ann. Math. Stat, vol. 27, pp.1065–1076, 1962.[10] L. Devroye and G. Lugosi, Combinatorial Methods inDensity Estimation, Springer, New York, 2001.[11] A. Doucet, N. de Freitas, and N. Gordon, Eds., Se-quential Monte Carlo Methods in Practice, Springer,New York, 2001.V - 784á à

Parzen Particle Filters

http://orbit.dtu.dk/ws/files/2756557/imm2870.pdf

Parzen Particle Filters

Abstract

Similar works

Full text

Available Versions

Online Research Database In Technology