Pierre Alquier

Analyse bayésienne de modèles de mélanges:
Approche par MCMC à sauts réversibles

Mémoire de Dea Paris 6 / GT Ensae

Encadré par Mme Dominique Picard
Remis le 2 juin 2003

1

Table des matières
1 Introduction : l’inférence bayésienne dans les modèles de mélanges
1.1 Les modèles de mélanges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 Définition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.2 Commentaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.3 Première approche pour l’estimation par méthodes bayésiennes . . . . . . . . . . .
1.2 Les méthodes MCMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Principe général . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2 Algorithme de Hastings-Metropolis . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.3 Algorithme de Gibbs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.4 Utilisation en pratique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Inférence bayésienne dans les modèles de mélanges . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Les modèles de mélanges comme modèles à données manquantes . . . . . . . . . .
1.3.2 Utilisation de l’algorithme de Gibbs . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.3 Utilisation de l’algorithme de Hastings-Metropolis . . . . . . . . . . . . . . . . . .
1.3.4 Problèmes rencontrés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
3
3
4
6
6
6
7
9
9
9
10
11
12

2 Détermination du nombre de composantes par MCMC à saut réversible
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Formulation du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 Objectif et problème principal . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Les algorithmes MCMC à sauts réversibles . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Présentation générale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Principe des RJMCMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 La méthode retenue par Richardson et Green (1997) . . . . . . . . . . . . . . . . . . . . .
2.3.1 Cadre général, notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Spécification du modèle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.3 Méthode de simulation : choix des types de mouvements . . . . . . . . . . . . . . .
2.4 Les résultats obtenus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.1 Les données . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 Interprétation des résultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.3 Sensibilité au choix des lois a priori . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.4 Applications à la classification bayésienne . . . . . . . . . . . . . . . . . . . . . . .

12
13
13
13
13
13
14
17
18
18
20
25
26
26
28
28

3 Critiques et extensions
3.1 Performances du simulateur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 Commentaires des auteurs sur les résultats obtenus . . . . . . . . . . . . . . . . . .
3.1.2 Une proposition d’amélioration de la méthode . . . . . . . . . . . . . . . . . . . . .
3.1.3 Une approche "concurrente" : méthode par diffusion avec sauts . . . . . . . . . . .
3.1.4 Comparaison simulateurs à temps continu / RJMCMC . . . . . . . . . . . . . . . .
3.2 Problème des échanges des labels des paramètres . . . . . . . . . . . . . . . . . . . . . . .
3.3 Généralisations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Composantes non normales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Cas multivarié . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 A priori non informatifs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Critère de convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 Discussion du choix des méthodes bayésiennes . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.1 Complexité des MCMC dans ce contexte . . . . . . . . . . . . . . . . . . . . . . . .
3.5.2 Sensibilité du choix de modèle aux lois a priori . . . . . . . . . . . . . . . . . . . .
3.5.3 Une comparaison numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29
29
29
30
31
34
34
35
35
35
35
36
36
36
36
37

2

4 Mise en pratique
4.1 Le programme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Choix effectués . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 Problèmes rencontrés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Résultats sur données simulées . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Présentation des données . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Résultats obtenus à k fixé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.3 Résultats des RJMCMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Résultats sur données réelles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Intérêt de travailler sur des données réelles . . . . . . . . . . . . . . . . . . . . . .
4.3.2 Présentation des données . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.3 Résultats et commentaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

37
37
37
37
37
37
37
38
39
39
40
40

1

Introduction : l’inférence bayésienne dans les modèles de mélanges

1.1

Les modèles de mélanges

1.1.1

Définition

Les modèles de mélanges ont été initialement introduits afin de représenter l’existence de différentsgroupes relativement homogènes dans une population donnée.
Supposons que l’on étudie une quantité d’intérêt, x, dans une population composée en fait de k groupes
relativement homogènes, en proportions p1 , ..., pk . Dans chaque groupe, la quantité d’intérêt est supposéesuivre une loi de densité fθi (x). Ainsi, la quantité x sera répartie sur l’ensemble de la population suivant
une loi de densité :g(x) =

kX

pi fθi (x)i=1

On appelle "mélange fini de distributions" une telle densité. Les différents fθi (x) sont appelés les"composantes du mélange".
La figure 1 représente la densité d’un mélange de deux lois normales.

Fig. 1 – Mélange de deux lois normales de moyennes distinctes :"0.7N (−1, 1) + 0.3N (2, 1)"

1.1.2

Commentaires

D’un point de vue pratique, les paramètres θi et pi sont presque toujours inconnus. Le but de ladémarche statistique sera donc de proposer une méthode d’estimation de ces paramètres, ce que nous
allons examiner d’un point de vue bayésien dans la fin de cette première partie. Il est un certain nombre decas où la valeur de k elle-même sera inconnue. Le but de ce mémoire est de présenter un certain nombrede méthodes, issues de la statistique bayésienne, permettant de considérer k comme un paramètre dumodèle et de l’estimer. Ceci est l’objet de la méthode présentée dans ce mémoire.
Remarquons tout d’abord que le modèle, présenté tel quel, n’est pas identifiable. En effet, considérons
un mélange simple à deux composantes. On a évidemment :4

pfθ1 + (1 − p)fθ2 = (1 − p)fθ2 + pfθ1

Ainsi les couples (p, θ1 , θ2 ) et (1 − p, θ2 , θ1 ) définissent la même loi. Il sera donc nécessaire d’imposer
certaines contraintes supplémentaires afin d’assurer l’identifiabilité du modèle, par exemple une contrainted’ordre sur les poids : p > 1 − p ou plus généralement p1 > ... > pk .
Notons que les modèles de mélanges peuvent être utilisés dans diverses applications dépassant le cadre
présenté précédemment, par exemple :- dans un problème quelconque d’inférence statistique ou le recours à une famille paramétrique simple
n’est pas satisfaisant, l’introduction de mélanges de distributions permet d’enrichir le modèle considéré.
Dans ce cas, les différentes composantes ne s’interprètent plus comme reflétant la présence de groupes
dans la population étudiée. Cf. Robert (1996).
- les modèles de mélange sont aussi vus comme un compromis entre certaines méthodes paramétriques
et non-paramétriques. Par exemple, pour estimer la densité de la loi d’un échantillon X1 , ..., Xn devariables indépendantes, une approche paramétrique simple consisterait à utiliser une loi normale N (µ, σ 2 )alors qu’une approche non-paramétrique possible est l’utilisation de l’estimateur du noyau :
fd
(x) =

n
1 X x − Xi
ϕ()
nhn i=1hn

où ϕ est la densité d’une loi N (0, 1). Ainsi, on approche la loi des observations par des lois normales,
dans un cas en imposant la présence d’une seule normale, dans l’autre cas en imposant la présence de nnormales. Ainsi, le modèle de mélanges peut être vu comme un compromis entre les deux méthodes en ce
qu’il utilisera k normales où 1 ≤ k ≤ n, où k peut être soit fixé par le statisticien, soit estimé comme nous
le verrons dans la seconde partie. Cf. Robert (1996) ou McLachlan & Peel (2000). L’article de Roeder
et Wasserman (1997) utilise effectivement des modèles de mélanges pour faire de l’estimation de densité
non paramétrique.
- l’utilisation des modèles de mélanges est particulièrement intéressante pour le problème de la classification. On peut supposer qu’il existe effectivement k classes dans la population étudiée, mais ne passavoir de quel groupe est issu un individu particulier. Si l’on souhaite "deviner" quel est le groupe (ou laclasse) d’un individu donné à partir de l’observation de x et que l’on adopte le modèle de mélanges :g(x) =

kX

pi fθi (x)i=1

alors un individu pour lequel x = X est classé de façon optimale dans la classe i∗ = arg max1≤i≤k pi fθi (X).
Donc, estimer les paramètres du modèle permet de proposer une classification des individus. De plus,
l’estimation de k apporte aux méthodes plus classiques une estimation du nombre possible de classes
dans la population.
- enfin, les modèles de mélanges sont utilisés lorsqu’on est en présence d’observations "aberrantes".
Supposons par exemple que l’on étudie le modèle N (θ, 1) (avec θ proche de 0) et que l’on ait une
observation aberrante, valant 10, dans l’échantillon. On peut alors utiliser le mélange "p1 N (θ, 1) + (1 −p1 )N (0, 100)" pour éviter que l’observation aberrante ne vienne entacher d’erreur l’estimation de θ. La
figure 2 illustre ceci.
1.1.3

Première approche pour l’estimation par méthodes bayésiennes

L’utilisation de méthodes bayésiennes pour estimer les paramètres d’un mélanges de distributions sejustifie tout d’abord par l’impossibilité, au moins en première approche, d’utiliser les méthodes classiquescomme le maximum de vraisemblance : en général, la vraisemblance d’un modèle de mélanges n’est pas
bornée...
D’un point de vue bayésien, la densité du modèle de mélange présenté précédemment (avec k fixé)devient :g(x|ξ) =

k
X
i=1

5

pi f (x|θi )Fig. 2 – Mélange de deux lois normales de même moyenne mais de variance différentes ("chapeau méxicain") : "0.5N (0, 0.5) + 0.5N (0, 2)"soit une vraisemblance, pour n observations x1 , ..., xn :
n
Y

j=1

g(xj |ξ) =

k
n X
Y

j=1 i=1

pi f (xj |θi )où ξ = (θ1 , ..., θk , p1 , ..., pk ). On spécifie aussi une loi a priori sur le paramètre ξ, π(ξ). Remarquons
qu’il est en général impossible de spécifier une loi a priori impropre sur le paramètre ξ (si on le fait, laloi a posteriori peut elle aussi ne pas être une mesure de probabilités !). Ceci interdit l’utilisation des lois
non-informatives usuelles (loi de Jeffreys). Le choix des lois a priori sera donc crucial dans un modèle demélange.
Dès lors, la loi posteriori est donnée par :π(ξ)π(ξ|x1 , ..., xn ) = Z

k
n X
Y

j=1 i=1

π(ξ)

k
n X
Y

j=1 i=1

pi f (xj |θi )pi f (xj |θi )dξle principal problème pour la mise en pratique des méthodes bayésiennes sur des modèles de mélanges
finis apparaît alors : le calcul de la loi a posteriori est souvent impossible sous forme exacte (notammentl’intégrale au dénominateur), ce qui interdit aussi le calcul explicite des estimateurs bayésiens, de plus lesapproximations numériques sont souvent trop longues en pratiques (présence de k n termes au numérateur
et au dénominateur).
En réalité, les problèmes posés ci-dessus n’ont pu être résolus qu’avec le développement des méthodes
de Monte-Carlo par Chaînes de Markov (MCMC).

6

1.2

Les méthodes MCMC

1.2.1

Principe général

Le but des algorithmes MCMC est de simuler une chaîne de Markov de loi de probabilité stationnaire
donnée de densité f , notée par exemple (xn )n∈N . Si la chaîne est simulée de façon à satisfaire les hypothèses
voulues (irréductibilité, apériodicité) on pourra appliquer le théorème ergodique.
Par exemple, pour g avec Ef |g(X)| < ∞ on a :n1Xg(xn ) = Ef (g(X))
lim
n→∞ n
i=1

Intuitivement, on souhaite donc que (xn ) "visite toutes les zones probables" sous f . Il faudra donc
simuler xn jusqu’à n = n2 où "n2 est grand". De plus, le choix de x0 va avoir une influence importante
sur les premières valeurs de la suite (xn ) produite. En particulier, si x0 est située dans une zone "peuprobable" selon f , il en sera probablement de même pour les premières valeurs x1 , x2 ,... Donc il conviendra
d’enlever les n1 premières réalisations de (xn ). n1 est appelé burn-in period. Ainsi, on estimera Ef (g(X))
par :
n2
1 X
g(xn )n i=n1

En pratique, le choix de n1 et n2 est basé sur les observations de (xn ) et de la convergence des quantitésétudiées.
1.2.2

Algorithme de Hastings-Metropolis

Etant donné une loi de densité f par rapport à une mesure dominante µ, l’algorithme de Hastings
a pour but de simuler la réalisation d’une chaîne de Markov ergodique de loi stationnaire de densité f ,et ceci même si f n’est connue qu’à une constante multiplicative près. Ceci pourrait notamment être lecas de la loi a posteriori dans le cas de l’inférence bayésienne dans les modèles de mélanges présentée
précédemment :f (ξ) = π(ξ|x1 , ..., xn ) =

k
n X
Y
1
pi f (xj |θi )π(ξ)
C
j=1 i=1

Pour ceci, il est nécessaire de disposer d’une probabilité de transition Q de densité q(x|y) par rapport
à µ, vérifiant :- on peut simuler rapidement une variable aléatoire de loi q(.|y),- soit q est connue sous forme analytique, soit q est symétrique, c’est-à-dire que ∀x, ∀y, q(x|y) = q(y|x).
L’algorithme est alors le suivant. On se donne x0 arbitraire et on passe de xn à xn+1 ainsi :- on simule une variable x′n selon la loi q(.|xn ),- on calcule la probabilité d’acceptation
f (x′n ) q(xn |x′n )a(xn , x′n ) = min 1,
f (xn ) q(x′n |xn )puis avec une probabilité a(xn , x′n ) on "accepte" x′n c’est-à-dire que l’on prend xn+1 = x′n , et avec une
probabilité 1 − a(xn , x′n ) on rejette x′n , on prend alors xn+1 = xn . Ceci est réalisé le plus simplement
possible en simulant la réalisation d’une variable uniforme sur [0, 1], un , puis en comparant un aveca(xn , x′n ) : acceptation si un < a(xn , x′n ), rejet sinon.
Plusieurs remarques sont alors à faire. Tout d’abord, il faut vérifier si l’algorithme décrit ci-dessus
est bien un "algorithme MCMC", c’est-à-dire si il produit bien une chaîne de Markov de probabilité
stationnaire de densité f . Si c’est le cas, les xn produits nous serviront ensuite à étudier les propriétés def à l’aide de théorèmes limites sur les chaînes de Markov, comme le théorème ergodique.
En fait, il est évident que (xn ) est une chaîne de Markov, de probabilité de transition P (x, y) =a(x, y)q(y|x) + (1 − a(x, y))δx (y). Ceci permet de démontrer que f est bien une densité invariante pour(xn ), en utilisant le fait que f satisfait la detailed balance condition :7

f (y)P (y, x) = f (x)P (x, y)Ceci se vérifie facilement :f (x) q(y|x))q(x|y) + (1 − a(x, y))δx (y)f (y)P (y, x) = f (y) min(1,f (y) q(x|y)... = min(f (y)q(x|y), f (x)q(y|x)) + f (y)(1 − a(x, y))I(y = x)... = f (x) min(f (y) q(x|y), 1)q(y|x) + f (x)(1 − a(y, x))δy (x)f (x) q(y|x)⇒ f (y)P (y, x) = f (x)P (x, y)Démontrons donc maintenant que f est invariante. Il suffit pour cela de démontrer que, quel que soit
A ensemble mesurable :ZZ
P (x, A)f (x)dx =f (x)dxAOr :Z

P (x, A)f (x)dx =Z Z

P (x, y)f (x)dydx =A

Z Z

IA (y)P (x, y)f (x)dydx

On utilise alors la detailed balance condition f (y)P (y, x) = f (x)P (x, y) et on obtient :Z
Z Z
P (x, A)f (x)dx =IA (y)P (y, x)f (y)dydx

Toutes les quantités dans l’intégrale double étant positives on peut appliquer le théorème de Fubini :ZZZ
P (x, A)f (x)dx = IA (y)f (y)P (y, x)dx dyR
Or P (., .) étant une probabilité de transition on a ∀y, P (y, x)dx = 1 donc :ZZ
P (x, A)f (x)dx =f (y)dyA

ce qui est exactement ce qu’il fallait démontrer.
Enfin, si ∀x, ∀y, q(x, y) > 0 la chaîne de Markov (xn ) est irréductible. Elle est apériodique car il est
toujours possible en général que xn+1 = xn . Elle est récurrente positive car f est la densité d’une loi deprobabilités. On est donc dans les conditions d’applications du théorème ergodique.
1.2.3

Algorithme de Gibbs

Un autre algorithme, dit algorithme de Gibbs, permet de simuler une chaîne de Markov ergodique de loi
stationnaire ν(θ1 , ..., θk ) donnée, à condition de savoir simuler toutes les lois de la forme : ν(θi |θ1 , ..., θi−1 , θi+1 , ..., θk ).
Comme nous le verrons rapidement, ce cadre est particulièrement adapté à celui de l’inférence bayésienne.
L’algorithme est le suivant :- on se donne (θ10 , ..., θk0 ) artibraire.
- à l’étape n, on dispose donc de (θ1n , ..., θkn ). On simule alors :θ1n+1 ∼ ν(θ1 |θ2n , ..., θkn )θ2n+1 ∼ ν(θ2 |θ1n+1 , θ3n , ..., θkn )8θ3n+1 ∼ ν(θ3 |θ1n+1 , θ2n+1 , θ4n , ..., θkn )jusqu’à :n+1)
θkn+1 ∼ ν(θk |θ1n+1 , ..., θk−1

Les valeurs simulées par l’algorithme de Gibbs seront utilisées exactement de la même façon que celles
simulées par la méthode de Hastings-Metropolis. En particulier il faudra définir la burn-in period n1 , etn2 .
Il faut enfin maintenant vérifier que ν est bien une loi stationnaire pour la chaîne de Markov ainsi
simulée. Tout d’abord on note que le noyau de transition de cette chaîne s’écrit :′K(θ, dθ′ ) = ν(θ1′ |θ2 , ..., θk )ν(θ2′ |θ1′ , θ3 , ..., θk )...ν(θk′ |θ1′ , ..., θk−1)dθ′

Vérifions maintenant que si θn ∼ ν alors θn+1 ∼ ν, c’est-à-dire que pour tout A ensemble mesurable
on a :Z
P(θn+1 ∈ A) =ν(θn+1 )dθn+1
A

Or on a :P(θn+1 ∈ A) =
... =
or :

Z

Z

IA (θn+1 )K(θn , dθn+1 )ν(θn )dθnn+1
IA (θn+1 )ν(θ1n+1 |θ2n , ..., θkn )...ν(θkn+1 |θ1n+1 , ..., θk−1
)ν(θn )dθn+1 dθn

ν(θn ) = ν(θ1n |θ2n , ..., θkn )ν 1 (θ2n , ..., θkn )où ν j est le loi marginale de toutes les composantes sauf la j-ième.
On peut dès lors sortir de l’intégrale en dθn dθn+1 l’intégrale en dθ1n qui vaut :Z
ν(θ1n |θ2n , ..., θkn )dθ1n = 1
dès lors on peut regrouper les facteurs :ν(θ1n+1 |θ2n , ..., θkn )ν 1 (θ2n , ..., θkn ) = ν(θ1n+1 , θ2n , ..., θkn )On est alors à :P(θn+1 ∈ A) =Z

IA (θn+1 )ν(θ2n+1 |θ1n+1 , θ3n , ..., θkn )...

n+1
)ν(θ1n+1 , θ2n , ..., θkn )dθn+1 d(θ2n , ..., θkn )...ν(θkn+1 |θ1n+1 , ..., θk−1

On recommence la même manipulation :ν(θ1n+1 , θ2n , ..., θkn ) = ν(θ2n |θ1n+1 , θ3n , ..., θkn )ν 2 (θ1n+1 , θ3n , ..., θkn )ce qui permet d’intégrer en θ2n , et ainsi de suite jusqu’à intégrer en θkn . On a alors :ZZ
ν(θn+1 )dθn+1
P(θn+1 ∈ A) = IA (θn+1 )ν(θ1n+1 , ..., θkn+1 )dθn+1 =A

ce qu’il fallait démontrer.

91.2.4

Utilisation en pratique

Notons tout d’abord qu’on peut montrer que l’algorithme de Gibbs peut se mettre sous la forme d’uncas particulier de l’algorithme de Hastings-Metropolis dans lequel les probabilités d’acceptation seraient
toujours égales à 1.
Cependant, on voit bien que ces deux algorithmes ne présentent pas les mêmes avantages pratiquentet seront utilisés dans des contextes différents. En particulier, l’algorithme de Hastings-Métropolis nen’impose que de savoir simuler selon q(.|.) qui est choisie par l’utilisateur. On va voir dans la suite quel’algorithme de Gibbs est particulièrement adapté au cas de l’analyse bayésienne des modèles de mélanges
à k fixe.
Notons aussi que certains algorithmes dits hybrides peuvent être utilisés. En effet, chaque étape del’algorithme de Gibbs :n+1 n, θi+1 , ..., θkn )θin+1 ∼ ν(θi |θ1n+1 , θ2n+1 , ..., θi−1

est en fait un noyau de transition pour lequel la loi cible ν est invariante. Si on ne sait pas simuler la loi
n+1 n, θi+1 , ..., θkn ) dans une de ces étapes on peut remplacer cette étape par l’utilisation
ν(θi |θ1n+1 , θ2n+1 , ..., θi−1d’une probabilité de transition annexe, q(.|.), et utiliser une méthode d’acceptation-rejet comme dansl’algorithme de Hastings-Metropolis. La loi ν sera toujours invariante pour la chaîne simulée. Il fautensuite vérifier au cas par cas que la convergence a bien lieu (en vérifiant l’irréductibilité, l’apériodicité,...). C’est une méthode de ce type qui sera utilisé pour l’analyse bayésienne des mélanges de distributions
lorsque k est inconnu.

1.3

Inférence bayésienne dans les modèles de mélanges

1.3.1

Les modèles de mélanges comme modèles à données manquantes

Afin de pouvoir appliquer les algorithmes présentés précédemment aux modèles de mélange, il est
nécessaire de modifier encore la façon dont est présentée le modèle :g(x|ξ) =kXpi f (x|θi )i=1

Pour cela on crée en fait une variable supplémentaire "artificielle", z, qui représente la composante dumélange dont est issu x. En fait, z=i où i est le groupe (ou la classe) de l’observation. On suppose doncque la variable étudiée est le couple (x, z) mais que le modèle est à "données manquantes", c’est-à-direque la partie z n’est pas observée. On a alors :P(z = i|p1 , ..., pk ) = pig(x|ξ, z) = f (x|θz )Selon le point de vue bayésien, les zj peuvent aussi être vus comme des paramètres et le modèle
devient alors un modèle hiérarchique, qui s’écrit :g(x|ξ, z) = f (x|θz )∀i ∈ [1, k], P(z = i|ξ) = pipour lequel il faut choisir une loi a priori π(ξ) = π(θ, p).
Pour un échantillon de n observations x1 , ..., xn on a donc :(x1 , ..., xn )|(z1 , ..., zn , ξ) ∼kYY

i=1 j/zj =if (xj |θi )Notons enfin que selon les auteurs les conventions sont légèrement différentes pour zj : parfois, zj =(0, ..., 0, 1, 0, ..., 0) avec le 1 en position i où i est la classe de xj .
10

1.3.2

Utilisation de l’algorithme de Gibbs

L’utilisation de l’algorithme de Gibbs nécessite, comme on vient de le voir, de pouvoir simuler certaines
lois. Ceci impose, dans notre cadre, de choisir des lois conjuguées pour les paramètres θi et pi . Les lois
utilisées ici sont celles proposées dans Robert (1996) et Diebolt et Robert (1994).
Si les f (.|θi ) appartiennent à une même famille exponentielle :f (x|θi ) = h(x)eθi x−ψ(θ)
alors une loi conjuguée possible pour θi est :π(θi |α, β) = Ceθi α−βψ(θi )et il est possible de trouver une loi conjugée pour les pi , une loi dite de Dirichlet D(α1 , ..., αk ) dedensité :π(p1 , ..., pk |α1 , ..., αk ) =kkXΓ(α1 + ... + αk ) Y αi −1
pi = 1)(pi
I(pi ≥ 0))I(
Γ(α1 )...Γ(αk ) i=1
i=1

par rapport à la mesure uniforme sur {(p1 , ..., pk ) ∈ [0, 1]k /p1 + ... + pk = 1}.
Rappelons que le modèle est sous la forme :P(zj = i) = pi

(x1 , ..., xn )|(z1 , ..., zn , ξ) ∼kYY

i=1 j/zj =if (xj |θi )où ξ = (θ1 , ..., θk , p1 , ..., pk ). Or :Yj/zj =i⇒Yf (xj |θi ) = 

j/zj =iYj/zj =iθh(xj ) e i

f (xj |θi ) = Ki (x)π(θj |αi +Pj/zj =i xj −ψ(θi )card({j/zj =i})X

xj , βi + card({j/zj = i}))j/zj =iP
En notant alors ni = card({j/zj = i}) et xi = n1i j/zj =i xj ceci se réécrit simplement :Yf (xj |θi ) = Ki (x)π(θi |αi + ni xi , βi + ni )j/zj =i

Donc :(x1 , ..., xn )|(z1 , ..., zn , ξ) ∼ K(x)

k
Y

i=1

π(θi |αi + ni xi , βi + ni )et finalement on obtient la loi a posteriori des paramètres sous la forme :
!
! k
k
Y
Y
ξ|(x1 , ..., xn , z1 , ..., zn ) ∼π(θi |αi + ni xi , βi + ni )piαi +nii=1

i=1

En particulier les pi et les θi , sachant les xj et les zj , sont indépendants et on trouve les lois nécessaires
pour la mise en place de l’algorithme de Gibbs :∀i, θi |(p1 , ..., pk , x1 , ..., xn , z1 , ..., zn , (θi′ , i′ 6= i)) ∼ π(θi |αi + ni xi , βi + ni )(p1 , ..., pk )|(x1 , ..., xn , z1 , ..., zn , θ1 , ..., θk )) ∼ D(α1 + n1 , ..., αk + nk )11

pi f (xj |θi )∀i, ∀j, P(zj = i|p1 , ..., pk , x1 , ..., xn , θ1 , ..., θk ) = Pk
t=1 pt f (xj |θt )Finalement, l’algorithme de Gibbs s’applique de la façon suivante. On choisit une valeur quelconqueξ 0 et des valeurs zj0 , et on calcule les valeurs initiales de ni et xi , notées n0i et x0i . Le passage de ξ n et desvaleurs zjn à ξ n+1 et zjn+1 se fait de la façon suivante :- on simule :∀i, θin+1 ∼ π(θi |αi + nni xni , βi + nni )- on simule :) ∼ D(α1 + nn1 , ..., αk + nnk ), ..., pn+1
(pn+1
1
k

- on simule : ∀j, zjn suivant la loi :- on met à jour les valeurs :pn+1 f (xj |θn+1 )P(zj = i) = Pk i n+1 i n+1
f (xj |θt )t=1 ptnn+1= card({j/zjn+1 = i})
i
=
xn+1
i

X

1

xj
nn+1
i
j/z n+1 =i
j

1.3.3

Utilisation de l’algorithme de Hastings-Metropolis

Il serait aussi possible d’utiliser la méthode de Hastings-Metropolis dans le cas des modèles de mélanges. En fait, dans la mesure où cet algorithme ne nécessite pas de savoir simuler la loi d’un paramètre
conditionnellement à tous les autres, l’utilisation de cet algorithme impose moins de restrictions sur leslois a priori.
Rappelons la formulation du problème sous forme hiérarchique :g(x|ξ, z) = f (x|θz )Notons ζ = (ξ, z), la loi a priori sur ζ, notée p(ζ) est donnée par :∀i ∈ [1, k], P(z = i|ξ) = piet la loi a priori π(ξ) = π(θ, p).
On peut faire les mêmes choix que dans la partie sur l’algorithme de Gibbs pour p(.), mais ceci n’est
pas une obligation. Rappelons cependant qu’il ne faut pas utiliser de lois a priori impropres.
Alors :
(x1 , ..., xn )|ζ ∼kYY

i=1 j:zj =if (xj |θi )p(ζ|x) ∝ p(ζ)p(x|ζ)
p(ζ|x) ∝ π(ξ) 

n
Y

j=1



pj 
12

k
Y
Y

i=1 j:zj =if (xj |θi )p(ζ|x) ∝ π(ξ)

kYY

i=1 j:zj =i

pj f (xj |θi )Un dernier choix à faire est celui de q(ζ ′ |ζ) ; en général q choisit d’abord un des paramètres
selonP
une certaine loi de probabilités puis propose une modification de ce paramètre. Ainsi, q = m∈M qm , où
qm ne modifiera qu’un seul paramètre. Notons que dans ce cas, qm n’est pas vraiment la densité d’une
probabilité de transition car en général, ∀ζ, qm (Z|ζ) < 1 où Z est l’ensemble où le vecteur de paramètres
ζ varie. En fait, afin que q soit une probabilité de transition on a :X∀ζ,qm (Z|ζ) = 1m∈M

voire ≤ 1 si on garde la possibilité de ne pas proposer de modification des paramètres à chaque étape.
Ici, par exemple, on pourrait prendre q1 (.|ζ) qui propose une modification de θ, q2 (.|ζ) qui propose une
modification de z et q3 (.|ζ) pour p, avec M = {1, 2, 3}.
Reste alors à calculer la probabilité d’acceptation d’une modification de ζ en ζ ′ :


k
Y
Y
p′j f (xj |θi′ )  q(ζ|ζ ′ )π(ξ ′ )
′
i=1 j:zj =i


′


a(ζ, ζ ) = min 1,

k
Y Y


 q(ζ ′ |ζ)π(ξ)
pj f (xj |θi ) 
i=1 j:zj =i

1.3.4

Problèmes rencontrés

Un premier problème apparaît rapidement lorsque l’on analyse les résultats, c’est-à-dire lorsqu’on
visualise les lois a posteriori des paramètres : celles-ci sont multimodales, ce qui est du au fait que lesparamètres de différentes composantes ont pu être échangés au cours des simulations ! En effet, il n’y apas de "paramétrage naturel" des modèles de mélanges : par exemple, si il y a deux composantes (k = 2),
laquelle devrait être rattachée à p1 et θ1 et laquelle à p1 et θ1 ? Pour remédier à ce problème, différentes
solutions ont été proposées, la plus simple étant d’imposer une contrainte d’ordre sur les θi . Par exemple,
si les différentes composantes sont toutes des lois normales, avec θi = (µi , σi2 ) alors on peut imposerµ1 < ... < µk . En fait, cette méthode ne résout pas complètement le problème. D’autres propositions sontétudiées dans la partie suivante.
Un autre problème important reste celui de la convergence des méthodes par MCMC. Afin que laconvergence soit assurée pour les quantités d’intérêt, il faut s’assurer que la chaîne de Markov simulée dansl’espace des paramètres a bien visité les zones de fortes probabilités selon la loi a posteriori. Notamment,
la façon dont on choisit les valeurs initiales est très importante. Si on les choisit mal, la période dedémarrage ou burn-in period peut-être bien plus longue que la taille de la chaîne que l’on a simulé, ce
qui signifie que les résultats obtenus n’auront aucun sens. Certains auteurs ont proposé d’initialiser lesalgorithmes en utilisant un algorithme de classification classique et rapide, par exemple la méthode des(k) centres mobiles.

2

Détermination du nombre de composantes par MCMC à sautréversible

Comme on l’a vu en introduction, dans certains cas les composantes f (.|θi ) ont une signification"physique", et représentent des "groupes" d’individus : k peut alors être connu à l’avance, mais ce n’est
pas toujours le cas ! Il devient alors intéressant de considérer k comme un paramètre à estimer. Il se peutmême que k soit pratiquement le seul paramètre d’intérêt du modèle !

13

2.1

Introduction

2.1.1

Formulation du problème

On s’intéresse donc maintenant à l’étude de modèles de mélanges pour lesquels le nombre de composantes est inconnu. Le modèle reste :kXpi f (x|θi )i=1

mais le rôle de k est transformé : k n’est plus une constante, mais un paramètre qu’il faudra lui aussi
estimer. On a donc maintenant :g(x|ξ, k) =kXpi f (x|θi )i=1

Dans une approche bayésienne, il faudra donc choisir une loi a priori sur le paramètre k. En général,celle-ci est une loi de Poisson de paramètre λ, ou parfois une loi uniforme sur un intervalle d’entiers,[1, 100] par exemple.
Le modèle proposé est donc un modèle hiérarchique. En effet :ξ|k ∼ π(.|k)
où π(.|k) est une loi sur l’espace dans lequel varie ξ = (θ1 , ..., θk , p1 , ..., pk ). On voit bien que les lois
π(.|k = 2) et π(.|k = 3) par exemple ne peuvent pas être les mêmes car elles ne sont pas définies sur desespaces de même dimension. Cette formulation sera rendue plus explicite dans la suite où on présente laméthode de Green (1997) pour aborder ce problème.
2.1.2

Objectif et problème principal

Il suffit de connaître la loi de (k, ξ)|x ou mieux les lois de k|x puis de ξ|k, x. En particulier la loi de
k|x sera cruciale pour le choix du nombre de composantes.
Le problème est alors le même que dans le cas de l’analyse bayésienne des modèles de mélanges à
k fixé : ces lois sont donnés par une formule pas exploitable en pratique, et il va falloir recourir à desméthodes par MCMC.

2.2

Les algorithmes MCMC à sauts réversibles

2.2.1

Présentation générale

Le problème rencontré ici reste donc le même que pour les modèles de mélanges à k fixé : la loi aposteriori en général pas être connue explicitement. Il faut donc utiliser des méthodes MCMC afin desimuler une chaîne de Markov de loi stationnaire cette loi a posteriori.
Seulement, cette loi est définie sur un espace particulier. Si Θk est l’espace du paramètre ξk =(θ1 , ..., θk , p1 , ..., pk ), alors la loi a posteriori sur t = (k, ξ) est définie sur :[[[Θ = [{1} × Θ1 ] [{2} × Θ2 ] [{3} × Θ3 ] ...
Ainsi, la chaîne de Markov simulée doit être capable de "sauter" d’un espace Θk dans un autre, alorsmême que ces espaces ne sont pas de même dimension.
Le premier à avoir réussi à mettre cette méthode en pratique est Peter Green, en adaptant l’algorithme
de Hastings-Metropolis à des sauts d’un espace Θk à un autre. Ainsi, les sauts peuvent être accepté, ou
rejetés, comme n’importe quel mouvement dans l’algorithme de Hastings-Metropolis. Afin d’assurer laconvergence de la méthode, les sauts doivent satisfaire une condition de réversibilité ce qui a conduit Greenà appeler la méthode "MCMC à sauts réversibles", ou reversible jump MCMC, abrégée en RJMCMC. Il a
tout d’abord présenté l’idée dans la discussion d’un article de Grenander et Miller (1994) et l’a développé
dans Reversible jump Markov chain Monte Carlo computation and Bayesian model determination (1995).
14

Dans cet article, cette méthode est utilisée dans différents types de problèmes bayésiens dans lesquels
la dimension du paramètre n’est pas connue a priori : choix des regresseurs dans une regression, segmentation d’image, évaluation de la tendance d’une série avec un nombre inconnu de changements detendance. Le problème de l’estimation d’un modèle de mélanges avec un nombre inconnu de composantesétait seulement cité comme exemple.
Richardson et Green, dans On Bayesian analysis of mixtures with an unknown number of components(1997), ont consacré un article complet à l’application de cette méthode aux modèles de mélanges. La
suite de cette partie est consacrée à la description de leur méthode et à la présentation de leurs résultats.
2.2.2

Principe des RJMCMC

Cette section présente la méthode des RJMCMC telle qu’elle est exposée par Green (1995).
L’algorithme proposé est donc une adaptation de celui de Hastings-Metropolis. Notons π(t) poursimplifier la loi a posteriori du paramètre t (qui contient k, les θi et les pi ), c’est-à-dire que l’on omet denoter le "sachant l’observation x" dans toute cette sous-partie 2.2.2. A chaque étape, différents types demodification des paramètres sont proposées, puis rejetées ou acceptées :- mise à jour des poids pi ,- mise à jour des paramètres θi ,- mise à jour des zj ,- changement du nombre de composantes k.
C’est ce dernier type de modification qui pose un problème. En fait, plusieurs changements peuventêtre proposés : augmenter ou diminuer k de 1, de 2...
Concentrons-nous pour le moment sur la probabilité de transition q dans ce dernier cas uniquement(modification de k). On note l’ensemble de toutes les modifications possibles M , et comme on l’a fait
précédemment on écrit :X
q(A|t) =qm (A|t)m∈Moù qm (A|t) est la probabilité de proposer l’ensemble A en partant de l’état t si on tente une modification
de type m. On rappelle que :X
qm (Θ|t) = 1m∈Mou ≤ 1 si on garde la possibilité de ne pas proposer de modifications à certaines étapes.
Notons que pour un état t donné il est toujours possible qu’il y ait un type de modification m telque : qm (Θ, t) = 0. C’est même toujours le cas en pratique : en effet, m peut être une modification desparamètres θ ou des poids p mais m peut aussi impliquer un changement du nombre de composantes k :par exemple 1 ↔ 2, 1 ↔ 3, 2 ↔ 3... Il est logique que si k = 1, la modification m = 2 ↔ 3 ne puisse pasêtre proposée...
Le but est maintenant de calculer la probabilité d’acceptation am (t, t′ ) qu’il faut associer à qm sil’on souhaite simuler une chaîne de Markov de probabilité stationnaire π. Pour cela, on peut se rappeler
que la detailed balance condition, si elle n’est pas une condition nécessaire, est une condition suffisantepour assurer ceci. Rappelons que dans un algorithme de Hastings-Metropolis simple, la probabilité detransition de la chaîne de Markov simulée était :P (t, B) = a(t, B)q(B|t) + (1 − a(t, B))IB (t)elle pourra s’écrire dans ce cas :P (t, B) =

X Z

m∈M

+IB (t)

X Z

m∈M

B

qm (dt′ |t)am (t, t′ )

′

Θ

′

qm (dt |t)[1 − am (t, t )]
15

!

!

+ IB (t)(1 −Xm∈M

qm (Θ|t))Le premier terme correspond à l’acceptation de la modification proposée par q, le second terme à son
rejet et le troisième terme au cas hypothétique où aucune modification n’aurait été proposée.
Pour simplifier, Green note :!
X
X Zs(t) = 1 −qm (Θ|t) +qm (dt′ |t)[1 − am (t, t′ )]
m∈M

Θ

m∈M

donc :X Z

P (t, B) =

m∈M

′

B

′

qm (dt |t)am (t, t )!

+ IB (t)s(t)La detailed balance condition sera assurée si, ∀A, ∀B mesurables dans Θ on a :ZZπ(dt)P (t, B) =π(dt)P (t, A)A

B

soit :

XZ
m

=

π(dt)
A

XZ
m

Z

π(dt′ )
B

′

B

Z

′

qm (dt |t)am (t, t ) +Z

qm (dt|t′ )am (t′ , t) +A

π(dt)s(t)A∩BZπ(dt′ )s(t′ )B∩ALes deux intégrales simples sont évidemment les mêmes. Donc une condition suffisante, et simple,
pour assurer la condition désirée est la réversibilité de chaque saut proposé, soit ∀m ∈ M, ∀A, ∀B :
Z
Z
Z
Z
′
′
′
π(dt)
qm (dt |t)am (t, t ) =π(dt )
qm (dt|t′ )am (t′ , t)
A

B

B

A

Or, supposons que la mesure π(dt)qm (dt′ |t) ait une densité fm (t, t′ ) par rapport à une mesure νmsymétrique sur l’espace produit Θ2 , on a alors :
Z
Z
Z Z
π(dt)
qm (dt′ |t)am (t, t′ ) =νm (dt, dt′ )fm (t, t′ )am (t, t′ )
A

B

A

B

et symétriquement :ZZ
Z Zπ(dt′ )
qm (dt|t′ )am (t′ , t) =νm (dt′ , dt)fm (t′ , t)am (t′ , t)
B

A

B

A

ainsi, la condition d’égalité s’écrit est assurée si am vérifie :∀m ∈ M, ∀(t, t′ ) ∈ Θ2 : am (t, t′ )fm (t, t′ ) = am (t′ , t)fm (t′ , t)
Pour ceci il suffit de prendre :
fm (t′ , t)am (t, t′ ) = min 1,
fm (t, t′ )que l’on peut réécrire formellement :π(dt′ )qm (dt|t′ )am (t, t′ ) = min 1,
π(dt)qm (dt′ |t)

En fait, cette forme peut être rendue plus explicite à condition de préciser la forme des qm .
On se place dans le cas où qm propose un saut de l’espace [{k1 } × Θk1 ] de dimension d1 à l’espace[{k2 } × Θk2 ] de dimension d2 > d1 , ou le saut de [{k2 } × Θk2 ] dans [{k1 } × Θk1 ]. Rappelons que t = (k, ξ)et que dans le couple (k1 , ξ) le ξ est de dimension d1 alors qu’il est de dimension d2 dans dans le couple
16

(k2 , ξ). On supposera que π(dt|k = k1 ) admet une densité fk (t) par rapport à la mesure de Lebesgue
dans Θk1 , on notera par commodité :π(dt) = π(k)π(dt|k) = π(k)fk (t)dt = π(t)dt
On note j((k1 , ξ1 )) la probabilité de choisir le type de saut m lorsqu’on est en (k1 , ξ1 ), et la probabilitéde choisir le saut m (en sens inverse) lorsqu’on est en (k2 , ξ2 ) sera notée j((k2 , ξ2 )). En fait :Zj((k1 , ξ1 )) =qm (dt′ |(k1 , ξ1 ))Θj((k2 , ξ2 )) =Z

qm (dt|(k2 , ξ2 ))Θ

Il faut maintenant expliciter la façon dont sera réalisé un saut. On l’a vu, intuitivement le problème
est que la dimension des deux espaces entre lesquels on saute n’est pas la même. La façon d’accomplir cesaut proposée par Green est la suivante :- on génère un vecteur aléatoire de dimension d = d2 − d1 , noté u(1) , selon une loi de densité q1 (u(1) )par rapport à la mesure de Lebesgue.
- on prend ensuite t′ = g(t, u(1) ) où g est une fonction déterministe. g doit vérifier :∀t, ∀t′ , ∃!u(1) : g(t, u(1) ) = t′∃ginv , ∀t, ∀u(1) : ginv (g(t, u(1) )) = t
Autrement dit, pour un saut de Θk1 dans Θk2 il faut "compléter" la dimension par l’introduction deu(1) qui est aléatoire et utiliser g, alors que pour un saut inverse de Θk2 dans Θk1 il faut utiliser ginv etla transformation proposée est déterministe.
Maintenant que qm est décrite il reste à déterminer la mesure symétrique νm et la densité fm . Greendéfinit νm par :∀A ⊂ {k1 } × Θk1 , ∀B ⊂ {k2 } × Θk2 mesurables :νm (A × B) = νm (B × A) = λ{(t, u(1) )/ξ1 ∈ A, g(t, u(1) ) ∈ B}

où λ est la mesure de Lebesgue sur Rd2 , et :∀A ⊂ Θ, ∀B ⊂ Θ mesurables :νm (A × B) = νm {(A ∩ Θk1 ) × (B ∩ Θk2 )} + νm {(A ∩ Θk2 ) × (B ∩ Θk1 )}Notons que νm est symétrique, ce qui était nécessaire. Alors on a :π(dt)qm (dt′ |t) = f (t, t′ )νm (dt, dt′ )avec f définie de la façon suivante : ∀t = (k1 , ξ1 ) ∈ Θk1 , ∀t′ = (k2 , ξ2 ) ∈ Θk2
f (t, t′ ) = π(t)j(t)q1 (u(1) )f (t′ , t) = π(t′ )j(t′ )∂g(t, u(1) )∂(t, u(1) )et f (x, x′ ) = 0 dans tous les autres cas. Dans ces cas on obtient donc une formule plus explicite pourla probabilité d’acceptation :∂g(t, u(1) )π(t′ )j(t′ )′a(t, t ) = min 1,
π(t)j(t)q1 (u(1) ) ∂(t, u(1) )La preuve de cette assertion n’est donnée explicitement ni dans Green (1995) ni dans Richardson etGreen (1997). En fait, ce résultat peut être vu de façon intuitive en utilisant la formule générique :17π(dt′ )qm (dt|t′ )am (t, t ) = min 1,
π(dt)qm (dt′ |t)′et la définition de qm .
Le résultat est néanmoins démontré rigoureusement à partir de la detailed balance condition dans
Waagepetersen et Sorensen (2000) ou dans Saint-Pierre et Garel (2001). Rappelons que cette conditions’écrivait :
Z
Z
Z
Z
π(dt)
qm (dt′ |t)am (t, t′ ) =π(dt′ )
qm (dt|t′ )am (t′ , t)
A

B

B

A

On a :

... =
... =
et :

Z

Z
π(t)dt
A

Z Z

π(dt)
A

Z

Z

qm (dt′ |t)am (t, t′ ) = ...

B

du(1) q1 (u(1) )j(t)am (t, g(t, u(1) ))IB (g(t, u(1) ))q1 (u(1) )j(t)am (t, g(t, u(1) ))IB (g(t, u(1) ))π(t)IA (t)du(1) dt
Z
... =

... =

Z

π(dt′ )B

Z

Z

qm (dt|t′ )am (t′ , t) = ...
Aπ(t′ )dt′
BZAj(t′ )am (t′ , t)dδginv (t′ ) (t)π(t′ )dt′ IB (t′ )j(t′ )am (t′ , ginv (t′ ))IA (ginv (t′ ))Il suffit alors d’obtenir le changement de variable t′ = g(t, u(1) ) pour obtenir :
... =

Z Z

π(g(t, u(1) ))IA (t)IB (g(t, u(1) ))j(g(t, u(1) ))am (g(t, u(1) ), t)∂g(t, u(1) )du(1) dt∂(t, u(1) )La relation d’égalité entre les deux résultats est effectivement satisfaite par la solution donnée par
Green :π(g(t, u(1) ))j(g(t, u(1) )) ∂g(t, u(1) )a(t, t′ ) = min 1,
π(t)j(t)q1 (u(1) )∂(t, u(1) )2.3

La méthode retenue par Richardson et Green (1997)

C’est donc dans On Bayesian analysis of mixtures with an unknown number of components (1997)
que Richardson et Green ont appliqué la première fois cette idée de RJMCMC pour estimer le nombre
de composantes dans un mélange de distributions. Cette sous-partie présente les choix qui ont été faits
par les auteurs dans cet article pour les points suivants :- le modèle estimé.
- les lois a priori.
- les types de mouvement (m ∈ M ) et les transitions associées (qm ), ainsi que les calculs des probabilitésd’acceptation associées (am ).
Notons que l’aspect théorique des RJMCMC ayant été présenté dans Green (1995), l’article de Richardson et Green détaille simplement une application de cette méthode pour un modèle de mélanges
particulier, en précisant cependant que la méthode peut être adaptée sans difficulté à d’autres modèles.
Il s’agit ici d’un mélange de lois normales unidimensionnelles.
18

Pour ce qui est du choix de M et des qm , les auteurs précisent que les choix ont été faits de façon
arbitraire, ou du moins qu’aucune optimisation des performances de l’algorithme n’a été menée à travers
ce choix : les résultats obtenus sont néanmoins satisfaisants. La recherche de qm "optimaux" a suscité
depuis la parution d’un certain nombre d’articles, comme ceci sera montré dans la suite.
2.3.1

Cadre général, notations

On adopte ici les notations retenues par Richardson et Green qui sont légèrement différentes. Les
observations sont notées yi pour i ∈ {1, ..., n}, et y = (y1 , ..., yn ). Le nombre de composantes est donc un
paramètre du modèle, toujours noté k et les paramètres de chaque composante sont toujours noté θj , j ∈{1, ..., k}, et θ = (θ1 , ..., θn ). Les variables d’allocations sont toujours notées zi ∈ {1, ..., k}, j ∈ {1, ..., n}.
Les poids sont maintenant notés wj pour j ∈ {1, ..., k}, et w = (w1 , ..., wk ).
yi |k, θ, w ∼

k
X
j=1

wj f (.|θj ) = f (yi |k, θ, w)∀j ∈ {1, ..., k}, P(zi = j|k, w) = wj
yi |k, θ, z ∼ f (.|θzi )Devant le nombre de lois considérées, une notation générique est adoptée : p(.). La loi jointe desparamètres et des variables et la suivante :p(k, w, z, θ, y) = p(k)p(w|k)p(z|w, k)p(θ|z, w, k)p(y|θ, z, w, k)p(z|w, k) a été complètement définie précédemment, il reste à préciser p(k), p(w|k), p(θ|z, w, k) =p(θ|k) et p(y|θ, z, w, k) = p(y|θ, z) c’est-à-dire f (.|θ).
En fait, un niveau hiérarchique supplémentaire est rajouté par les auteurs. Ceci permettra notammentde diminuer l’influence des lois a priori sur les résultats obtenus, ce point est discuté dans la présentation
des résultats, ou de les calibrer de façon adéquate. Les auteurs introduisent l’hyperparamètre λ qui sera
paramètre de la loi a priori de k, δ pour w et η pour θ. L’intérêt est le suivant : si les lois a posteriori
semblent trop sensible à des choix particuliers de η, par exemple, le niveau hiérarchique supérieur permetde munir η lui-même d’une loi a priori ce qui a pour effet de minimiser cette dépendance.
On obtient finalement la loi jointe des observations, paramètres et hyperparamètres :p(λ, δ, η, k, w, z, θ, y) = p(λ)p(δ)p(η)p(k|λ)p(w|k, δ)p(z|w, k)p(θ|k, η)p(y|θ, z)
On peut poser t = (λ, δ, η, k, w, z, θ) pour recoller aux notations des parties précédentes et alléger lesnotations ; on a :p(t|y) ∝ p(λ, δ, η, k, w, z, θ, y)
2.3.2

Spécification du modèle

Détaillons maintenant les choix retenus par les auteurs pour les dernières lois qui n’ont pas encore été
spécifiées.
Tout d’abord, il faut préciser la loi du paramètre k. Le choix retenu dans les autres applications des
RJMCMC par Green était une loi de Poisson de paramètre λ, notée ici P(λ). Cependant, afin de faciliter
la présentation, Green retient ici une loi uniforme sur {1, ..., kmax = λ} :k ∼ U ({1, ..., kmax })1
I{1,...,kmax } (k0 )kmax
Ensuite, comme on l’a mentionné précédemment, les observations suivent des lois normales unidimensionnelles :19
P(k = k0 ) =yi |θ, z ∼ N (µzi , σz2i )donc on peut détailler ici θj = (µj , σj2 ) ∈ R × R+ . Ainsi :f (yi |θ, z) = p−1e2πσz2i(yi −µz )2
i
2
2σz
i

Il faut donc maintenant préciser les lois a priori sur θj = (µj , σj2 ). On suppose en fait que tous les µj
et tous les σj2 sont indépendants et :µj |η, k ∼ N (ξ, κ−1 )σj2 |η, k ∼ IG(α, β)avec l’hyperparamètre η = (ξ, κ−1 , α, β), et IG est une loi gamma inverse c’est-à-dire que :σj−2 |η, k ∼ Γ(α, β)p(σj−2 |η, k) =−2
βα(σj−2 )α−1 e−βσj IR+ (σj−2 )Γ(α)

On obtient par le calcul :p(σj2 |η, k) =

− β2
βα
1
σ
−2
j I
e
R+ (σj )Γ(α) (σj2 )α+1

Une restriction cependant est à apporter : comme on l’a vu dans la première partie il est nécessaired’imposer des contraintes d’identifications sur les composantes. Ici, les contraintes retenues par les auteurssont :µ1 < ... < µk
Notons que, pour ce qui est des paramètres θ et p, le choix des lois a priori est assez semblable à celui
fait par Robert (1996) et présentée dans la première partie. La raison va être présentée assez rapidement :en fait, lors des simulations, Green et Richardson n’utilisent la méthode de Hastings-Metropolis que pourles changements de k. Afin de simplifier le problème, pour la modification des autres paramètres, c’est la
méthode de Gibbs qui a été retenue. Le choix de lois a priori fait ici permet d’appliquer l’algorithme deGibbs, ceci sera détaillé bien sûr par la suite.
Le choix retenu pour le vecteur des poids w = (w1 , ..., wk ) est une loi de Dirichlet symétrique deparamètre δ :w|δ, k ∼ D(δ, ..., δ)
déjà présentée dans la première partie.
Reste enfin à spécifier les lois des hyperparamètres. En fait, dans chaque application de la méthodeproposée par les auteurs, λ = kmax sera fixé à 30, ce qui n’aura aucune conséquence comme on le verra
par la suite.
Pour le choix des autres hyperparamètres, les auteurs précisent avoir voulu être "le moins informatif
possible", rappelons que l’utilisationde loi complètement non-informatives dans ce contexte est impossible.
Afin que la loi (a priori) des µj soit relativement plate sur l’intervalle où varient les données, les auteurs
prennent ξ égal au milieu de cet intervalle et κ un petit multiple de 1/R2 où R est la longueur de cetintervalle.
L’hyperparamètre δ est fixé, égal à 1 dans toutes les applications. Si on revient à la définition de laloi de Dirichlet (donnée dans la première partie), ceci signifie simplement que w est de loi uniforme surle simplexe de Rk .
20

Reste finalement à déterminer α et β. Après avoir mené de premières simulations à α et β fixés, lesauteurs se sont aperçus que la loi a posteriori obtenue sur les paramètres σj2 dépendait fortement desvaleurs de α et β.
Rappelons que :E(σj−2 |α, β, k) =αβV(σj−2 |α, β, k) =αβ2Afin de diminuer cette dépendance, les auteurs ont donc décidé d’introduire un niveau hiérarchique
supplémentaire en prenant en fait α fixe mais :β|g, h ∼ Γ(g, h)

Les auteurs remarquent qu’en fixant α > 1 > g on introduit une information a priori sur les σj2 enles contraignant à être du même ordre de grandeur, ce qui est en fait souhaitable si on souhaite avoir desrésultats cohérents. Les auteurs retiennent donc ce choix et prennent h égal à un petit multiple de 1/R2 .
Finalement, une fois introduit ce dernier niveau hiérarchique, les résultats obtenus qui seront présentés
plus loin se sont avérés assez peu sensibles aux choix particuliers faits pour g, h, α, κ.
Notons enfin que, pour toutes les mises en oeuvre sur données réelles que les auteurs font à la fin del’article, on a : κ = 1/R2 , α = 2, g = 0.2, h = 10/R2 , δ = 1 et kmax = 30.
2.3.3

Méthode de simulation : choix des types de mouvements

Les auteurs distinguent six types principaux de mouvements :(a) modification des poids w ;(b) modification des paramètres des lois normales, (µ, σ) ;(c) modification des variables d’allocation z ;(d) modification de l’hyperparamètre β ;(e) modification de k : séparer une composante en deux ou en réunir deux en une ;(f) modification de k : créer ou supprimer une composante vide.
Plutôt que de choisir un type de mouvement au hasard, les auteurs utilisent un algorithme qui passesystématiquement par les étapes (a), (b), (c), (d), (e) et (f) dans l’ordre à chaque étape ou sweep. En
fait, le choix de lois a priori qu’ils ont effectué leur permet de traiter les étapes (a) à (d) par la méthodeproposée par Diebolt et Robert (1994), c’est-à-dire par un algorithme de Gibbs (présenté en détail dansla première partie). En revanche, les étapes (e) et (f) ne peuvent être traitées que par un algorithme deHastings-Metropolis, c’est pour ces étapes que la méthode des RJMCMC développée par Green (1995)
se révèle utile.
Examinons donc tout d’abord les étapes (a) à (d) qui nécessitent de connaître la loi conditionnelle
de chaque paramètre, sachant tous les autres et sachant les observations. La loi jointe des paramètres,
hyperparamètres et observations était écrite sous la forme :p(t, y) = p(λ, δ, η, k, w, z, θ, y)p(t, y) = p(λ)p(δ)p(η)p(k|λ)p(w|k, δ)p(z|w, k)p(θ|k, η)p(y|θ, z)Comme on a retenu λ = kmax et δ = (g, h) constants on a p(λ) = p(δ) = 1, et on a η = (ξ, κ−1 , α, β)avec tous ces hyperparamètres constants sauf β qui suit une loi gamma. On obtient finalement :khg g−1 −hβ  Y δ−1 
1...
I{1≤k≤kmax }βewjp(λ, δ, η, k, w, z, θ, y) =
kmax
Γ(g)
j=1




−(µj −ξ)2
k
k
α
Y
Y
−1
2κ
−2
β
e
n
 ...
√
wj j  (σj−2 )α−1 e−βσj  
... 
−1
Γ(α)
2πκ
j=1
j=1
j=1


k
Y

21




−(yi −µz )2
i
2
n
2σz
Y
i
e


p
... 
 I{β>0,w1 >0,...,wk >0,w1 +...+wk =1,κ−1 >0,σ12 >0,...,σk2 >0}
2
2πσ
zi
i=1

où nj = card({i = 1, ..., n : zi = j}). Comme on impose en plus la contrainte µ1 < ... < µk , il faudrait
encore en principe multiplier cette expression par Iµ1 <...<µk k!
L’étape (a) consiste donc à modifier les paramètres de poids w. Le calcul de la loi conditionnellesachant tous les autres paramètres et les observations est direct à partir de la formule ci-dessus (on omet
les indicatrices pour simplifier l’écriture) :



k
k
Y
Y
n
wj j 
wjδ−1  p(w|...) ∝ 
j=1

j=1

p(w|...) ∝kY

δ−1+njwj

j=1⇒ w|... ∼ D(δ + n1 , ..., δ + nk )L’étape (b) consiste à modifier les paramètres des lois normales (µ, σ). On a d’après la formule de laloi jointe, pour le paramètre µ :!
−(yi −µz )2
n
k
i
Y
Y
−(µj −ξ)222σziee 2κ−1 p(µ|...) ∝ 
i=1

j=1


k
2
X (yi − µj )2
X
(µ
−
ξ)
1
j


+
p(µ|...) ∝ exp −
2 j=1
κ−1
2σj2
i:z =j




i





 
k
 1X
X
X
ξ
1
yi 
µ2j  1 +− 2µj  −1 +p(µ|...) ∝ exp −
2
−1
 2

κ
σ
κ
σ2
i:z =j j
i:z =j j
j=1
i



i



X

σ −2
yi + κξ j i/z =j

1
i


⇒ µj |... ∼ N 
,

−2
−2

σj nj + κσj nj + κ notée N (m, s2 ), sous la contrainte µ1 < ... < µk , ce qui signifie en fait que la loi jointe des µj s’écrit :Iµ1 <...<µk k!

k
Y

j=1

√



1
exp − 2 (µj − m)2
2s
2πs2
1

et pour le paramètre σ :−(yi −µz )2
i
2
n
2σz
Y
i
−2
e


p
Γ(α)(σj−2 )α−1 e−βσj  p(σ|...) ∝ 

2
2πσzi
i=1
j=1


p(σ|...) ∝

k
Y

j=1

k
Y



(σj−2 )α−1  q

1
2πσj2

nj




exp −βσj2 − σj222

X (yi − µj )2

i:zi =j2

p(σ|...) ∝

k
Y

j=1

(σj−2 )

nα−1− 2j

exp −σj−2 β +1 X

2 i:z =ji(yi − µj )2 

1
1 X⇒ σj−2 |... ∼ Γ(α + nj , β +(yi − µj )2 )2
2i/zi =jL’étape (c) consiste à modifier les variables d’allocation z. On a :wj −e⇒ P(zi = j|...) ∝σj(yi −µj )22σ 2
jEnfin l’étape (d) consiste à modifier l’hyperparamètre β. On a :kX
σj−2 
p(β|...) ∝ β g+kα exp −hβ − β
j=1⇒ β|... ∼ Γ(g + kα, h +kX

σj−2 )

j=1

Examinons maintenant les étapes (e) et (f). Ces étapes sont les deux qui font intervenir la méthodedes RJMCMC développée par Green (1995) : il faut ici détailler quelles seront les types de mouvements
qm , m ∈ M et calculer les probabilités d’acceptation. On peut se demander pourquoi les auteurs ont
retenu les deux étapes e et f : une seule de ces deux étapes aurait permis à l’algorithme de modifier leparamètre k. Nous reviendrons sur ce choix dans la suite car depuis, différents auteurs ont proposé desaméliorations de la méthode qui jouent sur la présence ou non de ces deux types d’étapes e et f .
L’étape (e) consiste à séparer une composante en deux ou à en regrouper deux en une seule. Donccette étape ne peut modifier k qu’en lui ajoutant ou lui soustrayant 1. Les types de mouvements possibles
seront donc :M = {1 ↔ 2, 2 ↔ 3, ..., kmax − 1 ↔ kmax }Lorsqu’on est en k, seuls les mouvements possibles seront donc k − 1 ↔ k et k ↔ k + 1, choisis
respectivement avec une probabilité dk (d pour death, "mort" d’une composante) et bk (b pour birth,"naissance" d’une composante). Bien sûr comme k ∈ {1, ..., kmax } on a bkmax = d1 = 0. Le choix des autres
dk et bk est arbitraire, les auteurs choisissent : dkmax = b1 = 1 et pour k ∈ {2, ..., kmax − 1}, bk = dk = 12 .
Détaillons tout d’abord comment deux composantes sont regroupées (si en k + 1 on a choisi la modification k ↔ k + 1). Deux composantes j1 et j2 sont choisies au hasard (de façon uniforme) telles queµj1 < µj2 et telles qu’il n’y ait pas d’autre µj dans l’intervalle [µj1 , µj2 ]. Ces deux composantes seront
regroupées, en une nouvelle composantes que l’on notera j ∗ .
Tout d’abord toutes les variables d’allocation sont modifiées : tous les zi valant j1 ou j2 vaudront
maintenant j ∗ . Ceci induit un changement d’indexation des composantes : par exemple si j1 = 4 et
j2 = 5 alors j ∗ = 4 et il faut renommer la composante 6 composante 5, la 7 devient 6 et ainsi de suite
jusqu’à la k + 1 qui devient k.
Les auteurs imposent ensuite que wj ∗ , µj ∗ et σj2∗ vérifient :w j ∗ = w j1 + w j2
w j ∗ µ j ∗ = w j1 µ j1 + w j2 µ j2
wj ∗ (µ2j ∗ + σj2∗ ) = wj1 (µ2j1 + σj21 ) + wj2 (µ2j2 + σj22 )ce qui définit complètement la nouvelle composante.
23

Le mouvement inverse, c’est-à-dire le passage de k à k + 1 doit vérifier lui aussi ces équations. Cettefois, comme la dimension du paramètre va augmenter, il est nécessaire de simuler des variables aléatoires
comme cela a été décrit dans la présentation générale des RJMCMC. En fait, la dimension augmente de
3 il faut donc simuler trois variables aléatoires indépendantes :u1 ∼ Be(2, 2), u2 ∼ Be(2, 2), u3 ∼ Be(1, 1)
où Be(p, q) avec p > 0, q > 0 dénote une loi Beta de densité :be(x|p, q) =où B(p, q) est la fonction Beta :xp−1 (1 − x)q−1
I]0,1[ (x)B(p, q)B(p, q) =Γ(p)Γ(q)Γ(p + q)Notons u = (u1 , u2 , u3 ). La loi de probabilité q suivant laquelle on génère u, nécessaire au calcul
des probabilités d’acceptation dans le cadre des RJMCMC, est donc complètement déterminée. On a sa
densité :q(u) = be(u1 |2, 2)be(u2 |2, 2)be(u3 |1, 1)

On choisit cette fois une composante j ∗ au hasard qu’on va séparer en deux : j1 et j2 en réindexant
les composantes suivantes. On pose :w j1 = u 1 w j ∗wj2 = (1 − u1 )wj ∗ainsi la contrainte wj ∗ = wj1 + wj2 est respectée, puis :r
w j2µ j1 = µ j ∗ − u 2 σ j ∗
w j1
r
w j1
µ j2 = µ j ∗ + u 2 σ j ∗w j2
ainsi la contrainte wj ∗ µj ∗ = wj1 µj1 + wj2 µj2 est respectée, et enfin :σj21 = u3 (1 − u22 )σj2∗

wj ∗
w j1

σj22 = (1 − u3 )(1 − u22 )σj2∗

wj ∗
w j2

ainsi la contrainte wj ∗ (µ2j ∗ + σj2∗ ) = wj1 (µ2j1 + σj21 ) + wj2 (µ2j2 + σj22 ) est respectée.
Ainsi, la fonction g qui permet de passer des paramètres t et du vecteur aléatoire u dans le modèle à
k composantes au paramètre t′ = g(t, u) dans le modèle à k + 1 composantes est entièrement déterminéepar les équations précédentes.
La dernière modification à faire pour achever la création des deux composantes est d’actualiser lesvariables d’allocation zi pour les i tels que zi = j ∗ . Elles sont tirées au hasard suivant la loi :wj −eP(zi = j|...) ∝σj(yi −µj )2
2σ 2
jpour j ∈ {j1 , j2 }.
Ainsi toutes les modifications possibles dans l’étape (e) ont été définies, reste encore à calculer les
probabilités d’acceptation. En fait, on a de façon générique (suivant les notations que l’on avait dans la24partie précédente, où am désignait la probabilité d’acceptation d’une modification, et qm est la loi selonlaquelle on propose la modification) :
p(dt′ )qm (dt|t′ )′am (t, t ) = min 1,
p(dt)qm (dt′ |t)
Donc on peut écrire : am (t, t′ ) = min(1, A), avec :A=
et donc :p(dt′ )qm (dt|t′ )p(dt)qm (dt′ |t)am (t′ , t) = min(1, A−1 )On a développé précédemment comment Green (1995) a obtenu pour les RJMCMC, dans le cas del’augmentation de la dimension du paramètre (donc ici la séparation de deux composantes) :A=p(t′ )j(t′ ) ∂g(t, u)p(t)j(t)q(u) ∂(t, u)où j(t) est la probabilité de proposer une modification de type m c’est-à-dire k ↔ k + 1 lorsqu’on est
en t, et est donc donnée par bk , de même j(t′ ) = dk
Ici on a :12
wjδ−1+lwjδ−1+lp(t′ )L′ p(k + 1)
1
2
...
=
(k + 1) δ−1+l
1 +l2
p(t)
L p(k)wj ∗B(δ, kδ)
r



κ
1
exp − κ (µj1 − ξ)2 + (µj2 − ξ)2 − (µj ∗ − ξ)2 ...
...
2π
2
!−α−1

σj21 σj22
βα− σj−2+ σj−2exp −β(σj−2
...
∗ )
2
2
1
Γ(α)
σj ∗

′

où l1 est le nombre d’observations allouées à j1 , l2 le nombre d’observations allouées à j2 et LL est le
ratio de la nouvelle vraisemblance sur l’ancienne, c’est-à-dire si on note respectivement zi et zi′ les variablesd’allocation dans les modèles à k + 1 et k composantes et de même (µi , σi ) (µ′i , σi′ ) les paramètres on a :−(yi −µ′z′ )2−(yi −µ′z′ )2
i
i
n exp
2σz′2′
2σz′2′
Y
Y exp
i
i
q
q
′2
′2
2πσz′
2πσz′
i=1
i:zi′ =j1 ,j2
L′
i
i



=
2 =
−(yi −µj ∗ )2−(yi −µzi )
L
n exp
n
exp
Y
Y
2σz2
2σj2∗
i
p
q
2
2πσ
zi
2πσj2∗
i=1
i:zi =j ∗et :j(t′ )
dk+1=
j(t)q(u)bk Palloc be(u1 |2, 2)be(u2 |2, 2)be(u3 |1, 1)

où Palloc est la probabilité de l’allocation entre j1 et j2 des variables regroupées sous j ∗ , c’est-à-dire
que :




(yi −µz′ )2
wz′
i
i
exp
−
2
σz ′
2σz′
Y 

i
i






Palloc = wj
wj2(yi −µj1 )2(yi −µj2 )2
1
i:zi =j ∗+ σj exp − 2σ2σj exp −
2σ 2
1

j1

25

2

j2

et :
wj ∗ |µj1 − µj2 |σj21 σj22∂g(ξ2 )=
u2 (1 − u22 )u3 (1 − u3 )σj ∗∂(ξ1 , u(1) )On obtient donc une expression exacte de A. Enfin, pour le mouvement en sens inverse (regroupement
de deux composantes) la probabilité d’acceptation est alors :min(1, A−1 )L’étape (f) enfin consiste à créer ou supprimer une composante vide. Ceci est en fait plus simple àdécrire. Ici, l’ensemble des modifications possibles est choisi similaire au précédent :M = {1 ↔ 2, 2 ↔ 3, ..., kmax − 1 ↔ kmax }En k, le mouvement k ↔ k + 1 (création d’une composante vide) est choisi avec la probabilité bk et
le mouvement k − 1 ↔ k (suppression d’une composante vide) avec la probabilité dk , les bk et dk étant
les mêmes que dans l’étape e.
Pour la création d’une composante, on introduit donc j ∗ avec wj ∗ ∼ Be(1, k), µj ∗ ∼ N (ξ, κ−1 ) et−2
σj ∗ ∼ Γ(α, β). Enfin, toutes les variables de poids doivent être renormalisées de façon à ce que leursomme soit 1, c’est-à-dire qu’on remplace wj par wj (1 − wj ∗ ) pour tout j 6= j ∗ .
Pour le mouvement inverse, on choisit au hasard une composante vide j ∗ (si il y en a une) de façonw
uniforme et on l’efface, il suffit de renormaliser les poids en remplaçant wj par P ′ j∗ w ′ .
j 6=j

La probabilité d’acceptation est toujours déterminée par am (t, t′ ) = min(1, A) oùA=jp(k + 1)1
(1 − wj ∗ )n+kδ−k (k + 1)...
wδ−1∗p(k) B(kδ, δ) j


1
dk+1
(1 − wj ∗ )k−1...
(k0 + 1)bk be(wj ∗ |1, k)

pour la création d’une composante et min(1, A−1 ) pour la destruction. Dans cette expression la première ligne correspond au rapport des lois jointes des observations et paramètres. On constate que lavraisemblance des observations n’apparaît pas du tout : en effet elle n’est pas modifiée dans un mouvement de ce type. La seconde ligne contient le jacobien de la transformation (entre crochets) et le rapport
k+1
du choix du type de modification particulier retenu : (k0d+1)b
où k0 est le nombre de composantes vides
k1
avant la création d’une nouvelle .
Ainsi tous les types de modifications sont clairement définis. Pour vérifier la validité de l’algorithme,
les auteurs concluent en remarquant que la chaîne simulée est clairement irréductible et apériodique.

2.4

Les résultats obtenus

La fin de l’article de Richardson et Green est consacrée à la mise en pratique de cet algorithme et à despropositions d’améliorations de la méthode. Rappelons notamment que les propositions de modifications
des paramètres dans les étapes (e) et (f) n’ont fait l’objet d’aucune optimisation selon les auteurs : ainsi lesperformances de l’algorithme devraient donc pouvoir être améliorées. L’article est suivi d’une discussion
dans lesquelles d’autres propositions ou critiques sont faites, certaines ayant trouvé des prolongements
dans d’autres articles par la suite. Ces commentaires et propositions seront l’objet de la partie 3. La finde la partie 2 est consacrée à la présentation des résultats obtenus par Richardson et Green.
1 Remarquons que dans l’article original de Richardson et Green, la puissance du dernier (1 − w ∗ ) n’est pas k − 1 mais
jk, il s’agit d’une erreur de calcul repérée par Tobias Rydén que les auteurs ont corrigé sur leur site web peu de temps après
la parution de l’article.

26

2.4.1Les données

L’algorithme proposé a donc été appliqué à trois ensembles de données différents :- le premier concerne des mesures de l’activité enzymatique dans le sang d’un groupe de 245 individus.
L’enzyme étudiée étant impliquée dans le métabolisme de substances carcinogéniques, la classification est
ici un objectif important car elle peut permettre de détecter différents groupes, génétiquement différents,
dans la population. On a pour cet ensemble de données : R = 2.86 et ξ = 0.45, donc κ = 0.122 et h = 1.22.
- le second concerne la mesure de l’acidité dans 155 lacs du Wisconsin. On a ici : R = 4.18 et ξ = 5.02,
donc κ = 0.057 et h = 0.573.
- le dernier, dit "galaxy data set", contient la mesure de vélocité de 82 galaxies. On a ici : R = 25.11
et ξ = 21.73, donc κ = 0.0016 et h = 0.016. Le "galaxy data set" figure en annexe.
Ces trois jeux de données avaient déjà été étudiés par d’autres auteurs dans le cadre de modèle demélanges, en particulier le "galaxy data set" (au moins trois autres articles). Ceci permet une première
comparaison des résultats obtenus avec des méthodes plus classiques.
2.4.2Interprétation des résultats

Des simulations selon la loi a posteriori ont été menées dans chacun des trois cas ci-dessus pas lesauteurs. Deux remarques sont à faire :- les auteurs ont réalisé 200000 itérations de l’algorithme décrit ci-dessus (ou 200000 sweeps), dont les
100000 premières on été considérées comme burn-in period. Autrement dit, n1 = 100000 et n2 = 200000.
Aucun critère quantitatif n’est proposé pour le choix de ces valeurs, Richardson et Green affirmentd’ailleurs qu’elles sont vraisemblablement supérieures à ce qui est nécessaire.
- dans les trois cas, le nombre de composantes n’a jamais dépassé 24, donc le choix de kmax = 30 n’aété d’aucune conséquence.
Pour ce qui est du paramètre k, l’interprétation des résultats est extrêmement simple puisqu’on obtientles probabilités a posteriori de chaque valeur possible.
Par exemple, pour le premier jeu de données (enzymes), on a :p(k=...|y)
1
2
3
4
5
6
7
8
9
≥ 10

Valeur obtenue
0.000
0.023
0.290
0.317
0.206
0.095
0.041
0.017
0.007
0.003

On peut représenter ces probabilités sur un histogramme, ce qui est fait par exemple dans Green(1995) (cf. figure 3).
Ainsi, le nombre de composantes est très vraisemblablement 3, 4 ou 5 (avec probabilité 0.813), lavaleur de k la plus probable est 4. Toutes les valeurs nécessaires pour le "choix de modèle" peuvent êtrecalculées à partir de ces valeurs, par exemple les facteurs de Bayes :p(k = k1 |y) p(k = k2 )p(k = k2 |y) p(k = k1 )Pour ce qui est de l’estimation de la densité obtenue, les auteurs proposent plusieurs approches :- on peut choisir, à l’aide des valeurs de p(k|y), une valeur de k particulière k0 . Ensuite, on ne considère
plus que les étapes de la simulation où k = k0 . On retient alors l’estimateur :Bk1 ,k2 =b w)
fb(yi |k = k0 , θ,
b =
27

k0
X
j=1

w
bj f (.|θbj )Fig. 3 – Loi a posteriori de k.
où w
b et θb sont des résumés des lois a posteriori des poids et des paramètres, par exemple les moyennes
des valeurs pour lesquelles k = k0 :n2X

Ikm =k0 θm

m=n1
θb =
n2
X

Ikm =k0

m=n1

Le problème est que, comme on le verra par la suite, les densités a posteriori sont souvent multimodales
et les résultats obtenus par cette méthode ne sont pas très satisfaisants.
- une deuxième approche consiste, toujours à k = k0 fixé, à estimer :E[f (.|k, w, θ)|k = k0 , y]en calculant l’estimateur :
n2
X

m=n1
fe(.|k = k0 , θ, w) =Ikm =k0

k0
X

wjm f (.|θjm )

j=1

n2
X

Ikm =k0

m=n1- enfin la dernière proposition de Richardson et Green consiste à sortir de l’optique paramétrique pourse concentrer sur l’estimation de la densité, en estimant :E[f (.|k, w, θ)|y]autrement dit, on n’impose plus une valeur fixe de k. L’estimateur proposé est le suivant :
kX
max

k =0
fe(.|k, θ, w) = 0

fe(.|k = k0 , θ, w)card{l ∈ {n1 , ..., n2 } : k l = k0 }28

n2 − n1

L’inconvénient des deux dernières méthodes est que la quantité approchée, par exemple E[f (.|k, w, θ)|k =k0 , y], n’est pas elle-même la densité d’un mélange.
Une proposition pour remédier à cela est faite par Christian Robert dans la discussion de l’article
de Richardson et Green : il propose de n’utiliser la méthode des RJMCMC que dans le but de choisir k(donc comme outil de sélection de modèles), puis une fois ce choix fait d’utiliser l’algorithme de Gibbs à
k = k0 fixé, comme cela a été décrit dans la première partie, pour étudier les paramètres.
2.4.3

Sensibilité au choix des lois a priori

Richardson et Green étudient ensuite en détail la dépendance des résultats obtenus aux lois a priorichoisies.
Ils étudient tout d’abord la sensibilité de la loi a posteriori de k aux lois a priori. Ils refont pour cela
les mêmes simulations en prenant différentes valeurs pour les hyperparamètres α, β, g, h et comparent esrésultats obtenus. La loi a posteriori de k s’étant avérée très sensible à a valeur du ration α/β, les auteurs
ont choisi comme on l’a vu un modèle hiérarchique avec β lui-même aléatoire. Une fois ce choix fait, la loi
a posteriori
particulières de α, g et h, tant que lap gde k s’est avérée assez insensible aux choix des valeursR(ordre de grandeur des σj ) restait comprise entre R5 et 20. Les choix pour les simulations
quantité hαR"définitives" ont été faits de façon à ce que cette valeur soit 10 . La loi a priori sur les µj , N (ξ, κ−1 )a une légère influence sur la loi a posteriori de k. Les auteurs observent qu’en réduisant légèrement lavariance κ−1 on favorise un plus grand nombre de composantes car on accepte plus facilement deux
composantes avec des moyennes très proches. Mais si κ−1 devient trop petit, on finit par diminuer lenombre de composantes car on refuse la présence de composantes dont la moyenne serait située près desvaleurs extrêmes des données.
C’est ensuite la dépendance de la loi a posteriori des paramètres (σ 2 , µ) aux lois a priori qui est
examinée par les auteurs. Tout d’abord, ils constatent que lorsque κ−1 diminue, les moyennes µ ont
tendance à être moins dispersées, mais ils remarquent aussi que cet effet est beaucoup plus sensible lorsqueβk est grand. Les σj2 sont, comme k, particulièrement sensibles aux valeurs de α
, là aussi particulièrement
lorsque k est grand. En revanche, une fois β rendu aléatoire, les résultats obtenus sur σ varient peu avec
les valeurs de α, g et h.
2.4.4

Applications à la classification bayésienne

Comme on l’a vu en introduction, un des intérêts des modèles de mélange est que les résultats obtenus
peuvent être utilisés pour faire de la classification, c’est-à-dire pour rattacher les individus à un groupeparticulier, ici à une composante particulière.
Ceci est examiné par Richardson et Green à la fin de l’article On Bayesian analysis of mixtures withan unknown number of components (1997).
On peut bien sûr classer sans grande difficulté les observations yi pour i = 1, ..., n : on dispose pourcela, grâce à la simulation MCMC, d’une approximation deP(zi = j|y, k = k0 )une fois que l’on a fixé le nombre de classes k = k0 . Le problème soulevé par les auteurs est que
lorsqu’on est en présence d’une nouvelle observation, disons y ∗ , on ne pourra obtenir une approximationde :P(z ∗ = j|y, y ∗ , k = k0 )qu’en recommençant complètement les simulations sur ce nouvel échantillon ! Ils proposent alors l’approximation (très intuitive) :Z Z∗∗
P(z = j|y, y , k = k0 ) =P(z ∗ = j|y, y ∗ , k = k0 , θ, w)p(θ, w|y, y ∗ , k = k0 )dθdw
... =

Z Z

P(z ∗ = j|y ∗ , k = k0 , θ, w)p(θ, w|y, y ∗ , k = k0 )dθdw
29

... ≈Z Z

P(z ∗ = j|y ∗ , k = k0 , θ, w)p(θ, w|y, k = k0 )dθdw

On peut estimer cette dernière intégrale par la quantité :!
(y ∗ − µj )2
q
exp −
2σj2
2πσj2


k0
X
w
(y ∗ − µJ )2
p J exp −
2σJ2
2πσJ2
J=1
wj

On classe enfin les observations dans la classe à laquelle leur probabilité d’appartenance est la plus
forte :
zbi = arg max P(zi = j|y, k = k0 )jzb∗ = arg max P(zi = j|y, y ∗ , k = k0 )j

3Critiques et extensions

Dans la dernière partie de l’article, les auteurs proposent quatre voies pour prolonger, améliorer ou
approfondir leur travail :- la première concerne les performances du simulateur. La question est de savoir si les types demouvements retenus sont optimaux.
- la deuxième concerne la présentation et l’interprétation des lois a posteriori. Pour les auteurs, denouvelles voies devraient être explorées, en particulier l’étude de lois a posteriori jointes de différents
paramètres.
- la troisième concerne les lois a priori. En particulier, il y a des cas où il pourrait être intéressant dene pas choisir des µj indépendants a priori.
- enfin la quatrième concerne le modèle lui-même. Tout les auteurs n’ont analysé que des mélanges delois normales. Pour eux, il faudrait étudier d’autres distributions. En particulier, dans certains cas une
densité unimodale peut-être très bien approchée par un mélange de deux lois normales avec des moyennes
très proches, ou par une seule loi plus "plate" que la loi normale. De même, on peut vouloir imposer denouvelles contraintes sur les paramètres (par exemple, certains modèles de mélange en génétique imposent
des contraintes sur les poids).
Dans la discussion et des articles parus depuis, certains de ces problèmes sont abordés : les performances du simulateur et la présence de lois non normales en particulier. En revanche, d’autres points
sont critiqués dans la discussion, en particulier le problème d’échange des labels des paramètres au cours
des simulations (déjà abordé pour le cas où k est fixe dans la première partie), l’absence d’un critère
quantitatif pour déterminer quand la convergence de la méthode est effectuée (valeurs de n1 et n2 ), ainsi
que tout simplement le choix de méthodes bayésiennes dans un tel cadre. On va détailler maintenantchacun de ces points.

3.1

Performances du simulateur

3.1.1

Commentaires des auteurs sur les résultats obtenus

Les auteurs examinent plusieurs critères. Tout d’abord, le simulateur a-t’il bien exploré toutes lesvaleurs probables de k ? Pour ceci, il peut être intéressant de regarder quelle est la proportion de mouvements de types (e) (séparation ou regroupement de deux composantes) qui ont été acceptés : elle varie
entre 8% et 14% selon les cas.
La présence de composantes vides pourrait biaiser l’inférence sur k (modifications de type f ), ceci
n’est pas le cas selon les auteurs car au cours des simulations de nombre moyen de composantes vides est
resté assez faible (selon les cas entre 0.1 et 0.57).
30

Le choix des valeurs initiales n’a que peu d’influence : les auteurs ont commencé toutes les simulations
avec k = 1, la première séparation de composantes a toujours été acceptée en moins de 5 étapes. Demême, en commençant avec k = 20 il fallait moins de 100 étapes pour descendre à k = 10.
Les auteurs remarquent cependant que l’inférence sur k, notamment le calcul de facteurs de Bayes,n’a de sens que pour les valeurs de k qui ont été visitées assez souvent par la chaîne de Markov simulée.
Pour ce qui est des valeurs des paramètres, les auteurs remarquent que leur simulateur pourrait mêmeavoir des performances supérieures à une MCMC fonctionnant à k fixé. Ils utilisent pour ceci un exemple
simple avec 200 "observations" :∀i ∈ {1, ..., 50}, xi ∼ N (2.5, 1)∀i ∈ {51, ..., 100}, xi ∼ N (4, 1)∀i ∈ {101, ..., 150}, xi = −xi−100∀i ∈ {151, ..., 200}, xi = −xi−100
et font "comme si" (x1 , ..., x200 était un échantillon de variables indépendantes et identiquement
distribuées provenant d’un mélange de distributions normales. Ils étudient alors le cas (volontairement
mal spécifié) où k = 3, par deux méthodes différentes :- par MCMC avec k = 3 fixe.
- par RJMCMC, en ne retenant que les étapes où k = 3.
Il y a en réalité quatre composantes dans ce "mélange". L’algorithme simple, une fois regroupée deux
composantes, en fait les deux de moyenne −4 et −2.5 ou les deux de moyennes 2.5 et 4, a énormément demal à changer de configuration. En revanche, l’algorithme utilisant les RJMCMC change fréquemmentde configuration puisqu’il peut fréquemment passer à 4 composantes et ensuite revenir à 3 en regroupantdeux composantes qui ne seront pas toujours les mêmes.
3.1.2

Une proposition d’amélioration de la méthode

Brooks, Giudici et Roberts (2000) remarquent que les simulateurs MCMC ont de très bonnes performances et que ce qui peut ralentir la convergence des méthodes utilisant les RJMCMC est la trop grande
proportion de rejet des propositions de mouvements où k varie. Leur article étudie donc de nouvelles
propositions de modifications dont le but est d’augmenter cette proportion.
La probabilité d’acceptation étant de la formeam (t, t′ ) = min(1, A)
am (t′ , t) = min(1, A−1 )
où
A=

∂t′
p(t′ )j(t′ )p(t)j(t)q(u) ∂(t, u)

il peut être judicieux de chercher des propositions de modifications pour lesquelles A est le plus proche
possible de 1. Les auteurs proposent différentes méthodes visant à augmenter la proportion de "sauts"d’une valeur de k à une autre acceptés. Les résultats obtenus améliorent sensiblement les performances
des simulateurs.

313.1.3

Une approche "concurrente" : méthode par diffusion avec sauts

Un algorithme différent a été proposé par David Phillips et Adrian Smith (1996) pour résoudre parune approche bayésienne le problème de la détermination du nombre de composantes dans un mélange.
Il repose que des propositions déjà faites par d’autres auteurs, notamment Grenander et Miller (1994).
Notons que l’idée quand à l’inférence est exactement la même que celle de Richardson et Green, seule lafaçon de simuler dans la loi a posteriori diffère.
On présente rapidement en 3.1.3 leur méthode, puis une généralisation proposée par M. Stephens(2000), et enfin en 3.1.4 la comparaison entre ces méthodes et les RJMCMC menée par Olivier Cappé,
Christian Robert et Tobias Rydén (2003).
Présentation de la méthode de Phillips et Smith Tout d’abord, comme pour les RJMCMC laméthode proposée ne se limite pas aux modèles de mélanges mais peut être utilisée pour l’inférence bayésienne dans n’importe quel cas où la dimension même du paramètre à estimer est inconnue. Des exemples
similaires à ceux utilisés par Green (1995) sont présentés : détection d’un nombre inconnu d’objets sur
une image, détection d’un nombre inconnu de changements de tendance sur une série temporelle, choix
du nombre de regresseurs dans un modèle linéaire, et l’estimation des paramètres d’un modèle de mélange
avec un nombre inconnu de composantes.
Le modèle Le modèle retenu est pratiquement le même que celui utilisé par Richardson et Green (1997),
avec quelques simplifications. La seule différence repose sur la façon dont est simulé un échantillon suivant
la loi a posteriori : au lieu d’utiliser une chaîne de Markov, qui est un processus à temps discret, Phillips
et Smith simulent un processus à temps continu dont la loi stationnaire est la loi a posteriori recherchée :une diffusion avec sauts, les sauts correspondant au changement du nombre de composantes. Une autre
différence notable avec Richardson et Green est que Smith et Phillips n’introduisent pas les variables
indicatrices de la classe des observations, zi .
Le modèle est le suivant :yi |k, θ, w ∼

k
X
j=1

wj f (.|θj ) = f (yi |k, θ, w)où θj = (µj , σj2 ) et f (.|θj ) est la densité d’une loi N (µj , σj2 ),µj ∼ N (µ, τ 2 )avec la contrainte µ1 < ... < µk ,σj2 ∼ IG(γ, δ)(w1 , ..., wk )|k ∼ D(1, ..., 1)

c’est-à-dire une loi uniforme sur le simplexe de Rk . Enfin, Phillips et Smith désiraient faire le choix
habituel k ∼ P(λ). En fait, comme k = 0 est interdit, ils retiennent :′∀k ′ ≥ 1, P(k = k ′ ) =λk
λ(e − 1)k ′ !

Contrairement à Richardson et Green, Phillips et Smith n’introduisent pas un niveau hiérarchique
supplémentaire, ici µ, τ , γ, δ et λ sont constantes.
Enfin on note :θ = (θ1 , ..., θk ) ∈ Θk32

Le simulateur Les auteurs cherchent donc à simuler selon la loi a posteriori p(θ, k|y) notée poursimplifier π(θ, k). Ils utilisent pour ceci une diffusion pouvant sauter d’un espace Θk dans un autre. On
la notera {(θ, k)(t) , t ∈ R+ }.
On étudie tout d’abord la dynamique des sauts. Aux questions déjà soulevées par Green s’ajoutent
de nouvelles dues au temps continu : quand proposer des sauts ?
Si à la date t on est en (θ, k), la probabilité de quitter Θk est donnée par l’intensité de sautq((θ, k), (φ, h)) où h 6= k. En fait, la probabilité (infinitésimale) de sauter entre les dates t et t + dt
de (θ, k) vers un nouvel état inclus dans A ⊂ Θh est donnée par :Z
dtq((θ, k), (φ, h))dφA

Les auteurs notent J (θ, k) l’ensemble des (φ, h) qui peuvent être atteint en un seul saut, appelé"ensemble de saut". Ces ensembles doivent être "réversibles", c’est-à-dire que :(φ, h) ∈ J (θ, k) ⇒ (θ, k) ∈ J (φ, h)Pour décider quand aura lieu un saut, il faut calculer l’intensité marginale de saut :XZq((θ, k), (φ, h))dφr(θ, k) =h(φ,h)∈J (θ,k)

Ainsi, la probabilité (infinitésimale) de quitter Θk entre les dates t et t + dt de (θ, k) est donnée par :r(θ, k)dt
Les temps de sauts sont simulés en générant aléatoirement e1 , e2 , ... indépendantes selon une loi expontentielle de paramètre 1, E(1), et en choisissant ti date du i-ème saut telle que :Z tir((θ, k)(t) )dt ≥ ei
ti−1

avec t0 = 0 par convention (≥ et non pas = car on ne pourra bien sûr pas simuler directement unprocessus à temps continu mais une version discrétisée de ce processus, le pas de temps étant ∆).
Si juste avant la date ti on est en (θ, k), il faut choisir le nouvel état suivant un noyau de transition
Q à déterminer.
Deux propositions sont faites pour Q par les auteurs, l’une s’inspirant de l’algorithme de Gibbs etl’autre reposant sur celui de Hastings-Metropolis. Seule cette deuxième solution est appliquée aux modèles
de mélanges :q((θ, k), (φ, h)) = Cph ePh (φ)
koù pk = (eλλ−1)k! est la probabilité a priori de k, et Ph (φ) est le logarithme de la densité a priori deφ = (θ, σ 2 ) sachant que l’on a h composantes :
!#
"

h
Y
(µ −µ)2
− δ2
1
δγ
1
σ
− j2τ 2j√e
Ph (φ) =2 )γ+1 e
2
Γ(γ)
(σ
2πτ
j
j=1
A chaque date de saut, on simule donc (φ, h) suivant la loi a priori restreinte à J (θ, k). On accepte
ensuite (φ, h) avec probabilité :min(1, eLh (φ)−Lk (θ) )où L est la vraisemblance, c’est-à-dire que :Lh (φ) =hn XY

i=1 j=1

33

wj f (yi |θj )Sinon, on rejette le saut c’est-à-dire qu’on reste à l’état (θ, k). On obtient :Q((θ, k), (φ, h)) =C min(1, eLh (φ)−Lk (θ) ph ePh (φ) )r(θ, k)Décrivons maintenant comment est réalisé un tel saut. Un seul type de saut est envisagé ici : la créationou la suppression d’une composante, c’est-à-dire que ce saut est proche de l’étape (f) de l’algorithme deRichardson et Green. Il n’y a pas ici de possibilité de regrouper ou séparer deux composantes (étape (e)chez Richardson et Green).
Pour la création d’une composante, tout d’abord on génère µ∗ ∼ N (µ, τ 2 ) puis σ ∗2 ∼ IG(γ, δ). Resteà définir la probabilité w∗ associé. Si µj < µ∗ < µj+1 on simule un nouveau triplé (wj , w∗ , wj+1 ) suivant
sa loi conditionnelle aux autres wj ′ fixés. Si µ∗ < µ1 on simule un nouveau couple (w∗ , w1 ) suivant sa loi
conditionnelle aux autres wj ′ , de même si µ∗ > µk . On réindexe ensuite les différentes composantes. Le
calcul explicite de r(θ, k) donne finalement la probabilité d’acceptation de ce saut.
Pour la destruction d’une composante, on choisit la composante j à effacer (de façon uniforme), on
simule un nouveau couple (wj−1 , wj+1 ) suivant sa loi conditionnelle aux autres wj ′ fixés, on réindexe
ensuite les différentes composantes et on calcule la probabilité d’acceptation.
Reste donc à déterminer la dynamique entre les sauts. Elle est donnée par l’équation différentielle
stochastique : ∀t ∈]ti−1 , ti [,


dt d
d
(t)
(t)
(t)
dθ =
(Lk (θ )) + (Pk (θ )) + dWtk
2 dθdθoù {Wtk , t ≥ 0} est un mouvement brownien de dimension dim(Θk ). Bien sûr, comme on ne simule
pas directement le processus à temps continu mais une version discrétisée (avec pas de temps ∆), on
utilisera l’équation discrétisée :


√ (t)
d
∆ d
(t)
(t)
(t+∆)
(t)
(Lk (θ )) + (Pk (θ )) + ∆zk
θ
=θ +
2 dθ
dθ
(t)

où les zk sont indépentants et identiquement distribués selon N (0, 1).
Les résultats Les auteurs utilisent appliquent ensuite leur méthode sur deux jeux de données, dont legalaxy data set utilisé par Richardson et Green.
On peut comparer ici les résultats des deux méthodes, tout en remarquant que ces comparaisons n’ont
que peu de sens puisque les lois simulées ici ne sont pas les mêmes : les deux articles utilisaient en effet
des lois a priori différentes donc les lois a posteriori ne peuvent être les mêmes. En particulier, Phillips
et Smith n’ont pas introduit un niveau hiérarchique supplémentaire comme Richardson et Green. Leurs
résultats dépendent donc plus fortement des lois a priori (Phillips et Smith ont retenu ici : µ = 20,
τ 2 = 100, δ = 0.5 et λ = 3).
p(k=...|y)
3
4
5
6
7
8
9
10

Richardson
Green
0.061
0.128
0.182
0.199
0.160
0.109
0.071
0.040

&

Phillips & Smith
0.025
0.394
0.321
0.217
0.043

Les auteurs reconnaissent en fait dès le début de l’article qu’aucune comparaison théorique de leurméthode avec les RJMCMC n’avait pas été menée lors de la parution de leur article.

34

Contribution de Stephens Comme on l’a vu, Phillips et Smith traitent en fait de nombreux exemples
dans leur article et ne passent que peu de temps sur celui des mélanges de distributions. Stephens (2000a)
propose une description complète de l’adaptation de la méthode aux modèles de mélanges.
Sa méthode est légèrement différente puisqu’il choisit une autre dynamique de saut : il en résulte que
les sauts ne sont jamais refusés. En revanche, si le processus saute dans un état très improbable il n’yrestera que très peu de temps avant de sauter de nouveau.
Pour Stephens, ceci est un point positif, car en visitant même les zones de faible probabilité de la loi
a posteriori, son simulateur aurait de meilleures performances que celui de Green et Richardson.
3.1.4

Comparaison simulateurs à temps continu / RJMCMC

Cappé, Robert et Rydén (2003) ont comparé les différents types de simulateurs. Ils remarquent tout
d’abord que dans l’approche de Green et Richardson par RJMCMC, deux types de sauts sont proposés :- regrouper ou séparer des composantes (type (e)),- créer ou détruire des composantes (type (f)),
alors que les simulateurs à temps continu, qu’il surnomment BDMCMC (birth and death MCMC), de
Phillips et Smith et de Stephens n’utilisent que des créations ou destructions de composantes. Remarquantque ce choix est arbitraire ils créent des simulateurs plus généraux à temps continu, qu’ils appellent
CTMCMC (continuous time MCMC) utilisant aussi des regroupements et séparations de composantes.
La partie théorique de leurs comparaisons montre qu’en fait, les CTMCMC ne sont pas "beaucoup
plus générales" que les RJMCMC puisque, sous certaines conditions, une CTMCMC peut être approchéepar une suite de RJMCMC dont le pas de temps tendrait vers 0.
La partie empirique compare quatre simulateurs (RJMCMC avec sauts de type (f) seulement, de type(e) et (f), BDMCMC et RJMCMC) et montre que si les résultats sont comparables, le classement obtenu
est le suivant :- le meilleur simulateur serait la RJMCMC avec sauts de type (f) seulement,- presque au même niveau, BDMCMC,- moins bons : RJMCMC avec sauts de type (e) et (f) et CTMCMC, cette dernière méthode étant enplus très coûteuse en temps de calcul.

3.2

Problème des échanges des labels des paramètres

Pour beaucoup d’auteurs ayant contribué à la discussion, il s’agit la du principal problème de cetteméthode. Matthew Stephens utilise un exemple simple pour montrer que la contrainte µ1 < ... < µk n’est
pas suffisante : il considère le mélange
0.33N (0, 1) + 0.33N (5, 1) + 0.34N (5, 4)
et utilise un algorithme de Gibbs simple pour obtenir un échantillon de la loi a posteriori. Ses résultats
montrent clairement ce que l’on pouvait attendre : la contrainte µ1 < µ2 < µ3 ne sépare pas correctement
les deux composantes ayant la même moyenne, et les labels de ces deux composantes sont échangés aucours des simulations. Ceci est problématique car la méthode retenue ne fonctionne donc pas correctement sur l’exemple du "chapeau méxicain" (figure 2, dans la première partie). Dans le cas du "chapeau
méxicain" on peut résoudre ce problème en prenant σ12 < σ22 < σ32 comme contrainte à la place de cellesur les composantes. Mais ceci ne résout pas le problème dans le cas de l’exemple proposé par Stephens,où deux composantes ont la même variance.
Gilles Celeux remarque que la contrainte µ1 < ... < µk modifie la forme de la loi a posteriori desmoyennes et peut biaiser leur estimation, en particulier dans le cas où deux composantes auraient desmoyennes très proches. Pour lui, il vaudrait mieux ne pas imposer de contraintes pendant les simulations puis de permuter ensuite les résultats des MCMC à l’aide d’un algorithme à définir. Cependant,
il remarque aussi que cette méthode pourrait être difficile à implémenter dans le cas où k est inconnu.
Il affirme en attendant que le problème de l’échange du label des composantes est une des difficultés del’analyse bayésienne des modèles de mélanges qui n’avait lors de la parution de l’article de Richardson etGreen pas trouvé de solution satisfaisante.
La proposition d’un algorithme permettant de permuter après les simulations a été étudiée depuis par
exemple par Celeux, Hurn et Robert (1997) ou Stephens (2000b) qui obtiennent des résultats satisfaisants.
35

3.3

Généralisations

3.3.1

Composantes non normales

Les auteurs remarquent que k peut avoir été surestimé du fait de la présence de composantes nonnormales. En effet, comme on l’a mentionné précédemment, dans certains cas une densité unimodale
peut-être bien approchée par un mélange de deux lois normales avec des moyennes très proches, ou parune seule loi, plus "plate" que la loi normale. Il faudrait donc être en mesure de généraliser cette méthodeà des lois autres que la loi normale. La même remarque est faite par Yung-Hsin Chien et Edward Georgedans la discussion : l’inférence sur k pourrait être faussée par l’hypothèse de normalité des composantes.
Remarquons que ceci n’est pas un problème si le seul but est l’estimation de densité, c’est-à-dire si
on considère les RJMCMC comme une alternative à la méthode du noyau par exemple. Mais dans le casoù le but premier est la classification et où les composantes ont une signification physique, ceci pose unproblème important.
3.3.2

Cas multivarié

Dès le début de l’article, Richardson et Green affirmaient que si leur exposé se limitait aux mélanges
de gaussiennes unidimensionnelles, leur méthode est plus générale et devrait pouvoir s’adapter au cas
multidimensionnel. Dans la discussion, D. M. Titterington et Mike West remarquent qu’une généralisation
au cas multivarié n’est pas évidente, et pourrait impliquer de nouveaux types de mouvements. Gilles
Celeux note, en particulier, que les contraintes d’identifiabilité des composantes, µ1 < ... < µk , ne sontpas directement généralisables au cas multivarié. De nombreux auteurs proposent d’introduire une mesurede la proximité de deux composantes afin de favoriser le regroupement de deux composantes "proches".
Cette généralisation a été étudiée par Guillaume Saint-Pierre, cf. Saint-Pierre et Garel (2001). Outreles problèmes de généralisation des contraintes, le principal problème rencontré est la difficulté de la
généralisation des formules concernant le saut de type (e) c’est-à-dire la séparation ou la réunion de deux
composantes :w j ∗ = w j1 + w j2
w j ∗ µ j ∗ = w j1 µ j1 + w j2 µ j2
wj ∗ (µ2j ∗ + σj2∗ ) = wj1 (µ2j1 + σj21 ) + wj2 (µ2j2 + σj22 )Si les deux premières de ces équations se généralisent directement, ce n’est pas le cas de la troisième,
puisqu’il faut maintenant passer à des matrices de variance-covariance. De plus, la complexité des calculs
des probabilités d’acceptation devient problématique dans ce cas.
3.3.3

A priori non informatifs

Même si l’utilisation d’a priori totalement non informatifs n’est pas possible dans ce contexte, plusieurs
auteurs remarquent qu’il serait intéressant d’utiliser des a priori partiellement non informatifs afin delimiter la dépendance des résultats au choix des lois a priori.
Par exemple Christian Robert remarque qu’il peut être avantageux de paramétrer différemment lemodèle de mélanges : la composante i est vue comme une perturbation locale des précédentes. Dans lecas normal on obtient :p0 N (θ0 , τ02 ) +

k−2
X
i=1

(1 − p0 )...(1 − pi−1 )pi N (θ0 + ... + τ0 ...τi−1 θi , τ02 ...τi2 )...

2... + (1 − p0 )...(1 − pk−2 )N (θ0 + ... + τ0 ...τk−2 θk−1 , τ02 ...τk−1)

Cette paramétrisation permet l’utilisation de lois a priori non informatives comme :π(θ0 , τ0 ) =
36

1
τ0

∀i, pi ∼ U [0, 1]∀i > 0, τi ∼ U [0, 1]∀i > 0, θi ∼ N (0, ζ 2 )D’autres propositions ont été faites. Roeder et Wasserman (1997) ont proposé des lois a priori partiellement non informatives pour les paramètres dans les modèles de mélange à k fixé qui devraient pouvoirs’adapter au cas où k varie. L’idée est de ne plus supposer l’indépendance des paramètres. Par exemple
les moyennes sont supposées être proches les unes des autres, groupées autour d’un "centre" A qui suit
lui une loi complètement non informative (et donc impropre).

3.4

Critère de convergence

De nombreux auteurs remarquent l’absence de critère quantitatif pour déterminer quand les simulations doivent s’arrêter, c’est-à-dire pour choisir n2 . Ceci est désigné comme la recherche d’un "diagnostic
de convergence". Selon S. Brooks, les 200000 étapes de simulations ne sont pas suffisantes, ce que reconnaissent Richardson et Green tout en remarquant que de nouvelles simulations, avec 500000 étapes,
n’apportaient pas de modification sensible des résultats.

3.5

Discussion du choix des méthodes bayésiennes

3.5.1

Complexité des MCMC dans ce contexte

Certains auteurs critiquent le choix de méthodes bayésiennes dans ce contexte. Par exemple, Murray
Aitkin remarque que la méthode des RJMCMC est beaucoup plus complexe, et coûteuse en temps decalcul, qu’une méthode classique. Celle-ci pourrait consister en l’utilisation du maximum de vraisemblance
pour différentes valeurs de k, puis d’une procédure de choix de modèle reposant sur la vraisemblance pourdéterminer k.
Le problème est en fait le même que pour les méthodes bayésiennes : les équations définissant l’estimateur du maximum de vraisemblance (EMV) ne peuvent être résolues analytiquement. L’EMV devra
donc être approché par l’algorithme EM, à k fixé. Le choix de k est ensuite classiquement réalisé parl’utilisation de critères comme l’AIC ou le BIC.
L’algorithme EM converge en général assez rapidement et cette méthode est finalement moins complexe que les RJMCMC. Cependant, deux problèmes (au moins) sont posés par cette méthode :- l’algorithme EM peut très bien converger vers un maximum local de la vraisemblance qui ne sera
pas le maximum global.
- il n’y a pas, dans le cas des modèles de mélanges, de justification théorique à l’utilisation de l’AIC
et du BIC, même si ils donnent en pratique des résultats convenables (quoique surestimant en général lenombre de composantes).
Cette critique sur la complexité est nuancée par K. Mengersen et A. George qui remarquent que laméthode n’est finalement pas très complexe à implémenter.
De même, d’autres auteurs estiment que la méthode est d’un grand intérêt pratique : Walter Gilks
présente dans la discussion un exemple de modèle en génétique qui pourrait nécessiter l’utilisation de
RJMCMC.
Cependant, cette critique reste valable dans le cas mutidimensionnel, comme ceci a été mentionné
précédemment, cf. Saint-Pierre et Garel (2001).
3.5.2

Sensibilité du choix de modèle aux lois a priori

Murray Aitkin remarque aussi qu’il faudrait étudier la dépendance de la loi a posteriori de k à la loi apriori choisie pour k. Dans le cas multidimensionnel, selon Saint-Pierre et Garel (2001), cette dépendance
est importante, mais pourrait être réduite par l’introduction d’un niveau hiérarchique supplémentaire,par exemple :37

k ∼ P(λ)

où λ ne serait plus une valeur fixe mais suivrait lui-même une loi donnée.
La sensibilité aux autres lois a priori a été étudiée en détail par Richardson et Green dans leur article.
On voit ici cependant l’intérêt d’être en mesure d’utiliser des a priori non informatifs (au moins
partiellement) et pour cela de prendre en compte les propositions (présentées précédemment) de Robertdans la discussion ou de Roeder et Wasserman (1997).
3.5.3

Une comparaison numérique

Geoffrey McLachlan et David Peel proposent une comparaison des résultats obtenus par Richardson
et Green avec ceux obtenus, sur les mêmes données, en utilisant les mêmes méthodes que celles proposéespar Aitkin (algorithme EM puis critère AIC). Les résultats obtenus, au moins sur k, sont assez similaires
à ceux de Richardson et Green, sauf pour le jeu de données sur l’acidité, pour lequel les méthodes
bayésiennes donnent k ∈ {2, ..., 6} alors que les méthodes classiques rejettent k > 3.

4

Mise en pratique

4.1

Le programme

4.1.1

Choix effectués

Le modèles est exactement le même que celui étudié par Richardson et Green. La seule différence est
que je n’ai pas imposé de contrainte d’ordre lors des simulations (ce qui impose une légère modification
dans la formule de la probabilité d’acceptation).
J’ai réalisé le programme sous le logiciel R. Le code source figure en annexe.
4.1.2

Problèmes rencontrés

La méthode n’est pas très difficile à programmer en général mais pose certaines problèmes : complexitédes formules des probabilités d’acceptation, stockage des résultats (dimension variable !)...
En revanche, la lenteur d’execution du programme obtenu a été un réel problème.

4.2

Résultats sur données simulées

4.2.1

Présentation des données

Je présente ici les résultats obtenus pour un fichier de 30 "observations" simulées suivant une loi :21
N (0, 1) + N (4, 1)
3
3
Ces simulations figurent en annexe.
Dans un premier temps, j’ai appliqué la méthode décrite dans la partie 1, c’est-à-dire que l’on garde
k fixé, égal à 2. Sinon, j’ai utilisé exactement la même structure de lois que Richardson et Green. Dans
un second temps, j’ai utilisé la méthode avec saut réversible (k "inconnu").
4.2.2

Résultats obtenus à k fixé

J’ai utilisé seulement 5000 itérations après la période de démarage (très courte, de l’ordre de quelquesétapes).
En prenant comme résumé des simulations les moyennes a posteriori, on obtient l’estimation suivantede la loi :0.6445N (0.1369, 12655) + 0.3555N (4.3909, 1.1495)
Des lois a posteriori sont représentées dans les figures 4 (moyenne de la première composante), 5(variance) et 6 (poids).
38

Fig. 4 – Loi a posteriori de µ1

Fig. 5 – Loi a posteriori de σ12
4.2.3

Résultats des RJMCMC

100000 simulations ont été réalisées dans ce cas bien que la convergence semble avoir été assez rapide.
Les résultats obtenus pour la loi a posteriori de k sont représentées dans la figure 7.
Ainsi, la méthode surestime le nombre de composantes dans ce cas. Examinons les résultats conditionnellement à k = 2 et k = 3. Pour exploiter les résultats, on distingue les composantes par la mêmeméthode que Richardson et Green (contrainte d’ordre sur les moyennes).
Pour k = 2, voici les résultats obtenus en estimant les paramètres par leurs espérances a posteriori(paramètre estimé / vraie valeur) :ComposanteMoyenneVariancePoids

1
0.146 / 01.264 / 1
0.648 / 0.667

2
4.409 / 41.110 / 1
0.352 / 0.333

Examinons maintenant les résultats obtenus pour k = 3 composantes. Ceci pourra peut-être nous
aider à comprendre pourquoi k = 3 est plus probable au vu des données que k = 2.
ComposanteMoyenneVariancePoids

1-0.585
1.330
0.362

2
1.331
1.415
0.341

3
4.687
1.037
0.297

Il faut faire très attention en interprétant ces résultats. On pourrait être tenté de croire, par exemple,
que la première composante dans k = 2 a été "coupée en deux". En fait il n’en est rien, la moyenne
39Fig. 6 – Loi a posteriori de w1

Fig. 7 – Loi a posteriori de k pour les données simulées
a posteriori n’a pas vraiment de signification pour la deuxième composante dans k = 3. Pour s’enconvaincre, il suffit de regarder les lois a posteriori des trois moyennes (sous la contrainte d’ordre) : cf.
figure 8.
La deuxième composante, dans le cas k = 2, est tantôt un "double" de la première, tantôt un double
de la troisième.
Afin d’empêcher ce phénomène de se produire trop souvent, on peut utilise une loi a priori sur kqui "pénalise" plus les nombres élevés de composantes, en revenant par exemple à la loi de Poisson deparamètre λ, initialement proposée par Richardson et Green. Avec λ = 1, on obtient déjà des résultats
beaucoup plus satisfaisants avec une probabilité a posteriori pour k = 2 de l’ordre de 0.7, et pour k = 3
de l’ordre de 0.25.
Ces résultats donnent ceci dit un facteur de Bayes proche de 1, très légèrement en faveur de 3 composantes.

4.3

Résultats sur données réelles

4.3.1

Intérêt de travailler sur des données réelles

En général, lorsqu’on étudie un modèle statistique on se donne une famille de lois possibles pour les
observations :C = {Pθ , θ ∈ Θ}

et on fait l’hypothèse que la vraie loi des observations P∗ ∈ C. Ici, on a fait l’hypothèse que les donnéesétaient issues d’un mélange de gaussiennes. Ceci était justifié dans le cas des données simulées suivant
cette loi.
40

Fig. 8 – Les lois a posteriori pour les trois moyennes dans le cas où k = 3 pour les données simulées.
Or, pour des données réelles ceci n’est qu’une approximation. En pratique, dans un modèle de mélanges, ceci se traduit par exemple par le fait que deux composantes normales peuvent être utilisées pourreprésenter une seule classe dont la loi est plus "plate" que la loi normale, ou par la présence d’une
observation atypique qui peut favoriser la création d’une nouvelle classe.
Notons en effet que :- dans les comparaisons entre méthodes classiques et RJMCMC, McLachlan et Peel par exemple
affirmaient que le RJMCMC a une tendance à surestimer le nombre de composantes.
- pour les données simulées avec k = 2, il y a déjà eu surestimation du nombre de composantes enutilisant l’a priori uniforme.
Il est donc intéressant de voir comment réagit la méthode étudiée face à un jeu de données réelles.
4.3.2

Présentation des données

On utilise donc ici le galaxy data set déjà étudié par Richardson et Green. Ces données dont disponibles
en annexe.
4.3.3

Résultats et commentaires

Ici, seulement 50000 simulations ont été réalisées.
La loi a posteriori de k donnée par les simulations est représentée sur la figure 8.

Fig. 9 – Loi a posteriori de k pour le galaxy data set
Le tableau suivant compare les résultats que j’ai obtenu à ceux obtenus par Richardson et Green avec
un a priori uniforme.

41p(k=...|y)
3
4
5
6
7
8
9
10

Richardson
Green
0.061
0.128
0.182
0.199
0.160
0.109
0.071
0.040

&

Mes simulations
0.058
0.272
0.328
0.225
0.092
0.020
0.004
0.000

Ainsi, par rapport aux estimations proposées par McLachlan et Peel dans la discussion, on obtientplus de composantes (pour eux le nombre de composantes était compris entre 2 et 6 maximum).
En revanche, en utilisant une loi de Poisson a priori (λ = 1) on obtient des résultats plus proches deceux obtenus par des méthodes classiques :- le k le plus probable a posteriori est 4 avec une probabilité de 0.56 ;- il y a une forte évidence en faveur de 3 à 5 composantes, 3 a une probabilité 0.28 et 5 de 0.14, ce
qui donne une probabilité totale de 0.98 pour cet intervalle.
- en particulier, la probabilité de k = 2 estimée est nulle, après la burn-in period la MCMC ne revientplus jamais à deux composantes.

42

Bibliographie
Brooks, S.P., Giudici, P. & Roberts, G.O., Efficient construction of reversible jump MCMC proposal
distributions, 2003
Cappé, Olivier, Robert, Christian P. & Rydén, Tobias, Reversible jump, birth-and-death, and moregeneral continuous time MCMC samplers, 2002
Celeux, Gilles, Hurn, Merrilee & Robert, Christian, Computational and Inferential difficulties with
mixture posterior distributions, Journal of the American Statistical Association, 1997
Diebolt, Jean & Robert, Christian P., Estimation of finite mixture distribution through bayesian sampling, Journal of the Royal Statistical Society, 1994
Green, Peter J., Reversible jump Markov chain Monte Carlo computation and Bayesian model determination, Biometrika, 1995
Green, Peter J. & Richardson, Sylvia, On Bayesian analysis of mixtures with an unknown number ofcomponents, with discussion, Journal of the Royal Statistical Society, 1997
Green, Peter J. & Richardson, Sylvia, Correction on "On Bayesian analysis of mixtures with an
unknown number of components", disponible sur le site web de Peter Green, 1997
Grenander, U. & Miller, M., Representation of knowledge in complex systems, with discussion, Journalof the Royal Statistical Society, 1994
McLachlan, Geoffrey & Peel, David, Finite mixture models, Wiley series in probability and mathematical statistics, 2000
Phillips, D. B. & Smith, A. F., Bayesian model comparison via jump diffusions, in Markov Chain
Monte Carlo in practice, Chapman and Hall, 1996
Robert, Christian P., Mixtures of distributions : inference and estimation, in Markov Chain Monte
Carlo in practice, Chapman and Hall, 1996
Robert, Christian P. & Casella, George, Monte Carlo Statistical Methods, Springer texts in statistics,
1999
Roeder, Kathryn & Wasserman, Larry, Practical bayesian density estimation using mixture of normals,
Journal of the American Statistical Association, 1997
Saint-Pierre, Guillaume & Garel, Bernard, Reversible Jump MCMC for multivariate gaussian mixture
models, Technical report, ENSEEIHT, 2001
Stephens, M., Bayesian analysis of mixture models with an unknown number of components - analternative to reversible jump methods, Annals of Statistics, 2000 (a)Stephens, M., Dealing with label-switching in mixture models, Journal of the Royal Statistical Society,2000 (b)Waagepetersen, Rasmus & Sorensen, Daniel, A tutorial on Reversible Jump MCMC with a view towardapplications in QTL-mapping , International Statistical Review, 2001

43