Nous nous intéressons au modèle de mélange suivant, comme en cours : \[ \mathcal{M}=\left\{\pi \phi\left(\cdot ; \mu_1, 1\right)+(1-\pi) \phi\left(\cdot ; \mu_2, 1\right): \mu_1 \in \mathbb{R}, \mu_2 \in \mathbb{R}\right\}, \] avec \(\phi\) la densité gaussienne sur \(\mathbb{R}\) et \(\pi \neq \frac{1}{2}\) connu. La loi a priori de \(\left(\mu_1, \mu_2\right)\) est donnée par: \(\mu_1 \sim \mathcal{N}\left(\delta_1, \frac{1}{\lambda}\right), \mu_2 \sim \mathcal{N}\left(\delta_2, \frac{1}{\lambda}\right), \delta_1, \delta_2 \in \mathbb{R}\) et \(\lambda>0\), avec \(\mu_1\) et \(\mu_2\) indépendantes.

Modèle

Pour cette partie, on se base sur les notes de cours disponibles ici: https://perso.lpsm.paris/~rebafka/#enseignement.

Notons les paramètres de notre mélange par \[\theta=(\underline \pi,\underline \mu,\underline\sigma^2)=\left(\pi_1,\pi_2,\mu_1,\mu_2,\sigma_1^2,\sigma_2^2\right)\] Dans l’exemple du TP on a même: \[\theta=(\underline \pi,\underline \mu,\underline\sigma^2)=\left(\pi_1,\pi_2,\mu_1,\mu_2,1,1\right)\]

où \(0<\pi_k<1\), \(\pi_1 + \pi_2=1\), \(\mu_k\in\mathbb R\) et \(\sigma_k^2>0\), \(k=1,\dots,K\).

Avec les notations bayesiennes usuelles, on a: \[\begin{align*} \theta&\sim\pi(\theta)\\ z_i|\theta&\sim \sum_{k=1}^K\pi_k\delta_{\{k\}},\quad i=1,2\\ x_i|z_i,\theta&\sim \mathcal N(\mu_{z_i},\sigma^2_{z_i}),\quad i=1,2. \end{align*}\]

Loi a priori

On utilise la loi a priori de \(\theta\) suivante : \[\pi(\theta)=\pi(\underline\pi)\pi(\underline\mu|\underline\sigma^2)\pi(\underline\sigma^2) =\pi(\underline\pi)\prod_{k=1}^K\pi(\mu_k|\sigma^2_k)\pi(\sigma_k^2) \] avec \[\underline\pi \sim \mathcal D(\gamma_1,\dots,\gamma_K),\qquad \mu_k|\sigma^2_k\sim\mathcal N\left(\alpha_k,\frac{1}{\lambda_k}\right)\] où \(\mathcal{D}(\gamma_1,\gamma_2)\) désigne la loi de Dirichlet de paramètres \(\gamma_1,\gamma_2\). L’ensemble des hyperparamètres est donc finalement \[ \psi =(\gamma_1,\gamma_2, \alpha_1,\alpha_2,\lambda_1,\lambda_2) \] où \(\gamma_k>0\), \(\alpha_k\in\mathbb R\) et \(\lambda_k>0\).

Loi a posteriori

La loi a posteriori de \(\theta\) sachant \(x=(x_1,…,x_n)\) est donnée par \[ \pi(\theta|x) \infty \prod_{i=1}^n (\sum_{k=1}^K\pi_{k}f_{\mathcal{N}(\mu_{k},\sigma^{2}_{k})}(x_i))\pi(\theta) \] On en déduit que la fonction de log-aposteriori est donnée par \[\log\pi(\theta|x)=\sum_{i=1}^n\log(\sum_{k=1}^K\pi_{k}f_{\mathcal{N}(\mu_{k},\sigma^{2}_{k})}(x_i))+\log\pi(\pi)+\sum_{k=1}^K\log(\pi(\mu_k|\sigma^{2}_{k}))+\sum_{k=1}^K\log(\pi(\sigma^{2}_{k}))+cst\]

où cst est une constante que l’on ne précisera pas

Bonus: Simulation des variances

Il n’est pas demandé dans le TP de simuler les variances \(\sigma_1, \sigma_2\) car elles sont fixes et identiques. Nous avons cependant implémenté cette option et voici les rajouts à faire pour notre modèle: \[\mu_k|\sigma^2_k\sim\mathcal N\left(\alpha_k,\frac{\sigma_k^2}{\lambda_k}\right),\qquad \sigma_k^2\sim \mathcal{IG}\left(\frac{\lambda_k}2,\frac{\beta_k}2\right)\]

où \(\mathcal{IG}(\alpha,\beta)\) désigne la loi Gamma inverse. Les hyperparamètres de notre modèle sont donc: \[ \psi =(\gamma_1,\gamma_2, \alpha_1,\alpha_2,\lambda_1,\lambda_2,\beta_1,\beta_2) \]

Simulation

Voici le pseudo-code et l’explication de notre algorithme.

L’algorithme prend en entrée les arguments suivants:

obs: un vecteur d’observations.
theta.init: un vecteur de paramètres de mélange initiaux.
psi:le vecteur des hyperparamètres pour simuler nos lois.
R: le nombre d’itérations à effectuer (par défaut 1000).
B: le nombre d’itérations de mise en chauffe à ignorer (par défaut 500).

[1] "Paramètres pour la simulation de notre mélange:"
$pi
[1] 0.4326397 0.5673603

$mu
[1] -1.5845629  0.9948023

$sig
[1] 2.9613410 0.9493493

On voit que notre algorithme a bien trouvé nos paramètres \(\mu_1, \mu_2, \sigma_1, \sigma_2\) car les valeurs chaine de markov generé par notre algorithme gravitent autour des vrais valeurs de notre méalnge. On remarque aussi que les valeurs simulées ont moyent de variance pour les paramètres \(\mu_1, \mu_2\) que pour \(\sigma_1, \sigma_2\).

Le but du Gibbs sampler est de créer un échantillon de la loi aposteriori \(\pi(\theta|x)\). Ensuite, cet échantillon de valeurs simulées \(\theta(t)\), \(t=1,…,R\) est utilisé pour évaluer des estimateurs bayésiens comme la moyenne aposteriori. Ainsi, la valeur de \(\mu_1\) est estimé par \[\hat{\mu_1}=\frac{1}{R}\sum_{t=1}^{R}\mu_1^{(t)}\] Ainsi, la valeur de \(\mu_2\) est estimé par \[\hat{\mu_2}=\frac{1}{R}\sum_{t=1}^{R}\mu_2^{(t)}\]

Notre algorithme a bien évalué le log aposteri, en utilisant l’échantillon de valeurs simulées \(\theta^{(t)}\) et des moyennes aposteriori \(\hat{\mu_1}\) et \(\hat{\mu_2}\). En effet, on voit que les valeurs sont proches de -1.6 pour \(\mu_1\) (valeur théorique: -1.58) et de 1 pour \(\mu_2\) (valeur théorique: 0.99).

Gibbs comme outil de simulation

Comme nous l’avons vu précédemment, l’algorithme de Gibbs peut être utilisé comme un outil statistique. On peut aussi noter son utilité pour des densités quand ses lois conditionnelles sont faciles à simuler. Un exemple simple d’utilisation de l’algorithme de Gibbs est la simulation d’une gaussienne 2D de matrice de covariance non diagonale. On s’interesse donc à une gaussienne 2D de coefficient de corrélation \(\rho\). L’algorithme de Gibbs est très simple dans ce cas, il suffit de simuler les lois conditionnelles en sachant que: \[p(y|x) \simeq \mathcal{N}(\rho x,1-\rho^2)\qquad p(x|y) \simeq \mathcal{N}(\rho y,1-\rho^2)\].

Voici une comparaison entre une simulation par Gibbs et par le package “mvtnorm” (qui utilise la méthode de Schuermann ).

Nous avons ici pris comme paramètres: \(\mu_1 = \mu_2 = 0\) et \(\rho = 0.5\). On voit que les données simulées par nos deux méthodes concordent.

Algorithme de Metropolis-Hastings

L’algorithme de Metropolis-Hastings est une classe de méthode de Monte Carlo à chaîne de Markov dont l’idée principale est de générer une chaîne de Markov de telle sorte que sa distribution stationnaire soit la distribution cible. (à long terme les deux échantillons vont donc se ressembler)

Il permet de produire des échantillons à partir de distributions qui pourraient autrement être difficiles à échantillonner, et peut être utilisé pour décider de certaines valeurs proposées, même si nous ne connaissons pas la loi a posteriori.

Il génère des transitions d’état candidates à partir d’une distribution alternative, et accepte ou rejette chaque candidat de manière probabiliste. En d’autres mots, il génère des valeurs d’échantillon de manière que plus on produit de valeurs, plus la distribution de ces valeurs se rapproche de la distribution recherchée. Ces valeurs sont produites de manière itérative, la probabilité de l’échantillon à l’étape suivante ne dépend que de l’échantillon à l’étape considérée, la suite d’échantillons est donc une chaîne de Markov.

A chaque itération, l’algorithme choisit un candidat pour le prochain échantillon qui n’est basé que sur l’échantillon courant. Ensuite, le candidat est soit accepté avec une certaine probabilité (et dans ce cas il sera utilisé à l’étape suivante), soit rejeté avec la probabilité complémentaire (et dans ce cas l’échantillon actuel sera réutilisé à l’étape suivante). La probabilité d’acceptation est déterminée en comparant les valeurs d’une fonction f proportionnelle à la distribution-cible pour l’échantillon actuel, à la valeur de f pour l’échantillon candidat. Après un « grand nombre » d’itérations, la valeur des échantillons « perdent la mémoire » de l’état initial et suivent la distribution recherchée.

Nous allons appliquer cet algorithme afin de simuler la loi a posteriori
pour le même modèle de mélange \(\mathcal{M} = \left\{ \pi\phi(.;\mu_1,1)+(1-\pi)\phi(.;\mu_2,1) \right\}\) tel que:

De même loi a priori \(\mu_1 \sim \mathcal{N}(\delta_1,\frac{1}{\lambda}),~\mu_2 \sim \mathcal{N}(\delta_2,\frac{1}{\lambda})\) telle que:

Avec comme loi de proposition une loi normale centrée en la valeur courante des paramètres et de matrice de covariance diagonale \(\begin{pmatrix} \zeta^2 & 0 \\ 0 & \zeta^2 \end{pmatrix}\) pour \(\zeta^2 > 0\)

\(\underline{Algorithme}\)

\(\quad Initialisation:\)
\(\qquad\) Choix de valeur de départ arbitraire \(\mu_1^{(0)}~et~\mu_2^{(0)}\)
\(\quad Iteration~(t-1) \rightarrow t:\)
\(\quad\) 1. Simuler :
\(\qquad \tilde\mu_1 \sim \mathcal{N}(\mu_1^{(t-1)},\zeta^2)\)
\(\qquad \tilde\mu_2 \sim \mathcal{N}(\mu_2^{(t-1)},\zeta^2)\)
\(\quad\) 2. Avec probabilité : \(\frac{p(\tilde\mu_1, ~\tilde\mu_2 ~\mid~ x)}{p(\mu_1^{(t-1)}, ~\mu_2^{(t-1)} ~\mid~x )} \wedge 1\)
\(\qquad\) Accepter \(\tilde\mu_1~et~\tilde\mu_2\) et faire \(\mu_1^{(t)} = \tilde\mu_1~et~\mu_2^{(t)} = \tilde\mu_2\)
\(\qquad\) Sinon rejeter \(\tilde\mu_1~et~\tilde\mu_2\) et faire \(\mu_1^{(t)} = \mu_1^{(t-1)}~et~\mu_2^{(t)} = \mu_2^{(t-1)}\)

\(\underline{Simulations}\)

Nous allons aplliquer l’algorithme en une dimension sur le mélange suivant : \(0.3~\mathcal{N}(-2,1)+0.7~\mathcal{N}(5,1)\)
avec deux valeurs de départ : -2 et 5
et quatre valeurs de \(\lambda\) : 10, 1, 0.1 et 0.001
correspondant aux quatres valeurs de variance : 0.1, 1, 10, 100

Nous commençons avec n=\(10^4\) échantillons de Monte Carlo

Nous remarquons que la probabilité d’acceptation est très élévée au debut, qu’elle diminue ensuite, et qu’elle recroît sur la derniere simulation.

Pour mieux comprendre l’évolution de la probabilité d’acceptation, nous allons afficher celle-ci en fonction de chaque itération avec les paramètres précédents. : \(0.3~\mathcal{N}(-2,1)+0.7~\mathcal{N}(5,1)\)
avec deux valeurs de départ : -2 et 5, les meme que les parametres du melange En prenanr un nombre d’iterations egale à 100

On remarque que la probabilité d’acceptation est faible à chaque iterations.mais cela ne signifie pas pour autant que l’algorithme ne converge pas vers les bons paramètres simulés que nous cherchions au départ. la prueve dans l’histogramme suivant :

La plupart des valeurs générées par notre algorithme se situent autour des vraies espérances.

Voyons maintenant ce que nous avons obtenons avec n=\(10^6\) échantillons de Monte Carlo

Nous remarquons bien que la difficulté de cet algorithme est le calibrage de la proportion d’acceptation des paramètres proposés: Si \(\zeta^2\) est trop grand, presque toutes les étapes de l’algorithme MH seront rejetées. Si \(\zeta^2\) est trop petit, presque toutes les étapes seront acceptées.

Autrement dit, plus la variance de la densité de proposition est grande, moins la probabilité d’acceptation sera élévée.

Ce qui est comprehensible si notre vraisemblance est assez compliquée pour avoir plusieurs maxima locaux. L’algorithme va parcourir la courbe jusqu’à ce qu’il trouve des valeurs de probabilité maximale. Si la variance est trop courte, il va rester coincé sur les maxima locaux, car il sera toujours en échantillonnage des valeurs à proximité : l’algorithme croira donc qu’il est bloqué sur la distribution cible. Si la variance est trop grande, une fois qu’il se retrouve coincé sur un maximum local, il va rejeter les valeurs jusqu’à ce qu’il trouve d’autres régions de probabilité maximale. Si on propose la valeur au MAP (ou une autre région pareil, de probabilité maximale locale qui est plus grande que les autres), avec une grande variance, il va finir par rejeter presque toutes les autres valeurs parce que la différence entre cette région et les autres sera trop grande!

Quelle que soit la variance, tant que la probabilité de sélectionner la valeur de cette région maximale globale est positive, notre chaîne va converger… Pour eviter ce problème, on peut proposer différentes variations au cours d’une période de “chauffe” pour chaque paramètre et viser un certain taux d’acceptation satisfaisant. Sachant qu’un “bon” taux d’acceptation dépendra forcement de la forme de notre distribution a posteriori.

C’est le principe du recuit simulé. Lorsque l’on aplatit une surface, il est plus facile de s’y déplacer, alors que si on l’aiguise, il devient plus difficile de le faire sauf autour des pics.

Etant donné une densité \(\pi(x)\), on peut définir une densité associée \(\pi_{\alpha}(x) ∝ \pi(x)^{\alpha}\) pour \(\alpha\) > 0 assez grand. Ces distributions partagent toutes les mêmes modes. Lorsque \(\alpha\) > 1, la surface de \(\pi \alpha\) est plus contrastée que celle de \(\pi\) : Les pics sont plus hauts et les vallées plus basses. À l’inverse, abaisser \(\alpha\) à des valeurs inférieures à 1 rend la surface plus lisse en abaissant les pics et en relevant les vallées.

Nous pouvoir voir l’effet cumulatif d’une petite variance pour la proposition de marche aléatoire et une diminution de la puissance \(\alpha\) sur les representations suivantes: Nous avons des échantillons de \(10^4\) points, commencés au voisinage du mode parasite pour les distributions cibles \(\pi \alpha\) lorsque \(\alpha\) = 1, 0,1, 0,01 (de haut en bas), et la proposition est une marche aléatoire avec une variance de 0,1. (avec (1,3) comme point de départ).

Pour la véritable distribution cible, le premier tracé, 10 000 itérations de l’algorithme ne sont pas du tout suffisantes pour supprimer l’attraction du mode inférieur. Lorsque \(\alpha\) = 0,1, on peut espérer que quelques milliers d’itérations supplémentaires pourraient amener la chaîne de Markov vers l’autre mode. Pour \(\alpha\) = 0,01, seules quelques itérations suffisent pour changer de mode, étant donné que la selle entre les deux modes n’est pas beaucoup plus basse que les modes eux-mêmes.

Par rapport à l’échantillonneur de Gibbs, les algorithmes de Metropolis-Hastings sont des algorithmes MCMC prêts à l’emploi, car ils peuvent être réglés sur un éventail de possibilités beaucoup plus large. Pour une échelle suffisamment grande, l’algorithme de Metropolis-Hastings surmonte les inconvénients de l’échantillonneur de Gibbs.

Comparaison des algorithmes

Dans cette section, nous allons comparer les performances des deux algorithmes que nous avons implémentés en utilisant un mélange de deux gaussiennes comme données d’entrée. Les deux gaussiennes ont une moyenne de 0 et 10 respectivement et une variance de 1. Elles sont également proportionnelles à 0,3 et 0,7 respectivement. Nous utiliserons les mêmes paramètres initiaux pour mesurer les performances de chaque algorithme. Voici à quoi ressemble notre mélange de gaussiennes.

En supposant que nous connaissions déjà assez bien notre loi de mélange, nous allons utiliser les espérances des lois à priori initiales avec les mêmes valeurs que le mélange (soit 0 et 10). Cela constituera un test de l’efficacité en fixant lambda à 1 et en effectuant le même nombre d’itérations (100 itérations).

Nous allons afficher un histogramme qui représente les différentes valeurs d’espérance simulées et leurs fréquences, ainsi que l’évolution de la log-vraisemblance à chaque itération.

test1

Nous commencons par Hasting-Metrpolis.

Maintenant, voyons Gibbs.

interpretation

Au cours des premières parties, nous avons constaté que les deux algorithmes ont réussi à simuler correctement les paramètres des espérances des deux composantes du mélange et à générer des valeurs assez proches à chaque itération. Cependant, il y a quelques différences que l’on peut remarquer.

Tout d’abord, en regardant les histogrammes, nous voyons que pour l’algorithme de Metropolis-Hastings, la plupart des valeurs simulées sont concentrées autour des vraies valeurs des espérances du mélange, proches de 0 et 10, mais parfois il génère des valeurs un peu éloignées de ces valeurs, d’une variation d’environ 0,5. En revanche, l’algorithme de Gibbs semble être un peu plus précis, car bien qu’il génère moins de valeurs qui reflètent la vraie valeur, les espérances simulées sont très proches de cette dernière, avec une variation ne dépassant pas 0,1.

Ensuite, en regardant l’évolution de la vraisemblance, nous voyons que Metropolis-Hastings maximise la vraisemblance à posteriori au fil des itérations. Il semble atteindre le maximum après environ trente itérations et se stabilise autour de lpost = -2400. Pour l’algorithme de Gibbs, il semble que celui-ci maximise rapidement la log-vraisemblance, mais qu’ensuite, au fil des itérations, il recommence un nouveau cycle de maximisation. Nous constatons que la valeur du maximum est assez différente pour Gibbs, qui atteint à chaque cycle une valeur de lpost = -5700 environ.

test2

Pour avoir une idée plus claire des différences entre ces deux algorithmes, nous allons les tester de nouveau sur le même mélange, mais cette fois en changeant les paramètres de la loi à priori initiale. Nous allons prendre cette fois mu1_0 = -10 et mu2_0 = 0, en gardant toujours lambda = 1. Cela nous permettra de voir comment les algorithmes se comportent lorsque nous partons des paramètres de départ plus éloignés des paramètres du mélange.

Pour Hasting-Metrpolis

Pour gibs

Nous constatons que lorsque nous initialisons les espérances de départ à des valeurs plus éloignées des “vraies” valeurs, l’algorithme de Metropolis-Hastings simule au début de ces itérations des valeurs proches des valeurs initiales, puis il commence petit à petit à trouver les bons paramètres du mélange, à partir de la 20ème itération, où la log-vraisemblance à posteriori est maximisée et se stabilise.

Pour l’algorithme de Gibbs, même en changeant les valeurs de la loi à priori, son comportement reste globalement le même. La fonction de la log-vraisemblance à posteriori se maximise très rapidement et les espérances générées sont très proches des “vrais” paramètres du mélange.

Bonus : Comparaison avec EM

Dans cette section, nous allons comparer les résultats de ces deux algorithmes basés sur les chaînes de Markov avec l’algorithme de Expectation-Maximization (EM), qui est basé sur une technique de maximisation de la vraisemblance en utilisant une boucle de maximisation de l’espérance (E-step) et de maximisation (M-step). Pour ce faire, nous allons toujours utiliser le même mélange de données et toujours avec les paramètres initiaux éloignés des “vrais” paramètres, toujours en faisant 100 itérations.

Nous allons afficher l’histogramme des valeurs et la log-vraisemblance à posteriori.

Nous constatons que l’algorithme EM semble être plus efficace que les deux premiers algorithmes. En regardant les histogrammes, nous voyons qu’il a estimé seulement trois mauvaises valeurs pour les deux espérances sur 100 itérations. De plus, nous pouvons voir clairement que la log-vraisemblance à posteriori atteint le maximum beaucoup plus rapidement, après seulement cinq itérations.

Sources :
- Marin and Robert (2014)