Cours
Cours
cours SOD333
Filtrage Bayésien
et Approximation Particulaire
François Le Gland
INRIA Rennes et IRMAR
http://www.irisa.fr/aspi/legland/ensta/
i
Objectif du cours
1 Introduction 1
1.1 Importance de l’information a priori . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Estimation bayésienne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Cadre gaussien . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2 Exemples 15
2.1 Recalage altimétrique de navigation inertielle . . . . . . . . . . . . . . . . . . . . 15
2.2 Suivi visuel par histogramme de couleur . . . . . . . . . . . . . . . . . . . . . . . 19
2.3 Poursuite d’une cible furtive (track–before-detect) . . . . . . . . . . . . . . . . . 22
2.4 Navigation en environnement intérieur . . . . . . . . . . . . . . . . . . . . . . . . 24
3 Filtrage de Kalman 29
3.1 Systèmes linéaires gaussiens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.2 Filtre de Kalman . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.3 Lisseur de Kalman . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
iii
iv TABLE DES MATIÈRES
7 Filtrage bayésien 73
7.1 Modèles de Markov cachés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
7.2 Chaı̂nes de Markov partiellement observées . . . . . . . . . . . . . . . . . . . . . 78
B Inégalités 179
Introduction
En toute généralité, le filtrage consiste à estimer l’état d’un système dynamique, c’est–à–dire
évoluant au cours du temps, à partir d’observations partielles, généralement bruitées.
Typiquement, on dispose d’une suite Y0 , Y1 , · · · , Yn d’observations, par exemple obtenues
après traitement préalable du signal recueilli au niveau des capteurs. Chaque observation Yn est
reliée à l’état inconnu Xn par une relation du type
Yn = hn (Xn ) + Vn , (1.1)
où Vn est un bruit, qui modélise l’erreur d’observation. On précisera plus loin dans ce cours la
notion de bruit, en terme de variables alátoires le plus souvent centrées (de moyenne nulle).
Une hypothèse assez commune est de supposer que les variables aléatoires V0 , V1 , · · · , Vn sont
indépendantes entre elles. A cause de cette hypothèse d’indépendance mutuelle des bruits d’ob-
servation, et a fortiori en absence de bruit, seule l’observation Yn participe à l’estimation de
l’état caché Xn , c’est–à–dire qu’on se trouve confronté à une succession de problèmes d’estima-
tion découplés : dans la relation (1.1), l’observation Yn est disponible (par définition) tandis que
ni l’état caché Xn ni le bruit Vn ne sont disponibles, et il faut arriver à retrouver (estimer) l’état
caché Xn au vu de l’observation Yn et malgré la présence du bruit Vn .
Tel qu’il est formulé, le problème de l’estimation de l’état caché Xn à partir des observations
Y0 , Y1 , · · · , Yn est en général mal–posé :
• dans le cas favorable où m = d, et même en absence de bruit, il n’est pas toujours possible
d’inverser la relation (1.1) qui peut très bien posséder plusieurs solutions distinctes,
1
2 CHAPITRE 1. INTRODUCTION
Pour lever l’indétermination, c’est–à–dire pour garantir l’existence d’une solution unique, et
pour résoudre le problème de cohérence temporelle, la solution classique consiste à utiliser des
informations supplémentaires sur la suite cachée, par exemple sous la forme de fonctions de
coût portant sur l’état initial ou sur les transitions entre deux états successifs. Par exemple, on
cherchera à minimiser le critère
∑
n ∑
n
J(x0:n ) = c0 (x0 ) + ck (xk−1 , xk ) + dk (xk ) ,
k=1 k=0
par rapport à la suite x0:n = (x0 , x1 , · · · , xn ), qui combine des fonctions de coût représentant
l’information a priori sur la solution avec des fonctions de coût d’une autre nature, qui peuvent
représenter par exemple un terme d’attache aux données, de la forme
hk (x) = 1
2 |Yk − hk (x)|2 ou bien hk (x) = 1
2 (Yk − hk (x))∗ Ik (Yk − hk (x)) ,
pour tout k = 0, 1, · · · , n, avec l’intreprétation que la suite recherchée doit également vérifier à
chaque instant l’équation d’observation en un sens approché. Plus généralement, ces fonctions de
coût peuvent juste représenter une contrainte (ou une propriété) que la suite recherchée devrait
vérifier (ou posséder). En absence d’information a priori , le critère se réduit simplement à
∑
n ∑
n
J(x0:n ) = 1
2 |Yk − hk (xk )| 2
ou bien J(x0:n ) = 1
2 (Yk − hk (xk ))∗ Ik (Yk − hk (xk )) ,
k=0 k=0
par rapport à l’état xk , pour tout k = 0, 1, · · · , n, avec les conséquences déjà évoquées en terme
d’indétermination et de possible incohérence temporelle. Un exemple classique de fonctions de
coût représentant l’information a priori est
c0 (x) = 1
2 |x − µ|2 ou bien c0 (x) = 1
2 (x − µ)∗ Σ−1
0 (x − µ) ,
ck (x, x′ ) = 1
2 |x′ − fk (x)|2 ou bien ck (x, x′ ) = 1
2 (x′ − fk (x))∗ Q−1 ′
k (x − fk (x)) ,
avec l’interprétation que l’état xk recherché doit être proche de fk (xk−1 ), ou de manière équiva-
lente que la transition (xk−1 , xk ) recherchée doit vérifier l’équation xk = fk (xk−1 ) dans un sens
approché, pour tout k = 1, · · · , n. On remarque que ces fonctions de coût sont (à une constante
additive près) de la forme
pour tout k = 1, · · · , n, où p0 (x) est la densité de probabilité initiale, et où pk (x′ | x) est la
densité de probabilité de transition, dans le modèle non–linéaire suivant avec bruits gaussiens
additifs
Xk = fk (Xk−1 ) + Wk avec Wk ∼ N(0, Qk ) ,
et avec condition initiale X0 ∼ N(µ, Σ). En effet (à une constante mutiplicative de normalisation
près)
P[X0 ∈ dx] ∝ exp{− 21 (x − µ)∗ Σ−10 (x − µ)} dx ∝ p0 (x) dx ,
et
pour tout k = 1, · · · , n. En toute généralité, si les relations (1.2) sont vérifiées pour une densité
de probabilité p0 (x) et pour des densités de probabilité de transition pk (x′ | x), pour tout
k = 1, · · · , n, alors le critère à minimiser peut s’écrire
∑
n ∑
n
J(x0:n ) = − log p0 (x0 ) − log pk (xk | xk−1 ) + dk (xk ) ,
k=1 k=0
∏
n ∑
n
exp{−J(x0:n )} = p0 (x0 ) pk (xk | xk−1 ) exp{− dk (xk )} ,
k=1 k=0
| {z }
p0:n (x0:n )
par rapport à la suite x0:n = (x0 , x1 , · · · , xn ). On remarque que p0:n (x0:n ) représente la densité
de probabilité conjointe des états successifs (X0 , X1 , · · · , Xn ) de la chaı̂ne de Markov caractérisée
par
∏
n ∑
n
exp{−J(x0:n )} dx0:n = p0 (x0 ) pk (xk | xk−1 ) exp{− dk (xk )} dx0:n . (1.3)
k=1 k=0
| {z }
p0:n (x0:n )
intégrales) du type
∫ ∫ ∫ ∫ ∑
n
··· f (x0:n ) exp{−J(x0:n )} dx0:n = ··· f (x0:n ) exp{− dk (xk )} p0:n (x0:n ) dx0:n
E E E E k=0
∑
n
= E[f (X0:n ) exp{− dk (Xk )} ] ,
k=0
pour des fonctions–test f définies sur l’espace des trajectoires En = E ×· · ·×E. Dans la pratique,
on verra comment résoudre ce problème de manière approchée, en simulant des échantillons de
variables aléatoires distribuées (approximativement) selon la distribution de Gibbs–Boltzmann
trajectorielle définie (à une constante multiplicative près) par (1.3).
Dans de nombreux cas, la prise en compte de l’information a priori peut se ramener au problème
statique suivant : étant donnés deux vecteurs aléatoires X et Y , qu’apporte le fait d’observer la
réalisation Y = y sur la connaissance que l’on a de X ?
Soit X et Y deux variables aléatoires à valeurs dans E et dans F respectivement, et soit
ϕ une application mesurable définie sur E à valeurs dans Rp . Par définition, un estimateur de
ϕ(X) à partir de l’observation de Y est un vecteur aléatoire ψ(Y ), où ψ est une application
mesurable définie sur F à valeurs dans Rp (par abus de notation, la variable aléatoire ψ(Y ) sera
également notée ψ).
dont la trace
entraı̂ne
b ) − ϕ(X)) (ϕ(Y
= E[ (ϕ(Y b ) − ϕ(X))∗ ] + E[ (ψ(Y ) − ϕ(Y
b )) (ψ(Y ) − ϕ(Y
b ))∗ ]
b )) (ϕ(Y
+ E[ (ψ(Y ) − ϕ(Y b ) − ϕ(X))∗ ] + E[ (ϕ(Y
b ) − ϕ(X)) (ψ(Y ) − ϕ(Y
b ))∗ ] ,
et on remarque que
b )) (ϕ(Y
E[ (ψ(Y ) − ϕ(Y b ) − ϕ(X))∗ ] =
∫ ∫
= b
(ψ(y) − ϕ(y)) b − ϕ(x))∗ P[X ∈ dx, Y ∈ dy]
(ϕ(y)
E F
∫ ∫
= b
(ψ(y) − ϕ(y)) b − ϕ(x))∗ P[X ∈ dx | Y = y] P[Y ∈ dy]
(ϕ(y)
E F
∫ ∫
{ }
= b
(ψ(y) − ϕ(y)) ∗ b − ϕ(x)) P[X ∈ dx | Y = y]
(ϕ(y) P[Y ∈ dy] = 0 ,
F E
b
par définition de ϕ(y). On a donc
b ) − ϕ(X)) (ϕ(Y
= E[ (ϕ(Y b ) − ϕ(X))∗ ] + E[ (ψ(Y ) − ϕ(Y
b )) (ψ(Y ) − ϕ(Y
b ))∗ ]
b ) − ϕ(X)) (ϕ(Y
≥ E[ (ϕ(Y b ) − ϕ(X))∗ ] ,
b
au sens des matrices symétriques, avec égalité pour ψ = ϕ. 2
Remarque 1.2 Compte tenu que le vecteur aléatoire (ϕ(Y b ) − ϕ(X)) est centré, la matrice
de corrélation d’erreur est aussi la matrice de covariance d’erreur, dans le cas particulier de
b
l’estimateur ϕ.
6 CHAPITRE 1. INTRODUCTION
On suppose que la distribution de probabilité jointe des vecteurs aléatoires X et Y possède une
densité
P[X ∈ dx, Y ∈ dy] = p(x, y) dx λ(dy) ,
sur Rm ×F , suffisamment régulière par rapport à la variable x ∈ Rm , avec les deux factorisations
alternatives
p(x, y) = p(x | y) p(y) = p(y | x) p(x) ,
On suppose que
∫ ∫ ∫ ∫ ∫
∂2 { ∂2 }
p(x, y) λ(dy) dx = p(x, y) λ(dy) dx = p′′ (x) dx = 0 .
Rm F ∂x2 Rm ∂x2 F Rm
∂2
J = −E[ log p(X, Y ) ] ,
∂x2
est inversible, alors la matrice de corrélation de l’erreur d’estimation est minorée (au sens des
matrices symétriques) par la relation suivante
M = E[ϕ′ (X)] ,
C ≥ M J −1 M ∗ . 2
8 CHAPITRE 1. INTRODUCTION
b ) − ϕ(X)) (ϕ(Y
E[ (ψ(Y ) − ϕ(X)) (ψ(Y ) − ϕ(X))∗ ] ≥ E[ (ϕ(Y b ) − ϕ(X))∗ ] ≥ M J −1 M ∗ ,
pour tout estimateur ψ, et la borne la plus à gauche est atteinte pour ψ = ϕ. b La borne donnée
par l’estimateur MMSE est donc plus fine que la borne de Cramér–Rao a posteriori, mais aussi
plus difficile à calculer : le plus souvent en effet on ne dispose pas de l’expression de l’estimateur
MMSE, mais l’expression des matrices J et M est assez facile à obtenir. La borne de Cramér–
Rao a posteriori peut même être assez grossière et atteinte par aucun estimateur, et on déduit
de l’encadrement ci–dessus que si la borne de Cramér–Rao a posteriori est atteinte, alors elle
est nécessairement atteinte pour l’estimateur MMSE ψ = ϕ. b
et
p(x) ∝ exp{− 21 (x − X̄)∗ Q−1
X (x − X̄) } ,
de sorte que
= 1
2 (y − h(x))∗ Q−1 1 ∗ −1
V (y − h(x)) + 2 (x − X̄) QX (x − X̄) + cste ,
et
∂2
− log p(x, y) = (h′ (x))∗ Q−1 ′ ∗ −1 ′′ −1
V h (x) − (y − h(x)) QV h (x) + QX ,
∂x2
d’où l’expression de la matrice d’information de Fisher
∂2
J = −E[ log p(X, Y )] = E[(h′ (X))∗ Q−1 ′ ∗ −1 ′′ −1
V h (X)] − E[V QV h (X)] + QX
∂x2
J = H ∗ Q−1 −1
V H + QX et J −1 = QX − QX H ∗ (H QX H ∗ + QV )−1 H QX ,
Dans le cas particulier des vecteurs aléatoires gaussiens, le résultat général obtenu ci–dessus
peut être précisé de la façon suivante.
et de matrice de covariance
R = QX − QXY Q−1
Y QY X ,
complément de Schur de la matrice QY dans la matrice–bloc QZ .
Remarque 1.8 Pour simuler un vecteur aléatoire gaussien de dimension m, de moyenne X(y) b
′ ′ ′
et de matrice de covariance R, il suffit de simuler un vecteur aléatoire gaussien Z = (X , Y ) de
dimension m + d, de même moyenne et de même matrice de covariance que Z = (X, Y ), et de
poser
ξ(y) = X ′ + QXY Q−1 ′
Y (y − Y ) .
On vérifie en effet que le vecteur aléatoire ξ(y) ainsi défini est gaussien, comme transformation
affine du vecteur aléatoire gaussien Z ′ , de moyenne
et de matrice de covariance
b
E[ (ξ(y) − X(y)) b
(ξ(y) − X(y)) ∗
]
−1
= E[ ((X ′ − X̄) − QXY Q−1 ′ ′ ′ ∗
Y (Y − Ȳ )) ((X − X̄) − QXY QY (Y − Ȳ )) ]
− QXY Q−1 ′ ′ ∗ −1 ′ ′ ∗ −1
Y E[ (Y − Ȳ ) (X − X̄) ] + QXY QY E[ (Y − Ȳ ) (Y − Ȳ ) ] QY QY X
= QX − QXY Q−1 −1 −1 −1
Y QY X − QXY QY QY X + QXY QY QY QY QY X
= QX − QXY Q−1
Y QY X = R ,
10 CHAPITRE 1. INTRODUCTION
b
ξ(y) − X(y) = X ′ + QXY Q−1 ′ −1
Y (y − Y ) − (X̄ + QXY QY (y − Ȳ ))
par différence.
0 ≤ R ≤ QX ,
au sens des matrices symétriques (la majoration est immédiate et la minoration résulte du
Lemme A.3), c’est–à–dire que l’utilisation de l’information supplémentaire Y = y, ne peut que
réduire l’incertitude que l’on a sur le vecteur aléatoire X. En outre, la matrice R ne dépend pas
de y, et peut donc être calculée avant même de disposer de la valeur prise par l’observation Y .
b = X(Y
Remarque 1.10 Soit X b ) l’estimateur du minimum de variance de X sachant Y .
Compte tenu que
b = X̄ + QXY Q−1 (Y − Ȳ ) ,
X Y
b = X̄ + QXY Q−1
X(y) Y (y − Ȳ )
( ) −1
QXY ′ QXY ′′ QY ′ 0 y ′ − Ȳ ′
= X̄ +
0 QY ′′ y ′′ − Ȳ ′′
= X̄ + QXY ′ Q−1 ′ ′ −1 ′′ ′′
Y ′ (y − Ȳ ) + QXY QY ′′ (y − Ȳ ) ,
′′
1.3. CADRE GAUSSIEN 11
et de matrice de covariance
R = QX − QXY Q−1
Y QY X
( ) −1
QXY ′ QXY ′′ QY ′ 0 QY ′ X
= QX −
0 QY ′′ QY ′′ X
= QX − QXY ′ Q−1 ′ ′′
−1
Y ′ QY X − QXY QY ′′ QY X .
′′
b ) = X̄ + QX H ∗ (H QX H ∗ + QV )−1 (Y − H X̄) ,
X(Y
R = QX − QX H ∗ (H QX H ∗ + QV )−1 H QX ,
ξ(Y ) = X ′ + QX H ∗ (H QX H ∗ + QV )−1 (Y − (H X ′ + V ′ )) .
Si en outre la matrice QX est inversible, alors il découle du Lemme A.1 d’inversion matricielle
que la matrice R est inversible, et
R−1 = H ∗ Q−1 −1
V H + QX = J ,
d’après l’expression obtenue dans l’Exemple 1.6 pour la matrice d’information de Fisher. Dans
ce cas particulier, la borne de Cramér–Rao a posteriori est donc atteinte, puisque
b ) − X) (X(Y
E[ (X(Y b ) − X)∗ ] = R = J −1 .
12 CHAPITRE 1. INTRODUCTION
Pour finir, on peut montrer directement la relation J = R−1 sans utiliser l’expression obtenue
dans l’Exemple 1.6. En effet, si la matrice R est inversible, ce qui est garanti dès que les matrices
QX et QV sont inversibles, alors on a
b
p(x | y) ∝ exp{− 12 (x − X(y))∗ −1 b
R (x − X(y))} ,
de sorte que
− log p(x | y) = 1 b
(x − X(y))∗ −1 b
R (x − X(y)) + cste ,
2
et
∂2
− log p(x | y) = R−1 ,
∂x2
et on retrouve bien l’expression de la matrice d’information de Fisher
∂2
J = −E[ log p(X | Y )] = R−1 .
∂x2
la matrice de covariance
− QXY Q−1 ∗ −1 ∗ −1
Y E[(Y − Ȳ ) (X − X̄) ] + QXY QY E[(Y − Ȳ ) (Y − Ȳ ) ] QY QY X
= QX − QXY Q−1
Y QY X = R ,
et la matrice de corrélation
QΞ Y = E[(Ξ − Ξ̄) (Y − Ȳ )∗ ]
par différence. En particulier, les vecteurs aléatoires gaussiens Ξ et Y sont décorrélés, donc
indépendants. Il suffit alors d’exprimer la fonction caractéristique de la distribution de probabi-
lité conditionnelle du vecteur aléatoire X = Ξ + QXY Q−1Y Y sachant Y
b ) − 1 u∗ R u} .
= exp{i u∗ X(Y 2
b ) et de
On reconnait la fonction caractéristique d’un vecteur aléatoire gaussien de moyenne X(Y
matrice de covariance R. 2
Conclusion On voit qu’il est important de disposer d’une information a priori sur l’état
inconnu Xn , par exemple de disposer d’une équation d’état décrivant l’évolution de Xn quand
n varie. On peut considérer deux types de modèles :
et dans chacun de ces deux cas, il est possible de résoudre exactement le problème de filtrage
de façon optimale, par la mise en œuvre :
Ces deux cas peuvent être vus comme des cas particuliers de modèles beaucoup plus généraux :
• les chaı̂nes de Markov à espace d’état quelconque (fini, dénombrable, continu, hybride,
etc.),
et dans ce cas il ne sera pas possible de résoudre exactement le problème de filtrage de façon
optimale, qui s’exprime pourtant très simplement en termes de distributions de Feynman–Kac, et
il faudra avoir recours à la mise en œuvre de méthodes de résolution approchées, en l’occurrence :
Exemples
Un avion survole une zone dont le relief est connu : la hauteur h(r) du relief en chaque point de
coordonnée horizontale r est connue, et enregistrée dans une carte numérique.
Dans la suite, la position horizontale de l’avion est notée r, la position verticale, ou altitude,
est notée z, et la vitesse horizontale est notée v. A l’instant 0, la position horizontale initiale
de l’avion est r0 , son altitude initiale est z0 et sa vitesse horizontale initiale est v0 . En réalité,
l’avion se déplace à l’altitude z = z0 constante et à la vitesse horizontale constante v = v0 .
Pour effectuer la navigation, c’est–à–dire pour permettre à l’avion d’estimer sa propre po-
sition horizontale rk et sa propre vitesse horizontale vk à chaque instant tk , on recueille (au
15
16 CHAPITRE 2. EXEMPLES
On introduit comme nouvelles variables d’état les erreurs d’estimation inertielle en position
horizontale δrk = rkINS − rk et en vitesse horizontale δvk = vkINS − vk , et le modèle d’état
2.1. RECALAGE ALTIMÉTRIQUE DE NAVIGATION INERTIELLE 17
1
0
11
0
0
1
0
1
0
1
1
00
1
1
00
1
0
1
01
0
1
01
0
1
0
11
0
0
1
0
1
0
1
1
00
1
0
1
01
0
1
0
11
0 hauteur au−dessus du terrain
0
1
0
1
0
1
1
00
1
1
00
1
0
1
0
11
0
0
11
0
0
position verticale 1
0
1
0
1
1
00
1
1
00
1
0
1
0
11
0
0
1
0
1
0
1
1
00
1
0
1
0
11
0
0
1
0
1
0
1
01
0
1
1
00
1
0 terrain
1
0
11
0
0
1
01
0
1
01
0 altitude du terrain
1
01
0
01
10
niveau zéro
2
un bruit blanc gaussien centré de variance σALT . Au même instant tk , le baro–altimètre fournit
BAR
une mesure bruitée zk de l’altitude de l’avion, c’est–à–dire
zkBAR = zk + wkBAR ,
où zk dénote l’altitude réelle de l’avion, et où la suite wkBAR est un bruit blanc gaussien centré de
2
variance σBAR . La hauteur du relief survolé à l’instant tk déduite à partir des mesures fournies
par le radio–altimètre et par le baro–altimètre est donc
hALT
k = zkBAR − dALT
k = h(rk ) + wkBAR − wkALT ,
et peut être reliée à l’erreur de position inertielle horizontale δrk par
hALT
k = h(rkINS − δrk ) + wkBAR − wkALT .
6600
6500
6400
6300
6200
6100
6000
5900
5800
5700
5600
0 10 20 30 40 50 60 70 80 90 100
En résumé, le modèle d’état utilisé pour le recalage altimétrique de navigation inertielle com-
prend :
• l’équation d’état
( ) ( ) ( ) ( )
δrk I2 ∆ I2 δrk−1 0
= +∆ ,
δvk 0 I2 δvk−1 wkINS
2
où la suite wkINS est un bruit blanc gaussien centré de variance σINS I2 ,
• la condition initiale
( ) ( )
δr0 σr20 I2 0
gaussienne, centrée, de matrice de covariance ,
δv0 0 σv20 I2
2.2. SUIVI VISUEL PAR HISTOGRAMME DE COULEUR 19
• et l’équation d’observation
hALT
k = h(rkINS − δrk ) + wkBAR − wkALT .
où la suite wkALT est un bruit blanc gaussien centré de variance σALT
2 , et où la suite wkBAR
2
est un bruit blanc gaussien centré de variance σBAR .
L’estimation inertielle horizontale rkINS fournie par la centrale inertielle, et la mesure hALT
k de
la hauteur du relief fournie par le radio–altimètre et par le baro–altimètre sont disponibles. La
fonction r 7→ h(r) n’est pas connue de façon analytique, mais définie point–par–point en allant
lire la carte numérique.
On souhaite réaliser un algorithme de suivi dans une séquence d’images numériques couleur. A
la lecture de la première image de la séquence, l’utilisateur sélectionne une zone de l’image, et
le suivi s’effectue de façon séquentielle sur l’ensemble de la séquence, voir Figure 2.5.
...
initialisation image 2 image 3 image 10
La méthode est construite sur l’algorithme SIR (souvent appelé algorithme condensation,
pour conditional density propagation, en vision par ordinateur). Elle repose sur l’hypothèse que
l’histogramme de couleur de la zone à suivre est constant le long de la séquence. Pour avoir plus
d’informations sur cette méthode de suivi visuel, on pourra lire [20].
On désigne sous le terme d’image numérique toute image (dessin, icône, photographie, etc.)
acquise, créée, traitée ou stockée sous forme binaire. On distingue généralement deux grandes
catégories d’images :
• les images vectorielles, dont la description informatique est composée d’objets géométriques
individuels (segments de droite, polygones, arcs de cercle, etc.), chacun définis par divers
attributs de forme, de position, de couleur, etc.
• les images matricielles, représentées par un tableau à deux dimensions dont chaque case
est un pixel (mot dérivé de l’anglais picture element, élément d’image). A chaque pixel
est associée une ou plusieurs valeurs décrivant son niveau de gris ou sa couleur.
20 CHAPITRE 2. EXEMPLES
Les images vectorielles sont utilisées essentiellement pour du graphisme ou en CAO. Lorsque
l’on s’interesse au traitement d’images et à la vision par ordinateur, la représentation utilisé est
la forme matricielle. Il existe plusieurs standards de codage de la couleur :
bitmap noir et blanc : en stockant un bit dans chaque case, il est possible de définir deux
couleurs (noir ou blanc).
bitmap 256 niveaux de gris : en stockant un octet dans chaque case, il est possible de définir
256 dégradés de gris allant du noir au blanc
palette de couleurs (colormap) : grâce à cette méthode, il est possible de définir une palette,
ou table des couleurs, contenant l’ensemble des couleurs pouvant être contenues dans
l’image, à chacune desquelles est associé un indice. Le nombre de bits réservé au codage
de chaque indice de la palette détermine le nombre de couleurs pouvant être utilisées. On
appelle ainsi image en couleurs indexées une image dont les couleurs sont codées selon
cette technique.
couleurs vraies (true color) : le codage de la couleur est réalisé sur trois octets, chaque
octet représentant la valeur d’une composante couleur par un entier de 0 à 255. Ces trois
valeurs codent généralement la couleur dans l’espace RVB (rouge, vert, bleu), mais d’autres
espaces de couleurs peuvent être utilisé. Le nombre de couleurs différentes pouvant être
ainsi représentées est de 256 x 256 x 256 possibilités, soit près de 16 millions de couleurs.
Une image numérique est avant tout un signal 2D. D’un point de vue mathématique, on considère
l’image comme une fonction de R × R dans Ω où le couplet d’entrée est une position spatiale
sur la grille des pixels, et où Ω est l’espace des valeurs de codage de la couleur (ou du niveau
de gris). Par extension, on parlera d’images en dimension 2D+t (t pour le temps) pour désigner
une séquence d’images numériques (ou vidéo numérique).
Remarque 2.1 Etant donné que l’écran effectue un balayage de gauche à droite et de haut en
bas, on désigne généralement par les coordonnées (0, 0) le pixel situé en haut à gauche de l’image,
ce qui signifie que les axes de l’image sont orientés de la façcon suivante : l’axe X est orienté de
gauche à droite, l’axe Y est orienté de haut en bas, contrairement aux notation conventionnelles
en mathématiques, où l’axe Y est orienté vers le haut.
Le but de cet algorithme est de suivre une région d’intérêt dans une séquence d’images. Cette
région est initialisée par l’utilisateur et sa forme est fixée a priori. On considèrera ici un rectangle,
paramétré par la position, en pixel, du centre du rectangle d = (x, y) et un paramètre d’échelle
s. Au pas de temps k (i.e. à l’image k), l’état du système à estimer sera donc Xk = (dk , sk ). Le
paramètre d’échelle permet de suivre un objet même si celui-ci avance ou s’éloigne dans l’axe
de la caméra (effet de zoom). A l’initialisation, l’utilisateur clique 4 points dans l’image, qui
vont définir le rectangle initial. Celui-ci est décrit par les coordonnées du point haut/gauche,
une largeur et une hauteur.
2.2. SUIVI VISUEL PAR HISTOGRAMME DE COULEUR 21
Équation d’état On s’intéresse à la situation où aucune information a priori n’est disponible
sur la nature de l’objet suivi. Dans ce cas, l’équation dynamique du système doit être peu
informative. On supposera donc un modèle à position constante
Xk = Xk−1 + Wk ,
Modèle de couleur La zone initiale à suivre est caractérisée par un histogramme de couleur.
Cet histogramme de référence est construit sur les N b couleurs les plus représentatives de cette
zone, comme montré sur la Figure 2.6. Cet histogramme de référence est noté q ∗ = {q ∗ (n) , n =
1, · · · , N b}, où q ∗ (n) représente le nombre normalisé de pixels de la zone initiale dont la couleur
∑
Nb
la plus proche est la couleur n. On a q ∗ (n) = 1. Pour plus d’informations sur les différents
n=1
espaces de couleur, on pourra se reporter à la page color space sous Wikipedia.
Comme décrit précédement, le but est de suivre une zone de l’image le long de la séquence,
sous l’hypothèse que son histogramme de couleur est invariant dans le temps. Au temps k,
l’histogramme de couleur qk (x) d’un état hypothèse x sera comparé au modèle de couleur de
référence q ∗ , et on définit la mesure de distance D entre ces deux histogrammes de couleur
normalisés
∑
Nb √
D(q ∗ , qk (x)) = ( 1 − q ∗ (n) qk (x, n) )1/2 ,
n=1
Pour favoriser les états hypothèses dont l’histogramme de couleur associé est proche de l’histo-
gramme de référence, on introduit la fonction de pondération
Une image radar est constituée par un tableau rectangulaire de p × p pixels, où l’intensité de
l’écho recueilli en un point est codée par un niveau de gris allant du plus foncé (écho de faible in-
tensité) au plus clair (écho de forte intensité). La même situation se rencontre avec un dispositive
opto–électronique, comme une caméra matricielle, où chaque pixel reçoit et affiche une intensité
lumineuse différente. En principe, si une cible est présente dans la scène 3D visée, elle apparaı̂tra
dans le plan–image sous la forme d’un pixel plus clair (ou d’un groupe de pixels adjacents plus
clairs) que les autres pixels de l’image, lesquels correspondent à l’écho d’objets secondaires de
moindre intensité et/ou à un bruit spatial, indépendant ou bien spatialement corrélé d’un pixel
à l’autre. Pour détecter (et localiser) la cible, il suffit en principe de rechercher dans l’image le
pixel (ou le groupe de pixels adjacents) le plus clair, c’est–à–dire de plus forte intensité. Au lieu
d’une recherche exhaustive, on utilise souvent une méthode de seuillage : rechercher les pixels
d’intensité supérieure à un seuil bien choisi, permet souvent d’obtenir directement le pixel de
plus forte intensité. En répétant cette opération pour chaque image successivement on peut ainsi
détecter d’abord, puis suivre, la cible dans une séquence d’images.
50 50
0 0
−50 −50
−100 −100
−100 −50 0 50 100 −100 −50 0 50 100
60 50
40 0
20 −50
0 −100
−1 0 1 2 −100 −50 0 50 100
Figure 2.7 – Image observée, position réelle, histogramme, détection (cible visible)
On s’intéresse ici au cas d’une cible furtive, caractérisée par un écho de très faible intensité,
c’est–à–dire d’une intensité du même ordre de grandeur que l’intensité caractéristique du bruit
présent dans l’image, voire même d’un ordre de grandeur inférieur. Dans ce cas, une méthode de
seuillage est inefficace : quel que soit le seuil choisi, rechercher les pixels d’intensité supérieure au
seuil ne permet plus d’isoler la cible au milieu du bruit. Même un opérateur humain est incapable,
sur une image isolée, de détecter la présence et la position de la cible. En revanche, un opérateur
2.3. POURSUITE D’UNE CIBLE FURTIVE (TRACK–BEFORE-DETECT) 23
humain est capable dans certains cas de suivre la cible dans une séquence d’images, comme une
succession de pixels (un dans chaque image de la séquence) animés d’un mouvement cohérent
au milieu de l’agitation incoordonnée des autres pixels. En quelque sorte, l’œil humain suit la
cible sans jamais la détecter vraiment : c’est ce genre de performance qu’il s’agit de reproduire
ici de manière algorithmique, connue sous le terme de track–before–detect, en s’appuyant sur un
modèle a priori pour le déplacement de la cible, qui favorise le mouvement cohérent de pixels
entre des images successives.
50 50
0 0
−50 −50
−100 −100
−100 −50 0 50 100 −100 −50 0 50 100
30 50
20 0
10 −50
0 −100
−2 −1 0 1 2 −100 −50 0 50 100
Figure 2.8 – Image observée, position réelle, histogramme, détection (cible furtive)
Chaque image peut se représenter comme un champ aléatoire (Yk (s) , s ∈ S) où l’indice s ∈ S
désigne le pixel ou de manière équivalente le site d’un réseau bi–dimensionnel. Par hypothèse,
l’intensité observée au pixel s ∈ S se décompose comme
L’intensité due au bruit seulement est modélisée comme un champ aléatoire gaussien (Bk (s) , s ∈
S) centré, de variance σB2 en tout pixel s ∈ S et décorrélé spatialement, c’est–à–dire
On définit le rapport signal à bruit (ou signal to noise ratio, SNR) en decibel, comme
I0
SNR = 20 log10 .
σB
On pose ici (par convention) I0 = 1 de sorte qu’un rapport signal à bruit de 20 dB correspond
à σB = 0.1 tandis qu’un rapport signal à bruit de 0 dB correspond à σB = 1.
La fonction de vraisemblance est donnée à une constante multiplicative près par l’expression,
en fonction de la variable r, de la densité du champ aléatoire observé (Yk (s) , s ∈ S) quand la
cible occupe la position r dans l’espace physique. On a donc par définition
1 ∑
gk (r) = exp{− 2 |Yk (s) − I(r, s)|2 } ,
2 σB
s∈S
et on remarque que
1 ∑ 1 ∑
gk (r) = exp{ 2 I(r, s) Yk (s) − 2 |I(r, s)|2 } ,
σB 2 σB
s∈C(r) s∈C(r)
à une constante multiplicative près, de sorte que le calcul porte seulement sur les 9 pixels du
voisinage C(r), et pas sur l’ensemble S de tous les pixels.
Le modèle a priori pour l’évolution de la cible est donné par le modèle d’état suivant
( ) ( ) ( ) ( )
rk I2 ∆ I2 rk−1 √ 0
= + σ ∆ ,
vk 0 I2 vk−1 wk
où wk est un vecteur aléatoire gaussien centré de matrice de covariance I2 , où la position initiale
r0 est distribuée uniformément dans l’espace physique défini ci–dessus, et où la vitesse initiale
v0 est distribuée uniformément dans le domaine délimité en module par vmin ≤ |v0 | ≤ vmax et
en orientation par [0, 2 π).
Sur chaque image, on peut rechercher le pixel de plus forte intensité observée, ou bien mettre
en œuvre une méthode de seuillage pour détecter les pixels d’intensité supérieure au seuil choisi.
On peut aussi extraire l’histogramme des intensités observées aux différents pixels de l’image.
Si le rapport signal à bruit est trop faible, alors une simple détection image par image s’avère
inefficace. On peut en revanche considérer les images successives comme des observations (ma-
tricielles), et mettre en œuvre un algorithme de filtrage pour effectuer directement le suivi.
Un utilisateur se déplace à l’intérieur d’un bâtiment dont le plan est disponible sous la forme
d’une carte numérique. L’utilisateur est caractérisé à l’instant tk
2.4. NAVIGATION EN ENVIRONNEMENT INTÉRIEUR 25
• par son orientation θk (un vecteur unitaire, ou un angle) par rapport à la direction de
référence correspondant au vecteur unitaire u = (1, 0), dirigé vers la droite sur la carte.
25
50
75
100
0 25 50 75 100 125 150
Le segment numéro k, joignant les positions rk et rk+1 occupées par l’utilisateur aux instants
tk et tk+1 respectivement, peut être caractérisé de la manière équivalente
• par sa longueur dk = |rk+1 − rk |, qui peut s’interpréter comme la distance parcourue par
l’utilisateur entre les instants tk et tk+1 ,
• et par son orientation θk (déja mentionnée), qui peut être définie de manière équivalente
par le vecteur unitaire uk = (rk+1 − rk )/dk ,
• une mesure α bk de la rotation effectuée par l’utilisateur à l’instant tk , avec une incertitude
caractérisée par un bruit gaussien additif de moyenne nulle et de variance σturn 2 ,
• et une mesure dbk de la distance parcourue par l’utilisateur entre les instants tk et tk+1 , avec
une incertitude caractérisée par un bruit gaussien additif de moyenne nulle et de variance
2
σwalk .
En d’autres termes
bk = αk + wkturn
α (modulo 2π) et dbk = dk + wkwalk , (2.1)
où wkwalk et wkturn sont deux variables aléatoires gaussiennes indépendantes, de moyenne nulle et
2
de variance σwalk 2 , respectivement.
et σturn
Les mesures bruitées db1 , · · · , dbnmax−1 et α
b1 , α
b2 , · · · , α
bnmax−1 (avec la convention α
b1 = 0) sont
recueillies par l’utilisateur le long de la trajectoire. À partir de ces mesures PNS incrémentales
bruitées, et à partir d’estimations de la position initiale r1 et de l’orientation initiale θ1 inconnues,
on peut essayer de reconstruire la position et l’orientation de l’utilisateur à chaque instant, par
intégration
θkPNS = θk−1
PNS
bk
+α (modulo 2π) et PNS
rk+1 = rkPNS + dbk u(θkPNS ) ,
où u(θ) = (cos θ, sin θ) désigne le vecteur unitaire associé à l’angle θ. La trajectoire estimée à
partir des mesures PNS seulement est représentée sur la Figure 2.10.
25
50
75
100
0 25 50 75 100 125 150
On remarque que la trajectoire estimée s’écarte de la trajectoire réelle, juste parce que les
erreurs sur les mesures PNS incrémentales s’accumulent au cours du temps.
2.4. NAVIGATION EN ENVIRONNEMENT INTÉRIEUR 27
Pour corriger la dérive de la trajectoire estimée à partir des mesures PNS seulement, l’idée
consiste à recueillir séparément des mesures fournies par d’autres capteurs. Dans la solution
proposée ici, à l’intérieur du bâtiment sont disposées des balises de ranging identiques, dont les
positions sont connues. Chaque balise est caractérisée par sa portée R, de sorte que
• tout utilisateur se trouvant à une distance inférieure à R par rapport à une balise est
détecté par cette balise,
25
50
75
100
0 25 50 75 100 125 150
En outre, si une balise détecte un utilisateur alors une mesure de la distance entre l’utilisateur
et cette balise est également disponible, avec une incertitude caractérisée par un bruit gaussien
2
additif de moyenne nulle et de variance σrange . Les éventuelles détections et mesures de distance
bruitées sont recueillies par l’utilisateur le long de la trajectoire, et sont disponibles pour le
reclalage de navigation.
28 CHAPITRE 2. EXEMPLES
• une fonction de vraisemblance associée à chaque balise active, c’est–à–dire à chaque balise
déclenchée par l’utilisateur,
Filtrage de Kalman
29
30 CHAPITRE 3. FILTRAGE DE KALMAN
et on suppose que
• même si l’état Xk−1 = x est connu exactement à l’instant (k − 1), on peut seulement dire
que l’état Xk à l’instant k est incertain, et distribué comme un vecteur aléatoire gaussien,
de moyenne Fk x + fk et de matrice de covariance QW k ,
• si l’état Xk−1 est incertain à l’instant (k − 1), et distribué comme un vecteur aléatoire
gaussien, de moyenne X̄k−1 et de matrice de covariance QX k−1 , alors cette incertitude se
propage à l’instant k : même en absence de bruit, c’est–à–dire même si Gk = 0, l’état Xk
à l’instant k est incertain, et distribué comme un vecteur aléatoire gaussien, de moyenne
Fk X̄k−1 + fk et de matrice de covariance Fk QX ∗
k−1 Fk .
Proposition 3.1 La suite {Zk = (Xk , Yk )} est une suite gaussienne à valeurs dans Rm+d .
Preuve. Comme sortie d’un système linéaire à entrées gaussiennes, la suite {Zk } est un proces-
sus aléatoire gaussien. En effet, pour tout instant n, le vecteur aléatoire (Z0 , Z1 , · · · , Zn ) peut
s’exprimer comme transformation affine du vecteur aléatoire (X0 , W1 , · · · , Wn , V0 , V1 , · · · , Vn )
qui par hypothèse est un vecteur aléatoire gaussien, donc le vecteur aléatoire (Z0 , Z1 , · · · , Zn )
est gaussien, comme transformation affine d’un vecteur aléatoire gaussien. 2
Remarque 3.2 Si les coefficients dépendent des observations passées, on parle de système
conditionnellement linéaire gaussien : on considère ainsi une suite d’états cachés {Xk } à va-
leurs dans Rm , vérifiant
où la suite {Wk } prend ses valeurs dans Rp , et une suite d’observations {Yk } à valeurs dans Rd ,
vérifiant
Yk = Hk (Y0:k−1 ) Xk + hk (Y0:k−1 ) + Vk ,
et on suppose que
Dans ce cas, la suite {Zk = (Xk , Yk )} n’est en général pas une suite gaussienne, mais on peut
vérifier que conditionnellement à Y0:k−1
Xk = Fk Xk−1 + fk + Wk , (3.3)
Yk = Hk Xk + hk + Vk , (3.4)
Y0:k = (Y0 , Y1 , · · · , Yk ) .
L’objectif est d’estimer de façon optimale et récursive le vecteur aléatoire Xk à partir de Y0:k .
Si on adopte le critère du minimum de variance, il s’agit d’après la Section 1.2 de calculer la
distribution de probabilité conditionnelle du vecteur aléatoire Xk sachant Y0:k . Comme le cadre
est gaussien, il suffit de calculer la moyenne et la matrice de covariance
bk = E[Xk | Y0:k ]
X et bk ) (Xk − X
Pk = E[(Xk − X bk )∗ | Y0:k ] .
b − = E[Xk | Y0:k−1 ]
X et b − ) (Xk − X
Pk− = E[(Xk − X b − )∗ | Y0:k−1 ] .
k k k
D’après la Remarque 1.9, les matrices de covariances conditionnelles Pk et Pk− ne dépendent pas
des observations, c’est–à–dire que
bk ) (Xk − X
Pk = E[(Xk − X bk )∗ ] et b − ) (Xk − X
Pk− = E[(Xk − X b − )∗ ] .
k k
Ik = Yk − E[Yk | Y0:k−1 ] ,
et d’après (3.4), on a
b − + hk ) ,
Ik = Yk − (Hk E[Xk | Y0:k−1 ] + hk + E[Vk | Y0:k−1 ]) = Yk − (Hk X k
Remarque 3.3 Par définition, toute fonction des variables (Y0 , · · · , Yk−1 , Yk ) peut s’exprimer
en fonction des variables (Y0 , · · · , Yk−1 , Ik ), et réciproquement. On en déduit que (Y0:k−1 , Ik )
contient exactement la même information que Y0:k .
Lemme 3.4 Le processus {Ik } est un processus gaussien à valeurs dans Rd , appelé processus
d’innovation. En particulier, le vecteur aléatoire Ik est gaussien, de moyenne nulle et de matrice
de covariance
QIk = Hk Pk− Hk∗ + QVk ,
b − , Ik ) est gaussien, de
et indépendant de Y0:k−1 . Plus généralement, le vecteur aléatoire (Xk − X k
moyenne nulle et de matrice de covariance
Pk− Pk− Hk∗
,
− − ∗ V
Hk Pk Hk Pk Hk + Qk
et indépendant de Y0:k−1 .
Preuve. D’après la Remarque 1.10, l’observation prédite E[Yk | Y0:k−1 ] dépend de façon af-
fine des observations passées (Y0 , Y1 , · · · , Yk−1 ), de sorte que l’innovation Ik dépend de façon
affine des observations (Y0 , Y1 , · · · , Yk ). On en déduit que le vecteur aléatoire (I0 , I1 , · · · , Ik ) est
gaussien, comme transformation affine d’un vecteur aléatoire gaussien.
Toujours d’après la Remarque 1.10, l’état prédit X b − = E[Xk | Y0:k−1 ] dépend de façon affine
k
des observations passées (Y0 , · · · , Yk−1 ), de sorte que le vecteur aléatoire (Y0 , · · · , Yk−1 , Xk −
Xb − , Ik ) dépend de façon affine du vecteur (Y0 , Y1 , · · · , Yk , Xk ) formé de l’état courant Xk et
k
des observations (Y0 , Y1 , · · · , Yk ). On en déduit que le vecteur aléatoire (Y0 , · · · , Yk−1 , Xk −
Xb − , Ik ) est gaussien, et donc a fortiori le vecteur aléatoire (Xk − X b − , Ik ) est gaussien, comme
k k
transformation affine d’un vecteur aléatoire gaussien. Compte tenu que
b − | Y0:k−1 ] = 0
E[Xk − X et E[Ik | Y0:k−1 ] = 0 ,
k
et on en déduit que
QIk = E[Ik Ik∗ ]
b − ) + Vk ) (Hk (Xk − X
= E[(Hk (Xk − X b − ) + Vk )∗ ]
k k
b − ) (Xk − X
= Hk E[(Xk − X b − )∗ ] H ∗ + E[Vk V ∗ ]
k k k k
b − )∗ ] H ∗ + Hk E[(Xk − X
+ E[Vk (Xk − X b −) V ∗]
k k k k
b − ) (Xk − X
= E[(Xk − X b − )∗ ] H ∗ + E[(Xk − X
b −) V ∗]
k k k k k
= Pk− Hk∗ .
b − ) est indépendant de
Dans cette dernière égalité, on a de nouveau utilisé le fait que (Xk − Xk
Vk , donc E[(Xk − Xb − ) V ∗ ] = 0. 2
k k
Remarque 3.5 Si la matrice de covariance QVk est inversible, alors a fortiori la matrice de
covariance QIk = Hk Pk− Hk∗ + QVk est inversible, pour tout instant k.
Remarque 3.6 Compte tenu que la distribution de probabilité conditionnelle du vecteur aléa-
toire Yk sachant Y0:k−1 est gaussienne, de moyenne Hk X b − + hk et de matrice de covariance QI ,
k k
et pourvu que la matrice QIk soit inversible, on obtient l’expression suivante
∏
n
Ln = b − + hk ))∗ (QI )−1 (Yk − (Hk X
exp{− 12 (Yk − (Hk X b − + hk )) }
k k k
k=0
∏
n
= exp{− 21 Ik∗ (QIk )−1 Ik } ,
k=0
Théorème 3.7 (Filtre de Kalman) On suppose que la matrice de covariance QVk est inver-
bk } et {Pk } vérifient les équations récurrentes
sible, pour tout instant k. Alors les suites {X
suivantes
b − = Fk X
X bk−1 + fk ,
k
et
X b − + Kk (Yk − (Hk X
bk = X b − + hk )) ,
k k
Pk = (I − Kk Hk ) Pk− ,
où la matrice
Kk = Pk− Hk∗ [Hk Pk− Hk∗ + QVk ]−1 ,
est appelée gain de Kalman, et avec les initialisations
b − = X̄0 = E[X0 ]
X et P0− = QX
0 = cov(X0 ) .
0
Remarque 3.9 On vérifie que la suite {Pk } ne dépend pas des observations : elle peut donc être
pré–calculée, en particulier dans le cas simple où les coefficients Fk = F , Hk = H, QWk = QW
et QVk = QV sont constants.
Preuve. On procède en plusieurs étapes. Le point central est la Proposition 1.7 qui sera
constamment utilisée.
b0 et P0 en fonction de X
Expression de X b − et P − :
0 0
Le vecteur aléatoire (X0 , Y0 ) est gaussien, de moyenne et de matrice de covariance données
par
b−
X P0− P0− H0∗
0
et ,
b−
H0 X 0 + h0 H0 P0− H0 P0− H0∗ + QV0
respectivement. D’après la Proposition 1.7, la distribution de probabilité conditionnelle du vec-
teur aléatoire X0 sachant Y0 est gaussienne, de moyenne
b0 = X
X b − + P − H0∗ [H0 P − H0∗ + QV0 ]−1 (Y0 − (H0 X
b − + h0 )) ,
0 0 0 0
et de matrice de covariance
P0 = P0− − P0− H0∗ [H0 P0− H0∗ + QV0 ]−1 H0 P0− .
3.2. FILTRE DE KALMAN 35
b − et P − en fonction de X
Expression de X bk−1 et Pk−1 :
k k
Le vecteur aléatoire (Xk , Y0 , · · · , Yk−1 ) est gaussien, et d’après la Proposition 1.7, la distri-
bution de probabilité conditionnelle du vecteur aléatoire Xk sachant Y0:k−1 est gaussienne, de
moyenne Xb − et de matrice de covariance P − . D’après l’équation (3.3), c’est–à–dire
k k
Xk = Fk Xk−1 + fk + Wk ,
on a
b − = Fk (Xk−1 − X
Xk − X bk−1 ) + Wk ,
k
de sorte que
b − ) (Xk − X
Pk− = E[(Xk − X b − )∗ ]
k k
bk−1 ) (Xk−1 − X
= Fk E[(Xk−1 − X bk−1 )∗ ] F ∗ + E[Wk W ∗ ]
k k
bk−1 )∗ ] F ∗ + Fk E[(Xk−1 − X
+ E[Wk (Xk−1 − X bk−1 ) W ∗ ]
k k
= Fk Pk−1 Fk∗ + QW
k .
bk et Pk en fonction de X
Expression de X b − et P − :
k k
bk = E[Xk | Y0:k ]
X
b − + E[Xk − X
= X b − | Y0:k ]
k k
b − + E[Xk − X
= X b − | Y0:k−1 , Ik ]
k k
b − + E[Xk − X
= X b − | Ik ] ,
k k
36 CHAPITRE 3. FILTRAGE DE KALMAN
de sorte que
bk ) (Xk − X
Pk = E[ (Xk − X bk )∗ ]
b − ) − E[Xk − X
= E[ ((Xk − X b − | Ik ]) ((Xk − X
b − ) − E[Xk − X
b − | Ik ])∗ ] .
k k k k
et
Pk = Pk− − Pk− Hk∗ [Hk Pk− Hk∗ + QVk ]−1 Hk Pk− ,
ce qui termine la démonstration. 2
bn = X
X bk−1 + Lk (X
bn − X
b −) ,
k−1 k k
n
Pk−1 = Pk−1 + Lk (Pkn − Pk− ) L∗k ,
on vérifie que la matrice Pk−1 − Lk Pk− L∗k est semi–définie positive, pour tout instant k. On en
déduit par récurrence arrière que la matrice Pkn (telle qu’elle est définie par l’équation rétrograde
de l’énoncé) est semi–définie positive, pour tout instant k. Par définition, Pnn = Pn , c’est–à–dire
que la relation est vraie au rank k = n. Si la relation est vraie au rang k, c’est–à–dire si la
matrice Pkn est semi–définie positive, alors nécessairement la matrice
n
Pk−1 = Pk−1 + Lk (Pkn − Pk− ) L∗k = (Pk−1 − Lk Pk− L∗k ) + Lk Pkn L∗k ,
aussi est semi–définie positive, c’est–à–dire que la relation est vraie au rang (k − 1).
Remarque 3.13 On vérifie par récurrence arrière que Pkn ≤ Pk , c’est–à–dire que la matrice de
covariance de l’erreur de lissage est plus petite (au sens des matrices symétriques) que la matrice
de covariance de l’erreur de filtrage, pour tout instant k. Par définition Pnn = Pn , c’est–à–dire
que la relation est vraie au rank k = n. Si la relation est vraie au rang k, c’est–à–dire si Pkn ≤ Pk ,
alors nécessairement Pkn ≤ Pk− compte tenu que Pk ≤ Pk− d’après la Remarque 3.8. En d’autres
termes, la différence (Pkn − Pk− ) est semi–définie négative, de sorte que la différence
n
Pk−1 − Pk−1 = Lk (Pkn − Pk− ) L∗k ,
n
aussi est semi–définie négative. En d’autres termes, Pk−1 ≤ Pk−1 , c’est–à–dire que la relation
est vraie au rang (k − 1).
Preuve. On remarque que le vecteur aléatoire Yk peut s’exprimer comme transformation affine
du vecteur aléatoire (Xk , Vk ), et donc a fortiori comme transformation affine du vecteur aléatoire
b − , Vk ). De même, le vecteur aléatoire Yk+p peut s’exprimer comme transformation
(Y0:k−1 , Xk − Xk
affine du vecteur aléatoire (Xk+p , Vk+p ), et par transitivité comme transformation affine du
vecteur aléatoire (Xk , Wk+1 , · · · , Wk+p , Vk+p ), et donc a fortiori comme transformation affine
du vecteur aléatoire (Y0:k−1 , Xk − X b − , Wk+1 , · · · , Wk+p , Vk+p ). On en déduit que le vecteur
k
aléatoire Y0:n = (Y0:k−1 , Yk , · · · , Yn ) peut s’exprimer comme transformation affine du vecteur
38 CHAPITRE 3. FILTRAGE DE KALMAN
n
Uk−1 b − , Zk+1:n ]
= E[Xk−1 | Y0:k−1 , Xk − X k
bk−1 + E[Xk−1 − X
= X bk−1 | Y0:k−1 , Xk − X
b − , Zk+1:n ]
k
bk−1 + E[Xk−1 − X
= X bk−1 | Y0:k−1 ] + E[Xk−1 − X
bk−1 | Xk − X
b −]
k
bk−1 | Zk+1:n ]
+ E[Xk−1 − X
bk−1 + E[Xk−1 − X
= X bk−1 | Xk − X
b −] ,
k
compte tenu que E[Xk−1 − X bk−1 | Y0:k−1 ] = 0 par définition, et où on a utilisé dans la dernière
b
égalité le fait que (Xk−1 − Xk−1 ) est indépendant de Zk+1:n , donc E[Xk−1 − X bk−1 | Zk+1:n ] = 0.
Par différence
Xk−1 − Uk−1
n bk−1 ) − (U n − X
= (Xk−1 − X bk−1 )
k−1
bk−1 ) − E[Xk−1 − X
= (Xk−1 − X bk−1 | Xk − X
b −] ,
k
de sorte que
E[(Xk−1 − Uk−1
n
) (Xk−1 − Uk−1
n
)∗ ]
bk−1 ) − E[Xk−1 − X
= E[ ((Xk−1 − X bk−1 | Xk − X
b − ])
k
bk−1 ) − E[Xk−1 − X
((Xk−1 − X bk−1 | Xk − X
b − ])∗ ] .
k
b − = Fk (Xk−1 − X
Xk − X bk−1 ) + Wk ,
k
3.3. LISSEUR DE KALMAN 39
de sorte que
bk−1 ) (Xk − X
E[(Xk−1 − X b − )∗ ]
k
bk−1 ) (Xk−1 − X
= E[(Xk−1 − X bk−1 )∗ ] F ∗ + E[(Xk−1 − X
bk−1 ) W ∗ ]
k k
= Pk−1 Fk∗ .
Dans cette dernière égalité, on a utilisé le fait que (Xk−1 − Xbk−1 ) et Wk sont indépendants, donc
b ∗
E[(Xk−1 − Xk−1 ) Wk ] = 0. On en déduit que le vecteur aléatoire gaussien (Xk−1 − X bk−1 , Xk − X
b −)
k
est de moyenne nulle et de matrice de covariance
Pk−1 Pk−1 Fk∗
.
−
Fk Pk−1 Pk
Par hypothèse la matrice Pk− est inversible, et d’après la Proposition 1.7 on a immédiatement
n
Uk−1 bk−1 + Pk−1 F ∗ (P − )−1 (Xk − X
=X b −) = X
bk−1 + Lk (Xk − X
b −) ,
k k k k
et
E[(Xk−1 − Uk−1
n
) (Xk−1 − Uk−1
n
)∗ ] = Pk−1 − Pk−1 Fk∗ (Pk− )−1 Fk Pk−1 = Pk−1 − Lk Pk− L∗k .
de sorte que
n
Pk−1 b n ) (Xk−1 − X
= E[ (Xk−1 − X b n )∗ ]
k−1 k−1
= E[ ((Xk−1 − Uk−1
n n
) + (Uk−1 b n )) ((Xk−1 − U n ) + (U n − X
−X b n ))∗ ]
k−1 k−1 k−1 k−1
= E[ (Xk−1 − Uk−1
n
) (Xk−1 − Uk−1
n
)∗ ] + E[ (Uk−1
n b n ) (U n − X
−X b n )∗ ]
k−1 k−1 k−1
+ E[ (Uk−1
n b n ) (Xk−1 − U n )∗ ] + E[ (Xk−1 − U n ) (U n − X
−X b n )∗ ]
k−1 k−1 k−1 k−1 k−1
• (Uk−1
n b n ) dépend de (Y0:k−1 , Xk − X
−X b − , Zk+1:n ),
k−1 k
• et E[Xk−1 − Uk−1
n b − , Zk+1:n ] = 0 par définition,
| Y0:k−1 , Xk − X k
b n ) et
Proposition 3.14 La matrice de corrélation Ckn entre les erreurs de lissage (Xk−1 − Xk−1
b n ) à deux instants successifs vérifie la relation suivante
(Xk − X k
b n ) (Xk − X
Ckn = E[ (Xk−1 − X b n )∗ ] = Lk P n .
k−1 k k
de sorte que
b n ) (Xk − X
Ckn = E[ (Xk−1 − X b n )∗ ]
k−1 k
= E[ (Xk−1 − Uk−1
n b n )∗ ] + Lk E[ (Xk − X
) (Xk − X b n ) (Xk − X
b n )∗ ]
k k k
= Lk Pkn .
b n ) = (Xk − X
• (Xk − X b − ) + (X
b− − X
b n ) dépend de (Y0:k−1 , Xk − X
b − , Zk+1:n ),
k k k k−1 k
• et E[Xk−1 − Uk−1
n b − , Zk+1:n ] = 0 par définition,
| Y0:k−1 , Xk − X k
bk + Pk F ∗ (P − )−1 (X
= X bn − Xb− )
k+1 k+1 k+1 k+1
bk + Pk rk ,
= X
3.3. LISSEUR DE KALMAN 41
et de même
−
Pkn = Pk + Lk+1 (Pk+1
n
− Pk+1 ) L∗k+1
∗ − − −
= Pk + Pk Fk+1 (Pk+1 )−1 (Pk+1
n
− Pk+1 ) (Pk+1 )−1 Fk+1 Pk
= Pk − Pk Πk Pk ,
pour tout k = 0, 1, · · · , n.
où les suites {rk } et {Πk } vérifient les équations récurrentes rétrogrades suivantes
b − + hk )) ,
rk− = (I − Kk Hk )∗ rk + Hk∗ Ξk (Yk − (Hk X k
Π− ∗ ∗
k = (I − Kk Hk ) Πk (I − Kk Hk ) + Hk Ξk Hk ,
et
rk−1 = Fk∗ rk− et Πnk−1 = Fk∗ Π−
k Fk ,
Pk = (I − Kk Hk ) Pk− = Pk− (I − Kk Hk )∗ ,
de sorte que
et par définition
Kk = Pk− Hk∗ Ξk ,
42 CHAPITRE 3. FILTRAGE DE KALMAN
de sorte que
(Pk− )−1 Kk = Hk∗ Ξk . (3.8)
D’après l’étape de correction du filtre de Kalman, on a
bk = X
X b − + Kk (Yk − (Hk X
b − + hk )) ,
k k
de sorte que
bn − X
X b− = X
bk − X
b − + Pk rk = Kk (Yk − (Hk X
b − + hk )) + Pk rk ,
k k k k
b − + hk )) ]
= Fk∗ (Pk− )−1 [ Pk rk + Kk (Yk − (Hk Xk
b − + hk )) ] ,
= Fk∗ [ (I − Kk Hk )∗ rk + Hk∗ Ξk (Yk − (Hk Xk
compte tenu des identités (3.7) et (3.8). D’après l’étape de correction du filtre de Kalman, on a
de sorte que
= Fk∗ [ (I − Kk Hk )∗ Πk (I − Kk Hk ) + Hk∗ Ξk Hk ] Fk ,
résolution de systèmes linéaires de dimension m de la forme Pk− x = b est requise, et passe par
exemple par la décomposition de Cholesky de la matrice Pk− . L’équation récurrente rétrograde
pour le calcul du lisseur X b n utilise les valeurs numériques du filtre X bk−1 et de la matrice de
k−1
covariance d’erreur de filtrage Pk−1 (à partir desquelles il est facile de reconstruire les valeurs
numériques du prédicteur X b − et de la matrice de covariance d’erreur de prédiction P − ). Ces
k k
valeurs numériques sont calculées dans la phase aller, et doivent donc être conservées en mémoire
pour être utilisées dans la phase retour. En revanche, cette équation récurrente rétrograde pour le
calcul du lisseur X b n n’utilise ni la valeur numérique de l’observation Yk ni celle de l’innovation
k−1
Ik = Yk − (Hk X b − + hk ).
k
Dans la formulation de Fraser–Potter, qui fait l’objet du Théorème 3.15, il n’y a pas de condition
nécessaire d’inversibilité pour la phase retour qui ne soit pas déjà nécessaire pour la phase aller.
Les expressions (3.6) pour le lisseur X b n et pour la matrice de covariance d’erreur de lissage P n
k k
utilisent les valeurs numériques du filtre X bk et de la matrice de covariance d’erreur de filtrage
Pk . Ces valeurs numériques sont calculées dans la phase aller, et doivent donc être conservées en
mémoire pour être utilisées dans la phase retour. L’équation récurrente rétrograde pour le calcul
de la variable rk utilise la valeur numérique de l’observation Yk ou de manière équivalente celle
de l’innovation Ik = Yk − (Hk X b − + hk ). Ces valeurs numériques sont calculées dans la phase
k
aller, et doivent donc être conservées en mémoire pour être utilisées dans la phase retour.
En conclusion :
• les deux formulations requièrent dans la phase aller une même condition d’inversibilité et
l’inversion de systèmes linéaires de dimension d,
• les deux formulations utilisent dans la phase retour les valeurs numériques du filtre et de
la matrice de covariance d’erreur de filtrage — ces valeurs numériques sont calculées dans
la phase aller, et doivent donc être conservées en mémoire pour être utilisées dans la phase
retour,
Remarque 3.16 Il est également possible d’obtenir une équation récurrente pour le lisseur,
dans le sens direct (et pas dans le sens rétrograde) et autonome (ne faisant pas intervenir ni le
filtre ni la matrice de covariance de l’erreur de filtrage). Par différence, on obtient
b n − Fk X
X bn = Xbk + Pk rk − Fk (X
bk−1 + Pk−1 rk−1 )
k k−1
bk − Fk X
= (X bk−1 ) + (Pk rk − Fk Pk−1 rk−1 ) .
44 CHAPITRE 3. FILTRAGE DE KALMAN
et on remarque que
b − + hk )) ]
Pk− rk− = Pk− [ (I − Kk Hk )∗ rk + Hk∗ Ξk (Yk − (Hk X k
b − + hk )) ,
= Pk rk + Kk (Yk − (Hk X k
de sorte que
bk − Fk X
(X bk−1 ) + (Pk rk − P − r− ) = 0 .
k k
D’autre part
= Pk rk − (Pk− − QW −
k ) rk
On en déduit que
X b n + (X
b n = Fk X bk−1 ) + (Pk rk − Fk Pk−1 rk−1 )
bk − Fk X
k k−1
b n + (X
= Fk X bk − Fk X
bk−1 ) + (Pk rk − P − r− ) + QW r−
k−1 k k k k
b n + QW r− ,
= Fk X k−1 k k
c’est–à–dire qu’on obtient une équation récurrente, dans le sens direct, et faisant seulement
intervenir la variable rk− .
Chapitre 4
où {Wk } prend ses valeurs dans Rp , et une suite d’observations {Yk } à valeurs dans Rd , vérifiant
Yk = hk (Xk ) + Vk , (4.2)
et on suppose que
• la suite {Vk } est un bruit blanc gaussien, de matrice de covariance QVk inversible,
• même si l’état Xk−1 = x est connu exactement à l’instant (k − 1), on peut seulement dire
que l’état Xk à l’instant k est incertain, et distribué comme un vecteur aléatoire gaussien,
de moyenne bk (x) et de matrice de covariance σk (x) σk∗ (x).
La plupart des propriétés obtenues à la Section 3.1 ne sont pas vraies pour le système décrit
par les équations (4.1) et (4.2). En particulier, le processus {Zk = (Xk , Yk )} n’est pas gaussien (ni
même conditionnellement gaussien), et les moments conditionnels de Xk sachant Y0:k ne peuvent
pas être calculés de manière simple. Deux approches pragmatiques sont présentées dans ce cha-
pitre, qui permettent d’obtenir des estimateurs sous–optimaux, c’est–à–dire qui n’atteignent
pas nécessairement le minimum de l’erreur quadratique moyenne, mais qui sont néanmoins très
largement utilisés en pratique. La première approche présentée à la Section 4.1 repose sur des
45
46 CHAPITRE 4. EXTENSIONS AUX SYSTÈMES NON–LINÉAIRES
Xk = bk (Xk−1 ) + σk (Xk−1 ) Wk ,
(4.3)
Yk = hk (Xk ) + Vk ,
et on suppose que les fonctions bk et hk sont dérivables. En linéarisant le système (4.3) autour
d’une suite déterministe donnée, ou bien autour de l’estimateur courant, on peut obtenir des
algorithmes sous–optimaux, qui sont décrits ci–dessous.
On se donne une suite (déterministe) {x̄k } à valeurs dans Rm , appelée trajectoire nominale (on
peut prendre par exemple x̄k comme une approximation de la moyenne de Xk ). La méthode
consiste à linéariser les fonctions bk et σk autour de x̄k−1 , c’est–à–dire
Le système non–linéaire (4.3) est alors remplacé par le système linéaire gaussien
Yk = HkL Xk + hLk + Vk ,
avec
FkL = b′k (x̄k−1 ) et fkL = −b′k (x̄k−1 ) x̄k−1 + bk (x̄k−1 ) ,
et avec
HkL = h′k (x̄k ) et hLk = −h′k (x̄k ) x̄k + hk (x̄k ) .
4.1. FILTRE DE KALMAN LINÉARISÉ, FILTRE DE KALMAN ÉTENDU 47
Ici, le vecteur aléatoire WkL = σk (x̄k−1 ) Wk est gaussien, centré et de matrice de covariance
QLk = σk (x̄k−1 ) σk∗ (x̄k−1 ). On applique alors exactement le filtre de Kalman à ce nouveau système,
et on obtient l’algorithme sous–optimal suivant
b − = bk (x̄k−1 ) + b′ (x̄k−1 ) (X
X bk−1 − x̄k−1 ) ,
k k
Pk− = b′k (x̄k−1 ) Pk−1 (b′k (x̄k−1 ))∗ + σk (x̄k−1 ) σk∗ (x̄k−1 ) ,
et
bk = X
X b − + Kk (Yk − (hk (x̄k ) + h′ (x̄k ) (X
b − − x̄k ))) ,
k k k
Kk = Pk− (h′k (x̄k ))∗ [ h′k (x̄k ) Pk− (h′k (x̄k ))∗ + QVk ]−1 .
bk = X
X b − + Kk (Yk − hk (X
b − )) .
k k
Au lieu de linéariser autour d’une trajectoire nominale déterministe {x̄k }, on peut utiliser l’es-
bk−1 , c’est–
timateur courant. La méthode consiste à linéariser les fonctions bk et σk autour de X
à–dire
bk−1 ) + b′ (X
bk (x) ≃ bk (X bk−1 ) (x − X
bk−1 ) et bk−1 ) ,
σk (x) ≃ σk (X
k
b − , c’est–à–dire
et à linéariser la fonction hk autour de X k
b − ) + h′ (X
hk (x) ≃ hk (X b − ) (x − X
b −) .
k k k k
Le système non–linéaire (4.3) est alors remplacé par le système conditionnellement linéaire
gaussien
Yk = HkL Xk + hLk + Vk ,
avec
bk−1 )
FkL = b′k (X et bk−1 ) X
fkL = −b′k (X bk−1 + bk (X
bk−1 ) ,
48 CHAPITRE 4. EXTENSIONS AUX SYSTÈMES NON–LINÉAIRES
et avec
b −)
HkL = h′k (X et b −) X
hLk = −h′k (X b − + hk (X
b −) ,
k k k k
et on remarque que
bk−1 + f L = bk (X
FkL X bk−1 ) et b − + hL = hk (X
HkL X b −) .
k k k k
et
bk = X
X b − + Kk (Yk − hk (X
b − )) ,
k k
b − )) P − ,
Pk = (I − Kk h′k (X k k
Remarque 4.1 Dans cet algorithme, la suite {Pk } dépend des observations, et ne peut donc
pas être pré–calculée.
Xk = bk (Xk−1 ) + σk (Xk−1 ) Wk ,
Yk = hk (Xk ) + Vk ,
et on ne suppose plus que les fonctions bk et hk sont dérivables, mais on suppose que les fonctions
bk , hk et σk et certaines fonctions associées, peuvent être intégrées par rapport à certaines
distributions de probabilité gaussiennes.
Au lieu de s’appuyer sur une linéarisation des fonctions autour de l’estimateur courant, on
se propose ici
4.2. FILTRE DE KALMAN UNSCENTED 49
Le premier point peut s’interpréter comme une projection, au sens de la distance de Kullback–
Leibler, sur la famille des distributions de probabilité gaussiennes.
▶ Le calcul des deux premiers moments (moyenne et matrice de covariance) de la distribution
de probabilité conditionnelle µ− k (dx) = P[Xk ∈ dx | Y0:k−1 ], c’est–à–dire le calcul de la moyenne
conditionnelle et de la matrice de covariance conditionnelle du vecteur aléatoire Xk sachant
Y0:k−1 , est facile. Par définition
b − = E[Xk | Y0:k−1 ]
X k
de sorte que
b − ) (Xk − X
Pk− = E[ (Xk − X b − )∗ | Y0:k−1 ]
k k
b − ) (bk (Xk−1 ) − X
= E[ (bk (Xk−1 ) − X b − )∗ | Y0:k−1 ]
k k
b − )∗ | Y0:k−1 ]
+ E[σk (Xk−1 ) Wk (bk (Xk−1 ) − X k
b − ) W ∗ σ ∗ (Xk−1 ) | Y0:k−1 ]
+ E[ (bk (Xk−1 ) − X k k k
b − )∗ | Y0:k−1 ] = 0 ,
= E[σk (Xk−1 ) E[Wk | Xk−1 , Y0:k−1 ] (bk (Xk−1 ) − X k
où on a encore utilisé dans la dernière égalité l’indépendance de (Y0 , · · · , Yk−1 , Xk−1 ) et de Wk ,
donc E[Wk | Xk−1 , Y0:k−1 ] = 0.
▶ En revanche, le calcul des deux premiers moments (moyenne et matrice de covariance) de
la distribution de probabilité conditionnelle µk (dx) = P[Xk ∈ dx | Y0:k ], c’est–à–dire le calcul
de la moyenne conditionnelle et de la matrice de covariance conditionnelle du vecteur aléatoire
Xk sachant Y0:k , n’est pas immédiat, et on commence par le calcul des deux premiers moments
(moyenne et matrice de covariance) de la distribution de probabilité conditionnelle jointe du
vecteur aléatoire (Xk , Yk ) sachant Y0:k−1 , qui est plus facile. On rappelle que
∫
b −
Xk = bk (x) µk−1 (dx) ,
Rm
On rappelle que
∫ ∫
Pk− = b − ) (bk (x) − X
(bk (x) − X b − )∗ µk−1 (dx) + σk (x) σk∗ (x) µk−1 (dx) ,
k k
Rm
de sorte que
= E[ (hk (Xk ) − Ybk− ) (hk (Xk ) − Ybk− )∗ | Y0:k−1 ] + E[Vk Vk∗ | Y0:k−1 ]
b − ) (Yk − Yb − )∗ | Y0:k−1 ]
Ck = E[ (Xk − Xk k
∫
= b − ) (hk (x) − Yb − )∗ µ− (dx) .
(x − X k k k
Rm
respectivement. Si la matrice QVk est inversible, alors a fortiori la matrice Ξk est inversible, et
d’après la Proposition 1.7 on obtient immédiatement les approximations suivantes
bk = X
X b − + Ck Ξ−1 (Yk − Yb − ) et Pk = Pk− − Ck Ξ−1 ∗
k k k k Ck ,
52 CHAPITRE 4. EXTENSIONS AUX SYSTÈMES NON–LINÉAIRES
Remarque 4.2 Si on suppose que les fonctions bk et hk sont dérivables, et qu’on utilise un
développement limité au premier ordre au voisinage de u = 0 dans les intégrales ci–dessus, on
retrouve les équations du filtre de Kalman étendu. L’idée ici est de ne pas linéariser , et de
calculer les intégrales en utilisant des formules de quadrature numérique.
c’est–à–dire que les deux premiers moments sont pris en compte exactement. Plus généralement
∫ ∑
+m
du
ϕ(u) exp{− 12 |u|2 } ≈ wi ϕ(ui ) ,
Rm (2π)m/2 i=−m
pour toute fonction ϕ définie sur Rm , c’est–à–dire que les σ–points (x−m , · · · , xm ) associés à la
distribution de probabilité gaussienne de vecteur moyenne µ et de matrice de covariance Σ, sont
définis par la relation xi = µ + Σ1/2 ui , soit
√ √
x0 = µ , xi = µ + Σ1/2 ei m + κ et x−i = µ − Σ1/2 ei m + κ ,
pour tout i = 1, · · · , m. On vérifie que
∑
+m ∑
+m ∑
m
∗
wi x i = µ et wi (xi − µ) (xi − µ) = Σ1/2 ei (Σ1/2 ei )∗ = Σ ,
i=−m i=−m i=1
c’est–à–dire que les deux premiers moments sont pris en compte exactement. Plus généralement
encore, soit X un vecteur aléatoire gaussien de vecteur moyenne µ et de matrice de covariance
Σ, et soit T une transformation non–linéaire définie sur Rm . Clairement
∫ ∑
+m
du
ϕ(T (µ + Σ1/2 u)) exp{− 21 |u|2 } ≈ wi ϕ(T (xi )) ,
(2π)m/2 i=−m
54 CHAPITRE 4. EXTENSIONS AUX SYSTÈMES NON–LINÉAIRES
pour toute fonction ϕ définie sur Rm , c’est–à–dire que les σ–points (x′−m , · · · , x′m ) associés au
vecteur aléatoire transformé X ′ = T (X), sont simplement obtenus par la relation x′i = T (xi ) à
partir des σ–points (x−m , · · · , xm ) associés au vecteur aléatoire X, soit
√ √
x′0 = T (µ) , x′i = T (µ + Σ1/2 ei m + κ) et x′−i = T (µ − Σ1/2 ei m + κ) ,
pour tout i = 1, · · · , m.
Avec ces formules de quadrature, on obtient l’algorithme de filtrage sous–optimal suivant.
b − et P − en fonction de X
Expression de X bk−1 et Pk−1 = Sk−1 S ∗ :
k k k−1
On introduit les σ–points
√ √
bk−1 ,
x0 = X bk−1 + Sk−1 ei
xi = X m+κ et bk−1 − Sk−1 ei
x−i = X m+κ ,
et la matrice de covariance
∑
+m ∑
+m
Pk− = b − ) (bk (xi ) − X
wi (bk (xi ) − X b − )∗ + wi σk (xi ) σk∗ (xi ) = Sk− (Sk− )∗ .
k k
i=−m i=−m
bk et Pk en fonction de X
Expression de X b − et P − = S − (S − )∗ :
k k k k
On introduit les σ–points
√ √
b− ,
x0 = X b − + S − ei
xi = X m+κ et b − − S − ei
x−i = X m+κ ,
k k k k k
la matrice de covariance
∑
+m
Ξk = wi (hk (xi ) − Ybk− ) (hk (xi ) − Ybk− )∗ + QVk ,
i=−m
et la matrice de corrélation
∑
+m
Ck = b − ) (hk (xi ) − Yb − )∗ ,
wi (xi − X k k
i=−m
et on pose
bk = X
X b − + Ck Ξ−1 (Yk − Yb − ) et Pk = Pk− − Ck Ξ−1 ∗ ∗
k k k k Ck = Sk Sk .
Chapitre 5
Il s’agit de la classe la plus générale de modèles d’état, et c’est aussi un cas particulier de la
classe plus générale des modèles de Markov cachés (pour lesquels l’espace d’état peut être très
général). On considère donc une suite d’états cachés {Xk } à valeurs dans Rm , vérifiant
Xk = fk (Xk−1 , Wk ) avec W k ∼ pW
k (dw) , (5.1)
avec des entrées bruitées {Wk } à valeurs dans Rp , pas nécessairement gaussiennes, et une condi-
tion initiale X0 ∼ η0 (dx) pas nécessairement gaussienne, et une suite d’observations {Yk } à
valeurs dans Rd , vérifiant
avec des bruits d’observation {Vk } additifs, à valeurs dans Rd , pas nécessairement gaussiens,
mais de distribution de probabilité qkV (v) dv absolument continue par rapport à la mesure de
Lebesgue dv. Les bruits blancs {Wk } et {Vk } sont indépendants entre eux et indépendants de
la condition initiale X0 . On ne suppose pas que les fonctions fk et hk sont dérivables. Pour la
suite, il sera suffisant de faire l’hypothèse suivante : pour tout instant k
• la distribution de probabilité du vecteur aléatoire Vk admet une densité qkV (v) qu’il est
facile d’évaluer pour tout v ∈ Rd .
Proposition 5.1 La suite {Xk } est une chaı̂ne de Markov à valeurs dans Rm , c’est–à–dire que
la distribution de probabilité conditionnelle par rapport au passé
55
56 CHAPITRE 5. AU–DELÀ DES SYSTÈMES LINÉAIRES GAUSSIENS
défini par ∫
Qk ϕ(x) = E[ ϕ(Xk ) | Xk−1 = x] = ϕ(fk (x, w)) pW
k (dw) ,
Rp
pour toute fonction test ϕ mesurable bornée, définie sur Rm .
pour toute fonction ϕ mesurable bornée définie sur Rm . Clairement, le résultat ne dépend que
de Xk−1 , c’est–à–dire que
et ∫
E[ ϕ(Xk ) | Xk−1 = x] = ϕ(fk (x, w)) pW
k (dw) . 2
Rp
Qk (x, dx′ ) = pW ′ ′
k (x − bk (x)) dx ,
c’est–à–dire que le noyau Qk (x, dx′ ) admet une densité. En effet, le changement de variable
x′ = bk (x) + w donne immédiatement
∫ ∫
Qk ϕ(x) = W
ϕ(bk (x) + w) pk (w) dw = ϕ(x′ ) pW ′ ′
k (x − bk (x)) dx ,
Rm Rm
Remarque 5.3 En général, le noyau Qk (x, dx′ ) n’admet pas de densité par rapport à la mesure
de Lebesgue sur Rm . En effet, conditionnellement à Xk−1 = x, le vecteur aléatoire Xk appartient
nécessairement au sous–ensemble
et dans le cas où p < m ce sous–ensemble M(x) est généralement, sous certaines hypothèses de
régularité, une sous–variété différentielle de dimension p dans l’espace Rm , c’est–à–dire un sous–
ensemble de mesure de Lebesgue nulle. La distribution de probabilité conditionnelle du vecteur
aléatoire Xk sachant Xk−1 = x ne peut donc pas avoir de densité par rapport à la mesure de
Lebesgue sur Rm .
5.1. SYSTÈMES NON–LINÉAIRES À BRUITS NON–GAUSSIENS 57
Proposition 5.4 La suite {Yk } vérifie l’hypothèse de canal sans mémoire, c’est–à–dire que pour
tout instant n
• conditionnellement aux états cachés X0:n les observations Y0:n sont mutuellement indépen-
dantes, ce qui se traduit par
∏
n
P[Y0:n ∈ dy0:n | X0:n ] = P[Yk ∈ dyk | X0:n ] ,
k=0
Exemple 5.5 Dans le cas particulier où le bruit additif Vk est un vecteur aléatoire gaussien
centré et de matrice de covariance identité, alors la probabilité d’émission
1
P[Yk ∈ dy | Xk = x] = exp{− 21 |y − hk (x)|2 } dy ,
(2 π)d/2
n ∫
∏
= ϕk (hk (Xk ) + vk ) qkV (vk ) dvk
k=0 Rd
n ∫
∏ ∏
n
= ϕk (yk ) qkV (yk − hk (Xk )) dyk = E[ϕk (Yk ) | Xk ] . 2
k=0 Rd | {z } k=0
P[Yk ∈ dyk | Xk ]
Yk = hk (rk , Xk ) + Vk ,
où la suite {rk } forme une chaı̂ne de Markov à valeurs dans un espace fini, correspondant à
différents régimes ou modes de fonctionnement.
Plus généralement, on peut aussi considérer un modèle de Markov caché où les états cachés
{Xk } forment une chaı̂ne de Markov à valeurs dans un espace E qui peut être très général, par
exemple un espace hybride continu / discret, un sous–ensemble défini par des contraintes, une
variété différentielle, un graphe, etc., de noyaux de transition
P[Xk ∈ dx′ | Xk−1 = x] = Qk (x, dx′ ) ,
et de distribution de probabilité initiale
P[X0 ∈ dx] = η0 (dx) ,
et où les observations {Yk } vérifient l’hypothèse de canal sans mémoire, c’est–à–dire que pour
tout instant n
5.2. MODÈLES DE MARKOV CACHÉS 59
• conditionnellement aux états cachés X0:n les observations Y0:n sont mutuellement indépen-
dantes, ce qui se traduit par
∏
n
P[Y0:n ∈ dy0:n | X0:n ] = P[Yk ∈ dyk | X0:n ] ,
k=0
où la mesure positive λFk (dy) définie sur F ne dépend pas de l’état caché x ∈ E, et par
abus de notation on définit la fonction de vraisemblance
gk (x) = gk (x, Yk ) ,
? ? ?
Yk−1 Yk Yk+1
où les flèches représentent la dépendance entre variables aléatoires. En d’autres termes, la dis-
tribution de probabilité conditionnelle jointe des observations Y0:n sachant les états cachés X0:n
vérifie
∏
n
P[Y0:n ∈ dy0:n | X0:n = x0:n ] = gk (xk , yk ) λF0 (dy0 ) · · · λFn (dyn ) .
k=0
Ce modèle peut paraı̂tre très abstrait à première vue, mais pour la suite il suffira que l’hypothèse
suivante soit vérifiée : pour tout instant k = 1, · · · , n
Malgré leur grande généralité, les modèles de Markov cachés ne permettent pas de prendre en
compte un certain nombre de systèmes non–linéaires à bruits non–gaussiens, qui correspondent
60 CHAPITRE 5. AU–DELÀ DES SYSTÈMES LINÉAIRES GAUSSIENS
à des situations d’intérêt pratique, par exemple les cas où les observations dépendent de la
transition de la chaı̂ne de Markov cachée
Xk = fk (Xk−1 , Wk ) ,
Yk = hk (Xk−1 , Xk ) + Vk ,
les modèles de Markov à variables latentes
Xk = fk (Xk−1 , Wk ) ,
Yk = hk (Yk−1 , Xk ) + Vk ,
où conditionnellement aux états cachés, les observations forment une chaı̂ne de Markov, ou bien
les systèmes d’état à bruits corrélés
Xk = fk (Xk−1 , Vk−1 , Wk ) ,
Yk = hk (Xk ) + Vk ,
où clairement le bruit Uk−1 = (Vk−1 , Wk ) dans l’équation d’état est corrélé au bruit d’observation
Vk−1 . Dans ce dernier exemple, une solution pragmatique consiste à reporter dans l’équation
d’état l’expression pour Vk−1 tirée de l’équation d’observation, de sorte que
Xk = fk (Xk−1 , Yk−1 − hk−1 (Xk−1 ), Wk ) ,
Yk = hk (Xk ) + Vk .
Les classes de modèles de plus en plus généraux présentés dans les deux prochaines sections,
permettent de prendre en compte ces situations.
Certains problèmes sont décrits par une chaı̂ne de Markov {Yk }, et pour disposer d’une plus
grande possibilité de modélisation on propose de faire dépendre les noyaux de transitions d’une
suite de variables aléatoires latentes {Xk }, formant elle–même une chaı̂ne de Markov. Cette
situation se rencontre par exemple dans les modèles à volatilité stochastique, et à la différence
de la situation précédente, l’estimation de la suite latente n’est pas un objectif en soi.
Dans ces modèles, les états cachés {Xk } forment une chaı̂ne de Markov à valeurs dans un
espace E, de noyaux de transition
P[Xk ∈ dx′ | Xk−1 = x] = Qk (x, dx′ ) ,
et de distribution de probabilité initiale
P[X0 ∈ dx] = η0 (dx) ,
et conditionnellement aux états cachés, les observations {Yk } forment une chaı̂ne de Markov à
valeurs dans F , c’est–à–dire que pour tout instant n
5.3. CHAÎNES DE MARKOV À PARAMÈTRES MARKOVIENS 61
• conditionnellement aux états cachés X0:n les observations Y0:n forment une chaı̂ne de
Markov, ce qui se traduit pour tout k = 1, · · · , n, par
où la mesure positive λF0 (dy) définie sur F ne dépend pas de l’état caché x ∈ E, et par
abus de notation on définit la fonction de vraisemblance
g0 (x) = g0 (x, Y0 ) ,
où la mesure positive λFk (y, dy ′ ) définie sur F ne dépend pas de l’état caché x′ ∈ E, et par
abus de notation on définit la fonction de vraisemblance
qui mesure l’adéquation d’un état quelconque x′ ∈ E avec les observations successives Yk−1
et Yk .
- Xk−1 - Xk - Xk+1 -
? ? ?
- Yk−1 - Yk - Yk+1 -
62 CHAPITRE 5. AU–DELÀ DES SYSTÈMES LINÉAIRES GAUSSIENS
où les flèches représentent la dépendance entre variables aléatoires. En d’autres termes, la dis-
tribution de probabilité conditionnelle jointe des observations Y0:n sachant les états cachés X0:n
vérifie
∏
n
= g0 (x0 , y0 ) λF0 (dy0 ) gk (xk , yk−1 , yk ) λFk (yk−1 , dyk )
k=1
∏
n ∏
n
= [ g0 (x0 , y0 ) gk (xk , yk−1 , yk ) ] λF0 (dy0 ) λFk (yk−1 , dyk ) .
k=1 k=1
Encore plus généralement, on peut considérer un modèle où les états cachés {Xk } ne forment
plus nécessairement une chaı̂ne de Markov, mais où conjointement états cachés et observations
{Zk } avec Zk = (Xk , Yk ) pour tout instant k = 0, 1, · · · , n, forment une chaı̂ne de Markov à
valeurs dans E × F , de distribution de probabilité initiale
où la mesure positive λF0 (dy) définie sur F , ne dépend pas de l’état caché x ∈ E, et de probabilités
de transition
où la mesure positive λFk (y, dy ′ ) définie sur F , dépend de l’observation précédente y ∈ F mais ne
dépend pas de la transition cachée (x, x′ ) ∈ E. En d’autres termes, la distribution de probabilité
jointe des états cachés X0:n et des observations Y0:n vérifie
∏
n
= γ0 (y0 , dx0 ) λF0 (dy0 ) Rk (yk−1 , yk , xk−1 , dxk ) λFk (yk−1 , dyk )
k=1
∏
n ∏
n
= [ γ0 (y0 , dx0 ) Rk (yk−1 , yk , xk−1 , dxk ) ] λF0 (dy0 ) λFk (yk−1 , dyk ) .
k=1 k=1
En toute généralité, les mesures positives γ0 (y, dx) et les noyaux positifs Rk (y, y ′ , x, dx′ ) peuvent
être factorisés comme
γ0 (y, dx) = g0imp (y, x) η0imp (y, dx) et Rk (y, y ′ , x, dx′ ) = gkimp (y, y ′ , x, x′ ) Qimp ′ ′
k (y, y , x, dx ) ,
• et d’une distribution de probabilité η0imp (y, dx) ou d’un noyau markovien Qimp ′ ′
k (y, y , x, dx ).
Une telle factorisation n’est évidemment pas unique, mais il existe toujours au moins la factori-
sation donnée par
γ0 (y, dx)
γ0 (y, dx) = γ0 (y, E) ,
γ (y, E)
| {z } | 0 {z }
gb0 (y) ηb0 (y, dx)
et
Rk (y, y ′ , x, dx′ )
Rk (y, y ′ , x, dx′ ) = Rk (y, y ′ , x, E) ,
Rk (y, y ′ , x, E)
| {z } | {z }
gbk (x, y, y ′ ) b ′
Qk (y, y , x, dx )′
Dans le cas particulier des modèles de Markov cachés, cette décomposition fait intervenir
• la probabilité d’émission
• et la transition de probabilité
b k (y ′ , x, dx′ ) ,
P[Xk ∈ dx′ | Xk−1 = x, Yk = y ′ ] = Q
R R R R
Yk−1 Yk Yk+1
Exemple 5.6 On considère un système non–linéaire avec des bruits gaussiens additifs et une
fonction d’observation linéaire
Xk = fk (Xk−1 ) + σk (Xk−1 ) Wk ,
Yk = Hk Xk + hk + Vk ,
où la condition initiale X0 est un vecteur aléatoire gaussien de moyenne X̄0 et de matrice de
covariance QX 0 , et où les suites {Wk } et {Vk } sont des bruits blancs gaussiens indépendants,
indépendants de la condition initiale X0 , de matrices de covariance identité et QVk respective-
ment, avec QVk inversible. Il résulte de la Proposition 1.7 que
X ′ = X + QX ∗ X ∗ V −1
0 H0 [H0 Q0 H0 + Q0 ] (Y0 − (H0 X + h0 + V )) .
5.4. CHAÎNES DE MARKOV PARTIELLEMENT OBSERVÉES 65
Xk = fk (Xk−1 ) + σk (Xk−1 ) Wk ,
Yk = Hk fk (Xk−1 ) + hk + Hk σk (Xk−1 ) Wk + Vk ,
d’où on déduit que conditionnellement à Xk−1 = x, le vecteur aléatoire (Xk , Yk ) est gaussien,
de moyenne et de matrice de covariance
fk (x) Σk (x) Σk (x) Hk∗
et ,
Hk fk (x) + hk ∗
Hk Σk (x) Hk Σk (x) Hk + Qk V
respectivement, avec Σk (x) = σk (x) σk∗ (x). Compte tenu que la matrice QVk est inversible, la
matrice Hk Σk (x) Hk∗ + QVk est inversible a fortiori, et il résulte de la Proposition 1.7 que
mk (y ′ , x) = fk (x) + Σk (x) Hk∗ [Hk Σk (x) Hk∗ + QVk ]−1 (y ′ − (Hk fk (x) + hk )) ,
Pk (x) = Σk (x) − Σk (x) Hk∗ [Hk Σk (x) Hk∗ + QVk ]−1 Hk Σk (x) ,
Pour évaluer la performance des algorithmes numériques de filtrage non–linéaire, y compris les
nombreuses variantes du filtrage particulaire, il est utile de disposer d’une borne inférieure sur
l’erreur commise par un estimateur quelconque de l’état caché. S’il s’agit d’estimer un paramètre
fixe, il est bien connu que la matrice d’information de Fisher associée au modèle statistique
permet d’obtenir une telle borne inférieure, sous le nom de borne de Cramér–Rao. Dans le cas
du filtrage bayésien, il s’agit d’estimer un paramètre aléatoire (et dynamique), à savoir la suite
des états cachés, pour lequel on dispose d’un modèle a priori : dans ce cadre bayésien, on peut
utiliser la notion de borne de Cramér–Rao a posteriori, pour laquelle des algorithmes de calcul
récursifs efficaces ont été obtenus.
On considère le modèle général d’une chaı̂ne de Markov partiellement observée, et on suppose
qu’il existe
et se ramener au problème statique considéré dans la Proposition 1.3 ci–dessus pour l’estimation
du vecteur aléatoire
ϕ(X0:n ) = Xn ,
sachant Y0:n .
Théorème 6.1 Sous les hypothèses de la Proposition 1.3, la matrice de corrélation de l’erreur
d’estimation (ψ(Y0:n ) − Xn ) est minorée par la relation suivante
67
68 CHAPITRE 6. BORNE DE CRAMÉR–RAO A POSTERIORI
avec
∂2
Dk− = −E[ log rk (Yk−1 , Yk , Xk−1 , Xk ) ] ,
∂x2k−1
∂2
Dk = −E[ log rk (Yk−1 , Yk , Xk−1 , Xk ) ] ,
∂xk−1 ∂xk
∂2
Dk+ = −E[ log rk (Yk−1 , Yk , Xk−1 , Xk ) ] .
∂x2k
Preuve. La densité jointe des vecteurs aléatoires X0:n et Y0:n est donnée par
∏
n
p0:n (x0:n , y0:n ) = r0 (x0 , y0 ) rk (yk−1 , yk , xk−1 , xk ) ,
k=1
d’où la log–densité
∑
n
log p0:n (x0:n , y0:n ) = log r0 (x0 , y0 ) + log rk (yk−1 , yk , xk−1 , xk )
k=1
On en déduit que
∂2 ∂2
∂x2 ∂x0:n−1 ∂xn
∂2 0:n−1
log p0:n (x0:n , y0:n ) =
log p0:n (x0:n , y0:n )
∂x20:n 2
∂
⋆
∂x2n
∂2 ∂2 ∂2
∂x2 ∂x0:n−2 ∂xn
0:n−2 ∂x0:n−2 ∂xn−1
[ log p0:n−1 (x0:n−1 , y0:n−1 )
∂2 ∂ 2
= ⋆
∂x2n−1 ∂xn−1 ∂xn
+ log rn (yn−1 , yn , xn−1 , xn ) ]
∂2
⋆ ⋆
∂x2n
∂2 ∂2
∂x2 0
0:n−2 ∂x0:n−2 ∂xn−1
= ∂2 log p0:n−1 (x0:n−1 , y0:n−1 )
⋆ 0
∂x2n−1
0 0 0
0 0 0
∂2 ∂2
0
+ ∂x2n−1 ∂xn−1 ∂xn log rn (yn−1 , yn , xn−1 , xn )
∂2
0 ⋆
∂x2n
An−1 Bn−1 0 0 0 0 An−1 Bn−1 0
J0:n =
⋆ Cn−1 0
+0 Dn− Dn
= ⋆ Cn−1 + Dn− Dn
, (6.2)
0 0 0 0 ⋆ +
Dn 0 ⋆ +
Dn
70 CHAPITRE 6. BORNE DE CRAMÉR–RAO A POSTERIORI
avec
∂2
Dn− = −E[ log rn (Yn−1 , Yn , Xn−1 , Xn ) ] ,
∂x2n−1
∂2
Dn = −E[ log rn (Yn−1 , Yn , Xn−1 , Xn ) ] ,
∂xn−1 ∂xn
∂2
Dn+ = −E[ log rn (Yn−1 , Yn , Xn−1 , Xn ) ] .
∂x2n
On remarque que ϕ′ (x0:n ) = (0 I) = Mn ne dépend pas de x0:n , et il résulte de la Proposi-
tion 1.3 que
−1
E[ (ψ(Y0:n ) − Xn ) (ψ(Y0:n ) − Xn )∗ ] ≥ Mn J0:n Mn∗ ,
et d’après le Lemme A.3 d’inversion matricielle et la Remarque A.4, on a
−1
( ) An Bn 0 ( ) ⋆ ⋆ 0
−1 ∗
Mn J0:n Mn = 0 I = 0 I = Jn−1 ,
Bn∗ Cn I ⋆ Jn−1 I
∗
où la matrice ∆n−1 = Cn−1 + Dn− − Bn−1 A−1 −
n−1 Bn−1 = Jn−1 + Dn est le complément de Schur
de la matrice An−1 dans la matrice–bloc An , de sorte que
r0 (y, x) = p0 (x) q0 (y | x) ,
et
rk (y, y ′ , x, x′ ) = pk (x′ | x) qk (y ′ | x′ ) ,
et dans ce cas
∂2
Dk− = −E[ log pk (Xk | Xk−1 ) ] ,
∂x2k−1
∂2
Dk = −E[ log pk (Xk | Xk−1 ) ] ,
∂xk−1 ∂xk
∂2 ∂2
Dk+ = −E[ log p k (X k | X k−1 ) ] − E[ log qk (Yk | Xk ) ] .
∂x2k ∂x2k
Exemple 6.2 Dans le cas particulier d’un système avec bruits gaussiens additifs, où
et où
Yk = hk (Xk ) + Vk avec Vk ∼ N(0, QVk ) ,
avec des matrices de covariance QW V
k et Qk inversibles, on obtient
−1
Dk+ = (QW
k ) + E[ [h′k (Xk )]∗ (QVk )−1 h′k (Xk ) ] .
72 CHAPITRE 6. BORNE DE CRAMÉR–RAO A POSTERIORI
Chapitre 7
Filtrage bayésien
L’objectif de ce chapitre est d’établir les équations du filtre non–linéaire, pour les systèmes
non–linéaires et non–gaussiens, ou plus généralement les équations du filtre bayésien, pour les
modèles de Markov cachés et les chaı̂nes de Markov partiellement observées. Il s’agit donc de
calculer la distribution de probabilité conditionnelle de la variable aléatoire Xk sachant Y0:k , et
la distribution de probabilité conditionnelle de la variable aléatoire Xk sachant Y0:k−1 , définies
par
µk (dx) = P[Xk ∈ dx | Y0:k ] et µ−
k (dx) = P[Xk ∈ dx | Y0:k−1 ] ,
respectivement.
P[X0:n ∈ dx0:n , Y0:n ∈ dy0:n ] = P[Y0:n ∈ dy0:n | X0:n = x0:n ] P[X0:n ∈ dx0:n ]
∏
n
= P[X0:n ∈ dx0:n ] gk (xk , yk ) λF0 (dy0 ) · · · λFn (dyn ) .
k=0
En intégrant par rapport aux variables x0:n , on obtient la distribution de probabilité jointe des
observations Y0:n , c’est–à–dire
∫ ∫ ∏
n
P[Y0:n ∈ dy0:n ] = ··· gk (xk , yk ) P[X0:n ∈ dx0:n ] λF0 (dy0 ) · · · λFn (dyn )
E E k=0
∏
n
= E[ gk (Xk , yk ) ] λF0 (dy0 ) · · · λFn (dyn ) .
k=0
73
74 CHAPITRE 7. FILTRAGE BAYÉSIEN
∏
n
= P[X0:n ∈ dx0:n ] gk (xk , yk ) λF0 (dy0 ) · · · λFn (dyn )
k=0
∏
n
= P[X0:n ∈ dx0:n | Y0:n = y0:n ] E[ gk (Xk , yk ) ] λF0 (dy0 ) · · · λFn (dyn ) ,
k=0
et on obtient
∏
n
gk (xk , yk ) P[X0:n ∈ dx0:n ]
k=0
P[X0:n ∈ dx0:n | Y0:n = y0:n ] = ,
∏
n
E[ gk (Xk , yk ) ]
k=0
pour toute suite y0:n d’observations. Pour toute fonction test f définie sur l’espace produit
E n+1 = E × · · · × E, on a
∫ ∫ ∏
n
··· f (x0:n ) gk (xk , yk ) P[X0:n ∈ dx0:n ]
E E k=0
E[f (X0:n ) | Y0:n = y0:n ] =
∏
n
E[ gk (Xk , yk ) ]
k=0
∏
n
E[ f (X0:n ) gk (Xk , yk ) ]
k=0
= ,
∏
n
E[ gk (Xk , yk ) ]
k=0
et on rappelle que
∏
n
P[Y0:n ∈ dy0:n ] = E[ gk (Xk , yk ) ] λF0 (dy0 ) · · · λFn (dyn ) ,
k=0
et comme ces identités sont vérifiées pour toute suite y0:n d’observations, on en déduit que la
distribution de probabilité conditionnelle jointe des états cachés X0:n sachant Y0:n est donnée
par
∏
n
E[f (X0:n ) gk (Xk ) ]
k=0
E[f (X0:n ) | Y0:n ] = , (7.1)
∏
n
E[ gk (Xk ) ]
k=0
7.1. MODÈLES DE MARKOV CACHÉS 75
où l’espérance porte seulement sur la suite des états cachés X0:n : les fonctions de vraisemblance
g0 (x), · · · , gn (x) sont définies par abus de notation comme
gk (x) = gk (x, Yk ) ,
et pour ϕ ≡ 1, on a
∏
n
⟨γn , 1⟩ = E[ gk (Xk ) ] = Ln .
k=0
De la même manière, la distribution de probabilité conditionnelle de l’état caché Xn sachant
Y0:n−1 est donnée par
∏
n−1
E[ϕ(Xn ) gk (Xk ) ]
⟨γn− , ϕ⟩
⟨µ−
n , ϕ⟩ = E[ϕ(Xn ) | Y0:n−1 ] = k=0
= ,
∏
n−1 ⟨γn− , 1⟩
E[ gk (Xk ) ]
k=0
∏
n−1
⟨γn− , ϕ⟩ = E[ϕ(Xn ) gk (Xk ) ] ,
k=0
et pour ϕ ≡ 1, on a
∏
n−1
⟨γn− , 1⟩ = E[ gk (Xk ) ] = ⟨γn−1 , 1⟩ .
k=0
Pour obtenir une équation récurrente permettant d’exprimer µk en fonction de µk−1 , il suffit
donc d’une équation récurrente permettant d’exprimer γk en fonction de γk−1 , puis de normaliser.
76 CHAPITRE 7. FILTRAGE BAYÉSIEN
Théorème 7.1 (Filtre bayésien) La suite {µk } vérifie l’équation récurrente suivante
prédiction correction
µk−1 −−−−−−−−−−→ µ−k = µ k−1 Qk −
− −−−−−−−−→ µk = gk · µ−
k ,
désigne l’action du noyau markovien Qk (x, dx′ ) sur la distribution de probabilité µk−1 (dx), et
la notation
gk µ−
gk · µ− = k
,
k
⟨µ−
k , g k⟩
et en particulier pour ϕ ≡ 1
⟨µk−1 Rk , 1⟩ = ⟨µ−
k , gk ⟩ ,
et en normalisant, on vérifie que
⟨µk−1 Rk , ϕ⟩ ⟨µ− , gk ϕ⟩
= k− = ⟨µk , ϕ⟩ .
⟨µk−1 Rk , 1⟩ ⟨µk , gk ⟩
Expression de µ−
n en fonction de µn−1 :
On remarque immédiatement que
∏
n−1
⟨γn− , 1⟩ = E[ gk (Xk ) ] = ⟨γn−1 , 1⟩ ,
k=0
∏
n−1
= E[ E[ϕ(Xn ) | X0:n−1 ] gk (Xk ) ]
k=0
∏
n−1
= E[ E[ϕ(Xn ) | Xn−1 ] gk (Xk ) ]
k=0
∏
n−1
= E[Qn ϕ(Xn−1 ) gk (Xk )] = ⟨γn−1 , Qn ϕ⟩ = ⟨γn−1 Qn , ϕ⟩ ,
k=0
pour toute fonction test ϕ définie sur E, où la dernière égalité exprime simplement que
∫
⟨γn−1 , Qn ϕ⟩ = γn−1 (dx) Qn ϕ(x)
E
∫ ∫ ∫ ∫
′ ′
= γn−1 (dx) [ Qn (x, dx ) ϕ(x ) ] = [ γn−1 (dx) Qn (x, dx′ ) ] ϕ(x′ )
E E E E
∫
= γn−1 Qn (dx′ ) ϕ(x′ ) = ⟨γn−1 Qn , ϕ⟩ .
E
Comme la fonction test ϕ est quelconque, on en déduit que
γn− = γn−1 Qn ,
et en normalisant, on obtient
γn− γn−1 Qn
µ−
n = − = = µn−1 Qn .
⟨γn , 1⟩ ⟨γ n−1 , 1⟩
78 CHAPITRE 7. FILTRAGE BAYÉSIEN
Expression de µn en fonction de µ−
n :
On a simplement
∏
n
⟨γn , ϕ⟩ = E[ϕ(Xn ) gk (Xk ) ]
k=0
∏
n−1
= E[ϕ(Xn ) gn (Xn ) gk (Xk ) ] = ⟨γn− , gn ϕ⟩ = ⟨gn γn− , ϕ⟩ ,
k=0
pour toute fonction test ϕ définie sur E, où la dernière égalité exprime simplement que
∫ ∫
⟨γn− , gn ϕ⟩ = [gn (x) ϕ(x)] γn− (dx) = ϕ(x) [gn (x) γn− (dx)] = ⟨gn γn− , ϕ⟩ .
E E
γn = gn γn− ,
et en normalisant, on obtient
γn gn γ − gn µ−
µn = = − n = − n ,
⟨γn , 1⟩ ⟨γn , gn ⟩ ⟨µn , gn ⟩
où la dernière égalité est obtenue en divisant numérateur et dénominateur par la constante de
normalisation ⟨γn− , 1⟩. 2
D’après la propriété de Markov, la distribution de probabilité jointe des états cachés X0:n et des
observations Y0:n vérifie
∏
n ∏
n
= [ γ0 (y0 , dx0 ) Rk (yk−1 , yk , xk−1 , dxk ) ] λF0 (dy0 ) λFk (yk−1 , dyk ) .
k=1 k=1
En intégrant par rapport aux variables x0:n , on obtient la distribution de probabilité jointe des
observations Y0:n , c’est–à–dire
P[Y0:n ∈ dy0:n ]
∫ ∫ ∏
n ∏
n
=[ ··· γ0 (y0 , dx0 ) Rk (yk−1 , yk , xk−1 , dxk ) ] λF0 (dy0 ) λFk (yk−1 , dyk ) .
E E k=1 k=1
7.2. CHAÎNES DE MARKOV PARTIELLEMENT OBSERVÉES 79
∏
n ∏
n
= [ γ0 (y0 , dx0 ) Rk (yk−1 , yk , xk−1 , dxk ) ] λF0 (dy0 ) λFk (yk−1 , dyk )
k=1 k=1
∫ ∫ ∏
n ∏
n
[ ··· γ0 (y0 , dx0 ) Rk (yk−1 , yk , xk−1 , dxk ) ] λF0 (dy0 ) λFk (yk−1 , dyk ) ,
E E k=1 k=1
et on obtient
∏
n
γ0 (y0 , dx0 ) Rk (yk−1 , yk , xk−1 , dxk )
k=1
P[X0:n ∈ dx0:n | Y0:n = y0:n ] = ∫ ∫ ,
∏
n
··· γ0 (y0 , dx0 ) Rk (yk−1 , yk , xk−1 , dxk )
E E k=1
pour toute suite y0:n d’observations. Pour toute fonction test f définie sur l’espace produit
E n+1 = E × · · · × E, on a
∫ ∫ ∏
n
··· γ0 (y0 , dx0 ) Rk (yk−1 , yk , xk−1 , dxk ) f (x0:n )
E E k=1
E[f (X0:n ) | Y0:n = y0:n ] = ∫ ∫ ,
∏
n
··· γ0 (y0 , dx0 ) Rk (yk−1 , yk , xk−1 , dxk )
E E k=1
et on rappelle que
P[Y0:n ∈ dy0:n ]
∫ ∫ ∏
n ∏
n
=[ ··· γ0 (y0 , dx0 ) Rk (yk−1 , yk , xk−1 , dxk ) ] λF0 (dy0 ) λFk (yk−1 , dyk ) .
E E k=1 k=1
et comme ces identités sont vérifiées pour toute suite y0:n d’observations, on en déduit que la
distribution de probabilité conditionnelle jointe des états cachés X0:n sachant Y0:n est donnée
par
∫ ∫ ∏
n
··· γ0 (dx0 ) Rk (xk−1 , dxk ) f (x0:n )
E E k=1
E[f (X0:n ) | Y0:n ] = ∫ ∫ , (7.2)
∏
n
··· γ0 (dx0 ) Rk (xk−1 , dxk )
E E k=1
80 CHAPITRE 7. FILTRAGE BAYÉSIEN
où la mesure positive γ0 (dx) et les noyaux positifs (non–normalisés) Rk (x, dx′ ) sont définis par
abus de notation comme
pour tout k = 1, · · · , n, et dépendent implicitement des observations Y0:n , mais celles–ci sont
considérées comme fixées dans les expressions ci–dessus. En particulier, la distribution de pro-
babilité conditionnelle de l’état caché Xn sachant Y0:n est donnée par
∫ ∫ ∏
n
··· γ0 (dx0 ) Rk (xk−1 , dxk ) ϕ(xn )
E E k=1 ⟨γn , ϕ⟩
⟨µn , ϕ⟩ = E[ϕ(Xn ) | Y0:n ] = ∫ ∫ = ,
∏
n
⟨γn , 1⟩
··· γ0 (dx0 ) Rk (xk−1 , dxk )
E E k=1
et pour ϕ ≡ 1, on a
∫ ∫ ∏
n
⟨γn , 1⟩ = ··· γ0 (dx0 ) Rk (xk−1 , dxk ) = Ln .
E E k=1
Pour obtenir une équation récurrente permettant d’exprimer µk en fonction de µk−1 , il suffit
donc d’une équation récurrente permettant d’exprimer γk en fonction de γk−1 , puis de normaliser.
Théorème 7.4 (Filtre bayésien) La suite {µk } vérifie l’équation récurrente suivante
µk−1 Rk
µk = ,
⟨µk−1 Rk , 1⟩
∏
n
Lk = ⟨µk−1 Rk , 1⟩ Lk−1 soit en itérant Ln = ⟨µk−1 Rk , 1⟩ ⟨γ0 , 1⟩ .
k=1
7.2. CHAÎNES DE MARKOV PARTIELLEMENT OBSERVÉES 81
Preuve. On a
∫ ∫ ∏
n
⟨γn , ϕ⟩ = ··· γ0 (dx0 ) Rk (xk−1 , dxk ) ϕ(xn )
E E k=1
∫ ∫ ∏
n−1 ∫
= ··· γ0 (dx0 ) Rk (xk−1 , dxk ) Rn (xn−1 , dxn ) ϕ(xn )
E E k=1 E
∫ ∫ ∏
n−1
= ··· γ0 (dx0 ) Rk (xk−1 , dxk ) Rn ϕ(xn−1 )
E E k=1
= ⟨γn−1 , Rn ϕ⟩ = ⟨γn−1 Rn , ϕ⟩ .
γn = γn−1 Rn ,
et en normalisant, on obtient
γn γn−1 Rn µn−1 Rn
µn = = = ,
⟨γn , 1⟩ ⟨γn−1 Rn , 1⟩ ⟨µn−1 Rn , 1⟩
où la dernière égalité est obtenue en divisant numérateur et dénominateur par la constante de
normalisation ⟨γn−1 , 1⟩. 2
On a déjà vu à la Section 5.4 que les mesures positives γ0 (y, dx) et les noyaux positifs
Rk (y, y ′ , x, dx′ ) peuvent être factorisés comme
γ0 (y, dx) = g0imp (y, x) η0imp (y, dx) et Rk (y, y ′ , x, dx′ ) = gkimp (y, y ′ , x, x′ ) Qimp ′ ′
k (y, y , x, dx ) ,
• et d’une distribution de probabilité η0imp (y, dx) ou d’un noyau markovien Qimp ′ ′
k (y, y , x, dx ),
et avec les abus de notation habituels, une telle décomposition implique que la mesure positive
γ0 (dx) et le noyau positif Rk (x, dx′ ) peuvent être factorisés comme
γ0 (dx) = g0imp (x) η0imp (dx) et Rk (x, dx′ ) = gkimp (x, x′ ) Qimp ′
k (x, dx ) , (7.3)
• d’une fonction de pondération positive, éventuellement aléatoire, g0imp (x) = g0imp (Y0 , x) ou
gkimp (x, x′ ) = gkimp (Yk−1 , Yk , x, x′ ),
• et d’une distribution de probabilité, éventuellement aléatoire, η0imp (dx) = η0imp (Y0 , dx) ou
d’un noyau markovien, éventuellement aléatoire, Qimp ′ imp ′
k (x, dx ) = Qk (Yk−1 , Yk , x, dx ).
82 CHAPITRE 7. FILTRAGE BAYÉSIEN
Cette décomposition est évidemment non unique : dans le cas particulier des modèles de Markov
cachés, le premier exemple d’une telle décomposition est donné naturellement par la définition
peut être interprétée pour tout état x ∈ E comme une mesure quantitative du recouvrement
entre l’application x′ 7→ gk (x′ ) et la distribution de probabilité Qk (x, dx′ ). En pratique, la
décomposition d’importance doit être telle que
• il est facile de simuler une variable aléatoire selon la distribution de probabilité η0imp (dx),
• il est facile de simuler pour tout x ∈ E une variable aléatoire selon la distribution de
probabilité Qimp ′
k (x, dx ),
quand bien même l’expression analytique du noyau positif Rk (x, dx′ ) serait inconnue, ou telle-
ment compliquée qu’il serait impossible en pratique de calculer des intégrales telles que
∫ ∫
′ ′ ′
Rk ϕ(x) = Rk (x, dx ) ϕ(x ) ou µ Rk (dx ) = µ(dx) Rk (x, dx′ ) .
E E
Exemple 7.5 Cette situation favorable se rencontre par exemple pour le système non–linéaire,
présenté dans l’Exemple 5.6, avec des bruits gaussiens additifs et une fonction d’observation
linéaire
Xk = fk (Xk−1 ) + σk (Xk−1 ) Wk ,
(7.4)
Yk = Hk Xk + hk + Vk ,
où les suites {Wk } et {Vk } sont des bruits blancs gaussiens indépendants, indépendants de
la condition initiale X0 , de matrices de covariance identité et QVk respectivement (avec QVk
inversible) à l’instant k. Dans ce cas en effet
7.2. CHAÎNES DE MARKOV PARTIELLEMENT OBSERVÉES 83
• il est facile de simuler un vecteur aléatoire X ′ selon la distribution de probabilité ηb0 (dx) :
il suffit de simuler deux vecteurs aléatoires gaussiens indépendants X et V , de moyenne
X̄0 et 0 et de matrice de covariance QX V
0 et Q0 respectivement, et de poser
X ′ = X + QX ∗ X ∗ V −1
0 H0 [H0 Q0 H0 + Q0 ] (Y0 − (H0 X + h0 + V )) ,
X = fk (x) + σk (x) W ,
et
X ′ = X + Σk (x) Hk∗ [Hk Σk (x) Hk∗ + QVk ]−1 (Yk − (Hk X + hk + V )) .
L’équation du filtre bayésien a été obtenue très simplement, mais il est en général impossible
de la résoudre, sauf dans le cas particulier des systèmes linéaires gaussiens, où elle se ramène aux
équations du filtre de Kalman, présentées au Chapitre 3. Il faut donc avoir recours à une approxi-
mation numérique, et on présente ci–dessous une approximation de type Monte Carlo, appelée
filtre particulaire, qui a connu un développement spectaculaire au cours des dernières années,
et qui est maintenant largement répendu, en particulier dans les applications en localisation,
navigation ou poursuite de mobiles, aussi bien dans le domaine militaire (aéronef, sous–marin,
bâtiment de surface, missile, drone, etc.), que dans le domaine civil, avec des applications en
robotique mobile ou en communications sans–fil.
84 CHAPITRE 7. FILTRAGE BAYÉSIEN
Chapitre 8
Généralisation : distributions de
Feynman–Kac
et où gk (x) sont des fonctions mesurables bornées (strictement positives) données, appelées
fonctions de sélection ou fonctions de fitness, pour tout k = 0, 1, · · · , n. L’hypothèse minimale,
faute de quoi le problème n’est pas bien posé, est que ⟨γn , 1⟩ > 0, ce qui est assuré par exemple
si les fonctions de sélection sont strictement positives.
∏
k−1
= E[ϕ(Xk ) gk (Xk ) gp (Xp ) ] = ⟨γk− , gk ϕ⟩ = ⟨gk γk− , ϕ⟩ ,
p=0
85
86 CHAPITRE 8. GÉNÉRALISATION : DISTRIBUTIONS DE FEYNMAN–KAC
∏
k−1
⟨γk− , ϕ⟩ = E[ϕ(Xk ) gp (Xp ) ]
p=0
∏
k−1
= E[ E[ϕ(Xk ) | X0:k−1 ] gp (Xp ) ]
p=0
∏
k−1
= E[ E[ϕ(Xk ) | Xk−1 ] gp (Xp ) ]
p=0
∏
k−1
= E[ Qk ϕ(Xk−1 ) gp (Xp ) ] = ⟨γk−1 , Qk ϕ⟩ = ⟨γk−1 Qk , ϕ⟩ ,
p=0
pour toute fonction mesurable bornée ϕ, de sorte que la distribution non–normalisée vérifie la
relation de récurrence linéaire
de sorte que la distribution normalisée vérifie la relation de récurrence non–linéaire décrite par
le schéma suivant
mutation pondération
µk−1 −−−−−−−−−−→ ηk = µk−1 Qk −−−−−−−−−−−−→ µk = gk · ηk ,
∏
n ∏
n
⟨γn , 1⟩ = E[ gk (Xk ) ] = ⟨ηk , gk ⟩ ,
k=0 k=0
c’est–à–dire que l’espérance d’un produit est remplacée par un produit d’espérances.
On remarque que la distribution non–normalisée vérifie aussi la relation de récurrence linéaire
c’est–à–dire que l’identité (8.5) est vérifiée pour k = (n − 1). D’autre part, si l’identité (8.5) est
vérifiée à l’étape k, alors
gk−1 (x) (Rk:n−1 Qn )(x, dx′ ) = gk−1 (x) (Rk (Rk+1:n−1 Qn ))(x, dx′ )
∫
= gk−1 (x) Rk (x, dx′′ ) (Rk+1:n−1 Qn )(x′′ , dx′ )
E
∫
= gk−1 (x) Qk (x, dx′′ ) gk (x′′ ) (Rk+1:n−1 Qn )(x′′ , dx′ )
E
∫
−
= gk−1 (x) Qk (x, dx′′ ) Rk+1:n (x′′ , dx′ )
E
∫
= Rk− (x, dx′′ ) Rk+1:n
−
(x′′ , dx′ )
E
−
= Rk:n (x, dx′ ) ,
c’est–à–dire que l’identité (8.5) est vérifiée à l’étape (k − 1). En particulier, il résulte de (8.5)
que
−
gk Rk+1:n−1 1 = gk Rk+1:n−1 Qn 1 = Rk+1:n 1, (8.6)
pour tout k = 0, 1, · · · , (n − 1).
η0 (dx) = r0 (x) η00 (dx) et Qk (x, dx′ ) = rk (x, x′ ) Q0k (x, dx′ ) , (8.7)
∏
n
= E[ϕ(Xn0 ) gk0 (Xk−1
0
, Xk0 ) ] ,
k=0
88 CHAPITRE 8. GÉNÉRALISATION : DISTRIBUTIONS DE FEYNMAN–KAC
pour toute fonction mesurable bornée ϕ, où la suite {Xk0 , k = 0, 1, · · · , n} est une chaı̂ne de
Markov, caractérisée par
pour toute fonction mesurable bornée ϕ, où γ0 est une mesure positive donnée et où {Rk , k =
1, · · · , n} sont des noyaux positifs (non–normalisés) donnés. Ce modèle général inclut comme
cas particulier le modèle (8.1), avec
γ0 (dx) = η0 (dx) g0 (x) et Rk (x, dx′ ) = Qk (x, dx′ ) gk (x′ ) .
∫ ∫ ∏
k−1 ∫
= ··· γ0 (dx0 ) Rp (xp−1 , dxp ) Rk (xk−1 , dxk ) ϕ(xk )
E E p=1 E
∫ ∫ ∏
k−1
= ··· γ0 (dx0 ) Rp (xp−1 , dxp ) Rk ϕ(xk−1 )
E E p=1
= ⟨γk−1 , Rk ϕ⟩ = ⟨γk−1 Rk , ϕ⟩ ,
pour toute fonction mesurable bornée ϕ, de sorte que la distribution non–normalisée vérifie la
relation de récurrence linéaire
γk = γk−1 Rk = µk−1 Rk ⟨γk−1 , 1⟩ . (8.9)
8.2. MODÈLE (APPARAMMENT) PLUS GÉNÉRAL 89
de sorte que la distribution normalisée vérifie la relation de récurrence non–linéaire décrite par
le schéma suivant
µk−1 Rk
µk−1 −−−−−−−−−→ µk = ,
⟨µk−1 Rk , 1⟩
c’est–à–dire qu’une intégrale multiple est remplacée par un produit d’intégrales doubles.
γ0 (dx) = g0imp (x) η0imp (dx) et Rk (x, dx′ ) = gkimp (x, x′ ) Qimp ′
k (x, dx ) , (8.11)
et la représentation probabiliste
∏
n
⟨γn , ϕ⟩ = E[ϕ(Xnimp ) gkimp (Xk−1
imp
, Xkimp ) ] , (8.12)
k=0
pour toute fonction mesurable bornée ϕ, où la suite {Xkimp , k = 0, 1, · · · , n} est une chaı̂ne de
Markov caractérisée par
avec la convention g0imp (x, x′ ) = g0imp (x′ ) pour k = 0. On remarque que chaque fonction de
sélection dépend de la transition courante (et pas seulement de l’état courant) de la chaı̂ne de
Markov.
La décomposition (8.11) est évidemment non unique : dans le cas particulier du modèle
considéré à la Section 8.1, le premier exemple d’une telle décomposition est donné naturellement
par la définition
dépend seulement de l’état de départ de la transition courante, et peut être interprétée pour tout
état x ∈ E comme une mesure quantitative du recouvrement entre l’application x′ 7→ gk (x′ ) et
la distribution de probabilité Qk (x, dx′ ), pour tout k = 1, · · · , n, et plus généralement, il existe
une décomposition
γ0 (dx) = g0 (x) r0 (x) η00 (dx) et Rk (x, dx′ ) = gk (x′ ) rk (x, x′ ) Q0k (x, dx′ ) ,
| {z } | {z }
′
g00 (x) 0
gk (x, x )
• il est facile de simuler une variable aléatoire selon la distribution de probabilité η0imp (dx),
• il est facile de simuler pour tout x ∈ E une variable aléatoire selon la distribution de
probabilité Qimp ′
k (x, dx ),
8.2. MODÈLE (APPARAMMENT) PLUS GÉNÉRAL 91
quand bien même l’expression analytique du noyau positif Rk (x, dx′ ) serait inconnue, ou telle-
ment compliquée qu’il serait impossible en pratique de calculer des intégrales telles que
∫ ∫
Rk ϕ(x) = Rk (x, dx′ ) ϕ(x′ ) où µ Rk (dx′ ) = µ(dx) Rk (x, dx′ ) .
E E
pour tout k = 1, · · · , n, décomposition qui devient triviale dans le cas particulier où la fonction
de sélection g0 (x) = cste est constante, et où la fonction de sélection gk (x, x′ ) = gk (x) ne
dépend que de l’état de départ de la transition courante, pour tout k = 1, · · · , n, on obtient la
représentation probabiliste équivalente
∏
n
⟨γn , ϕ⟩ = E[ϕ(Xn ) gk (Xk−1 , Xk ) ]
k=0
∫ ∫ ∏
n ∏
n
= ··· g0 (x0 ) η0 (dx0 ) gk (xk−1 , xk ) Qk (xk−1 , dxk ) ϕ(xn )
E E k=1 k=1
∫ ∫ ∏
n ∏
n
= ⟨η0 , g0 ⟩ ··· ηb0 (dx0 ) gbk (xk−1 ) b k (xk−1 , dxk ) ϕ(xn )
Q
E E k=1 k=1
∫ ∫ ∏
n−1 ∏
n
= ⟨η0 , g0 ⟩ ··· ηb0 (dx0 ) gbk+1 (xk ) b k (xk−1 , dxk ) ϕ(xn )
Q
E E k=0 k=1
∫ ∫ ∏
n−1 ∏
n
= ⟨η0 , g0 ⟩ ··· η0opt (dx0 ) gkopt (xk ) Qopt
k (xk−1 , dxk ) ϕ(xn )
E E k=0 k=1
∏
n−1
= ⟨η0 , g0 ⟩ E[ϕ(Xnopt ) gkopt (Xkopt ) ]
k=0
= ⟨η0 , g0 ⟩ ⟨γnopt− , ϕ⟩ ,
où la suite {Xkopt , k = 0, 1, · · · , n} est une chaı̂ne de Markov, caractérisée par
pour tout k = 1, · · · , n,
et où les fonctions de sélection gkopt (x′ ) sont définies par (8.14), c’est–à–dire
∫
opt ′ ′
gk (x ) = gbk+1 (x ) = gk+1 (x′ , x′′ ) Qk+1 (x′ , dx′′ ) , (8.17)
E
pour tout k = 0, 1, · · · , (n − 1). En pratique, cette décomposition n’est vraiment utile que si
• il est facile de simuler une variable aléatoire selon la distribution de probabilité ηb0 (dx),
• il est facile de simuler pour tout x ∈ E une variable aléatoire selon la distribution de
probabilité Qb k (x, dx′ ),
⟨γn , ϕ⟩ ⟨γnopt− , ϕ⟩
⟨µn , ϕ⟩ = = opt− = ⟨ηnopt , ϕ⟩ ,
⟨γn , 1⟩ ⟨γn , 1⟩
c’est–à–dire que la constante de normalisation ⟨γn , 1⟩ pour le modèle d’origine peut s’interpréter
opt
en terme de la constante de normalisation ⟨γn−1 , 1⟩ à l’instant précédent pour le modèle dit
8.3. MODÈLE À VALEURS TRANSITIONS OU TRAJECTOIRES 93
pour tout k = 1, · · · , n, et on déduit de (8.5) l’identité suivante entre noyaux positifs (non–
normalisés)
gkopt (Rk+1:n−1
opt
Qopt opt−
n ) = Rk+1:n = Rk+1:n , (8.20)
valide pour tout k = (n − 1), · · · , 1, 0.
Remarque 8.1 A titre de vérification, et sans repasser par l’identité (8.5), on peut montrer de
manière directe l’identité (8.20) par récurrence arrière. D’après (8.19) on a immédiatement
opt ′ ′
gn−1 (x) Qopt
n (x, dx ) = Rn (x, dx ) ,
c’est–à–dire que l’identité (8.20) est vérifiée pour k = (n − 1). D’autre part, si l’identité (8.20)
est vérifiée à l’étape k, alors d’après (8.19) on a
opt opt ′ opt opt opt ′
gk−1 (x) (Rk:n−1 Qopt opt
n )(x, dx ) = gk−1 (x) (Rk (Rk+1:n−1 Qn ))(x, dx )
∫
opt
= gk−1 (x) Rkopt (x, dx′′ ) (Rk+1:n−1
opt
Qopt ′′ ′
n )(x , dx )
E
∫
opt ′′ opt ′′ ′′ ′
= gk−1 (x) Qopt opt opt
k (x, dx ) gk (x ) (Rk+1:n−1 Qn )(x , dx )
E
∫
opt ′′ ′′ ′
= gk−1 (x) Qopt
k (x, dx ) Rk+1:n (x , dx )
E
∫
= Rk (x, dx′′ ) Rk+1:n (x′′ , dx′ )
E
c’est–à–dire que l’identité (8.20) est vérifiée à l’étape (k − 1). En particulier, il résulte de (8.20)
que
gkopt Rk+1:n−1
opt
1 = gkopt Rk+1:n−1
opt
Qopt
n 1 = Rk+1:n 1 , (8.21)
pour tout k = (n − 1), · · · , 1, 0.
Tous les modèles présentés jusqu’ici semblent pouvoir être vus comme des cas particuliers du
modèle
∏ n
⟨γn , ϕ⟩
⟨γn , ϕ⟩ = E[ϕ(Xn ) gk (Xk−1 , Xk ) ] et ⟨µn , ϕ⟩ = , (8.22)
⟨γn , 1⟩
k=0
pour toute fonction mesurable bornée ϕ, où {Xk , k = 0, 1, · · · , n} est une chaı̂ne de Markov
caractérisée par
et où gk (x, x′ ) sont des fonctions mesurables bornées (strictement positives) données, appelées
fonctions de sélection ou fonctions de fitness, pour tout k = 0, 1, · · · , n, avec la convention
g0 (x, x′ ) = g0 (x′ ) pour k = 0.
On remarque que chaque fonction de sélection dépend de la transition courante de la chaı̂ne
de Markov, ce qui inclus comme cas particuliers le cas où chaque fonction de sélection dépend
seulement de l’état d’arrivée de la transition courante, comme dans le modèle (8.1), et le cas où
chaque fonction de sélection dépend seulement de l’état de départ de la transition courante.
Le modèle (8.22) peut être vu comme un cas particulier du modèle (8.8), avec
∏
n
⟨γn , ϕ⟩ = E[ϕ(Xn ) gk (Xk−1 , Xk ) ]
k=0
∫ ∫ ∏
n ∏
n
= ··· ϕ(xn ) gk (xk−1 , xk ) η0 (dx0 ) Qk (xk−1 , dxk )
E E k=0 k=1
∫ ∫ ∏
n
= ··· ϕ(xn ) γ0 (dx0 ) Rk (xk−1 , dxk ) ,
E E k=1
pour toute fonction mesurable bornée ϕ, et on déduit de (8.9) que la distribution non–normalisée
vérifie la relation de récurrence linéaire
de sorte que
∫ ∫
⟨µ Rk , ϕ⟩ = [ gk (x, x′ ) µ(dx) Qk (x, dx′ ) ] ϕ(x′ )
E E
∫ ∫
= ϕ(x′ ) gk (x, x′ ) µ(dx) Qk (x, dx′ )
E E
= ⟨µ ⊗ Qk , gk ϕ ◦ π⟩ ,
µ Rk = (gk (µ ⊗ Qk )) ◦ π −1 .
∏
n ∏
n
⟨γn , 1⟩ = E[ gk (Xk−1 , Xk ) ] = ⟨η0 , g0 ⟩ ⟨µk−1 ⊗ Qk , gk ⟩ ,
k=0 k=1
c’est–à–dire que l’espérance d’un produit est remplacée par un produit d’espérances.
Remarque 8.2 Pour générer une variable aléatoire (X, X ′ ) distribuée selon (µ ⊗ Qk )(dx, dx′ ),
il suffit de générer d’abord une variable aléatoire X distribuée selon µ(dx), et de générer ensuite
une variable aléatoire X ′ distribuée selon Qk (X, dx′ ).
96 CHAPITRE 8. GÉNÉRALISATION : DISTRIBUTIONS DE FEYNMAN–KAC
= ⟨η0:n , g0:n ϕ ◦ π⟩
On remarque que
⟨γn , ϕ⟩ ⟨η0:n , g0:n ϕ ◦ π⟩
⟨µn , ϕ⟩ = = ,
⟨γn , 1⟩ ⟨η0:n , g0:n ⟩
pour toute fonction mesurable bornée ϕ, de sorte que la distribution normalisée µn s’exprime en
terme de la distribution de Gibbs–Boltzmann
g0:n η0:n
µ0:n = g0:n · η0:n = ,
⟨η0:n , g0:n ⟩
définie sur l’espace trajectoriel E n+1 = E × · · · × E, comme µn = µ0:n ◦ π −1 .
pour tout k = 1, · · · , n,
c’est–à–dire que l’état de départ de la nouvelle transition coı̈ncide avec l’état d’arrivée de la
transition précédente et l’état d’arrivée de la nouvelle transition est distribué à partir de l’état
de départ selon le noyau de transition du modèle (8.22), et on considère la distribution non–
normalisée et la distribution normalisée associée, définies par
∏
n
⟨γntr , f ⟩
⟨γntr , f ⟩ = E[f (Xntr ) gk (Xktr ) ] et ⟨µtr
n , f⟩ = , (8.27)
⟨γntr , 1⟩
k=0
pour toute fonction mesurable bornée f définie sur l’ensemble produit E tr = E × E. On vérifie
que
∏
n ∏
n
⟨γn , ϕ⟩ = E[ϕ(Xn ) gk (Xk−1 , Xk ) ] = E[ϕ ◦ π(Xntr ) gk (Xktr ) ] = ⟨γntr , ϕ ◦ π⟩ ,
k=0 k=0
c’est–à–dire que l’espérance d’un produit est remplacée par un produit d’éspérances.
Remarque 8.3 A titre de vérification et sans repasser par les représentations probabilistes as-
sociées (8.22) et (8.27), on peut montrer de manière directe la consistance des deux suites définies
par les relations de récurrence (8.24) et (8.28) respectivement, c’est–à–dire vérifier par récurrence
98 CHAPITRE 8. GÉNÉRALISATION : DISTRIBUTIONS DE FEYNMAN–KAC
que γktr ◦ π −1 = γk , pour tout k = 1, · · · , n. On suppose que l’hypothèse de récurrence est vraie
au rang (k − 1), c’est–à–dire que la distribution marginale (non–normalisée) de γk−1 tr coincide
tr
avec la distribution non–normalisée γk−1 , ou en d’autres termes γk−1 (E, dx2 ) = γk−1 (dx2 ), et
en utilisant la relation de récurrence (8.28) on remarque que
⟨γktr , ϕ ◦ π⟩ = ⟨γk−1
tr
k , ϕ ◦ π gk ⟩
Qtr
∫ ∫ ∫ ∫
= tr
γk−1 (dx1 , dx2 ) δx2 (dx′1 ) Qk (x′1 , dx′2 ) ϕ(x′2 ) gk (x′1 , x′2 )
E E E E
∫ ∫ ∫
= tr
γk−1 (dx1 , dx2 ) Qk (x2 , dx′2 ) ϕ(x′2 ) gk (x2 , x′2 )
E E E
∫ ∫
= tr
γk−1 (E, dx2 ) Qk (x2 , dx′2 ) ϕ(x′2 ) gk (x2 , x′2 )
E E
∫ ∫
= γk−1 (dx2 ) Qk (x2 , dx′2 ) ϕ(x′2 ) gk (x2 , x′2 )
E E
∫ ∫
= γk−1 (dx2 ) Qk (x2 , dx′2 ) ϕ(x′2 ) gk (x2 , x′2 )
E E
∫ ∫
= (γk−1 ⊗ Qk )(dx2 , dx′2 ) ϕ(x′2 ) gk (x2 , x′2 )
E E
= ⟨γk−1 ⊗ Qk , gk ϕ ◦ π⟩
en utilisant la relation de récurrence (8.24), c’est–à–dire que l’hypothèse de récurrence est vraie
au rang k.
▶ Chaı̂ne de Markov à valeurs trajectoires Plus généralement encore, tous les modèles
présentés jusqu’ici peuvent être vus comme des cas particuliers du modèle trajectoriel suivant.
On définit la variable aléatoire Xk• = X0:k = (X0 , · · · , Xk ) à valeurs dans l’espace produit
E k+1 = E ×· · ·×E dépendant du temps, pour tout k = 0, 1, · · · , n. Clairement la suite {Xk• , k =
0, 1, · · · , n} est une chaı̂ne de Markov, caractérisée par
Q•k (x0 , · · · , xk−1 , dx′0 , · · · , dx′k ) = δx0 (dx′0 ) · · · δxk−1 (dx′k−1 ) Qk (xk−1 , dx′k ) ,
pour tout k = 1, · · · , n,
8.3. MODÈLE À VALEURS TRANSITIONS OU TRAJECTOIRES 99
pour toute fonction mesurable bornée f définie sur l’ensemble produit E n+1 = E × · · · × E, où
les fonctions de sélection gk• (x0:k ) sont définies par
∏
n ∏
n
⟨γn , ϕ⟩ = E[ϕ(Xn ) gk (Xk−1 , Xk ) ] = E[ϕ ◦ π(Xn• ) gk• (Xk• ) ] = ⟨γn• , ϕ ◦ π⟩ ,
k=0 k=0
Pour une distribution de probabilité µ donnée, il s’agit d’approcher numériquement, par des
méthodes de Monte Carlo, l’intégrale ou de manière équivalente l’espérance mathématique
∫
⟨µ, ϕ⟩ = ϕ(x) µ(dx) = E[ϕ(X)] , (9.1)
E
où la variable aléatoire X a pour distribution de probabilité µ, pour toute fonction mesurable
bornée ϕ. Dans toute la suite, la notation S N (µ) désigne la distribution de probabilité empirique
1 ∑
N
S N (µ) = δ i
N ξ
i=1
1 ∑
N
⟨S N (µ), ϕ⟩ = ϕ(ξ i ) ,
N
i=1
1 ∑
N
⟨µ, ϕ⟩ ≈ ⟨S (µ), ϕ⟩ =
N
ϕ(ξ i ) ,
N
i=1
1 ∑
N
µ ≈ S (µ) =
N
δ i ,
N ξ
i=1
101
102 CHAPITRE 9. MÉTHODES DE MONTE CARLO
Théorème 9.1 La variable aléatoire ⟨S N (µ), ϕ⟩ est un estimateur non–biaisé de ⟨µ, ϕ⟩, et les
moments de l’erreur d’estimation vérifient
1
E| ⟨S N (µ) − µ, ϕ⟩ |2 = var(ϕ, µ) ,
N
et pour tout réel p ≥ 2, il existe une constante positive cp > 0 telle que
cp
{ E| ⟨S N (µ) − µ, ϕ⟩ |p }1/p ≤ √ ⟨µ, |ϕ − ⟨µ, ϕ⟩|p ⟩1/p ,
N
pour toute fonction mesurable bornée ϕ.
1 ∑
N
E| ⟨S (µ) − µ, ϕ⟩ | = E|
N 2
[ϕ(ξ i ) − ⟨µ, ϕ⟩ ] |2
N
i=1
1 ∑
N
1
= 2
E|ϕ(ξ i ) − ⟨µ, ϕ⟩|2 = ⟨µ, |ϕ − ⟨µ, ϕ⟩|2 ⟩ ,
N N
i=1
pour toute fonction mesurable bornée ϕ. Plus généralement, pour tout réel p ≥ 2
1 ∑
N
E| ⟨S N (µ) − µ, ϕ⟩ |p = E| [ϕ(ξ i ) − ⟨µ, ϕ⟩ ] |p
N
i=1
Bp 1 ∑
N
Bp
≤ p/2 E|ϕ(ξ i ) − ⟨µ, ϕ⟩|p = p/2 ⟨µ, |ϕ − ⟨µ, ϕ⟩|p ⟩ ,
N N N
i=1
Une approche traditionnelle pour calculer l’intégrale (9.1) est l’échantillonnage pondéré, ou
importance sampling, dans laquelle une nouvelle distribution de probabilité ν ≫ µ est utilisée,
qui domine la distribution de probabilité µ, c’est–à–dire qu’il existe une densité (ou dérivée de
Radon–Nikodym) dµ/dν telle que
∫ ∫
dµ dµ dµ
⟨µ, ϕ⟩ = ϕ(x) µ(dx) = ϕ(x) (x) ν(dx) = ⟨ν, ϕ ⟩ = E[ϕ(Ξ) (Ξ)] ,
E E dν dν dν
où la variable aléatoire Ξ a pour distribution de probabilité ν, pour toute fonction mesurable
bornée ϕ. S’il est facile
1 ∑
N
dµ dµ dµ
⟨µ, ϕ⟩ = ⟨ν, ϕ ⟩ ≈ ⟨S (ν), ϕ
N
⟩= ϕ(xi ) (xi ) ,
dν dν N dν
i=1
où les variables aléatoires (x1 , · · · , xN ) sont i.i.d., de distribution de probabilité commune ν,
c’est–à–dire que
1 ∑ dµ
N
dµ dµ N
µ= ν ≈ µN = S (ν) = (xi ) δxi .
dν dν N dν
i=1
1 ∑ dµ
N
⟨µN , 1⟩ = (xi ) ,
N dν
i=1
n’est pas nécessairement égale à 1, de sorte que l’approximation µN n’est pas nécessairement
normalisée. En revanche, il résulte immédiatement du Théorème 9.1 que la variable aléatoire
dµ dµ
⟨µN , ϕ⟩ = ⟨S N (ν), ϕ ⟩ est un estimateur non–biaisé de ⟨ν, ϕ ⟩ = ⟨µ, ϕ⟩, et la variance (non–
dν dν
asymptotique) de cet estimateur est
dµ 2 1 dµ
E|⟨µN − µ, ϕ⟩|2 = E⟨S N (ν) − ν, ϕ ⟩ = var(ϕ , ν) .
dν N dν
104 CHAPITRE 9. MÉTHODES DE MONTE CARLO
pour toute fonction mesurable bornée ϕ, où les variables aléatoires (x1 , · · · , xN ) sont i.i.d., de
distribution de probabilité commune ν. Il résulte de la décomposition suivante
µN
µ′N − µ = − µ = (µN − µ) − ⟨µN − µ, 1⟩ µ′N ,
⟨µN , 1⟩
que
⟨µ′N − µ, ϕ⟩ = ⟨µN − µ, ϕ⟩ − ⟨µN − µ, 1⟩ ⟨µ′N , ϕ⟩ ,
de sorte que
{ E|⟨µ′N − µ, ϕ⟩|2 }1/2 ≤ { E|⟨µN − µ, ϕ⟩|2 }1/2 + { E|⟨µN − µ, 1⟩|2 }1/2 ∥ϕ∥ ,
d’après l’inégalité (triangulaire) de Minkowski, pour toute fonction mesurable bornée ϕ. Ce type
d’approximation sera étudié en détail à la Section 9.2.
On remarque que
∫ ∫
dµ dµ dµ
var(ϕ , ν) = (ϕ(x) (x)) ν(dx) − (
2
ϕ(x) (x) ν(dx) )2
dν E dν E dν
∫
dµ
= (ϕ(x) (x))2 ν(dx) − ⟨µ, ϕ⟩2 ,
E dν
et ∫ ∫
dµ dµ
(ϕ(x) (x)) ν(dx) ≥ (
2
|ϕ(x)| (x) ν(dx) )2 = ⟨µ, |ϕ|⟩2 ,
E dν E dν
d’après l’inégalité de Jensen, d’où la borne inférieure suivante
dµ
var(ϕ , ν) ≥ ⟨µ, |ϕ|⟩2 − ⟨µ, ϕ⟩2 ≥ 0 ,
dν
indépendante du choix de la distribution de probabilité d’importance ν. On remarque que si la
fonction ϕ garde un signe constant, alors la borne inférieure est nulle.
Parmi tous les choix possibles pour la distribution de probabilité d’importance ν, il existe
un choix qui minimise la variance, c’est–à–dire que la borne inférieure est atteinte, même si ce
choix est en pratique inaccessible car il nécessite de connaı̂tre la constante de normalisation
∫
⟨µ, |ϕ|⟩ = |ϕ(x)| µ(dx) = E|ϕ(X)| ,
E
dont le calcul présente le même degré de difficulté que le calcul de l’intégrale (9.1) elle–même !
En effet, si on introduit
|ϕ| µ
ν∗ = |ϕ| · µ = ,
⟨µ, |ϕ|⟩
9.1. ÉCHANTILLONNAGE PONDÉRÉ 105
Remarque 9.5 Il est certainement possible de simuler une variable aléatoire distribuée selon
ν∗ , même si la constante de normalisation ⟨µ, |ϕ|⟩ est inconnue, en utilisant l’une ou l’autre des
méthodes proposées à la Section 9.2, mais la connaissance explicite de la constante de norma-
lisation est absolument nécessaire pour évaluer la fonction d’importance optimale dµ/dν∗ . En
pratique, le choix optimal n’est donc simplement pas utilisable. Cependant, on peut espérer que
des algorithmes adaptatifs, qui apprennent (de manière approchée) la distribution de probabilité
d’importance optimale ν∗ , produiront des estimateurs dont la variance approchera la variance
minimale.
1 ∑
N
dµ
⟨µM
N , ϕ⟩ = ϕ(xi ) M (xi ) ,
N dν
i=1
où les variables aléatoires (x1 , · · · , xN ) sont i.i.d. de distribution de probabilité commune ν M ,
alors la variance de cet estimateur est
1 dµ 1 dµ
E|⟨µM
N − µ, ϕ⟩| =
2
var(ϕ M , ν M ) = [var(ϕ , ν∗ ) + ⟨µ, |ϕ|2 ⟩ χ2 (ν∗ , ν M )] ,
N dν N dν∗
et si la fonction ϕ garde un signe constant, alors
1 2
E|⟨µM
N − µ, ϕ⟩| =
2
χ (ν∗ , ν M ) ⟨µ, |ϕ|2 ⟩ ,
N
dµ
compte tenu que la variance minimale var(ϕ , ν∗ ) est nulle dans ce cas, c’est–à–dire que
dν∗
l’erreur d’échantillonnage Monte Carlo et l’erreur d’approximation due à l’apprentissage de la
distribution de probabilité d’importance optimale ν∗ se multiplient au lieu de s’additionner !
Il peut arriver que la distribution de probabilité d’importance soit seulement connue à une
constante multiplicative près, par exemple dans la situation suivante. Si la distribution de pro-
babilité µ est de la forme
gη ⟨η, g ϕ⟩
µ=g·η = c’est–à–dire ⟨µ, ϕ⟩ = ,
⟨η, g⟩ ⟨η, g⟩
9.2. SIMULATION SELON UNE DISTRIBUTION DE GIBBS–BOLTZMANN 107
ou de manière équivalente
⟨γ, ϕ⟩
⟨γ, ϕ⟩ = ⟨η, g ϕ⟩ = E[g(Ξ) ϕ(Ξ)] et ⟨µ, ϕ⟩ = , (9.3)
⟨γ, 1⟩
où la variable aléatoire Ξ a pour distribution de probabilité η, pour toute fonction mesurable
bornée ϕ, et s’il est facile
T = inf{k ≥ 1 : g(Ξk ) ≥ M Uk } ,
qui représente le nombre d’itérations de l’algorithme d’acceptation / rejet nécessaires pour pro-
duire une variable aléatoire de distribution de probabilité µ = g · η.
• de distribution de probabilité µ = g · η,
⟨η, g⟩
• et de loi géométrique de paramètre ,
M
respectivement.
∏
n−1
E[ϕ(X) 1(T = n) ] = E[ϕ(Ξn ) 1(g(Ξ ) ≥ M U ) 1(g(Ξ ) < M U ) ]
n n k k
k=1
∏
n−1
= E[ϕ(Ξn ) 1(g(Ξ ) ≥ M U ) ] P[g(Ξk ) < M Uk ]
n n
k=1
total d’itérations, c’est–à–dire le nombre total de variables aléaoires simulées, qui est une me-
sure quantitative du temps de calcul de l’algorithme d’acceptation / rejet pour générer un
N –échantillon. D’après la loi des grands nombres
1 ∑ i
N
M
T −→ E[T ] = ,
N ⟨η, g⟩
i=1
∑N
M
( T i )1/2 ⟨S N (µ) − µ, ϕ⟩ =⇒ N(0, var(ϕ, µ)) ,
⟨η, g⟩
i=1
en distribution quand N ↑ ∞.
Remarque 9.7 La variance asymptotique est d’autant plus petite que le rapport
sup g(x)
M M x∈E
= ,
⟨η, g⟩ sup g(x) ⟨η, g⟩
x∈E
est petit (proche de 1). Le premier facteur est d’autant plus petit (proche de 1) que la borne M est
proche du supremum, c’est donc une caractéristique de la méthode, tandis que le second facteur
est d’autant plus petit (proche de 1) que le recouvrement entre la distribution de probabilité η
et la fonction g est grand, c’est–à–dire que la distribution de probabilité η est concentrée autour
des points où la fonction g prend ses plus grandes valeurs, c’est donc une caractéristique du
modèle lui–même.
1 ∑
N
⟨γ, ϕ⟩ = ⟨η, g ϕ⟩ ≈ ⟨S (η), g ϕ⟩ =
N
g(ξ i ) ϕ(ξ i ) ,
N
i=1
et
∑
N
g(ξ i ) ϕ(ξ i )
1 ∑
N
i=1
⟨µ, ϕ⟩ = ⟨g · η, ϕ⟩ ≈ ⟨g · S N (η), ϕ⟩ = pourvu que g(ξ i ) > 0 ,
∑
N N
i=1
g(ξ i )
i=1
1 ∑
N
γ ≈ γ N = g S N (η) = g(ξ i ) δ i
N ξ
i=1
110 CHAPITRE 9. MÉTHODES DE MONTE CARLO
et
∑
N
g(ξ i ) ∑
N
1 ∑
N
µ ≈ µN = g · S N (η) = δ i = wi δ i pourvu que g(ξ i ) > 0 ,
i=1
∑
N ξ
i=1
ξ N
i=1
g(ξ j )
j=1
g(ξ i )
wi = pour tout i = 1, · · · , N .
∑
N
g(ξ j )
j=1
Théorème 9.8 La variable aléatoire ⟨γ N , ϕ⟩ est un estimateur non–biaisé de ⟨γ, ϕ⟩, et les mo-
ments de l’erreur d’estimation vérifient
⟨γ N − γ, ϕ⟩ ⟨S N (η) − η, g ϕ⟩
= ,
⟨γ, 1⟩ ⟨η, g⟩
Théorème 9.10 La variable aléatoire ⟨µN , ϕ⟩ est un estimateur biaisé de ⟨µ, ϕ⟩, avec
E[ ⟨µN , ϕ⟩ ] = ⟨µ, ϕ⟩ + O(1/N ) , (9.10)
et les moments de l’erreur d’estimation vérifient
1 ⟨η, g 2 |ϕ − ⟨µ, ϕ⟩ |2 ⟩ 1/2
{ E| ⟨µN − µ, ϕ⟩ |2 }1/2 = √ ( ) + O(1/N ) , (9.11)
N ⟨η, g⟩2
et pour tout réel p ≥ 2
cp ⟨η, g p |ϕ − ⟨µ, ϕ⟩ |p ⟩ 1/p
{ E| ⟨µN − µ, ϕ⟩ |p }1/p ≤ √ ( ) + O(1/N ) ,
N ⟨η, g⟩p
pour toute fonction mesurable bornée ϕ.
⟨γ N − γ, ϕ⟩ ⟨γ N − γ, 1⟩ ⟨γ N − γ, 1⟩
= − ⟨µ, ϕ⟩ − ⟨µN − µ, ϕ⟩
⟨γ, 1⟩ ⟨γ, 1⟩ ⟨γ, 1⟩
⟨γ N − γ, ϕ − ⟨µ, ϕ⟩⟩ ⟨γ N , 1⟩
= − ⟨µN − µ, ϕ⟩ ( − 1) ,
⟨γ, 1⟩ ⟨γ, 1⟩
112 CHAPITRE 9. MÉTHODES DE MONTE CARLO
⟨γ N − γ, ϕ − ⟨µ, ϕ⟩⟩
⟨µN − µ, ϕ⟩ =
⟨γ, 1⟩
⟨γ N − γ, ϕ − ⟨µ, ϕ⟩⟩ ⟨γ N , 1⟩ ⟨γ N , 1⟩
−[ − ⟨µN − µ, ϕ⟩ ( − 1) ] ( − 1)
⟨γ, 1⟩ ⟨γ, 1⟩ ⟨γ, 1⟩
⟨γ N , 1⟩
+ ⟨µN − µ, ϕ⟩ ( − 1)2 .
⟨γ, 1⟩
⟨γ N − γ, ϕ − ⟨µ, ϕ⟩⟩
E =0,
⟨γ, 1⟩
de sorte que
⟨γ N − γ, ϕ − ⟨µ, ϕ⟩⟩ ⟨γ N , 1⟩ ⟨γ N , 1⟩
= − E[ ( − 1)] + E[⟨µN − µ, ϕ⟩ ( − 1)2 ] ,
⟨γ, 1⟩ ⟨γ, 1⟩ ⟨γ, 1⟩
⟨γ N , 1⟩
+ osc(ϕ) E| − 1|2 ,
⟨γ, 1⟩
où les deux termes dans la majoration sont d’ordre 1/N d’après (9.6) et (9.8).
9.2. SIMULATION SELON UNE DISTRIBUTION DE GIBBS–BOLTZMANN 113
Pour l’étude du moment d’ordre 2, en utilisant l’identité (9.8) et l’inégalité triangulaire, puis
l’inégalité de Hölder et la majoration grossière (9.12), on obtient
1 ⟨η, g 2 |ϕ − ⟨µ, ϕ⟩ |2 ⟩ 1/2
| {E|⟨µN − µ, ϕ⟩|2 }1/2 − √ ( ) |
N ⟨η, g⟩2
⟨γ N − γ, ϕ − ⟨µ, ϕ⟩⟩ ⟨γ N , 1⟩ ⟨γ N , 1⟩
≤ {E| ( − 1) |2 }1/2 + {E|⟨µN − µ, ϕ⟩ ( − 1)2 |2 }1/2
⟨γ, 1⟩ ⟨γ, 1⟩ ⟨γ, 1⟩
⟨γ N − γ, ϕ − ⟨µ, ϕ⟩⟩ ⟨γ N , 1⟩
+ {E| ( − 1) |p }1/p
⟨γ, 1⟩ ⟨γ, 1⟩
⟨γ N , 1⟩
+ {E|⟨µN − µ, ϕ⟩ ( − 1)2 |p }1/p
⟨γ, 1⟩
cp ⟨η, g p |ϕ − ⟨µ, ϕ⟩|p ⟩ 1/p
≤ √ ( )
N ⟨η, g⟩p
⟨γ N , 1⟩
+ osc(ϕ) {E| − 1|2p }1/p ,
⟨γ, 1⟩
où les deux derniers termes dans la majoration sont d’ordre 1/N d’après (9.7) et (9.9). 2
Théorème 9.11
√ ⟨γ N , 1⟩ √
N[ − 1] =⇒ N(0, V ) et N ⟨µN − µ, ϕ⟩ =⇒ N(0, v(ϕ)) ,
⟨γ, 1⟩
en distribution quand N ↑ ∞, pour toute fonction mesurable bornée ϕ, avec l’expression suivante
pour la variance asymptotique
⟨η, g 2 ⟩ ⟨η, g 2 |ϕ − ⟨µ, ϕ⟩|2 ⟩
V = −1 et v(ϕ) = ,
⟨η, g⟩2 ⟨η, g⟩2
114 CHAPITRE 9. MÉTHODES DE MONTE CARLO
respectivement.
Remarque 9.12 On vérifie que la variance asymptotique V coı̈ncide avec la variance non–
asymptotique donnée en (9.6), et que la variance asymptotique v(ϕ) coı̈ncide avec le terme
dominant de l’erreur quadratique moyenne non–asymptotique donnée en (9.11) ou de manière
équivalente avec le terme dominant de la variance non–asymptotique, compte tenu que le biais
donné en (9.10) est asymptotiquement négligeable.
Preuve. Il résulte du théorème central limite (dans sa version classique, pour des variables
indépendantes identiquement distribuées), que
√ ⟨γ N − γ, ϕ⟩ √ ⟨S N (η) − η, g ϕ⟩ var(g ϕ, η)
N = N =⇒ N(0, ),
⟨γ, 1⟩ ⟨η, g⟩ ⟨η, g⟩2
en distribution quand N ↑ ∞, et en particulier pour ϕ ≡ 1
√ ⟨γ N , 1⟩ √ ⟨γ N − γ, 1⟩ var(g, η)
N[ − 1] = N =⇒ N(0, ),
⟨γ, 1⟩ ⟨γ, 1⟩ ⟨η, g⟩2
en distribution quand N ↑ ∞. On remarque aussi que
1 ∑
N
⟨γ N , 1⟩ = g(ξ i ) −→ ⟨η, g⟩ = ⟨γ, 1⟩ ,
N
i=1
de sorte que
var(g (ϕ − ⟨µ, ϕ⟩), η) = ⟨η, g 2 |ϕ − ⟨µ, ϕ⟩|2 ⟩ . 2
Remarque 9.13 Le principe de l’échantillonnage pondéré est illustré sur la Figure 9.2. On
constate en particulier que l’approximation sera d’autant meilleure que la distribution d’impor-
tance η (densité a priori) et la fonction d’importance g (vraisemblance) se recouvrent mutuelle-
ment, de telle sorte que la fonction d’importance prend des valeurs significatives sur l’échantillon
généré. La décomposition d’importance sera au contraire mal–posée si les valeurs significatives
de la fonction d’importance sont obtenues dans les queues de la distribution d’importance, au-
quel cas la fonction d’importance prend des valeurs négligeables sur l’échantillon généré. Un
9.2. SIMULATION SELON UNE DISTRIBUTION DE GIBBS–BOLTZMANN 115
∫ sup g(x)
g(x′ )
η(dx′ )
x∈E
ou par le rapport inverse r= .
E sup g(x) ⟨η, g⟩
x∈E
116 CHAPITRE 9. MÉTHODES DE MONTE CARLO
1.2
0.8
0.6
0.4
0.2
0
−5 −4 −3 −2 −1 0 1 2 3 4 5
1.2
0.8
0.6
0.4
0.2
0
−5 −4 −3 −2 −1 0 1 2 3 4 5
Figure 9.1 – Densité a priori, échantillon (en haut) et histogramme associé à l’échantillon (en
bas)
9.2. SIMULATION SELON UNE DISTRIBUTION DE GIBBS–BOLTZMANN 117
0.8
0.6
0.4
0.2
0
−5 −4 −3 −2 −1 0 1 2 3 4 5
0.8
0.6
0.4
0.2
0
−5 −4 −3 −2 −1 0 1 2 3 4 5
∑
M ∑
M
η= wi m i avec wi = 1 ,
i=1 i=1
alors il est facile, en principe, de simuler une variable aléatoire selon la distribution de proba-
bilité η. Il suffit en effet de simuler d’abord une variable aléatoire I à valeurs dans l’ensemble
fini {1, · · · , M } et distribuée selon les poids (w1 , · · · , wM ), c’est–à–dire P[I = i] = wi pour
tout i = 1, · · · , M , puis de générer une variable aléatoire distribuée selon mI . La probabilité
de sélectionner une composante du mélange sera d’autant plus grande que le poids de cette
composante est grand. La question qui reste, et qui sera traitée à la Section 9.4, est donc de
savoir simuler une variable aléatoire à valeurs dans l’ensemble fini {1, · · · , M }.
1 ∑
N
S N (η) = δξi . (9.13)
N
i=1
9.3. ÉCHANTILLONNAGE ET APPROXIMATION D’UN MÉLANGE FINI 119
Les poids sont ici exploités pour sélectionner (avec remise) les composantes du mélange les mieux
loties, avec l’effet attendu que les composantes de plus forts poids seront sélectionnées plusieurs
fois, tandis que les composantes de moins forts poids pourront même être éliminées et ne plus
être représentées du tout dans l’approximation. Le nombre de fois que la i–ème composante
du mélange est sélectionée, ou de manière équivalente son nombre Ni de représentants dans
l’approximation, sera d’autant grand que le poids wi de cette composante est grand, et on peut
montrer que le vecteur aléatoire (N1 , · · · , NM ) suit une loi multinomiale. La question qui reste,
et qui sera traitée à la Section 9.4, est donc de savoir simuler un N –échantillon à valeurs dans
l’ensemble fini {1, · · · , M }, plus efficacement qu’en répétant N fois la simulation d’une seule
variable aléatoire.
Intuitivement, si tous les poids sont égaux à (ou proches de) 1/M , c’est–à–dire si la répartition
des poids de mélange est proche de l’équidistribution, alors il est inutile voire même contre–
productif de sélectionner les composantes du mélange, avec le risque de favoriser certaines com-
posantes au détriment des autres composantes, alors qu’en principe toutes les composantes ont
la même importance.
Théorème 9.14 La variable aléatoire ⟨S N (η), ϕ⟩ est un estimateur non–biaisé de ⟨η, ϕ⟩, et les
moments de l’erreur d’estimation vérifient
1
E| ⟨S N (η) − η, ϕ⟩ |2 = var(ϕ, η)
N
et pour tout réel p ≥ 2, il existe une constante positive cp > 0 telle que
cp
{ E| ⟨S N (η) − η, ϕ⟩|p }1/p ≤ √ ⟨η, |ϕ − ⟨η, ϕ⟩|p ⟩1/p ,
N
pour toute fonction mesurable bornée ϕ.
∑
M ∑
M ∑
M
= wi var(ϕ, mi ) + [ wi |⟨mi , ϕ⟩| − |
2
wi ⟨mi , ϕ⟩ |2 ] ,
i=1
| i=1 {z i=1 }
WM
où le terme WM représente la variance des moyennes intra–composantes affectées du poids de
chaque composante.
120 CHAPITRE 9. MÉTHODES DE MONTE CARLO
▶ Conservation des poids A l’opposé, on peut décider de conserver les poids et de simuler
un représentant pour chaque composante du mélange (ce qui impose que N est nécessairement
égal au nombre M de composantes du mélange initial), et poser
∑
M
ηM = wi δ ξ , (9.14)
i
i=1
Théorème 9.17 La variable aléatoire ⟨ηM , ϕ⟩ est un estimateur non–biaisé de ⟨η, ϕ⟩, et les
moments de l’erreur d’estimation vérifient
∑
M ∑
M
E| ⟨ηM − η, ϕ⟩ | = (
2
wi2 ) [ wi□ var(ϕ, mi ) ] ,
i=1 i=1
et pour tout réel p ≥ 2, il existe une constante positive cp > 0 telle que
∑M ∑
M
{ E| ⟨ηM − η, ϕ⟩|p }1/p ≤ cp ( wi2 )1/2 { wi□ ⟨mi , |ϕ − ⟨mi , ϕ⟩|p ⟩ }1/p ,
i=1 i=1
∑
M
wi□ = wi2 / [ wj2 ] pour tout i = 1, · · · , M .
j=1
∑
M
E| ⟨ηM − η, ϕ⟩ |2 = E| wi [ϕ(ξ i ) − ⟨mi , ϕ⟩ ] |2
i=1
∑
M
= wi2 E|ϕ(ξ i ) − ⟨mi , ϕ⟩|2
i=1
∑M ∑
M
= ( wi2 ) [ wi□ var(ϕ, mi ) ] ,
i=1 i=1
pour toute fonction mesurable bornée ϕ. Plus généralement, pour tout réel p ≥ 2, il résulte de
l’inégalité de Marcinkiewicz–Zygmund (B.2) que
∑
M
E| ⟨ηM − η, ϕ⟩ |p = E| wi [ϕ(ξ i ) − ⟨mi , ϕ⟩] |p
i=1
∑M ∑
M
≤ Bp ( 2 p/2
wi ) [ wi□ E| ϕ(ξ i ) − ⟨mi , ϕ⟩|p ]
i=1 i=1
∑M ∑
M
= Bp ( wi2 )p/2 [ wi□ ⟨mi , |ϕ − ⟨mi , ϕ⟩|p ⟩ ] ,
i=1 i=1
1 ∑ ∑ 1 ∑
M M M
Vred = [ wi ⟨mi , |ϕ|2 ⟩ − | wi ⟨mi , ϕ⟩|2 ] ≥ wi var(ϕ, mi ) ,
N N
i=1 i=1 i=1
∑
M
Vnored = wi2 var(ϕ, mi ) ,
i=1
de l’estimateur (9.14). A l’équidistribution, c’est–à–dire si tous les poids sont égaux entre eux
(et égaux à 1/M ), alors
1 ∑
M
Vred ≥ 2 var(ϕ, mi ) = Vnored ,
M
i=1
122 CHAPITRE 9. MÉTHODES DE MONTE CARLO
ce qui confirme l’intuition que redistribuer est contre–productif dans ce cas extrême. A l’inverse,
si la distribution des poids est complètement dégénérée, c’est–à–dire si tous les poids sont nuls
sauf le poids wa = 1 pour la composante a du mélange, alors
1 1
Vred = [ ⟨ma , |ϕ|2 ⟩ − |⟨ma , ϕ⟩|2 ] = var(ϕ, ma ) ≤ var(ϕ, ma ) = Vnored ,
N N
ce qui confirme l’intuition que redistribuer est certainement pertinent dans cet autre cas extrême.
N w i = N i + qi avec 0 ≤ qi < 1 .
∑
M ∑
M ∑
M ∑
M
N= N wi = (Ni + qi ) = Ni + N0 avec N0 = qi ,
i=1 i=1 i=1 i=1
et
∑
M ∑
M
Ni ∑
M
qi ∑
M
Ni N0
η= wi m i = mi + mi = mi + m0 ,
N N N N
i=1 i=1 i=1 i=1
avec
∑
M
qi
m0 = mi ,
N0
i=1
on déduit que (N − N0 ) représentants ont déjà été affectés à l’issue de cette première passe, et
il reste donc N0 représentants à affecter de manière à approcher la distribution de probabilité
résiduelle convenablement renormalisée m0 . L’approximation proposée consiste
1 ∑∑ 1 ∑
M Ni N0
ηN = δ i,j + δ 0,j , (9.15)
N ξ N ξ
i=1 j=1 j=1
9.3. ÉCHANTILLONNAGE ET APPROXIMATION D’UN MÉLANGE FINI 123
c’est–à–dire que
1 ∑∑ 1 ∑
M Ni N0
⟨ηN , ϕ⟩ = i,j
ϕ(ξ ) + ϕ(ξ 0,j ) ,
N N
i=1 j=1 j=1
et par différence
1 ∑∑ 1 ∑
M Ni N0
⟨ηN − η, ϕ⟩ = [ϕ(ξ ) − ⟨mi , ϕ⟩] +
i,j
[ϕ(ξ 0,j ) − ⟨m0 , ϕ⟩] ,
N N
i=1 j=1 j=1
Théorème 9.20 La variable aléatoire ⟨ηN , ϕ⟩ est un estimateur non–biaisé de ⟨η, ϕ⟩, et les
moments de l’erreur d’estimation vérifient
1 ∑ Ni
M
N0
E| ⟨ηN − η, ϕ⟩ |2 = [ var(ϕ, mi ) + var(ϕ, m0 ) ] ,
N N N
i=1
et pour tout réel p ≥ 2, il existe une constante positive cp > 0 telle que
c p ∑ Ni
M
N0
{ E| ⟨ηN − η, ϕ⟩|p }1/p ≤ √ [ ⟨mi , |ϕ − ⟨mi , ϕ⟩|p ⟩ + ⟨m0 , |ϕ − ⟨m0 , ϕ⟩|p ⟩ ]1/p ,
N i=1 N N
∑
M
qi ∑
M
qi
var(ϕ, m0 ) = ⟨mi , |ϕ|2 ⟩ − | ⟨mi , ϕ⟩ |2 ,
N0 N0
i=1 i=1
124 CHAPITRE 9. MÉTHODES DE MONTE CARLO
∑
M
N0 ∑ qi
M
= wi var(ϕ, mi ) − [⟨mi , |ϕ|2 ⟩ − |⟨mi , ϕ⟩ |2 ]
N N0
i=1 i=1
N0 ∑ qi ∑
M M
qi
+ [ ⟨mi , |ϕ|2 ⟩ − | ⟨mi , ϕ⟩ |2 ]
N N0 N0
i=1 i=1
∑
M
N0 ∑ qi
M ∑
M
qi
= wi var(ϕ, mi ) + [ |⟨mi , ϕ⟩|2 − | ⟨mi , ϕ⟩ |2 ] ,
N N0 N0
i=1
| i=1 {z i=1 }
WM
où le terme WM représente la variance des moyennes intra–composantes affectées du poids rési-
duel de chaque composante.
∑
M
Ni ∑
M
Ni
= ⟨mi , |ϕ| ⟩ −
2
|⟨mi , ϕ⟩ |2
N N
i=1 i=1
N 0 ∑ qi N 0 ∑ qi
M M
+ ⟨mi , |ϕ|2 ⟩ − | ⟨mi , ϕ⟩ |2
N N0 N N0
i=1 i=1
∑
M ∑
M
Ni N 0 ∑ qi
M
= wi ⟨mi , |ϕ| ⟩ − [
2
|⟨mi , ϕ⟩ | + 2
| ⟨mi , ϕ⟩ |2 ] ,
N N N0
i=1 i=1 i=1
tandis que
∑
M ∑
M
var(ϕ, η) = wi ⟨mi , |ϕ| ⟩ − |
2
wi ⟨mi , ϕ⟩ |2
i=1 i=1
∑
M ∑
M
Ni N 0 ∑ qi
M
= wi ⟨mi , |ϕ|2 ⟩ − | ⟨mi , ϕ⟩ + ⟨mi , ϕ⟩ |2 .
N N N0
i=1 i=1 i=1
1 ∑∑ 1 ∑
M Ni N0
E| ⟨ηN − η, ϕ⟩ |2 = E| [ϕ(ξ i,j ) − ⟨mi , ϕ⟩] + [ϕ(ξ 0,j ) − ⟨m0 , ϕ⟩] |2
N N
i=1 j=1 j=1
1 ∑∑ 1 ∑
M Ni N0
= E|ϕ(ξ i,j
) − ⟨m i , ϕ⟩|2
+ E|ϕ(ξ 0,j ) − ⟨m0 , ϕ⟩|2
N2 N2
i=1 j=1 j=1
1 ∑ Ni
M
N0
= [ var(ϕ, mi ) + var(ϕ, m0 ) ] ,
N N N
i=1
pour toute fonction mesurable bornée ϕ. Plus généralement, pour tout réel p ≥ 2, il résulte de
l’inégalité Marcinkiewicz–Zygmund (B.1) que
1 ∑∑ 1 ∑
M Ni N0
E| ⟨ηN − η, ϕ⟩ | = E|
p
[ϕ(ξ ) − ⟨mi , ϕ⟩] +
i,j
[ϕ(ξ 0,j ) − ⟨m0 , ϕ⟩] |p
N N
i=1 j=1 j=1
Bp 1 ∑ ∑ ∑
M Ni N0
≤ p/2 [ E| ϕ(ξ i,j ) − ⟨mi , ϕ⟩|p + E| ϕ(ξ 0,j ) − ⟨m0 , ϕ⟩|p ]
N N
i=1 j=1 j=1
Bp ∑ Ni
M
N0
= p/2
[ ⟨mi , |ϕ − ⟨mi , ϕ⟩|p ⟩ + ⟨m0 , |ϕ − ⟨m0 , ϕ⟩|p ⟩ ] ,
N N N
i=1
il n’est véritablement intéressant de sélectionner les composantes du mélange que si les poids
(w1 , · · · , wM ) sont très déséquilibrés. Plusieurs critères ont été proposés pour mesurer l’écart
à l’équidistribution, et pour décider de conserver les poids ou de les utiliser pour échantilloner
selon l’une ou l’autre des implémentations proposées ci–dessus, par exemple la distance du
126 CHAPITRE 9. MÉTHODES DE MONTE CARLO
∑
M
pi ∑
M
pi
2
χ (p, q) = qi ( − 1)2 et K(p, q) = pi log ,
qi qi
i=1 i=1
respectivement.
1 ∑ 1 ∑ ∑
M M M
M
0≤ (M wi − 1) =
2
(M wi ) − 1 = M
2
wi2 − 1 = −1 ,
M M Meff
i=1 i=1 i=1
∑
M
1 ≤ Meff = 1 / [ wi2 ] ≤ M ,
i=1
∑
M
Ent = − wi log wi ≤ log M ,
i=1
Remarque 9.24 Les résultats obtenus au Théorème 9.17 pour l’estimateur (9.14) peuvent être
ré–interprétés en terme de la taille effective de l’échantillon. En effet
1 ∑ □
M
E| ⟨ηM − η, ϕ⟩ |2 = wi var(ϕ, mi ) ,
Meff
i=1
9.4. ÉCHANTILLONNAGE SELON UNE DISTRIBUTION À SUPPORT FINI 127
et pour tout réel p ≥ 2, il existe une constante positive cp > 0 telle que
cp ∑ M
{ E| ⟨ηM − η, ϕ⟩|p }1/p ≤ √ { w□ ⟨mi , |ϕ − ⟨mi , ϕ⟩|p ⟩ }1/p ,
Meff i=1 i
Remarque 9.25 On peut aussi introduire l’approximation particulaire adaptative définie par
M
∑
avec ξi ∼ mi pour tout i = 1 · · · M
wi δ ξ
i
i=1
si Meff > cred M ,
ηM =
1 ∑
M
avec ξi ∼ η pour tout i = 1 · · · M
δξ
M i
i=1
si Meff ≤ cred M ,
avec l’expression suivante pour la variance de l’erreur d’estimation
1 ∑ □
M
wi var(ϕ, mi ) si Meff > cred M ,
Meff
E| ⟨ηM − η, ϕ⟩ |2 = i=1
1 var(ϕ, η) si Meff ≤ cred M .
M
La question qui reste est donc de savoir simuler une variable aléatoire I à valeurs dans l’en-
semble fini {1, · · · , M }, ou bien de savoir simuler un N –échantillon à valeurs dans l’ensemble
fini {1, · · · , M }, et on commence par introduire un découpage de l’intervalle [0, 1] en M seg-
ments adjacents de longueurs respectives w1 , · · · , wM . Ces M segments sont délimités par les
probabilités cumulées
s0 = 0 et si = w1 + · · · + wi ,
pour tout i = 1, · · · , M et on vérifie que sM = 1.
La méthode la plus directe est la méthode d’inversion, qui consiste à générer une variable
aléatoire U uniforme sur [0, 1] : si U appartient au j–ème segment, i.e. si
sj−1 < U ≤ sj ,
alors on pose I = j. Une recherche binaire en O(log2 M ) opérations permet d’obtenir ce résultat,
et il suffit donc de N O(log2 M ) opérations pour générer un N –échantillon à valeurs dans l’en-
semble fini {1, · · · , M } et distribué selon les poids (w1 , · · · , wM ).
Au lieu de répéter N fois l’opération de
128 CHAPITRE 9. MÉTHODES DE MONTE CARLO
0 ≤ U(1) ≤ · · · ≤ U(N ) ≤ 1 et 0 = s0 ≤ s1 ≤ · · · ≤ sM = 1 .
? ? ? ? ? ? ?
Avec cette méthode, il suffit donc de O(N log2 N )+O(M +N ) opérations pour générer un N –
échantillon à valeurs dans l’ensemble fini {1, · · · , M } et distribué selon les poids (w1 , · · · , wM ).
Une méthode plus efficace, qui évite l’étape préalable de ré–ordonner les variables aléatoires
uniformes, consiste à générer directement une N –statistique d’ordre uniforme, c’est–à–dire un
vecteur aléatoire (V1 , · · · , VN ) distribué comme le vecteur aléatoire (U(1) , · · · , U(N ) ) obtenu en
ré–ordonnant un N –échantillon (U1 , · · · , UN ) de variables aléatoires uniformes sur [0, 1]. L’un ou
l’autre des deux résultats suivants permet d’effectuer cette tâche en O(N ) opérations, et il suffit
donc de O(N ) + O(M + N ) opérations pour générer un N –échantillon à valeurs dans l’ensemble
fini {1, · · · , M } et distribué selon les poids (w1 , · · · , wM ).
1/i
Vi = Vi+1 Ui pour tout i = N − 1, · · · , 1.
Le vecteur aléatoire (V1 , · · · , VN ) est distribué comme le vecteur aléatoire (U(1) , · · · , U(N ) ) obtenu
en ré–ordonnant (U1 , · · · , UN ).
9.4. ÉCHANTILLONNAGE SELON UNE DISTRIBUTION À SUPPORT FINI 129
V1 ≤ · · · ≤ Vi ≤ Vi+1 ≤ · · · ≤ VN ,
pour toute fonction mesurable bornée ϕ définie sur l’intervalle [0, 1], de sorte que
∫ 1
= ϕ(Vi+1 x1/i ) P[Ui ∈ dx | Vi+1 , · · · , VN ]
0
∫ 1
= ϕ(Vi+1 x1/i ) dx
0
∫ 1
i v i−1
= ϕ(v) 1(0 ≤ v ≤ V ) i dv ,
0 i+1 V
i+1
pour toute fonction mesurable bornée ϕ définie sur l’intervalle [0, 1], de sorte que
i v i−1
P[Vi ∈ dv | Vi+1 , · · · , VN ] = 1(0 ≤ v ≤ V ) i dv ,
i+1 V
i+1
∏
N −1
P[V1 ∈ dv1 , · · · , VN ∈ dvN ] = P[VN ∈ dvN ] P[Vi ∈ dvi | Vi+1 = vi+1 , · · · , VN = vN ]
i=1
∏
N −1
N −1 i v i−1
= 1(0 ≤ v ≤ 1) N vN dvN 1(0 ≤ v ≤ v ) ii dvi
N
i=1
i i+1 vi+1
Si
Si = E1 + · · · + Ei et Vi = pour tout i = 1, · · · , N ,
SN +1
130 CHAPITRE 9. MÉTHODES DE MONTE CARLO
E1
ou bien par récurrence : V1 = et
SN +1
Ei
Vi = Vi−1 + pour tout i = 2, · · · , N .
SN +1
Le vecteur aléatoire (V1 , · · · , VN ) est distribué comme le vecteur aléatoire (U(1) , · · · , U(N ) ) obtenu
en ré–ordonnant un N –échantillon (U1 , · · · , UN ) de variables aléatoires uniformes sur [0, 1].
V1 ≤ · · · ≤ Vi ≤ Vi+1 ≤ · · · ≤ VN ,
Compte tenu que (S1 , · · · , Si ) contient exactement la même information que (E1 , · · · , Ei ), et
que la variable aléatoire Ei+1 est indépendante de (E1 , · · · , Ei ), on a
pour toute fonction ϕ mesurable bornée définie sur l’intervalle [0, ∞), de sorte que
∏
N
P[S1 ∈ ds1 , · · · , SN +1 ∈ dsN +1 ] = P[S1 ∈ ds1 ] P[Si+1 ∈ dsi | S1 = s1 , · · · , Si = si ]
i=1
∏
N
= 1(s ≥ 0) e−s1 ds1 1(s e−(si+1 −si ) dsi
i+1 ≥ si )
i=1
et
S1 SN
E[f (V1 , · · · , VN )] = E[f ( ,··· , )]
SN +1 SN +1
∫ ∞ ∫ ∞
s1 sN
= ··· f( ,··· , ) P[S1 ∈ ds1 , · · · , SN +1 ∈ dsN +1 ]
0 0 sN +1 sN +1
∫ ∞ ∫ ∞
s1 sN
= ··· f( ,··· , ) 1(0 ≤ s ≤ · · · ≤ s ≤ s e−sN +1 ds1 · · · dsN dsN +1
0 0 sN +1 sN +1 1 N N +1 )
∫ ∞ ∫ ∞
−sN +1
= ··· f (v1 , · · · , vN ) 1(0 ≤ v ≤ · · · ≤ v ≤ 1) sN
N +1 e dv1 · · · dvN dsN +1
0 0 1 N
∫ ∞ ∫ ∞ ∫ ∞
−sN +1
= ··· f (v1 , · · · , vN ) 1(0 ≤ v ≤ · · · ≤ v ≤ 1) dv1 · · · dvN sN
N +1 e dsN +1
0 0 1 N 0
∫ ∞ ∫ ∞
= ··· f (v1 , · · · , vN ) N ! 1(0 ≤ v ≤ · · · ≤ v ≤ 1) dv1 · · · dvN ,
0 0 1 N
pour toute fonction f mesurable bornée définie sur l’ensemble produit [0, ∞)N , compte tenu que
∫ ∞
sN e−s ds = N ! .
0
On en déduit que
Approximations particulaires
On se place dans le cadre général décrit au Chapitre 8, où différents modèles ont été considérés,
avec différents points de vue. Il ressort de la discussion que le modèle (8.22), où chaque fonction
de sélection dépend de la transition courante de la chaı̂ne de Markov, semble suffisamment
général pour inclure comme cas particuliers la plupart des modèles présentés jusqu’ici, mais il
ressort aussi de la discussion que le modèle (8.1) apparamment plus simple, où chaque fonction
de sélection dépend seulement de l’état courant (c’est–à–dire de l’état d’arrivée de la transition
courante) de la chaı̂ne de Markov, contient en fait le modèle (8.22) comme cas particulier,
pourvu qu’on change de point de vue et qu’on adopte le modèle (8.27) à valeurs transitions.
Cette remarque sera abondamment exploitée dans ce chapitre et dans les chapitres suivants, et
on considérera indifférement ces deux modèles et les différents points de vue associés.
Il s’agit ici d’approcher numériquement, par des méthodes de Monte Carlo, la distribution
non–normalisée et la distribution normalisée associée, définis
respectivement, selon le modèle utilisé. Le premier point de vue conduit aux algorithmes d’é-
chantillonnage pondéré (SIS, pour sequential importance sampling) qui sont des algorithmes de
Monte Carlo classiques (sans interaction), assez inefficaces, et le second point de vue conduit
aux algorithmes d’échantillonnage / ré–échantillonnage (SIR, pour sampling with importance
resampling) qui sont des algorithmes de Monte Carlo avec interaction, beaucoup plus efficaces.
133
134 CHAPITRE 10. APPROXIMATIONS PARTICULAIRES
et la fonction de poids
∏
n
g0:n (x0:n ) = gk (xk−1 , xk ) ,
k=0
et si on définit
∫ ∫
⟨η0:n , g0:n f ⟩ = E[f (X0:n ) g0:n (X0:n ) ] = ··· f (x0:n ) g0:n (x0:n ) η0:n (dx0:n ) ,
E E
pour toute fonction mesurable bornée f définie sur l’espace produit E × · · · × E = E n+1 , alors
on a
∏
n
⟨γn , ϕ⟩ = E[ϕ(Xn ) gk (Xk ) ] = E[ϕ ◦ π(X0:n ) g0:n (X0:n )] = ⟨η0:n , g0:n ϕ ◦ π⟩ ,
k=0
• pour tout k = 1, · · · , n, la variable aléatoire ξki est distribuée selon mik = Qk (ξk−1
i , ·),
∑
N
−1
µn ≈ µN
n = µ0:n ◦ π
N
= wni δ .
ξni
i=1
La simulation des trajectoires et le calcul des poids peuvent être décrits de la façon récursive
suivante
g0 (ξ0i )
w0i = ,
∑
N
g0 (ξ0j )
j=1
i
wk−1 i
gk (ξk−1 , ξki )
wki = .
∑
N
j j
wk−1 gk (ξk−1 , ξkj )
j=1
Exemple 10.1 Dans le cas particulier du système non–linéaire à bruits non–gaussiens décrit
par (5.1) et (5.2), simuler une variable aléatoire X selon la distribution de probabilité Qk (x, dx′ )
signifie simplement simuler une variable aléatoire W selon la distribution de probabilité pW k (dw),
et poser X = fk (x, W ), et évaluer la fonction de vraisemblance gk (x′ ) signifie simplement évaluer
qkV (Yk − hk (x′ )), d’où l’algorithme suivant
j=1
et on pose
i
wk−1 qkV (Yk − hk (ξki ))
wki = .
∑
N
j
wk−1 qkV (Yk − hk (ξkj ))
j=1
On remarque que les poids dépendent des trajectoires simulées, et sont en effet d’autant
plus élevés pour les trajectoires en adéquation avec les observations, mais en revanche ces poids
ne sont pas utilisés pour simuler les trajectoires : en poussant les choses à l’extrême, on peut
donc dire que les trajectoires sont simulées en aveugle, et l’algorithme se contente de pondérer
différemment les différentes trajectoires.
Comme ces différentes trajectoires sont en nombre fini, il est de moins en moins raisonnable,
au fur et à mesure que le temps passe, d’espérer qu’un nombre suffisant d’entre ces trajectoires
puisse être assez proche de la vraie trajectoire. Comme les poids s’accumulent au cours du temps
le long de chaque trajectoire, la situation typique est de voir une seule trajectoire recueillir un
poids beaucoup plus fort que toutes les autres, et ceci juste parce qu’au cours de son histoire
passée elle s’est trouvée plus souvent proche de la vraie trajectoire, quand bien même elle s’en
trouverait très éloignée à l’instant présent.
Ces phénomènes de dégénerescence des poids et d’importance excessive du passé sont bien
connus, et diverses solutions ont été proposées pour y remédier
• simuler les trajectoires selon un mécanisme qui prenne mieux en compte les observations,
au lieu de simuler les trajectoires en aveugle,
• multiplier les trajectoires de poids le plus fort, et éliminer les trajectoires de poids le plus
faible, en introduisant une étape de ré–échantillonage,
et il est également possible de combiner ces solutions, ce qui fait l’objet des algorithmes présentés
ci–dessous.
Il est plus facile ici de considérer d’abord le modèle (8.1) apparamment plus simple, où chaque
fonction de sélection dépend seulement de l’état courant (c’est–à–dire de l’état d’arrivée de la
transition courante) de la chaı̂ne de Markov, et de considérer ensuite comme un cas particulier le
modèle (8.22) apparamment plus général, où chaque fonction de sélection dépend de la transition
courante de la chaı̂ne de Markov, ce qui est possible pourvu qu’on change de point de vue et
qu’on adopte le modèle (8.27) à valeurs transitions.
Au lieu de simuler d’abord N trajectoires indépendantes de la chaı̂ne de Markov et d’évaluer
séparément les poids associés à chaque trajectoire simulée, le principe consiste à rechercher une
approximation
1 ∑ ∑ ∑
N N N
ηk ≈ ηkN = δ i et µk ≈ µN
k = wki δ avec wki = 1 , (10.1)
N ξk ξki
i=1 i=1 i=1
10.2. ÉCHANTILLONNAGE / RÉ–ÉCHANTILLONNAGE (SIR) 137
mutation pondération
µk−1 −−−−−−−−−−→ ηk = µk−1 Qk −−−−−−−−−−−−→ µk = gk · ηk , (10.2)
1 ∑
N
η0N = S N (η0 ) = δ i ,
N ξ0
i=1
où indépendamment pour tout i = 1, · · · , N , la variable aléatoire ξ0i est distribuée selon η0 .
Pour tout k = 0, 1, · · · , n, il est immédiat à partir de la définition (10.1) que
∑
N
gk (ξki ) ∑ N
µN
k = gk · ηkN = δ i = wki δ i ,
i=1
∑
N
j
ξk
i=1
ξk
gk (ξk )
j=1
gk (ξki )
wki = pour tout i = 1, · · · , N , (10.3)
∑
N
gk (ξkj )
j=1
et
1 ∑
N
⟨ηkN , gk ⟩ = gk (ξki ) ,
N
i=1
∑
N
µN
k−1 Qk = i
wk−1 mik où mik (dx′ ) = Qk (ξk−1
i
, dx′ ) ,
i=1
sous la forme désirée, et où indépendamment pour tout i = 1, · · · , N la variable aléatoire ξki
est distribuée selon le mélange fini µNk−1 Qk . On en déduit l’expression suivante pour l’erreur
d’approximation
1 ∑
N
⟨ηk − µk−1 Qk , ϕ⟩ =
N N
[ ϕ(ξki ) − ⟨µN
k−1 Qk , ϕ⟩ ] ,
N
i=1
pour toute fonction mesurable bornée ϕ, où Fk−1 N désigne la tribu engendrée par le système de
particules jusqu’à la (k − 1)–ème génération.
∑
N ∑
N ∑
N ∑
N
i
N= N wk−1 = (Nki + qki ) = Nki + Nk0 avec Nk0 = qki ,
i=1 i=1 i=1 i=1
et
∑
N ∑
N
Ni ∑
N
qi ∑
N
Ni Nk0 0
µN
k−1 Qk =
i
wk−1 mik = k
mik + k
mik = k
mik + mk ,
N N N N
i=1 i=1 i=1 i=1
avec
∑N
qki
m0k = i
0 mk ,
i=1
N k
10.2. ÉCHANTILLONNAGE / RÉ–ÉCHANTILLONNAGE (SIR) 139
on déduit que (N − Nk0 ) descendants ont déjà été affectés à l’issue de cette première passe, et
il reste donc Nk0 descendants à affecter de manière à approcher la distribution de probabilité
résiduelle convenablement renormalisée m0k . L’approche proposée consiste
0,Nk0
• à simuler un Nk0 –échantillon (ξk0,1 , · · · , ξk ) distribué selon le mélange fini m0k ,
i,Nki
• à simuler indépendamment pour tout i = 1, · · · , N un Nki –échantillon (ξki,1 , · · · , ξk )
distribué selon mik ,
toutes les variables aléatoires étant simulées de manière indépendantes, d’où l’approximation
i 0
∑ N Nk Nk
1 ∑∑ 1 ∑ 1 ∑
N N
µN
k−1 Qk =
i
wk−1 mik ≈ ηkN = δ i,j + δ 0,j = δ i , (10.5)
N ξk N ξk N ξk
i=1 i=1 j=1 j=1 i=1
et par différence
i 0
N Nk Nk
1 ∑∑ 1 ∑
⟨ηkN − µN
k−1 Qk , ϕ⟩ = [ϕ(ξki,j ) − ⟨mik , ϕ⟩] + [ϕ(ξk0,j ) − ⟨m0k , ϕ⟩] ,
N N
i=1 j=1 j=1
pour toute fonction mesurable bornée ϕ. Il résulte du Théorème 9.20 et de la Remarque 9.21
que
1 ∑ Nki
N
N0
E[ | ⟨ηk − µk−1 Qk , ϕ⟩ | | Fk−1 ] =
N N 2 N
[ var(ϕ, mik ) + k var(ϕ, m0k ) ] ,
N N N
i=1
et pour tout réel p ≥ 2
cp
{ E[ | ⟨ηkN − µN
k−1 Qk , ϕ⟩| | Fk−1 ] }
p N 1/p
≤ √ osc(ϕ) , (10.6)
N
pour toute fonction mesurable bornée ϕ, où Fk−1 N désigne la tribu engendrée par le système de
particules jusqu’à la (k − 1)–ème génération.
Résumé Dans toute cette classe d’algorithmes d’approximation particulaire, et quelle que soit
l’implémentation retenue pour mette en œuvre l’étape de sélection, l’évolution de la population
de particules et la mise–à–jour des poids sont ici couplées et peuvent être décrites de la façon
récursive suivante
– un individu ξbk−1
i 1 , · · · , ξN )
est sélectionné au sein de la population de particules (ξk−1 k−1
en fonction des poids (wk−1 , · · · , wk−1 ) et à l’aide de l’un des mécanismes de sélection
1 N
proposés,
– la variable aléatoire ξ i est simulée selon la distribution de probabilité Qk (ξbi , dx′ ),
k k−1
et on pose
gk (ξki )
wki = .
∑
N
gk (ξkj )
j=1
Au lieu de s’accumuler le long de chaque trajectoire comme dans le cas de l’algorithme SIS,
les poids sont ici utilisés pour redistribuer les particules, c’est–à–dire multiplier les particules de
plus fort poids et éliminer les particules de plus faible poids. Le gain escompté en ne conservant à
chaque pas de temps que les particules les plus pertinentes, est de concentrer ainsi les particules,
c’est–à–dire la puissance de calcul disponible, dans les régions d’intérêt de l’ensemble E.
Exemple 10.2 Dans le cas particulier du système non–linéaire à bruits non–gaussiens décrit
par (5.1) et (5.2), simuler une variable aléatoire X selon la distribution de probabilité Qk (x, dx′ )
signifie simplement simuler une variable aléatoire W selon la distribution de probabilité pWk (dw),
′
et poser X = fk (x, W ), et évaluer la fonction de vraisemblance gk (x ) signifie simplement évaluer
qkV (Yk − hk (x′ )), d’où l’algorithme suivant
j=1
– un individu ξbk−1
i 1 , · · · , ξN )
est sélectionné au sein de la population de particules (ξk−1 k−1
en fonction des poids (wk−11 , · · · , w N ) et à l’aide de l’un des mécanismes de sélection
k−1
proposés,
10.2. ÉCHANTILLONNAGE / RÉ–ÉCHANTILLONNAGE (SIR) 141
et on pose
qkV (Yk − hk (ξki ))
wki = .
∑
N
qkV (Yk − hk (ξkj ))
j=1
Dans le cas plus général du modèle (8.8) et pour une décomposition d’importance (8.11)
donnée, avec la représentation probabiliste (8.12) associée, ou bien pour le modèle (8.22) où la
décomposition d’importance est donnée de manière explicite dans la représentation probabiliste,
chaque fonction de sélection dépend de la transition courante de la chaı̂ne de Markov, mais il
suffit de changer de point de vue et d’adopter le modèle (8.27) à valeurs transitions, où chaque
fonction de sélection dépend seulement de l’état courant, puis de ré–exprimer dans ce cadre les
algorithmes proposés ci–dessus pour le modèle (8.1) apparamment plus simple.
Le principe consiste donc à rechercher une approximation
1 ∑ ∑
N N
ηktr ≈ ηkN,tr = δ i,1 i,2 et N,tr
k ≈ µk
µtr = wki δ , (10.7)
N (ξk , ξk ) (ξki,1 , ξki,2 )
i=1 i=1
1 ∑
N
η0N N
= S (η0 ) = δ i ,
N ξ0
i=1
où indépendamment pour tout i = 1, · · · , N , la variable aléatoire ξ0i est distribuée selon η0 .
Pour tout k = 0, 1, · · · , n, il est immédiat à partir de la définition (10.1) que
∑
N
gk (ξki,1 , ξki,2 ) ∑
N
µN,tr = gk · ηkN,tr = δ = wki δ ,
k
∑
N (ξki,1 , ξki,2 ) (ξki,1 , ξki,2 )
i=1 i=1
gk (ξkj,1 , ξkj,2 )
j=1
142 CHAPITRE 10. APPROXIMATIONS PARTICULAIRES
gk (ξki,1 , ξki,2 )
wki = pour tout i = 1, · · · , N , (10.9)
∑
N
gk (ξkj,1 , ξkj,2 )
j=1
et
1 ∑
N
⟨ηkN,tr , gk ⟩ = gk (ξki,1 , ξki,2 ) ,
N
i=1
∑
N
µN,tr
k−1 Qtr
k = i
wk−1 mik où mik (dx′1 , dx′2 ) = δ i,2 (dx′1 ) Qk (x′1 , dx′2 ) ,
ξk−1
i=1
g0 (ξ0i,2 )
w0i = .
∑
N
j,2
g0 (ξ0 )
j=1
et on pose
Finalement, l’approximation particulaire pour le modèle (8.22) est donnée sous la forme
∑
N
µk ≈ µN N,tr
k = µk ◦ π −1 = wki δ i,2 ,
ξk
i=1
−1
k ◦ π . On remarque aussitôt que toutes les étapes de cette
compte tenu de la relation µk = µtr
classe d’algorithmes d’approximation particulaire peuvent s’exprimer en terme de l’état d’arrivée
seulement des particules–transitions, ce qui donne
∑
N
µk ≈ µN = wki δ ,
k ξki
i=1
où l’évolution de la population de particules et la mise–à–jour des poids peuvent être décrites
de la façon récursive suivante
g0 (ξ0i )
w0i = .
∑
N
g0 (ξ0j )
j=1
– un individu ξbk−1
i 1 , · · · , ξN )
est sélectionné au sein de la population de particules (ξk−1 k−1
en fonction des poids (wk−1 , · · · , wk−1 ) et à l’aide de l’un des mécanismes de sélection
1 N
proposés,
– la variable aléatoire ξ i est simulée selon la distribution de probabilité Qk (ξbk−1 , dx′ ),
k
et on pose
gk (ξbk−1
i , ξ i)
k
wki = .
∑
N
gk (ξbk−1
j
, ξkj )
j=1
1 , · · · , ξN )
En résumé, les particules (ξk−1 k−1
1 , · · · , w N ) (étape de sélection),
• sont sélectionnées selon leurs poids respectifs (wk−1 k−1
144 CHAPITRE 10. APPROXIMATIONS PARTICULAIRES
Remarque 10.3 Dans le cas où les fonctions de sélection gk (x, x′ ) = gk (x′ ) ne dépendent que
de l’état d’arrivée de la transition, pour tout k = 0, 1, · · · , n, on retrouve évidemment comme
cas particulier le schéma d’approximation déjà décrit plus haut.
w0i = 1/N .
– un individu ξbk−1
i 1 , · · · , ξN )
est sélectionné au sein de la population de particules (ξk−1 k−1
en fonction des poids (wk−1 , · · · , wk−1 ) et à l’aide de l’un des mécanismes de sélection
1 N
proposés,
– la variable aléatoire ξ i est simulée selon la distribution de probabilité Qb k (ξbi , dx′ ),
k k−1
et on pose
gbk (ξbk−1
i )
wki = .
∑
N
gbk (ξbk−1
j
)
j=1
On constate néanmoins que les poids (wk1 , · · · , wkN ) qui servent à sélectionner les individus au sein
de la nouvelle population (ξk1 , · · · , ξkN ) ne dépendent en fait que de la population (ξbk−1
1 , · · · , ξbN )
k−1
et sont donc disponibles avant même que la nouvelle population ne soit générée. Il est plus efficace
dans ce cas d’effectuer la sélection plus tôt, et il suffit ici d’adopter le modèle (8.18).
Le principe consiste donc à rechercher une approximation
1 ∑ ∑ ∑
N N N
ηkopt ≈ ηkN,opt = δ i et µopt N,opt
k ≈ µk = wki δ avec wki = 1 ,
N ξk ξki
i=1 i=1 i=1
1 ∑
N
ηnopt ≈ ηnN,opt = δ i ,
N ξn
i=1
10.2. ÉCHANTILLONNAGE / RÉ–ÉCHANTILLONNAGE (SIR) 145
mutation pondération
µopt −−−−−−−−−→ ηkopt = µopt
k−1 −
opt
k−1 Qk − −−−−−−−−−−−→ µopt opt opt
k = gk · ηk , (10.10)
avec la condition initiale µopt 0 = g0opt · η0opt , où la notation · désigne le produit projectif. Ici,
la distribution initiale η0 (dx) est définie en (8.15), les noyaux markoviens Qopt
opt ′
k (x, dx ) sont
définis en (8.16) pour k = 1, · · · , n, et les fonctions de sélection gkopt (x′ ) sont définis en (8.17)
pour k = 0, · · · , (n − 1).
L’évolution de la population de particules et la mise–à–jour des poids peuvent être décrites
de la façon récursive suivante
g0opt (ξ0i )
w0i = .
∑
N
g0opt (ξ0j )
j=1
– un individu ξbk−1
i 1 , · · · , ξN )
est sélectionné au sein de la population de particules (ξk−1 k−1
en fonction des poids (wk−1 , · · · , wk−1 ) et à l’aide de l’un des mécanismes de sélection
1 N
proposés,
– la variable aléatoire ξ i est simulée selon la distribution de probabilité Qopt (ξbi , dx′ ),
k k k−1
et on pose
gkopt (ξki )
wki = .
∑
N
gkopt (ξkj )
j=1
– un individu ξbn−1
i 1 , · · · , ξN )
est sélectionné au sein de la population de particules (ξn−1 n−1
en fonction des poids (wn−1 1 , · · · , w N ) et à l’aide de l’un des mécanismes de sélection
n−1
proposés,
– la variable aléatoire ξ i est simulée selon la distribution de probabilité Qopt bi , dx′ ).
n (ξ
n n−1
En ré–organisant différemment les calculs, et en utilisant les définitions (8.15), (8.16) et (8.17),
on voit que l’évolution de la population de particules et la mise–à–jour des poids peuvent être
aussi décrites de la façon récursive suivante
146 CHAPITRE 10. APPROXIMATIONS PARTICULAIRES
– un individu ξbk−1
i 1 , · · · , ξN )
est sélectionné au sein de la population de particules (ξk−1 k−1
en fonction des poids (wk−1 1 , · · · , w N ) et à l’aide de l’un des mécanismes de sélection
k−1
proposés,
– la variable aléatoire ξki est simulée selon la distribution de probabilité Q b k (ξbi , dx′ ).
k−1
Exemple 10.4 Dans le cas particulier du système non–linéaire avec des bruits gaussiens additifs
et une fonction d’observation linéaire décrit par (7.4), on obtient l’algorithme suivant
i q(Yk − Hk fk (ξk−1
i ) − hk , Hk Σk (ξk−1
i ) Hk∗ + QVk )
wk−1 = ,
∑
N
j
q(Yk − Hk fk (ξk−1 j
) − hk , Hk Σk (ξk−1 ) Hk∗ + QVk )
j=1
∗
où par définition Σk (x) = σk (x) QW k σk (x) pour tout x ∈ E,
– on sélectionne un individu ξbk−1
i 1 , · · · , ξN )
au sein de la population de particules (ξk−1 k−1
en fonction des poids (wk−1 , · · · , wk−1 ) et à l’aide de l’un des mécanismes de sélection
1 N
proposés,
– on simule deux vecteurs aléatoires gaussiens indépendants Wki et Vki , centrés et de
matrice de covariance QW V
k et Qk respectivement,
– on pose
i
ξk|k−1 = fk (ξbk−1
i
) + Wki ,
et
ξki = ξk|k−1
i
+ Σk (ξbk−1
i
) Hk∗ [Hk Σk (ξbk−1
i
) Hk∗ + QVk ]−1 (Yk − (Hk ξk|k−1
i
+ hk + Vki )) .
10.2. ÉCHANTILLONNAGE / RÉ–ÉCHANTILLONNAGE (SIR) 147
g0 (ξ0i )
w0i = .
∑
N
g0 (ξ0j )
j=1
∑
N
i
Neff = 1 / [ (wk−1 )2 ] ,
i=1
proposés,
– la variable aléatoire ξ i est simulée selon la distribution de probabilité Qk (ξbi , dx′ ),
k k−1
et on pose
gk (ξbk−1
i , ξi )
k
wki = ,
∑
N
gk (ξbk−1
j
, ξkj )
j=1
sinon, si Neff > cred N , alors indépendamment pour tout i = 1, · · · , N , la variable aléatoire
ξki est simulée selon la distribution de probabilité Qk (ξk−1
i , dx′ ), et on pose
i
wk−1 i
gk (ξk−1 , ξki )
wki = .
∑
N
j j
wk−1 gk (ξk−1 , ξkj )
j=1
148 CHAPITRE 10. APPROXIMATIONS PARTICULAIRES
end for
[calcul des poids normalisés]
for i = 1 · · · N do
wki ∝ gk (ξki )
end for
end loop
150 CHAPITRE 10. APPROXIMATIONS PARTICULAIRES
end for
[calcul des poids normalisés]
for i = 1 · · · N do
wki ∝ gk (ξbk−1
i , ξi )
k
end for
end loop
10.2. ÉCHANTILLONNAGE / RÉ–ÉCHANTILLONNAGE (SIR) 151
end for
[sélection]
for i = 1 · · · N (indépendemment) do
choisir ξbk−1i au sein de la population (ξk−11 , · · · , ξ N ) en fonction des poids respectifs
k−1
(wk−11 , · · · , wN )
k−1
end for
[propagation]
for i = 1 · · · N (indépendemment) do
ξki ∼ Qb k (ξbi , dx′ )
k−1
end for
end loop
152 CHAPITRE 10. APPROXIMATIONS PARTICULAIRES
Chapitre 11
Estimation d’erreur
de sorte que
γkN γ0N
= gk · ηkN = µN
k et = g0 · η0N = µN
0 .
⟨γkN , 1⟩ ⟨γ0N , 1⟩
Pour k = 0, on a par différence
γ0N − γ0 = g0 (η0N − η0 ) ,
de sorte que
⟨γ0N − γ0 , ϕ⟩ = ⟨η0N − η0 , g0 ϕ⟩ , (11.2)
pour toute fonction ϕ mesurable bornée. Pour tout k = 1, · · · , n, on a par différence
153
154 CHAPITRE 11. ESTIMATION D’ERREUR
de sorte que
⟨γkN − γk , ϕ⟩ = ⟨γk−1
N
− γk−1 , Qk (gk ϕ)⟩ + ⟨ηkN − µN
k−1 Qk , gk ϕ⟩ ⟨γk−1 , 1⟩ ,
N
(11.3)
pour toute fonction ϕ mesurable bornée. On constate que l’erreur d’approximation au rang k
évaluée pour la fonction ϕ, peut s’exprimer à l’aide
La décomposition (11.3) est à la base des démonstrations par récurrence de la convergence dans
Lp et de la normalité asymptotique des approximations particulaires.
Au vu de la relation de récurrence (8.3), l’hypothèse minimale ⟨γn , 1⟩ > 0 est équivalente
à supposer que ⟨ηk , gk ⟩ > 0 pour tout k = 0, 1, · · · , n. Cette condition est trivialement vérifiée
si les fonctions gk sont strictement positives et ne s’annulent donc en aucun point, pour tout
k = 0, 1, · · · , n. En revanche, si la fonction gk peut s’annuler en certains points, et même si
⟨ηk , gk ⟩ > 0, il peut quand même arriver que gk (ξki ) = 0 pour tout i = 1, · · · , N , c’est–à–dire
que toutes les particules sont affectées d’un poids nul, auquel cas ⟨ηkN , gk ⟩ = 0 et ⟨γkN , 1⟩ = 0,
de sorte que l’approximation µN k n’est pas définie. Soit τ
N le temps d’extinction du système de
τ N = inf{k ≥ 1 : ⟨γkN , 1⟩ = 0}
1 ∑
N
= inf{k ≥ 1 : ⟨ηkN , gk ⟩ = gk (ξki ) = 0}
N
i=1
définie. Pour tout k = 1, · · · , n, sur l’ensemble {τ N > k − 1}, les distributions γk−1N et µNk−1
sont bien définies, donc les distributions ηk = S (µk−1 Qk ) et γk = gk ηk ⟨γk−1 , 1⟩ aussi sont
N N N N N N
bien définies, mais rien n’empêche que la constante de normalisation ⟨ηkN , gk ⟩ soit nulle, de sorte
que la distribution µN k = gk · ηk n’est pas nécessairement définie. Sur l’ensemble {τ
N N > k}
Le temps d’extinction est bien sûr infini dans le cas où les fonctions gk sont strictement
positives et ne s’annulent en aucun point, pour tout k = 0, 1, · · · , n.
Les résultats suivants seront démontrés dans ce chapitre et dans le chapitre suivant : une
borne non–asymptotique
cn,p
sup {E| ⟨µN
n − µn , ϕ⟩ | }
p 1/p
≤√ ,
ϕ : ∥ϕ∥=1 N
11.1. PROBABILITÉ D’EXTINCTION 155
On donne d’abord une majoration exponentielle de la probabilité d’extinction P[τ N ≤ n]. Cette
partie peut donc être sautée dans le cas où les fonctions gk sont strictement positives et ne
s’annulent en aucun point, pour tout k = 0, 1, · · · , n, puisque le temps d’extinction est infini
dans ce cas.
Proposition 11.1 Il existe des constantes positives an > 0 et bn > 0 telles que
P[τ N ≤ n] ≤ an exp{−bn N } ,
et
⟨γkN − γk , 1⟩
FkN = P[ | |> 1
et τ N > k−1 ] .
⟨γk , 1⟩ 2
Clairement, pour tout k = 0, 1, · · · , n, l’application c 7→ EkN (c) est décroissante et FkN ≤ EkN ( 12 ).
⟨γ0N − γ0 , 1⟩
Pour k = 0, si ⟨γ0N , 1⟩ = 0 alors nécessairement | | > 12 , de sorte que
⟨γ0 , 1⟩
⟨γ0N − γ0 , 1⟩
P[τ N = 0] = P[ ⟨γ0N , 1⟩ = 0 ] ≤ P[ | |> 1
] = F0N .
⟨γ0 , 1⟩ 2
156 CHAPITRE 11. ESTIMATION D’ERREUR
Pour tout k = 1, · · · , n, les bons ensembles {τ N > k} ⊆ {τ N > k − 1} sont emboı̂tés, et sur
⟨γ N − γk , 1⟩
l’ensemble {τ N > k−1}, si ⟨γkN , 1⟩ = 0 alors nécessairement | k | > 12 , de sorte que
⟨γk , 1⟩
⟨γkN − γk , 1⟩
P[τ N = k] = P[ ⟨γkN , 1⟩ = 0 et τ N > k−1 ] ≤ P[ | |> 1
et τ N > k−1 ] = FkN .
⟨γk , 1⟩ 2
On en déduit que
∑
n ∑
n
P[τ N
≤ n] = P[τ N
= k] ≤ FkN .
k=0 k=0
⟨γ0N − γ0 , ϕ⟩
P[ | | > c ] = P[ |⟨η0N − η0 , g0 ϕ⟩| > c ⟨η0 , g0 ⟩ ] .
⟨γ0 , 1⟩
où les v.a. (ξ0i , · · · , ξ0N ) sont indépendantes de distribution commune η0 , et on vérifie que
et
1 ∑ 1 ∑
N N
⟨η0N − η0 , g0 ϕ⟩ = [g0 (ξ0i ) ϕ(ξ0i ) − ⟨η0 , g0 ϕ⟩] = (Xi − E(Xi )) .
N N
i=1 i=1
2 c2
≤ 2 exp{− N} .
r02
Il en résulte que
⟨γ0N − γ0 , ϕ⟩ 2 c2
P[ | | > c ] ≤ 2 exp{− 2 N } ,
⟨γ0 , 1⟩ r0
et en prenant le supremum par rapport aux fonctions mesurables positives ϕ telles que ∥ϕ∥ = 1,
on obtient
2 c2
E0N (c) ≤ 2 exp{− 2 N } .
r0
11.1. PROBABILITÉ D’EXTINCTION 157
⟨γkN − γk , ϕ⟩ = ⟨γk−1
N
− γk−1 , Qk (gk ϕ)⟩ + ⟨ηkN − µN
k−1 Qk , gk ϕ⟩ ⟨γk−1 , 1⟩ ,
N
(11.3)
valide sur l’ensemble {τ N > k−1} pour toute fonction ϕ mesurable bornée, et on en déduit que
⟨γk−1
N −γ
k−1 , Qk (gk ϕ)⟩ ⟨η N − µN
k−1 Qk , gk ϕ⟩ ⟨γ N − γk−1 , 1⟩
≤ | |+| k | ( 1 + | k−1 |) .
⟨γk−1 , 1⟩ ⟨ηk , gk ⟩ ⟨ηk , gk ⟩ ⟨γk−1 , 1⟩
Si sur l’ensemble {τ N > k−1} (et a fortiori sur l’ensemble {τ N > k})
⟨γk−1
N −γ
k−1 , Qk (gk ϕ)⟩ ⟨γk−1
N −γ
k−1 , 1⟩
| |≤ 1
c ⟨ηk , gk ⟩ et | |≤ 1
,
⟨γk−1 , 1⟩ 2 ⟨γk−1 , 1⟩ 2
et
| ⟨ηkN − µN
k−1 Qk , gk ϕ⟩ | ≤
1
3 c ⟨ηk , gk ⟩ ,
⟨γkN − γk , ϕ⟩
| |≤c.
⟨γk , 1⟩
On en déduit que
⟨γkN − γk , ϕ⟩
EkN (c) = sup P[ | | > c et τ N > k ]
ϕ≥0 : ∥ϕ∥=1 ⟨γk , 1⟩
⟨γk−1
N −γ
k−1 , Qk (gk ϕ)⟩
≤ sup P[ | |> 1
c ⟨ηk , gk ⟩ et τ N > k−1 ]
ϕ≥0 : ∥ϕ∥=1 ⟨γk−1 , 1⟩ 2
(11.6)
⟨γk−1
N −γ
k−1 ,1⟩
+ P[ | ⟨γk−1 ,1⟩ |> 1
2 et τ N > k−1 ]
+ sup P[ | ⟨ηkN − µN
k−1 Qk , gk ϕ⟩ | >
1
3 c ⟨ηk , gk ⟩ et τ N > k−1 ] .
ϕ≥0 : ∥ϕ∥=1
N ,
Dans le second membre de (11.6), le deuxième terme s’interprète immédiatement comme Fk−1
et on se propose d’étudier successivement le premier et le troisième terme.
⟨γk−1
N −γ
k−1 , Qk (gk ϕ)⟩
P[ | |> 1
c ⟨ηk , gk ⟩ et τ N > k−1 ]
⟨γk−1 , 1⟩ 2
⟨γk−1
N −γ
k−1 , Qk (gk ϕ)⟩ ⟨ηk , gk ⟩
= P[ | |> 1
c et τ N > k−1 ]
⟨γk−1 , 1⟩ sup Qk (gk ϕ)(x) 2 sup Qk (gk ϕ)(x)
x∈E x∈E
⟨γk−1
N −γ
k−1 , Qk (gk ϕ)⟩ ⟨ηk , gk ⟩
≤ P[ | |> 1
c et τ N > k−1 ]
⟨γk−1 , 1⟩ sup Qk (gk ϕ)(x) 2 sup gk (x)
x∈E x∈E
⟨γk−1
N − γk−1 , ϕ⟩ c c
≤ sup P[ | |> 1
et τ N > k−1 ] = Ek−1
N
( 12 ) ,
ϕ≥0 : ∥ϕ∥=1 ⟨γk−1 , 1⟩ 2 rk rk
de sorte que
⟨γk−1
N −γ
k−1 , Qk (gk ϕ)⟩ c
sup P[ | |> 1
c ⟨ηk , gk ⟩ et τ N > k−1 ] ≤ Ek−1
N
( 12 ),
ϕ≥0 : ∥ϕ∥=1 ⟨γk−1 , 1⟩ 2 rk
ce qui fournit une majoration du premier terme figurant dans le second membre de (11.6).
Pour toute fonction mesurable positive ϕ telle que ∥ϕ∥ = 1, on définit
E(Xi | Fk−1
N
) = ⟨µN
k−1 Qk , gk ϕ⟩ ,
et
1 ∑
N
⟨ηkN − µN
k−1 Qk , gk ϕ⟩ = [ gk (ξki ) ϕ(ξki ) − ⟨µN
k−1 Qk , gk ϕ⟩ ]
N
i=1
1 ∑
N
= (Xi − E(Xi | Fk−1
N
)) .
N
i=1
c2
≤ 2 exp{− 2
9 N} .
rk2
11.1. PROBABILITÉ D’EXTINCTION 159
On en déduit que
P[ |⟨ηkN − µN
k−1 Qk , gk ϕ⟩ | >
1
3 c ⟨ηk , gk ⟩ et τ N > k−1 | Fk−1
N
]
=1 P[ |⟨ηkN − µN
k−1 Qk , gk ϕ⟩ | >
1
c ⟨ηk , gk ⟩ | Fk−1
N
]
(τ N > k−1) 3
c2
≤ 2 exp{− 2
9 N} ,
rk2
de sorte que, en prenant l’espérance
c2
P[ |⟨ηkN − µN
k−1 Qk , gk ϕ⟩ | >
1
3 c ⟨ηk , gk ⟩ et τ N > k−1 ] ≤ 2 exp{− 2
9 N} ,
rk2
et finalement, en prenant le supremum par rapport aux fonctions mesurables positives ϕ telles
que ∥ϕ∥ = 1, on obtient
c2
sup P[ |⟨ηkN − µN
k−1 Qk , gk ϕ⟩ | >
1
3 c ⟨ηk , gk ⟩ et τ N > k−1 ] ≤ 2 exp{− 2
9 N} ,
ϕ≥0 : ∥ϕ∥=1 rk2
ce qui fournit une majoration du troisième terme figurant dans le second membre de (11.6).
En reportant ces majorations dans le second membre de (11.6), on obtient
c c2
EkN (c) ≤ Ek−1
N
( 12 N
) + Fk−1 + 2 exp{− 2
9 N} ,
rk rk2
pour tout k = 1, · · · , n, avec la condition initiale
2 c2
E0N (c) ≤ 2 exp{− N} .
r02
On peut montrer par récurrence que
EkN (c) ≤ ek max(exp{−dk c2 N }, exp{−fk N }) = ek exp{− min(dk c2 , fk ) N }) ,
où ek > 0, dk > 0 et fk > 0 sont des réels positifs. En particulier pour c = 12 , on a
∑
n
≤ ( ek ) exp{− min min( 41 dk , fk ) N } ,
k=0,1,··· ,n
k=0
⟨γnN , 1⟩
{ E[ 1 N | − 1 |p ] }1/p ≤ znN,p , (11.7)
(τ > n) ⟨γn , 1⟩
et
sup { E[ 1 | ⟨µN
n − µn , ϕ⟩ | ] }
p 1/p
≤ 2 znN,p , (11.8)
ϕ : ∥ϕ∥=1 (τ N > n)
2 cp N,p 2 cp 2 cp
zkN,p ≤ rk (1 + √ ) zk−1 +√ et z0N,p ≤ √ . (11.9)
N N N
de sorte que
⟨γnN − γn , ϕ⟩ ⟨γ N , 1⟩
1 N | ⟨µN
n − µn , ϕ⟩ | ≤ 1(τ N > n) [ | | + ∥ϕ∥ | n − 1|] ,
(τ > n) ⟨γn , 1⟩ ⟨γn , 1⟩
⟨γnN , 1⟩ ⟨γ N − γn , ϕ⟩ p 1/p
{ E[ 1 N | − 1 |p ] }1/p ≤ sup { E[ 1 N | n | ]} ,
(τ > n) ⟨γn , 1⟩ ϕ : ∥ϕ∥=1 (τ > n) ⟨γn , 1⟩
et
⟨γnN − γn , ϕ⟩ p 1/p
sup { E[ 1 N | ⟨µN
n − µn , ϕ⟩ | ] }
p 1/p
≤2 sup { E[ 1 N | | ]} ,
ϕ : ∥ϕ∥=1 (τ > n) ϕ : ∥ϕ∥=1 (τ > n) ⟨γn , 1⟩
11.2. ESTIMATION D’ERREUR DANS LP 161
d’après l’inégalité (triangulaire) de Minkovski. Pour démonter le Théorème 11.2 il suffit donc de
prouver que la suite définie par
⟨γkN − γk , ϕ⟩ p 1/p
zkN,p = sup { E[ 1 | | ]} avec k = {τ
AN N
> k} ,
ϕ : ∥ϕ∥=1 AN
k ⟨γk , 1⟩
de sorte que
2 cp
z0N,p ≤ √ r0 .
N
Pour tout k = 1, · · · , n, on rappelle la décomposition
⟨γkN − γk , ϕ⟩ = ⟨γk−1
N
− γk−1 , Qk (gk ϕ)⟩ + ⟨ηkN − µN
k−1 Qk , gk ϕ⟩ ⟨γk−1 , 1⟩ ,
N
(11.3)
⟨γkN − γk , ϕ⟩ p 1/p
sup { E[1 | | ]}
ϕ : ∥ϕ∥=1 AN
k ⟨γk , 1⟩
⟨γk−1
N −γ
k−1 , Qk (gk ϕ)⟩ p
≤ sup { E[1 | | ] }1/p
ϕ : ∥ϕ∥=1 AN
k−1 ⟨γk , 1⟩
k−1 Qk , gk ϕ⟩ ⟨γk−1 , 1⟩ p
⟨ηkN − µN N
+ sup { E[ 1 | | ] }1/p ,
ϕ : ∥ϕ∥=1 AN
k−1 ⟨γk , 1⟩
{ E[ 1 | ⟨γk−1
N
− γk−1 , Qk (gk ϕ)⟩ |p ] }1/p
AN
k−1
pour tout x ∈ E et pour toute fonction mesurable bornée ϕ, et en divisant par ⟨γk , 1⟩ =
⟨ηk , gk ⟩ ⟨γk−1 , 1⟩ on obtient
⟨γk−1
N −γ
k−1 , Qk (gk ϕ)⟩ p
sup { E[ 1 | | ] }1/p
ϕ : ∥ϕ∥=1 AN
k−1 ⟨γk , 1⟩
sup gk (x)
x∈E ⟨γk−1
N −γ
k−1 , ϕ⟩ p
≤ sup { E[ 1 | | ] }1/p .
⟨ηk , gk ⟩ ϕ : ∥ϕ∥=1 AN
k−1 ⟨γk−1 , 1⟩
D’autre part
2 cp
{ E[ | ⟨ηkN − µN
k−1 Qk , gk ϕ⟩ | | Fk−1 ] }
p N 1/p
≤ √ sup gk (x) ∥ϕ∥ ,
N x∈E
{ E[ 1 | ⟨ηkN − µN
k−1 Qk , gk ϕ⟩ ⟨γk−1 , 1⟩ | ] }
N p 1/p
AN
k−1
2 cp
≤ √ sup gk (x) ∥ϕ∥ { E[ 1 N ⟨γk−1
N
, 1⟩p ] }1/p ,
N x∈E Ak−1
pour toute fonction mesurable bornée ϕ, et en divisant par ⟨γk , 1⟩ = ⟨ηk , gk ⟩ ⟨γk−1 , 1⟩ on obtient
⟨ηkN − µN
k−1 Qk , gk ϕ⟩ ⟨γk−1 , 1⟩ p
N
sup { E[ 1 | | ] }1/p
ϕ : ∥ϕ∥=1 AN
k−1 ⟨γk , 1⟩
sup gk (x)
2 cp x∈E ⟨γk−1
N , 1⟩
≤√ { E[ 1 | |p ] }1/p ,
N ⟨ηk , gk ⟩ AN
k−1 ⟨γk−1 , 1⟩
et on remarque que
⟨γk−1
N , 1⟩ ⟨γ N − γk−1 , 1⟩ p 1/p
{ E[ 1 | |p ] }1/p ≤ 1 + { E[ 1 N | k−1 | ]}
AN
k−1 ⟨γk−1 , 1⟩ Ak−1 ⟨γk−1 , 1⟩
N −γ
⟨γk−1 k−1 , ϕ⟩ p
≤ 1+ sup { E[ 1 | | ] }1/p ,
ϕ : ∥ϕ∥=1 AN
k−1 ⟨γk−1 , 1⟩
et la fonction de poids
∏
n
g0:n (x0:n ) = gk (xk−1 , xk ) ,
k=0
et si on définit
∫ ∫
⟨η0:n , g0:n f ⟩ = E[f (X0:n ) g0:n (X0:n ) ] = ··· f (x0:n ) g0:n (x0:n ) η0:n (dx0:n ) ,
E E
pour toute fonction mesurable bornée f définie sur l’espace produit E × · · · × E = E n+1 , alors
on peut réécrire le flot linéaire comme une intégrale et appliquer la méthode d’échantillonnage
pondéré vu à la Section 9.2. En effet, dans le cas particulier f = ϕ ◦ π où la fonction f ne dépend
que de la dernière variable, c’est–à–dire prend la forme suivante
on a
∏
n
⟨γn , ϕ⟩ = E[ϕ(Xn ) gk (Xk ) ] = E[ϕ ◦ π(X0:n ) g0:n (X0:n )] = ⟨η0:n , g0:n ϕ ◦ π⟩ ,
k=0
163
164 CHAPITRE 12. TCL POUR LES APPROXIMATIONS PARTICULAIRES
1 ∑
N
⟨γnN , ϕ⟩ = ⟨S (η0:n ), g0:n ϕ ◦ π⟩ =
N i
g0:n (ξ0:n ) ϕ ◦ π(ξ0:n
i
),
N
i=1
de sorte que
⟨γnN − γn , ϕ⟩ ⟨S N (η0:n ) − η0:n , g0:n ϕ ◦ π⟩
= ,
⟨γn , 1⟩ ⟨η0:n , g0:n ⟩
pour toute fonction mesurable bornée ϕ définie sur E, et il suffit d’appliquer le Théorème 9.11.
Théorème 12.1
√ ⟨γ N , 1⟩ √
N[ n − 1 ] =⇒ N(0, Vn ) et N ⟨µN
n − µn , ϕ⟩ =⇒ N(0, vn (ϕ)) ,
⟨γn , 1⟩
en distribution quand N ↑ ∞, pour toute fonction mesurable bornée ϕ, avec l’expression suivante
pour la variance asymptotique
⟨η0:n , g0:n
2 ⟩ ⟨η0:n , g0:n
2 |ϕ ◦ π − ⟨µ , ϕ⟩ |2 ⟩
n
Vn = −1 et vn (ϕ) = ,
⟨η0:n , g0:n ⟩2 ⟨η0:n , g0:n ⟩2
respectivement.
(12.2)
et γ0N = g0 S N (η 0) = g0 η0N ,
pour la distribution non–normalisée, alors il est facile de voir que
de sorte que
γkN γ0N
= gk · ηkN = µN
k et = g0 · η0N = µN
0 ,
⟨γkN , 1⟩ ⟨γ0N , 1⟩
c’est–à–dire que (12.2) correspond exactement à l’algorithme SIR avec ré–échantillonnage mul-
tinomial, et en itérant (12.3)
∏
n
⟨γn , 1⟩ =
N
⟨ηkN , gk ⟩ .
k=0
12.2. ÉCHANTILLONNAGE / RÉ–ÉCHANTILLONNAGE (SIR) 165
Théorème 12.2 Pour l’approximation particulaire du modèle (8.1), avec redistribution multi-
nomiale
√ ⟨γ N , 1⟩ √
N[ n − 1 ] =⇒ N(0, Vn ) et N ⟨µN
n − µn , ϕ⟩ =⇒ N(0, vn (ϕ)) ,
⟨γn , 1⟩
en distribution quand N ↑ ∞, pour toute fonction mesurable bornée ϕ, avec l’expression suivante
pour la variance asymptotique
∑
n
⟨ηk , (gk Rk+1:n 1)2 ⟩ ∑
n
⟨ηk , |gk Rk+1:n (ϕ − ⟨µn , ϕ⟩) |2 ⟩
Vn = [ − 1] et vn (ϕ) = ,
⟨ηk , gk Rk+1:n 1⟩2 ⟨ηk , gk Rk+1:n 1⟩2
k=0 k=0
respectivement, où
∏
n
Rk+1:n ϕ(x) = Rk+1 · · · Rn ϕ(x) = E[ϕ(Xn ) gp (Xp ) | Xk = x] ,
p=k+1
Corollaire 12.3 Pour l’approximation particulaire du modèle (8.1), avec redistribution multi-
nomiale √
N ⟨ηnN − ηn , ϕ⟩ =⇒ N(0, vn− (ϕ)) ,
en distribution quand N ↑ ∞, pour toute fonction mesurable bornée ϕ, avec l’expression suivante
pour la variance asymptotique
∑
n −
⟨ηk , |Rk+1:n (ϕ − ⟨ηn , ϕ⟩) |2 ⟩
vn− (ϕ) = − ,
k=0
⟨ηk , Rk+1:n 1⟩2
respectivement, où
∏
n−1
− −
Rk+1:n ϕ(x) = Rk+1 · · · Rn− ϕ(x) = E[ϕ(Xn ) gp (Xp ) | Xk = x] ,
p=k
−
pour tout k = 0, 1, · · · , n, avec la convention Rn+1:n ϕ(x) = ϕ(x).
pour toute fonction ϕ mesurable bornée. On vérifie que la v.a. ZN ′′ est mesurable par rapport
en distribution quand N ↑ ∞, pour toute fonction mesurable bornée ϕ, avec l’expression suivante
pour la variance asymptotique
∑
n−1
⟨ηk , |gk Rk+1:n−1 (Qn ϕ − ⟨µn−1 , Qn ϕ⟩) |2 ⟩
vn−1 (Qn ϕ) =
⟨ηk , gk Rk+1:n−1 1⟩2
k=0
∑
n−1 −
⟨ηk , |Rk+1:n (ϕ − ⟨ηn , ϕ⟩) |2 ⟩
= − ,
k=0
⟨ηk , Rk+1:n 1⟩2
et
∑
n−1 −
⟨ηk , |Rk+1:n (ϕ − ⟨ηn , ϕ⟩) |2 ⟩
= − + ⟨ηn , |ϕ − ⟨ηn , ϕ⟩ |2 ⟩
k=0
⟨ηk , Rk+1:n 1⟩2
∑
n −
⟨ηk , |Rk+1:n (ϕ − ⟨ηn , ϕ⟩) |2 ⟩
= − . 2
k=0
⟨ηk , Rk+1:n 1⟩2
pour toute fonction mesurable bornée ϕ, et compte tenu que ⟨γnN , 1⟩ −→ ⟨γn , 1⟩ en probabilité
quand N ↑ ∞, le Théorème 12.2 pour la distribution normalisée découle de (12.4) et du lemme
de Slutsky, avec
∑
n
var(gk Rk+1:n (ϕ − ⟨µn , ϕ⟩), ηk )
vn (ϕ) = Vn (ϕ − ⟨µn , ϕ⟩) = ,
⟨ηk , gk Rk+1:n 1⟩2
k=0
et on vérifie que
de sorte que
pour tout k = 0, 1, · · · , n.
∑
k
var(gp Rp+1:k ϕ, ηp ) ∑
k−1
var(gp Rp+1:k−1 Rk ϕ, ηp ) var(gk ϕ, ηk )
Vk (ϕ) = = + ,
⟨ηp , gp Rp+1:k 1⟩2 ⟨ηp , gp Rp+1:k−1 1⟩ ⟨ηk , gk ⟩
2 2 ⟨ηk , gk ⟩2
p=0 p=0
pour tout k = 1, · · · , n.
168 CHAPITRE 12. TCL POUR LES APPROXIMATIONS PARTICULAIRES
′′
√ ⟨γk−1
N −γ
k−1 , Rk ϕ⟩ Vk−1 (Rk ϕ)
ZN = N =⇒ N(0, ),
⟨γk , 1⟩ ⟨ηk , gk ⟩2
en distribution quand N ↑ ∞, compte tenu que ⟨γk , 1⟩ = ⟨ηk , gk ⟩ ⟨γk−1 , 1⟩. D’autre part, on
remarque que
′
√ ⟨ηkN − µN
k−1 Qk , gk ϕ⟩ ⟨γk−1 , 1⟩
N
ZN = N
⟨γk , 1⟩
√ ⟨ηkN − µN
k−1 Qk , gk ϕ⟩ ⟨γk−1 , 1⟩
N
= N
⟨ηk , gk ⟩ ⟨γk−1 , 1⟩
où conditionnellement par rapport à la tribu Fk−1 N engendrée par le système de particules jus-
qu’à la (k − 1)–ème génération, les variables aléatoires ξk1 , · · · , ξkN sont i.i.d. de distribution de
probabilité commune µN k−1 Qk . On vérifie que la variable aléatoire
var(gk ϕ, µNk−1 Qk )
s2i,N = E[ |Xi,N |2 | Fk−1
N
]= ,
⟨ηk , gk ⟩2
et bornée
sup gk (x)
x∈E
|Xi,N | ≤ 2 ∥ϕ∥ ,
⟨ηk , gk ⟩
pour tout i = 1, · · · , N . Clairement
1 ∑ 2
N
var(gk ϕ, µNk−1 Qk ) var(gk ϕ, ηk ) ⟨γk−1
N , 1⟩
s2N = si,N = −→ et θN = −→ 1 ,
N ⟨ηk , gk ⟩2 ⟨ηk , gk ⟩2 ⟨γk−1 , 1⟩
i=1
′ var(gk ϕ, ηk )
E[exp{j u ZN } | Fk−1
N
] −→ exp{− 12 u2 },
⟨ηk , gk ⟩2
en distribution quand N ↑ ∞, pour toute fonction mesurable bornée ϕ, avec l’expression suivante
pour la variance asymptotique
opt
∑
n
⟨µk , (Rk+1:n 1)2 ⟩
Vn−1 = [ − 1] ,
⟨µk , Rk+1:n 1⟩2
k=0
170 CHAPITRE 12. TCL POUR LES APPROXIMATIONS PARTICULAIRES
et
∑
n
⟨µk , |Rk+1:n (ϕ − ⟨µn , ϕ⟩) |2 ⟩
vnopt− (ϕ) = ,
⟨µk , Rk+1:n 1⟩2
k=0
respectivement, où
∏
n
Rk+1:n ϕ(x) = Rk+1 · · · Rn ϕ(x) = E[ϕ(Xn ) gp (Xp ) | Xk = x] ,
p=k+1
opt
∑
n−1
⟨ηkopt , (gkopt Rk+1:n−1
opt
1)2 ⟩ ∑
n
⟨µk , (Rk+1:n 1)2 ⟩
Vn−1 = [ − 1] = [ − 1] ,
⟨ηkopt , gkopt Rk+1:n−1
opt
, 1⟩2 ⟨µk , Rk+1:n 1⟩2
k=0 k=0
et
∑
n
⟨ηkopt , |Rk+1:n
opt−
(ϕ − ⟨ηnopt , ϕ⟩) |2 ⟩ ∑
n
⟨µk , |Rk+1:n (ϕ − ⟨µn , ϕ⟩) |2 ⟩
vnopt− (ϕ) = = ,
⟨ηkopt , Rk+1:n
opt−
1⟩2 ⟨µk , Rk+1:n 1⟩2
k=0 k=0
opt
∑
n
⟨µk , (Rk+1:n 1)2 ⟩ ∑
n
⟨ηk , gk ⟩ ⟨ηk , gk (Rk+1:n 1)2 ⟩
Vn−1 = [ − 1 ] = [ − 1] ,
⟨µk , Rk+1:n 1⟩2 ⟨ηk , gk Rk+1:n 1⟩2
k=0 k=0
à comparer avec
∑
n
⟨ηk , (gk Rk+1:n 1)2 ⟩
Vn = [ − 1] ,
⟨ηk , gk Rk+1:n 1⟩2
k=0
pour les constantes de normalisation, et
∑
n
⟨µk , |Rk+1:n (ϕ − ⟨µn , ϕ⟩) |2 ⟩ ∑
n
⟨ηk , gk ⟩ ⟨ηk , gk |Rk+1:n (ϕ − ⟨µn , ϕ⟩) |2 ⟩
vnopt− (ϕ) = = ,
⟨µk , Rk+1:n 1⟩2 ⟨ηk , gk Rk+1:n 1⟩2
k=0 k=0
12.2. ÉCHANTILLONNAGE / RÉ–ÉCHANTILLONNAGE (SIR) 171
à comparer avec
∑
n
⟨ηk , |gk Rk+1:n (ϕ − ⟨µn , ϕ⟩) |2 ⟩
vn (ϕ) = ,
⟨ηk , gk Rk+1:n 1⟩2
k=0
pour les distributions normalisées.
Dans le cas plus général du modèle (8.8) et pour une décomposition d’importance (8.11)
donnée, avec la représentation probabiliste (8.12) associée, ou bien pour le modèle (8.22) où la
décomposition d’importance est donnée de manière explicite dans la représentation probabiliste,
chaque fonction de sélection dépend de la transition courante de la chaı̂ne de Markov, mais il
suffit de changer de point de vue et d’adopter le modèle (8.27) à valeurs transitions, où chaque
fonction de sélection dépend seulement de l’état courant, puis de ré–exprimer dans ce cadre
le Théorème 12.2 établi ci–dessus pour le modèle (8.1) apparamment plus simple. On obtient
ainsi un premier résultat intermédiaire, qu’il suffit ensuite de ré–interpréter en terme de la
décomposition d’importance donnée. On introduit le noyau positif
pour tout k = 1, · · · , n.
Théorème 12.9 Pour l’approximation particulaire du modèle (8.22), avec redistribution mul-
tinomiale
√ ⟨γ N , 1⟩ √
N[ n − 1 ] =⇒ N(0, Vn ) et N ⟨µN
n − µn , ϕ⟩ =⇒ N(0, vn (ϕ)) ,
⟨γn , 1⟩
en distribution quand N ↑ ∞, pour toute fonction mesurable bornée ϕ, avec l’expression suivante
pour la variance asymptotique
∑
n
⟨µk−1 Rk□ , | Rk+1:n 1 |2 ⟩ ∑
n
⟨µk−1 R□ , | Rk+1:n (ϕ − ⟨µn , ϕ⟩) |2 ⟩
Vn = [ − 1] et vn (ϕ) = k
,
⟨µk−1 Rk , Rk+1:n 1⟩2 ⟨µk−1 Rk , Rk+1:n 1⟩2
k=0 k=0
respectivement, où
∏
n
Rk+1:n ϕ(x) = Rk+1 · · · Rn ϕ(x) = E[ϕ(Xn ) gp (Xp−1 , Xp ) | Xk = x] ,
p=k+1
en distribution quand N ↑ ∞, pour toute fonction mesurable bornée F définie sur l’ensemble
produit E × E, avec l’expression suivante pour la variance asymptotique
∑
n
⟨η tr , (gk Rk+1:n
tr 1)2 ⟩ ∑
n
⟨ηktr , |gk Rk+1:n
tr
n , F ⟩) | ⟩
(F − ⟨µtr 2
Vntr = [ ktr − 1] et vntr (F ) = ,
⟨ηk , gk Rk+1:n
tr 1⟩2 ⟨ηktr , gk Rk+1:n
tr 1⟩2
k=0 k=0
172 CHAPITRE 12. TCL POUR LES APPROXIMATIONS PARTICULAIRES
respectivement, où
∏
n
tr
Rk+1:n F (x1 , x2 ) = E[F (Xntr ) gp (Xptr ) | Xktr = (x1 , x2 )] ,
p=k+1
∫ ∫ ∫ ∫
⟨ηktr , F ⟩ = F (x′1 , x′2 ) µtr tr ′ ′
k−1 (dx1 , dx2 ) Qk (x1 , x2 , dx1 , dx2 )
E E E E
∫ ∫ ∫ ∫
= F (x′1 , x′2 ) µtr ′ ′ ′
k−1 (dx1 , dx2 ) δx2 (dx1 ) Qk (x1 , dx2 )
E E E E
∫ ∫ ∫
′ ′
= µtr
k−1 (dx1 , dx2 ) F (x2 , x2 ) Qk (x2 , dx2 )
E E E
∫ ∫
′ ′
= µtr
k−1 (E, dx2 ) F (x2 , x2 ) Qk (x2 , dx2 )
E E
∫ ∫
= µk−1 (dx2 ) F (x2 , x′2 ) Qk (x2 , dx′2 ) ,
E E
toute pour toute fonction mesurable bornée F définie sur l’ensemble produit E × E, et compte
tenu que µtr
k−1 (E, dx2 ) = µk−1 (dx2 ), c’est–à–dire que
ηktr = µk−1 ⊗ Qk .
∏
n
tr
Rk+1:n (ϕ ◦ π)(x1 , x2 ) = E[ϕ(Xn ) gp (Xp−1 , Xp ) | Xk−1 = x1 , Xk = x2 ]
p=k+1
∏
n
= E[ϕ(Xn ) gp (Xp−1 , Xp ) | Xk = x2 ]
p=k+1
= Rk+1:n ϕ(x2 ) ,
tr
Rk+1:n (ϕ ◦ π) = (Rk+1:n ϕ) ◦ π ,
12.2. ÉCHANTILLONNAGE / RÉ–ÉCHANTILLONNAGE (SIR) 173
−1 = µ , de sorte que
n ◦π
compte tenu que le résultat ne dépend que de x2 . On rappelle que µtr n
⟨µtr
n , ϕ ◦ π⟩ = ⟨µn , ϕ⟩ pour toute fonction de la forme F = ϕ ◦ π. On en déduit alors que
∫ ∫
= µk−1 (dx) Qk (x, dx′ ) |gk (x, x′ )|2 | Rk+1:n ϕ(x′ ) − Rk+1:n 1(x′ ) ⟨µn , ϕ⟩)) |2
E E
∫ ∫
= µk−1 (dx) Rk□ (x, dx′ ) | Rk+1:n ϕ(x′ ) − Rk+1:n 1(x′ ) ⟨µn , ϕ⟩)) |2
E E
et
∫ ∫
⟨ηktr , (gk Rk+1:n
tr
1)2 ⟩ = µk−1 (dx) Qk (x, dx′ ) |gk (x, x′ )|2 | Rk+1:n 1(x′ ) |2
E E
∫ ∫
= µk−1 (dx) Rk□ (x, dx′ ) | Rk+1:n 1(x′ ) |2
E E
et finalement
∫ ∫
⟨ηktr , gk Rk+1:n
tr
1⟩ = µk−1 (dx) Qk (x, dx′ ) gk (x, x′ ) Rk+1:n 1(x′ )
E E
∫ ∫
= µk−1 (dx) Rk (x, dx′ ) Rk+1:n 1(x′ )
E E
= ⟨µk−1 Rk , Rk+1:n 1⟩ ,
Inversion matricielle
H Q H∗ + R et H ∗ R−1 H + Q−1
− Q H ∗ (H Q H ∗ + R)−1 H
=I ,
= Q H ∗ (H Q H ∗ + R)−1 R . 2
175
176 ANNEXE A. INVERSION MATRICIELLE
Remarque A.2 Cette formule permet de remplacer l’inversion de la matrice (H ∗ R−1 H +Q−1 )
de dimension m, par l’inversion de la matrice (H Q H ∗ + R) de dimension d, avec d ≤ m en
général. En particulier, dans le cas où d = 1, la matrice H = h∗ est un vecteur ligne, la matrice
R = r est un scalaire, et la formule devient
h h∗ Q h h∗ Q
( + Q−1 )−1 = Q − .
r r + h∗ Q h
Remarque A.4 Si la matrice A est inversible, alors la matrice ∆ = D − C A−1 B est appelée
complément de Schur de la matrice A dans la matrice–bloc M , la matrice M est inversible si et
seulement si la matrice ∆ est inversible, et
⋆ ⋆
M −1 = .
⋆ ∆ −1
et on remarque que
−1
I 0 ∆ 0 I ⋆ ∆−1 0 I ⋆ ∆−1 ⋆
= = .
⋆ I 0 ⋆ 0 I ⋆ ⋆ 0 I ⋆ ⋆
de sorte que ( )
u∗ −u∗ B D−1 A B u
= u∗ ∆ u ,
B∗ D −D−1 B ∗ u
pour tout vecteur u, ce qui permet de conclure. 2
178 ANNEXE A. INVERSION MATRICIELLE
Annexe B
Inégalités
On regroupe dans cette annexe plusieurs résultats non–asymptotiques sur les sommes de va-
riables aléatoires indépendantes mais pas nécessairement identiquement distribuées : inégalité
de Khintchine, inégalité de Marcinkiewicz–Zygmund pour les moments d’ordre p ≥ 1, inégalité
exponentielle de Hoeffding pour les probabilités de déviation.
On appelle suite de Rademacher une suite de variables aléatoires indépendantes prenant les
valeurs −1 ou +1 avec probabilité 12 .
Proposition B.1 (Inégalité de Khintchine) Pour tout réel p ≥ 0, il existe une constante
positive Ap > 0 telle que
∑
N ∑N
E| εi ci | ≤ Ap (
p
c2i )p/2 ,
i=1 i=1
pour toute suite (c1 , · · · , cN ) de réels et pour toute suite de Rademacher (ε1 , · · · , εN ).
∑
N
Par homogénéité, on peut supposer que c2i = 1 sans perte de généralité.
i=1
Si l’inégalité est vraie pour un entier p ≥ 0, alors pour tout réel 0 ≤ q ≤ p
∑
N ∑
N
E| εi ci | ≤ {E|
q
εi ci |p }q/p ≤ Aq/p
p ,
i=1 i=1
d’après l’inégalité de Jensen, compte tenu que l’application x 7→ |x|q/p est concave, et il suffit de
montrer l’inégalité pour tout entier p ≥ 1, le cas p = 0 étant trivial.
Preuve. On remarque que
1
ex ≥ e|x| 1(x ≥ 0) ≥ |x|p 1(x ≥ 0) ,
p!
pour tout entier p ≥ 1, de sorte que
1 1 1
|x|p = |x|p 1(x ≥ 0) + |x|p 1(x ≤ 0) ≤ ex + e−x ,
p! p! p!
179
180 ANNEXE B. INÉGALITÉS
1 ∑
N ∑N ∑
N ∑N
E| εi ci |p ≤ E exp{ εi ci } + E exp{− εi ci } = 2 E exp{ εi ci } .
p!
i=1 i=1 i=1 i=1
Finalement
∑
N ∏
N ∑
N
√
E exp{ εi ci } = ( 12 eci + 12 e−ci ) ≤ exp{ 12 c2i } = e,
i=1 i=1 i=1
1
compte tenu de l’inégalité 12 (ex +e−x ) ≤ e 2
x2
, valide pour tout réel x, et l’inégalité de Khintchine
√
est démontrée avec Ap = 2 e p!. 2
∑
N ∑N
E| Xi | ≤ Bp E(
p
Xi2 )p/2 ,
i=1 i=1
1 ∑ 1 ∑ 2 p/2
N N
Bp
E| Xi |p ≤ p/2 E( Xi ) ,
N N N
i=1 i=1
1 ∑ 2 p/2 1 ∑
N N
( Xi ) ≤ |Xi |p ,
N N
i=1 i=1
1 ∑ Bp 1 ∑
N N
E| Xi |p ≤ p/2 ( E|Xi |p ) . (B.1)
N N N
i=1 i=1
Remarque B.4 Plus généralement, pour tout vecteur de probabilité (w1 , · · · , wN ), et quitte à
remplacer Xi par wi Xi pour tout i = 1, · · · , N , on obtient
∑
N ∑N ∑N ∑
N
E| wi Xi | ≤ Bp E(
p 2 2 p/2
wi X i ) = B p ( wi ) E(
2 p/2
wi□ Xi2 )p/2 ,
i=1 i=1 i=1 i=1
wi2
wi□ = ,
∑
N
wj2
j=1
181
où (ε1 , · · · , εN ) est une suite de Rademacher, où la suite (X1′ , · · · , XN ′ ) a la même distribution
′
que la suite (X1 , · · · , XN ), et où les suites (X1 , · · · , XN ), (X1 , · · · , XN′ ) et (ε , · · · , ε ) sont
1 N
mutuellement indépendantes. Pour tout i = 1, · · · , N et compte tenu que les variables aléatoires
(Xi − Xi′ ) et (Xi′ − Xi ) ont la même distribution, on vérifie que
E[ϕ(εi (Xi − Xi′ ))] = E( E[ϕ(εi (Xi − Xi′ )) | Xi , Xi′ ] )
∑
N ∑
N ∑
N
≤ 2p−1 ( E| εi Xi |p + E| εi Xi′ |p ) = 2p E| εi Xi |p ,
i=1 i=1 i=1
182 ANNEXE B. INÉGALITÉS
et
∑
N ∑
N ∑N
E| εi Xi | = E( E[ |
p
εi Xi | | X1 , · · · , XN ] ) ≤ Ap E(
p
Xi2 )p/2 ,
i=1 i=1 i=1
∑
N
sym p sym p
∑N
E| Xi | = E|SN | ≤
p p
E|SN | = E|RN | = 2 Ap E(
p
Xi2 )p/2 ,
i=1 i=1
Lemme B.5 Soit X une variable aléatoire réelle, de moyenne nulle et à valeurs bornées, c’est–
à–dire que a ≤ X ≤ b. Alors
x−a b−x
x= b+ a,
b−a b−a
et de la convexité de la fonction exponentielle que
x−a b−x
exp{s x} ≤ exp{s b} + exp{s a} ,
b−a b−a
pour tout a ≤ x ≤ b, de sorte que
−a b
E[exp{s X}] ≤ exp{s b} + exp{s a} ,
b−a b−a
compte tenu que E(X) = 0. On pose
−a b
p= de sorte que 1−p= ,
b−a b−a
et il vient
s a = −s p (b − a) = −p u et s b = s (1 − p) (b − a) = (1 − p) u ,
ce qui définit
ϕ(u) = −p u + log(p exp{u} + 1 − p) .
183
1 ∑
N
2 N c2
P[ | (Xi − E(Xi ))| ≥ c ] ≤ 2 exp{− },
1 ∑
N N
i=1
(bi − ai )2
N
i=1
Preuve. On utilise la méthode de majoration de Chernoff : pour tout réel positif λ > 0, il
résulte de l’inégalité de Markov, de l’indépendance des variables aléatoires (X1 , · · · , XN ) et du
Lemme B.5, que
∑
N ∑
N
P[ (Xi − E(Xi )) ≥ c ] = P[ exp{λ (Xi − E(Xi ))} ≥ exp{λ c} ]
i=1 i=1
∑
N
≤ exp{−λ c} E[exp{λ (Xi − E(Xi ))}]
i=1
∏
N
≤ exp{−λ c} E[exp{λ (Xi − E(Xi ))}]
i=1
∏
N
≤ exp{−λ c} exp{ 81 λ2 (bi − ai )2 }
i=1
∑
N
≤ exp{−λ c + 81 λ2 (bi − ai )2 } ,
i=1
184 ANNEXE B. INÉGALITÉS
∑
N ∑
N
2 c2
P[ (Xi − E(Xi )) ≥ c ] ≤ min exp{−λ c + 8 λ
1 2
(bi − ai )2 } = exp{− },
i=1
λ>0
i=1
∑
N
(bi − ai ) 2
i=1
4c
obtenue pour la valeur λ = . En posant Xi′ = −Xi et compte tenu que
∑
N
(bi − ai )2
i=1
pour tout i = 1, · · · , N (de sorte que l’oscillation est invariante par changement de signe), il
vient
∑
N ∑
N
2 c2
P[ (Xi − E(Xi )) ≤ −c ] = P[ (Xi′ − E(Xi′ )) ≥ c ] ≤ exp{− },
i=1 i=1
∑
N
(bi − ai )2
i=1
∑
N ∑
N ∑
N
P[ | (Xi − E(Xi ))| ≥ c ] = P[ (Xi − E(Xi )) ≥ c ou (Xi − E(Xi )) ≤ −c ]
i=1 i=1 i=1
∑
N ∑
N
= P[ (Xi − E(Xi )) ≥ c ] + P[ (Xi − E(Xi )) ≤ −c ]
i=1 i=1
2 c2
≤ 2 exp{− }. 2
∑
N
(bi − ai )2
i=1
Annexe C
On regroupe dans cette annexe quelques généralisations du théorème central limite, dont la
version la plus classique concerne la somme de variables aléatoires i.i.d. (indépendantes et iden-
tiquement distribuées). On commence par rappeler des majorations bien connues.
Lemme C.1 On a
Preuve. On définit
G1 (x) = e−x − (1 − x) ,
pour tout x ≥ 0, et on vérifie que
pour tout x ≥ 0, de sorte que la fonction G1 est croissante sur [0, ∞) et G1 (x) ≥ G1 (0) = 0, ce
qui prouve la minoration. On définit ensuite
G2 (x) = e−x − (1 − x) − 21 x2 ,
pour tout x ≥ 0, de sorte que la fonction G2 est décroissante sur [0, ∞) et G2 (x) ≤ G2 (0) = 0,
ce qui prouve la majoration.
185
186 ANNEXE C. THÉORÈME CENTRAL LIMITE CONDITIONNEL
On définit
′
F (λ) = ej (λ x+(1−λ) x ) ,
pour tout 0 ≤ λ ≤ 1, où x et x′ sont fixés. On vérifie que
′
F ′ (λ) = j (x − x′ ) ej (λ x+(1−λ) x ) ,
et on en déduit que
∫ 1 ∫ 1
j x′ ′ ′ ′
e jx
−e = F (1) − F (0) = F (λ) dλ = j (x − x ) ej (λ x +(1−λ) x) dλ ,
0 0
de sorte que ∫ 1
j x′ ′ ′
|e jx
−e | = |x − x | | ej (λ x +(1−λ) x) dλ| ≤ |x − x′ | .
0
On définit de même
′
F (λ) = e−(λ x+(1−λ) x ) ,
pour tout 0 ≤ λ ≤ 1, où x, x′ ≥ 0 sont fixés. On vérifie que λ x + (1 − λ) x′ ≥ 0 et
′
F ′ (λ) = −(x − x′ ) e−(λ x+(1−λ) x ) ,
et on en déduit que
∫ 1 ∫ 1
′ ′
e−x − e−x = F (1) − F (0) = F ′ (λ) dλ = −(x − x′ ) e−(λ x+(1−λ) x ) dλ ,
0 0
de sorte que
∫ 1
−x −x′ ′ ′
|e −e | = |x − x | | e−(λ x+(1−λ) x ) dλ| ≤ |x − x′ | . 2
0
il vient
∫ 1
E[ exp{j u X}] = 1 − u s − u 1
2
2 2 2
(1 − λ) E[ |X|2 [exp{j λ u X} − 1] ] dλ ,
0
de sorte que
∫ 1
1
|E[ exp{j u X}] − (1 − 12 u2 s2 )| ≤ E[ |X|2 | exp{j λ u X} − 1| ] dλ .
u2 0
187
On définit
Z(u) = |X|2 | exp{j λ u X} − 1| ≤ 2 |X|2 ,
pour tout réel u, de sorte que la famille (Z(u) , u ∈ R) est uniformément intégrable, et on vérifie
que Z(u) converge vers 0 presque sûrement quand u → 0. Il suffit alors d’appliquer le théorème
de convergence dominée de Lebesgue pour conclure. 2
Lemme C.3 Soit X une variable aléatoire centrée, de variance s2 . Pour tout réel positif c > 0
et pour tout réel u
≤ 1
6 c s2 |u|3 + E[ 1(|X| > c) |X|2 ] u2 .
de sorte que
∫ 1
(1 − λ) |ej λ x − 1| dλ = min( 61 |x|, 1) .
0
il vient
∫ 1
E[ exp{j u X}] = 1 − u s − u 1
2
2 2 2
(1 − λ) E[ |X|2 [exp{j λ u X} − 1] ] dλ ,
0
de sorte que
∫ 1
|E[ exp{j u X}] − (1 − u s )| ≤ u 1
2
2 2 2
(1 − λ) E[ |X|2 | exp{j λ u X} − 1| ] dλ
0
∫ 1
= u E[ |X|
2 2
(1 − λ) | exp{j λ u X} − 1| dλ ]
0
1
|E[ exp{j u X}] − (1 − 12 u2 s2 )| ≤ E[ |X|2 min( 16 |u X|, 1) ] ,
u2
fournit une autre preuve du Lemme C.2. On définit
pour tout réel u, de sorte que la famille (Z(u) , u ∈ R) est uniformément intégrable, et on vérifie
que Z(u) converge vers 0 presque sûrement quand u → 0. Il suffit alors d’appliquer le théorème
de convergence dominée de Lebesgue pour conclure.
Théorème C.5 Soit (X1 , · · · , XN ) des variables aléatoires i.i.d. centrées et de variance s2 .
Alors
1 ∑
N
√ Xi =⇒ N(0, s2 ) ,
N i=1
en distribution, quand N ↑ ∞.
∑
N ∏
N
E[exp{j u Xi }] = E[exp{j u Xi } = ( E[exp{j u X}] )N ,
i=1 i=1
où la variable aléatoire X est distribuée comme chacune des variables aléatoires (X1 , · · · , XN ).
D’après les majorations classiques rappelées dans les Lemmes C.2 et C.1-(i), on a
et
| exp{− 21 u2 s2 } − (1 − 21 u2 s2 ) | ≤ 1
8 u4 s4 ,
En utilisant la majoration
|aN − bN | ≤ N |a − b| ,
on obtient
∑
N
| E[exp{j u Xi }] − exp{− 12 N u2 s2 } | = | ( E[exp{j u X} )N − ( exp{− 12 u2 s2 } )N |
i=1
≤ N | E[exp{j u X}] − exp{− 12 u2 s2 } |
≤ N R(u) u2 + 81 N s4 u4 ,
v
et en posant u = √ on obtient
N
v ∑
N
v s4 v 4
| E[exp{j √ Xi }] − exp{− 12 s2 v 2 } | ≤ R( √ ) v 2 + 1
8 .
N i=1 N N
v ∑
N
E[exp{j √ Xi }] −→ exp{− 12 s2 v 2 } ,
N i=1
quand N ↑ ∞. 2
Théorème C.6 Soit (X1,N , · · · , XN,N ) des variables aléatoires indépendantes et centrées. On
pose
s2i,N = E|Xi,N |2 et pour tout c > 0 Fi,N (c) = E[ 1(|X > c) |Xi,N | ] ,
2
i,N |
pour tout i = 1, · · · , N . Si
1 ∑ 2
N
s2N = si,N −→ s2 ,
N
i=1
1 ∑ √ ∑
N N
Xi,N
FN (ε) = Fi,N (ε N )) = E[ 1 Xi,N | √ |2 ] −→ 0 , (C.1)
N (| √ | > ε) N
i=1 i=1
N
quand N ↑ ∞, alors
1 ∑
N
√ Xi,N =⇒ N(0, s2 ) ,
N i=1
en distribution, quand N ↑ ∞.
190 ANNEXE C. THÉORÈME CENTRAL LIMITE CONDITIONNEL
D’après les majorations classiques rappelées dans les Lemmes C.3 et C.1-(i), on a
et
| exp{− 12 u2 s2i,N } − (1 − 12 u2 s2i,N ) | ≤ 1
8 u4 s4i,N ,
et d’après l’inégalité triangulaire, on a
En utilisant la majoration
∑
N ∑
N
|a1 · · · an − b1 · · · bN | = | a1 · · · ai−1 (ai − bi ) bi+1 · · · bN | ≤ |ai − bi | ,
i=1 i=1
∏
N ∏
N
=| E[exp{j u Xi,N }] − exp{− 12 u2 s2i,N } |
i=1 i=1
∑
N
≤ | E[exp{j u Xi,N }] − exp{− 12 u2 s2i,N } |
i=1
∑
N ∑
N ∑
N
≤ 1
6 c s2i,N |u|3 + Fi,N (c) u2 + 1
8 s4i,N u4 .
i=1 i=1 i=1
On remarque que
∑
N
s2i,N = E[ 1(|X |Xi,N |2 ] + E[ 1(|X | > c) |Xi,N |2 ] ≤ c2 + Fi,N (c) ≤ c2 + Fi,N (c) ,
i,N | ≤ c) i,N
i=1
On obtient ainsi
∑
N
| E[exp{j u Xi,N }] − exp{− 12 N u2 s2N } |
i=1
∑
N ∑
N
≤ 1
6 c N s2N |u|3 + Fi,N (c) u2 + 1
8 (c2 + Fi,N (c)) N s2N u4 ,
i=1 i=1
v √
et en posant u = √ et c = ε N on obtient
N
v ∑
N
| E[exp{j √ Xi,N }] − exp{− 12 v 2 s2N } |
N i=1 (C.2)
≤ 1
6 ε s2N |v|3 + FN (ε) v 2 + 1
8 (ε2 + FN (ε)) s2N v 4 4 ,
v ∑
N
| E[exp{j √ Xi,N }] − exp{− 12 v 2 s2 } |
N i=1
≤ 1
6 ε s2N |v|3 + FN (ε) v 2 + 1
8 (ε2 + FN (ε)) s2N v 4 + 21 v 2 | s2N − s2 | .
v ∑
N
lim sup | E[exp{j √ Xi,N }] − exp{− 12 v 2 s2 } | ≤ 1
6 ε s2 |v|3 + 1
8 ε2 s 2 v 4 ,
N ↑∞ N i=1
v ∑
N
E[exp{j √ Xi,N }] −→ exp{− 21 v 2 s2 } ,
N i=1
quand N ↑ ∞. 2
Théorème C.7 On suppose que conditionnellement par rapport à FN , les variables aléatoires
(X1,N , · · · , XN,N ) sont indépendantes et centrées, et on pose
s2i,N = E[ |Xi,N |2 | FN ] et pour tout c > 0 Fi,N (c) = E[ 1(|X > c) |Xi,N | | FN ] ,
2
i,N |
1 ∑ 2
N
s2N = si,N −→ s2 et θN −→ 1 ,
N
i=1
192 ANNEXE C. THÉORÈME CENTRAL LIMITE CONDITIONNEL
1 ∑ √ ∑
N N
Xi,N
FN (ε) = Fi,N (ε N )) = E[ 1 Xi,N | √ |2 | FN ] −→ 0 , (C.3)
N (| √ | > ε) N
i=1 i=1
N
en probabilité quand N ↑ ∞, alors pour tout réel v
v θN ∑
N
E[exp{j √ Xi,N } | FN ] −→ exp{− 12 v 2 s2 } , (C.4)
N i=1
en probabilité quand N ↑ ∞.
∑
N ∏
N
E[exp{j u Xi,N } | FN ] = E[exp{j u Xi,N } | FN ] .
i=1 i=1
v θN ∑
N
| E[exp{j √ Xi,N } | FN ] − exp{− 12 v 2 θN s }|
2 2
N i=1
≤ 1
6 ε s2N θN
3
|v|3 + FN (ε) θN
2 2
v + 1
8 (ε2 + FN (ε)) s2N θN
4 4
v .
| exp{− 12 v 2 θN sN } − exp{− 12 v 2 s2 } | ≤
2 2 1
2 v 2 | θN sN − s2 | ,
2 2
v θN ∑
N
| E[exp{j √ Xi,N } | FN ] − exp{− 12 v 2 s2 } |
N i=1
≤ 1
6 ε s2N θN
3
|v|3 + FN (ε) θN
2 2
v + 1
8 (ε2 + FN (ε)) s2N θN v + 21 v 2 | θN
4 4
sN − s2 | .
2 2
∆N (ε) = 1
6 ε s2N θN
3
|v|3 + FN (ε) θN
2 2
v + 1
8 (ε2 + FN (ε)) s2N θN v + 12 v 2 | θN
4 4
sN − s2 |
2 2
−→ ∆(ε) = 1
6 ε s2 |v|3 + 1
8 ε2 s2 v 4 ,
en probabilité quand N ↑ ∞. Soit η > 0 fixé. On rappelle que le réel v est aussi fixé. Il existe
alors ε = ε(η) > 0 tel que
∆(ε) = 61 ε |v|3 + 18 ε2 v 4 < 12 η .
C.2. TCL POUR DES VARIABLES ALÉATOIRES INDÉPENDANTES 193
v θN ∑
N
| E[exp{j √ Xi,N } | FN ] − exp{− 12 v 2 s2 } |
N i=1
On en déduit que
v θN ∑
N
P[ | E[exp{j √ Xi,N } | FN ] − exp{− 12 v 2 s2 } | > η] ≤ P[ |∆N (ε) − ∆(ε)| > 1
2 η] ,
N i=1
de sorte que
v θN ∑
N
E[exp{j √ Xi,N } | FN ] −→ exp{− 21 v 2 s2 } ,
N i=1
en probabilité quand N ↑ ∞. 2
Remarque C.8 Si les variables aléatoires (X1,N , · · · , XN,N ) sont bornées, i.e. si |Xi,N | ≤ K
pour tout i = 1, · · · , N , alors
∑
N
Xi,N K2 ∑
N
Xi,N
FN (ε) = E[ 1 Xi,N | √ |2 | FN ] ≤ P[ | √ | > ε | FN ] ,
(| √ | > ε) N N N
i=1 i=1
N
et on remarque que
Xi,N √
P[ | √ | > ε | FN ] ≤ 1 ,
N (K > ε N )
pour tout i = 1, · · · , N , de sorte que
FN (ε) ≤ K 2 1 √ −→ 0 ,
(K > ε N )
Remarque C.9 Si les variables aléatoires (X1,N , · · · , XN,N ) vérifient la condition de Lyapunov
conditionnelle : pour un certain δ > 0
∑
N
Xi,N
E[ | √ |2+δ | FN ] −→ 0 , (C.5)
i=1
N
en probabilité quand N ↑ ∞, il est facile de montrer directement le résultat (C.4), sans devoir
passer par l’intermédiaire du Théorème C.7.
La première généralisation concerne un théorème central limite conditionnel pour une somme
de variables aléatoires i.i.d., et son application à la convergence en distribution de la somme de
deux variables aléatoires
• une somme de variables aléatoires indépendantes mesurables par rapport à une sous–tribu,
• et une somme de variables aléatoires conditionnellement indépendantes par rapport à la
sous–tribu.
distribution vers une variable aléatoire gaussienne centrée, de variance V ′ , au sens où pour tout
u fixé
′
E[exp{j u ZN } | FN ] −→ exp{− 21 u2 V ′ } ,
en probabilité (et dans L1 , par le théorème de convergence dominée de Lebesgue) quand N ↑ ∞,
et si la variable aléatoire ZN′′ est mesurable par rapport à F , et converge en distribution vers
N
une variable aléatoire gaussienne centrée, de variance V ′′ , i.e. si pour tout u fixé
′′
E[exp{j u ZN }] −→ exp{− 12 u2 V ′′ } ,
quand N ↑ ∞, alors la variable aléatoire ZN = ZN ′ + Z ′′ converge en distribution vers une
N
variable aléatoire gaussienne centrée, de variance V = V ′ + V ′′ , quand N ↑ ∞.
+ exp{− 21 u2 V ′ } [ E[ exp{j u ZN
′′
}] − exp{− 21 u2 V ′′ }] ,
C.3. TCL CONDITIONNEL 195
| E[exp{j u ZN }] − exp{− 12 u2 V } |
′
≤ E| E[exp{j u ZN } | FN ] − exp{− 21 u2 V ′ } |
′′
+ | E[exp{j u ZN }] − exp{− 12 u2 V ′′ } | ,
[1] Brian D. O. Anderson and John B. Moore. Optimal filtering. Prentice–Hall Information
and System Sciences Series. Prentice–Hall, Englewood Cliffs, NJ, 1979.
[2] M. Sanjeev Arulampalam, Simon Maskell, Neil J. Gordon, and Tim Clapp. A tutorial on
particle filters for online nonlinear / non–Gaussian Bayesian tracking. IEEE Transactions
on Signal Processing, SP–50(2 (Special issue on Monte Carlo Methods for Statistical Signal
Processing)) :174–188, February 2002.
[3] Nathalie Bartoli and Pierre Del Moral. Simulation et algorithmes stochastiques. Cépaduès,
Toulouse, 2001.
[4] Niclas Bergman. Posterior Cramér–Rao bounds for sequential estimation. In Arnaud Dou-
cet, Nando de Freitas, and Neil Gordon, editors, Sequential Monte Carlo methods in practice,
Statistics for Engineering and Information Science, chapter 15, pages 321–338. Springer–
Verlag, New York, 2001.
[5] Olivier Cappé, Simon J. Godsill, and Éric Moulines. An overview of existing methods and
recent advances in sequential Monte Carlo. Proceedings of the IEEE, 95(5 (Special issue on
Large–Scale Dynamic Systems)) :899–924, May 2007.
[6] Olivier Cappé, Éric Moulines, and Tobias Rydén. Inference in hidden Markov models.
Springer Series in Statistics. Springer–Verlag, New York, 2005.
[7] Dan Crişan and Arnaud Doucet. A survey of convergence results on particle filtering
methods for practitioners. IEEE Transactions on Signal Processing, 50(3) :736–746, March
2002.
[8] Pierre Del Moral. Feynman–Kac formulae. Genealogical and interacting particle systems
with applications. Probability and its Applications. Springer–Verlag, New York, 2004.
[9] Luc Devroye. Non–uniform random variate generation. Springer–Verlag, New York, 1986.
[10] Randal Douc, Olivier Cappé, and Éric Moulines. Comparison of resampling schemes for
particle filtering. In Proceedings of the 4th Symposium on Image and Signal Processing and
Analysis, Zagreb 2005, pages 64–69. IEEE–SPS, September 2005.
[11] Randal Douc and Éric Moulines. Limit theorems for weighted samples with applications
to sequential Monte Carlo methods. The Annals of Statistics, 36(5) :2344–2376, October
2008.
197
198 BIBLIOGRAPHIE
[12] Randal Douc, Éric Moulines, and David S. Stoffer. Nonlinear time series : Theory, methods
and applications with R examples. Texts in Statistical Science. Chapman & Hall / CRC
Press, Boca Raton, 2014.
[13] Arnaud Doucet and Christophe Andrieu. Particle filters for partially observed Gaussian
state space models. Journal of the Royal Statistical Society, Series B, 64(4) :827–836,
December 2002.
[14] Arnaud Doucet, Nando de Freitas, and Neil Gordon, editors. Sequential Monte Carlo
methods in practice. Statistics for Engineering and Information Science. Springer–Verlag,
New York, 2001.
[15] Arnaud Doucet, Simon J. Godsill, and Christophe Andrieu. On sequential Monte Carlo
sampling methods for Bayesian filtering. Statistics and Computing, 10(3) :197–208, July
2000.
[16] Fredrik Gustafsson, Fredrik Gunnarsson, Niclas Bergman, Urban Forssell, Jonas Jansson,
Rickard Karlsson, and Per-Johan Nordlund. Particle filters for positioning, navigation, and
tracking. IEEE Transactions on Signal Processing, SP–50(2 (Special issue on Monte Carlo
Methods for Statistical Signal Processing)) :425–437, February 2002.
[17] Allan Gut. Probability : A graduate course. Springer Texts in Statistics. Springer–Verlag,
New York, 2005.
[18] Hans R. Künsch. Recursive Monte Carlo filters : Algorithms and theoretical analysis. The
Annals of Statistics, 33(5) :1983–2021, October 2005.
[19] Jun S. Liu. Monte Carlo strategies in scientific computing. Springer Series in Statistics.
Springer–Verlag, New York, 2001.
[20] Patrick Pérez, Carine Hue, Jaco Vermaak, and Michel Gangnet. Color–based probabilistic
tracking. In Anders Heyden, Gunnar Sparr, Mads Nielsen, and Peter Johansen, editors,
Proceedings of the 7th European Conference on Computer Vision (ECCV’02), Copenhagen
2002, volume 2350 of Lecture Notes in Computer Science, pages 661–675. Springer–Verlag,
Berlin, June 2002.
[21] Dinh-Tuan Pham. Stochastic methods for sequential data assimilation in strongly nonlinear
systems. Monthly Weather Review, 129(5) :1194–1207, May 2001.
[22] Branko Ristić, M. Sanjeev Arulampalam, and Neil J. Gordon. Beyond the Kalman filter :
Particle filters for tracking applications. Artech House, Boston, 2004.
[23] Christian P. Robert and George Casella. Monte Carlo statistical methods. Springer Texts
in Statistics. Springer–Verlag, New York, 2nd edition, 2004.
[24] Thomas Schön, Fredrik Gustafsson, and Per-Johan Nordlund. Marginalized particle filters
for mixed linear / nonlinear state–space models. IEEE Transactions on Signal Processing,
SP–53(7) :2279–2289, July 2005.
[25] Sebastian Thrun, Wolfram Burgard, and Dieter Fox. Probabilistic robotics. Intelligent
Robotics and Autonomous Agents. The MIT Press, Cambridge, MA, 2005.
BIBLIOGRAPHIE 199
[26] Petr Tichavský, Carlos H. Muravchik, and Arye Nehorai. Posterior Cramér–Rao bounds for
discrete–time nonlinear filtering. IEEE Transactions on Signal Processing, SP–46(5) :1386–
1396, May 1998.
[27] Miroslav Šimandl, Jakub Královec, and Petr Tichavský. Filtering, predictive and smoo-
thing Cramér–Rao bounds for discrete–time nonlinear dynamic systems. Automatica,
37(11) :1703–1716, November 2001.