Modèles à structure latente

TP2 Modèles de mélange, Model-Based Clustering, algorithme EM

Modèles de mélange

Densité d’une loi de mélange

Nous allons tracer plusieurs densités de loi de mélange gaussiens en dimension 1, en faisant varier leurs paramètres.

Simulations d’échantillons

\(\bullet~\underline{Exemple~1}\)
Nous allons simuler des échantillons de taille 200
Pour la loi de densité \(f(.;(\frac{1}{3},\frac{1}{6},\frac{1}{2},0,5,10,1,1,4)) = \frac{1}{3}\phi(.;0,1) + \frac{1}{6}\phi(.;5,1) + \frac{1}{2}\phi(.;10,4)\)

\(\bullet~\underline{Exemple~2}\)
Nous allons maintenant simuler des échantillons de taille 1 000
Pour la loi de densité \(f(.;(\frac{1}{3},\frac{1}{3},\frac{1}{3},-2,0,2,1,1,1)) = \frac{1}{3}\phi(.;-2,1) + \frac{1}{3}\phi(.;0,1) + \frac{1}{3}\phi(.;2,1)\)

\(\bullet~\underline{Exemple~3}\)
Nous allons, pour finir, simuler des échantillons de taille 1 000
Pour la loi de densité \(f(.;(\frac{3}{4},\frac{1}{8},\frac{1}{8},10,0,-10,1,2,3)) = \frac{3}{4}\phi(.;10,1) + \frac{1}{8}\phi(.;0,2) + \frac{1}{8}\phi(.;-10,3)\)

Conservation des vraies composantes

\(\bullet~\underline{Exemple~1}~~~f(.;(\frac{1}{3},\frac{1}{6},\frac{1}{2},0,5,10,1,1,4)) = \frac{1}{3}\phi(.;0,1) + \frac{1}{6}\phi(.;5,1) + \frac{1}{2}\phi(.;10,4)\)

Nous regardons maintenant ce que l’on obtient si nous multiplions ces composantes par un certain \(\lambda\), correspondant a la multiplication de chaque composante par leurs probabilités.

\(\bullet~\underline{Exemple~2}~~~f(.;(\frac{1}{3},\frac{1}{3},\frac{1}{3},-2,0,2,1,1,1)) = \frac{1}{3}\phi(.;-2,1) + \frac{1}{3}\phi(.;0,1) + \frac{1}{3}\phi(.;2,1)\)

\(\bullet~\underline{Exemple~3}~~~f(.;(\frac{3}{4},\frac{1}{8},\frac{1}{8},10,0,-10,1,2,3)) = \frac{3}{4}\phi(.;10,1) + \frac{1}{8}\phi(.;0,2) + \frac{1}{8}\phi(.;-10,3)\)

Comme la premiere composante a un “gros poids” face aux deux autres composantes, la densité du mélange est plus haute en la moyenne de celle-ci.

Comparaisons par la règle du MAP

Nous allons comparer les «vraies» classes et celles obtenues, par la règle du MAP calculée avec les vrais paramètres du mélange.

\(\bullet~\underline{Exemple~1}\)

MAP Estimate

Parameter | MAP_Estimate
------------------------
comp1     |         0.38
comp2     |         5.18
comp3     |         9.89
mixture   |         0.33

Nous pouvons voir que le maximum de vraissemblance pour la mixture est proche de celle de la premiere composante.

\(\bullet~\underline{Exemple~2}\)

MAP Estimate

Parameter  | MAP_Estimate
-------------------------
comp1bis   |        -1.88
comp2bis   |         0.12
comp3bis   |         1.64
mixturebis |         0.10

Nous pouvons voir que le maximum de vraissemblance pour la mixture est proche de celle de la deuxieme composante.

\(\bullet~\underline{Exemple~3}\)

MAP Estimate

Parameter  | MAP_Estimate
-------------------------
comp1ter   |         9.70
comp2ter   |        -0.18
comp3ter   |        -9.72
mixtureter |        -0.06

Nous pouvons voir que le maximum de vraissemblance pour la mixture est proche de celle de la deuxieme composante.

Vraissemblance

Nous avons ici illustrer la likelihood en fonction de \(\sigma_1\).
A noter que notre likelihood est modifiée, en effectuant une transformation exponentielle, pour eviter des valeurs trop basses numériquement.
On voit que le maximum est bien atteint à la bonne valeur de \(\sigma_1\), c’est à dire \(\sigma_1 = 1\) qui est la bonne valeur théorique.

Situation en dimension 2

Comme demandé, nous choisissons un mélange de deux composantes bien séparées dans le modèle [\(p\) _ \(L_k\) _ \(B_k\)].

Pour cela, nous prennons:
 - Avec comme proportion \(p_1=1/4\), une loi normale centrée reduite \(\mathcal{N}(0,1)\), pour les deux dimensions
 - Avec comme proportion \(p_2=3/4\), une loi normale \(\mathcal{N}(10,2)\), pour les deux dimensions

On note:
 - \(p\) la proportion de la deuxième composante, soit \(p=3/4\)
 - \(\lambda A\) la décomposition de sa matrice de covariances

1. Loi “originale”

Voici les isodensités de cette loi, ainsi qu’un échantillon de taille 1000:

Nous pouvons très bien voir les deux classes (les deux lois bi-dimensionnelles) sur ce plot 3D, qui montre la densité du mélange en deux dimensions.
Nous remarquons que les isodensités sont plutôt de forme cylindrique (surtout celle au sommet des gaussiennes).

2. On double \(\lambda\)

On représente quelques isodensités, ainsi qu’un échantillon de taille 1 000 de la loi obtenue:

Nous remarquons que les isodensités sont maintenant plutôt de forme elliptique.

3. On divise la deuxieme proportion par 2

On représente quelques isodensités, ainsi qu’un échantillon de taille 1 000 de la loi obtenue (avec la valeur d’origine de \(\lambda\)):

Nous remarquons comme une rotation de 180 degrès, nous avons de belles isodensités bien rondes.

Algorithme EM, dimension 1

Implémentation de l’algorithme

1. Algorithme EM dans ce cadre

Initialisation: paramètre \(\theta\)
Répéter:
-Étape E :
Pour chaque point de données \(x_i\), calculer l’attente de la variable latente \(z_i\) en fonction des paramètres du modèle actuel \(\theta\): \(q_i(z_i) = p(z_i \mid x_i, \theta)\)
-Étape M :
Mettre à jour les paramètres du modèle \(\theta\) en utilisant l’étape de maximisation:
\(\theta = \arg\max_{\theta} \sum_{i} \sum_{z_i} q_i(z_i) \log p(x_i, z_i \mid \theta)\)
Jusqu’à: Convergence

2. Codage de l’algorithme EM

Nous allons implementer l’algorithme EM dans notre cas, qui est un melange de gaussienne en dimension 1, d’abord pour un melange de deux composantes, puis on verra pour des composantes quelconques

Nous allons prendre une mixture gaussienne de deux composantes avec les parametres suivants :
mean = c(1,1), sd = c(0,10) et proportion = c(0.4,0.6)

Nous avons fixé les parametres initiaux arbitrairement:

[1] 0.392005 0.607995 2.046043 8.690435 3.718406 5.529920

La sortie est de l’ordre sivant : (p1,p2,mu1,mu2,sigma1,sigma2)

On remarque que l’estimation de l’algorithme n’est pas loin des ‘vraies’ valeurs des parametres du melange, mais un peu loin de la porportion.

Implementation:
Maintenant nous allons implementer la fonction qui prend en parametre un nombre de composantes K

Nous avons essayé de faire le maximum pour que l’algorithme soit le plus optimal possible, il s’arrete quand l’erreur de convergence de la difference entre les nouvelles valeurs du log de vraissemblance et l’ancienne ne depassent pas 0.0001, sinon apres 1000 iterations.

\(\underline{Exemple~avec~deux~composants}\)

Nous allons tester l’implementation EM sur un echantillion d’un melange de deux composantes de taille 1000, avec une proportion egale, une gaussienne centre reduite pour la premiere composante, et une moyenne de 20 et variance de 10 pour la deuxieme.

[1] "p1_vrai: 0.5             p1_estimé: 0.462807705520327"

[1] "p2_vrai: 0.5             p2_estimé: 0.537192294479673"

[1] "mu1_vrai:0               mu1_estimé:22.1747686327561"

[1] "mu2_vrai:20              mu2_estimé: 0.0501854938347355"

[1] "sigma1_vrai:1            sigma1_estimé: 8.13312558501076"

[1] "sigma2_vrai:10           sigma2_estimé: 0.918765566593624"

Nous pouvons voir que pour cet exemple, le EM arrive a se debrouiller pas mal, et il retrouve les parametres et la proportion des deux classes.

\(\bullet\) Fonction qui renvoie la densité graphiquement et compare avec la “vraie” densité des données

Maintenant nous allons implenter une fonction qui affiche la densité du melange ajusté avec les parametres estimés, et pour bien visualiser le resultat, nous allons aussi afficher la densité du ‘vrai’ melange

Comme nous le voyons sur le plot, la densité estimée avec le EM, colle presque parfaitement et couvre les deux ‘vraies’ classes de depart.

\(\bullet\) La classification des observations par MAP

MAP Estimate: 0.12

MAP Estimate: -0.03

Nous pouvons voir qu’avec la classification par maximum a posteriori nous donne un resultat assez similaire pour les deux melanges gaussiens, une coupure autour de la moyenne en 0, car la variance est faible pour la premiere composante par raport à la deuxieme, et donc le maximum se trouve tres probalement vers la moyenne de cette classe.

Application sur des données simulées

1. Suivi de l’évolution de la log-vraisemblance au cours des itérations

Dans la fonction EM que nous avons implementé, nous avons mis en place un parametre qui stocke l’evolution du maximum de vraissemblance au cours des iterations, nous allons aficher cela par la suite : Pour 100 iterations, le log-lik trouve son maximum apres 2 iterations et se stabilise tres rapidement.

2. Algorithme avec plusieurs paramètres initiaux différents

Exemple avec un melange de trois gaussiennes : mean = c(-5,5,12), sd = c(1,4,8), et proportion = c(0.2,0.3,0.5)

[1] "p1_vrai: 0.3             p1_estimé: 0.3"

[1] "p2_vrai: 0.3             p2_estimé: 0.280000000028213"

[1] "p3_vrai: 0.4             p3_estimé: 0.419999999971787"

[1] "mu1_vrai:-20               mu1_estimé:-20.3475446147184"

[1] "mu2_vrai:0              mu2_estimé: -0.580803629290619"

[1] "mu3_vrai:20              mu3_estimé: 20.101847904426"

[1] "sigma1_vrai:1            sigma1_estimé: 1.10058143882416"

[1] "sigma2_vrai:2           sigma2_estimé: 2.01444044069603"

[1] "sigma3_vrai:2           sigma3_estimé: 2.27244424589851"

3. Densité ajustée

Notre implementation arrive a faire presque les memes classes que la densité originale!
Les estimations etaient vraiment bonnes

4. Comparaison des classifications obtenues par MAP

MAP Estimate: -19.87

MAP Estimate: -20.41

Quand nous avons agmenté le nombre de classes, le log-lik augmente cette fois au fur et a mesures pour qu’il retrouve son max au bout de 15 iterations

\(\underline{Exemple~2}\)

Exemple avec un melange de quatres gaussiennes :
mean = c(-20,-10,0,10), sd = c(1,2,1,2), et proportion = c(0.2,0.2,0.3,0.3)

Nous prenons toujours des moyennes eloignées pour chaque classe pour qu’on puisse les distinguer sur les representations graphiques

[1] "p1_vrai: 0.2             p1_estimé: 0.337557580899849"

[1] "p2_vrai: 0.2             p2_estimé: 0.204691962082917"

[1] "p3_vrai: 0.3             p3_estimé: 0.252271641019615"

[1] "p4_vrai: 0.3             p2_estimé: 0.20547881599762"

[1] "mu1_vrai:-20               mu1_estimé:9.87296588852944"

[1] "mu2_vrai:-10              mu2_estimé: -19.9574281235813"

[1] "mu3_vrai:0              mu3_estimé: -0.18287481880476"

[1] "mu4_vrai:10              mu4_estimé: -9.79855914932671"

[1] "sigma1_vrai:1            sigma1_estimé: 2.00076380868719"

[1] "sigma2_vrai:2           sigma2_estimé: 0.895164841059247"

[1] "sigma3_vrai:1           sigma3_estimé: 0.955463140384335"

[1] "sigma4_vrai:2           sigma4_estimé: 2.13915575322573"

\(\bullet\) Representation graphique du deuxieme exemple

Meme avec 4 composantes, l’algo EM arrive trés bien a classer les observations dans leurs bonnes classes

5. Trouver un cas de croissance avec palier de la log-vraisemblance lors des itérations de EM.

Cette fois avec 4 classes, la vraissemblance prends plus de temps pour trouver son maximum, ce qui donne une sorte de fonction croissante

Algorithme EM, dimension quelconque (Mixmod)

Rmixmod sur des données simulées

Nous allons comparer graphiquement la densité du mélange estimé à la vraie densité de l’échantillon

Voici notre mélange gaussien, avec le même nombre d’échantillons pour chaque gaussienne:
\(\text{comp}_1 \simeq \mathcal{N}(0,1)\)
\(\text{comp}_2 \simeq \mathcal{N}(5,1)\)
\(\text{comp}_3 \simeq \mathcal{N}(10,4)\)

Maintenant nous allons utiliser Rmixmod et son algorithme EM pour ce mélange gaussien.

[1] 1

mixmodCluster nous donne les clusters suivants:
Cluster 1 :
-proportion = 0.1555
-means = 4.4874
-variances = 1.5376
Cluster 2
-proportion = 0.3477
-means = -0.0805
-variances = 1.1103
Cluster 3
-proportion = 0.4968
-means = 9.8012
-variances = 19.2738

C’est à dire, il nous donne comme mélange:
\(0.1555\mathcal{N}(4.4874 ,1.5376 )+0.3477\mathcal{N}(-0.0805 ,1.1103 )+ 0.4968\mathcal{N}(9.8012,19.2738)\)

On voit donc que l’algorithme donne des valeurs assez proches, bien que non exactes.
Cependant, c’est assez normal car nous faisons un mélange avec un échantillon total n=200.

\(\bullet\) Illustration de Rmixmod avec des gaussiennes 2D

Illustration du package Rmixmod pour un modèle de mélange 2D avec les vecteurs d’ésperance suivants:
\[\mu_1=(0, 0), \mu_2= (3, 3) , \mu_3=(1.5, 6)\] et les matrices de covariances suivantes:
\[\sigma_1 = \text{Diag}(1, 0, 0, 1), \sigma_2 = \text{Diag}(1, 0.5, 0.5, 1), \sigma_3 = \text{Diag}(1, 0.2, 0.2, 1) \]
et de proportions respectives: \((0,3, 0,2, 0,5)\).

[1] 1

[1] 2

Output du modèle:
-nbCluster = 3
-model name = Gaussian_pk_Lk_C
-criterion = BIC(7646.5941)
-likelihood = -3778.3967
- proportions: (0.4925, 0.1956, 0.3120)
- moyennes: (1.4542, 6.0012 ), (3.0599, 3.1260 ), (0.0061, 0.0712 )
variances :
| 1.0257 0.1967 | | 0.8607 0.1650 | | 0.9987 0.1915 |
| 0.1967 1.0832 | | 0.1650 0.9089 | | 0.1915 1.0546 |

On voit que l’algorithme a trouvé de très bonnes approximations des proportions et esperances.
De même les resultats obtenus pour les matrices de covariances sont très bons, avec une erreur plus prononcé pour le deuxième cluster.

\(\bullet\) Illustration de Rmixmod avec des gaussiennes 3D

Pouvant utiliser l’algorithme EM de Rmixmod en diemnsion quelconque, on fera une dernière illustration en 3d.
En effet, une illustration en dimension supérieure ne serait pas très pertinente visuellement.
Notre illustration ce fait avec quatre gaussiennes 3d, tous de matrice de convariance égale à l’identité, et avec les esperances suivantes: \[ \mu_1=(0, 0, 0),\mu_2=(2, 2, 2), \mu_3=(-2, -2, -2), \mu_4=(1, 0, -2)\]

[1] 1

[1] 2

[1] 3

Output de l’algorithme:
-nbCluster = 4
-model name = Gaussian_pk_Lk_C
-criterion = BIC(10644.7285)
-likelihood = -5239.4712
- moyennes: (0.9661,-0.0617,-1.8577 ), (2.0548,2.0647,1.9462 ), (0.0430,0.2588,0.0972)
- covariances: toutes quasi-diagonales (1,1,1)

On peut observer que les matrices de covariance diagonales sont relativement bien respectés, mais avec des variances inexactes (notamment le cluster 3) que les moyennes ont été très bien approximés, de même pour le nombre de cluster qui est exacte.
On notera néanmoins que cet exemple jouet est relativement facile et bien choisi pour l’algorithme.

Rmixmod avec nbCluster = 2:8

Revenons à notre problème initiale pour Rmixmod, c’est à dire le mélange de gaussienne 1D suivant:
Même proportion pour chaque gaussienne et:
\(\text{comp}_1 \simeq \mathcal{N}(0,1)\) \(\text{comp}_2 \simeq \mathcal{N}(5,1)\) \(\text{comp}_3 \simeq \mathcal{N}(10,4)\)

Jouons avec les arguments d’entrée de Rmixmod

\(\underline{Exemple~1}\)

Nous simulons les trois échantillons gaussiens suivant:
  - un échantillon de taille 125, pour la loi de densité \(\mathcal{N}(0.1,0.2)\)
  - un échantillon de taille 200, pour la loi de densité \(\mathcal{N}(0.8,0.2)\)
  - un échantillon de taille 160, pour la loi de densité \(\mathcal{N}(1.4,0.4)\)
avec les proportions respectives suivantes: \((0.258, 0.412, 0.330)\).

\(\bullet\) Variation de l’argument nbCluster

[1] 1

Les trois gaussiennes trouvées sont: \[ \mathcal{N}(0.7774 ,0.0620), \mathcal{N}(0.0856 ,0.0389), \mathcal{N}(1.4066 ,0.1759)\] avec comme proportions respectives: \(0.4315, 0.2446, 0.1759\)
La likelihood est de: -403.8474.
Le criètre BIC est de: 857.1679
Les données trouvées sont très proches des données simulées.

[1] 1

Bien evidemment, le critère BIC est plus élevé avec un nombre de clusters fixé égale à 5 et la likelihood est pratiquement la même:
La likelihood est de: -401.4026
Le criètre BIC est de: 889.3833
Ce résultat est cleui qu’on attendait.
En effet, on a fait exprès de demander un nombre de gaussiennes plus grand que le vértiable nombre.

\(\bullet\) Variation de l’argument models

Nous fixons désormais le nombre de cluster a 3.

[1] 1

Quand on change sur les modeles Gaussiennes de proportions libres, on trouve que le modele est Gaussian_pk_Lk_Bk et le BIC diminue de 811.8644 a 781.5433 ca veut dire que ce modele est plus precis.

[1] 1

Dans ce cas, on trouve que le modele est Gaussian_p_Lk_Ck.
Le critère BIC maintenant est de 800.3081.

[1] 1

(Les gaussiennes sont en proportions égales)
On trouve que dans ce cas, le modele est Gaussian_p_Lk_D_Ak_D, et le Bic est de 812.7693

[1] 1

(Les gaussiennes sont en proportions égales)
Dans ce cas, on trouve que le modele est Gaussian_p_Lk_C et le BIC est de 770.9846 (Le Bic plus bas).
Cela nous dit que ce modele est le meilleur pour predire ces donnees.

[1] 1

En donnant en argument 1:8 à mixmodCluster, c’est à dire en lui indiquant que nous voulons un modèle de mélange avec entre 1 et 8 cluster, l’algorithme essaye de trouver le nombre de cluster ressemblant le plus à notre mélange.
Le critère d’information bayésienne BIC est utiliser par défaut par l’algorithme.
L’utilisation de ce critère permet à l’algorithme de ne pas surestimer le nombre de cluster, en effet BIC pénalise un trop gros nombre nombre de clusters.

\(\underline{Exemple~2}\)
Notre illustration de Rmixmod est tiré d’un modèle de mélange simple avec trois gaussiennes. Nous allons maintenant donner un problème un peu plus difficile à Rmixmod, et voir ce que l’algorithme donne comme résultat, avec sept gaussiennes.

Nous simulons les septs échantillons gaussiens suivant:
  - un échantillon de taille 60, pour la loi de densité \(\mathcal{N}(1,0.2)\)
  - un échantillon de taille 60, pour la loi de densité \(\mathcal{N}(0,0.3)\)
  - un échantillon de taille 80, pour la loi de densité \(\mathcal{N}(-0.3,0.5)\)
  - un échantillon de taille 100, pour la loi de densité \(\mathcal{N}(0.7,0.7)\)
  - un échantillon de taille 125, pour la loi de densité \(\mathcal{N}(-0.1,0.2)\)
  - un échantillon de taille 100, pour la loi de densité \(\mathcal{N}(1.6,0.1)\)
  - un échantillon de taille 120, pour la loi de densité \(\mathcal{N}(-2,0.4)\)

[1] 1

Output du modèle:
-nbCluster = 4
-model name = Gaussian_pk_Lk_C
-criterion = BIC(1743.6639)
-likelihood = -836.2511
Cluster 1
-proportion = 0.1872
-means = -2.0836
-variances = 0.1568
Cluster 2
-proportion = 0.2799
-means = -0.1348
-variances = 0.0457
Cluster 3
-proportion = 0.3942
-means = 0.5322
-variances = 0.5031
Cluster 4
-proportion = 0.1387
-means = 1.6112
-variances = 0.0065

Dans ce cas là, on voit que l’algorithme nous trouve un modèle avec 4 gaussiennes.
C’est surement dû au faite que notre mélange comporte des distributions relativement proches, et que le critère BIC donne une pénalité sur le nombre de clusters.

Regardons maintenant le même algorithme en changeant les valeurs possibles de Nb_Cluster

[1] 1

Dans ce cas là le modèle donne:
-nbCluster = 7
-criterion = BIC(1782.7787)
-likelihood = -826.6969
Cluster 1
-proportion = 0.0917
-means = -1.7996
-variances = 0.0756
Cluster 2
-proportion = 0.1703
-means = 0.9642
-variances = 0.1026
Cluster 3
-proportion = 0.0949
- means = -2.3672
-variances = 0.0659
Cluster 4
-proportion = 0.1355
-means = -0.2156
-variances = 0.1596
Cluster 5
-proportion = 0.2848
-means = -0.0801
-variances = 0.0467
Cluster 6
-proportion = 0.0758
-means = 0.6028
-variances = 0.7774
Cluster 7
-proportion = 0.1471
-means = 1.6133
-variances = 0.0068

Le critère BIC et la likelihood ont forcément légèrement augmenté.
Le modèle n’est pas très performant mais il faut noté que les echantillons pour chaque gaussiennes etaient bas et celles-ci se chevauchent beaucoup. Cela montre en quelque sorte les limites de l’algorithme EM, qui est evidemment très dépendant du nombre de données et des densités en entrée.

Choix du nombre de classes

1. Mélange unidimensionnel

Nous simulons des échantillons de taille 200, dans des proportions égales.
Nos deux gaussiennes sont: \(\mathcal{N}(10,1)\) et \(\mathcal{N}(2,1)\).

\(\bullet~\underline{Critère~BIC}\)

[1] 1

\(\bullet~\underline{Critère~ICL}\)

[1] 1

Ce critère nous donne une estimation de la qualité du GMM en termes de prédiction des données dont on dispose. Plus le BIC est faible, plus le modèle est capable de prédire les données qu’on a. Afin d’éviter l’overfitting, cette technique pénalise les modèles avec un grand nombre de clusters.

2. Données multidimensionnelles “seeds”

Nous allons utiliser Rmixmod pour sélectionner le nombre de classes et la forme du modèle par BIC d’une part, par ICL d’autre part.

\(\bullet~\underline{Critère~BIC}\)

[1] 1 2 3 4 5 6 7

\(\bullet~\underline{Critère~ICL}\)

[1] 1 2 3 4 5 6 7

3. Sélection unique du nombre de classes

\(\bullet~\underline{Forme~pkLI}\)

[1] 1 2 3 4 5 6 7

\(\bullet~\underline{Forme~pkLkC}\)

[1] 1 2 3 4 5 6 7

Avec le modele pk_LI, le BIC est de 2682.7555 est trop eleve que le modele Gaussian_pk_Lk_C avec son BIC est de -1644.7596.
Cela nous signifie que le modele pk_Lk_C est mieux que le modele pk_LI

Comparaison des modeles de Gaussiennes

Par ce graphique, nous comparons tous les modèles de mélange gaussien disponible dans Rmixmod par rapport au critère BIC et le nombre de clusters k. Bien que les graphiques semblent à première vue difficiles à lire, ils révèlent trois groupes différents de modèles, qui sont : (“General”,“Spherical”,“Diagonal”). Par exemple, pour “Pk_LI”, le BIC est de 2682 et ce modèle est de la famille “Spherical”, le BIC de pK_LKC est de -1644 et il appartient à la famille “General”.

Le but de cette partie est de sélectionner le meilleur modèle ayant le critère BIC le plus petit pour notre algorithme EM pour bien choisir un bon nombre de composantes. Les deux modèles étudiées sont : “pk_Lk_C” et “pk_LI”.

Finalement, nous voyons qu’avec le modele pk_LI, le bic est de 2682, ce qui est bien plus élevé qu’avec le modèle pk_Lk_C qui lui a un score BIC de -1644. Cela signifie donc que selon le critère BIC, le modèle sélectionné est pk_Lk_C.

Question bonus

Pour illustrer un exemple où Rmixmod pense qu’on a affaire a un mélange gaussien bien qu’il s’agit en réalité d’une seulle gaussienne, il suffit de réduire le nombre de simulation de notre gaussienne et de lui donner une une matrice de covariance avec de grandes valeurs tout en ayant une gaussienne en grande dimension. Notre exempble ce fait avec la gaussienne 5d suivante:
\(\mathcal{N}(\mu, \sigma)\) avec \(\mu = (0,0,0,0,0)\) et \(\sigma = 100*Id_5\)

****************************************
*** INPUT:
****************************************
* nbCluster =  1 2 3 4 5 6 7 8 9 10 
* criterion =  BIC 
****************************************
*** MIXMOD Models:
* list =  Gaussian_pk_Lk_C 
* This list includes only models with free proportions.
****************************************
* data (limited to a 10x10 matrix) =
      V1     V2     V3      V4     V5    
 [1,] -3.9   -18.19 6.592   4.596  16.17 
 [2,] -18.56 -2.868 17.5    1.164  13.84 
 [3,] 5.742  1.365  9.142   -18.01 -3.399
 [4,] 6.063  13.41  7.673   1.937  11.41 
 [5,] 0.1386 -11.05 -0.2516 -1.637 3.701 
 [6,] -3.808 6.53   20.61   -17.97 5.841 
 [7,] -7.228 -6.292 -18.16  -2.593 3.346 
 [8,] -14.27 19.39  -7.595  -22.79 -1.141
 [9,] 23.52  15.96  12.77   7.89   4.615 
[10,] -4.381 -15.08 -22.23  -11.79 -17.83
* ... ...
****************************************
*** MIXMOD Strategy:
* algorithm            =  EM 
* number of tries      =  1 
* number of iterations =  200 
* epsilon              =  0.001 
*** Initialization strategy:
* algorithm            =  smallEM 
* number of tries      =  10 
* number of iterations =  5 
* epsilon              =  0.001 
* seed                 =  NULL 
****************************************


****************************************
*** BEST MODEL OUTPUT:
*** According to the BIC criterion
****************************************
* nbCluster   =  1 
* model name  =  Gaussian_pk_Lk_C 
* criterion   =  BIC(3796.9187)
* likelihood  =  -1852.4076 
****************************************
*** Cluster 1 
* proportion =  1.0000 
* means      =  -1.3219 -0.1873 0.2007 -0.7230 0.3885 
* variances  = |   102.0197    -2.0238    -8.9963    11.8388    -2.8029 |
               |    -2.0238   103.1458    -8.9114     5.0213    19.2153 |
               |    -8.9963    -8.9114   117.4786    -3.6591    20.3573 |
               |    11.8388     5.0213    -3.6591    91.7639     3.3772 |
               |    -2.8029    19.2153    20.3573     3.3772    84.9032 |
****************************************

En effet, le modèle trouve un nombre de clusters égale à 9. C’est à dire, l’algorithme croit avoir un mélange de 9 gaussiennes. Sachant que Rmixmod utilse un critère “BIC”, “ICL” ou “NEC” il est difficile de lui faire penser qu’on a un mélange au lieu d’une seule gaussienne. Les deux manières d’accentuer ce phénomène sont: la réduction de la taille d’échantillon (ici n= 100) L’augmentation de la dimension des gaussiennes (ici d=5)

Sources :
- Pour les tracés en 3D: https://plotly.com/
- La reference pour rmixmod: https://cran.r-project.org/web/packages/Rmixmod/Rmixmod.pdf