Modèle de mélange

Densité de mélange

Considérons un mélange de $K$ lois gaussiennes dont la densité est donnée par $$f_\theta(x)=\sum_{k=1}^K\pi_kf_{\mathcal N(\mu_k,\sigma^2_k)}(x),$$ 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$.

Pour gérer les paramètres d'un mélange gaussien sous R nous allons travailler avec des listes. Exemple : Les paramètres du mélange suivant : $$f_\theta(x)= 0.2*f_{\mathcal N(2,4)}(x)+0.2*f_{\mathcal N(-5,1)}(x)+0.6*f_{\mathcal N(2,0.01)}(x)$$ sont stocké dans la liste suivante :

In [ ]:
theta = list(pi=c(.2,.2,.6), mu=c(2,-5,2), sig=c(2,1,.1))
theta

Remarquons que les trois vecteurs dans la liste theta ont la même longueur et que les $\pi_k$ et les $\sigma_k$ vérifient les conditions precisées ci-dessus.

Question 1

Ecrire une fonction drnormmix qui prend en argument un vecteur x et une liste de paramètres theta d'un mélange gaussien. Cette fonction évalue la densité du mélange gaussien associé aux paramètres theta en les éléments de x.

In [ ]:
dnormmix <- function(x,theta){
    
    return(dens)
}

Question 2

En changeant les paramètres d'un mélange gaussien, essayez de reproduire les six graphiques suivants, qui sont tous des mélanges gaussiens.

image

In [ ]:






















Question 3
    

  • Considérons les mélanges de loi uniforme $U[0,\lambda]$ dont la densité est donnée par 
$$f_\theta(x)=\sum_{k=1}^K\pi_kf_{U[0,\lambda_k]}(x),$$
où $0<\pi_k<1$, $\sum_{k=1}^K\pi_k=1$ et $\lambda_k>0$, $k=1,\dots,K$. Ecrire une fonction dunifmix pour évaluer la densité du mélange. Quelles formes de densité peut-on obtenir avec un mélange uniforme ?
    • 

  • Mêmes questions pour les mélanges de loi Gamma $\Gamma(\alpha,\beta)$ dont la densité est donnée par 
$$f_\theta(x)=\sum_{k=1}^K\pi_kf_{\Gamma(\alpha_k,\beta_k)}(x),$$
où $0<\pi_k<1$, $\sum_{k=1}^K\pi_k=1$ et $\alpha_k>0, \beta_k>0$ pour $k=1,\dots,K$.  Comparer aux formes des mélanges gaussiens.
    • 














In [ ]:



    
























Simulation de données
Question 4
    

  • Ecrire une fonction rnormmix pour générer des réalisations d'un mélange gaussien d'ordre $K$. Les arguments de la fonction sont la taille d'échantillon n et une liste theta avec les paramètres du mélange. Ne pas faire de boucle sur n.
    • 

  • Tester votre fonction en générant un grand échantillon d'un mélange gaussien dont on compare l'histogramme à sa densité par un graphique.
    • 














In [ ]:



    


rnormmix <- function(n,theta){

    return(obs)
}























Fonction de log-vraisemblance
Dans l'optique du calcul de l'estimateur de maximum de vraisemblance, considérons la fonction de log-vraisemblance d'un modèle de mélange  défini par
$$\ell(\theta)=\sum_{i=1}^n\log f_\theta(x_i),$$
où $f_\theta$ est la densité du mélange et  $(x_1,\dots,x_n)$ un échantillon du mélange $f_\theta$.


Question 5
Ecrire une fonction lvraisnorm qui prend en argument un échantillon observé obs et une liste  param de paramètres d'un mélange gaussien. La fonction renvoie la valeur de la log-vraisemblance du mélange gaussien associé.













In [ ]:



    


lvraisnorm <- function(param,obs){
    
    return(logvrais)
}























L'EMV est définit comme le point $\hat\theta$ où la log-vraisemblance $\ell(\theta)$ atteint son maximum. Dans des modèles de mélange, $\theta$ n'est jamais un scalaire, ce qui rend la visualisation de $\ell(\theta)$ compliquée.


Dans la question suivante nous allons considérer le mélange gaussien d'ordre 2 de densité
$$f_\theta(x) = p f_{\mathcal N(0,1)}(x) + (1-p) f_{\mathcal N(2.5,1)}(x).$$
Afin de faciliter la visualisation de la log-vraisemblance, nous supposons que les écarts-types  $\sigma_1=\sigma_2=1$ ainsi que les proportions $p$ et $1-p$ sont connus, et seulement les  moyennes $\mu_1$ et $\mu_2$ sont inconnus et à estimer à partir des données.
Dans ce cas, notons la fonction de log-vraisemblance  par
$$\ell(\mu_1,\mu_2)=
\sum_{i=1}^n\log\left\{ p f_{\mathcal N(\mu_1,1)}(x_i) + (1-p) f_{\mathcal N(\mu_2,1)}(x_i)\right\}.$$


Nous allons utiliser la fonction image pour la représentation graphique de $\ell(\mu_1,\mu_2)$. La syntaxe est la suivante :







image(x, y, z, zlim, xlim, ylim, col = heat.colors(12),
      add = FALSE, xaxs = "i", yaxs = "i", xlab, ylab,
      breaks, oldstyle = FALSE, useRaster, ...)












où x et y sont les grilles de l'abscisse et de l'ordonnée (des vecteurs de valeurs strictement croissantes de longueur r et s, resp.) et z est une matrice de taille r$\times$s contenant les valeurs de la fonction à tracer. Pour les autres options consulter les pages d'aide de R.


Exemple : Illustration de la fonction $(x,y)\mapsto y^x$ :













In [ ]:



    


x <- seq(0,2,by=.03)
y <- seq(0,1,by=.03)
M <- matrix(NA,nrow=length(x),ncol=length(y))
for (k in 1:length(x))
  M[k,] <- y^x[k]
image(x,y,M,col=grey(0:255/255))























On peut ajouter des lignes de niveau par la fonction contour. Ajouter aux instructions ci-dessus la commande :


contour(x,y,M,add=T,nlevels=10,col='yellow')




















Question 6
    

  • Utiliser la fonction image pour visualiser la fonction de log-vraisemblance 
$\ell(\mu_1,\mu_2)$ pour un échnatillon simulé de taille 1000 du mélange gaussien donné ci-dessus avec $p=0.5$. Constater que la log-vraisemblance n'est pas concave. Combien de maxima locaux admet-elle ? Commenter.
    • 

  • Même question lorsque $p=0.7$. 
    • 

  • Supposons que l'on applique la méthode de Newton-Raphson pour maximiser la fonction $(\mu_1,\mu_2)\mapsto \ell(\mu_1,\mu_2)$. Va-t-elle toujours converger vers l'EMV quand $p=0.5$ ? Et  quand $p=0.7$ ? 
    • 














In [ ]: