Filtrage Bayésien et Approximation Particulaire

140 Pages·2008·0.6 MB·French

´
Ecole Nationale Sup´erieure de Techniques Avanc´ees
Module B7 : Commande des Syst`emes
Filtrage Bay´esien
et Approximation Particulaire
version du 7 septembre 2008
Fran¸cois Le Gland
IRISA / INRIA Rennes
http://www.irisa.fr/aspi/legland/ensta/

i
Objectif du cours
Le filtrage consiste `a estimer de fa¸con r´ecursive un ´etat cach´e au vu d’observations. Le
domaine d’application principal est la localisation, la navigation et la poursuite de mobiles, dans
ledomainemilitaire,maisaussienrobotiquemobile,envisionparordinateur,encommunications
sans–fil (GSM en ext´erieur, WiFi en indoor), ou` il s’agit de combiner : un mod`ele a priori de
d´eplacement du mobile, des mesures issues de capteurs, et´eventuellemnent une base de mesures
de r´ef´erences, disponibles par exemples sous la forme d’une carte num´erique (mod`ele num´erique
de terrain, carte de couverture, etc.).
Le probl`eme de filtrage poss`ede une solution explicite, appel´ee filtre de Kalman, dans le cas
particulier des syst`emes lin´eaires gaussiens. Dans le cas plus g´en´eral des mod`eles de Markov
cach´es, des m´ethodes de simulations efficaces sont apparues r´ecemment, sous le nom de filtrage
particulaire. L’objectif de ce cours est de pr´esenter diff´erents algorithmes de filtrage particulaire,
de les mettre en œuvre dans le cadre de travaux pratiques, et de d´emontrer quelques r´esultats de
convergence en utilisant le cadre tr`es g´en´eral de l’approximation particulaire des distributions
de Feynman–Kac.
ii
Table des mati`eres
1 Introduction 1
1.1 Importance de l’information a priori . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Estimation bay´esienne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Cadre gaussien . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2 Syst`emes lin´eaires gaussiens 15
3 Filtrage de Kalman 19
3.1 Filtre de Kalman . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2 Filtre de Kalman lin´earis´e, filtre de Kalman ´etendu . . . . . . . . . . . . . . . . . 24
4 Au del`a des syst`emes lin´eaires gaussiens 29
4.1 Syst`emes non–lin´eaires `a bruits non–gaussiens . . . . . . . . . . . . . . . . . . . . 29
4.2 Mod`eles de Markov cach´es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.3 Chaˆınes de Markov `a param`etres markoviens . . . . . . . . . . . . . . . . . . . . 33
4.4 Chaˆınes de Markov partiellement observ´ees . . . . . . . . . . . . . . . . . . . . . 35
5 Borne de Cram´er–Rao a posteriori 39
6 Filtrage bay´esien 45
6.1 Mod`eles de Markov cach´es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
6.2 Chaˆınes de Markov partiellement observ´ees . . . . . . . . . . . . . . . . . . . . . 50
7 G´en´eralisation : distributions de Feynman–Kac 55
7.1 Mod`ele de base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
7.2 Mod`ele (apparamment) plus g´en´eral . . . . . . . . . . . . . . . . . . . . . . . . . 57
7.3 Mod`ele `a valeurs transitions ou trajectoires . . . . . . . . . . . . . . . . . . . . . 60
iii
iv TABLE DES MATIE`RES
8 M´ethodes de Monte Carlo 69
8.1 E´chantillonnage pond´er´e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
8.2 Simulation selon une distribution de Gibbs–Boltzmann . . . . . . . . . . . . . . . 74
8.2.A Acceptation / rejet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
8.2.B E´chantillonnage pond´er´e . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
8.3 Simulation selon un m´elange fini . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
8.4 E´chantillonnage selon une distribution `a support fini . . . . . . . . . . . . . . . . 90
9 Approximations particulaires 95
9.1 E´chantillonnage pond´er´e (SIS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
9.2 E´chantillonnage / r´e–´echantillonnage (SIR) . . . . . . . . . . . . . . . . . . . . . 98
9.3 R´e–´echantillonnage adaptatif . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
10 Estimation d’erreur 107
11 TCL pour les approximations particulaires 111
11.1 E´chantillonnage pond´er´e (SIS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
11.2 E´chantillonnage / r´e–´echantillonnage (SIR) . . . . . . . . . . . . . . . . . . . . . 112
A Inversion matricielle 119
B In´egalit´es 123
C TCL pour les tableaux triangulaires 127
Chapitre 1
Introduction
En toute g´en´eralit´e, le filtrage consiste `a estimer l’´etat d’un syst`eme dynamique, c’est–`a–dire
´evoluant au cours du temps, `a partir d’observations partielles, g´en´eralement bruit´ees.
Typiquement, on dispose d’une suite Y ,Y ,··· ,Y d’observations, par exemple obtenues
0 1 n
apr`es traitement pr´ealable du signal recueilli au niveau des capteurs. Chaque observation Y est
n
reli´ee `a l’´etat inconnu X par une relation du type
n
Y = h (X )+V , (1.1)
n n n n
ou` V est un bruit, qui mod´elise l’erreur d’observation. On pr´ecisera plus loin dans ce cours la
n
notion de bruit, en terme de variables al´atoires le plus souvent centr´ees (de moyenne nulle).
1.1 Importance de l’information a priori
Une hypoth`ese assez commune est de supposer que les variables al´eatoires V ,V ,··· ,V
0 1 n
sont ind´ependantes entre elles. A cause de cette hypoth`ese d’ind´ependance mutuelle des bruits
d’observation, et a fortiori en absence de bruit, seule l’observation Y participe `a l’estimation
n
de l’´etat cach´e X , c’est–`a–dire qu’on se trouve confront´e `a une succession de probl`emes d’esti-
n
mation d´ecoupl´es : dans la relation (1.1), l’observation Y est disponible (par d´efinition) tandis
n
que ni l’´etat cach´e X ni le bruit V ne sont disponibles, et il faut arriver `a retrouver (estimer)
n n
l’´etat cach´e X au vu de l’observation Y et malgr´e la pr´esence du bruit V .
n n n
Tel qu’il est formul´e, le probl`eme de l’estimation de l’´etat cach´e X `a partir des observations
n
Y ,Y ,··· ,Y est en g´en´eral mal–pos´e :
0 1 n
• en g´en´eral, la dimension m de la variable cach´ee est plus grande que la dimension d de
l’observation : mˆeme en absence de bruit, on ne peut pas inverser la relation (1.1) qui
poss`ede plus d’inconnues que d’´equations,
• dans le cas favorable ou` m = d, et mˆeme en absence de bruit, il n’est pas toujours possible
d’inverser la relation (1.1) qui peut tr`es bien poss´eder plusieurs solutions distinctes,
1
2 CHAPITRE 1. INTRODUCTION
• la situation est ´evidemment encore plus compliqu´ee en pr´esence de bruit : `a cause du
ph´enom`ene de d´ecouplage cit´e plus haut, la suite X ,X ,··· ,X reconstitu´ee peut ne pas
0 1 n
ˆetre pertinente en tant que trajectoire, mˆeme si chacune des estimations est pertinente
s´epar´ement.
Pour lever l’ind´etermination, c’est–`a–dire pour garantir l’existence d’une solution unique, et
pour r´esoudre le probl`eme de coh´erence temporelle, la solution classique consiste `a utiliser des
informations suppl´ementaires sur la suite cach´ee, par exemple sous la forme de fonctions de
couˆt portant sur l’´etat initial ou sur les transitions entre deux ´etats successifs. Par exemple, on
cherchera `a minimiser le crit`ere
n n
X X
J(x ) = c (x )+ c (x ,x )+ d (x ) ,
0:n 0 0 k k−1 k k k
k=1 k=0
par rapport `a la suite x = (x ,x ,··· ,x ), qui combine des fonctions de couˆt repr´esentant
0;n 0 1 n
l’information a priori sur la solution avec des fonctions de couˆt d’une autre nature, qui peuvent
repr´esenter par exemple un terme d’attache aux donn´ees, de la forme
h (x) = 1 |Y −h (x)|2 ou bien h (x) = 1 (Y −h (x))∗I (Y −h (x)) ,
k 2 k k k 2 k k k k k
pour tout k = 0,1,··· ,n, avec l’intrepr´etation que la suite recherch´ee doit ´egalement v´erifier `a
chaqueinstantl’´equationd’observationenunsensapproch´e.Plusg´en´eralement,cesfonctionsde
couˆt peuvent juste repr´esenter une contrainte (ou une propri´et´e) que la suite recherch´ee devrait
v´erifier (ou poss´eder). En absence d’information a priori, le crit`ere se r´eduit simplement `a
n n
X X
J(x ) = 1 |Y −h (x )|2 ou bien J(x ) = 1 (Y −h (x ))∗I (Y −h (x )) ,
0:n 2 k k k 0:n 2 k k k k k k k
k=0 k=0
ce qui revient en absence de couplage `a minimiser s´epar´ement le crit`ere
1 |Y −h (x )|2 ou bien 1 (Y −h (x ))∗R−1(Y −h (x )) ,
2 k k k 2 k k k k k k k
par rapport `a l’´etat x , pour tout k = 0,1,··· ,n, avec les cons´equences d´ej`a ´evoqu´ees en terme
k
d’ind´etermination et de possible incoh´erence temporelle. Un exemple classique de fonctions de
couˆt repr´esentant l’information a priori est
c (x) = 1 |x−µ|2 ou bien c (x) = 1 (x−µ)∗Σ−1(x−µ) ,
0 2 0 2 0
avec l’interpr´etation que l’´etat initial x recherch´e doit ˆetre proche de µ, et
0
c (x,x0) = 1 |x0−f (x)|2 ou bien c (x,x0) = 1 (x0−f (x))∗Q−1(x0−f (x)) ,
k 2 k k 2 k k k
avec l’interpr´etation que l’´etat x recherch´e doit ˆetre proche de f (x ), ou de mani`ere ´equiva-
k k k−1
lente que la transition (x ,x ) recherch´ee doit v´erifier l’´equation x = f (x ) dans un sens
k−1 k k k k−1
approch´e, pour tout k = 1,··· ,n. On remarque que ces fonctions de couˆt sont (`a une constante
additive pr`es) de la forme
c (x) = −logp (x) et c (x,x0) = −logp (x0 | x) , (1.2)
0 0 k k
1.1. IMPORTANCE DE L’INFORMATION A PRIORI 3
pour tout k = 1,··· ,n, ou` p (x) est la densit´e de probabilit´e initiale, et ou` p (x0 | x) est la
0 k
densit´e de probabilit´e de transition, dans le mod`ele non–lin´eaire suivant avec bruits gaussiens
additifs
X = f (X )+W avec W ∼ N(0,Q ) ,
k k k−1 k k k
etavecconditioninitialeX ∼ N(µ,Σ).Eneffet(`auneconstantemutiplicativedenormalisation
0
pr`es)
P[X ∈ dx] ∝ exp{−1 (x−µ)∗Σ−1(x−µ)}dx ∝ p (x)dx ,
0 2 0 0
et
P[X ∈ dx0 | X = x] ∝ exp{−1 (x0−f (x))∗Q−1(x0−f (x))}dx0 ∝ p (x0 | x)dx0 ,
k k−1 2 k k k k
pour tout k = 1,··· ,n. En toute g´en´eralit´e, si les relations (1.2) sont v´erifi´ees pour une densit´e
de probabilit´e p (x) et pour des densit´es de probabilit´e de transition p (x0 | x), pour tout
0 k
k = 1,··· ,n, alors le crit`ere `a minimiser peut s’´ecrire
n n
X X
J(x ) = −logp (x )− logp (x | x )+ d (x ) ,
0:n 0 0 k k k−1 k k
k=1 k=0
ce qui revient `a maximiser
n n
Y X
exp{−J(x )} = p (x ) p (x | x ) exp{− d (x )} ,
0:n 0 0 k k k−1 k k
k=1 k=0
| {z }
p (x )
0:n 0:n
par rapport `a la suite x = (x ,x ,··· ,x ). On remarque que p (x ) repr´esente la densit´e
0:n 0 1 n 0:n 0:n
deprobabilit´econjointedes´etatssuccessifs(X ,X ,··· ,X )delachaˆınedeMarkovcaract´eris´ee
0 1 n
par
• la densit´e de probabilit´e initiale p (x ),
0 0
• et les densit´es de probabilit´e de transition p (x0 | x), pour tout k = 1,··· ,n.
k
Comme alternative au point de vue de l’optimisation d´eterministe d´evelopp´e jusqu’ici, on adop-
tera dans ce cours un point de vue d’estimation bay´esienne, c’est–`a–dire qu’on remplacera le
probl`eme de minimisation d´eterministe, avec prise en compte de l’information a priori en terme
de fonctions de couˆt, par le probl`eme du calcul de la distribution de Gibbs–Boltzmann d´efinie
(`a une constante multiplicative pr`es) sur l’espace des trajectoires E = E ×···×E par
n
n n
Y X
exp{−J(x )} dx = p (x ) p (x | x ) exp{− d (x )} dx . (1.3)
0:n 0:n 0 0 k k k−1 k k 0:n
k=1 k=0
| {z }
p (x )
0:n 0:n
En d’autres termes, on remplacera le probl`eme de calculer le mode, c’est–`a–dire la trajectoire
x = (x ,x ,··· ,x ) de plus forte densit´e, par le probl`eme de calculer des esp´erances (ou des
0:n 0 1 n
4 CHAPITRE 1. INTRODUCTION
int´egrales) du type
Z Z Z Z n
X
··· f(x ) exp{−J(x )} dx = ··· f(x ) exp{− d (x )} p (x ) dx
0:n 0:n 0:n 0:n k k 0:n 0:n 0:n
E E E E
k=0
n
X
= E[f(X ) exp{− d (X )}] ,
0:n k k
k=0
pourdesfonctions–testf d´efiniessurl’espacedestrajectoiresE = E×···×E.Danslapratique,
n
on verra comment r´esoudre ce probl`eme de mani`ere approch´ee, en simulant des ´echantillons de
variables al´eatoires distribu´ees (approximativement) selon la distribution de Gibbs–Boltzmann
trajectorielle d´efinie (`a une constante multiplicative pr`es) par (1.3).
1.2 Estimation bay´esienne
Dans de nombreux cas, la prise en compte de l’information a priori peut se ramener au
probl`eme statique suivant : ´etant donn´es deux vecteurs al´eatoires X et Y, qu’apporte le fait
d’observer la r´ealisation Y = y sur la connaissance que l’on a de X?
Soit X et Y deux variables al´eatoires `a valeurs dans E et dans F respectivement, et soit
φ une application mesurable d´efinie sur E `a valeurs dans Rp. Par d´efinition, un estimateur de
φ(X) `a partir de l’observation de Y est un vecteur al´eatoire ψ(Y), ou` ψ est une application
mesurable d´efinie sur F `a valeurs dans Rp (par abus de notation, la variable al´eatoire ψ(Y) sera
´egalement not´ee ψ).
I Estimateur MMSE Soit ψ un estimateur de φ(X) sachant Y. Naturellement ψ = ψ(Y)
n’est pas ´egal `a φ(X) : une mesure de l’´ecart entre l’estimateur et la vraie valeur est fournie par
la matrice (de dimension p×p) de corr´elation d’erreur
E[(ψ(Y)−φ(X))(ψ(Y)−φ(X))∗] , (1.4)
dont la trace
trace E[(ψ(Y)−φ(X))(ψ(Y)−φ(X))∗] = E|ψ(Y)−φ(X)|2 ,
est l’erreur quadratique moyenne. L’estimateur du minimum d’erreur quadratique moyenne
(MMSE, pour minimum mean–square error) de φ(X) sachant Y est un estimateur φbtel que
E[(φb(Y)−φ(X))(φb(Y)−φ(X))∗] ≤ E[(ψ(Y)−φ(X))(ψ(Y)−φ(X))∗] ,
au sens des matrices sym´etriques, pour tout autre estimateur ψ.
LaProposition1.1ci–dessousmontrequecetestimateurestobtenu`al’aidedeladistribution
deprobabilit´econditionnelledeX sachantY = y,d´efinie`apartirdeladistributiondeprobabilit´e
jointe de (X,Y) par la d´ecomposition
P[X ∈ dx,Y ∈ dy] = P[X ∈ dx | Y = y] P[Y ∈ dy] . (1.5)