Analysis of a sequential Monte Carlo method for ... - CiteSeerX

ARTICLE IN PRESS Signal Processing 90 (2010) 1609–1622

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Analysis of a sequential Monte Carlo method for optimization in dynamical systems Joaquı´n Mı´guez Department of Signal Theory & Communications, Universidad Carlos III de Madrid, Avenida de la Universidad 30, 28911 Legane´s, Madrid, Spain

a r t i c l e in fo

abstract

Article history: Received 16 April 2009 Received in revised form 31 August 2009 Accepted 3 November 2009 Available online 12 November 2009

We investigate a recently proposed sequential Monte Carlo methodology for recursively tracking the minima of a cost function that evolves with time. These methods, subsequently referred to as sequential Monte Carlo minimization (SMCM) procedures, have an algorithmic structure similar to particle filters: they involve the generation of random paths in the space of the signal of interest (SoI), the stochastic selection of the fittest paths and the ranking of the survivors according to their cost. In this paper, we propose an extension of the original SMCM methodology (that makes it applicable to a broader class of cost functions) and introduce an asymptotic-convergence analysis. Our analytical results are based on simple induction arguments and show how the SoI-estimates computed by a SMCM algorithm converge, in probability, to a sequence of minimizers of the cost function. We illustrate these results by means of two computer simulation examples. & 2009 Elsevier B.V. All rights reserved.

Keywords: Sequential Monte Carlo Stochastic optimization Nonlinear dynamics Nonlinear tracking Dynamic optimization

1. Introduction We address the extension and analysis of a recently proposed family of algorithms termed ‘‘cost reference particle filters’’ (CRPFs) that can be used for nonlinear tracking in dynamical systems. The original method was proposed in [21,22], with some theoretical developments in [19] and applications to target tracking in [20,3,9]. Although the new methodology was introduced as an algorithmic generalization of conventional particle filters (and, indeed, the standard particle filtering techniques can be interpreted as CRPFs from a purely procedural viewpoint [19]), a CRPF can be plainly described as a sequential Monte Carlo algorithm that tracks the minima of a cost function that evolves with time. To be precise, assume that we wish to sequentially estimate a discrete-time random signal of interest (SoI), fxt gt2N , from a series of associated observations, fyt gt2N . The estimation criterion is the minimization of a given cost function that depends

E-mail address: [email protected] 0165-1684/$ - see front matter & 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2009.11.007

on both sequences, denoted Ct ðx1 ; . . . ; xt ; y1 ; . . . ; yt Þ, i.e., our desired estimate of the SoI is the sequence x^ 1 ; . . . ; x^ t that minimizes Ct for the observed data. A CRPF generates random paths in the space of fxt gt2N , then uses a stochastic mechanism to select the most promising ones (and discard the others) and, finally, ranks the surviving paths according to their cost. The resemblance to particle filters is apparent, but the ultimate goal of the algorithm is the minimization of the prescribed cost and, for this reason, in this paper we adopt the more accurate term ‘‘sequential Monte Carlo minimization’’ (SMCM) for procedures of this type. The main contribution of this paper is the analysis of the convergence of SMCM methods, which was missing in [21,19]. We prove, using simple induction arguments, that an adequately constructed SMCM algorithm yields estimates of the SoI that converge asymptotically (with time and the number of samples in the signal space), and in probability, to a certain sequence of minimizers of the prescribed cost. Additionally, we provide a description of the SMCM class of algorithms that is concise but more general than the statements in [21,19]. Namely, it

ARTICLE IN PRESS 1610

J. Mı´guez / Signal Processing 90 (2010) 1609–1622

encompasses a larger class of cost functions1 that can be used in the algorithm, including many instances that depend on the local dynamics of the random SoI. This is a significant departure from the more restrictive assumptions in [21,19]. We illustrate our results by way of two simulation examples. In the first one, we consider a typical nonlinear dynamical system that has been studied by several authors in the context of nonlinear filtering (see [16] and references therein). We address the tracking of the state of this system by applying the SMCM methodology with a cost function that depends explicitly on the time evolution of the state (this was not strictly feasible with the methods described in [21,19]) but differs significantly from the actual system dynamics. The design of this example is fully consistent with the assumptions of our analysis and, hence, we also use it to illustrate how convergence is achieved. In the second example we study a more complex model that describes the dynamics of a chaotic CO2 laser [18]. The state of the laser includes five dynamic variables and there is only one observable (the laser intensity). We apply the SMCM methodology using four different cost functions and compare the resulting performance with a standard particle filter. In particular, we show that the standard particle filter is very sensitive to mismatches in the assumed form of the observations, while the SMCM algorithms provide a more stable performance. The rest of the paper is organized as follows. We start with a brief outline of our notation in Section 2. The general description of the SMCM class of algorithms is given in Section 3. Section 4 is devoted to the asymptotic convergence analysis. The numerical examples are presented in Section 5 and, finally, Section 6 contains some concluding remarks.

2. Notation We write boldface lower-case letters for column vectors, e.g., x. Sets are denoted with upper-case calligraphic letters, e.g., I , while lower-case greek letters are used to represent probability density functions (pdf’s) or probability mass functions (pmf’s), e.g., p. The joint pdf/pmf, p, of two random vectors x and y is written as pðx; yÞ, and the corresponding conditional pdf/pmf is denoted as pðxjyÞ. A random sample, xðkÞ , drawn from the distribution characterized by the pdf/pmf p, is denoted as xðkÞ # pðxÞ. Probability measures are denoted as P. We may write Pp to make explicit that the associated pdf/pmf is p. For sequences of measures indexed by discrete time t 2 N, we write Ptp , implying that the associated pdf/pmf is pt . The Borel s$algebra in the n-dimensional space Rn is denoted as Bn , while the Lebesgue measure of a measurable set A 2 Bn is distinctly denoted as ‘ðAÞ. 1 It is not the goal of this paper, however, to discuss the degree of suitability of one cost function or the other. This is, in general, dependent on the SoI and the observations and it is the task of the user to choose the type of cost according to her knowledge of the problem at hand.

The SoI at time t 2 N is denoted as xt all through the paper (except in Section 5.2 where we also use zn to denote xt at t ¼ 10n). Notation xt1 :t2 is shorthand for the set fxt1 ; . . . ; xt2 g. The sequences x0 ; x1 ; . . . ; xt ; . . . and x1 ; x2 ; . . . ; xt ; . . . are compactly denoted as fxt gN& and fxt gN , respectively (N& ¼ N [ f0g is the set of non-negative integers). 3. Sequential Monte Carlo minimization 3.1. Cost functions We are interested in tracking the state of a dynamical system that generates a sequence of vector-valued observations (often called ‘‘measurements’’) fyt 2 Rny gN , where ny Z1 is the dimension of the observation space and t denotes discrete time. The SoI is the system state, denoted by the sequence fxt 2 Rnx gN& , where nx Z 1 is the dimension of the signal space. If one represents the dynamics of the state signal, fxt gN& , and the dependence of the measurements, fyt gN , on fxt gN& , by means of a Markov state-space model, then the estimation of xt given the observations y1:t is the typical kind of problem addressed by particle filters and other stochastic filtering techniques [11,24,1]. In this paper we adopt a different perspective. We assume an estimation criterion that consists in the minimization of a sequence of cost functions that involve both the observations and the states. Specifically, at each time t we have the ability to evaluate the cost Ct ðx0:t0 ; y1:t Þ, where y1:t are the available measurements and x0:t0 2 Rnx ' ( ( ( ' Rnx is any trial sequence in the pathspace of the SoI. Our ultimate goal is to find the sequence of states that minimizes the cost at each time step. To be specific, the cost at time t is a function of ðt þ1Þnx þ tny real arguments, namely Ct : Rnx ' ( ( ( ' Rnx ' Rny ' ( ( ( ' Rny -½0; 1Þ;

x0 ; . . . ; xt ; y1 ; . . . ; yt *Ct ðx0:t ; y1:t Þ;

ð1Þ

and we aim at solving the optimization problems S t ðy1:t Þ9argminCt ðx0:t ; y1:t Þ; x0:t

t ¼ 0; 1; 2; . . . ;

ð2Þ

where S t ðy1:t Þ denotes the set of sequences in the pathspace of the SoI that minimize the cost Ct for the given observations y1:t . Note that: (a) S 0 does not depend on any measurement (not available at that time) and (b) the cost functions can present many global minima, hence S t ðy1:t Þ is not necessarily a singleton. We wish to search for elements of S t ðy1:t Þ using a sequential and recursive algorithm and, accordingly, we restrict the cost functions to be considered here to be those that can be computed recursively. Specifically, if we let Cð0Þ be shorthand for C0 ðx0 Þ and, in general, Cðt$1Þ 9Ct$1 ðx0:t$1 ; y1:t$1 Þ, then we assume that CðtÞ can be recursively computed from Cðt$1Þ as CðtÞ ¼ Ct ðx0:t ; y1:t Þ9Hðct ðxt$1:t ; yt Þ; Cðt$1Þ Þ;

ð3Þ

where ct : Rnx ' Rnx ' Rny -½0; 1Þ is a marginal cost function (that depends only on the t-th observation) and

ARTICLE IN PRESS J. Mı´guez / Signal Processing 90 (2010) 1609–1622

function H : ½0; 1Þ ' ½0; 1Þ-½0; 1Þ incorporates the marginal cost into the overall cost. Typical costs can be additive, e.g., Hðct ðxt$1:t ; yt Þ; Cðt$1 ÞÞ ¼ Cðt$1Þ þct ðxt$1:t ; yt Þ, multiplicative, e.g., Hðct ðxt$1:t ; yt Þ; Cðt$1Þ Þ ¼ Cðt$1Þ ct ðxt$1:t ; yt Þ, or present other forms. The marginal cost function, ct , may have multiple minima. We do not require continuity (unless explicitly stated) and only the ability to compute ct ðxt$1:t ; yt Þ and update CðtÞ given Cðt$1Þ is taken for granted. 3.2. Algorithm The minimization of Ct ðx0:t ; y1:t Þ can be attempted using an algorithm similar to the standard particle filter of [14]. In order to describe the proposed methodology we need to introduce additional notation. Specifically:

+ Let Mt denote the number of samples (also called

+ + + + +

particles) drawn at time t. The k-th particle at this time, nx drawn denoted xðkÞ t , is a point in the state space R from a probability distribution (to be introduced below). ðkÞ ðkÞ Let CðtÞ 9Hðct ðxðkÞ t$1:t ; yt Þ; Cðt$1Þ Þ be shorthand for the cost of the k-th stream of particles, xðkÞ 0:t , given the collected measurements y1:t . Let X 0 2 Bnx denote a finite measurable set, hence ‘ðX 0 Þ o1, termed the prior set. Let Bt be a pmf defined over the set f1; . . . ; Mt$1 g. We refer to it as the selection pmf. Let s : ½0; 1Þ-½0; 1Þ be termed the selection function. Let fr~ t ðxt jxt$1 ÞgN be a family of propagation pdf’s, used to generate new particles in the state space.

The recursive SMCM algorithm investigated in this paper can be outlined as follows. (1) Initialization: Draw M0 random samples from the prior set X 0 . The initial costs take a non-negative constant ðkÞ ¼ Cð0Þ Z0 for all k. value, Cð0Þ ðkÞ Mt (2) Recursive step: Let Ot 9fxðkÞ 0:t ; CðtÞ gk ¼ 1 be the collection of particles and costs at time t. (a) Selection: Draw Mt þ 1 random indices from f1; . . . ; Mt g according to the pmf Bt þ 1 , i.e., ðkÞ

i

# Bt þ 1 ðiÞ;

k ¼ 1; . . . ; Mt þ 1 ; i 2 f1; . . . ; Mt g:

ð4Þ

1611

particles with the lowest cost, although averaging techniques can also be used [21,19]. The selection pmf, Bt , and the selection function, sð(Þ, may be assigned different forms by the user of the algorithm. The selection step above is the counterpart of resampling [12,10] in particle filters, hence its goal is to propagate the most promising particles while discarding high-cost ones. An adequate combination of sð(Þ and Bt þ 1 ð(Þ can yield any of the resampling procedures commonly used in particle filtering, as well as other schemes. For instance, if we define normalized weights that decrease exponentially with the costs, ðiÞ g, and choose a constant selection function, wtðiÞ pexpf$CðtÞ sðcÞ ¼ 1=M, then we reproduce exactly the resampling scheme of the standard particle filter with a likelihood of the exponential family [14]. Ref. [19] contains a study of selection methods. The notation is slightly different from the one in this paper, but it is straightforward (and actually convenient) to rewrite the algorithms in [19] using functions sð(Þ and Bt þ 1 ð(Þ. For example, the selection scheme of [19, Section 4] is achieved with sðcÞ ¼ c and ðiÞ ðkÞ Bt þ 1 ðiÞpðCðtÞ $ mink2f1;...;Mt g CðtÞ þ Mt$1 Þ$3 . For the purpose of our convergence analysis in Section 4, we will assume that every particle at time t has a non-zero probability to be selected for propagation at time t þ 1, i.e., that Bt þ 1 ðiÞ 4 0 for all i 2 f1; . . . ; Mt g. No particular form of sð(Þ will be assumed. An important feature of the proposed formulation of the method is the definition of the marginal cost component, ct . We allow for the marginal cost to vary with time and depend on the pair xt$1:t , while in [21,19] this term can only depend on xt and has a fixed functional form. The family of functions Ct that can now be represented is larger and, in particular, it includes costs that depend on the dynamics of fxt gN& . For (an extreme) example, if ct ðxt$1:t ; yt Þ ¼ Jxt $ yt J þ105 Jxt $ xt$1 J (where J ( J denotes the norm of a vector), then any particle xðkÞ t that departs significantly from xðkÞ t$1 will be heavily penalized when computing the cost. Obviously, the notation ct ðxt$1:t ; yt Þ also allows for simpler functions that depend only on xt and not on xt$1 .

3.3. Particle filtering and SMCM

~ ~ xðkÞ t þ 1 # r t þ 1 ðxt þ 1 jx t Þ;

ð5Þ

A particle filter (PF) is a Monte Carlo algorithm for the recursive approximation of the sequence of posterior probability distributions of the state of a Markov dynamic model [12,11,6,8]. In this section we briefly review the class of PFs derived from the sequential importance resampling methodology [12] and relate it to the SMCM framework of this paper. Consider the state-space model described by

ðkÞ ~ ðkÞ CðtðkÞþ 1Þ ¼ Hðct þ 1 ðx~ t ; xðkÞ t þ 1 ; yt þ 1 Þ; C ðtÞ Þ;

ð6Þ

x0 # tðx0 Þ;

ðkÞ x~ 0:t

ðkÞ C~ ðtÞ

ðkÞ

ði Þ Set ¼ x0:t and ðkÞ ~ ðkÞ Mt þ 1 ~ O t 9fx~ 0:t ; C ðtÞ gk ¼ 1 .

ðkÞ

ði Þ ¼ sðCðtÞ Þ

in order to build

(b) Propagation: Draw new particles and compute new costs, ðkÞ

Mt þ 1 ðkÞ Ot þ 1 ¼ fxðkÞ 0:t þ 1 ; Cðt þ 1Þ gk ¼ 1 ,

in order to build ~ ðkÞ ðkÞ xðkÞ 0:t þ 1 ¼ fx 0:t ; xt þ 1 g.

where

Estimation using the particles in Ot þ 1 ¼ fxðkÞ 0:t þ 1 ; Mt þ 1 can be as simple as choosing the sequence of CðtðkÞþ 1Þ gk ¼ 1

xt # tðxt jxt$1 Þ

and

yt # lðyt jxt Þ;

ð7Þ

where tðx0 Þ is the prior pdf of the state signal fxt gN& , tðxt jxt$1 Þ is the transition pdf that determines the system

dynamics and lðyt jxt Þ is the conditional pdf of the observations yt (i.e., the likelihood of xt ). The posterior pdf of x0:t given the observations y1:t (denoted p) can be

ARTICLE IN PRESS J. Mı´guez / Signal Processing 90 (2010) 1609–1622

1612

recursively decomposed, using Bayes’ theorem, as

pðx0:t jy1:t Þplðyt jxt Þtðxt jxt$1 Þpðx0:t$1 jy1:t$1 Þ:

ð8Þ

It is possible to use (8) and the importance sampling principle (IS) [7] to recursively generate properly weighted samples in the state space [12]. Specifically, consider a proposal pdf cðx0:t jy1:t Þ that can be factorized as

cðx0:t jy1:t Þpcðxt jxt$1 ; yt Þcðx0:t$1 jy1:t$1 Þ:

ð9Þ

Then, we can draw samples in the state space sequentially, i.e., given the stream xðiÞ 0:t$1 we can draw ðiÞ xðiÞ t # cðxt jxt$1 ; yt Þ

ð10Þ xðiÞ 0:t

ðiÞ ¼ fx0:t$1 ; xtðiÞ g.

and append the new sample to obtain We can also compute proper importance weights recursively, namely wtðiÞ p

ðiÞ pðx0:t jy1:t Þ

ðiÞ cðx0:t jy1:t Þ

p

ðiÞ lðyt jxtðiÞ ÞtðxðiÞ t jxt$1 Þ ðiÞ wt$1 : ðiÞ ðiÞ cðxt jxt$1 ; yt Þ

ð11Þ

Eqs. (10) and (11) yield the sequential importance sampling (SIS) algorithm [12]. Each weighted sample ðiÞ ; wðiÞ to as a particle and fx0:t t g, for i ¼ 1; . . . ; Mt , is referred PMt ðiÞ the weights are normalized, i ¼ 1 wt ¼ 1. A practical problem of this procedure is the so-called degeneracy of the weights [12], that renders the algorithm useless after a few time steps. This difficulty is overcome by introducing resampling steps [4,12,10,19]. In its conceptually simpler form, resampling consists in drawing Mt times from the discrete distribution that assigns probability mass wðiÞ t to the i-th particle, and then assign uniform weights to the resampled particles (this procedure is known as multinomial resampling [10]). The SIS algorithm combined with resampling steps is referred to as sequential importance resampling (SIR). The standard PF introduced in [14], also termed bootstrap filter (BF), is obtained from the general SIR scheme by choosing cðxt jxt$1 ; yt Þ ¼ tðxt jxt$1 Þ and resampling at every time step. This yields a very simple algorithm that can be outlined as follows (assume Mt ¼ M for all t). xðiÞ 0

w0ðiÞ

(1) Initialization: Draw # tðx0 Þ and set ¼ 1=M, for i ¼ 1; . . . ; M. (2) Recursive step: Given fxtðiÞ ; wtðiÞ gM i ¼ 1, ðiÞ (a) resample to obtain a new set of particles fx~ t gM i¼1 with equal weights, ðiÞ ðiÞ ðiÞ (b) draw new samples xt þ 1 # tðxt þ 1 jx~ t Þ, i ¼ 1; . . . ; M, and ðiÞ (c) assign new weights, wðiÞ t þ 1 plðyt þ 1 jxt þ 1 Þ. Note that we do not keep the complete streams of ðiÞ samples, x0:t , in this simple algorithm. Most integrals with respect to (w.r.t.) the filtering pdf pðxt jy1:t Þ can be easily approximated by summations using the particles generated by the BF. Specifically, if f is some integrable function of the state xt then we can write Z M X f ðxt Þpðxt jy1:t Þdxt , f ðxtðiÞ ÞwtðiÞ : ð12Þ i¼1

See [17,1] for a survey of convergence results of PFs.

Any SIR algorithm can be easily expressed within the SMCM framework. In particular, a generic SIR procedure with proposal pdf cðxt jxt$1 ; yt Þ, multinomial resampling and M particles can be interpreted as a SMCM algorithm with Mt ¼ M,

rt ðxt jxt$1 Þ ¼ cðxt jxt$1 ; yt Þ;

ð13Þ

Hða; bÞ ¼ a þb;

ð14Þ

and ct ðxt$1:t ; yt Þ ¼ logðcðxt jxt$1 ; yt ÞÞ $ logðlðyt jxt ÞÞ $ logðtðxt jxt$1 ÞÞ:

ð15Þ

The corresponding importance weights can be recovered ðiÞ as wðiÞ t pexpf$CðtÞ g and, as a consequence, the multinomial resampling is implemented within the SMCM framework with BtðiÞþ 1 pwtðiÞ and sðcÞ ¼ 1. The two algorithms (SIR and SMCM) become identical if we draw initial samples at time t ¼ 0 from the same distribution, tðx0 Þ, and resample/ select in the same time steps. 4. Asymptotic convergence 4.1. Cost minimizers Given the statement of the SMCM algorithm in Section 3, the straightforward question to pose is whether the algorithm can produce a sequence of state estimates at a given time t with a cost which is arbitrarily close to the minimum one. In particular, let xo0:tjt ¼ fxo0jt ; xo1jt ; . . . ; xotjt g 2 S t ðy1:t Þ

ð16Þ

be a sequence of states (possibly not unique) that minimizes the overall cost for the sequence of observations y1:t and let t x^ 0:tjt ¼ fx^ 0jt ; x^ 1jt ; . . . ; x^ tjt g 2 S M t ðy1:t Þ;

ð17Þ

t SM t ðy1:t Þ

where is the set of approximate minimizers generated by the SMCM algorithm, i.e., t SM t ðy1:t Þ9arg

min

M

x0:t 2fxðkÞ g t 0:t k ¼ 1

Ct ðx0:t ; y1:t Þ:

ð18Þ

Can we guarantee that x^ 0:tjt -xo0:tjt , for some x^ 0:tjt 2 o t SM t ðy1:t Þ and x0:tjt 2 S t ðy1:t Þ? This is a hard problem in general, unless we impose further constraints on the cost function Ct . Assume that T observations, y1:T , are to be collected. In general, the optimal value of the SoI at time t o T, xotjT , cannot be computed before time T. Since any ðMt Þ at time t SMCM algorithm draws the samples xð1Þ t ; . . . ; xt and does not modify them at any later time, it is difficult to ensure the convergence of the sequence x^ 0:TjT except for cost functions with suitable structures (e.g., when the overall cost is additive and the marginal cost at time t depends only on xt and not on xt$1 ). The original question on the convergence of x^ 0:t can be simplified, but still retain a good deal of practical interest, if we restrict our attention to the latest realization of the system state, xt . Thus, we can ask whether a suitably designed SMCM algorithm can ensure that x^ tjt -xotjt in some statistical sense. Our answer is, to some extent, positive. Indeed, we can show that, by adequately

ARTICLE IN PRESS J. Mı´guez / Signal Processing 90 (2010) 1609–1622

choosing the family of propagation densities fr~ t gN , a SMCM method can produce particles arbitrarily close to the minimizer xotjt at time t with high probability.

drawn from a mixture density, xtðiÞ # rt ðxt Þ9

4.2. Assumptions and notations Let us consider an arbitrary but fixed series fxotjt gN& , where xotjt is the ðt þ 1Þ$th vector in some sequence xo0:tjt 2 S t ðy1:t Þ (note that there may exist several possible sequences of minimizers to choose from, since the sets fS t ðy1:t ÞgN are not necessarily singletons). All through this section, as well as in the proofs in Appendix A, we assume that the sequence of observations, fyt gN , is unknown but fixed. This is a rather common assumption in the analysis of particle filters (see, e.g., [5,17,1]) and, in our case, it implies that the minimizers fxotjt gN are deterministic, as they depend only on the observations and the form of the cost function. In the subsequent analysis, randomness is only due to the selection and propagation steps of the SMCM algorithm, not to the observations. We equip the set Rnx with a proper distance dð(; (Þ and introduce a sequence of functions at : Rnx -Rnx , t 2 N& , in the state space. The role of the latter functions is to provide a (possibly rough) approximation to the dynamics of the sequence of minimizers fxotjt gN . Specifically, we assume that the at ’s comply with the following statements. Assumption 1. For any set A - Rnx , let

at ðAÞ9fx 2 Rnx : x ¼ aðx0 Þ for some x0 2 Ag

ð19Þ

denote the image of A under the function at . We assume that, if A 2 Bnx , then at ðAÞ 2 Bnx for all t. Assumption 2. There exists a finite constant A 2 R þ such that At 9dðat ðxotjt Þ; xotþ 1jt þ 1 Þ o A for all t. Assumption 3. For any pair x; x0 2 Rnx and all t 2 N& , if dðx; x0 Þ o 1 then dðat ðxÞ; at ðx0 ÞÞ o1. Next, we use the newly introduced functions fat gN& to construct the propagation densities of the SMCM algorithm. In particular, we assume that ðkÞ r~ t ðxt jxðkÞ t$1 Þ ¼ gðxt ; at$1 ðxt$1 Þ; rt Þ;

ð20Þ

where gðx; c; rÞ denotes a pdf which is

+ strictly positive in the open ball with center c 2 Rn + +

x

and radius r 40, denoted Bðc; rÞ9fx 2 Rnx : dðx; cÞ o rg, i.e., gðx; c; rÞ 40 8x 2 Bðc; rÞ, uniformly continuous in Bðc; rÞ and R null outside Bðc; rÞ, i.e., Bðc;rÞ gðx; c; rÞdx ¼ 1.

Many densities comply with these assumptions, e.g., the uniform pdf in Bðc; rÞ, denoted UðBðc; rÞÞ, or a Gaussian density truncated outside the set Bðc; rÞ. Then, we assume that the selection pmf, Bt ðkÞ, k 2 f1; . . . ; Mt$1 g, is strictly positive (implying that all particles at time t $ 1 have a non-zero probability to be propagated to time t) and formally combine the selection and propagation stages of the algorithm into a single step. The particles are then

1613

M t$1 X

k¼1

Bt ðkÞgðxt ; at ðxðkÞ t$1 Þ; rt Þ;

i ¼ 1; . . . ; Mt ð21Þ

for some prescribed rt 40. The sequence of radii frt gN must be selected by the user and we will show that, if every rt is finite but sufficiently large, then the SMCM algorithm is guaranteed to converge in probability. 4.3. Asymptotic convergence in probability We will show that a SMCM algorithm can be designed in such a way that, for sufficiently large natural numbers fMt gN , the state-space samples x&tjt 9arg

min

M

t x2fxðkÞ t gk ¼ 1

dðx; xotjt Þ

ð22Þ

converge to xotjt , in probability, when t-1. The main convergence result in this paper is stated below. Proposition 1. If we consider

+ a fixed sequence of minimizers fxotjt gN , + a sequence of functions fat gN that satisfy Assumptions 1, &

+ +

2 and 3 for the sequence fxotjt gN , a sequence of propagation pdf’s of the form in Eq. (21) (this implies strictly positive selection pmf’s, Bt ðiÞ 40 for all t and all i 2 f1; . . . ; Mt$1 g), and a sequence of sufficiently large finite radii frt gN ,

then limt-1 limMt -1 dðx&tjt ; xotjt Þ ¼ 0 in probability. To be specific, convergence of the latter limit in probability means that for any arbitrarily small e; d 4 0 there exist te 2 N and a sequence of natural numbers fMt;e;d o1gN& , such that Pfdðx&tjt ; xotjt Þ o eg 41 $ d for all t 4te and all sequences fMt gN& that satisfy the inequalities Mt 4 Mt;e;d for t ¼ 0; 1; 2; . . .. Notation Pfdðx&tjt ; xotjt Þ o eg is shorthand for PfA& g, where

ðMt Þ > o . 2 Rnx ' ( ( ( ' Rnx : dðxðkÞ A& 9f½xð1Þ t ; . . . ; xt t ; xtjt Þ o e

for some k 2 f1; . . . ; Mt gg:

ð23Þ

The proof of Proposition 1 is given in Appendix A and a detailed construction of PfA& g is provided in Appendix B. Formally, Proposition 1 states that it is possible to design SMCM algorithms that generate, with high probability, particles arbitrarily close to the minimizers xotjt whenever t, rt and Mt are sufficiently large. Indeed, we can prove that Pfdðx&tjt ; xotjt Þ o eg-1 exponentially fast as Mt -1 (see Eq. (A.23) in Appendix A for details). Intuitively, the proof of Proposition 1 shows that the algorithm generates a sequence of particle sets, t Lt ¼ fxtðiÞ gM , t ¼ 0; 1; 2; . . ., that converges toward, and i¼1 eventually ‘‘capture’’, the minimizers xotjt with high probability. This illustrated in Fig. 1. The particle sets L1 ; . . . ; L4 are depicted with solid black lines, while the corresponding minimizers, xo1j1 ; . . . ; xo4j4 , are denoted by stars. At time t ¼ 2; 3; 4, the image of Lt$1 , at ðLt$1 Þ, is shaded and delimited by a dashed line. The Figure is ‘‘read’’ from left to right. At time t ¼ 1, the cloud of

ARTICLE IN PRESS J. Mı´guez / Signal Processing 90 (2010) 1609–1622

1614

Λ1 Λ2 a1(Λ1)

Λ3 o * x1/1

x

Λ4

a2(Λ2)

o * 2/2

o * 3/3

x

a3(Λ3)

xo * 4/4

time (t=1)

(t=2)

(t=3)

(t=4)

Fig. 1. Illustrative depiction of the convergence of a SMCM algorithm according to Proposition 1. See Appendix A for a formal argument.

particles, L1 , is far away from the minimizer, xo1j1 , but the algorithm is constructed in such a way that, at time t ¼ 2:

+ the distance from x2j2 to the border of a1 ðL1 Þ is bounded,

+ the particle cloud L2 is a proper superset of the image a1 ðL1 Þ, i.e., a1 ðL1 Þ - L2 , and

+ the distance between xo2j2 and the boundary of L2 is

smaller than the distance from xo1j1 to the boundary of L1 .

At time t ¼ 3 the distance between the particle cloud L3 and xo3j3 is shorter than the distance between L2 and xo2j2 and, finally, at time t ¼ 4, xo4j4 lies ‘‘within’’ L4 . Remarkably, this kind of convergence can be achieved without a precise knowledge of the system dynamics (i.e., we do not use an explicit state-space model of the system). Instead, we need a sequence of functions fat gN& that comply with Assumptions 1–3. Such functions provide information on the dynamics of the minimizers, indeed, but this is very rough, since it reduces to a bound on the departure of xotjt from at$1 ðxot$1jt$1 Þ. Note that Fig. 1 is just a piece of artwork for the illustration of the intuitive ideas underlying the proof of Proposition 1. See Appendix A for the complete formal argument. Fig. 3 in Section 5.1 shows how the particle clouds (in a 1-dimensional space) actually converge to a sequence of exact minimizers for an actual SMCM algorithm. If we constrain the marginal cost function to depend only on the current state, i.e., ct ðxt$1:t ; yt Þ ¼ ct ðxt ; yt Þ, and to be continuous, we can obtain a straightforward corollary of Proposition 1. Let us consider two sequences, xot 2 S t ðyt Þ9arg minn ct ðxt ; yt Þ x t 2R

x

ð24Þ

Claim 1. Let us consider

+ a fixed sequence of marginal-cost minimizers fxot gN , + +

where xot 2 S t ðyt Þ, a sequence of functions at that satisfy Assumptions 1, 2 and 3 for the sequence fxot gN , and a sequence of propagation pdf’s of the form in Eq. (21).

If ct ðxt ; yt Þ is uniformly continuous at xt ¼ xot for all t 2 N and the finite radii frt gN are sufficiently large, then limt-1 limMt -1 dðx&t ; xot Þ ¼ 0 in probability. Specifically, the Claim states that for any arbitrarily small e; d 40 there exist te 2 N and a sequence fMt;e;d o1gN& such that Pfdðx&t ; xot Þ o eg 41 $ d for all t 4 te and all sequences fMt gt2N& satisfying Mt 4 Mt;e;d . A proof can be constructed that follows exactly the same steps as the proof of Proposition 1. The continuity of ct ðxt ; yt Þ is simply used to show that a sample sufficiently close to xot also has an ‘‘almost minimal’’ marginal cost. 5. Examples 5.1. A nonlinear one-dimensional system As a first example, we apply a SMCM algorithm to the minimization of the cost function ! "2 x2 Ct ðx0:t ; y1:t Þ ¼ ct ðxt$1:t ; yt Þ ¼ ð1 $ eÞ yt $ t þ eðxt $ at$1 ðxt$1 ÞÞ2 ; 20

ð26Þ

and t x&t 9S M t ðy t Þ9arg

for t ¼ 1; 2; . . ., that denote a sequence of marginal-cost minimizers and the approximate minimizers generated by the SMCM method, respectively. The following result states that x&t -xot in probability.

min

M

t xt 2fxðkÞ t gk ¼ 1

ct ðxt ; yt Þ

ð25Þ

t 2 N, where yt ; xt 2 R are the observations and the SoI, respectively, e Z0 is a non-negative weight and we set at ðxÞ ¼ x þ 8cosð1:2tÞ. When compared to the generic

ARTICLE IN PRESS J. Mı´guez / Signal Processing 90 (2010) 1609–1622

definition of (3), Ct has the form Hða; bÞ ¼ a, i.e., at each time step only the marginal cost is used to assess the quality of a particle. We have selected this specific structure because in this way we can obtain the minimizers xotjt analytically and, therefore, it is easier to illustrate the convergence of the SMCM algorithm in the way predicted either by Proposition 1 (when e 4 0) or Claim 1 (when e ¼ 0). To be specific, if we take e ¼ 0 then we readily find that ( 0 when yt r 0; o pffiffiffiffiffiffiffiffiffiffi ð27Þ xt ¼ 7 20yt when yt Z 0; hence there is a unique minimizer when the observation yt is negative and two optimal (equally valid) solutions when yt is positive. If we choose a positive weight, e 4 0, then the solution is similar, but it also involves xt$1 , namely xotjt ¼ xot

xot$1jt ¼ xotjt $ 8cosð1:2tÞ:

and

ð28Þ

This very simple system enables us to illustrate the convergence of both dðx&tjt ; xotjt Þ ¼ jx&tjt $ xotjt j and dðx&t ; xot Þ ¼ jx&t $ xot j. For the simulations, we generate the SoI and the observations from the nonlinear dynamic system x0 # N ðx0 ; 0; 1Þ; xt ¼ 0:5xt$1 þ yt ¼

ð29Þ

25xt$1 þ 8cosð1:2tÞ þ vt ; 1 þ x2t$1

x2t þ mt ; 20

ð30Þ ð31Þ

where xt 2 R is the system state at time t 2 N& , yt 2 R is the corresponding measurement and the noise processes vt # Nð0; s2m Þ and mt # Nð0; s2m Þ are zero-mean Gaussian with variances s2v ¼ 1 and s2m ¼ 2, respectively. Model (29)–(31) has been studied by several authors in the stochastic filtering literature [23,14,16]. In our case, it is just a mechanism to produce the signals in a way that departs explicitly from the selected sequence of functions at (compare the latter with the dynamic Eq. (30)). We have applied the simple SMCM algorithm below, with a fixed number of particles M ¼ Mt , to the sequential optimization of the cost function (26). (1) Initialization: Draw M i.i.d. samples from a Gaussian distribution centered at x ¼ 500, xðkÞ 0 # N ðx0 ; 500; 1Þ, k ¼ 1; . . . ; M, and take C0ðkÞ ¼ 0 for all k. (2) Recursive step: The selection pdf of choice at time t þ 1 is [21] ! " 1 $2 Bt þ 1 ðiÞp RtðiÞþ 1 $ minRðjÞ ; ð32Þ tþ1 þ M j where RtðiÞþ 1 ¼ ðyt þ 1 $ at ðxtðiÞ ÞÞ2 is a ‘‘risk’’ that we associate to the i-th particle. Therefore, indices are drawn as iðkÞ # Bt þ 1 ðiÞ, for k ¼ 1; . . . ; M and i 2 f1; . . . ; Mg. Taking the identity selection function, sðcÞ ¼ c, yields the set of selected particles ðkÞ ðkÞ ðkÞ fx~ ðkÞ ¼ xði Þ ; C~ ¼ C ði Þ gM . The propagation pdf t

t

t

t

k¼1

1615

of choice is ~ ðkÞ ~ ðkÞ ~ xðkÞ t þ 1 # r ðxt þ 1 jx t Þ ¼ TN ðxt þ 1 ; at ðx t Þ; 10Þ;

k ¼ 1; . . . ; M;

ð33Þ

2

where TN ðx; m; s Þ denotes a Gaussian density with mean m and variance s2 truncated outside the interval ðm $ 10s; m þ 10sÞ. The cost of the k-th particle is ðkÞ CtðkÞ ¼ ct ðx~ ðkÞ t ; xt þ 1 ; yt þ 1 Þ. Let us note that the SMCM algorithm described above strictly complies with the assumptions of Proposition 1 and Claim 1. Indeed,

+ the sequence of minimizers is deterministic (once +

+ +

the observations fyt gt for the simulation have been generated); the truncated Gaussian densities TN ðx; m; s2 Þ with support in Bðm; 10sÞ ¼ ðm $ 10s; m þ10sÞ comply with the definition of a density gðx; m; 10sÞ as specified in Section 4.2; the functions at ðxÞ ¼ x þ 8cosð1:2tÞ satisfy Assumptions 1, 2 and 3; and the radius rt ¼ 10s (for all t) has been found to be large enough to attain convergence in all simulations.

In particular, Assumptions 1 and 3 are straightforward to check, using the continuity of at , while, for Assumption 2, we note just that the upper bound A can be explicitly computed for each simulation (since the sequence of observations has a finite length). For our first simulation, we have independently generated 500 sequences, y0:T ðjÞ, with T ¼ 200 and j ¼ 1; 2; . . . ; 500, and computed the corresponding sequences of minimizers xotjt ðjÞ, for t ¼ 0; . . . ; T and j ¼ 1; . . . ; 500. Then, we have applied the SMCM algorithm with M ¼ 10; 000 particles to compute the sequence of approximate minimizers x&tjt ðjÞ, j ¼ 1; . . . ; 500, for every sequence of observations. Fig. 2 (left) shows the distance between the true and approximate minimizers as a function of time and P500 o 1 averaged over the 500 simulations, i.e., 500 j ¼ 1 jxtjt ðjÞ & $xtjt ðjÞj, as a solid line. It is clearly seen that the distance is very large during the initial time steps (notice the mismatch between the prior density of x0 and the initialization of the SMCM algorithm, which has been introduced with the purpose of illustrating the convergence of the algorithm with time explicitly). However, the proposed procedure converges steadily and, for t 4 40 time steps, x&tjt is already a good approximation of xotjt , with a mean distance of , 4 ' 10$4 . In the same plot, we also depict the averaged distance between xotjt and the minimum-cost particle. When e ¼ 0, the minimum-cost particle and x&tjt coincide (i.e., if e ¼ 0 then x&tjt ¼ x&t as defined in Eq. (25) and used in Claim 1), but this is not necessarily the case when e 40. Specifically, the dashed line in the figure shows this averaged distance when e ¼ 0:1. The reason for this mismatch is that the second term in the cost function (26) depends on the sample x&t$1jt , which is not updated by the SMCM algorithm at time t. This fact illustrates the

ARTICLE IN PRESS J. Mı´guez / Signal Processing 90 (2010) 1609–1622

1616

1000

1

100 0.1

10 1

0.01

0.1 0.01

0.001

0.001 0.0001

20

60

100

140

180

t (time)

0.0001

0

2000

4000

6000

8000

10000

M (no. of particles)

Fig. 2. Left: Average distance (a) between xotjt and x&tjt (solid line) and (b) between xotjt and the minimum-cost particle for a weight e ¼ 0:1 (dashed line), both of them as a function of time. Right: Average distance (a) between xotjt and x&tjt (solid line and hollow circles, ‘3’ Þ and (b) between xotjt and the minimum-cost particle for weight e ¼ 0:1 (dashed line and black squares, ‘’’ Þ as a function of the number of particles M.

particles and minimizers

400 350 300 250 200 150 100 50 0 10

20

30

40 50 t (time)

60

70

80

Fig. 3. Convergence of the particle sets. At time t ¼ 10, the set of samples M L10 ¼ fxðiÞ t gi ¼ 1 is far from the corresponding minimizers. The distance consistently reduces as time evolves and, after t ¼ 60, at least one of the minimizer at each time step lies among the samples in Lt .

need for a careful choice of the costs in practical problems where the aim is the estimation of the system state. Be aware, nevertheless, that this result is not contradictory with Proposition 1. Indeed, the simulation shows that x&tjt -xotjt both with e ¼ 0 and e 40, but in the latter case x&tjt is not the minimum-cost particle. While convergence with time is clearly observed in Fig. 2 (left), the accuracy of the steady-state distance between the approximate and true minimizers depends on the value of M. Fig. 2 (right) illustrates the convergence of the approximate minimizers with the number of particles, M. It shows, with a solid line connecting hollow circles, the distance between xotjt and x&tjt averaged over the last 100 time steps of 200 independent simulations. It can be seen how the distance decreases steadily with M. The average distance between the minimum-cost particle, for e ¼ 0:1, and the true minimizers is depicted with a dashed line connecting black squares. This distance also decreases with M, but apparently converges to a positive ‘‘floor’’

value of 0.1. Again, this is due to the second term in the cost function (26) and the fact that the sample x&t$1jt is inherited from time t $ 1 and not updated when yt is received. We conclude this example with a plot that shows how the particle sets Lt ¼ fxtðiÞ gM i ¼ 1 eventually ‘‘capture’’ the minimizers xotjt as time evolves. Only for this simulation, we set M ¼ 500. Fig. 3 is a numerical counterpart of the (merely illustrative) diagram in Fig. 1. It shows the minimizers xotjt , represented by dark-colored ‘*’ signs, for t ¼ 10; 20; . . . ; 80, together with the particle sets (lightcolored ‘+’ signs). Initially, the samples xtðiÞ are far away from the xotjt but, as time evolves, this distance is reduced. At times t ¼ 60; 70 and 80, at least one minimizer lies among the particles in Lt . 5.2. A chaotic CO2 laser model As a second example, we consider a more challenging system that describes the dynamics of a CO2 laser with modulated losses in a chaotic regime [18]. A deterministic model of the laser dynamics is given by the set of differential equations x_ 1 ðsÞ ¼ Q0 x1 ðsÞx2 ðsÞ $ Q1 ðsÞx1 ðsÞ; x_ 2 ðsÞ ¼ $ g1 x2 ðsÞ $ 2Q0 x1 ðsÞx2 ðsÞ þGx3 ðsÞ þ x4 ðsÞ þ P; x_ 3 ðsÞ ¼ $ g1 x3 ðsÞ þ Gx2 ðsÞ þx5 ðsÞ þP; x_ 4 ðsÞ ¼ $ g2 x4 ðsÞ þ Zx2 ðsÞ þ Gx5 ðsÞ þ ZP; x_ 5 ðsÞ ¼ $ g2 x5 ðsÞ þ Zx3 ðsÞ þ Gx4 ðsÞ þ ZP;

ð34Þ

where s 2 R denotes continuous time, x_ i ðsÞ ¼ dxi =ds denotes the time derivative of the i-th state variable, Q1 ðsÞ ¼ Q0 ½1 þ a sin2 ðB0 þ FðsÞÞ. and FðsÞ ¼ H sinð2pfsÞ is an external sinusoidal forcing signal with frequency f (Hz). In the above equations, x1 ðsÞ represents the laser output intensity, x2 ðsÞ is the population inversion between the two resonant levels, and x3 ðsÞ, x4 ðsÞ and x5 ðsÞ account for molecular exchanges between the two levels resonant with the radiation field and the other rotational levels of

ARTICLE IN PRESS J. Mı´guez / Signal Processing 90 (2010) 1609–1622

Table 1 Fixed parameter values for the CO2 laser model in chaotic regime.

a B0

H

4 0.403 0.060283

g1

g2

f

Q0

G

P

Z

1 7

29.95 10.0643 1.0643 0.05 0.0263 10

1617

mn # Nð0; s2m Þ is a zero-mean Gaussian perturbation with variance s2m ¼ 10$9 , z1;n 9x1;10n is the laser intensity variable and n 2 N. If we let zn ¼ ½z1;n ; . . . ; z5;n .> (with the obvious notation zi;n ¼ xi;10n ), then we can rewrite the statespace model used for the simulation of the CO2 laser as ðstate eq:Þ zn ¼ g10 ðzn$1 ; v10ðn$1Þ þ 1:10n Þ;

the same vibrational band. See [18] for physical details, including the meaning of the fixed parameters Q0 , a, B0 , H, g1 , g2 , G, Z and P. The values selected for the simulation ensure chaotic behavior and are shown in Table 1. We address the problem of tracking the latent variables fxi ðsÞg5i ¼ 1 in (34). In order to simulate the evolution of the system, we transform the continuoustime deterministic model (34) into a discrete-time dynamic random system. Time discretization is most easily carried out by Euler’s method and we assume a white Gaussian perturbation of the latent variables, i.e., x1;t ¼ x1;t$1 þ

Q0 x1;t$1 x2;t$1 $ Q1;t x1;t$1 ; T $1

$ðg1 þ 2Q0 x1;t$1 Þx2;t$1 þGx3;t$1 þx4;t$1 þP x2;t ¼ x2;t$1 þ þv2;t ; T $1

x3;t ¼ x3;t$1 þ

$g1 x3;t$1 þ Gx2;t$1 þ x5;t$1 þ P þ v3;t ; T $1

x4;t ¼ x4;t$1 þ

$g2 x4;t$1 þ Zx2;t$1 þ Gx5;t$1 þZP þ v4;t ; T $1

x5;t ¼ x5;t$1 þ

$g2 x5;t$1 þ Zx3;t$1 þ Gx4;t$1 þZP þ v5;t ; T $1

t 2 N, where T ¼ 5 ' 10$3 time units is the time step used to approximate the derivatives by differences, xi;t ¼ xi ðs ¼ tTÞ and vi;t # Nð0; s2v Þ, Q1;t ¼ Q1 ðs ¼ tTÞ, i ¼ 2; 3; 4; 5, are Gaussian noise terms with zero mean and variance s2v ¼ 10$6 . Note that x1;t is deterministic given fxi;t$1 g5i ¼ 1 (otherwise the features of the dynamics are changed abruptly compared to the deterministic continuous-time case). The state of system (35) consists of xt ¼ ½x1;t ; . . . ; x5;t .> 2 R5 and we start the simulation with the state value x0 ¼ ½8:0221 ' 10$03 ; 1:5817; 1:6770; 1:6384 ' 10; 1:6722 ' 10.> :

ð36Þ For convenience, we represent model (35) as a vector function g, i.e., xt ¼ gðxt$1 ; vt Þ, where vt ¼ ½0; v2;t ; . . . ; v5;t .> is the noise vector. Furthermore, we use the recursive notation gk ðxt$1 ; vt:t þ k$1 Þ9gðgk$1 ðxt$1 ; vt:t þ k$2 Þ; vt þ k$1 Þ;

ðobservation eq:Þ yn ¼ z1;n þmn :

ð39Þ

Ca;n ðz0:n ; y1:n Þ ¼ ca;n ðzn ; yn Þ ¼ jyn $ z1;n j;

ð40Þ

Cb;n ðz0:n ; y1:n Þ ¼ 12Cb;n$1 ðz0:n$1 ; y1:n$1 Þ þca;n ðzn ; yn Þ;

ð41Þ

8$ $2 > < $$1 $ z1;n $$ yn $ Cc;n ðz1:n ; y1:n Þ ¼ cc;n ðzn ; yn Þ ¼ $ > :1

ð42Þ

In this example we study the performance of the SMCM algorithms built from four different cost functions

8 > :1 $4

ð35Þ

ð38Þ

if yn Z U; if yn o U; ð43Þ

where n ¼ 1; 2; . . . and U ¼ 1:6 ' 10 is a threshold value. Note that functions Ca;n and Cc;n do not depend on the history of the SoI but only on the last value, i.e., they are of the form Hða; bÞ ¼ a, while functions Cb;n and Cd;n are constructed from the same marginal costs as Ca;n and Cc;n , respectively, but they depend on the past values of the SoI. The reason for the use of a threshold value, U, in the definition of Cc;n and Cd;n above is that the laser intensity can take very small values that are completely masked by the observational noise. In this case, the observation is not useful for the algorithm. The generic SMCM algorithm with fixed M ¼ Mt ¼ 500 is constructed as follows. (1) Initialization: Produce a sample of size M by ðkÞ setting zðkÞ 1;0 ¼ x1;0 and drawing zi;0 # Uðxi;0 $ 1; xi;0 þ 1Þ, i ¼ 2; . . . ; 5, k ¼ 1; . . . ; M. Assign initial costs, C0ðkÞ for k ¼ 1; 2; . . . ; M. This is different for the different ðkÞ ðkÞ ðkÞ ðkÞ ¼ Cb;0 ¼ 0 and Cc;0 ¼ Cd;0 ¼ 1. functions, namely, Ca;0 ðkÞ ðkÞ M (2) Recursive step: Assume On ¼ fzn ; Cn gk ¼ 1 is available. The selection pmf at time n þ 1 is ðkÞ $2 Bn þ 1 ðiÞpðRðiÞ ; n þ 1 $ minRn þ 1 þ 1=MÞ k

ð44Þ

10 ðkÞ where RðkÞ n þ 1 ¼ jyn þ 1 $ ½g ðzn ; 0Þ.1 j and ½u.i denotes

the i-th element of vector u. Note that RðkÞ n þ 1 can be interpreted as a prediction of the cost of the k-th particle. Indices are drawn as iðkÞ # Bn þ 1 ðiÞ, for k ¼ 1; . . . ; M and i 2 f1; . . . ; Mg, and the selected particles are ðkÞ ðkÞ ðkÞ fz~ ðkÞ ¼ znði Þ ; C~ ¼ Cnði Þ gM . We propagate the partin

n

k¼1

cles using Eq. (38) to obtain zðkÞ n þ 1 , k ¼ 1; . . . ; M, and then compute the corresponding costs CnðkÞþ 1 .

Note that in this algorithm the composition of functions g10 plays the role of an in the general scheme.

ARTICLE IN PRESS J. Mı´guez / Signal Processing 90 (2010) 1609–1622

1618

1.75

observed true BF SMCM (Cc)

1 0.01

dynamic variable no. 2

laser output intensity (mV)

100

0.0001 1e-06 1e-08 1e-10

true BF SMCM (Cc)

1.7 1.65 1.6 1.55 1.5 1.45

0

10

20

30

40

50

0

10

time (ms)

20

30

40

50

time (ms)

Fig. 4. Left: The laser output intensity (z1;n ), the noisy observations (jyn j) and its estimates using the BF and the SMCM algorithm with cost function Cc;n . Right: The dynamic variable z2;n together with its estimates obtained using the BF and the SMCM algorithm with cost function Cc;n .

Table 2 Normalized mean square error for each dynamic variable using the BF and the SMCM algorithms. ENMSE 1 BF SMCM (a) SMCM (b) SMCM (c) SMCM (d)

ENMSE 2

ENMSE 3

ENMSE 4

ENMSE 5

1:4 ' 10$3

7:8 ' 10$6

4:1 ' 10$5

3:7 ' 10$6

3:6 ' 10$5

1:7 ' 10$3

7:9 ' 10$5

8:8 ' 10$5

3:7 ' 10$5

6:5 ' 10$5

$3

3:1 ' 10

$3

2:3 ' 10

1:9 ' 10$2

$5

8:0 ' 10

$5

3:3 ' 10

4:4 ' 10$5

$5

8:8 ' 10

$5

8:8 ' 10

9:5 ' 10$5

4:2 ' 10

$5

1:4 ' 10

$5

2:5 ' 10$5

6:4 ' 10$5 6:4 ' 10$5 7:1 ' 10$5

Linear observations: yn ¼ z1;n þ mn . Table 3 Normalized mean square error for each dynamic variable using the BF and the SMCM algorithms.

BF SMCM (a) SMCM (b) SMCM (c) SMCM (d)

ENMSE 1

ENMSE 2

ENMSE 3

ENMSE 4

ENMSE 5

5:1 ' 10$1

1:7 ' 10$4

5:6 ' 10$5

1:1 ' 10$4

5:0 ' 10$5

3:1 ' 10$3

8:1 ' 10$5

9:0 ' 10$5

4:1 ' 10$5

6:6 ' 10$5

1:2 ' 10$3 5:0 ' 10$3

$2

4:6 ' 10

8:5 ' 10$5 3:5 ' 10$5 $5

5:8 ' 10

9:4 ' 10$5 8:8 ' 10$5 $5

8:8 ' 10

4:5 ' 10$5 1:5 ' 10$5 3:4 ' 10

$5

7:1 ' 10$5 6:5 ' 10$5 6:5 ' 10$5

Nonlinear observations: yn ¼ maxf10$4 ; z1;n þ mn g.

The propagation density does not match the assumptions of Proposition 1, though, as it does not have a finite support. In order to assess the tracking performance of the proposed method we have also applied the standard bootstrap filter (BF) described in Section 3.3 to the statespace model given by Eqs. (38) and (39). Fig. 4 illustrates the dynamics of the laser model and the performance of the SMCM algorithm with cost function Cc;n and the BF. Fig. 4 (left) shows a typical realization of the output intensity variable, z1;n , the absolute value of the associated observations, jyn j, and the estimates from the two Monte Carlo algorithms. It can be seen that the Gaussian noise completely hides the variable z1;n for relatively long periods of time. Only some spikes can be actually observed with accuracy, but tracking can be achieved, nonetheless, both with the BF and the SMCM algorithm. Valid estimates can also be computed for the latent variables z2;n ; . . . ; z5;n . In

particular, Fig. 4 (right) shows the signal z2;n and its estimates. We have numerically approximated the normalized mean square error (NMSE) for each variable zi;n , i ¼ 1; . . . ; 5. Let z^ i;n ðkÞ denote the estimate of variable zi;n ðkÞ obtained in the k-th simulation trial. For the SMCM algorithms, z^ i;n ðkÞ is the i-th component of the particle with the least cost, while for the BF the estimate results from the weighted mean of the particles. If we assume that Ns independent simulations have been carried out, each one consisting of a sequence of Ny observations, then the error is approximated as ¼ ENMSE i

Ny Ns X jzi;n ðkÞ $ z^ i;n ðkÞj2 1 X : Ny Ns k ¼ 1 n ¼ 1 PNy jzi;l ðkÞj2 l¼1

ð45Þ

, i ¼ 1; . . . ; 5, computed Table 2 shows the values ENMSE i from Ns ¼ 200 independent simulations, for each algorithm. With M ¼ 500 particles, all methods

ARTICLE IN PRESS J. Mı´guez / Signal Processing 90 (2010) 1609–1622

guarantee an adequate tracking of all variables and attain a similar performance (the error in the estimation of z1;n is only significantly higher for the SMCM algorithm with cost function Cd;n ). From the perspective of the SMCM methodology, it is interesting to look into the relative performance when there is some discrepancy between the assumed model and the BF. In particular, we have simulated nonlinear observations of the form yn ¼ maxf10$4 ; z1;n þ mn g, with mn # Nð0; 10$9 Þ as before.2 We have applied the five algorithms without modification (in particular, the BF is still designed assuming linear observations) and evalu, through another 200 ated the average error values, ENMSE i independent simulation trials. Our results, shown in Table 3, illustrate the superior robustness of the SMCM algorithms, which attain error values very similar to than those shown in Table 2, whereas the performance of the BF degrades considerably for z1;n ; z2;n and z4;n .

1619

Acknowledgements This work has been partially supported by the Ministry of Science and Innovation of Spain (Project MONIN, Ref. TEC-2006-13514-C02-01/TCM, and program Consolider-Ingenio 2010, Project CSD2008-00010 COMONSENS) and the Autonomous Community of Madrid (Project PROMULTIDIS-CM, Ref. S-0505/TIC/0233).

Appendix A. Proof of Proposition 1 In this appendix we address the formal proof of Proposition 1. We start with two auxiliary results, that involve the introduction of some new notation, and then proceed with the main proof. The latter is based on a simple induction argument and a classical inequality for the sum of bounded random variables. A.1. Auxiliary results

6. Conclusions We have revisited the sequential Monte Carlo minimization method originally introduced in [21,22,19]. As a result, we have:

+ proposed an extended description of the methodology +

that allows to tackle the optimization of a larger class of cost functions and analyzed the asymptotic convergence of SMCM algorithms both with time and the number of particles.

Our results show that properly designed SMCM algorithms can generate, with high probability, arbitrarily accurate estimates of a sequence of cost minimizers. The proof, shown in Appendix A, is lengthy but based on elementary results and a simple induction argument (an intuitive sketch of the argument has been provided in Section 4.3). One remarkable consequence of the analysis is that the convergence of SMCM algorithms depends very weakly on the dynamics of the signal of interest. Indeed, virtually any function that provides a rough approximation of the evolution of the cost minimizer at time t into the minimizer at time t þ 1 enables the design of a convergent algorithm. These results have been illustrated by way of a simple example involving a nonlinear one dimensional system. Finally, we have also considered the problem of tracking the latent physical variables in the dynamic model of a CO2 laser. In this case, we have compared the tracking performance of four SMCM algorithms (built around different cost functions) and a standard particle filter. We have shown that, at least for this example, the SMCM algorithm is more robust to variations in the form of the observations than the particle filter. 2 This model makes practical sense because experimental observations also show such a ‘‘floor’’ effect, i.e., very small values of the intensity variable x1;t are not practically observable [18].

We first use the sequence of functions fat gN& , to establish a connection between the dynamics of the minimizers, fxotjt gN& , and the dynamics of an adequately defined sequence of sets. In particular, let fX t gN& be a sequence of sets in Bnx and let us assume that they satisfy the following assumptions. t Assumption 4. For all t 2 N& , X t ¼ [ni ¼ 1 Oi;t , where nt is a (finite) natural number and O1;t ; . . . ; On;t are bounded open cells in Rnx .

Let A; B be proper subsets of Rnx and let dðx; x0 Þ denote a proper distance between points x; x0 2 Rnx . We define the distance between A and B as dðA; BÞ9 inf dðy; xÞ; x2A;y2B

ðA:1Þ

2Ag. and the complement of A as A9fy 2 Rnx : y= Assumption 5. For all t 2 N we assume that dðX t ; at$1 ðX t$1 ÞÞ 4 0. Equivalently, at$1 ðX t$1 Þ is a proper subset (not equal to) X t . The class of sets that complies with Assumption 4 is very large. In particular, note that for any measurable set E 2 Bnx , such that ‘ðEÞ o1, and for any e 40, there exists a sequence of bounded open cells fOi gni¼ 1 , with n o 1, such that the measure of the symmetric difference between the set E and its approximation E n ¼ [ni¼ 1 Oi is less than e [2, Theorem 15.10], i.e., ‘ðE DE n Þ ¼ ‘ððE $ E n Þ[ ðE n $ EÞÞ o e. Next, we introduce some notation needed to express the properties of the sequence fX t gN& relative to the sequence of minimizers fxotjt gN& . The distance between a point x 2 Rnx and a set X 2 Rnx is defined as dðx; X Þ9inf x0 2X dðx; x0 Þ. In our case, it is of interest to upper bound the separation between a minimizer xotþ 1jt þ 1 and the set X t þ 1 , given the previous distance dðX t ; xotjt Þ. We start with a bound on the departure between two arbitrary points x; x0 2 Rnx when they are propagated

ARTICLE IN PRESS J. Mı´guez / Signal Processing 90 (2010) 1609–1622

1620

through the function at , namely 0 Dsup t þ 1 9 sup dðat ðxÞ; at ðx ÞÞ; x;x0 2Rnx

subject to dðx; x0 Þ r dðX t ; xotjt Þ;

ðA:2Þ while the corresponding upper bound for the distance increment is sup Dsup t þ 1 9Dt þ 1

$

dðX t ; xotjt Þ:

ðA:3Þ

Note that, because of Assumption 3, if dðX t ; xotjt Þ o 1 then the supremum in (A.2) is bounded, i.e., Dsup t þ 1 o 1. From

Assumption 5, we can also define the sequence of positive real numbers dout t þ 1 9dðX t þ 1 ; at ðX t ÞÞ 4 0;

t 2 N& ;

ðA:4Þ

and, taking together (A.3), (A.4) and Assumption 2, we introduce the sequence sup Kt þ 1 9dout t þ 1 $ At $ Dt þ 1 ;

t 2 N& ;

ðA:5Þ

with maximum lower bound denoted as K9inf t2N Kt . It can be shown (see the proof of Lemma 1 below) that Kt þ 1 is the worst-case reduction of the distance dðX t ; xotjt Þ that we achieve at time t þ 1, i.e., Kt þ 1 rdðX t ; xotjt Þ$ dðX t þ 1 ; xotþ 1jt þ 1 Þ. Lemma 1. If the following assumptions hold for some sequence fxotjt gN : (a) K Z0. (b) There exists L 2 ð0; þ 1Þ such that dðX 1 ; xo1j1 Þ rL.

(c) Given L, for all e 4 0 there exists te 2 N such that, for all P t 4te , tn ¼ 2 Kn 4 L $ e.

Then, for all t 4 te , dðX t ; xotjt Þ o e. Moreover, if L4 dðX 1 ; xo1j1 Þ þ e then dðX t ; xotjt Þ ¼ 0 for all t 4 te . Proof. From Assumption 5 in the definition of sequence fX t gN& , we deduce an initial upper bound of dðX t ; xotjt Þ, t Z 2,

dðX t ; xotjt Þ rmaxf0; dðat$1 ðX t$1 Þ; xotjt Þ $ dout t g;

ðA:6Þ

which using the triangular inequality, dðat$1 ðX t$1 Þ; xotjt Þ r dðat$1 ðX t$1 Þ; at$1 ðxot$1jt$1 ÞÞ þ dðat$1 ðxot$1jt$1 Þ; xotjt Þ;

ðA:7Þ and Assumption 2 yields dðX t ; xotjt Þ r maxf0; dðat$1 ðX t$1 Þ; at$1 ðxot$1jt$1 ÞÞ þ At$1 $ dout t g:

From the definition of

Dsup t

and

dðat$1 ðX t$1 Þ; at$1 ðxot$1jt$1 ÞÞ r Dsup t

Dsup t ,

¼

ðA:8Þ

we have

o Dsup t þ dðX t$1 ; xt$1jt$1 Þ;

ðA:9Þ

hence, substituting (A.9) into (A.8), we arrive at the inequality out dðX t ; xotjt Þ rmaxf0; dðX t$1 ; xot$1jt$1 Þ þ Dsup t þ At$1 $ dt g

ðA:10Þ

and, since we have previously defined At$1 , we readily obtain dðX t ; xotjt Þ rmaxf0; dðX t$1 ; xot$1jt$1 Þ $ Kt g:

Kt ¼ dout t

$ Dsup t $ ðA:11Þ

The assumption Kt ZK Z 0, for t Z 2, enables us to recursively apply (A.11) to find the relationship ( ) t X o o ðA:12Þ Kn ; dðX t ; xtjt Þ rmax 0; dðX 1 ; x1j1 Þ $ n¼2

valid for any t Z 2. Since we have assumed that dðX 1 ; xo1j1 Þ rL and also that for all e 40 there exists te 2 Pt N such that n ¼ 2 Kn 4 L $ e whenever t 4 te , then it follows from (A.12) that dðX t ; xotjt Þ o e whenever t 4 te . Pt o Moreover, if, additionally, n ¼ 2 Kn 4 L $ e 4 dðX 1 ; x1j1 Þ we trivially arrive at dðX t ; xotjt Þ ¼ 0. & Intuitively, Lemma 1 states that if the distance between the first minimizer xo1j1 and the set X 1 is finite, and the P P sum tn ¼ 2 Kn ¼ tn ¼ 2 jKn j is large enough, then the sets X t will eventually get arbitrarily close, or even include, the corresponding minimizers xotjt . The following corollary of Lemma 1 provides a more direct way of applying this result. Lemma 2. If dðX 1 ; xo1j1 Þ rJ o1 and K 40 for some sequence fxotjt gN , then there exists tJ 2 N such that, for all t 4 tJ , dðX t ; xotjt Þ ¼ 0. Pt Proof. Since n ¼ 2 Kn Z ðt $ 1ÞK and K 40, then for any L ¼ J þ 2e (where e 40) there exists tJ 2 N such that ðt $ 1ÞK 4 L ¼ J þ 2e 4 dðX 1 ; xo1j1 Þ þ e, for all t 4tJ . Hence, we can apply Lemma 1 to obtain dðX t ; xotjt Þ ¼ 0 for all t 4 tJ . & P Therefore, if the sum tn ¼ 2 Kn is not bounded, as t-1, the sets fX t gN will eventually include the minimizers xotjt . P We note that limt-1 tn ¼ 2 Kn ¼ 1 does not necessarily imply that limt-1 ‘ðX t Þ ¼ 1. This can be avoided by a proper choice of ‘‘contractive’’ functions at such that ‘ðat ðX t ÞÞ o‘ðX t Þ. A.2. Main proof The proof of Proposition 1 consists of two parts: (1) We show how a SMCM algorithm can generate a sequence of sets fX t 2 Bnx gN& that complies with the assumptions of Lemma 2 and such that we can draw samples arbitrarily close to any point in any set X t . (2) We apply Lemma 2 to show that, for sufficiently large t 2 N, dðX t ; xotjt Þ can be made as small as we wish and, therefore, we can draw samples from X t arbitrarily close to xotjt with high probability. A.2.1. Construction of the sequence fX t gN& Let e 40 and d 40 be arbitrarily small but fixed positive real numbers and let X 0 ¼ [ni 0¼ 1 Oi;0 , with arbitrary n0 o 1 and bounded open cells Oi;0 2 Bnx . Since ‘ðOi;0 Þ o 1 for all 1 r ir n0 and, as a consequence, ‘ðX 0 Þ o 1, it is possible to build a proper pdf uniformly continuous and strictly positive in X 0 , denoted r0 . If we draw a set of M0 independent and identically distributed 0 , we can invoke samples from r0 , denoted O0 ¼ fx0ðiÞ gM i¼1 the weak law of large numbers (see, e.g., [13, Chapter 7]) to claim that, for any x 2 X 0 , there exists M0;e;d 2 N such

ARTICLE IN PRESS J. Mı´guez / Signal Processing 90 (2010) 1609–1622

that, for all M0 4 M0;e;d , n o e for some k 2 f1; . . . ; M0 g 4 1 $ d: P dðx; xðkÞ 0 Þo 2

The expected value of Zt is non-zero, namely ðA:13Þ

Note that the event ‘‘dðx; xðkÞ 0 Þ o e=2 for some k 2 f1; . . . ; M0 g’’ implies the event ‘‘mink2f1;...;M0 g dðx; xðkÞ 0 Þo e=2’’ . Next, we build the sequence fX t gN recursively. Our induction hypothesis is that X t$1 is a finite Borel set that ðiÞ Mt$1 gi ¼ 1 drawn from a pdf contains all discrete points fxt$1 rt$1 , where Mt$1 is a natural number (hence, Mt$1 o 1). This is actually weaker than Assumption 4 but is sufficient for our argument. Let us define 9 sup dðx; xðkÞ t$1 Þ: x2X t$1 ;k2f1;...;Mt$1 g

Z

ðA:14Þ

Since ‘ðX t$1 Þ o 1 by assumption, Z o1 and the set Mt$1 fxðkÞ t$1 gk ¼ 1 is a finite Z$net of X t$1 (recall that Mt$1 can be arbitrarily large but finite). As a consequence, it turns out that at$1 ðX t$1 Þ D

M t$1 [

k¼1

at$1 ðBðxðkÞ t$1 ; ZÞÞ;

ðA:15Þ

where Bðx; ZÞ9fx0 2 Rnx : dðx; x0 Þ o Zg 2 Bnx denotes an open ball centered at x with radius Z 40. Moreover, if we construct X t9

M t$1 [

k¼1

Bðat$1 ðxðkÞ t$1 Þ; rt Þ;

ðA:16Þ

with sufficiently large rt , then we can write at$1 ðX t$1 Þ D

M t$1 [

k¼1

at$1 ðBðxðkÞ t$1 ; ZÞÞ - X t :

ðA:17Þ

sup

Z

sup

x;x0 2X t$1 ;

0

dðat$1 ðxÞ; at$1 ðx ÞÞ

x;x0 2Rnx ;dðx;x0 Þ o 2Z

Mt 1 X E½zðkÞ . ¼ em ; Mt k ¼ 1 t

ðA:21Þ

where em is independent of Mt . For conciseness, let Q ðe; Mt Þ denote the event 0 ‘‘dðxðkÞ t ; x Þ o e=2 for some k 2 f1; . . . ; Mt g’’ . Obviously, Q ðe; Mt Þ is equivalent to the event ‘‘Zt 40’’ and we can write PfQ ðe; Mt Þg ¼ PfZt 40g 4 PfjZt $ em jo xg

for any small constant 0 o x o em . Therefore,

ðA:22Þ 2

PfQ ðe; Mt Þg 4 1 $ PfjZt $ em jZ xg Z1 $ 2expf$2Mt x g;

ðA:23Þ

where we have applied one of Hoeffding’s inequalities [15, Theorem 1] to obtain the bound. Since limMt -1 2 expf$2Mt x g ¼ 0, for any d 4 0 we can find Mt;e;d such that, whenever Mt 4Mt;e;d , the inequality 1 $ 2expf$ 2 2Mt x g 41 $ d is satisfied. A.2.2. Application of Lemma 2 In order to apply Lemma 2 to the sequence fX t gN& we ¼ dðX t ; at$1 ðX t$1 ÞÞ be large need to guarantee that dout t enough. Specifically, the inequality sup dout t 4 At$1 þ Dt

ðA:24Þ

must be satisfied for all t. We recall that A 4 At$1 ¼ is defined dðxotjt ; at$1 ðxot$1jt$1 ÞÞ by Assumption 2 and Dsup t by (A.3) and (A.2). In our case, from the definition of dout t , must be chosen to (A.16) and (A.18) we obtain that dout t satisfy the inequality sup x;x0 2Rnx ;

dðat$1 ðxÞ; at$1 ðx0 ÞÞ:

dðx;x0 Þ o 2Z

ðA:25Þ

Since it is always possible to choose rt such that

dðat$1 ðxÞ; at$1 ðx0 ÞÞ

ðA:18Þ

dðx;x0 Þ o 2Z

(note that the suprema are finite because of Assumption 3). Given Mt$1 o1, (A.16) and (A.18) ensure that (A.17) holds and, therefore, X t complies with Assumptions 4 and 5. We naturally associate the propagation pdf rt defined in Eq. (21) to the set X t . Recall that, since we have assumed Bt ðkÞ 40 for all k 2 f1; . . . ; Mt$1 g (which is possible because Mt$1 o 1), it follows that rt is uniformly continuous and strictly positive in X t . Let xðkÞ t # rt ðxt Þ, k ¼ 1; . . . ; Mt . We need to prove that we can draw samples arbitrarily close to an arbitrary point x0 2 X t . First, we note that, from the construction of X t and the choice of rt we obtain n % e &o 0 9em 4 0 ðA:19Þ P xðkÞ t 2 B x; 2 (note that ‘ðX t Þ o 1, since Mt$1 o 1 and rt o 1). We also introduce the auxiliary random variables 8 % & < 1 if xðkÞ 2 B x0 ; e ; t ðkÞ 2 zt 9 ðA:20Þ : 0 otherwise; and the (random) sample mean Zt 9ð1=Mt Þ

E½Zt . ¼

dout t Z rt $

In particular, (A.17) holds if we take rt 4

1621

PMt

ðkÞ k ¼ 1 zt r 1.

rt 4

sup

dðx;x0 Þ o 2Z

x;x0 2Rnx ;

dðat$1 ðxÞ; at$1 ðx0 ÞÞ þ Dsup t þ At$1 ; ðA:26Þ

then it is also always possible to make the inequality (A.24) hold for all t and, as a consequence, K ¼ inf t2N dout t $ At$1 $ Dsup t 40. Therefore, we can apply Lemma 2 to prove the existence of te 2 N such that, for all t 4te , dðX t ; xotjt Þ o

e

2

:

ðA:27Þ

Also, from the definition of the sequences fX t gN& and frt gN& , there exists Mt;e;d such that, for all Mt 4 Mt;e;d and all x 2 X t , PfdðxðkÞ t ; xÞ o

e

2

for some k 2 f1; . . . ; Mt gg 4 1 $ d:

ðA:28Þ

Taking together (A.27) and (A.28), we arrive at n o e o P dðxðkÞ ¼ e; for some k 2 f1; . . . ; Mt g 41 $ d: t ; xtjt Þ o2 2 ðA:29Þ

o Since the event ‘‘dðxðkÞ t ; xtjt Þ o e for some k 2 f1; . . . ; Mt g’’ o implies the event ‘‘mink2f1;...;Mt g dðxðkÞ t ; xtjt Þ o e’’ , the proof is complete.

ARTICLE IN PRESS J. Mı´guez / Signal Processing 90 (2010) 1609–1622

1622

Appendix B. The probability measure for the random vector x&t Recall that the sequence of observations y1:t is fixed, hence the randomness in the SMCM method comes exclusively from the sampling process using the sequence of pdf’s rt defined in Eq. (21). Let us put together the samples at time t in order to >

>

tÞ ; . . . ; xðM .> xst ¼ ½xð1Þ t t Q ðiÞ M 2R . At time t ¼ 0, xs0 # rs0 ðxs0 Þ ¼ i ¼0 1 r0 ðx0 Þ while, Q t ðiÞ at time t and given xst$1 , xst # rst ðxst Þ ¼ M i ¼ 1 rt ðxt Þ QMt PMt$1 ðiÞ ðkÞ s ¼ i ¼ 1 k ¼ 1 r~ t ðxt jxt$1 ÞBt ðkÞ. The process fxt gN is firstorder Markovian and, as a consequence, we can write down the joint pdf for the sequence xs0:t , namely Q rs0:t ðxs0:t Þ ¼ tn ¼ 0 rsn ðxsn Þ. Therefore, we construct the associated joint probability measure as Z P0:t rs0:t ðxs0:t Þdxs0 ; . . . ; dxst ; ðB:1Þ rs fAg ¼

construct

the

random

vector

nx Mt

A

where A 2 Bnx ðM0 þ ((( þ Mt Þ , and the marginal measure of xst is ( Z 'Z rs0:t ðxs0:t Þdxs0 ; . . . ; dxst$1 dxst : Ptrs fDg ¼ D

Rnx ðM0 þ ((( þ Mt$1 Þ

‘‘dðx&tjt ; xotjt Þ o

ðB:2Þ

e’’ in Proposition 1 is the set A& 2 The event nx Mt defined by Eq. (23), hence, Pfdðx&tjt ; xotjt Þ o eg ¼ Ptrs fA& g. B References [1] A. Bain, D. Crisan, Fundamentals of Stochastic Filtering, Springer, Berlin, 2008. [2] R.G. Bartle, The Elements of Integration and Lebesgue Measure. Wiley Classics Library, Wiley & sons, New York, 1995. [3] M.F. Bugallo, S. Xu, P.M. Djuric´, Performance comparison of EKF and particle filtering methods for maneuvering targets, Digital Signal Processing 17 (October 2007) 774–786. [4] J. Carpenter, P. Clifford, P. Fearnhead, Improved particle filter for nonlinear problems, IEE Proceedings—Radar, Sonar and Navigation 146 (1) (February 1999) 2–7. [5] D. Crisan, Particle filters—a theoretical perspective, in: A. Doucet, N. de Freitas, N. Gordon (Eds.), Sequential Monte Carlo Methods in Practice, Springer, Berlin, 2001, pp. 17–42 (Chapter 2). [6] D. Crisan, A. Doucet, A survey of convergence results on particle filtering, IEEE Transactions on Signal Processing 50 (3) (March 2002) 736–746.

[7] M.H. DeGroot, M.J. Schervish, Probability and Statistics, third ed., Addison-Wesley, New York, 2002. [8] P.M. Djuric´, J.H. Kotecha, J. Zhang, Y. Huang, T. Ghirmai, M.F. Bugallo, J. Mı´guez, Particle filtering, IEEE Signal Processing Magazine 20 (5) (September 2003) 19–38. [9] P.M. Djuric´, M. Vemula, M.F. Bugallo, Target tracking by particle filtering in binary sensor networks, IEEE Transactions on Signal Processing 56 (6) (June 2008) 2229–2238. [10] R. Douc, O. Cappe´, E. Moulines, Comparison of resampling schemes for particle filtering, in: Proceedings of the Fourth International Symposium on Image and Signal Processing and Analysis, September 2005, pp. 64–69. [11] A. Doucet, N. de Freitas, N. Gordon (Eds.), Sequential Monte Carlo Methods in Practice, Springer, New York, USA, 2001. [12] A. Doucet, S. Godsill, C. Andrieu, On sequential Monte Carlo sampling methods for Bayesian filtering, Statistics and Computing 10 (3) (2000) 197–208. [13] B.V. Gnedenko, The Theory of Probability, sixth ed., Gordon and Breach Science Publishers, 1997. [14] N. Gordon, D. Salmond, A.F.M. Smith, Novel approach to nonlinear and non-Gaussian Bayesian state estimation, IEE Proceedings— F 140 (2) (1993) 107–113. [15] W. Hoeffding, Probability inequalities of sums of bounded random variables, Journal of the American Statistical Association 58 (301) (March 1963) 13–30. [16] G. Kitagawa, S. Sato, Monte Carlo smoothing and self-organising state-space model, in: A. Doucet, N. de Freitas, N. Gordon (Eds.), Sequential Monte Carlo Methods in Practice, Springer, Berlin, 2001, pp. 177–196 (Chapter 9). ¨ [17] H.R. Kunsch, Recursive Monte Carlo filters: algorithms and theoretical bounds, The Annals of Statistics 33 (5) (2005) 1983–2021. ˜o, E. Allaria, M.A.F. Sanjua´n, R. Meucci, F.T. Arecchi, Coupling [18] I.P. Marin scheme for complete synchronization of periodically forced chaotic CO2 lasers, Physical Review E 70 (036208) (2004). [19] J. Mı´guez, Analysis of selection methods for cost-reference particle filtering with applications to maneuvering target tracking and dynamic optimization, Digital Signal Processing 17 (2007) 787–807. [20] J. Mı´guez, A. Arte´s-Rodrı´guez, Particle filtering algorithms for tracking a maneuvering target using a network of wireless dynamic sensors, EURASIP Journal on Advances in Signal Processing 2006 (2006) 1–16 (Article ID 83042, doi:10.1155/ASP/2006/83042). [21] J. Mı´guez, M.F. Bugallo, P.M. Djuric´, A new class of particle filters for random dynamical systems with unknown statistics, EURASIP Journal on Advances in Signal Processing 2004 (15) (November 2004) 2278–2294. [22] J. Mı´guez, M.F. Bugallo, P.M. Djuric´. Erratum to ‘‘a new class of particle filters for random dynamic systems with unknown statistics’’, EURASIP Journal on Advances in Signal Processing 2006 (2006) 1–16 (Article ID 78708, doi:10.1155/ASP/2006/78708). [23] M. Netto, L. Gimeno, M. Mendes, On the optimal and suboptimal nonlinear filtering problem for discrete-time systems, IEEE Transactions on Automatic Control 23 (6) (1978) 1062–1067. [24] B. Ristic, S. Arulampalam, N. Gordon, Beyond the Kalman Filter: Particle Filters for Tracking Applications, Artech House, Boston, 2004.