Fichier PDF

Partage, hébergement, conversion et archivage facile de documents au format PDF

Partager un fichier Mes fichiers Convertir un fichier Boite à outils PDF Recherche PDF Aide Contact



PolyCoursDB .pdf



Nom original: PolyCoursDB.pdf

Ce document au format PDF 1.5 a été généré par LaTeX with hyperref package / pdfTeX-1.40.12, et a été envoyé sur fichier-pdf.fr le 31/07/2013 à 20:13, depuis l'adresse IP 142.137.x.x. La présente page de téléchargement du fichier a été vue 1952 fois.
Taille du document: 2.3 Mo (75 pages).
Confidentialité: fichier public




Télécharger le fichier (PDF)









Aperçu du document


Master Recherche Optique, Photonique, Signal et Image
Option Signal-Image

Introduction aux Modèles Markoviens
pour le Signal et l’Image

Stéphane DERRODE
stephane.derrode@centrale-marseille.fr

Version 1.3 - Janvier 2012

2

Notes bibliographiques

Ce cours présente une introduction aux modèles markoviens pour le traitement numérique du signal et de
l’image. Les modèles statistiques exposés ici sont classiques (modèle aveugle, chaîne de Markov, champs
de Markov), et quelques passages de ce document s’inspire largement de cours ou de rapports techniques,
certains disponibles sur Internet, dont voici les principaux :
? Olivier Cappé, Modèles de mélange et modèles de Markov cachés pour le traitement automatique de la
parole, juin 2000, http://www.tsi.enst.fr/~cappe/cours/tap.pdf.
? Wojciech Pieczynski et Alain Hillion, Bases mathématiques pour le traitement des images de télédétection,
Master Statistique, Université Paris VI, mars 2006.
? François Le Gland : Introduction au filtrage en temps discret. Filtrage de Kalman et Modèles de Markov cachés, Université de Rennes 1, Master Recherche Électronique, spécialité SISEA (Signal, Image,
Systèmes Embarqués, Automatique), http://www.irisa.fr/aspi/legland/rennes-1/
? Polycopié de cours de M. Sigelle and F. Tupin [ST99] intitulé Champs de Markov en traitement d’images,
perso.telecom-paristech.fr/~tupin/cours/polymrf.pdf.
D’autres part, un certain nombre d’articles et de livres qui traitent en détail les sujets abordés dans ce cours
peuvent être consultés :
? Article de L. R. Rabiner [Rab89] intitulé A tutorial on hidden Markov models and selected applications
in speech recognition, http://www.cs.ubc.ca/~murphyk/Bayes/rabiner.pdf.
? Article de W. Pieczynski [Pie03], intitulé Modèles de Markov en traitement d’images, http://www-public.
int-evry.fr/~pieczyn/A31.pdf.
? Rapport technique de J. Bilmes [Bil97] intitulé A gentle tutorial on the EM algorithm and its application
to parameter estimation for Gaussian mixture and hidden Markov models, http://www.icsi.berkeley.edu/
ftp/global/pub/techreports/1997/tr-97-021.pdf.
? Livre de O. Cappé, E. Moulines et T. Rydén [CMR05], intitulé Inference in hidden Markov models.
Mais bien sûr, des centaines d’autres références sont possibles. . . !

3

4

Table des matières
1 Introduction et motivations

8

2 Décision bayésienne

11

2.1

Position du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.2

Stratégie bayésienne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.3

Exemple : cas gaussien . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

3 Modèle de mélange - cas indépendant

17

3.1

Modèle de mélange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.2

L’algorithme EM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

3.2.1

Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

3.2.2

Propriétés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

EM et mélange gaussien . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

3.3.1

Quantité intermédiaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

3.3.2

Formules de ré-estimation des paramètres . . . . . . . . . . . . . . . . . . . . . . . . .

25

3.4

Simulations et exemples en segmentation d’images . . . . . . . . . . . . . . . . . . . . . . . .

27

3.5

Synoptique de l’algorithme complet de mélange aveugle . . . . . . . . . . . . . . . . . . . . .

28

3.3

4 Chaînes de Markov cachées
4.1

4.2

31

Le modèle de chaîne de Markov cachée . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

4.1.1

Loi de X a priori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

4.1.2

Loi de (X, Y ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

4.1.3

Probabilités « forward-backward » . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

4.1.4

Loi de X a posteriori

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

Décision bayésienne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

4.2.1

39

Critère du MPM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5

4.2.2

Critère du MAP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

L’algorithme EM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

4.3.1

Quantité intermédiaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

4.3.2

Ré-estimation des paramètres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

4.4

Simulations et exemples en segmentation d’images . . . . . . . . . . . . . . . . . . . . . . . .

44

4.5

Synoptique de l’algorithme des CMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

4.3

5 Champs de Markov cachés

49

Loi de X a priori : champs de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

5.1.1

Champs de Gibbs, champs de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

5.1.2

Echantillonneurs de Gibbs et de Metropolis . . . . . . . . . . . . . . . . . . . . . . . .

51

5.1.2.1

Echantillonneur de Gibbs . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

5.1.2.2

L’algorithme de Metropolis . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

Quelques MRF fondamentaux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

5.1.3.1

Modèle d’Ising . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

5.1.3.2

Modèle de Potts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

5.1.3.3

Modèle markovien gaussien . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

5.2

Loi a posteriori et loi conjointe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

5.3

Décision bayésienne . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

5.3.1

Critère du MPM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

5.3.2

Critère TPM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

58

5.3.3

Critère du MAP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

5.3.4

Algorithme ICM (approximation du MAP) . . . . . . . . . . . . . . . . . . . . . . . .

59

Estimation des paramètres du modèle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

5.4.1

Estimation avec échantillon d’apprentissage . . . . . . . . . . . . . . . . . . . . . . . .

60

5.4.2

Estimation sans échantillon d’apprentissage . . . . . . . . . . . . . . . . . . . . . . . .

61

Exemples en segmentation d’images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

5.5.1

Segmentation d’une image radar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

5.5.2

Comparaison de méthodes de segmentation . . . . . . . . . . . . . . . . . . . . . . . .

63

5.1

5.1.3

5.4

5.5

6 Filtrage de Kalman et extensions

67

6.1

Système linéaire gaussien . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

6.2

Filtre de Kalman . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

6

6.3

6.2.1

Prédiction/correction

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

6.2.2

Théorème de Kalman-Bucy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

6.2.3

Démonstration du théorème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

Filtre de Kalman étendu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

Bibliographie

73

7

Chapitre 1

Introduction et motivations
L’étude probabiliste des phénomènes s’introduit naturellement lorsqu’il existe une incertitude sur la mesure
décrivant un phénomène. Considérons deux phénomènes « mesurés » par deux réels x et y. Lorsqu’on cherche à
étudier des liens entre les deux phénomènes il existe, en dehors de la théorie des probabilités, deux possibilités :
soit un lien déterministe y = f (x), soit aucun lien. Le calcul des probabilités permet d’introduire une infinité
de « liens intermédiaires » : les deux phénomènes peuvent être plus au moins liés. Lorsqu’on « observe »
x, on dispose sur y d’une certaine information, sans pour autant pouvoir le calculer explicitement. Cette
information est modélisée par une « mesure de probabilité » notée P Y |x : pour tout A ⊂ R, P Y |x [A] est un
nombre dans [0, 1] donnant la « probabilité » pour que y soit dans A. y est ainsi considéré comme réalisation
d’une variable aléatoire Y et P Y |x est la loi de Y conditionnelle à x. Si on veut faire des raisonnements
généraux, valables pour tout x, on est amené à le considérer également comme une réalisation d’une variable
aléatoire X dont le comportement est décrit par la loi de probabilité PX . On arrive ainsi à la loi de probabilité


du couple (X, Y ), donnée par PX et la famille P Y |x , x ∈ R, modélisant les « liens stochastiques » entre
les deux phénomènes.
Ce type de modélisations peut être utile dans le traitement de certains problèmes se posant en imagerie.
D’une façon générale X contient l’information que l’on recherche mais n’est pas directement observable (on
dit généralement qu’elle est « cachée »). On observe, ou mesure, Y = y et on cherche à retrouver, ou à
« estimer », la réalisation cachée x. Considérons, à titre d’exemple, le problème suivant : on cherche à savoir,
à partir d’une image satellite, si un certain pixel de l’image représente de la forêt ou de l’eau. X prend
ainsi ses valeurs dans un ensemble de deux éléments Ω = {ω1 = "eau", ω2 = "forêt"} et l’observation Y = y
est donnée par un nombre représentant un niveau de gris. La loi de X, appelée « a priori », est donc une
probabilité sur Ω et modélise, de façon générale, la connaissance que l’on a sur le phénomène modélisé par
X « a priori », i.e. sans aucune mesure. Si on sait, dans notre exemple, que l’image a été prise dans une
région où il y a trois fois plus d’eau que de forêt, on posera PX (ω1 ) = 0.75 et PX (ω2 ) = 0.25. Les deux lois
conditionnelles P Y |ω1 et P Y |ω2 modélisent plusieurs phénomènes différents. Les classes « eau » et « forêt »
ne produisent pas une mesure unique (présence de vagues, « variabilité naturelle » de la forêt. . . ), d’où une
variation stochastique de la mesure Y = y. À cette variation peuvent s’ajouter divers « bruits », dus à la
transmission, l’acquisition. . .
L’étape suivante, après la définition de la loi du couple PX,Y modélisant les interactions stochastiques entre
les « mesures » décrivant les phénomènes, est la définition d’une règle de calcul de la réalisation cachée
X = x à partir de la réalisation observée Y = y. Une fois retenue, la règle, ou « stratégie », y = s(x)
est déterministe. Mais, contrairement au cas de lien déterministe entre y et x, on dispose généralement de
plusieurs choix possibles. Le choix est opéré à partir d’un « critère de qualité », ce dernier étant fonction
des résultats que l’on attend de s. Si on reprend l’exemple ci-dessus on peut considérer comme critère de
qualité la « probabilité de se tromper » qui peut, pour s donnée, être calculée à partir de PX,Y . La stratégie

8

(a)

(b)

Figure 1.1 – Deux exemples : (a) Signal illustrant l’impact de la nicotine sur le rythme cardiaque. (b) Image
satellitaire (SPOT) de l’île de la réunion.
s qui optimise ce critère est la stratégie bayésienne avec la fonction de perte « 0 − 1 ». Si pour une raison
quelconque on considère que les deux types d’erreurs « décider qu’il y a de l’eau alors qu’il y a de la forêt »
et « décider qu’il y a de la forêt alors qu’il y a de l’eau » ne sont pas de gravité égale, ce qui est fréquent
dans les problèmes de détection, on modifie la fonction de perte et on arrive à une stratégie s différente de
la précédente. Dans la pratique le calcul explicite optimisant un critère donné est parfois impossible, on est
alors amené à rechercher une stratégie s « sous-optimale ». Finalement, lorsque PX,Y est donnée, on choisit
un critère correspondant à la nature du problème que l’on veut résoudre et on cherche s optimisant ce critère.
En fait, PX,Y est rarement connue avec précision dans la pratique et on doit, dans une phase préalable
d’« estimation », rechercher des paramètres dont la connaissance est suffisante pour la détermination de s.
Le cas le plus fréquent est celui où on a une idée sur la forme générale de PX,Y et on considère qu’elle fait
partie d’une famille PX,Y,θ , θ ∈ Θ. On cherche alors à estimer θ. En reprenant notre exemple, supposons que
les lois de Y conditionnelles à ω1 et ω2 sont gaussiennes, notons f1 et f2 les densités correspondantes. Nous
sommes dans le cas « paramétrique », θ a six composantes :
? les lois « a priori » : π(k) = PX=ωk , k ∈ Ω ;
? les moyennes et écarts-type définissant f1 et f2 : θ k = {µk , σk }, k ∈ Ω.
Dans les cas « paramétrique » et « non paramétrique », on distingue deux sous-cas :
1. estimation avec « échantillon d’apprentissage » (supervisée) et
2. estimation « sans échantillon d’apprentissage » (non supervisée).
Le premier est celui où on dispose d’un échantillon x01 , x02 , . . . , x0M de réalisations de X. Les réalisations de
Y étant toujours observables, on estime alors les paramètres nécessaires à la détermination de s à partir
0
de (x01 , y10 ), . . . , (x0N , yN
). Dans notre exemple simple le cas « avec échantillon d’apprentissage » est celui
où on dispose dans l’image de M endroits où la nature du terrain ("eau" ou "forêt") est connue. Notons
x = {x1 , . . . , xM } l’échantillon « observé », x1 = {x1 , . . . , xP } et x2 = {x1 , . . . , xQ } les sous-échantillons
P
"eau" et "forêt" (P + Q = M ). Les lois a priori π(1) et π(2) peuvent être estimées par les fréquences M
et
Q
M , et les paramètres θ k , k ∈ Ω des gaussiennes par les moyennes et écarts-type empiriques calculés à partir
des deux sous-échantillons x1 et x2 .
Le deuxième est celui où les paramètres nécessaires à la détermination de s doivent être estimés à partir de
l’échantillon y de Y , dont la loi est un mélange de lois gaussiennes : on arrive au problème statistique général
de reconnaissance de mélange. Le cas le plus général, mais aussi le plus difficile à traiter, est celui de
9

l’estimation non paramétrique sans échantillon d’apprentissage.
Finalement la démarche générale, que l’on retrouve dans le traitement de nombreux problèmes (en traitement
de signal, économie, médecine, . . . ) est la suivante :
? on définit la forme générale de PX,Y ;
? on définit un critère de qualité de l’« estimation » de X = x (caché) à partir de Y = y (observé).
? on recherche une « stratégie » s optimisant le critère choisi.
? on estime les paramètres nécessaires à la mise en œuvre de s.
L’objectif du cours est de présenter l’utilisation des modèles de Markov, dans le cadre du traitement numérique du signal. Les modèles présentés sont très généraux, et peuvent s’appliquer directement dans d’autres
domaines comme l’économie et la finance, la génomique, la reconnaissance de formes. . . Dans ce cours, nous
ne nous intéresserons qu’aux modèles à temps discret (indice n discret) et à états discrets (X ∈ Ω ⊂ Z),
qui correspondent bien aux signaux numériques tels que ceux extraits de la parole [HCF+ 06], ou aux images
considérées comme des matrices de niveaux de gris. Ces deux domaines ont d’ailleurs largement contribué au
développement des modèles de chaîne de Markov cachée et de champs de Markov caché, que nous détaillons
plus particulièrement dans ce cours.
La démarche générale présentée ci-dessus conduit à plusieurs modélisations stochastiques qui différent selon
les hypothèses de dépendance que l’on considère pour modéliser les liens stochastiques entre les données
observées et les données cachées. Le cours commence par introduire le principe de la stratégie bayésienne
de la décision (chapitre 1). Puis, nous continuons en posant le problème de l’estimation d’un mélange fini
dans le cas aveugle, c’est-à-dire celui où les échantillons sont supposés indépendants les uns des autres
(variables i.i.d.). Nous détaillons l’algorithme itératif EM (Expectation-Maximization) et les formules de
ré-estimation dans le cas gaussien (chapitre 2). Cette hypothèse d’indépendance pouvant apparaître trop
restrictive dans bon nombre d’applications (séries chronologiques, images. . . ), nous poursuivons avec deux
modélisations qui permettent de considérer des liens plus riches entre les données : les chaînes de Markov
(1D) et champs de Markov (2D), présentés dans les chapitres 3 et 4 respectivement. Pour ces deux
modèles, nous détaillons les algorithmes d’estimation EM associés (cas gaussien), permettant d’aboutir à des
méthodes de classification non supervisées. Les algorithmes sont illustrés dans le contexte de la segmentation
d’images, en particulier sur les images radar pour lesquelles ces méthodes markoviennes s’avèrent robustes
et particulièrement performantes.

10

Chapitre 2

Décision bayésienne
Nous présentons brièvement dans ce chapitre le principe de la classification bayésienne, qui est traitée dans
bon nombre d’ouvrages classiques. Une présentation orientée « reconnaissance de formes » peut être consultée
dans [DHS01] (chapitres 2 et 3 essentiellement) et dans [TK06] (chapitre 2 essentiellement). Les démarches
considérées ici sont très générales et s’appliquent à un grand nombre de problèmes en dehors de la classification des signaux et de la segmentation statistique d’images.

2.1

Position du problème

Selon le schéma général, on observe une réalisation d’une variable aléatoire Y (une valeur y de R) et on
souhaite « estimer » la valeur cachée k du paramètre. L’ensemble des paramètres Ω sera supposé fini, Ω =
{1, . . . , K}, ses éléments appelés « classes » et tout estimateur « stratégie de classification » 1 .
Supposons maintenant que nous nous trouvons devant un problème de classification des données et que nous
connaissons la fréquence d’apparition des classes. Par exemple, on classe les individus en classe « homme »
et « femme », uniquement à partir de leur poids. On sait a priori (ce qui signifie ici « avant l’observation »)
que la population que nous devons classer contient deux tiers d’hommes et un tiers de femmes. Une telle
connaissance a priori peut être modélisée par une probabilité (dite a priori) sur Ω. Cette probabilité peut
alors être considérée comme la loi d’une variable aléatoire X et les p (. |X = k ) apparaissent comme les lois
de Y conditionnelles à X.
Finalement, la loi a priori p (X = k) = p (k) = π(k) sur Ω et les lois conditionnelles p (y |X = k ) = fk (y)
sur R définissent une probabilité p (y, k) = π(k) fk (y) sur R × Ω 2 , dite loi du couple ou loi conjointe. La loi
p (y) est appelée densité mélange ou, plus simplement, mélange
p (y) =

K
X

p (y, k) =

k=1

K
X

π(k) fk (y).

(2.1)

k=1

La probabilité conditionnelle p (k |y ) (sachant y ∈ R) sur Ω, dite loi a posteriori s’écrit
p (k |y ) =

p (y, k)
π(k) fk (y)
= K
.
p (y)
X
π(i) fi (y)

(2.2)

i=1

1. Par abus de notation, et lorsque cela ne peut engendrer de confusion, nous écrirons p (x) à la place de p (X = x), x ∈ Ω
et p (y) à la place de p (Y ∈ dy), y ∈ R.
2. par rapport à la mesure δ ⊗ ν, où δ est la mesure de comptage et ν la mesure de Lebesgue sur R.

11

Intuitivement, la différence entre la probabilité a priori p (k) et la probabilité a posteriori p (k |y ) sur Ω
illustre l’apport de l’information (sur l’identité de la classe non observable) contenue dans l’observation (a
priori signifie « avant » l’observation, et a posteriori signifie « après » l’observation). On retrouve le fait que
si les variables sont indépendantes, l’observation de l’une d’entre elles n’apporte aucune connaissance sur le
comportement de l’autre et donc ces deux probabilités sont égales.

2.2

Stratégie bayésienne

Considérons une probabilité sur Ω × R qui est une loi d’un couple de variables aléatoires (X, Y ). Ainsi
(x, y) ∈ Ω × R étant une réalisation de (X, Y ), le problème de la classification devient celui de l’estimation
de la réalisation inobservable de la variable X à partir de la variable observable Y .
Considérons une stratégie de classification sˆ : R −→ Ω. Pour chaque réalisation (x, y) = (X, Y ), sˆ peut
donner la bonne réponse, i.e. sˆ(y) = x , ou se tromper, i.e. sˆ(y) 6= x. Supposons que les différentes erreurs
ne sont pas de gravité équivalente. On le modélise en définissant une application L : Ω × Ω −→ R+ dite
fonction de perte :
(
0
si i = j,
L(i, j) =
(2.3)
λi,j sinon.
Le nombre réel λi,j modélise la gravité de l’erreur « on a choisi la classe i alors que la vraie classe est
j ». Insistons sur le fait que la « perte » modélisée par L ne fait pas partie de la modélisation probabiliste
considérée. Par ailleurs, à une erreur donnée, deux utilisateurs peuvent avoir des intérêts différents, et donc
les pertes qu’ils associent à une même erreur peuvent être différentes. Ainsi que nous allons le voir dans la
suite, la possibilité de l’utilisation des fonctions de perte différentes introduit une grande généralité - et une
grande souplesse - des modèles probabilistes utilisés à des fins de classification des données.
À stratégie sˆ et fonction de perte L données, comment mesurer la qualité de sˆ ? Supposons que l’on a N
observations indépendantes y = {y1 , . . . , yN }, chacune correspondant à une classe inconnue, à classer. En
notant x = {x1 , . . . , xN } les classes correspondantes, la perte globale est
L (ˆ
s(y1 ), x1 ) + . . . + L (ˆ
s(yN ), xN ) .
On cherche à minimiser cette perte globale, ce qui revient à minimiser son quotient par N . Par la loi des
grands nombres, ce dernier tend vers :
L (ˆ
s(y1 ), x1 ) + . . . + L (ˆ
s(yN ), xN )
−→ E [L (ˆ
s(Y ), X)] .
N →+∞
N
On constate qu’à « long terme », la qualité d’une stratégie sˆ est mesurée par E [L (ˆ
s(Y ), X)], qui est appelée
« perte moyenne ». La stratégie bayésienne sˆB est celle parmi toutes les stratégies pour laquelle la perte
moyenne est minimale :
E [L (ˆ
sB (Y ), X)] = min E [L (ˆ
s(Y ), X)].
(2.4)


La qualité de sˆB est ainsi appréhendée via la loi des grands nombres et on ne peut rien dire pour une seule
observation (ou même un petit nombre).
Montrons que la stratégie bayésienne associée à la fonction de perte définie par l’équation (2.3) est
"
#
K
K
X
X

sB (y) = k] ⇐⇒ ∀j ∈ Ω,
λk,i p (i |y ) ≤
λj,i p (i |y ) ,
i=1

i=1

soit encore
sˆB (y) = arg min
j∈Ω

K
X
i=1

12

λj,i p (i |y ).

(2.5)

Démonstration : En appliquant la formule de Fubini 3 à E [L (ˆ
s(Y ), X)], on peut écrire :




E [L (ˆ
s(Y ), X)] = E E [ L (ˆ
s(Y ), X)| Y ] .
|
{z
}
φ(y)

4

Nous obtenons :
φ(y) =

K
X

L (ˆ
s(y), i) p (i |y ).

i=1

PK
L’élément sˆ(y) = k, qui minimise φ(y), minimise la quantité i=1 λj,i p (i |y ), ce qui donne (2.5). Notons en
R
effet que sˆB ainsi déterminée minimise bien E [L (ˆ
s(Y ), X)] car on a E [φ(Y )] = R φ(y) p (y) dy, et donc la
minimisation de φ en tout point minimise bien l’intégrale (car p(y) > 0).
Remarque 1 : Pour calculer la perte moyenne ξ (qui est minimale pour la stratégie bayésienne) associée
à la stratégie sˆ et à la fonction de perte L, on utilise toujours le résultat de la note 3 (en conditionnant
par X) et celui de la note 4 :

Z
φ(y) p (y) dy =

ξ = E [L (ˆ
s(Y ), X)] =

Z X
K

R

π(i)fi (y) L (ˆ
s(y), i) dy.

R i=1

Nous disposons ainsi de la stratégie qui assure, à long terme, d’avoir une perte minimale et, de plus, il est
possible de calculer sa valeur (cf. exemple ci-après).
Remarque 2 : Ainsi la stratégie bayésienne dépend des λi,j que l’on choisit de façon subjective. Si on
souhaite détecter une classe donnée avec une précision , on peut calculer les coefficients λi,j de façon
à ce que la stratégie bayésienne correspondante vérifie cette condition. Ce type de possibilités montre la
puissance de la modélisation en question.

Exemple : Soit Ω = {1, . . . , K} et la fonction de perte L0−1 définie par :
(
0 si i = j
L0−1 (i, j) =
1 sinon

(2.6)

L0−1 (ˆ
s(y), k) désigne alors la valeur, au point (k, y), de la fonction indicatrice du sous-ensemble de Ω × R
sur lequel sˆ se trompe et donc E [L0−1 (ˆ
s(Y ), X)] représente la probabilité pour que sˆ se trompe. Ainsi dans
ce cas la stratégie bayésienne sˆB définie par
sˆB (y) = k si ∀j ∈ Ω,

p (k |y ) ≥ p (j |y ) ,

(2.7)

qui est un cas particulier de (2.5), est celle pour laquelle la probabilité de se tromper est minimale 5 . Sachant
qu’en vertu de la loi des grands nombres la probabilité d’un événement peut être vue comme la fréquence de
son apparition lorsque le phénomène se reproduit un grand nombre de fois de façon indépendante, la stratégie
définie ci-dessus est celle qui produira, lorsqu’on l’utilisera dans un grand nombre de cas indépendants, la
plus petite proportion d’erreurs.
Ainsi sˆB (y) consiste, dans ce cas, à associer à chaque y ∈ R l’élément de Ω dont la probabilité a posteriori,
i.e. conditionnelle à Y = y, est maximale. Cette règle de décision est aussi appelée celle du « maximum de
vraisemblance a posteriori ». Notons que les probabilités a posteriori de (2.7) peuvent être remplacées par
les « fonctions discriminantes » π(j)fj (y), et la stratégie sˆB (y) s’écrire
sˆB (y) = k si ∀j ∈ Ω,

π(k)fk (y) ≥ π(j)fj (y)

(2.8)

3. Pour deux variables aléatoires réelles U , V et une fonction quelconque Ψ, la formule de Fubini est E [Ψ(U, V )] =
E [E [ Ψ(U, V )| U ]] = E [E [ Ψ(U, V )| V ]].
Z
h(u) p (u |v ) du.

4. En utilisant la version discrète du résultat classique suivant : E [ h(U )| V = v] =
R

5. On le démontre en exprimant l’eq. (2.5) dans ce cas particulier de 2 classes : nous sélectionnons la classe 1 si λ1,1 p (1 |y ) +
λ1,2 p (2 |y ) = p (2 |y ) est plus petit que λ2,1 p (1 |y ) + λ2,2 p (2 |y ) = p (1 |y ), d’où le résultat.

13

Figure 2.1 – Dessin de deux densités gaussiennes de paramètres θ1 = {100, 6} et θ2 = {110, 3}.
Cette dernière écriture est intéressante pour son interprétation graphique (cf. question 1 de la section 2.3).
La perte minimale définie dans la remarque précédente s’écrit dans le cas de la fonction de perte L0−1
Z
p (y) − max π(i)fi (y) dy.

ξ=

i∈Ω

R

(2.9)

ou, plus simplement encore dans le cas de 2 classes :
Z
ξ=

min π(i)fi (y) dy.

R i∈Ω

(2.10)

Ce résultat sera interprété dans la question 2 de la section 2.3.
Remarque 3 : Pour faciliter la lecture, nous avons pour l’instant considéré le cas scalaire, c’est-à-dire y ∈ R
(D = 1). Les résultats énoncés s’étendent sans difficulté au cas vectoriel où les observations sont vectorielles,
c’est-à-dire y ∈ RD (D > 1).

2.3

Exemple : cas gaussien

Dans le contexte gaussien considéré ici, notons µk et Σk les vecteurs moyens et les matrices de covariance
n
o
D × D correspondantes, θk = µk , Σk . On rappelle que fk (.) est décrite par l’expression analytique
suivante :



1
1
− 21
−1
t
fk y =
exp − (y − µk ) Σk (y − µk ) .
(2.11)
D |Σk |
2
(2π) 2
Lorsque D = 1, nous retrouvons bien évidemment la densité d’une gaussienne mono-dimensionnelle θk =
{µk , σk } :
"

2 #
1
1 y − µk

fk (y) =
exp −
.
(2.12)
2
σk
σk 2π
À titre d’exemple, la figure 2.1 montre deux gaussiennes qui serviront pour les exercices suivants.
Questions 1. Considérons le cas de deux gaussiennes de paramètres θ1 = {0, σ} et θ2 = {a, σ/2} (a réel) et
de proportions π(1) = 1/3 et π(2) = 2/3. Calculer de manière analytique les seuils de décision bayésienne,
c’est-à-dire les Υ valeurs {τ1 , . . . , τΥ } qui séparent les deux classes sur R. Calculer les valeurs numériques
pour les lois dont les valeurs des paramètres sont données dans la figure 2.1.
14

(a)

(b)

Figure 2.2 – (a) Mélange des deux lois gaussiennes de la figure 2.1 dans des proportions données par
π(1) = 1/3 et π(2) = 2/3. (b) En couleur cyan (resp. magenta) apparaît la courbe de « π(k) fk (.) »
maximum (resp. minimum).

Réponse.



1
1

exp −
2
3 σ 2π
1
exp
2

2
y
σ

"

=

1
4

exp −
2
3 σ 2π

=

4

3y 2 − 8ay + 4a2

=

4σ 2 ln 2

3y 2 − 8ay + 4a2 − 4σ 2 ln 2

=

0.

"

2(y − a)
σ

2


2
y
σ



y−a

2 #

σ
2

#



Le discriminant ∆ = 64a2 − 12 (4a2 − 4σ 2 ln 2) = 16 a2 + 3σ 2 ln 2 est toujours positif et les racines réelles
(c’est à dire les seuils de décision) sont données par τ1 =


8a− ∆
6

et τ2 =


8a+ ∆
.
6

En utilisant les valeurs numériques de la figure 2.1, nous obtenons ∆ = 2797.8, τ1 = 104.5 et τ2 = 122.1,
ce que l’on peut vérifier sur le graphe (b) de la figure 2.2. Ainsi, les valeurs de y comprises dans ]τ1 ; τ2 [
sont associées à la classe 2, alors que les valeurs en dehors de cet intervalle sont associées à la classe 1.
Nous avons autant de chance de commettre une erreur ou de ne pas en commettre en décidant, au niveau
des seuils (y = τ1 et y = τ2 ), de classer la donnée dans l’une ou l’autre classe.
La graphe (a) de la figure 2.2 représente le mélange des deux classes, cf. eq. (2.1). La courbe cyan du
graphe (b) représente, en chaque valeur y, la plus forte valeur π(k)fk (y) pour k ∈ Ω, alors que la courbe
magenta représente le minimum.

Questions 2. Dans le cas de la fonction de perte L0−1 définie par l’éq. (2.6), calculer de manière analytique
la perte moyenne ξ, cf. eq. (2.9), en utilisant les paramètres de la question 1. Calculer les valeurs numériques
pour les valeurs des paramètres données dans la figure 2.1.
Réponse. La courbe magenta modélise la probabilité de se tromper et la courbe la cyan la probabilité de
ne pas commettre d’erreur. La perte moyenne est donc représentée par la surface sous la courbe magenta,
et est calculée en intégrant cette fonction.

15

Calcul :

τ1

Z

Z

τ2

Z

π(2)f2 (y) dy +

ξ=
−∞

π(2)f2 (y) dy .

τ1

|

{z

}

A

+∞

π(1)f1 (y) dy +
τ2

|

{z

}

B

|

{z
C

}

Nous obtenons pour le terme A :
A

Z

τ1





y−a 2
dy
2
σ
"

2 #
Z τ1
2 (y − a)
2
2
√ √
exp −
dy.
σ
3σ 2 π −∞
4

3 σ 2π

=

=

1
exp −
2
−∞





2
2
En posant z =
(y − a) (ainsi dz =
dy), nous avons
σ
σ
1 2
A= √
3 π
2
En notant erf(x) = √
π

Z

Z

τ1 −a


exp −z 2 dz.





−∞

x

exp −z 2 dz avec limx→∞ erf (x) = 1, nous avons





0

A=

1
3





2
(τ1 − a)
σ

1 + erf


.

Par un calcul similaire, nous obtenons pour B et C :
B

=

C

=









τ2
τ1


− erf
σ 2
σ 2



2
1
1 − erf
(τ2 − a)
.
3
σ
1
6



erf

,

En utilisant les valeurs numériques, nous obtenons A = 0.023, B = 0.075 et C = 1.71 10−5 , ce qui donne
finalement une perte moyenne de ξ = 0.098.

Questions 3. Calculer les seuils de décision bayésienne dans le cas suivant : µ2 = µ1 et σ2 = 2σ1 . De même
avec deux lois exponentielles de paramètres λ1 et λ2 (on rappelle que la densité d’une loi exponentielle s’écrit,
pour y > 0, f (y) = λ exp [−λy]).
Questions 4. Si on remplace la fonction de perte L0−1 en



0 si i = j
L0 (i, j) =

1

2


1

si i = 1 et j = 2

(2.13)

si i = 2 et j = 1

qu’obtenons-nous (les erreurs n’ont pas le même poids) ? Même question pour L1−0 (on favorise les erreurs).
Comment pondérer les pertes entre les classes de manière à avoir le même pourcentage d’erreur pour chacune
des deux classes (i.e. α1 = α2 dans la question 2) ?
Questions 5. Exercice 3 de l’examen de février 2010.

16

Chapitre 3

Modèle de mélange - cas indépendant
n
o
On considère ici un signal numérique se présentant sous forme d’un vecteur y = y 1 , . . . , y n , . . . , y N de N
t
vecteurs de paramètres y n = yn1 , . . . , ynD de dimension D, observés à intervalle régulier. Pour exemple, on
peut citer
? les pixels d’une image à niveaux de gris (D = 1) ou multi-spectrale (D = 3, 4 ou 5),
? les températures et pressions dans le cadre d’applications météorologiques (D = 2),
? les coefficients cepstraux, dans le cadre du traitement de la parole (D = 10 ou 20), . . .
On trouve sur Internet des démonstrateurs sur les modèles de mélange et l’algorithme EM, notamment http:
//wiki.stat.ucla.edu/socr/index.php/SOCR_EduMaterials_ModelerActivities_MixtureModel_1.

3.1

Modèle de mélange

Le modèle de mélange probabiliste (probabilistic mixture model) [TSM85] consiste à supposer que les
vecteurs observés y n sont des réalisations de variables aléatoires Y n mutuellement indépendantes 1 , qui
suivent toutes une même loi f (.) ayant la forme suivante :

K
X

f yn =
π(k) fk y n .

(3.1)

k=1

PK
Chaque fk (.) est une densité de probabilité et les π(k) sont des scalaires positifs tels que k=1 π(k) = 1.
Nous allons voir que ce modèle peut être interprété en supposant que les données observées sont réparties
dans K classes (appelées aussi composantes du mélange).

1. Soient A1 , . . . , Am , m événements. On dit qu’ils sont mutuellement indépendants si pour toute famille J d’indices dans
[1, . . . , m], nous avons :

!
P

\

Aj

=

j∈J

Y
j∈J

Ils sont en particulier deux à deux indépendants (réciproque fausse).

17

P (Aj ).

(a)

(b)

(c)

(d)

Figure 3.1 – Échantillon (a) et indicatrices (b) de N = 20 valeurs générées à l’aide d’un mélange de K = 3
composantes gaussiennes : (c) et (d). Les proportions théoriques (ayant servies à la simulation) sont données
par π(1) = 0.1 (rouge), π(2) = 0.55 (vert) et π(3) = 0.35 (bleu). Les paramètres θk = {µk , σk } des 3
gaussiennes sont donnés par (−2, 4), (−0.8, 1) et (1, 2).

Remarque 1 : On parle de modèle de mélange paramétrique lorsque les densités fk (.) des composantes
sont issues de familles paramétriques telles que les lois gaussienne, gamma, . . . On notera θk les paramètres
de fk (.) et Θ = {θ1 , . . . , θK }.
À titre d’illustration la courbe (d) de la figure 3.1 montre le résultat du mélange de trois lois gaussiennes
1D. Un exemple de mélange de trois lois gaussiennes 2D est illustré dans la figure 3.2. Les vecteurs moyens,
cf. (2.11) page 14, sont donnés par
µ1 = (−3, −1)t ,

µ2 = (0, 0)t ,

µ3 = (2, 0)t ,

et les matrices de covariance par













1
0.5
2
−0.8
1
0.1
, Σ2 =
, Σ3 =
0.5
2
−0.8
2
0.1
3
On distingue les trois formes ellipsoïdales (qui n’ont d’ailleurs pas la même orientation, ce qui traduit le
fait que les matrices de covariances sont différentes pour chaque composante du mélange). Il est clair que
selon la manière dont les vecteurs moyens des composantes diffèrent, l’individualisation des composantes
sera plus où moins marquée.
Σ1 =

18

Figure 3.2 – Représentation des courbes isoprobabilité d’un mélange de trois lois gaussiennes 2D selon les
mêmes proportions que dans la figure 3.1.

Remarque 2 : Comment réalise t’on un tirage aléatoire selon un mélange pondéré de densités ? Cela se
fait en deux étapes
1. D’abord on tire un numéro k de classe selon les probabilités {π(i)}i∈Ω .
2. Ensuite, on tire selon la densité sélectionnée fk (.).
Pour illustrer ce résultat, nous avons conduit l’expérience suivante, cf. figure 3.3. Nous avons ajouté un
bruit gaussien à chacune des classes de l’image de cible (a) dont la classe noire représente environ 2/3 des
pixels de l’image, et la classe blanche 1/3. Nous avons procédé de la manière suivante. Nous avons bruité
chacun des pixels noirs avec la loi gaussienne f2 de paramètres {110, 3}. Les pixels de la classe blanche ont
quand à eux été bruités avec la loi f1 de paramètres {100, 6}. L’image bruitée et son histogramme sont
présentés dans les figures (b) et (c). Nous avons alors appliqué la décision bayésienne avec l’ensemble des
paramètres de la simulation, et obtenu l’image classée (d). L’erreur de classification obtenue en comptant
le nombre de pixels différents entre (a) et (d) s’élève à 9.76%. Ce résultat est une bonne approximation
du calcul théorique (rappel : ξ = 0.098), ce qui est conforme au principe de la loi des grands nombres car
l’image a pour dimensions 128 × 128 = 16384 pixels.
Il est remarquable de constater que le pourcentage d’erreur dans la classe noire est de α1 = 3.36% (ce qui
correspond bien à une erreur de 2.26% par rapport au total des pixels, à rapprocher de A + C), alors que
celui dans la classe blanche est de α2 = 22.60% (ce qui correspond bien à une erreur de 7.40% par rapport
au total des pixels, à rapprocher de B). Ainsi les erreurs commises entres les deux classes ne sont pas les
mêmes : une classe est bien plus erronée que la seconde.
À titre d’illustration, la figure 3.4 montre le résultat d’un tirage de 4000 échantillons issus du mélange
2D décrit dans la remarque ci-dessus. Le poids de chaque composante du mélange π(k) se traduit par la
proportion statistique de vecteurs venant de chacune des composantes (par application de la loi des grands
nombres).
La figure 3.5 montre le résultat d’un tirage de 128 × 128 échantillons avec les mêmes paramètres que
ceux utilisés pour obtenir l’image 3.3(b). Ces deux images sont donc issues du même mélange et ont
approximativement le même histogramme. Par contre, elles présentent des aspects très différents. . .

Comme le laisse présager ce qui précède, la variable indicatrice Xn est une donnée constitutive du problème
qui présente l’inconvénient de ne pouvoir être observée en pratique : on observe des réalisations du vecteur
aléatoire Y n sans savoir de manière certaine quelle est la classe du mélange associée à chaque observation.
Au sens de l’algorithme EM, la variable Xn constitue une donnée latente, c’est-à-dire fortement suggérée
par le problème considéré (on parle également de donnée non-observée ou manquante). Nous verrons que
l’introduction de ces données non-observées permet de résoudre de manière élégante un problème d’estimation
relativement complexe.
Ce modèle est décrit par le terme « aveugle » dans la mesure où il ne prend pas en compte l’influence des
autres échantillons. Cette hypothèse est acceptable dans certaines applications, mais la prise en compte de
19

(a)

(b)

(c)

(d)

Figure 3.3 – Décision bayésienne sur une image bruitée avec paramètres connus. (a) image originale (π(2) '
0.33 et π(1) ' 0.67). (b) image bruitée avec les paramètres donnés dans le texte et (c) son histogramme
normalisé, à rapprocher du mélange des deux lois gaussiennes de la figure 2.2(a). (d) image classée obtenue
par décision bayésienne.

20

Figure 3.4 – Tirage aléatoire de 4000 échantillons issus du mélange de trois gaussiennes 2D présenté dans
la remarque ci-dessus. La couleur de chaque point montre la classe k à laquelle l’échantillon appartient.

Figure 3.5 – Image construite par tirage selon le même mélange que celui utilisé pour obtenir la figure 3.3(b).
l’aspect séquentiel (c’est à dire l’ordre dans lequel les données sont lues, comme les séries chronologiques ou
spatiales) s’avère parfois primordial. Les modèles markoviens que nous aborderons dans les chapitres 4 et 5
permettront d’infléchir cette hypothèse.

3.2

L’algorithme EM

L’algorithme « Expectation-Maximization (EM) » [DLR77, MK96, Bil97] est une méthode générale pour
trouver l’estimée du maximum de vraisemblance d’un ensemble de paramètres (noté Θ) d’une distribution
donnée à partir d’un échantillon. On sait en effet que cette stratégie d’estimation conduit à des estimateurs
asymptotiquement efficace, c’est-à-dire « optimaux » lorsque le nombre de données observées devient important. Le problème est que la fonction de vraisemblance L prend ici la forme relativement complexe suivante :
N X
K


Y
L Θ|y = p y|Θ =
π(k) fk y n ,

(3.2)

n=1 k=1

n
o
où y = y 1 , . . . , y n , . . . , y N désigne l’ensemble des N vecteurs dont on dispose pour estimer les paramètres
du modèle, et Θ regroupe les paramètres du modèle à estimer, c’est-à-dire les proportions π(k), les moyennes
µk et les matrices de covariance Σk pour k ∈ Ω = {1, . . . , K} d’un mélange gaussien. Dans cette équation,
les densités fk (.) doivent être remplacées par leur expression donnée à l’équation (2.11).
On conçoit aisément qu’il y ait quelques difficultés à maximiser la fonction de vraisemblance définie par (3.2)
21

par rapport aux paramètres Θ du modèle. Outre une expression analytique plutôt complexe, cette fonction
de vraisemblance présente le défaut de ne pas être convexe, c’est-à-dire qu’elle n’admet a priori pas de maximum unique [Min83]. Par conséquent même l’utilisation d’algorithmes d’optimisation classiques (gradient,
Newton) est ici délicate car il faudrait auparavant cerner le domaine dans lequel on désire rechercher la
valeur optimale de Θ, chose qui est plutôt malaisée lorsque Θ est un paramètre de dimension très élevée
comportant des données non-homogènes (par exemple les poids π(k) et les matrices de covariance Σk ).

3.2.1

Principe

L’algorithme EM apporte une solution extrêmement générale à ce type de problèmes pour lesquels le
modèle statistique considéré peut être complété en faisant appel à des données latentes. Dans le cas du
modèle de mélange nous avons vu que la variable non observée, qu’il est judicieux de faire intervenir, est
l’indicatrice de la composante du mélange xn associée à chaque vecteur observé y n . Si ces données latentes
pouvaient être observées, la solution du problème serait beaucoup plus simple.
L’idée à la base de l’algorithme EM consiste à raisonner sur les données complètes (données observées et
données latentes) tout en prenant en compte le fait que l’information disponible sur les données latentes
ne peut venir que des données observées. Le principe de l’algorithme EM est de compenser les données
manquantes en les remplaçant par leur moyenne. Ceci se traduit par l’algorithme itératif suivant :
1. estimer les données manquantes étant donné la valeur courante des paramètres,
2. estimer les nouveaux paramètres étant donné les estimées des données manquantes courantes,
3. réitérer les étapes 1. et 2. jusqu’à convergence.
Plus précisément, chaque itération ` se décompose en deux étapes [DLR77] :


ë Expectation Calcul de la fonction auxiliaire Q Θ|Θ(`)


h
i

Q Θ|Θ(`) = E ln p y, x |Θ y, Θ(`) .

?
?
?
?

(3.3)

Θ désigne les vrais paramètres,
Θ(`) désigne la valeur estimée des paramètres du modèle à la `e itération de l’algorithme,
y et x désignent respectivement l’ensemble des données observées et des données latentes associées,


p y, x |Θ = L Θ|x, y désigne la vraisemblance conjointe des données observées et latentes.

ë Maximisation de la fonction auxiliaire


Θ(`+1) = arg max Q Θ|Θ(`) ,
Θ

(3.4)

permettant d’obtenir une nouvelle estimée Θ(`+1) .

3.2.2

Propriétés

Propriété 1 : La
et la plus importante propriété de l’algorithme EM est le fait que la suite des
n première
o
(`)
valeurs estimées Θ
est construite de façon à ce que la vraisemblance des données observées augmente
22

à chaque itération de l’algorithme. En effet


X





= ln p y Θ(`)
ln p y Θ(`)
p x y, Θ(`)
x

|

{z
}





p x, y Θ(`)
X

=
p x y, Θ(`) ln
(`)
p
x
Θ
y,
x



X



X



p x y, Θ(`) ln p x, y Θ(`) −
p x y, Θ(`) ln p x y, Θ(`)
=
=1

x

x



i

i
h





= E ln p x, y Θ(`) y, Θ(`) − E ln p x y, Θ(`) y, Θ(`)
|
{z
} |
{z
}
h



Q(Θ(`) |Θ(`) )

(3.5)

R(Θ(`) |Θ(`) )

et
ln p y |Θ



=

ln

X

p x, y |Θ



x




(`)
p x y, Θ


ln 
p x, y |Θ

p x y, Θ(`)
x

 

p x, y |Θ
y, Θ(`) 
ln E 
p x y, Θ(`)


=

=

X

Or d’après l’inégalité de Jensen 2 , sachant que la fonction ln est concave et non convexe




p x, y |Θ
y, Θ(`) 
ln p y |Θ
≥ E  ln
(`)
p x y, Θ
h
h

i

i



≥ E ln p x, y |Θ y, Θ(`) − E ln p x y, Θ(`) y, Θ(`)
{z
} |
{z
}
|
Q(Θ|Θ(`) )

(3.6)

R(Θ(`) |Θ(`) )



Si Θ(`+1) = arg maxΘ Q Θ|Θ(`) , alors Q(Θ(`+1) |Θ(`) ) ≥ Q(Θ(`) |Θ(`) ), et donc depuis (3.5) et (3.6)






ln p y Θ(`+1) − ln p y Θ(`)



Q(Θ(`+1) |Θ(`) ) − R(Θ(`) |Θ(`) ) − Q(Θ(`) |Θ(`) ) + R(Θ(`) |Θ(`) )



Q(Θ(`+1) |Θ(`) ) − Q(Θ(`) |Θ(`) ) ≥ 0

Ainsi la vraisemblance croit toujours, ce qui constitue une condition suffisante pour assurer la convergence
de EM. La première vertu de l’algorithme EM est donc de permettre de construire une suite d’estimateurs
des paramètres du modèle pour laquelle la vraisemblance croît.
Bien sûr, ce qui rend l’algorithme EM intéressant en pratique est le fait que les deux équations (3.3) et (3.4)
vont avoir une forme analytique explicite, et ce pour une classe très large de modèles statistiques. Ce point
est détaillé pour les modèles de mélange au paragraphe 3.3, puis pour les modèles de Markov cachés conditionnellement gaussiens aux chapitres 4 et 5.
Propriété 2 : La seconde propriété importante de l’algorithme EM et que si l’on différencie la quantité
2. qui dit que pour une fonction φ convexe, i.e. ∀t ∈ [0, 1], φ(tx+(1−t)y) ≤ tφ(x)+(1−t)φ(y), nous avons φ (E [X]) ≤ E [φ(X)]

23

intermédiaire par rapport au paramètre Θ, il vient


h
i



(`)
(`)


∂Q Θ|Θ
∂E ln p x y, Θ y, Θ


∂ ln p y |Θ


=
+




∂Θ
∂Θ
∂Θ

(`)
Θ=Θ
|
{z
}
Θ=Θ(`)
A

Intéressons nous au terme terme A :
h
i


∂E ln p x y, Θ y, Θ(`)
∂Θ

Or, lorsque Θ = Θ(`) , A peut s’écrire :
"

#
X ∂p x y, Θ
A|Θ=Θ(`) =


∂Θ
x

.

(3.7)

Θ=Θ(`)






∂ X
=
ln p x y, Θ p x y, Θ(`) 
∂Θ x
i
X ∂ h


=
ln p x y, Θ p x y, Θ(`)
∂Θ
x




(`)
X p x y, Θ
∂p x y, Θ

.


=
∂Θ
p x y, Θ
x

Θ=Θ(`)






∂ X

=
p x y, Θ 
∂Θ x


Θ=Θ(`)


∂1
=
= 0.
∂Θ Θ=Θ(`)

L’algorithme EM possède donc une seconde vertu, en ce qu’il permet de calculer le gradient de la fonction
d’objectif (la vraisemblance) aux points Θ(`) . Cette propriété est très importante pour la convergence de
l’algorithme puisque qu’elle implique que les points stables de l’algorithme (c’est-à-dire des points tels que
Q (Θ|Θ∗ ) soit maximum en Θ∗ ) sont des points stationnaires de la vraisemblance, c’est à dire pour lesquels
∂ ln p(y|Θ )

= 0.
∂Θ


Θ=Θ

Avec quelques hypothèses concernant la régularité de la fonction intermédiaire Q(.) et la structure topologique
de l’espace des paramètres Θ, l’éq. (3.7) permet de montrer que l’algorithme EM ne peut converger que
vers des points stationnaires de la vraisemblance. En pratique, la vraisemblance n’étant pas une fonction
convexe des paramètres Θ, c’est-à-dire pouvant présenter plusieurs maximums locaux, le comportement de
l’algorithme EM dépend fortement de la valeur initiale Θ(0) depuis laquelle on le fait démarrer.
Variantes de EM : L’algorithme EM, bien que très performant et souvent simple à mettre en œuvre, pose
quand même parfois quelques problèmes qui ont donné lieu à des développements complémentaires. On
distingue trois variantes qui permettent de palier (au moins partiellement) à certaines difficultés connues
de EM :
ý GEM ( Generalized EM) GEM a été proposé en même temps qu’EM par Dempster et al. [DLR77]
qui ont prouvé que pour assurer la convergence vers un maximum local de vraisemblance, il n’est pas
nécessaire de maximiser Q à chaque étape mais qu’une simple amélioration de Q est suffisante.
ý CEM ( Classification EM) L’algorithme EM se positionne dans une optique estimation, c’est-à-dire
qu’on cherche à maximiser la vraisemblance du paramètre, sans considération de la classification faite a
posteriori en utilisant la règle de Bayes. L’approche classification, proposée par Celeux et Govaert [CG91]
consiste à optimiser, non pas la vraisemblance du paramètre, mais directement la vraisemblance complétée.
ý SEM ( Stochastic EM) Afin de réduire le risque de tomber dans un maximum local de vraisemblance,
Celeux et Diebolt [CD85] proposent d’intercaler une étape stochastique de classification entre les étapes E
et M.

3.3

EM et mélange gaussien



Nous commençons par expliciter le calcul de la quantité Q Θ|Θ(`) (valable pour n’importe quel type de
loi), puis nous donnons les formules de ré-estimation des paramètres dans le cas d’un mélange gaussien.
24

3.3.1

Quantité intermédiaire

Pour le modèle de mélange, nous avons déjà rencontré la densité conjointe des données observées et des
données latentes. On en déduit que la log-vraisemblance des données complètes s’écrit
N X
K


X
ln p y, x |Θ =
I(Xn =k) .
ln π(k) fk y n
n=1 k=1

Dans cette équation le logarithme est passé à l’intérieur de la sommation sur l’indice k uniquement car cette
somme se réduit à un seul terme du fait de la présence des fonctions indicatrices. La quantité intermédiaire
de l’algorithme EM s’écrit donc
N X
K

X
i

h

Q Θ|Θ(`) =
E I(Xn =k) y, Θ(`)
ln π(k) fk y n

(3.8)

n=1 k=1




avec fk y n = fk y n |µk , Σk .
Le dernier terme de cette équation s’écrit
h
i

E I(Xn =k) y, Θ(`)

=

K
X





I(Xn =j) p Xn = j y, Θ(`)

j=1

=





p Xn = k y, Θ(`)

=

γn(`) (k).
(`)

À l’itération `, le calcul de la quantité intermédiaire de l’algorithme EM se réduit donc au calcul de γn (k)
pour 1 ≤ n ≤ N et 1 ≤ k ≤ K. Ces quantités peuvent être calculées simplement par application de la
formule de Bayes












p Xn = k Θ(`) p y n Xn = k, Θ(`)
p Xn = k Θ(`) p y n Xn = k, Θ(`)


γn(`) (k) =
= K
,




X
p y n Θ(`)
(`)

p Xn = j Θ
p y n Xn = j, Θ(`)
j=1

que l’on peut écrire, étant entendu que les quantités du membre de droite sont calculées à partir des paramètres estimés à l’itération (`),

(`)
π (`) (k) fk y n
γn(`) (k) = K
.
(3.9)

X
(`)
(`)
π (j) fj
yn
j=1

En insérant ces valeurs dans (3.8), nous pouvons écrire
N X
K

X


Q Θ|Θ(`) =
γn(`) (k) ln π(k) fk y n .

(3.10)

n=1 k=1

3.3.2

Formules de ré-estimation des paramètres

C’est la maximisation de l’eq. (3.8) par rapport aux paramètres π(k), µk et Σk qui fournit les nouvelles
valeurs estimées pour l’itération ` + 1. Le problème est que cette quantité dépend de paramètres vectoriels
(µk ), voire matriciels (Σk ). Pour simplifier, nous proposons dans la suite les démonstrations concernant le
cas de mélanges scalaires.
25

ý Moyenne : En dérivant (3.10) par rapport à µk , nous écrivons


h
i
(`)
N ∂ γn
∂Q Θ|Θ(`)
(k) ln (fk (yn ))
X
=
∂µk
∂µk
n=1
h
i
(`)
N ∂ γn
(k) ln (fk (yn )) ∂ [ln (f (y ))]
X
k
n
=

[ln
(f
(y
))]
∂µ
k
n
k
n=1
N
X

=

yn − µk
,
σk2

γn(`) (k)

n=1

ce qui donne, en annulant la dérivée
N
X
(`+1)

µk

=

γn(`) (k) yn

n=1
N
X

.

(3.11)

γn(`) (k)

n=1

ý Variance : En dérivant (3.10) par rapport à σk2 , nous écrivons
i
h


(`)
N ∂ γn
(k) ln (fk (yn ))
∂Q Θ|Θ(`)
X
=
∂σk2
∂σk2
n=1
i
h
(`)
N ∂ γn
(k) ln (fk (yn )) ∂ [ln (f (y ))]
X
k
n
=
2

[ln
(f
(y
))]
∂σ
k
n
k
n=1


N
2
X
(yn − µk )
1
=
γn(`) (k)
− 2 ,
4
2σk
2σk
n=1
ce qui donne, en annulant la dérivée
N
X
(`+1) 2

(σk

) =

γn(`) (k)



(`+1)

yn − µk

2

n=1
N
X

.

(3.12)

γn(`) (k)

n=1

ý Le cas des poids des composantes de mélange π(k) est plus délicat car il faut prendre en compte
PK
la contrainte :
π(j) = 1 par le biais de la technique des multiplicateurs de Lagrange. On note


j=1


PK
(`)
P Θ|Θ
= Q Θ|Θ(`) + 1 − j=1 π(j) λk


∂P Θ|Θ(`)

h
i
(`)
N ∂ γn
(k) ln (π(k))
X

=

∂π(k)

∂π(k)

n=1
N
X

=

γn(`) (k)

n=1

d’où, en annulant la dérivée π (`+1) (k) =

1
λk

λk =

PN

N X
K
X

1
− λk ,
π(k)

(`)

n=1

− λk

γn (k). Or étant donné

γn(`) (j) =

n=1 j=1

N
X
n=1

26

1 = N,

PK

j=1

π(j) = 1, nous avons

(a)

(b)

Figure 3.6 – Deux images de classes utilisées pour l’initialisation des paramètres. (a) obtenue à partir d’un
initialisation au hasard ; (b) obtenue par une initialisation avec l’algorithme des k-moyennes.

d’où le résultat final
π (`+1) (k) =

N
1 X (`)
γ (k).
N n=1 n

(3.13)

On peut retenir cet algorithme en se souvenant que la probabilité a priori est re-estimée par la moyenne des
probabilités a posteriori, les moyennes sont re-estimées par les moyennes des observations « pondérées par
les probabilités a posteriori », et les variances sont re-estimées par les variances empiriques en considérant
des observations « pondérées par les probabilités a posteriori ».

3.4

Simulations et exemples en segmentation d’images

La segmentation supervisée de l’image 3.3(b) par décision bayésienne produit un taux d’erreur de 9.76%.
Voyons maintenant les taux d’erreur obtenus lorsque l’on estime les paramètres avec EM, à partir d’un
initialisation de ces paramètres selon les deux cas de figures suivants :
CAS 1 : initialisation au hasard. Nous avons généré une image de deux classes en procédant par tirage
aléatoire uniforme et équiprobable. L’image ainsi obtenue, observable sur la figure 3.6(a), produit un
taux d’erreur de 49.26% par rapport à l’image originale. En utilisant cette image de classes et l’image
bruitée, et des estimateurs à partir des données complètes pour estimer Θ(0) , nous avons obtenu, après


` = 300 itérations de EM, les estimées suivantes : {ˆ
π (1) = 0.68; π
ˆ (2) = 0.32}, µ
ˆ1 = 99.75; σ
ˆ12 = 34.706


et µ
ˆ2 = 109.97; σ
ˆ22 = 9.09 (à comparer avec les vrais paramètres). Le résultat de segmentation, obtenu
par décision bayésienne sur la base de ces estimées, est observable sur la figure 3.7(a).
CAS 2 : initialisation par l’algorithmes des k-moyennes. Ici, nous avons segmenté l’image avec l’algorithme des k-moyennes. L’image ainsi obtenue, observable sur la figure 3.6(b), produit un taux d’erreur de 12.87% par rapport à l’image originale. En utilisant cette image de classes et l’image bruitée, et des estimateurs à partir des données complètes pour estimer Θ(0) , nous avons obtenu, après


` = 300 itérations de EM, les estimées suivantes : {ˆ
π (1) = 0.32; π
ˆ (2) = 0.68}, µ
ˆ1 = 99.76; σ
ˆ12 = 34.71


et µ
ˆ2 = 109.97; σ
ˆ22 = 9.09 (à comparer avec les vrais paramètres). Le résultat de segmentation, obtenu
par décision bayésienne sur la base de ces estimées, est observable sur la figure 3.7(b).
Ces deux résultats montrent la grande robustesse de EM qui, à partir de deux points de départ Θ(0) très
différents, conduit a deux estimées Θ(100) très proches des vrais paramètres. Ce bon résultat est obtenu car
la fonction de vraisemblance de cet exemple est très lisse et présente peu de maxima locaux. Il est quand
même important de souligner le nombre d’itérations nécessaires dans le premier cas (cf. courbes présentées
dans la remarque ci-dessous).
27

(a) τ = 9.51%

(b) τ = 9.46%

Figure 3.7 – Segmentation par décision bayésienne suite à une estimation EM. (a) cas 1 ; (b) cas 2.
Remarque 3 : À titre d’illustration, la figure 3.8 dresse l’évolution des paramètres du mélange de la
figure 3.3(b) et de la quantité intermédiaire Q(.) de l’eq. 3.10 au cours de l’algorithme EM, à partir d’une
initialisation au hasard (cas 1). La convergence est lente mais illustre bien la capacité de EM à retrouver
les paramètres cachés, même très éloignés !
Remarque 4 : Pour illustrer le concept dans le cas d’un mélange 2D, nous avons segmenté une image de
galaxie constituée de 2 bandes spectrales, cf. figure 3.9. Les résultats de segmentation pour 2, 3 et 4 classes
sont présentés dans la figure 3.10.

3.5

Synoptique de l’algorithme complet de mélange aveugle

L’algorithme de la figure 3.11 reprend l’ensemble des étapes nécessaires pour réaliser la segmentation non
supervisée d’un signal numérique à l’aide du modèle de mélange aveugle développé dans ce chapitre.
Questions 1. Exercice 2 de l’examen de février 2011.

28

(a) Proportions π(k)

(b) Moyennes µk

(c) Variances σk2

(d) Q

Figure 3.8 – Tracés de l’évolution de l’estimation des paramètres du mélange de la figure 3.3(b) par EM.
Les traits horizontaux représentent les vrais valeurs vers lesquelles les courbes devraient converger.

(a) Bande J

(b) Bande R

Figure 3.9 – Les deux bandes d’une image bi-spectrale de galaxie.

29

(a) 2 classes

(b) 3 classes

(c) 4 classes

Figure 3.10 – Segmentation de l’image 3.9 en 2, 3 et 4 classes.

Require: y et K
1. Initialisation : Donner une valeur initiale aux paramètres.
Segmenter y −→ x(0) .

Estimer les paramètres sur les données complètes y, x(0) −→ Θ(0) .
2. Estimation EM : Trouver Θ(L) à partir de Θ(0) .
for ` = 1 to L do {À partir des paramètres Θ(`) }
(`)
Calculer les probabilités a posteriori : γn (.).
(`+1)
(`+1)
Estimer des paramètres du bruit : µk
et Σk
.
Estimer des probabilités a priori : π (`+1) (.) .
end for
3. Segmentation : Appliquer la décision bayésienne.
Figure 3.11 – Segmentation non supervisée par mélange aveugle gaussien.

30

Chapitre 4

Chaînes de Markov cachées
Comme nous avons pu le voir précédemment, les résultats des deux chapitres précédents (Décision bayésienne
et EM sur données indépendantes) peuvent s’appliquer directement à la segmentation d’images en procédant
« pixel par pixel ». Soit S = [1, N ] l’ensemble des pixels, qui se présente généralement sous forme d’une grille
rectangulaire. Pour chaque observation y n (un nombre pour une image numérique) sur le pixel n ∈ S, on
applique la règle de l’eq. (2.5), page 12 qui donne sˆB (y n ) = xn . Nous avons vu que cette utilisation de la
stratégie Bayésienne assure une optimalité de la classification. Cependant, la classification est faite « point
par point » (on ne tient pas compte des « voisins ») et l’optimalité concerne uniquement ce type de démarche.
Intuitivement, lorsque l’on classe un pixel n, regarder des pixels voisins de n devrait apporter de l’information
supplémentaire. Notons bien que nous ne connaissons pas les classes des pixels voisins, mais uniquement les
observations faites sur ces pixels. Tenir compte des valeurs des autres pixels revient à supposer que les
différentes variables aléatoires associées aux pixels ne sont pas nécessairement indépendantes. L’objectif des
modélisations markoviennes est d’introduire des modèles permettant de tenir compte de cette dépendance. La
modélisation par champs de Markov est bien adaptée à la structure 2D des images (on parle de « dépendance
spatiale ») ; nous aborderons cette modélisation dans le chapitre suivant. La modélisation par chaîne de
Markov, présentée dans ce chapitre, est bien adaptée aux signaux numériques 1D et correspond bien souvent
à une « dépendance temporelle » (série chronologique d’observations). Les références sur ce modèle à consulter
sont en priorité [MZ97, Rab89].
Remarque 1 : Il est possible d’utiliser une modélisation par chaîne de Markov pour segmenter les images.
La démarche consiste à « déplier l’image » en un vecteur 1D. Ce dépliement peut être réalisé ligne par
ligne, colonne par colonne, ou bien en zig-zag à partir d’un coin (technique utilisée dans les méthodes de
codage par DCT). Il est aussi possible d’utiliser les « Space Filling Curves » [Sag94]. Une « courbe qui
remplit l’espace » est une courbe qui parcourt un espace multi-dimensionnel (inclus dans Z2 dans notre
cas) en passant une et une seule fois par tous les pixels de cet espace. Parmi celles-ci, la courbe de Hilbert
se construit de manière récursive, ce qui lui confère une structure auto-similaire, cf. figure 4.1. Elle possède
une propriété importante : le passage d’un point de la courbe à son voisin (sur la courbe) ne s’accompagne
jamais d’un déplacement supérieur à un pixel dans l’image. Ainsi la courbe reste longtemps dans une zone
restreinte de l’image.
Ces parcours sont inversibles, ce qui permet, une fois les traitements effectués sur le signal 1D, de reconstruire l’image 2D. Ainsi, quand nous parlerons d’image, il faudra comprendre l’agencement spatial des
pixels issu du dépliement de l’image en une chaîne grâce au parcours de Hilbert. Précisons enfin que ce
parcours a été étendu par W. Skarbeck [Ska92] pour prendre en compte des images dont les nombres de
lignes L et de colonnes C sont paires et non plus uniquement des puissances de deux comme dans la courbe
initiale.

On considère donc deux processus stochastiques :

31

1 2 3 4 5 6 7 8
1
2
3
4
5
6
7
8
y
y1

yn−1 yn yn+1

yN

Figure 4.1 – Construction du parcours de Hilbert pour une image 8 × 8. À chaque pixel de coordonnées
(l, c) est associée une abscisse n indiquant le numéro d’ordre dans le parcours de la courbe, depuis la position
initiale (1, 1) → n = 1 jusqu’à la position finale (1, 8) → n = 64.

? le processus observé Y = {Y 1 , . . . , Y N } 1 où chaque variable vecteur aléatoire Y n est à valeur dans RD ,
et
? le processus caché X = {X1 , . . . , XN } où chaque variable aléatoire est à valeur dans Ω = [1, K].
Les réalisations de X sont inobservables (justifiant l’appellation « caché ») et le problème de la segmentation
est celui de l’estimation de X = x à partir de l’observation Y = y (i.e. le signal numérique à classer). On
parle aussi de la restauration de X. Ce chapitre suit la démarche générale présentée dans l’introduction. Nous
abordons d’abord la définition de la loi du couple de processus (X, Y ) à partir de la loi a priori p (X) et la


famille P Y |x des lois de Y conditionnelles à X = x. Ceci nous permet d’établir que, conditionnellement à
Y , X est une chaîne de Markov non homogène. Puis, nous explicitons deux méthodes de décision bayésienne,
i.e. le « MPM » et le « MAP », pour lesquels nous donnons des procédures exactes de calcul (c’est à dire
sans avoir recours à des simulations). Nous abordons enfin le problème de l’estimation des paramètres du
modèles par EM [DLR77, MK96, Bil97], en établissant les équations équivalentes au modèle de mélange
indépendant présenté dans le chapitre précédent. La dernière section est consacrée aux exercices et aux
simulations numériques.

4.1

Le modèle de chaîne de Markov cachée

Le processus aléatoire X est une variable aléatoire discrète prenant ses valeurs dans ΩN . Le cardinal de ce
dernier ensemble est trop important pour qu’on puisse calculer la probabilité de chaque élément (nombre de
réalisations possibles est K N ) ; il est ainsi impossible de traiter le problème de la classification bayésienne
directement. Exemple pour une image : N = 256 × 256 et K = 4.
Un des moyens pour contourner cette difficulté est de supposer l’indépendance des variables Xn (chapitre
précédent), mais sa justification intuitive reste difficile dans de nombreuses applications. Le modèle des
chaînes de Markov permet de contourner cette difficulté. Ainsi, le modèle de chaîne de Markov cachée est
plus riche que le modèle de mélange dans le sens où il permet de rendre compte des interactions temporelles
en substituant à l’hypothèse de V.A. i.i.d. celle d’une évolution markovienne.
1. La notation Y 1:n désignera la séquence de V.A.





Y 1, . . . , Y n .

32

4.1.1

Loi de X a priori

Une série de variables aléatoires X = {X1 , . . . , XN } prenant leurs valeurs dans l’ensemble des états (ou
classes) Ω = [1, K] est une chaîne de Markov si elle vérifie pour tout n ≥ 1 :


p xn x1:n−1 = p (xn |xn−1 ) .

(4.1)

La loi du processus X est alors déterminée par la loi de X1 , dite « loi initiale », π(k) = p (X1 = k) 2 et
la suite de matrices de transitions a(n) (k, l) = p (xn+1 = l |xn = k ). Nous supposerons dans la suite que la
chaîne de Markov est homogène (indépendante de l’indice temporel) et qu’en particulier
c(n) (k, l) = p (xn = k, xn+1 = l)
ne dépendent pas de n (on note alors c(k, l)). La loi initiale π(k) = p (X1 = k) est alors donnée par
π(k) =

K
X

c(k, l).

l=1

La matrice de transitions A, dont les éléments s’écrivent a(k, l) = p (X2 = l |X1 = k ), est alors également
indépendante de n et donnée par
c(k, l)
a(k, l) =
.
π(k)
La loi de probabilité du vecteur aléatoire X s’écrit donc
p (x)

=

p (x1 ) p (x2:N |x1 )

= p (x1 ) p (x2 |x1 ) p (x3:N |x1 , x2 )
= p (x1 ) p (x2 |x1 ) p (x3:N |x2 )
= ...
= p (x1 ) p (x2 |x1 ) . . . p (xn |xn−1 ) .
que nous pouvons finalement écrire

p (x) = π(x1 )

N
Y

a (xn−1 , xn ).

(4.2)

n=2

Simulations 1 : Les courbes (a) et (b) de la figure 4.2 montrent deux simulations de chaînes de Markov
à K = 3 classes, avec les mêmes probabilités initiales données par
π(1) = 0.1,

π(2) = 0.55,

π(3) = 0.35,

et les matrices de transitions respectives suivantes :



0.7
A1 = 0.3
0.2

0.1
0.6
0.3



0.2
0.1 ,
0.5



0.2
A2 = 0.4
0.3

0.4
0.3
0.5



0.4
0.3
0.2

La simulation de x se fait par tirage successif : xn est le résultat du tirage selon les probabilités
{a(xn−1 , k)}k∈Ω a .
a. Pour x1 , on tire selon les probabilités a priori.
2. Les probabilités initiales ont le même rôle que les probabilités a priori du chapitre précédent.

33

(a) CM - A1

(b) CM - A2

Figure 4.2 – Les courbes montrent deux simulations de chaînes de Markov à K = 3 classes, avec les
probabilités initiales et les matrices de transition données dans le texte.

4.1.2

Loi de (X, Y )

Il nous faut maintenant définir les lois de Y conditionnelles à X. Pour cela, nous allons émettre deux
hypothèses
H1 Conditionnellement aux états {Xn }, les observations {Yn } sont mutuellement indépendantes :
N


Y
p y n |x .
p y |x =

(4.3)

n=1

H2 Chaque observation Yn ne dépend que de l’état Xn au même instant :




p y n |x = p y n |xn .

(4.4)

Ainsi la loi du couple (X, Y ) est définie par :
p x, y



= p (x)

N


Y
p y n |xn ,
n=1



= p (x1 ) p y 1 |x1

N
Y



p (xn |xn−1 ) p y n |xn ,

n=2

soit finalement
N
Y


p x, y = π(x1 ) fx1 y 1
a(xn−1 , xn ) fxn y n .

(4.5)

n=2




Les densités conditionnelles fxn y n = p y n |xn sont également appelées lois d’attache aux données ou
lois d’émission selon que le modèle est utilisé en traitement du signal ou en traitement des images. Elles
représentent notre connaissance sur la dispersion des niveaux de gris associée à chaque classe (bruit, variation
naturelle des valeurs. . . ). Ces K densités sont généralement considérées gaussiennes mais d’autres formes de
distributions (paramétrique ou non) restent possibles. De manière locale, les hypothèses permettent d’écrire
directement à partir de (4.3)





p y n , y n+1 |xn , xn+1 = fxn y n fxn+1 y n+1 ,
(4.6)
34

Figure 4.3 – Représentation de la structure de dépendance d’une CMC à bruit indépendant.
et donc





p xn , xn+1 , y n , y n+1 = c(xn , xn+1 ) fxn y n fxn+1 y n+1 .
On note également que



p xn+1 , y n+1 xn , y n
=


=





p y n , y n+1 , xn+1 |xn
p yn , y n+1 |xn , xn+1




=
p (xn+1 |xn )
p y n |xn
p y n |xn


p y n+1 |xn+1 p (xn+1 |xn ) ,

en utilisant l’eq. (4.6), soit






p xn+1 , y n+1 xn , y n = a(xn , xn+1 ) fxn+1 y n+1 .

(4.7)

Ce dernier résultat aboutit au schéma de la figure 4.3 qui représente la structure de dépendance entre les
processus X et Y selon les deux hypothèses énoncées ci-dessus.
De la même manière, on montre que








p xn+1 , y n+1 x1:n , y 1:n = p xn+1 , y n+1 xn , y n ,
c’est à dire que le processus joint Z = (X, Y ) (zn = (xn , yn )) est markovien homogène, tout comme X, à la
différence que son espace d’états n’est pas fini (il est mixte entre Ω et R). Ainsi, il est aisé de montrer que
le processus Y |X est une chaîne de Markov homogène.
Par contre, le processus observé Y seul n’est pas markovien puisque



p y n y 1:n−1

=

K



X

p y n , Xn = k y 1:n−1
k=1

=

K




X


y 1:n−1
p y n Xn = k,
p Xn = k y 1:n−1

k=1

=

K
X





fk y n p Xn = k y 1:n−1 .

k=1





Par analogie avec le filtre le Kalman, le terme p Xn = k y 1:n−1 apparait comme une prédiction de l’état
suivant connaissant les observations
passées.
Ainsi
la loi de Y n conditionnellement à son passé est un modèle




de mélange dont les poids p Xn = k y 1:n−1 dépendent du passé complet du signal (et pas seulement de
Y n−1 ). Nous verrons ce qu’il est de la markovianité de la loi de X a posteriori (X |Y ) dans la section 4.1.4,
page 38.
35

(a) Bruitage de 4.2(a)

(b) Bruitage de 4.2(b)

(c) Histogramme (a)

(d) Histogramme (b)

Figure 4.4 – (a) et (b) Bruitage des deux chaînes de la fig. 4.2 avec les mêmes densités gaussiennes que
pour les simulations de la fig. 3.1, page 18. Les figures (c) et (d) montrent les histogrammes dans les deux
cas mais avec des simulations de 5000 échantillons.

Simulations 2 : Les courbes (a) et (b) de la fig. 4.4 montrent le résultat du bruitage des deux chaînes de
Markov de la fig. 4.2. Une fois la chaîne de Markov x simulées, la valeur bruitée yn de chaque état xn est
obtenue en effectuant un tirage selon la loi de densité fxn (.). On parle de bruit indépendant. Les courbes
(c) et (d) montrent les histogrammes pour ces deux cas mais avec des simulations de 5000 échantillons. Il
est intéressant de remarquer que les deux mélanges sont sensiblement différents entre eux, et également
différent du mélange avec les mêmes paramètres mais sous hypothèse d’indépendance (cf. fig. 3.1(d),
page 18). Ceci montre l’influence de la markovianité sur le mélange.

4.1.3

Probabilités « forward-backward »

Nous définissons dans ce paragraphe les probabilités « forward-backward » ou « avant-arrière » qui jouent
un rôle crucial, aussi bien au niveau de la décision (section 4.2) que de l’estimation des paramètres par
EM (section 4.3). Nous définissons directement la version conditionnelle des probabilités « avant-arrière »,
proposée par P. A. Devijver [Dev85] pour pallier aux problème numériques d’« underflow » des probabilités
« avant-arrière » originales proposées par Baum et al. [BPSW70].
36

ý Notons αn (.) les probabilités « avant »






1


αn (k) = p Xn = k y 1:n =
p Xn = k, y n y 1:n−1 ,
Sn




avec Sn = p y n y 1:n−1 , n > 1 le facteur de normalisation. Elles se calculent par récursion
? n=1:


p X1 = k, y 1
π(k) fk y 1

.
= K
∀k ∈ Ω, α1 (k) = p X1 = k y 1 =
X
p y1

π(l) fl y 1

l=1

? 1≤n<N −1 :
∀k ∈ Ω,

αn+1 (k)

=
=
=

=

=

1





p Xn+1 = k, y n+1 y 1:n



Sn+1





1
p y n+1 Xn+1 = k, y 1:n
p Xn+1 = k y 1:n


Sn+1
1
Sn+1



p y n+1 |Xn+1 = k

K
X



p Xn = l, Xn+1 = k y 1:n



l=1

1
fk y n+1
Sn+1



1
fk y n+1
Sn+1



K
X







l=1
K
X





p Xn+1 = k Xn = l, y 1:n
p Xn = l y 1:n
a (l, k) αn (l).

l=1





On voit apparaître, en seconde ligne, le terme de prédiction, qui s’écrit donc p Xn+1 = k y 1:n =
PK
l=1 a (l, k) αn (l).
ý Notons βn (.) les probabilités « arrière »





|X
=
k
p y n+1:N |Xn = k
p
y
n
1
n+1:N



=

.
βn (k) =


S
n+1
p y n+1:N y 1:n
p y n+2:N y 1:n+1

Elles se calculent par récursion
? n=N :
∀k ∈ Ω,

βN (k) = 1.

? 1≤n<N :
∀k ∈ Ω,

βn (k)

=

=

=

=

K

X

1

1

Sn+1




p y n+2:N y 1:n+1

l=1

1

1

K
X

Sn+1




p y n+2:N y 1:n+1

l=1

1
Sn+1
1
Sn+1

p y n+1:N , Xn+1 = l |Xn = k





K
X




a(k, l) p y n+1 |
Xn =
k, Xn+1


a(k, l) fl y n+1



βn+1 (l).

l=1

37



p (Xn+1 = l |Xn = k ) p y n+1 , y n+2:N |Xn = k, Xn+1 = l

l=1
K
X






p y n+2:N |
Xn =
k, Xn+1 = l



=l

p y n+2:N y 1:n+1



4.1.4

Loi de X a posteriori

Intéressons-nous maintenant à la loi de X a posteriori et calculons
a
˜n (k, l) = p Xn+1




p Xn+1 = l, Xn = k y

= l Xn = k, y = PK
.

j=1 p Xn+1 = j, Xn = k y

Or
X


p Xn+1 = l, Xn = k y =


p x y =

xn+2 ,...,xN ∈Ω

1

p y

X

(4.8)


p x, y .

xn+2 ,...,xN ∈Ω

Soit, en utilisant l’expression de la loi du couple (4.5),





X
a(k, l)fl y n+1
a(l, xn+2 )fxn+2 y n+2 . . . a(xN −1 , xN )fxN y N
a
˜n (k, l) =

xn+2 ,...,xN ∈Ω




a(k, xn+1 )fxn+1 y n+1 . . . a(xN −1 , xN )fxN y N

X

.

xn+1 ,...,xN ∈Ω

Enfin, à partir de la définition récursive des probabilités arrières de l’eq. (4.8), le produit des termes
Sn+2 . . . SN au numérateur et au dénominateur s’annulant, nous obtenons


a(k, l) fl y n+1 βn+1 (l)
,
a
˜n (k, l) = K


X
a(k, m) fm y n+1 βn+1 (m)
m=1

soit finalement
a
˜n (k, l) =



a(k, l) fl y n+1 βn+1 (l)
Sn+1 βn (k)

.

(4.9)

Ainsi, la loi de X a posteriori est celle d’une chaîne de Markov non homogène dont la matrice de transition
est définie par les probabilités de l’eq. (4.9). La loi initiale par donnée par, cf. (4.13)

(4.10)
p X1 = k y = γ1 (k) = α1 (k) β1 (k).
Remarque 2 : L’intérêt de préciser la loi de X a posteriori se justifie lorsqu’on utilise les méthodes
d’estimation SEM [CD85] et ICE [Pie94], qui peuvent être utilisées à la place de l’algorithme EM. Celles-ci
font en effet appel à des réalisations de X a posteriori, simulables grâce aux équations (4.10) et (4.9).

4.2

Décision bayésienne

L’estimation de X à partir de Y est réalisée en appliquant une règle de décision bayésienne. Comme présenté
dans le premier chapitre, cette règle est construite à partir d’une fonction L qui mesure la perte entre x et
ˆ . Comme on le sait, l’estimateur bayésien de X correspondant à la perte L est celui qui
son estimé sˆ(y) = x
minimise le coût moyen :


ˆ = arg min E L (X = x, x
ˆ )| y
sˆ(y) = x
ˆ ∈ΩN
x

=

arg min

ˆ ∈ΩN
x

X


ˆ ) p X = x y .
L (X = x, x

x∈ΩN

En traitement du signal et des images, deux fonctions de pertes sont le plus souvent utilisées, aboutissant aux
estimateurs bien connus suivants : le MPM et le MAP. Comme nous allons maintenant le voir, les chaînes de
Markov permettent le calcul explicite des solutions bayésiennes de classification, contrairement au modèle
des champs de Markov qui fait appel à des méthodes d’approximation de type Monte-Carlo (MPM) ou recuit
simulé (MAP).
38

4.2.1

Critère du MPM

La fonction de coût de l’estimateur du Maximum a Posteriori des Marges (MPM) est proportionnelle
au nombre de pixels mal classés
N
X
ˆ) =
L1 (x, x
1xn 6=xˆn ,
n=1

Or


ˆ )| y =
E L (X = x, x

X


ˆ ) p X = x y
L (X = x, x

x∈ΩN

=

N
X X


L (xn , x
ˆ n ) p x y

x∈ΩN n=1

=

N
X
X


L (xn , x
ˆ n ) p x y

n=1 x∈ΩN

=

N X
X

X


L (xn , x
ˆn ) p xn , xn y

n=1 xn ∈Ω xn ∈ΩN −1

=

N X
X


L (xn , x
ˆ n ) p xn y ,

n=1 xn ∈Ω

soit finalement

N
X



ˆ )| y =
E L (xn , x
ˆn )| y .
E L (X = x, x

(4.11)

n=1

L’espérance globale se décompose en l’espérance en chaque site et donc la décision peut se prendre en chaque
site. L’estimateur correspondant a donc pour expression
X

PM
x
ˆM
= arg min
L (xn , x
ˆn ) p Xn = xn y
n
x
ˆn ∈Ω

xn ∈Ω

X


p Xn = xn y

=

arg min

=


arg min 1 − p Xn = x
ˆ n y

x
ˆn ∈Ω
xn 6=x
ˆn

x
ˆn ∈Ω

Soit finalement
∀n ∈ [1 · · · N ],


PM
x
ˆM
= arg max p Xn = x
ˆ n y .
n
x
ˆn ∈Ω

(4.12)


Il est maintenant facile d’exprimer γn (k) = p Xn = k y , appelées « probabilités marginales a posteriori » 3
en fonction des probabilités avant et arrière selon :




p Xn = k, y n+1:N y 1:n



γn (k) =

p y n+1:N y 1:n





βn (k)




= p y n+1:N Xn = k, y 1:n p Xn = k y 1:n
p y n+1:N |Xn = k




y 1:n
p y n+1:N Xn = k,


.
= αn (k) βn (k)
p y n+1:N |Xn = k
3. Ce sont bien les mêmes que celles présentées dans le chapitre précédent.

39

Soit finalement
γn (k) = αn (k) βn (k) .

(4.13)

Nous en profitons pour définir l’expression des probabilités conjointes a posteriori ξn (k, i)

p Xn = k, Xn+1 = i y = a
˜n (k, i) γn (k), qui nous serviront par la suite et qui s’écrivent :

ξn (k, i) =

4.2.2



αn (k) βn+1 (i) a(k, i) fi y n+1
Sn+1

.

=

(4.14)

Critère du MAP

La stratégie Maximum A Posteriori (MAP) consiste à pénaliser de la même façon toutes les configurations estimées qui sont différentes de x. La fonction de coût correspondante est donnée par
ˆ ) = 1{x6=xˆ } ,
L2 (x, x
et l’estimateur correspondant a pour expression
ˆ M AP
x

=
=
=

arg min

ˆ ∈ΩN
x

arg min

ˆ ∈ΩN
x

X


ˆ ) p X = x y
L (x, x

x∈ΩN

X


p X = x y

ˆ
x6=x


ˆ y .
arg min 1 − p X = x
ˆ ∈ΩN
x

Soit finalement



ˆ M AP = arg max p X = x
ˆ y = arg max p y, X = x
ˆ = arg max p y |X = x
ˆ p (X = x
ˆ) .
x
ˆ ∈ΩN
x

ˆ ∈ΩN
x

ˆ ∈ΩN
x

(4.15)

Le MAP a pour principal inconvénient de pénaliser aveuglément les erreurs, indépendamment du nombre de
sites erronés.
La question que l’on se pose est la suivante : « Étant donné Θ, comment calculer de manière
n efficace la
o
séquence d’états x = {x1 , . . . , xN }, pour laquelle la probabilité d’une séquence observée y = y 1 , . . . , y N
est maximale ? » L’algorithme de Viterbi [For73, Rab89] (méthode de « programmation dynamique ») donne
la solution de ce problème d’estimation. La présentation de l’algorithme que nous donnons s’inspire largement
de celle présentée dans [TK06], pages 309-311.
La figure 4.5 montre un diagramme où chaque point d’une colonne représente l’une des K classes. Les colonnes
successives correspondent aux observations successives. Les flèches désignent les transitions d’une classe vers
une autre, comme les observations sont obtenues séquentiellement. Ainsi, chaque xn correspond à un chemin
spécifique de transitions successives. Chaque transition d’une classe i vers une classe j est caractérisée par
la probabilité p (j |i ).
La tâche qui consiste à maximiser l’eq. (4.15) peut s’écrire de la manière suivante. Pour chercher le chemin
optimal, nous devons associer un coût à chacune des transitions, cf. eq. (4.7), page 35 :





d (xn−1 = i, xn = j) = p xn = j, y n xn−1 = i, y n−1 = a(i, j) fj y n .

Nous noterons, de manière abusive, d (x0 , x1 = j) = π(j) fj y 1 le coût initial.
40

Figure 4.5 – Diagramme en trellis pour l’algorithme de Viterbi.
On cherche à calculer le coût sur un chemin du diagramme noté c. En désignant xc(n) la classe du ne élément
du chemin c, le coût total d’un chemin c s’écrit :
P =

N
Y


d xc(n−1) , xc(n) ,

n=1

en utilisant l’eq. (4.5). Nous chercherons plutôt a optimiser la quantité équivalente suivante :
D = ln P =

N
X

N

X
ln d xc(n−1) , xc(n) ≡
d xc(n−1) , xc(n) .

n=1

n=1

Le coût d’atteindre l’état xc(n) à l’échantillon n <= N , s’écrit
n
X

D xc(n) =
d xc(m−1) , xc(m) .
m=1

Le coût maximum d’atteindre l’état k au ne élément du chemin se décompose récursivement de la manière
suivante :


δn (k) = max δn−1 (j) + d xc(n−1) = j, xc(n) = k ,
(4.16)
j∈Ω

avec ∀k ∈ Ω, δ0 (k) = 0. Il est maintenant immédiat de voir que le chemin optimal se termine par la classe :
x?N = arg max δn (k).
k∈Ω

(4.17)

En revenant à la figure 4.5, on s’aperçoit qu’il y a K transitions possibles pour chaque nœud xn . La relation
récursive de l’eq. (4.16) suggère que, dans la recherche d’un maximum, nous n’avons besoin de conserver
41

1. Récurrence
Pour n = 1



δ (k) = π(k) f y
1
k
1
ψ1 (k) = 0;

Pour 2 ≤ n ≤ N



δn (k) = max δn−1 (j) + d xc(n−1) = j, xc(n)=k
j∈Ω


ψn (k) = arg max δn−1 (j) + d xc(n−1) = j, xc(n)=k
j∈Ω

2. Décodage
Pour n = N
x?N = arg max δN (k)
k∈Ω

Pour 1 ≤ n < N
x?n = ψn+1 x?n+1



Figure 4.6 – L’algorithme de Viterbi.

qu’une seule de ces transitions pour chaque nœud : celle qui donne le coût maximal δn (k). Ainsi, à chaque
étape il y a uniquement K chemins possibles. Ainsi le nombre d’opérations est K pour chaque nœud, K 2
pour tous les nœuds, et donc N K 2 au total. Ce dernier nombre est à comparer avec la méthode brutale qui
comporte N K N étapes. L’algorithme est repris entièrement dans la figure 4.6.
Remarque 3 : L’estimateur MPM minimise, en moyenne, le nombre d’observations mal classées, ce qui
correspond au taux d’erreur de classification, contrairement au MAP qui pénalise de la même façon une
erreur sur un pixel et une erreur sur plusieurs pixels.

4.3

L’algorithme EM



Comme dans le cas indépendant, nous commençons par expliciter la quantité Q Θ|Θ(`) , puis nous donnons
les formules de re-estimation (appelées « formules de ré-estimation de Baum-Welch ») de tous les paramètres
du modèle des chaînes de Markov cachées.

4.3.1

Quantité intermédiaire

En procédant comme au paragraphe 3.3.1, on écrit le logarithme de la vraisemblance complétée selon :
"
ln p y, x |Θ



=



ln π(x1 )fx1 y 1

N
Y



a(xn−1 , xn )fxn y n



#

n=2

=

ln [π(x1 )] +

N
X

N
h
i X
ln fxn y n +
ln [a(xn−1 , xn )].

n=1

n=2

42

Ainsi la quantité intermédiaire s’écrit

h
i


= E ln p y, x |Θ y, Θ(`) ,
Q Θ|Θ(`)
K
X

h
i

ln [π(k)] E I(X1 =k) y, Θ(`)

k=1

=

+

K X
N
X

i
h i h

ln fk y n E I(Xn =k) y, Θ(`)

k=1 n=1

+

K X
K X
N
X

i
h

ln [a(k, i)] E I(Xn−1 =k,Xn =i) y, Θ(`)

k=1 i=1 n=2

soit
K

X
(`)
Q Θ|Θ(`) =
γ1 (k) ln [π(k)]

+

K X
N
X

h i
γn(`) (k) ln fk y n

k=1 n=1

k=1

+

K X
K X
N
X

ξn(`) (k, i) ln [a(k, i)],

(4.18)

k=1 i=1 n=2
(`)

(`)

en notant respectivement γn (.) et ξn (., .) les quantités des équations (4.13) et (4.14), lorsqu’elles sont
calculées à partir des probabilités αn (.), βn (.), a(., .) et f. (.) obtenues à l’itération `.

4.3.2

Ré-estimation des paramètres

La forme de la quantité intermédiaire de l’algorithme EM est donc très proche de celle obtenue pour le
modèle de mélange en (3.10). La partie maximisation par rapport aux différents paramètres du modèles fait
intervenir trois quantités différentes.
ý Pour les paramètres des densités gaussiennes Υ, c’est le premier terme de (4.18) qui intervient. Dans
la mesure où celui-ci est exactement l’analogue de celui obtenu pour le mélange en (3.10), les équations
de ré-estimation des vecteurs moyens et des matrices de covariances sont donc exactement identiques aux
équations (3.11) et (3.12) :
N
X

µ(`+1)
=
k

γn(`) (k) y n

n=1
N
X

,

(4.19)

γn(`) (k)

n=1

et
N
X
(`+1)

Σk

=

γn(`) (k)



y n − µ(`+1)
k



y n − µ(`+1)
k

n=1
N
X

0
.

(4.20)

γn(`) (k)

n=1

ý Intéressons aux paramètres de la loi de X (Ψ). Pour la distribution initiale π, c’est le troisième terme
PK
de (4.18) qu’il est nécessaire de considérer. Compte tenu de la contrainte j=1 π(j) = 1, la résolution se fait
comme dans le cas du modèle de mélange, et on obtient une formule analogue à (3.13) :
π (`+1) (k) =

N
1 X (`)
γ (k).
N n=1 n

43

(4.21)

Enfin, pour la matrice de transition, chaque ligne vérifie une contrainte de normalisation qui conduit à
introduire autant de multiplicateurs de Lagrange. En procédant comme dans le cas des mélanges (pour les
poids du mélange), on trouve
N
−1
X

a(`+1) (k, i) =

ξn(`) (k, i)

n=1
N
−1
X

.

(4.22)

γn(`) (k)

n=1

Comme dans le cas d’un mélange aveugle, l’algorithme EM permet d’obtenir des formules explicites de
ré-estimation des paramètres, ce qui en fait l’un de ses principaux intérêts. En contre-partie, l’algorithme
converge vers un maximum local de la vraisemblance, ne garantissant pas l’optimalité des paramètres trouvés.
D’autres méthodes d’estimation existent, telles que SEM [CD85] (Stochastic EM ) et ICE [Pie94] (Iterated
Conditional Estimation). Ces méthodes estiment certains des paramètres en utilisant des estimateurs à partir
des données complètes, en simulant X conditionnellement Y .

4.4

Simulations et exemples en segmentation d’images

Nous reprenons exactement les mêmes données que celles de la section 3.4, page 27, et nous présentons des
expériences numériques similaires. Ainsi, dans cette expérience, les observations sont des scalaires et les lois
d’attache aux données 1D.
À partir de l’image bruitée de la figure 3.3(b), nous avons estimé les paramètres de la chaîne de Markov
cachée avec EM. Pour cela, nous avons initialisé les paramètres Θ(0) soit (CAS 1) avec une initialisation
des classes au hasard (cf. fig. 3.6(a)), soit (CAS 2) avec une initialisation des classes par l’algorithme des
k-moyennes (cf. fig. 3.6(b)) :
CAS 1 - Initialisation au hasard En utilisant cette image de classes et l’image bruitée,
et des estimateurs à partir des données complètes pour estimer Θ(0) , nous avons obtenu, après ` = 300 itérations de EM, les estimées suivantes : {ˆ
π (1) = 0.67; π
ˆ (2) = 0.33},



a(1, 1) = 0.99, a
ˆ(1, 2) = 0.01, a
ˆ(2, 1) = 0.02, a
ˆ(2, 2) = 0.98},
µ
ˆ1 = 95.12; σ
ˆ12 = 41.25
et


2
µ
ˆ2 = 109.61; σ
ˆ2 = 66.71 (à comparer avec les vrais paramètres). Les résultats de segmentation
MPM et MAP, obtenus par décision bayésienne sur la base de ces estimées, sont observables sur les
figures (a) et (b) de 4.7.
CAS 2 - Initialisation par K-moyennes En utilisant cette image de classes et l’image bruitée, et des estimateurs à partir des données complètes, nous avons obtenu, après ` = 300 itérations de EM, les estimées
suivantes : {ˆ
π (1) = 0.67; π
ˆ (2) = 0.33}, {ˆ
a(1, 1) = 0.99, a
ˆ(1, 2) = 0.01, a
ˆ(2, 1) = 0.02, a
ˆ(2, 2) = 0.98},




µ
ˆ1 = 95.13; σ
ˆ12 = 41.40 et µ
ˆ2 = 109.68; σ
ˆ22 = 66.23 (à comparer avec les vrais paramètres). Les résultats de segmentation MPM et MAP, obtenus par décision bayésienne sur la base de ces estimées, sont
observables sur les figures (c) et (d) de 4.7.
Remarque 4 : À titre d’illustration, la figure 4.8 dresse l’évolution des paramètres du mélange de la
figure 3.3(b), page 20 et de la log-vraisemblance au cours de l’algorithme EM, à partir d’une initialisation
au hasard (cas 1). Il est remarquable de constater que la convergence numériques des paramètres est bien
plus rapide que dans le cas indépendant.
Remarque 5 : Pour illustrer le concept dans le cas d’un mélange 2D, nous avons segmenté la même image
de galaxie que celle utilisée pour le cas de mélange indépendant (cf. figure 3.9, page 29). Les résultats
de segmentation (critère MPM) pour 2, 3 et 4 classes sont présentés dans la figure 4.9. Ils offrent une
segmentation plus homogène (moins de pixels isolés) que la méthode de mélange aveugle.

44

(a) MPM - τ = 1.90%

(b) MAP - τ = 1.95%

(c) MPM - τ = 1.90%

(d) MAP - τ = 1.95%

Figure 4.7 – Segmentation par décision bayésienne suite à une estimation du mélange par EM/CMC. (a)
et (b) cas 1 ; (c) et (d) cas 2.

4.5

Synoptique de l’algorithme des CMC

L’algorithme de la figure 4.10 reprend l’ensemble des étapes nécessaires pour réaliser la segmentation non
supervisée d’un signal numérique à l’aide du modèle de chaîne de Markov cachée. Pour les images, une étape
préalable consiste à convertir l’image 2D en un signal 1D, avec le parcours de Peano (cf. fig 4.1), et à réaliser
l’opération inverse sur le signal segmenté, pour reconstituer l’image segmentée.

45

(a) Proportions π(k)

(b) Coeff. a(1, 1)

(c) Moyennes µk

(d) Variances σk2

(e) Log-vraisemblance

Figure 4.8 – Tracés de l’évolution de l’estimation des paramètres du mélange de la figure 3.3(b) par
EM/CMC. Les lignes représentent les vrais valeurs.

46

(a) 2 classes

(b) 3 classes

(c) 4 classes

Figure 4.9 – Segmentation de l’image 3.9 en 2, 3 et 4 classes (critère MPM).

Require: y et K
1. Initialisation : Donner une valeur initiale aux paramètres.
Segmenter y −→ x(0) .

Estimer les paramètres sur les données complètes y, x(0) −→ Θ(0) .
2. Estimation EM : Trouver Θ(L) à partir de Θ(0) .
for ` = 1 to L do {À partir des paramètres Θ(`) }
Calculer les probabilités « avant-arrière » : αn (.)(`) et βn (.)(`) .
(`)
(`)
Calculer les probabilités a posteriori : γn (.) et ξn (., .).
(`+1)
(`+1)
Estimer des paramètres du bruit : µk
et Σk
.
(`+1)
(`+1)
Estimer des paramètres de Markov : π
(.) et a
(., .).
end for
3. Segmentation : Appliquer une décision bayésienne.
(L)
MPM : directement à partir de γn (.).
MAP : algorithme de Viterbi (utilisant Θ(L) ).

Figure 4.10 – Segmentation non supervisée par EM/CMC gaussien.

47

48

Chapitre 5

Champs de Markov cachés
Tout comme nous l’avons fait dans le chapitre précédent, nous allons tenir compte des valeurs des pixels
voisins, c’est à dire que nous allons supposer que les différentes variables aléatoires associées aux pixels ne
sont pas nécessairement indépendantes. Si le voisinage se définit par l’échantillon précédent dans le modèle
par chaîne, il est défini sur une fenêtre 2D, dans le modèle par champs. Ce modèle est donc particulièrement
bien adapté au traitement des images et a suscité beaucoup d’intérêts depuis les travaux précurseurs de D.
Geman [GG84], J. Besag [Bes86], B. Chalmond [Cha89] et J. Marroquin [MMT87] dans la fin des années 80.
Une synthèse intéressante sur les champs de Markov cachés et leurs applications en traitement des images
est proposée dans le livre de S. Z. Li [Li01] et dans l’article de P. Pérez [P´98].
On considère donc toujours deux processus stochastiques : le processus observé Y et le processus caché X.
Le problème que l’on se pose reste l’estimation d’une réalisation de X à partir de l’observation y. Dans ce
modèle, les processus ne sont plus des vecteurs aléatoires, mais des champs aléatoires. Pour X :

 X1,1 ,
X=
...,

XL,1 ,


. . . , X1,C 
...,
...

. . . , XL,C

Pour éviter les indices de ligne et de colonne, et lorsqu’il n’y aura pas de confusion possible, la notation
suivante sera utilisée Xs , s ∈ S, S désignant l’ensemble des sites du champs. Le nombre total de variables
aléatoires est égale à N = L × C, c’est à dire le nombre de pixels de l’image que l’on souhaite traiter.
En nous référant toujours à la démarche proposée en introduction, nous établissons les lois de X a priori
et a posteriori en introduisant les hypothèses nécessaires pour appliquer des règles de décisions bayésiennes
(MAP et MPM) et estimer l’ensemble des paramètres du modèle.

5.1
5.1.1

Loi de X a priori : champs de Markov
Champs de Gibbs, champs de Markov

Le champs X est markovien s’il vérifie l’hypothèse suivante :
Hypothèse 1 : Pour tout s ∈ S, p (Xs |Xt , t 6= s ) = p (Xs |Xt , t ∈ Vs ) ou Vs est un « voisinage » de s dont
la forme géométrique est indépendante de s 1 .
1. On ne considère généralement que le 4-voisinage et le 8-voisinage.

49


Documents similaires


Fichier PDF polycoursdb
Fichier PDF foreground background segmentation using temporal and spatial
Fichier PDF projet dona gombert lecocq
Fichier PDF 17schaeffer et al
Fichier PDF martingales
Fichier PDF rapport projet 2a giraud remi


Sur le même sujet..