Dans ce TP vous allez mettre en oeuvre l'algorithme EM pour un mélange gaussien de densité $$f_\theta(x)=\sum_{k=1}^K\pi_kf_{\mathcal N(\mu_k,\sigma^2_k)}(x),\qquad x\in\mathbb R,$$ où $0<\pi_k<1$, $\sum_{k=1}^K\pi_k=1$, $\mu_k\in\mathbb R$ et $\sigma_k>0$, $k=1,\dots,K$.
Notons les paramèters du modèle par $$\theta=\left(\pi_1,\dots,\pi_K,\mu_1,\dots,\mu_K,\sigma_1,\dots,\sigma_K\right)$$ Comme au dernier notebook, sous R nous allons stocker les paramètres dans une liste :
theta = list(pi=c(3,2,5)/10, mu=c(4,-5,0), sig=c(1.3,1,1.5))
theta
Soient $\mathbf x=(x_1,\dots,x_n)$ un échantillon du mélange gaussien de densité $f_\theta$.
L'algorithme EM est un algorithme itératif pour approcher l'EMV en tirant profit de la structure latente du modèle. Il consiste à calculer succéssivmeent pour $t=0,1,\dots$ \begin{align} \theta^{(t+1)}=\arg\max\theta Q(\theta|\theta^{(t)}), \end{align} où $Q(\theta|\theta')$ est l'espérance conditionnelle de la fonction de log-vraisemblance des données complètes sous le paramètre $\theta'$, et $\theta^{(0)}$ est un point initial choisi par l'utilisateur. En cours, nous avon montré que $\theta^{(t+1)}$ est donné par \begin{align} \pi_k^{(t+1)}&= \frac1n\sum{i=1}^n\alpha{k|i}^{(t)},\quad k=1,\dots,K \ \mu_k^{(t+1)}& =\frac1{n\pi_k^{(t+1)}}{\sum{i=1}^nxi\alpha{k|i}^{(t)}},\quad k=1,\dots,K \qquad\qquad\qquad\qquad(*)\ \sigmak^{(t+1)} &= \sqrt{\frac1{n\pi_k^{(t+1)}}{ \sum{i=1}^n\alpha_{k|i}^{(t)}(x_i-\mu_k^{(t+1)})^2}},\quad k=1,\dots,K \end{align}
avec \begin{align} \alpha{k|i}^{(t)} = \mathbb P{\theta^{(t)}}(Z=k|X=xi) = \frac{\pi_k^{(t)}f{\mathcal N(\muk^{(t)},\sigma_k^{2(t)})}(x_i)}{f{\theta^{(t)}}(x_i)},\quad k=1,\dots,K,i=1,\dots,n . \end{align}
update.theta <- function(obs, theta){
return(theta.new)
}
Ecrire une fonction algoEM qui effectue l'algorithme EM. Ses arguments sont un vecteur d'observations obs et une liste avec les valeurs initiales des paramètres. Cette fonction retourne une liste avec les paramètres estimés par l'algorithme EM et le nombre d'itérations effectuées. Plus précisément, la fonction algoEM itère des mises à jours de paramètres par les formules en $(*)$ jusqu'à convergence. On utilise comme critère d'arrêt $$\left|\frac{\ell(\theta^{(t+1)})-\ell(\theta^{(t)})}{\ell(\theta^{(t)})}\right|<\varepsilon,$$ où $\ell(\theta)$ désigne la fonction de log-vraisemblance et avec $\varepsilon=10^{-3}$ comme seuil. Pour éviter que l'algorithme tourne indéfiniment, on s'arrête au plus tard après $R=200$ itérations. On peut déclarer $\varepsilon$ et $R$ comme des arguments de la fonction avec des valeurs par défaut.
Testez votre fonction sur des données simulées en initialisant avec les vraies valeurs des paramètres. Visualiser la densité du mélange estimée par l'algorithme EM en la superposant à l'histogramme des données. Rajouter également la densité du mélange avec les vraies valeurs des paramètres.
algoEM <- function(obs, theta.init){
resultat <- list(emv = theta.new, nb.iter = it)
return(resultat)
}
Essayer différents points initiaux. Est-ce que l'algorithme converge toujours vers le même point ? Quel effet sur le nombre d'itérations. Comme points initiaux essayer
Commenter.
On a vu qu'il est surtout important de bien initialiser les moyennes $\mu_k$. Or, les observations $x_i$ d'un mélange gaussien prennent des valeurs près des moyennes $\mu_k$. On peut en déduire l'initialisation suivante :
Afin d'ameliorer la performance, on peut appliquer l'algorithme EM avec plusierus (disons $L$) points initiaux différents et retenir parmi les $L$ solutions $\hat\theta_1,\dots,\hat\theta_L$ celle avec la plus grande valeur de la fonction de log-vraisemblance. Autrement dit, la solution finale $\hat\theta$ est donnée par $\hat\theta\in\{\hat\theta_1,\dots,\hat\theta_L\}$ tel que $$\ell(\hat\theta) \geq\ell(\hat\theta_l),\qquad l=1,\dots,L$$
Ecrire une fonction EMgauss qui performe la procédure décrite ci-dessus. Cette fonction prend en argument l'échantillon obs, le nombre de composants K et le nombre de points initiaux souhaités L.
EMgauss <- function(obs,K,L){
return(theta)
}
On reprend les données réelles vues en cours :
Les données sont disponibles dans le fichier donneesmelange.RData que vous pouvez charger par l'instruction suivante :
load('donneesmelange.RData')
table(ailes)
table(chlor)
Pour chaque jeu de données :