Considèrons un mélange gaussien d'ordre $K$ dans un contexte de statistique bayésienne.
Notons les paramètres d'un mélange gaussien par $$\theta=(\underline \pi,\underline \mu,\underline\sigma^2)=\left(\pi_1,\dots,\pi_K,\mu_1,\dots,\mu_K,\sigma_1^2,\dots,\sigma_K^2\right),$$ où $0<\pi_k<1$, $\sum_{k=1}^K\pi_k=1$, $\mu_k\in\mathbb R$ et $\sigma_k^2>0$, $k=1,\dots,K$.
On considère le modèle suivant: \begin{align} \theta&\sim\pi(\theta)\ zi|\theta&\sim \sum{k=1}^K\pik\delta{{k}},\quad i=1,\dots,n\ xi|z_i,\theta&\sim \mathcal N(\mu{zi},\sigma^2{z_i}),\quad i=1,\dots,n. \end{align}
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{\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 de densité $$f(x)=\frac{\beta^\alpha}{\Gamma(\alpha)}x^{-\alpha-1}\exp\left\{-\frac\beta x\right\},\qquad x>0.$$ L'ensemble des hyperparamètres est $$ \psi =(\gamma_1,\dots,\gamma_K, \alpha_1,\dots,\alpha_K,\lambda_1,\dots,\lambda_K,\beta_1,\dots,\beta_K), $$ où $\gamma_k>0$, $\alpha_k\in\mathbb R$, $\lambda_k>0$ et $\beta_k>0$.
Pour simuler de la loi de Dirichlet on utilisera la fonction rdirichlet(n,gamma) du package gtools que l'on charge par la commande :
library(gtools)
Pour simuler de la loi Gamma inverse, on note la propriété suivante : $$X\sim\mathcal{IG}(\alpha,\beta)\Longleftrightarrow \frac1X\sim\Gamma(\alpha,\beta).$$
psi <- list(gamma=c(1,1,1), alpha=c(1,1,1), lambda=c(2,2,2), beta=c(1,1,1))
rthetaB <- function(psi){
return(theta)
}
rnormmixB <- function(n,psi){
return(list(obs=obs,theta=theta))
}
La loi a posterior de $\theta$ sachant $\mathbf x=(x_1,\dots,x_n)$ est donnée par $$\pi(\theta|\mathbf x)\propto \prod_{i=1}^n\left(\sum_{k=1}^K\pi_kf_{\mathcal N(\mu_k,\sigma_k^2)}(x_i) \right)\times\pi(\theta) $$ En en déduit que la fonction de log-aposteriori est donnée par $$\log \pi(\theta|\mathbf x) = \sum_{i=1}^n\log\left(\sum_{k=1}^K\pi_kf_{\mathcal N(\mu_k,\sigma_k^2)}(x_i) \right)+\log(\pi(\underline\pi)) +\sum_{k=1}^K\log(\pi(\mu_k|\sigma^2_k))+ \sum_{k=1}^K\log(\pi(\sigma_k^2))+\text{cst}, $$ où $\text{cst}$ est une constante que l'on ne précisera pas.
Remarquons que pour $X\sim\mathcal{IG}(\alpha,\beta)$ et $Y\sim\Gamma(\alpha,\beta)$, les densités associées vérifient $$f_X(s) = \frac1{s^2}f_Y\left(\frac1s\right),\qquad s>0.$$
psi <- list(gamma=c(3,3), alpha=c(0,4), lambda=c(4,4), beta=c(2,2))
logaposteriori <- function(theta,x,psi){
return(l)
}
Le Gibbs sampler exploite la structure des variables latentes $\mathbf z=(z_1,\dots,z_n)$. Plus précisément, l'algorithme consiste à simuler successivement des lois conditionnelles suivantes : \begin{align} \pi(\mathbf z|\mathbf x,\theta)&=\prod{i=1}^n\pi( z_i|x_i,\theta)\qquad\text{ avec }\qquad\ \pi( z_i|x_i,\theta)\propto \pi{zi}f{\mathcal N(\mu{z_i},\sigma{zi}^2)}(x_i) \ \pi(\theta|\mathbf z,\mathbf x)&\propto \pi(\underline \pi|\mathbf z,\mathbf x) \times\pi(\underline\mu,\underline\sigma^2|\mathbf z,\mathbf x)\qquad\text{ avec }\ &\quad\underline \pi|(\mathbf z,\mathbf x)\sim\mathcal D(\gamma_1+n_1,\dots,\gamma_K+n_k)\qquad \text{où}\qquad n_k=#{i:z_i=k}\ &\quad\pi(\underline\mu,\underline\sigma^2|\mathbf z,\mathbf x) \propto\prod{k=1}^K\pi(\muk,\sigma^2_k|\mathbf z,\mathbf x) =\prod{k=1}^K\left( \pi(\muk|\sigma_k)\pi(\sigma_k)\prod{i:zi=k}f{\mathcal N(\mu_k,\sigma_k^2)}(x_i)\right) \end{align} Plus précisément, on peut montrer que \begin{align} \muk|\sigma_k,\mathbf z,\mathbf x &\sim\mathcal N\left(\frac{\sum{i:zi=k}x_i+\alpha_k\lambda_k}{n_k+\lambda_k},\frac{\sigma_k^2}{n_k+\lambda_k}\right)\ \sigma_k^2|\mathbf z,\mathbf x &\sim\mathcal{IG}\left(\frac{n_k+\lambda_k}2,\frac12\left(\sum{i:zi=k}x_i^2+\lambda_k\alpha_k^2+\beta_k-\frac{(\sum{i:z_i=k}x_i+\alpha_k\lambda_k)^2}{n_k+\lambda_k} \right) \right) \end{align}
rz <- function(theta,obs){
return(z)
}
rtheta <- function(z,obs,psi){
return(theta)
}
gibbs <- function(obs,theta.init,psi,R=1000,B=500){
return(list(theta=theta,z=z))
}
Le but du Gibbs sampler est de créer un échantillon de la loi aposteriori $\pi(\theta|\underline x)$. Ensuite, cet échantillon de valeurs simulées $\theta^{(t)},t=1,\dots,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=\frac1R\sum_{t=1}^R\mu_1^{(t)}.$$
Un outil simple pour apprécier la convergence de l'algorithme consiste à surveiller l'évolution de la moyenne des valeurs simulées. Autrement dit, on considère les tracés du type $$r\mapsto\frac1r\sum_{t=1}^r\mu_1^{(t)}\qquad\text{ pour }r=1,\dots,R.$$
Le Gibbs sampler fournit également un échantillon simulé des étiquette $z_i^{(t)}$. On peut étudier pour chaque observation $x_i$ la proportion avec laquelle le Gibbs sampler l'affecte au premier (deuxième,...) groupe, c'est-à-dire la proportion $$p_{i,k}=\#\{t:z_i^{(t)}=k\}/R\qquad\text{ pour }k=1,\dots,K,\quad i=1,\dots,n.$$ Dans le cas d'un mélange gaussien, on considère que l'algorithme a convergé si les tracés des estimateurs du type $r\mapsto\frac1r\sum_{t=1}^r\mu_1^{(t)}$ se sont stabilisés et si l'affection des observations $x_i$ a un groupe est stable.
n <- 2000
psi <- list(gamma=c(3,3), alpha=c(0,4), lambda=c(4,4), beta=c(2,2))
set.seed(1711)
rthetaB(psi)
res <- rnormmixB(n,psi)
sim <- gibbs(res$obs,res$theta,psi)
tracer.estim <- function(sim){
return(estim)
}
class.zi <- function(sim){
return(class)
}
Pour appliquer l'échantillonneur de Gibbs, il faut choisir des valeurs des hyperparamètres. Nous allons faire le choix suivant pour un mélange de deux gaussiennes :
Et pour initialiser $\theta$, on utilisera :
Appliquer le Gibbs sampler aux données ailes et chlor disponible dans le fichier donneesmelange.RData. Contrôler la convergence de l'algorithme à l'aide des tracés des estimateurs bayésiens en fonction du nombre d'itérations. Représenter graphiquement les probabilités $p_{i,1}$ des observations d'appartenir au premier groupe du mélange. Comparer les résultats à ceux obtenus par l'algorithme EM.