Keywords: monitoring, mechanical system, modal parameter, recursive maximum likelihood, estimation, tracking, particle lter, linear tangent particle lter surveillance, syst me m canique, param tre modal, maximum de vraisemblance r cursif, estimation, poursuite, ltre particulaire, ltre particulaire lin aire tangent

Easy (FLITE).

I R I S A I N S T I T U T D E R E C H E R C H E E N I N F O R M A T I Q U E E T S Y S T È M E S A L É A T O I R E S
1 Introduction

The purpose of this work is to obtain as much intuition as possible, through numerical experiments in a simple case where exact solutions are explicitly available, about the particle approximation of nite signed measures. A prototypical example of a nite signed measure is the derivative, w.r.t. a parameter of the model, of some probability distributions related with a hidden Markov chain. This includes prior, prediction, ltering probability distributions, etc. Two points of view are considered here, to feel the quality of the approximation, at least in a qualitative manner : (i) how accurate is the particle approximation of the nite signed measure, in view of an histogram representation of the weighted particle system ? (ii) considering the log likelihood function and the score function, how close is the approximate expression provided by the particle approximation to the exact expression ? These two questions seem closely related, however the numerical experiments presented in this work show that one of the two particle approximation schemes fails to satisfy the rst criteria (quality of the approximation of the nite signed measure), and that both schemes satisfy the second criteria (quality of the approximation of the statistics).

The initial motivation for this work was provided by an application to monitoring the integrity of structural and mechanical systems. Detecting and localizing damages for monitoring the integrity of structural and mechanical systems is a topic of growing interest, due to the aging of many engineering constructions and machines and to increased safety norms. Automatic global vibration based monitoring techniques turn out to be useful alternatives to visual inspections or local non destructive (e.g. ultrasonic) evaluations performed manually.

Health monitoring techniques based on processing vibration measurements basically handle two types of characteristics: the structural parameters (mass, sti ness, exibility, damping) and the modal parameters (modal frequencies, and associated damping values and mode shapes), see [START_REF] Hans G Nther Natke | Model Aided Diagnosis of Mechanical Systems : Fundamentals, Detection, Localization, and Assessment[END_REF]9,[START_REF] Mevel | Merging sensor data from multiple measurement setups for nonstationary subspace based modal analysis[END_REF]. A central question for monitoring is to compute changes in those characteristics and to assess their signi cance. For the frequencies, crucial issues are then: how to compute the changes, to assess that the changes are signi cant, to handle correlations among individual changes. A related issue is how to compare the changes in the frequencies obtained from experimental data with the sensitivity of modal parameters obtained from an analytical model. Furthermore, it has been widely acknowledged that changes in frequencies bear useful information for damage detection.

Our contribution in this work is to design a particle ltering method to track the modal parameters. Particle ltering techniques are a set of powerful and versatile simulation based methods to perform optimal state estimation in nonlinear non Gaussian state space models, and we consider here an approach combining particle ltering and gradient algorithm to perform recursive maximum likelihood parameter estimation and tracking. In the next section, the modeling issues are introduced and some key parameterizations are discussed. Section 3 details the particle approximation mechanisms. Section 4 is devoted to particle approximation of nite signed measures, with numerical experiments. In Section 5, the particle implementation of the recursive maximum likelihood (RML) algorithm is described.

Mechanical model 2.1 Dynamical model and structural parameters

It is assumed that the behavior of the mechanical system can be described by a stationary linear dynamical system, and that, in the frequency range of interest, the input forces can be modeled as a non-stationary white noise. This results in :

M Z(t) + C Ż(t) + KZ(t) = ν(t) Y (t) = LZ(t) (1)
where t denotes continuous time, M , C and K are the mass, damping and sti ness matrices respectively, the (high dimensional) vector Z collects the displacements of the degrees of freedom of the structure, the PI n 1604 external (non measured) force ν is modeled as a non stationary white noise with time varying covariance matrix Q ν (t), measurements are collected in the (often, low dimensional) vector Y , and the matrix L indicates which components of the state vector are actually measured, i.e. where the sensors are located. The modes or eigenfrequencies, denoted generically by µ, are solutions of :

det(µ 2 M + µ C + K) = 0 .
(2)

The frequency f and the damping coe cient d are recovered from a continuous eigenvalue µ through :

f = ℑ(µ) 2 π
and d = -ℜ(µ) .

Some comments are in order on parameterizations of interest for damage detection and localization. Since a local damage in the structure reduces the sti ness and increases the damping, many damage detection techniques have been proposed which monitor the sti ness matrix K. Monitoring its inverse K -1 , namely the exibility matrix, has proven more tractable and computationally feasible [START_REF] Hans G Nther Natke | Model Aided Diagnosis of Mechanical Systems : Fundamentals, Detection, Localization, and Assessment[END_REF][START_REF] Ewins | Modal testing : Theory, Practice and Application[END_REF]5]. In some cases, other structural parameterizations such as volumic mass and Young elasticity modulus may be preferable [START_REF] Farrar | Vibration based structural damage identi cation[END_REF][START_REF] Hans G Nther Natke | Model Aided Diagnosis of Mechanical Systems : Fundamentals, Detection, Localization, and Assessment[END_REF]. Also, several methods in the literature are based on a transmissibility matrix [START_REF] Mevel | Merging sensor data from multiple measurement setups for nonstationary subspace based modal analysis[END_REF][START_REF] Pandey | Damage detection in structures using changes in exibility[END_REF], which involves the processing of input output data. However, in the case of non measured input excitation, processing output only data is mandatory [START_REF] Guyader | A particle implementation of the recursive MLE for partially observed di usions[END_REF]3]. On the other hand, a reduced sti ness and an increased damping result in decreased natural frequencies. Thus, monitoring the modal parameters is relevant.

State space model and parameterization

Sampling model (1) at rate 1/∆ yields the discrete time model in state space form 12,[START_REF] Hermans | Modal testing and analysis of structures under operational conditions : industrial applications[END_REF] :

X k+1 = F X k + W k Y k = H X k (4)
where the state and the output are :

X k = Z(k ∆) Ż(k ∆)
and

Y k = Y (k ∆) , (5)
the state transition and observation matrices are :

F = exp(A ∆)
where

A = 0 I -M -1 K -M -1 C and H = L 0 . (6)
The measurement equation in (4) with H as in (6) implicitly assumes that the available sensors measure the (relative) displacements of the degrees of freedom themselves. The nature of the sensors used only in uences the observation matrix H. In (4), the unmeasured state noise {W k , k ≥ 0} is assumed to be a zero mean Gaussian white noise, with covariance matrix :

Q k △ = E[W k W * k ] = (k+1) ∆ k ∆ exp(A s) Q(s) exp(A * s) ds and
Q(s) = 0 0 0 M -1 Q ν (s) M - * .
The state X and the observed output Y have dimensions 2m and r respectively, with r (often much) smaller than 2m in practice.

Irisa

In view of the key assumption that it is on the other hand easy to simulate r.v.'s with probability distribution Q(x, dx ′ ), the idea is to approximate the predictor µ k|k-1 with the empirical probability distribution associated with an N sample, i.e.

µ k|k-1 ≈ µ N k|k-1 = 1 N N i=1 δ ξ i k|k-1 .
This approximation is completely characterized by the set {ξ i k|k-1 , i = 1, • • • , N } of particles, and the algorithm is completely described by the mechanism which builds

{ξ i k+1|k , i = 1, • • • , N } from {ξ i k|k-1 , i = 1, • • • , N }.
This mechanism is as follows : (i) the correction step is applied exactly to µ N k|k-1 , which results in

µ N k = Ψ k • µ N k|k-1 = N i=1 Ψ k (ξ i k|k-1 ) δ ξ i k|k-1 N j=1 Ψ k (ξ j k|k-1 ) = N i=1 ω i k δ ξ i k|k-1 ,
i.e. particles {ξ i k|k-1 , i = 1, • • • , N } are now weighted, with weights {ω i k , i = 1, • • • , N } which are more heavy for those particles which are more consistent with the current observation Y k , Irisa (ii) instead of trying to compute Q µ N k , the following particle approximation

µ N k+1|k = S N (Q µ N k ) = 1 N N i=1 δ ξ i k+1|k ,
is used, where the N sample {ξ i k+1|k , i = 1, • • • , N } has precisely the probability distribution Q µ N k , i.e. N independent r.v.'s are simulated with common probability distribution Q µ N k , which can be achieved in the following manner : independently for any

i = 1, • • • , N ξ i k ∼ µ N k (dx)
⇐= easy, since the probability distribution µ N k is discrete,

ξ i k+1|k ∼ Q(ξ i k , dx ′ )
⇐= easy, by assumption.

3.3 Linear tangent kernel / extended kernel, etc.

If the transition kernel Q(x, dx ′ ) depends on a parameter, then the lter µ k depends also on the parameter, and one would like to compute the linear tangent lter w k , i.e. the derivative of the lter µ k w.r.t. the parameter. To this end, one needs rst to study the linear tangent kernel Γ(x, dx ′ ), i.e. the derivative of the transition kernel Q(x, dx ′ ) w.r.t. the parameter, and the following assumption is made Assumption AC : The following probabilistic representation holds for the linear tangent kernel

Γ(x, dx ′ ) Γ φ(x) = E Γ(x, dx ′ ) φ(x ′ ) = E[φ(X k+1 ) Ξ k+1 | X k = x] ,
where

{ (X k , Ξ k ) , k ≥ 0} is a Markov chain taking values in the product space E × F , such that P[X k+1 ∈ dx ′ , Ξ k+1 ∈ ds ′ | X k = x, Ξ k = s] = P[X k+1 ∈ dx ′ , Ξ k+1 ∈ ds ′ | X k = x] = K(x, dx ′ , ds ′ ) .
The following assumption, which extends the similar assumption introduced in Section 3.1, is made

It is easy to simulate a r.v. (X, Ξ) with probability distribution K(x, dx ′ , ds ′ ), even though the analytical expression of the kernel K(x, dx ′ , ds ′ ) is not known, or is so complicated that it is pratically impossible to compute such integrals as

Γ φ(x) = E×F s ′ φ(x ′ ) K(x, dx ′ , ds ′ ) or Γ µ(dx ′ ) = E×F µ(dx) s ′ K(x, dx ′ , ds ′ ) . Example In our model, the Markov chain {X k , k ≥ 0} taking values in E = R d , is de ned by X k+1 = F X k + W k ,
where only the matrix F depends on the parameter, and where {W k , k ≥ 0} is a sequence of independent r.v.'s taking values in R d with probability distribution p(w) dw (coe cients are assumed time independent for simplicity). For any x ∈ R d , the transition kernel Q(x, dx ′ ) is given by

Q(x, dx ′ ) = p(x ′ -F x) dx ′ ,
and one can show directly that

Γ(x, dx ′ ) = -p ′ (x ′ -F x) ∂F ∂θ x dx ′ = -p ′ p (x ′ -F x) ∂F ∂θ x Q(x, dx ′ ) .

PI n 1604

It follows that

Γ φ(x) = E Γ(x, dx ′ ) φ(x ′ ) = E φ(x ′ ) -p ′ p (x ′ -F x) ∂F ∂θ x Q(x, dx ′ ) = E[φ(X k+1 ) -p ′ p (X k+1 -F X k ) ∂F ∂θ X k | X k = x] ,
i.e. Assumption AC is satis ed, with

Ξ k+1 = -p ′ p (W k ) ∂F ∂θ X k .
Notice that in the above example, the r.v. Ξ k+1 depends only on (X k , W k ), in which case it does not seem necessary to simulate Ξ k+1 in addition to W k . This apparently very particular situation is actually very general, as the following result shows. Lemma 3.1. Under Assumption AC

Γ(x, dx ′ ) = I(x, x ′ ) Q(x, dx ′ ) , with I(x, x ′ ) = E[Ξ k+1 | X k = x, X k+1 = x ′ ] ,
for any x, x ′ ∈ E. Proof. For any probability distribution µ on E, and any pair B, B ′ of Borel subsets of

E E µ [Ξ k+1 1 (X k+1 ∈ B ′ , X k ∈ B) ] = B×B ′ E[Ξ k+1 | X k = x, X k+1 = x ′ ] µ(dx) Q(x, dx ′ ) ,
and

E µ [Ξ k+1 1 (X k+1 ∈ B ′ ) 1 (X k ∈ B) ] = B E[Ξ k+1 1 (X k+1 ∈ B ′ ) | X k = x] µ(dx) = B Γ(x, B ′ ) µ(dx) , hence taking B = E yields E Γ(x, B ′ ) µ(dx) = E { B ′ E[Ξ k+1 | X k = x, X k+1 = x ′ ] Q(x, dx ′ ) } µ(dx) ,
and since the probability distribution µ is arbitrary, it holds

Γ(x, B ′ ) = B ′ E[Ξ k+1 | X k = x, X k+1 = x ′ ] Q(x, dx ′ ) ,
which proves the result.

By de nition

Γ φ(x) = E×F φ(x ′ ) s ′ K(x, dx ′ , ds ′ ) and Q φ(x) = E×F φ(x ′ ) K(x, dx ′ , ds ′ ) , hence Γ(x, dx ′ ) = F s ′ K(x, dx ′ , ds ′ ) and Q(x, dx ′ ) = F K(x, dx ′ , ds ′ ) .

Irisa

On the product space E × E × F , de ne the projection π 0 : (x, x ′ , s ′ ) -→ x on the ( rst) space E, the projection π : (x, x ′ , s ′ ) -→ x ′ on the (second) space E and the projection π F : (x, x ′ , s ′ ) -→ s ′ on the auxiliary space F . For any probability distribution µ on the space E, the probability distribution µ ⊗ K is de ned on the product space E × E × F by

(µ ⊗ K)(dx, dx ′ , ds ′ ) = µ(dx) K(x, dx ′ , ds ′ ) .

It follows that

Q µ(dx ′ ) = E×F µ(dx) K(x, dx ′ , ds ′ ) = (µ ⊗ K) • π -1 (dx ′ ) ,
and

Γ µ(dx ′ ) = E×F µ(dx) s ′ K(x, dx ′ , ds ′ ) = E×F π F (x, x ′ , s ′ ) (µ ⊗ K)(dx, dx ′ , ds ′ ) = (π F (µ ⊗ K)) • π -1 (dx ′ ) ,
and if the nite signed measure w is absolutely continuous w.r.t. µ, then

Q w(dx ′ ) = E×F w(dx) K(x, dx ′ , ds ′ ) = E×F dw dµ (x) µ(dx) K(x, dx ′ , ds ′ ) = E×F ( dw dµ • π 0 )(x, x ′ , s ′ ) (µ ⊗ K)(dx, dx ′ , ds ′ ) = (( dw dµ • π 0 ) (µ ⊗ K)) • π -1 (dx ′ ) , i.e. Q µ = (µ ⊗ K) • π -1 and Γ µ = (π F (µ ⊗ K)) • π -1 ,
and

( w ≪ µ =⇒ Q w = (( dw dµ • π 0 ) (µ ⊗ K)) • π -1 ) .
Lemma 3.2. Under Assumption AC, Γ µ ≪ Q µ for any probability distribution µ on E, with Radon

Nikodym derivative (which depends on µ)

d(Γ µ) d(Q µ) (x ′ ) = E µ [Ξ n+1 | X n+1 = x ′ ] .
Proof. For any Borel subset

B ′ of E, it holds Γ µ(B ′ ) = E µ(dx) Γ(x, B ′ ) = E µ(dx) E[Ξ k+1 1 (X k+1 ∈ B ′ ) | X k = x] = E µ [Ξ k+1 1 (X k+1 ∈ B ′ ) ] = B ′ E µ [Ξ k+1 | X k+1 = x ′ ] Q µ(dx ′ ) ,
hence Γ µ is a signed measure, absolutely continuous w.r.t. Q µ, with Radon Nikodym derivative

d(Γ µ) d(Q µ) (x ′ ) = E µ [Ξ k+1 | X k+1 = x ′ ]
PI n 1604

For completeness, the following elementary property is recalled Lemma 3.3. If the nite signed measure w is absolutely continuous w.r.t. the probability distribution µ,

then Q w ≪ Q µ, with Radon Nikodym derivative d(Q w) d(Q µ) (x ′ ) = E µ [ dw dµ (X k ) | X k+1 = x ′ ] .
Proof. For any Borel subset

B ′ of E, it holds Q w(B ′ ) = E w(dx) Q(x, B ′ ) = E dw dµ (x) µ(dx) Q(x, B ′ ) = E µ [ dw dµ (X k ) 1 (X k+1 ∈ B ′ ) ] = B ′ E µ [ dw dµ (X k ) | X k+1 = x ′ ] Q µ(dx ′ ) ,
hence Q w is a signed measure, absolutely continuous w.r.t. Q µ, with Radon Nikodym derivative

d(Q w) d(Q µ) (x ′ ) = E µ [ dw dµ (X k ) | X k+1 = x ′ ]
The explicit expression of the Radon Nikodym derivatives will not be used in the sequel : only the qualitative properties

Γ µ ≪ Q µ and ( w ≪ µ =⇒ Q w ≪ Q µ ) ,
will be used.

By de nition

F k (µ) w = Ψ k w µ, Ψ k - w, Ψ k µ, Ψ k Ψ k µ µ, Ψ k ,
is the derivative at point µ and in the direction w, of the mapping µ -→ Ψ k • µ.

then F k (µ) w ≪ Ψ k • µ, with Radon Nikodym derivative d(F k (µ) w) d(Ψ k • µ) (x) = dw dµ (x) -Ψ k • µ, dw dµ .

Particle approximation of some nite signed measures

With the notations of the previous section, it easily seen that the probability distribution Q µ and the nite signed measures Γ µ and Q w can be put in the general form (r (µ ⊗ K)) • π -1 for some appropriate choice of the weight function r, namely r ≡ 1, r = π F and r = dw dµ • π 0 respectively. The weighted particle approximation of a nite signed measure of the general form r (µ ⊗ K) is de ned by

r (µ ⊗ K) ≈ r S N (µ ⊗ K) = 1 N N i=1 r(ξ i 0 , ξ i , Ξ i ) δ (ξ i 0 , ξ i , Ξ i ) ,

Irisa

Particle based Methods for Parameter Estimation and Tracking 9 where the N sample { ξ i 0 , ξ i , Ξ i , i = 1, • • • , N } has precisely the probability distribution µ ⊗ K, i.e. one simulate N independent r.v.'s with common probability distribution µ ⊗ K, which can be achieved in the following manner : independently for any

i = 1, • • • , N ξ i 0 ∼ µ(dx)
and

(ξ i , Ξ i ) ∼ K(ξ i 0 , dx ′ , ds ′ ) ,
and the corresponding particle approximation for the marginal measure (r

(K ⊗ µ)) • π -1 is de ned by (r (µ ⊗ K)) • π -1 ≈ (r S N (µ ⊗ K)) • π -1 = 1 N N i=1 r(ξ i 0 , ξ i , Ξ i ) δ ξ i .
In particular for the weight functions r ≡ 1, r = π F and r = dw dµ

• π 0 , it holds Q µ = (µ ⊗ K) • π -1 ≈ 1 N N i=1 δ ξ i , Γ µ = (π F (µ ⊗ K)) • π -1 ≈ 1 N N i=1 Ξ i δ ξ i ,
and

Q w = (( dw dµ • π 0 ) (µ ⊗ K)) • π -1 ≈ 1 N N i=1 dw dµ (ξ i 0 ) δ ξ i ,
respectively. For any test function φ de ned on E, it holds

E| 1 N N i=1 r(ξ i 0 , ξ i , Ξ i ) φ(ξ i ) -(r (µ ⊗ K)) • π -1 , φ | ≤ 1 √ N { E×E×F |φ(x ′ )| 2 |r(x, x ′ , s ′ )| 2 µ(dx) K(x, dx ′ , ds ′ ) } 1/2 ≤ 1 √ N { E×E×F |r(x, x ′ , s ′ )| 2 µ(dx) K(x, dx ′ , ds ′ ) } 1/2 φ ,
and in particular for the weight functions r ≡ 1, r = π F and r = dw dµ

• π 0 , it holds sup φ =1 E| 1 N N i=1 φ(ξ i ) -Q µ, φ | ≤ 1 √ N , sup φ =1 E| 1 N N i=1 Ξ i φ(ξ i ) -Γ µ, φ | ≤ 1 √ N {sup x∈E E×F |s ′ | 2 K(x, dx ′ , ds ′ ) } 1/2 ,
and

sup φ =1 E| 1 N N i=1 dw dµ (ξ i 0 ) φ(ξ i ) -Q w, φ | ≤ 1 √ N { E | dw dµ (x)| 2 µ(dx) } 1/2 ,
respectively.

PI n 1604 3.5 Joint particle approximation of the lter and the linear tangent lter

Recall that the evolution of the sequence {µ k , k ≥ 0} taking values in the space of probability distributions on E, is described by the following two steps

µ k-1 prediction ---------→ µ k|k-1 = Q µ k-1 correction ---------→ µ k = Ψ k • µ k|k-1 .
As we want to estimate H, Ψ k also depends on θ. If w k denotes at each time instant the linear tangent lter, i.e. the derivative of the lter µ k w.r.t. the parameter, then the evolution of the sequence {w k , k ≥ 0} taking values in the linear tangent space to the space of probability distributions on E, i.e. taking values in the space of nite signed measures on E with zero total mass, is described by the following two steps, which are linear tangent versions of the prediction step and correction step respectively

w k-1 linear tangent prediction ------------→ w k|k-1 = Q w k-1 + Γ µ k-1 linear tangent correction ------------→ w k = F k (µ k|k-1 ) w k|k-1 + G k (µ k|k-1 ) .
where

G k (µ) = ∂Ψ k ∂θ µ µ, Ψ k - µ, ∂Ψ k ∂θ µ, Ψ k Ψ k µ µ, Ψ k = [ ∂ log Ψ k ∂θ -Ψ k • µ, ∂ log Ψ k ∂θ ] Ψ k • µ .
Under Assumption AC, it is easily seen by induction, and using Lemmas 3.2, 3.3 and 3.4, that at each time instant w k|k-1 ≪ µ k|k-1 and w k ≪ µ k . In view of this absolute continuity property, and of the key assumption that it is easy to simulate r.v.'s with probability distribution K(x, dx ′ , ds ′ ), the idea is to jointly approximate the predictor µ k|k-1 and its derivative w k|k-1 w.r.t. the parameter with the empirical probability distribution and with a weighted empirical distribution associated with the same and unique N sample, i.e.

µ k|k-1 ≈ µ N k|k-1 = 1 N N i=1 δ ξ i k|k-1
and

w k|k-1 ≈ w N k|k-1 = 1 N N i=1 ρ i k|k-1 δ ξ i k|k-1
.

With this de nition w N k|k-1 ≪ µ N k|k-1 , with Radon Nikodym derivative

r N k|k-1 (x) = dw N k|k-1 dµ N k|k-1 (x) = 1 |I N k|k-1 (x)| i∈I N k|k-1 (x) ρ i k|k-1 ,
where

I N k|k-1 (x) = {i = 1, • • • , N : ξ i k|k-1 = x}, for any x in the support supp µ N k|k-1 of the discrete probability distribution µ N k|k-1 . Indeed µ N k|k-1 = 1 N x∈supp µ N k|k-1 |I N k|k-1 (x)| δ x ,
and

w N k|k-1 = 1 N x∈supp µ N k|k-1 [ i∈I N k|k-1 (x) ρ i k|k-1 ] δ x .
Notice that in most cases, the particle locations {ξ i k|k-1 , i = 1, • • • , N } turn out to be all distinct, and the much simpler relation

r N k|k-1 (ξ i k|k-1 ) = ρ i k|k-1 ,

Irisa holds for any

i = 1, • • • , N .
This approximation is completely characterized by the set {ξ i k|k-1 , ρ i k|k-1 , i = 1, • • • , N } of particles and weights, and the algorithm is completely described by the mechanism which builds

{ξ i k+1|k , ρ i k+1|k , i = 1, • • • , N } from {ξ i k|k-1 , ρ i k|k-1 , i = 1, • • • , N }.
This mechanism is as follows : (i) the correction step is applied exactly to µ N k|k-1 , which results in

µ N k = Ψ k • µ N k|k-1 = N i=1 Ψ k (ξ i k|k-1 ) δ ξ i k|k-1 N j=1 Ψ k (ξ j k|k-1 ) = N i=1 ω i k δ ξ i k|k-1 ,
as previously, and the linear tangent correction step is applied exactly to w N k|k-1 , which results in

w N k = F k (µ N k|k-1 ) w N k|k-1 + G k (µ N k|k-1 ) = [ r N k|k-1 -Ψ k • µ N k|k-1 , r N k|k-1 ] Ψ k • µ N k|k-1 + [ ∂ log Ψ k ∂θ -Ψ k • µ N k|k-1 , ∂ log Ψ k ∂θ ] Ψ k • µ N k|k-1 = [ r N k|k-1 + ∂ log Ψ k ∂θ -µ N k , r N k|k-1 + ∂ log Ψ k ∂θ ] µ N k ,
(ii) instead of trying to compute Q µ N k , the following particle approximation

µ N k+1|k = S N (Q µ N k ) = 1 N N i=1 δ ξ i k+1|k ,
is used as previously, instead of trying to compute

Q w N k = Q((r N k|k-1 + ∂ log Ψ k ∂θ ) µ N k ) -µ N k , r N k|k-1 + ∂ log Ψ k ∂θ Q µ N k = ( ((r N k|k-1 + ∂ log Ψ k ∂θ ) • π 0 ) (µ N k ⊗ Q) ) • π -1 -µ N k , r N k|k-1 + ∂ log Ψ k ∂θ Q µ N k ,
the following weighted particle approximation

( ((r N k|k-1 + ∂ log Ψ k ∂θ ) • π 0 ) S N (µ N k ⊗ Q) ) • π -1 -S N (µ N k ), r N k|k-1 + ∂ log Ψ k ∂θ S N (Q µ N k ) = 1 N N i=1 [ r N k|k-1 (ξ i k ) + ∂ log Ψ k ∂θ (ξ i k ) - 1 N N j=1 [ r N k|k-1 (ξ j k ) + ∂ log Ψ k ∂θ (ξ j k ) ] ] δ ξ i k+1|k ,
is used, and instead of trying to compute

Γ µ N k = (π F (µ N k ⊗ K)) • π -1 ,
the following weighted particle approximation

(π F S N (µ N k ⊗ K)) • π -1 = 1 N N i=1 Ξ i k+1 δ ξ i k+1|k ,
is used, hence nally the weighted particle approximation

w N k+1|k = 1 N N i=1 [ r N k|k-1 (ξ i k ) + ∂ log Ψ k ∂θ (ξ i k ) - 1 N N j=1 [ r N k|k-1 (ξ j k ) + ∂ log Ψ k ∂θ (ξ j k ) ] + Ξ i k+1 ] ρ i k+1|k δ ξ i k+1|k , PI n 1604
where the

N sample { ξ i k , ξ i k+1|k , Ξ i k+1 , i = 1, • • • , N } has precisely the probability distribution µ N k ⊗K,
i.e. N independent r.v.'s are simulated with common probability distribution µ N k ⊗ K, which can be achieved in the following manner : independently for any

i = 1, • • • , N ξ i k ∼ µ N k (dx)
⇐= easy, since the probability distribution µ N k is discrete,

(ξ i k+1|k , Ξ i k+1 ) ∼ K(ξ i
k , dx ′ , ds ′ ) ⇐= easy, by assumption.

The proposed particle approximation of the linear tangent optimal lter is especially attractive, since it uses the same particle system already used in the approximation of the optimal lter : only one dimensional weights are needed in addition. Other closely related particle approximation schemes will be presented in Sections 4.1, 4.2, 4.3 and 4.4 on some simple models where comparison with explicit exact expressions provided by Kalman ltering is possible.

Di erent redistribution schemes

In practice, after a few iterations all but one particle will have negligible weights, and a large computational e ort is devoted to updating particles whose contribution is almost zero. To avoid this phenomenon of degeneracy of particle weights, we resort to resampling. The basic idea of resampling, and more generally of redistribution, is to eliminate particles with small weights and to replicate particles with large weights.

Resampling

In the selection step, we want to generate a N sample {ξ

i k+1 , i = 1, • • • , N } with discrete probability distribution N i=1 ω i k δ ξ i k .
One of the most direct method is based on an inversion method and consists in generating a uniform r.v. U on [0, 1], and if

ω 1 k + • • • ω j k ≤ U < ω 1 k + • • • ω j+1 k then choose ξ j k .
The new weigths are reset to 1/N , hence the following algorithm :

Resampling

Let {ξ i k , ω i k , i = 1, • •
• , N } be the particle system we want to resample from. Construct cumulative distribution function (CDF) :

c 1 = 0. For i = 2, • • • , N , set c i = c i-1 + ω i k . For i = 1, • • • , N , U ∼ U [0, 1].
Find α i so that c αi ≤ U < c αi+1 . The new particle system is given by : For

i = 1, • • • , N , set ξ i k+1 = ξ αi k and ω i k+1 = 1/N .
Some heuristic measures have been proposed to decide when it is necessary to resample. In this work, we choose to resample at each time iteration.

Irisa

Particle based Methods for Parameter Estimation and Tracking 13 3.6.2 Comb method A more determistic way to resample is to multiply ⌊N ω i k ⌋ times the particle ξ i k . Then, we use a comb method to select

N i=1 (N ω i k -⌊N ω i k ⌋
) particles as described in the following algorithm :

Comb method

Let {ξ i k , ω i k , i = 1, • •
• , N } be the particle system we want to redistribute from. Let σ be a random permutation on {1, • • • , N }.

For i = 1, • • • , N , replace i with σ(i). Construct CDF : Set α 1 = ω σ(1) k - 1 N ⌊N ω σ(1) k ⌋. For i = 2, • • • , N , compute α i = α i-1 + ω σ(i) k - 1 N ⌊N ω σ(i) k ⌋.
The new particle system is given by :

For i = 1, • • • , N -1, replicate the particle ξ σ(i) k a number of times equal to 1 + ⌊N ω σ(i) k ⌋ if α i < i/N ≤ α i+1 , and to ⌊N ω σ(i) k
⌋ otherwise, and keep the particle

ξ σ(N ) k . For i = 1, • • • , N , set ω i k+1 = 1/N .
Notice that the last particle is always selected in this algorithm, and the random permutation is introduced so as to avoid this undesirable systematic e ect.

An alternative comb method

Another deterministic algorithm has been proposed :

PI n 1604 Alternative comb method Let {ξ i k , ω i k , i = 1, • • • , N } be the particle system we want to redistribute from. Let σ be a random permutation on {1, • • • , N }. For i = 1, • • • , N , replace i with σ(i).

Construct CDF :

Set

α 1 = ω σ(1) k . For i = 2, • • • , N , compute α i = α i-1 + ω σ(i) k . For i = 1, • • • , N , let β i denote the nearest integer to N α i . Construct : Set γ 1 = β 1 . For i = 2, • • • , N , set γ i = β i -β i-1
. The new particle system is given by :

For i = 1, • • • , N , replicate the particle ξ σ(i) k
a number of times equal to γ i , and set

ω i k+1 = 1/N .

Conclusions

We have seen several methods to redistribute our particle system. The rst one seems very natural but is slower than the two others. In practice, we will use a deterministic method since the particle approximations produced by all three methods are very similar, as the following experimental results show :

Particle approximation with 500 particles The purpose of this section is to obtain as much intuition as possible, through numerical experiments in simple cases where exact solutions are explicitly available, about the particle approximation of nite signed measures. A prototypical example of a nite signed measure is the derivative, w.r.t. a parameter of the model, of some probability distribution related with a Markov chain. This includes prior, prediction, ltering probability distributions, and the like. Two points of view are considered here, to feel the quality of the approximation, at least in a qualitative manner : (i) how accurate is the particle approximation of the nite signed measure, in view of an histogram representation of the weighted particle system ? (ii) considering the log likelihood function and the score function, how close is the approximate expression provided by the particle approximation to the exact expression ? These two questions seem closely related, however the numerical experiments presented in this section show that one of the two particle approximation schemes fails to satisfy the rst criteria (quality of the approximation of the nite signed measure), and that both schemes satisfy the second criteria (quality of the approximation of the statistics). Explanations about the graphical outputs are provided at the end of this section.

Part A : AR(1) model

The model considered in this example is the simplest scalar AR(1) model

X k+1 = a X k + σ W k , X 0 ∼ N( X0 , σ 2 0 ) , i.e. µ 0 (dx) = P[X 0 ∈ dx] = 1 √ 2 π σ 0 exp{- 1 2 σ 2 0 (x -X0 ) 2 } p 0 (x) dx ,
where {W k , k ≥ 0} is a standard Gaussian white noise sequence, and where only a is considered as an unknown parameter, i.e. the data X0 and σ 2 are known and σ 2 0 could depend on a. It follows from the model that

X k+1 | X k = x ∼ N(a x, σ 2 ) , i.e. Q(x, dx ′ ) = P[X k+1 ∈ dx ′ | X k = x] = 1 √ 2 π σ exp{- 1 2 σ 2 (x ′ -a x) 2 } q(x, x ′ ) dx ′ .
The goal here is to design a particle approximation scheme for the derivative, with respect to the parameter a, of the probability distribution of the state.

Exact expressions

Obviously, in such a simple model, an exact expression can easily be obtained. Indeed

X k ∼ N( Xk , σ 2 k ) , i.e. µ k (dx) = P[X k ∈ dx] = 1 √ 2 π σ k exp{- 1 2 σ 2 k (x -Xk ) 2 } p k (x) dx ,
where

Xk+1 = a Xk and σ 2 k+1 = a 2 σ 2 k + σ 2 .
The log density of the r.v.

X k is log p k (x) = cste - 1 2 log σ 2 k - 1 2 σ 2 k (x -Xk ) 2 , hence the logarithmic derivative w.r.t. the parameter a ∂ log p k ∂a (x) = - 1 2 σ 2 k ∂σ 2 k ∂a + 1 2 σ 4 k ∂σ 2 k ∂a (x -Xk ) 2 + 1 σ 2 k ∂ Xk ∂a (x -Xk ) ,
and

w k (dx) = ∂µ k ∂a (dx) = [ 1 2 σ 4 k ∂σ 2 k ∂a [ (x -Xk ) 2 -σ 2 k ] + 1 σ 2 k ∂ Xk ∂a (x -Xk ) ] µ k (dx) ,

Irisa

where

∂ Xk+1 ∂a = Xk + a ∂ Xk ∂a and ∂σ 2 k+1 ∂a = 2 a σ 2 k + a 2 ∂σ 2 k ∂a .
In the special case where X0 = 0, it holds Xk = 0, ∂ Xk ∂a = 0, hence

w k (dx) = 1 2 σ 4 k ∂σ 2 k ∂a (x 2 -σ 2 k ) µ k (dx) ,
and the support of w + k , the positive part of w k , is the set {x :

dw k dµ k (x) > 0} = {x : |x| > σ k }.
In the stationary case where |a| < 1 and X0 = 0, it holds

σ 2 k = σ 2 0 = σ 2 1 -a 2 and ∂σ 2 0 ∂a = ∂σ 2 k ∂a = 2 a σ 2 (1 -a 2 ) 2 , hence µ k (dx) = µ 0 (dx)
and

w k (dx) = w 0 (dx) = a σ 2 (x 2 -σ 2 0 ) µ 0 (dx) .

Preliminary computations

The transition log density is

log q(x, x ′ ) = cste - 1 2 log σ 2 - 1 2 σ 2 (x ′ -a x) 2 , hence the logarithmic derivative w.r.t. the parameter a ∂ log q ∂a (x, x ′ ) = x σ 2 (x ′ -a x) ,
and

Γ(x, dx ′ ) = ∂Q ∂a (x, dx ′ ) = x σ 2 (x ′ -a x) Q(x, dx ′ ) .
Di erentiating with respect to the parameter a, throughout the recursion

µ k+1 (dx ′ ) = Q µ k (dx ′ ) = ∞ -∞ µ k (dx) Q(x, dx ′ ) , (10)
yields

w k+1 (dx ′ ) = ∂µ k+1 ∂a (dx ′ ) = ∞ -∞ [ ∂µ k ∂a (dx) Q(x, dx ′ ) + µ k (dx) ∂Q ∂a (x, dx ′ ) ] ,
hence

w k+1 (dx ′ ) = Q w k (dx ′ ) + Γ µ k (dx ′ ) = ∞ -∞ [ w k (dx) + x σ 2 (x ′ -a x) µ k (dx) ] Q(x, dx ′ ) . (11)
Since

w k ≪ µ k , it holds w k+1 (dx ′ ) = Q w k (dx ′ ) + Γ µ k (dx ′ ) = ∞ -∞ [ dw k dµ k (x) + x σ 2 (x ′ -a x) ] µ k (dx) Q(x, dx ′ ) ,
i.e. w k+1 can be interpreted as the marginal of the nite signed measure w k,k+1 de ned on the product space R × R by

w k,k+1 (dx, dx ′ ) = w k (dx) Q(x, dx ′ ) + µ k (dx) Γ(x, dx ′ ) = [ dw k dµ k (x) + x σ 2 (x ′ -a x) ] µ k (dx) Q(x, dx ′ ) . (12)
PI n 1604 4.1.3 A particle approximation scheme Assuming that the following particle approximations

µ k ≈ µ N k = 1 N N i=1 δ ξ i k and w k ≈ w N k = 1 N N i=1 ρ i k δ ξ i k ,
are available at time index k, and plugging these approximations into equations ( 10) and ( 11), yields

Q µ N k (dx ′ ) = 1 N N i=1 Q(ξ i k , dx ′ ) , (13)
and

Q w N k (dx ′ ) + Γ µ N k (dx ′ ) = 1 N N i=1 [ ρ i k + ξ i k σ 2 (x ′ -a ξ i k ) ] Q(ξ i k , dx ′ ) , (14)
which can be interpreted as the marginals of the nite signed

measures m = (m i , i = 1, • • • , N ) and s = (s i , i = 1, • • • , N ) de ned on the product space E N = {1, • • • , N } × R by s i (dx ′ ) = [ ρ i k + ξ i k σ 2 (x ′ -a ξ i k ) ] r i (x ′ ) 1 N Q(ξ i k , dx ′ ) m i (dx ′ )
, respectively. Introducing the particle approximation

m ≈ m N = 1 N N i=1 δ (τ i k , ξ i k+1 )
and the weighted particle approximation

s ≈ s N = 1 N N i=1 r τ i k (ξ i k+1 ) δ (τ i k , ξ i k+1 ) ,
where independently for any i = 1, • • • , N , the pair (τ i k , ξ i k+1 ) is jointly distributed according to the probability distribution m = (m i , i = 1, • • • , N ), or using systematic sampling alternatively

τ i k = i and ξ i k+1 ∼ Q(ξ τ i k k , dx ′ ) ,
and marginalizing, yields

µ k+1 ≈ µ N k+1 = 1 N N i=1 δ ξ i k+1 ,
and

w k+1 ≈ w N k+1 = 1 N N i=1 r i (ξ i k+1 ) δ ξ i k+1 = 1 N N i=1 ρ i k+1 δ ξ i k+1 ,
where independently for any

i = 1, • • • , N ξ i k+1 ∼ Q(ξ i k , dx ′ ) .

Irisa

Particle based Methods for Parameter Estimation and Tracking 19 Notice that the weight

ρ i k+1 = r i (ξ i k+1 ) = ρ i k + ξ i k σ 2 (ξ i k+1 -a ξ i k ) ,
does not depend only on the position of ξ i k+1 , but on the position of the pair

(ξ i k , ξ i k+1
). An approximation of the support of w + k+1 , the positive part of w k+1 , is provided by the set {ξ i k+1 , i ∈ I N,+ k+1 } of particles with positive weight, where

I N,+ k+1 = {i = 1, • • • , N : ρ i k + ξ i k σ 2 (ξ i k+1 -a ξ i k ) > 0} .
4.1. 4 An alternate particle approximation scheme Equations ( 13) and ( 14) read

Q µ N k (dx ′ ) = 1 N N i=1 Q(ξ i k , dx ′ ) = [ 1 N N i=1 q(ξ i k , x ′ ) ] dx ′ ,
and

Q w N k (dx ′ ) + Γ µ N k (dx ′ ) = 1 N N i=1 [ ρ i k + ξ i k σ 2 (x ′ -a ξ i k ) ] Q(ξ i k , dx ′ ) = [ 1 N N i=1 [ ρ i k + ξ i k σ 2 (x ′ -a ξ i k ) ] q(ξ i k , x ′ ) ] dx ′ = r N k+1 (x ′ ) Q µ N k (dx ′ ) ,
where

r N k+1 (x ′ ) = N j=1 [ ρ j k + ξ j k σ 2 (x ′ -a ξ j k ) ] q(ξ j k , x ′ ) N j=1 q(ξ j k , x ′ ) = N j=1 [ ρ j k + ξ j k σ 2 (x ′ -a ξ j k ) ] exp{-1 2 (x ′ -a ξ j k ) 2 } N j=1 exp{-1 2 (x ′ -a ξ j k ) 2 } .
Using systematic sampling, yields the particle approximation

Q µ N k ≈ 1 N N i=1 δ ξ i k+1 ,
where independently for

any i = 1, • • • , N ξ i k+1 ∼ Q(ξ i k , dx ′ ) ,
and the weighted particle approximation

Q w N k + Γ µ N k ≈ w N k+1 = 1 N N i=1 r N k+1 (ξ i k+1 ) δ ξ i k+1 = 1 N N i=1 ρ i k+1 δ ξ i k+1 , PI n 1604
where the weight

ρ i k+1 = r N k+1 (ξ i k+1 ) = N j=1 [ ρ j k + ξ j k σ 2 (ξ i k+1 -a ξ j k ) ] exp{-1 2 (ξ i k+1 -a ξ j k ) 2 } N j=1 exp{-1 2 (ξ i k+1 -a ξ j k ) 2 }
, depends only on the position of ξ i k+1 . An approximation of the support of w + k+1 , the positive part of w k+1 , is provided by the set {ξ i k+1 , i ∈ I N,+ k+1 } of particles with positive weight, where

I N,+ k+1 = {i = 1, • • • , N : r N k+1 (ξ i k+1 ) > 0} .

Part B : HMM situation

Consider now the HMM situation, where the state, still described by the same scalar AR(1) model as above, is not observed directly and where observations are available instead, which are related to the hidden state by

Y k = c X k + s V k ,
where {V k , k ≥ 0} is a another standard Gaussian white noise sequence, independent of {W k , k ≥ 0}, and where only a is considered as an unknown parameter, i.e. the data X0 , σ 0 , σ, c and s are known. It follows from the model that

Y k | X k = x ∼ N(c x, s 2 ) , i.e. P[Y k ∈ dy | X k = x] = 1 √ 2 π s exp{- 1 2 σ 2 (y -c x) 2 } g(x, y) dy ,
and let

Ψ k (x) = g(x, Y k ) ,
denote the likelihood function. The goal here is to design a particle approximation scheme for the derivative, with respect to the parameter a, of the conditional probability distribution of the hidden state, given the observations.

Exact expressions (Kalman lter)

In such a simple model, an exact expression can easily be obtained, via the prediction / correction steps of the Kalman lter framework. Indeed, the prediction step reads

X k | Y 1 , • • • , Y k-1 ∼ N( Xk|k-1 , σ 2 k|k-1 ) , i.e. µ k|k-1 (dx) = P[X k ∈ dx | Y 1 , • • • , Y k-1 ] = 1 √ 2 π σ k|k-1 exp{- 1 2 σ 2 k|k-1 (x -Xk|k-1 ) 2 } p k|k-1 (x) dx ,
where

Xk|k-1 = a Xk-1 and σ 2 k|k-1 = a 2 σ 2 k-1 + σ 2 .

Irisa

Particle based Methods for Parameter Estimation and Tracking 21 Proceeding exactly as above yields

w k|k-1 (dx) = ∂µ k|k-1 ∂a (dx) = [ 1 2 σ 4 k|k-1 ∂σ 2 k|k-1 ∂a [ (x -Xk|k-1 ) 2 -σ 2 k|k-1 ] + 1 σ 2 k|k-1 ∂ Xk|k-1 ∂a (x -Xk|k-1 ) ] µ k|k-1 (dx) ,
where

∂ Xk|k-1 ∂a = Xk-1 + a ∂ Xk-1 ∂a and ∂σ 2 k|k-1 ∂a = 2 a σ 2 k-1 + a 2 ∂σ 2 k-1 ∂a .
The correction step reads

X k | Y 1 , • • • , Y k ∼ N( Xk , σ 2 k ) ,
i.e.

µ k (dx) = P[X k ∈ dx | Y 1 , • • • , Y k ] = 1 √ 2 π σ k exp{- 1 2 σ 2 k (x -Xk ) 2 } p k (x) dx ,
where the Kalman gain is de ned by

K k = c σ 2 k|k-1 c 2 σ 2 k|k-1 + s 2 hence 1 -K k c = s 2 c 2 σ 2 k|k-1 + s 2 ,
and where

Xk = Xk|k-1 + K k (Y k -c Xk|k-1 ) = Xk|k-1 + c σ 2 k|k-1 c 2 σ 2 k|k-1 + s 2 (Y k -c Xk|k-1 ) ,
and

σ 2 k = (1 -K k c) σ 2 k|k-1 = s 2 σ 2 k|k-1 c 2 σ 2 k|k-1 + s 2 .
Therefore

w k (dx) = ∂µ k ∂a (dx) = [ 1 2 σ 4 k ∂σ 2 k ∂a [ (x -Xk ) 2 -σ 2 k ] + 1 σ 2 k ∂ Xk ∂a (x -Xk ) ] µ k (dx) ,
where

∂ Xk ∂a = s 2 c 2 σ 2 k|k-1 + s 2 ∂ Xk|k-1 ∂a + c s 2 (c 2 σ 2 k|k-1 + s 2 ) 2 ∂σ 2 k|k-1 ∂a (Y k -c Xk|k-1 ) ,
and

∂σ 2 k ∂a = s 4 (c 2 σ 2 k|k-1 + s 2 ) 2 ∂σ 2 k|k-1 ∂a .
From the decomposition

Y k = c Xk|k-1 + I k , PI n 1604
where the innovation

I k = Y k -c Xk|k-1 = c (X k -Xk|k-1 ) + V k ∼ N(0, s 2 k )
with

s 2 k = c 2 σ 2 k|k-1 + s 2 , is independent of the past observations Y 1 , • • • , Y k-1 , it follows that Y k | Y 1 , • • • , Y k-1 ∼ N(c Xk|k-1 , s 2 k ) ,
i.e.

P[Y k ∈ dy | Y 1 , • • • , Y k-1 ] = 1 √ 2 π s k exp{- 1 2 s 2 k (y -c Xk|k-1 ) 2 } g k (y) dy .
Using the straightforward identity

P[Y 1 ∈ dy 1 , • • • , Y n ∈ dy n ] = n k=1 P[Y k ∈ dy k | Y 1 = y 1 , • • • , Y k-1 = y k-1 ] = n k=1 g k (y k ) dy k ,
yields the following expression for the (suitably normalized) log likelihood function

ℓ n = 1 n log n k=1 g k (Y k ) = 1 n n k=1 log g k (Y k ) = - 1 n n k=1 [ 1 2 log(2πs 2 k ) + 1 2 s 2 k (Y k -c Xk|k-1 ) 2 ] ,
and for the score function

∂ℓ n ∂a = 1 n n k=1 [ - 1 2 s 2 k ∂s 2 k ∂a + 1 2 s 4 k ∂s 2 k ∂a (Y k -c Xk|k-1 ) 2 + c s 2 k ∂ Xk|k-1 ∂a (Y k -c Xk|k-1 ) ] ,
where

∂s 2 k ∂a = c 2 ∂σ 2 k|k-1 ∂a .
Notice that in full generality

P[Y k ∈ dy | Y 1 , • • • , Y k-1 ] = ∞ -∞ P[Y k ∈ dy, X k ∈ dx | Y 1 , • • • , Y k-1 ] = ∞ -∞ P[Y k ∈ dy | X k = x] P[X k ∈ dx | Y 1 , • • • , Y k-1 ] = [ ∞ -∞ g(x, y) µ k|k-1 (dx) g k (y)
] dy , hence the following equivalent expression for the log likelihood function

ℓ n = 1 n n k=1 log g k (Y k ) = 1 n n k=1 log ∞ -∞ Ψ k (x) µ k|k-1 (dx) , (15)
and for the score function

∂ℓ n ∂a = 1 n n k=1 ∞ -∞ Ψ k (x) ∂µ k|k-1 ∂a (dx) ∞ -∞ Ψ k (x) µ k|k-1 (dx) = 1 n n k=1 ∞ -∞ Ψ k (x) w k|k-1 (dx) ∞ -∞ Ψ k (x) µ k|k-1 (dx) . ( 16
) Irisa 4.2.

Preliminary computations

Di erentiating with respect to the parameter a, throughout the recursions

µ k (dx) = (Ψ k • µ k|k-1 )(dx) = Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) , (17)
and

µ k+1|k (dx) = Q µ k (dx ′ ) = ∞ -∞ µ k (dx) Q(x, dx ′ ) , (18)
yields

w k (dx) = ∂µ k ∂a (dx) = Ψ k (x) ∂µ k|k-1 ∂a (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) - ∞ -∞ Ψ k (x ′ ) ∂µ k|k-1 ∂a (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) , hence w k (dx) = (F k (µ k|k-1 )w k|k-1 )(dx) = Ψ k (x) w k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) - ∞ -∞ Ψ k (x ′ ) w k|k-1 (dx ′ ) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) , (19)
and

w k+1|k (dx ′ ) = Q w k (dx ′ ) + Γ µ k (dx ′ ) = ∞ -∞ [ w k (dx) + x σ 2 (x ′ -a x) µ k (dx) ] Q(x, dx ′ ) . (20)
Since

w k|k-1 ≪ µ k|k-1 , it holds F k (µ k|k-1 )w k|k-1 ≪ Ψ k • µ k|k-1 , with Radon Nikodym derivative dw k dµ k (x) = d(F k (µ k|k-1 )w k|k-1 ) d(Ψ k • µ k|k-1 ) (x) = dw k|k-1 dµ k|k-1 (x) - ∞ -∞ dw k|k-1 dµ k|k-1 (x ′ ) (Ψ k • µ k|k-1 )(dx ′ ) ,
and

w k+1|k (dx ′ ) = Q w k (dx ′ ) + Γ µ k (dx ′ ) = ∞ -∞ [ dw k dµ k (x) + x σ 2 (x ′ -a x) ] µ k (dx) Q(x, dx ′ ) .

A particle approximation scheme

Assuming that the following particle approximations

µ k|k-1 ≈ µ N k|k-1 = 1 N N i=1 δ ξ i k|k-1
and

w k|k-1 ≈ w N k|k-1 = 1 N N i=1 ρ i k|k-1 δ ξ i k|k-1 , (21)
are available at time index k, and plugging these approximations into equations ( 17) and [START_REF] Juang | Applied System Identi cation[END_REF], yields

µ N k = Ψ k • µ N k|k-1 = N i=1 Ψ k (ξ i k|k-1 ) δ ξ i k|k-1 N j=1 Ψ k (ξ j k|k-1 ) = N i=1 ω i k δ ξ i k|k-1 , PI n 1604
and

w N k = F k (µ N k|k-1 )w N k|k-1 = N i=1 ρ i k|k-1 Ψ k (ξ i k|k-1 ) δ ξ i k|k-1 N j=1 Ψ k (ξ j k|k-1 ) - N j=1 ρ j k|k-1 Ψ k (ξ j k|k-1 ) N j=1 Ψ k (ξ j k|k-1 ) N i=1 Ψ k (ξ i k|k-1 ) δ ξ i k|k-1 N j=1 Ψ k (ξ j k|k-1 ) = N i=1 ρ i k|k-1 ω i k δ ξ i k|k-1 -[ N j=1 ρ j k|k-1 ω j k ] N i=i ω i k δ ξ i k|k-1 = N i=1 [ ρ i k|k-1 -[ N j=1 ρ j k|k-1 ω j k ] ] ω i k δ ξ i k|k-1
.

Plugging these expressions into equations ( 18) and ( 20), yields

Q µ N k (dx ′ ) = N i=1 ω i k Q(ξ i k|k-1 , dx ′ ) , (22)
and

Q w N k (dx ′ ) + Γ µ N k (dx ′ ) = N i=1 [ ρ i k|k-1 -[ N j=1 ρ j k|k-1 ω j k ] + ξ i k|k-1 σ 2 (x ′ -a ξ i k|k-1 ) ] ω i k Q(ξ i k|k-1 , dx ′ ) , (23)
which can be interpreted as the marginals of the nite signed measures

m = (m i , i = 1, • • • , N ) and s = (s i , i = 1, • • • , N ) de ned on the product space E N = {1, • • • , N } × R by s i (dx ′ ) = [ ρ i k|k-1 + ξ i k|k-1 σ 2 (x ′ -a ξ i k|k-1 ) r i (x ′ ) -[ N j=1 ρ j k|k-1 ω j k ] ] ω i k Q(ξ i k|k-1 , dx ′ ) m i (dx ′ ) . Notice that N j=1 ∞ -∞ r j (x ′ ) m j (dx ′ ) = N j=1 ρ j k|k-1 ω j k ,
is just a normalizing constant. Introducing the particle approximation

m ≈ m N = 1 N N i=1 δ (τ i k , ξ i k+1|k )
, and the weighted particle approximation

s ≈ s N = 1 N N i=1 [ r τ i k (ξ i k+1 ) -[ 1 N N j=1 r τ j k (ξ j k+1 ) ] ] δ (τ i k , ξ i k+1|k )
,

where independently for any

i = 1, • • • , N , the pair (τ i k , ξ i k+1|k
) is jointly distributed according to the probability distribution m = (m i , i = 1, • • • , N ), i.e.

τ i k ∼ (ω j k , j = 1, • • • , N )
and

ξ i k+1|k ∼ Q(ξ τ i k k|k-1 , dx ′ ) ,

Irisa

and marginalizing, yields

µ k+1|k ≈ µ N k+1|k = 1 N N i=1 δ ξ i k+1|k . w k+1|k ≈ w N k+1|k = 1 N N i=1 [ r τ i k (ξ i k+1|k ) -[ 1 N N j=1 r τ j k (ξ j k+1|k ) ] ] δ ξ i k+1 = 1 N N i=1 ρ i k+1|k δ ξ i k+1|k .
Notice that the weight

ρ i k+1|k = r τ i k (ξ i k+1|k ) -[ 1 N N j=1 r τ j k (ξ j k+1|k ) ] ,
where

r τ i k (ξ i k+1|k ) = ρ τ i k k|k-1 + ξ τ i k k|k-1 σ 2 (ξ i k+1|k -a ξ τ i k k|k-1 ) ,
does not depend only on the position ξ i k+1|k , but on the position of the pair (ξ

τ i k k|k-1 , ξ i k+1|k ).
Statistics Plugging any joint particle approximation of the form (21) into equations ( 15) and ( 16), yields the following approximation

ℓ n ≈ 1 n n k=1 log[ 1 N N i=1 Ψ k (ξ i k|k-1 ) ] ,
for the log likelihood function, and the following approximation

∂ℓ n ∂θ ≈ 1 n n k=1       1 N N i=1 ρ i k|k-1 Ψ k (ξ i k|k-1 ) 1 N N i=1 Ψ k (ξ i k|k-1 )       = 1 n n k=1 [ N i=1 ρ i k|k-1 ω i k ] ,
for the score function, respectively. 22) and ( 23) read

Numerical results

X0

Q µ N k (dx ′ ) = N i=1 ω i k Q(ξ i k|k-1 , dx ′ ) = [ N i=1 ω i k q(ξ i k|k-1 , x ′ ) ] dx ′ ,
and

Q w N k (dx ′ ) + Γ µ N k (dx ′ ) = N i=1 [ ρ i k|k-1 -[ N j=1 ρ j k|k-1 ω j k ] + ξ i k|k-1 σ 2 (x ′ -a ξ i k|k-1 ) ] ω i k Q(ξ i k|k-1 , dx ′ ) = [ N i=1 [ ρ i k|k-1 -[ N j=1 ρ j k|k-1 ω j k ] + ξ i k|k-1 σ 2 (x ′ -a ξ i k|k-1 ) ] ω i k q(ξ i k|k-1 , x ′ ) ] dx ′ = r N k+1 (x ′ ) Q µ N k (dx ′ ) ,
where

r N k+1 (x ′ ) = N j=1 [ ρ j k|k-1 + ξ j k|k-1 σ 2 (x ′ -a ξ j k|k-1 ) ] ω j k q(ξ j k|k-1 , x ′ ) N j=1 ω j k q(ξ j k|k-1 , x ′ ) -[ N j=1 ρ j k|k-1 ω j k ] = N j=1 [ ρ j k|k-1 + ξ j k|k-1 σ 2 (x ′ -a ξ j k|k-1 ) ] ω j k exp{-1 2 (x ′ -a ξ j k|k-1 ) 2 } N j=1 ω j k exp{-1 2 (x ′ -a ξ j k|k-1 ) 2 } -[ N j=1 ρ j k|k-1 ω j k ] .
Resampling yields the particle approximation

Q µ N k ≈ µ N k+1|k = 1 N N i=1 δ ξ i k+1|k ,
where independently for any

i = 1, • • • , N ξ i k+1|k ∼ Q µ N k (dx ′ ) ,
which can be achieved for instance by taking

ξ i k ∼ µ N k (dx)
and

ξ i k+1|k ∼ Q( ξ i k , dx ′ ) ,
or even more explicitly by taking

τ i k ∼ (ω j k , j = 1, • • • , N ) , ξ i k = ξ τ i k k|k-1 and ξ i k+1|k ∼ Q( ξ i k , dx ′ ) ,
and the weighted particle approximation

Q w N k + Γ µ N k ≈ w N k+1|k = 1 N N i=1 r N k+1 (ξ i k+1|k ) δ ξ i k+1|k = 1 N N i=1 ρ i k+1|k δ ξ i k+1|k ,

Irisa

Particle based Methods for Parameter Estimation and Tracking 29 where the weight

ρ i k+1|k = r N k+1 (ξ i k+1|k ) = N j=1 [ ρ j k|k-1 + ξ j k|k-1 σ 2 (ξ i k+1|k -a ξ j k|k-1 ) ] ω j k exp{-1 2 (ξ i k+1|k -a ξ j k|k-1 ) } N j=1 ω j k exp{-1 2 (ξ i k+1|k -a ξ j k|k-1 ) 2 } -[ N j=1 ρ j k|k-1 ω j k ] ,
depends only on the position of ξ i k+1|k .

Numerical results

Part C : 2 dimensional case

We consider now the following model :

X k+1 = F X k + W k Y k = H X k + V k , (24)
where {W k , k ≥ 0} and {V k , k ≥ 0} are standard Gaussian independent white noise sequences, with covariance matrix Σ and S respectively. Here

F = f 11 f 12 f 21 f 22
is considered as an unknown parameter and X 0 ∼ N( X0 , Σ 0 ). It follows from the model that

X k+1 | X k = x ∼ N(F x, Σ) , i.e. Q(x, dx ′ ) = P[X k+1 ∈ dx ′ | X k = x] = 1 2 π √ det Σ exp{- 1 2 (x ′ -F x) * Σ -1 (x ′ -F x)} q(x, x ′ ) dx ′ ,
and

Y k | X k = x ∼ N(H x, S) , Irisa i.e. P[Y k ∈ dy | X k = x] = 1 2 π √ det S exp{- 1 2 (y -H x) * S -1 (y -H x)} g(x, y) dy ,
and let

Ψ k (x) = g(x, Y k ) ,
denote the likelihood function. The goal here is to design a particle approximation scheme for the derivative, with respect to the parameters f 11 , f 12 , f 21 , f 22 , of the conditional probability distribution of the hidden state, given the observations.

Exact expressions (Kalman lter)

In such a simple model, an exact expression can easily be obtained, via the prediction / correction steps of the Kalman lter framework. Indeed, the prediction step reads

X k | Y 1 , • • • , Y k-1 ∼ N( Xk|k-1 , Σ k|k-1 ) ,
i.e.

µ k|k-1 (dx) = P[X k ∈ dx | Y 1 , • • • , Y k-1 ] = 1 2 π det Σ k|k-1 exp{- 1 2 (x -Xk|k-1 ) * Σ -1 k|k-1 (x -Xk|k-1 )} p k|k-1 (x) dx ,
where

Xk|k-1 = F Xk-1
and

Σ k|k-1 = F Σ k-1 F * + Σ .
The log density of the r.v.

X k | Y 1 , • • • , Y k-1 is log p k|k-1 (x) = cste -1 2 log det Σ k|k-1 -1 2 (x -Xk|k-1 ) * Σ -1 k|k-1 (x -Xk|k-1 ) ,
hence the logarithmic derivative w.r.t. the parameter θ where θ ∈ {f 11 , f 12 , f 21 , f 22 }

∂ log p k|k-1 ∂θ (x) = 1 2   - ∂ det Σ k|k-1 ∂θ det Σ k|k-1 + ( ∂ Xk|k-1 ∂θ ) * Σ -1 k|k-1 (x -Xk|k-1 ) -(x -Xk|k-1 ) * ∂Σ -1 k|k-1 ∂θ (x -Xk|k-1 ) + (x -Xk|k-1 ) * Σ -1 k|k-1 ∂ Xk|k-1 ∂θ    ,
and

w k|k-1 (dx) = ∂µ k|k-1 ∂θ (dx) = ∂ log p k|k-1 ∂θ (x) µ k|k-1 (dx) ,
where

∂ Xk|k-1 ∂θ = ∂F ∂θ Xk-1 + F ∂ Xk-1 ∂θ , ∂Σ k|k-1 ∂θ = ∂F ∂θ Σ k-1 F * + F ∂Σ k-1 ∂θ F * + F Σ k-1 ( ∂F ∂θ ) * , ∂ det Σ k|k-1 ∂θ det Σ k|k-1 = trace[ ∂Σ k|k-1 ∂θ Σ -1 k|k-1 ] ,
PI n 1604 and

∂Σ -1 k|k-1 ∂θ = -Σ -1 k|k-1 ∂Σ k|k-1 ∂θ Σ -1 k|k-1 .
The last equality comes from the di erentiation of Σ k|k-1 Σ -1 k|k-1 = I. We have used the following result :

Proposition. If U is an invertible d × d matrix, then

∂ det U ∂θ det U = trace[ ∂U ∂θ U -1 ] . Proof. Writing U = ( u 1 , • • • , u d ), we have u 1 ∧ • • • ∧ u d = det U ( e 1 ∧ • • • ∧ e d ) ,
where

( e 1 , • • • , e d )
is the canonical basis. Di erentiating w.r.t. the parameter yields

∂ ∂θ ( u 1 ∧ • • • ∧ u d ) = ∂ det U ∂θ ( e 1 ∧ • • • ∧ e d ) = d i=1 u 1 ∧ • • • ∧ u i-1 ∧ ∂ u i ∂θ ∧ u i+1 ∧ • • • ∧ u d . Let ∂ u i ∂θ = n j=1 v ij u j , hence ∂ ∂θ ( u 1 ∧ • • • ∧ u d ) = d i,j=1 v ij ( u 1 ∧ • • • ∧ u i-1 ∧ u i ∧ u j ∧ u i+1 ∧ • • • ∧ u d = [ n i=1 v ii ] det U ( e 1 ∧ • • • ∧ e d ) ,
and it follows that

∂ det U ∂θ = [ n i=1 v ii ] det U .
Writing

u i = d k=1
u ij e k , and e k = d j=1 u kj u j , with U -1 = (u kj ) 1≤j,k≤d . we have

∂ u i ∂θ = d k=1 ∂u ik ∂θ e k = d j=1 [ d k=1 ∂u ik ∂θ u kj ] u j , hence v ij = d k=1 ∂u ij ∂θ u kj and n i=1 v ii = d i,k=1 ∂u ik ∂θ u ki = trace[ ∂U ∂θ U -1 ] .
The correction step reads

X k | Y 1 , • • • , Y k ∼ N( Xk , Σ k ) , Irisa i.e. µ k (dx) = P[X k ∈ dx | Y 1 , • • • , Y k ] = 1 2 π √ det Σ k exp{- 1 2 (x -Xk ) * Σ -1 k (x -Xk )} p k (x) dx ,
where the Kalman gain is de ned by

K k = Σ k|k-1 H * H Σ k|k-1 H * + Σ -1
and where

Xk = Xk|k-1 + K k (Y k -H Xk|k-1 )
and

Σ k = (I -K k H) Σ k|k-1 .
Therefore

∂ log p k ∂θ (x) = 1 2   - ∂ det Σ k ∂θ det Σ k + ( ∂ Xk ∂θ ) * Σ -1 k (x -Xk ) -(x -Xk ) * ∂Σ -1 k ∂θ (x -Xk ) + (x -Xk ) * Σ -1 k ∂ Xk ∂θ    ,
and

w k (dx) = ∂µ k ∂θ (dx) = ∂ log p k ∂θ (x) µ k (dx) ,
where

∂ Xk ∂θ = (I -K k H) ∂ Xk|k-1 ∂θ + ∂K k ∂θ (Y k -H Xk|k-1 ) , ∂Σ -1 k ∂θ = -Σ -1 k ∂Σ k ∂θ Σ -1 k , ∂ det Σ k ∂θ det Σ k = trace[ ∂Σ k ∂θ Σ -1 k ] , ∂Σ k ∂θ = (I -K k H) ∂Σ k|k-1 ∂θ ,
and

∂K k ∂θ = -[ ∂Σ k ∂θ Σ -1 k|k-1 + Σ k ∂Σ k|k-1 ∂θ ] H -1 .
The last equality is obvious since Σ k Σ -1 k|k-1 = I -K k H, and we have used the following result :

Proposition. If Q and R are two symmetric, square and positive de nite matrices, then

(H * R -1 H + Q -1 ) -1 = Q -Q H * (H * Q H + R) -1 H Q .

PI n 1604

Proof. We have

H * Q H + R ≥ R
and

H * R -1 H + Q -1 ≥ Q -1 hence H * Q H + R and H * R -1 H + Q -1 are invertible. It is now easy to check that [Q -Q H * (H * Q H + R) -1 H Q] [H * R -1 H + Q -1 ] = Q H * R -1 H + I -Q H * (H * Q H + R) -1 (H * Q H + R -R) R -1 H -Q H * (H * Q H + R) -1 H = I .
We have

Σ k = Σ k|k-1 -Σ k|k-1 H * (H Σ k|k-1 H * + S) -1 H Σ k|k-1 , hence Σ -1 k = H * S -1 H + Σ -1 k|k-1 ,
and

∂Σ -1 k ∂θ = ∂Σ -1 k|k-1 ∂θ = -Σ -1 k|k-1 ∂Σ k|k-1 ∂θ Σ -1 k|k-1 .
It follows from

Σ k Σ -1 k|k-1 = I -K k H , that ∂Σ k ∂θ = (I -K k H) ∂Σ k|k-1 ∂θ (I -K k H) * .
From the decomposition

Y k = H Xk|k-1 + I k ,
where the innovation

I k = Y k -H Xk|k-1 = H (X k -Xk|k-1 ) + V k ∼ N(0, S k ) with S k = H Σ k|k-1 H * + S , is independent of the past observations Y 1 , • • • , Y k-1 , it follows that Y k | Y 1 , • • • , Y k-1 ∼ N(H Xk|k-1 , S k ) ,
i.e.

P[Y k ∈ dy | Y 1 , • • • , Y k-1 ] = 1 2 π √ det S k exp{-1 2 (y -H Xk|k-1 ) * S -1 k (y -H Xk|k-1 )} g k (y)
dy .

Using the straightforward identity

P[Y 1 ∈ dy 1 , • • • , Y n ∈ dy n ] = n k=1 P[Y k ∈ dy k | Y 1 = y 1 , • • • , Y k-1 = y k-1 ] = n k=1 g k (y k ) dy k ,

Irisa

Particle based Methods for Parameter Estimation and Tracking 35 yields the following expression for the (suitably normalized) log likelihood function

ℓ n = 1 n log n k=1 g k (Y k ) = 1 n n k=1 log g k (Y k ) = cste - 1 n n k=1 [ 1 2 log det S k + 1 2 (y -H Xk|k-1 ) * S -1 k (y -H Xk|k-1 ) ] ,
and for the score function

∂ℓ n ∂θ = 1 2 n n k=1   - ∂ det S k ∂θ det S k + (H ∂ Xk|k-1 ∂θ ) * S -1 k (y -H Xk|k-1 ) -(y -H Xk|k-1 ) * ∂S -1 k ∂θ (y -H Xk|k-1 ) + (y -H Xk|k-1 ) * S -1 k H ∂ Xk|k-1 ∂θ    ,
where

∂S k ∂θ = H ∂Σ k|k-1 ∂θ H * and ∂S -1 k ∂θ = -S -1 k ∂S k ∂θ S -1 k .
Notice that in full generality

P[Y k ∈ dy | Y 1 , • • • , Y k-1 ] = ∞ -∞ P[Y k ∈ dy, X k ∈ dx | Y 1 , • • • , Y k-1 ] = ∞ -∞ P[Y k ∈ dy | X k = x] P[X k ∈ dx | Y 1 , • • • , Y k-1 ] = [ ∞ -∞ g(x, y) µ k|k-1 (dx) g k (y)
] dy , hence the following equivalent expression for the log likelihood function

ℓ n = 1 n n k=1 log g k (Y k ) = 1 n n k=1 log ∞ -∞ Ψ k (x) µ k|k-1 (dx) , (25)
and for the score function

∂ℓ n ∂θ = 1 n n k=1 ∞ -∞ Ψ k (x) ∂µ k|k-1 ∂θ (dx) ∞ -∞ Ψ k (x) µ k|k-1 (dx) = 1 n n k=1 ∞ -∞ Ψ k (x) w k|k-1 (dx) ∞ -∞ Ψ k (x) µ k|k-1 (dx)
.

(26)

Preliminary computations

We have to design a particle approximation scheme for the derivative w.r.t. f 11 , f 12 , f 21 , f 22 . The transition log density is

log q(x, x ′ ) = cste -1 2 log det Σ -1 2 (x ′ -F x) * Σ -1 (x ′ -F x) , with Σ -1 = Σ 11 Σ 12 Σ 12 Σ 22 , x = x 1 x 2 and x ′ = x ′ 1 x ′ 2
, hence the logarithmic derivative w.r.t. the parameter θ is

I(x, x ′ ) = ∂ log q ∂θ (x, x ′ ) = 1 2 [ ( ∂F ∂θ x) * Σ -1 (x ′ -F x) + (x ′ -F x) * Σ -1 ∂F ∂θ x ] ,
PI n 1604 and w.r.t. the parameters f 11 , f 12 , f 21 , f 22

∂ log q ∂f 11 (x, x ′ ) = x 1 Σ 11 (x ′ 1 -f 11 x 1 -f 12 x 2 ) + x 1 Σ 12 (x ′ 2 -f 21 x 1 -f 22 x 2 ) , ∂ log q ∂f 12 (x, x ′ ) = x 2 Σ 11 (x ′ 1 -f 11 x 1 -f 12 x 2 ) + x 2 Σ 12 (x ′ 2 -f 21 x 1 -f 22 x 2 ) , ∂ log q ∂f 21 (x, x ′ ) = x 1 Σ 12 (x ′ 1 -f 11 x 1 -f 12 x 2 ) + x 1 Σ 22 (x ′ 2 -f 21 x 1 -f 22 x 2 ) , ∂ log q ∂f 22 (x, x ′ ) = x 2 Σ 12 (x ′ 1 -f 11 x 1 -f 12 x 2 ) + x 2 Σ 22 (x ′ 2 -f 21 x 1 -f 22 x 2 ) ,
and

Γ(x, dx ′ ) = ∂Q ∂θ (x, dx ′ ) = I(x, x ′ ) Q(x, dx ′ )
where θ ∈ {f 11 , f 12 , f 21 , f 22 }. Di erentiating with respect to the parameter θ, throughout the recursions

µ k (dx) = (Ψ k • µ k|k-1 )(dx) = Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) , (27)
and

µ k+1|k (dx) = Q µ k (dx ′ ) = ∞ -∞ µ k (dx) Q(x, dx ′ ) , (28)
yields

w k (dx) = ∂µ k ∂θ (dx) = Ψ k (x) ∂µ k|k-1 ∂θ (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) - ∞ -∞ Ψ k (x ′ ) ∂µ k|k-1 ∂θ (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) , hence w k (dx) = (F k (µ k|k-1 )w k|k-1 )(dx) = Ψ k (x) w k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) - ∞ -∞ Ψ k (x ′ ) w k|k-1 (dx ′ ) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) , (29)
and

w k+1|k (dx ′ ) = Q w k (dx ′ ) + Γ µ k (dx ′ ) = ∞ -∞ [ w k (dx) + I(x, x ′ ) µ k (dx) ] Q(x, dx ′ ) . (30)
Since

w k|k-1 ≪ µ k|k-1 , it holds F k (µ k|k-1 )w k|k-1 ≪ Ψ k • µ k|k-1 , with Radon Nikodym derivative dw k dµ k (x) = d(F k (µ k|k-1 )w k|k-1 ) d(Ψ k • µ k|k-1 ) (x) = dw k|k-1 dµ k|k-1 (x) - ∞ -∞ dw k|k-1 dµ k|k-1 (x ′ ) (Ψ k • µ k|k-1 )(dx ′ ) ,
and

w k+1|k (dx ′ ) = Q w k (dx ′ ) + Γ µ k (dx ′ ) = ∞ -∞ [ dw k dµ k (x) + I(x, x ′ ) ] µ k (dx) Q(x, dx ′ ) .
Irisa 4.3.3 A particle approximation scheme Assuming that the following particle approximations

µ k|k-1 ≈ µ N k|k-1 = 1 N N i=1 δ ξ i k|k-1 and w k|k-1 ≈ w N k|k-1 = 1 N N i=1 ρ i k|k-1 δ ξ i k|k-1 , (31)
are available at time index k, and plugging these approximations into equations ( 27) and (29), yields

µ N k = Ψ k • µ N k|k-1 = N i=1 Ψ k (ξ i k|k-1 ) δ ξ i k|k-1 N j=1 Ψ k (ξ j k|k-1 ) = N i=1 ω i k δ ξ i k|k-1
, and

w N k = F k (µ N k|k-1 )w N k|k-1 = N i=1 ρ i k|k-1 Ψ k (ξ i k|k-1 ) δ ξ i k|k-1 N j=1 Ψ k (ξ j k|k-1 ) - N j=1 ρ j k|k-1 Ψ k (ξ j k|k-1 ) N j=1 Ψ k (ξ j k|k-1 ) N i=1 Ψ k (ξ i k|k-1 ) δ ξ i k|k-1 N j=1 Ψ k (ξ j k|k-1 ) = N i=1 ρ i k|k-1 ω i k δ ξ i k|k-1 -[ N j=1 ρ j k|k-1 ω j k ] N i=i ω i k δ ξ i k|k-1 = N i=1 [ ρ i k|k-1 -[ N j=1 ρ j k|k-1 ω j k ] ] ω i k δ ξ i k|k-1
.

Plugging these expressions into equations ( 28) and (30), yields

Q µ N k (dx ′ ) = N i=1 ω i k Q(ξ i k|k-1 , dx ′ ) , (32)
and

Q w N k (dx ′ ) + Γ µ N k (dx ′ ) = N i=1 [ ρ i k|k-1 -[ N j=1 ρ j k|k-1 ω j k ] + I(ξ i k|k-1 , x ′ ) ] ω i k Q(ξ i k|k-1 , dx ′ ) , (33)
which can be interpreted as the marginals of the nite signed measures

m = (m i , i = 1, • • • , N ) and s = (s i , i = 1, • • • , N ) de ned on the product space E N = {1, • • • , N } × R by s i (dx ′ ) = [ ρ i k|k-1 + I(ξ i k|k-1 , x ′ ) r i (x ′ ) -[ N j=1 ρ j k|k-1 ω j k ] ] ω i k Q(ξ i k|k-1 , dx ′ ) m i (dx ′ ) . Notice that N j=1 ∞ -∞ r j (x ′ ) m j (dx ′ ) = N j=1 ρ j k|k-1 ω j k ,
is just a normalizing constant. Introducing the particle approximation

m ≈ m N = 1 N N i=1 δ (τ i k , ξ i k+1|k )
, PI n 1604 and the weighted particle approximation

s ≈ s N = 1 N N i=1 [ r τ i k (ξ i k+1|k ) -[ 1 N N j=1 r τ j k (ξ j k+1|k ) ] ] δ (τ i k , ξ i k+1|k )
, where independently for any

i = 1, • • • , N , the pair (τ i k , ξ i k+1|k
) is jointly distributed according to the probability distribution m = (m i , i = 1, • • • , N ), i.e.

τ i k ∼ (ω j k , j = 1, • • • , N )
and

ξ i k+1|k ∼ Q(ξ τ i k k|k-1 , dx ′ ) ,
and marginalizing, yields

µ k+1|k ≈ µ N k+1|k = 1 N N i=1 δ ξ i k+1|k . w k+1|k ≈ w N k+1|k = 1 N N i=1 [ r τ i k (ξ i k+1|k ) -[ 1 N N j=1 r τ j k (ξ j k+1|k ) ] ] δ ξ i k+1|k = 1 N N i=1 ρ i k+1|k δ ξ i k+1|k .
Notice that the weight

ρ i k+1|k = r τ i k (ξ i k+1|k ) -[ 1 N N j=1 r τ j k (ξ j k+1|k ) ] ,
where

r τ i k (ξ i k+1|k ) = ρ τ i k k|k-1 + I(ξ τ i k k|k-1 , ξ i k+1|k ) ,
does not depend only on the position ξ i k+1|k , but on the position of the pair (ξ τ i k k|k-1 , ξ i k+1|k ).

Numerical results

X0

Q µ N k (dx ′ ) = N i=1 ω i k Q(ξ i k|k-1 , dx ′ ) = [ N i=1 ω i k q(ξ i k|k-1 , x ′ ) ] dx ′ ,
and

Q w N k (dx ′ ) + Γ µ N k (dx ′ ) = N i=1 [ ρ i k|k-1 -[ N j=1 ρ j k|k-1 ω j k ] + I(ξ i k|k-1 , x ′ ) ] ω i k Q(ξ i k|k-1 , dx ′ ) = [ N i=1 [ ρ i k|k-1 -[ N j=1 ρ j k|k-1 ω j k ] + I(ξ i k|k-1 , x ′ ) ] ω i k q(ξ i k|k-1 , x ′ ) ] dx ′ = r N k+1 (x ′ ) Q µ N k (dx ′ ) ,
where

r N k+1 (x ′ ) = N j=1 [ ρ j k|k-1 + I(ξ j k|k-1 , x ′ ) ] ω j k q(ξ j k|k-1 , x ′ ) N j=1 ω j k q(ξ j k|k-1 , x ′ ) -[ N j=1 ρ j k|k-1 ω j k ] = N j=1 [ ρ j k|k-1 + I(ξ j k|k-1 , x ′ ) ] ω j k exp{-1 2 (x ′ -F ξ j k|k-1 ) * Σ -1 (x ′ -F ξ j k|k-1 )} N j=1 ω j k exp{-1 2 (x ′ -F ξ j k|k-1 ) * Σ -1 (x ′ -F ξ j k|k-1 )} -[ N j=1 ρ j k|k-1 ω j k ] .
Resampling yields the particle approximation

Q µ N k ≈ µ N k+1|k = 1 N N i=1 δ ξ i k+1|k ,
where independently for any

i = 1, • • • , N ξ i k+1|k ∼ Q µ N k (dx ′ ) ,
which can be achieved for instance by taking

ξ i k ∼ µ N k (dx)
and

ξ i k+1|k ∼ Q( ξ i k , dx ′ ) ,
or even more explicitly by taking

τ i k ∼ (ω j k , j = 1, • • • , N ) , ξ i k = ξ τ i k k|k-1 and ξ i k+1|k ∼ Q( ξ i k , dx ′ ) ,
and the weighted particle approximation

Q w N k + Γ µ N k ≈ w N k+1|k = 1 N N i=1 r N k+1 (ξ i k+1|k ) δ ξ i k+1|k = 1 N N i=1 ρ i k+1|k δ ξ i k+1|k ,

Irisa

where the weight

ρ i k+1|k = r N k+1 (ξ i k+1|k ) = N j=1 [ ρ j k|k-1 + I(ξ j k|k-1 , ξ i k+1|k ) ] ω j k exp{-1 2 (ξ i k+1|k -F ξ j k|k-1 ) * Σ -1 (ξ i k+1|k -F ξ j k|k-1 )} N j=1 ω j k exp{-1 2 (ξ i k+1|k -F ξ j k|k-1 ) * Σ -1 (ξ i k+1|k -F ξ j k|k-1 )} -[ N j=1 ρ j k|k-1 ω j k ] ,
depends only on the position of ξ i k+1|k .

Numerical results In this last section, we consider a more general model :

X0 Σ 0 F Σ H S N (# particles)
X k+1 = F X k + W k Y k = H X k + V k , (34)
where {W k , k ≥ 0} and {V k , k ≥ 0} are standard Gaussian independent white noise sequences, with covariance matrix Σ and S respectively. Here

F =    f 11 • • • f 1m ′ . . . . . . . . . f m ′ 1 • • • f m ′ m ′   
and

H =    h 11 • • • h 1m • • • 0 . . . . . . . . . . . . . . . . . . h r1 • • • h rm • • • 0   
are considered as unknown parameters and X 0 ∼ N( X0 , Σ 0 ). As X is a vector with 2m rows, we set m ′ = 2m.

It follows from the model that

X k+1 | X k = x ∼ N(F x, Σ) , i.e. Q(x, dx ′ ) = P[X k+1 ∈ dx ′ | X k = x] = 1 (2 π) m ′ /2 √ det Σ exp{- 1 2 (x ′ -F x) * Σ -1 (x ′ -F x)} q(x, x ′ ) dx ′ ,
and

Y k | X k = x ∼ N(H x, S) ,
i.e.

P[Y k ∈ dy | X k = x] = 1 (2 π) r/2 √ det S exp{- 1 2 (y -H x) * S -1 (y -H x)} g(x, y) dy ,

Irisa

and let

Ψ k (x) = g(x, Y k ) ,
denote the likelihood function. The goal here is to design a particle approximation scheme for the derivative, with respect to the parameters In such a simple model, an exact expression can easily be obtained, via the prediction / correction steps of the Kalman lter framework. Indeed, the prediction step reads

f 11 , • • • , f m ′ m ′ , h
X k | Y 1 , • • • , Y k-1 ∼ N( Xk|k-1 , Σ k|k-1 ) ,
i.e.

µ k|k-1 (dx) = P[X k ∈ dx | Y 1 , • • • , Y k-1 ] = 1 (2 π) m ′ /2 det Σ k|k-1 exp{- 1 2 (x -Xk|k-1 ) * Σ -1 k|k-1 (x -Xk|k-1 )} p k|k-1 (x) dx ,
where

Xk|k-1 = F Xk-1
and

Σ k|k-1 = F Σ k-1 F * + Σ .
The log density of the r.v.

X k | Y 1 , • • • , Y k-1 is log p k|k-1 (x) = cste -1 2 log det Σ k|k-1 -1 2 (x -Xk|k-1 ) * Σ -1 k|k-1 (x -Xk|k-1 ) ,
hence the logarithmic derivative w.r.t. the parameter θ where θ ∈ {f 11 ,

• • • , f rm ′ , h 11 , • • • , h rm } ∂ log p k|k-1 ∂θ (x) = 1 2   - ∂ det Σ k|k-1 ∂θ det Σ k|k-1 + ( ∂ Xk|k-1 ∂θ ) * Σ -1 k|k-1 (x -Xk|k-1 ) -(x -Xk|k-1 ) * ∂Σ -1 k|k-1 ∂θ (x -Xk|k-1 ) + (x -Xk|k-1 ) * Σ -1 k|k-1 ∂ Xk|k-1 ∂θ    ,
and

w k|k-1 (dx) = ∂µ k|k-1 ∂θ (dx) = ∂ log p k|k-1 ∂θ (x) µ k|k-1 (dx) ,
where

∂ Xk|k-1 ∂θ = ∂F ∂θ Xk-1 + F ∂ Xk-1 ∂θ , ∂Σ k|k-1 ∂θ = ∂F ∂θ Σ k-1 F * + F ∂Σ k-1 ∂θ F * + F Σ k-1 ( ∂F ∂θ ) * , ∂ det Σ k|k-1 ∂θ det Σ k|k-1 = trace [ ∂Σ k|k-1 ∂θ Σ -1 k|k-1 ] ,
PI n 1604 and

∂Σ -1 k|k-1 ∂θ = -Σ -1 k|k-1 ∂Σ k|k-1 ∂θ Σ -1 k|k-1 .
The last equality comes from the di erentiation of Σ k|k-1 Σ -1 k|k-1 = I.

The correction step reads

X k | Y 1 , • • • , Y k ∼ N( Xk , Σ k ) ,
i.e.

µ k (dx) = P[X k ∈ dx | Y 1 , • • • , Y k ] = 1 (2 π ) m ′ /2 √ det Σ k exp{- 1 2 (x -Xk ) * Σ -1 k (x -Xk )} p k (x) dx ,
where the Kalman gain is de ned by

K k = Σ k|k-1 H * H Σ k|k-1 H * + Σ -1
and where

Xk = Xk|k-1 + K k (Y k -H Xk|k-1 )
and

Σ k = (I -K k H) Σ k|k-1 . Therefore ∂ log p k ∂θ (x) = 1 2   - ∂ det Σ k ∂θ det Σ k + ( ∂ Xk ∂θ ) * Σ -1 k (x -Xk ) -(x -Xk ) * ∂Σ -1 k ∂θ (x -Xk ) + (x -Xk ) * Σ -1 k ∂ Xk ∂θ    ,
and

w k (dx) = ∂µ k ∂θ (dx) = ∂ log p k ∂θ (x) µ k (dx) ,
where

∂ Xk ∂θ = (I -K k H) ∂ Xk|k-1 ∂θ + ∂K k ∂θ (Y k -H Xk|k-1 ) -K k ∂H ∂θ Xk|k-1 , ∂Σ -1 k ∂θ = -Σ -1 k ∂Σ k ∂θ Σ -1 k , ∂ det Σ k ∂θ det Σ k = trace [ ∂Σ k ∂θ Σ -1 k ] , ∂Σ k ∂θ = (I -K k H) ∂Σ k|k-1 ∂θ -[ ∂K k ∂θ H + K k ∂H ∂θ ] Σ k|k-1 ,

Irisa

and

∂K k ∂θ = Σ k|k-1 ( ∂H ∂θ ) * [H Σ k|k-1 H * + Σ] -1 + ∂Σ k|k-1 ∂θ H * [H Σ k|k-1 H * + Σ] -1 -Σ k|k-1 H * [H Σ k|k-1 H * + Σ] -1 ] ∂H ∂θ Σ k|k-1 H * + H ∂Σ k|k-1 ∂θ H * + H Σ k|k-1 ( ∂H ∂θ ) * [H Σ k|k-1 H * + Σ] -1 .
From the decomposition

Y k = H Xk|k-1 + I k ,
where the innovation

I k = Y k -H Xk|k-1 = H (X k -Xk|k-1 ) + V k ∼ N(0, S k ) with S k = H Σ k|k-1 H * + S , is independent of the past observations Y 1 , • • • , Y k-1 , it follows that Y k | Y 1 , • • • , Y k-1 ∼ N(H Xk|k-1 , S k ) , i.e. P[Y k ∈ dy | Y 1 , • • • , Y k-1 ] = 1 (2 π) r/2 √ det S k exp{-1 2 (y -H Xk|k-1 ) * S -1 k (y -H Xk|k-1 )} g k (y)
dy .

Using the straightforward identity

P[Y 1 ∈ dy 1 , • • • , Y k ∈ dy k ] = n k=1 P[Y k ∈ dy k | Y 1 = y 1 , • • • , Y k-1 = y k-1 ] = n k=1 g k (y k ) dy k ,
yields the following expression for the (suitably normalized) log likelihood function

ℓ n = 1 n log n k=1 g k (Y k ) = 1 n n k=1 log g k (Y k ) = cste - 1 2n n k=1 [ log(det S k ) + (y -H Xk|k-1 ) * S -1 k (y -H Xk|k-1 ) ] ,
and for the score function

∂ℓ n ∂θ = 1 2n n k=1   - ∂ det S k ∂θ det S k + (H ∂ Xk|k-1 ∂θ + ∂H ∂θ Xk|k-1 ) * S -1 k (y -H Xk|k-1 ) -(y -H Xk|k-1 ) * ∂S -1 k ∂θ (y -H Xk|k-1 ) + (y -H Xk|k-1 ) * S -1 k (H ∂ Xk|k-1 ∂θ + ∂H ∂θ Xk|k-1 )    ,
where

∂S k ∂θ = ∂H ∂θ Σ k|k-1 H * + H ∂Σ k|k-1 ∂θ H * + H Σ k|k-1 ( ∂H ∂θ ) * and ∂S -1 k ∂θ = -S -1 k ∂S k ∂θ S -1 k .

PI n 1604

Notice that in full generality

P[Y k ∈ dy | Y 1 , • • • , Y k-1 ] = ∞ -∞ P[Y k ∈ dy, X k ∈ dx | Y 1 , • • • , Y k-1 ] = ∞ -∞ P[Y k ∈ dy | X k = x] P[X k ∈ dx | Y 1 , • • • , Y k-1 ] = [ ∞ -∞ g(x, y) µ k|k-1 (dx) g k (y)
] dy , hence the following equivalent expression for the log likelihood function

ℓ n = 1 n n k=1 log g k (Y k ) = 1 n n k=1 log ∞ -∞ Ψ k (x) µ k|k-1 (dx) , (35)
and for the score function

∂ℓ n ∂θ = 1 n n k=1 ∞ -∞ Ψ k (x) ∂µ k|k-1 ∂θ (dx) + ∂Ψ k (x) ∂θ µ k|k-1 (dx) ∞ -∞ Ψ k (x) µ k|k-1 (dx) = 1 n n k=1 ∞ -∞ Ψ k (x) w k|k-1 (dx) + ∂Ψ k (x) ∂θ µ k|k-1 (dx) ∞ -∞ Ψ k (x) µ k|k-1 (dx)
.

Preliminary computations

We have to design a particle approximation scheme for the derivative w.r.t.

f 11 , • • • , f m ′ m ′ , h 11 , • • • , h rm .
The transition log density is

log q(x, x ′ ) = cste -1 2 log(det Σ) -1 2 (x ′ -F x) * Σ -1 (x ′ -F x) ,
hence the logarithmic derivative w.r.t. the parameter θ is

I(x, x ′ ) = ∂ log q ∂θ (x, x ′ ) = 1 2 [ ( ∂F ∂θ x) * Σ -1 (x ′ -F x) + (x ′ -F x) * Σ -1 ∂F ∂θ x ] ,
and

Γ(x, dx ′ ) = ∂Q ∂θ (x, dx ′ ) = I(x, x ′ ) Q(x, dx ′ ) where θ ∈ {f 11 , • • • , f m ′ m ′ , h 11 , • • • , h rm }.
Di erentiating with respect to the parameter θ, throughout the recursions

µ k (dx) = (Ψ k • µ k|k-1 )(dx) = Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) , (36)
and

µ k+1|k (dx) = (Q µ k )(dx ′ ) = ∞ -∞ µ k (dx) Q(x, dx ′ ) , (37)
Irisa yields

w k (dx) = ∂µ k ∂θ (dx) = Ψ k (x) ∂µ k|k-1 ∂θ (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) - ∞ -∞ Ψ k (x ′ ) ∂µ k|k-1 ∂θ (dx ′ ) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) + ∂ log Ψ k ∂θ (x) Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) - ∞ -∞ ∂ log Ψ k ∂θ (x ′ ) Ψ k (x ′ ) µ k|k-1 (dx ′ ) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) , hence, with S k (x) = ∂ log Ψ k ∂θ (x) w k (dx) = (F k (µ k|k-1 )w k|k-1 )(dx) + (G k (µ k|k-1 ))(dx) = Ψ k (x) w k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) - ∞ -∞ Ψ k (x ′ ) w k|k-1 (dx ′ ) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) + S k (x) Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) - ∞ -∞ S k (x ′ ) Ψ k (x ′ ) µ k|k-1 (dx ′ ) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) Ψ k (x) µ k|k-1 (dx) ∞ -∞ Ψ k (x ′ ) µ k|k-1 (dx ′ ) . and w k+1|k (dx ′ ) = Q w k (dx ′ ) + Γ µ k (dx ′ ) = ∞ -∞ [ w k (dx) + I(x, x ′ ) µ k (dx) ] Q(x, dx ′ ) . (38)
Since

w k|k-1 ≪ µ k|k-1 , it holds F k (µ k|k-1 )w k|k-1 ≪ Ψ k • µ k|k-1 , and G k (µ k|k-1 ) ≪ Ψ k • µ k|k-1 with Radon Nikodym derivative dw k dµ k (x) = dw k|k-1 dµ k|k-1 (x) + S k (x) - ∞ -∞ [ dw k|k-1 dµ k|k-1 (x ′ ) + S k (x ′ ) ] (Ψ k • µ k|k-1 )(dx ′ ) ,
and

w k+1|k (dx ′ ) = Q w k (dx ′ ) + Γ µ k (dx ′ ) = ∞ -∞ [ dw k dµ k (x) + I(x, x ′ ) ] µ k (dx) Q(x, dx ′ ) .

A particle approximation scheme

Assuming that the following particle approximations

µ k|k-1 ≈ µ N k|k-1 = 1 N N i=1 δ ξ i k|k-1 and w k|k-1 ≈ w N k|k-1 = 1 N N i=1 ρ i k|k-1 δ ξ i k|k-1 , (39)
are available at time index k, and plugging these approximations into equations ( 36) and (38), yields

µ N k = Ψ k • µ N k|k-1 = N i=1 Ψ k (ξ i k|k-1 ) δ ξ i k|k-1 N j=1 Ψ k (ξ j k|k-1 ) = N i=1 ω i k δ ξ i k|k-1 , PI n 1604
(iv) the fourth graphic shows the exact derivative of the density w.r.t. the parameter a, displayed in blue solid line, and an histogram of the weighted particle approximation associated with the weighted particle system {(ξ i k|k-1 , ρ i k|k-1 ) , i = 1, • • • , N }, i.e. for each subinterval A of length |A|, a black solid bar of (either positive or negative) height

1 |A| [ 1 N N i=1 ρ i k|k-1 1 (ξ i k|k-1 ∈ A) ] ,
is drawn (either upwards or downwards) at the center of the interval A.

Even if the second graphic shows a good match between the blue solid line and the black histogram bars, which implies in particular a good match between the region where the green solid line lies above the zero axis and the region where the red solid bars are drawn upwards, it could nevertheless happen that some particles with a negative weight fall within this region, at the expense of having their contribution compensated by the contribution of particles with a positive weight, in such a way that the balance over each subinterval of this region has the right sign and magnitude : to detect this situation, which denotes an une cient allocation of signed weights to particles, the third graphic shows the exact derivative of the density w.r.t. the parameter a, displayed in blue solid line, and separately an histogram of the positive and negative parts of the weighted particle approximation associated with the weighted particle system {(ξ i k|k-1 , ρ i k|k-1 ) , i = 1, • • • , N }, i.e. for each subinterval A of length |A|, a red solid bar of positive height

1 |A| [ 1 N N i=1 ρ i k|k-1 1 (ρ i k|k-1 ≥ 0) 1 (ξ i k|k-1 ∈ A) ] ,
is drawn upwards at the center of the interval A, and a green solid bar of negative height

1 |A| [ 1 N N i=1 ρ i k|k-1 1 (ρ i k|k-1 ≤ 0) 1 (ξ i k|k-1 ∈ A) ] ,
is drawn downwards at the center of the interval A, (v) the fth graphic shows the weighted particle system {(ξ i k|k-1 , ρ i k|k-1 ) , i = 1, • • • , N }, i.e. an (either green or red) solid bar of (either positive or negative) height ρ i k|k-1 , is drawn (either upwards or downwards) at each particle location ξ i k|k-1 . In Part C, we have graphics of the score function for each parameter and the log likelihood function with the same conventions as above.

Conclusions

In this section, we have seen two di erent methods to approximate the derivative of the lter w.r.t. the parameter. The second particle approximation scheme provides better results with less particles but the rst scheme is faster and very easy to compute. As our goal is to design a particle approximation scheme to study real data with a high sampling rate, we will use the rst scheme.

Moreover, the simulations show that the deterministic resampling step described before leads to very good approximations of the optimal log likelihood and score functions.

Particle implementation of the RML algorithm

The algorithm presented in this section is a recursive version of the gradient algorithm usually referred to as recursive maximum likelihood (RML). It requires essentially the ability to compute the optimal lter and the derivative of this lter with respect to the parameter of interest. We apply our method to approximate the derivative of the lter to perform RML. If the parameter θ is actually time varying but one does not have a dynamic model for its evolution, a standard approach to track this parameter consists of using the recursive algorithm presented in Section 5.1 using a xed step size γ instead of a decreasing sequence γ k . Selecting the step size is a di cult problem. If γ is too large, the statistical uctuations around the parameter are too large. If γ is too small, the algorithm loses its tracking ability. In this section, the true parameter symbolized by a solid blue line is approximated by the estimator θ k provided by the RML in solid black line and the averaged estimator θ k in red blue line.

Numerical results

We have used a reprojection method in order to ensure the convergence of our algorithm. In fact, the parameter θ = (θ 1 , • • • , θ k ) ∈ Θ, where Θ is a compact subset of R k . The updating step does not ensure that θ k+1 ∈ Θ even if θ k ∈ Θ. A standard approach in stochastic approximation to prevent divergence consists of reprojecting θ k+1 inside Θ whenever the value obtained, say θk+1 , does not belong to Θ. One has Θ = k i=1 [θ i,min , θ i,max ] and the reprojection procedure simply consists of setting θ i,k+1 = θ i,min if θ i,k+1 < θ i,min and θ i,k+1 = θ i,max if θ i,k+1 > θ i,max . For further results and asymptotic theory, see 8].

Conclusions

We have seen the e ciency of particle approximation through several examples. Indeed, we are able to estimate and track a multidimensional parameter using good approximations of optimal lters and its derivatives. In the estimation problem, the true parameter is well approximated by the estimator using particle methods even in the two dimensional case where the four coe cients of a square matrix are estimated at the same time. The tracking problem has also been studied, with good results.

4

  Particle approximation of nite signed measures, with numerical experiments

4

 4 An alternate particle approximation scheme Equations (

4

 4 An alternate particle approximation scheme Equations (32) and (33) read

2

 2 Particle lter and parameter tracking

  for Parameter Estimation and Tracking 57

  [START_REF] Douc | Asymptotics of the maximum likelihood estimator for general hidden Markov models[END_REF]

  , • • • , h rm , of the conditional probability distribution of the hidden state, given the observations.

	4.4.1 Exact expressions (Kalman lter)

PI n 1604

and

S k (ξ j k|k-1 ) Ψ k (ξ j k|k-1 )

(ρ j k|k-1 + S k (ξ j k|k-1 ))

.

Plugging these expressions into equations ( 37) and (38), yields

and

which can be interpreted as the marginals of the nite signed measures m = (m i , i = 1, • • • , N ) and s = (s i , i = 1, • • • , N ) de ned on the product space

is just a normalizing constant. Introducing the particle approximation

, and the weighted particle approximation

,

Irisa

Particle based Methods for Parameter Estimation and Tracking

49

where independently for any

and

and marginalizing, yields

Notice that the weight

where

does not depend only on the position ξ i k+1|k .

An alternate particle approximation scheme

Equations ( 40) and (41) read

and

where

PI n 1604

Resampling yields the particle approximation

where independently for any

which can be achieved for instance by taking

and

or even more explicitly by taking

and the weighted particle approximation

where the weight

depends only on the position of ξ i k+1 .

About the graphical outputs

The gures presented in Part B, where the HMM situation is considered, correspond to the particle approximation of the prediction density and its derivative w.r.t. the parameter a, and to the re sampled particle approximation of the correction density and its derivative w.r.t. the parameter a.

(i) the rst graphic shows the exact trajectory of the log likelihood function, displayed in black solid line, as provided by the Kalman lter, and the approximate trajectory, displayed in blue solid line, as provided by the particle lter approximation, (ii) similarly, the second graphic shows the exact trajectory of the score function, displayed in black solid line, as provided by the Kalman lter, and the approximate trajectory, displayed in blue solid line, as provided by the particle lter approximation, (iii) the third graphic shows the exact density, displayed in blue solid line, and an histogram of the particle approximation associated with the particle system {ξ i k|k-1 , i = 1, • • • , N }, i.e. for each subinterval A of length |A|, a black solid bar of height

is drawn upwards at the center of the interval A, Irisa

Particle lter and parameter estimation

In this section, the parameter is denoted by θ and dependence w.r.t. the parameter appears explicitly in the notation for the transition kernel Q θ (x, dx ′ ), and for the linear tangent kernel K θ (x, dx ′ , ds ′ ). It is well known that in such a parametric model, the log likelihood function for the estimation of the parameter θ can be written as

and the score function, i.e. the derivative of the log likelihood function w.r.t. the parameter, can be written as

where the lter {µ θ k , k ≥ 0} and the linear tangent lter {w θ k k ≥ 0} satisfy

and

respectively. It is natural to consider the RML algorithm to identify the parametric model, which is de ned by the following relation

where typically γ k ≃ k -2/3 , and the averaged estimator (which achieves the minimum variance of the estimation error) is obtained by post processing

i.e.

where the adaptive lter { µ k , k ≥ 0} and the adaptive linear tangent lter { w k , , k ≥ 0} satisfy the same equations as the lter and the linear tangent lter respectively, in which the value of the parameter is adapted at each time instant according to equation (42), i.e.

and

Irisa

Particle based Methods for Parameter Estimation and Tracking 53 respectively. The particle implementation of the RML equation ( 42) is

and

and the corresponding algorithm is described below.

Particle implementation of the RML algorithm

end for

parameter update]

The mathematical analysis of the asymptotic properties of the estimator θ N k as k → ∞ and N → ∞ is far beyond the scope of this work, and would rely on joint stability properties of the lter and the linear tangent lter, which is a very di cult question. Even the asymptotic properties of the estimator θ k as k → ∞ are di cult to prove, unless some mixing assumption holds for the transition kernels Q θ (x, dx ′ ) and the linear tangent kernels Γ θ (x, dx ′ ), which practically implies that the state space E should be compact, see e.g. 10] where only the non recursive ML is studied.