--- jupytext: text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.16.0 kernelspec: name: python3 display_name: Python 3 --- # Statistiques multivariées La plupart des jeux de données réels comportent des dizaines, voire des centaines de variables mesurées simultanément sur les mêmes individus. Analyser ces variables une par une revient à ignorer l'information contenue dans leurs relations. Les **statistiques multivariées** exploitent précisément ces relations : elles décrivent la structure d'un nuage de points en grande dimension, détectent des groupes naturels, réduisent la dimensionnalité tout en préservant l'information essentielle. Ce chapitre couvre les outils fondamentaux : la distribution normale multivariée, l'analyse en composantes principales (ACP), l'analyse factorielle des correspondances (AFC), les méthodes de clustering (k-means, clustering hiérarchique, DBSCAN), et un aperçu des méthodes de réduction de dimension modernes. ```{code-cell} python3 :tags: [hide-input] import numpy as np import pandas as pd import matplotlib.pyplot as plt import matplotlib.patches as mpatches from matplotlib.patches import Ellipse import seaborn as sns from scipy import stats from scipy.cluster.hierarchy import dendrogram, linkage, fcluster from scipy.spatial.distance import pdist from sklearn.decomposition import PCA from sklearn.preprocessing import StandardScaler from sklearn.cluster import KMeans, DBSCAN from sklearn.metrics import silhouette_score, silhouette_samples from sklearn.datasets import load_iris, load_wine, make_blobs, make_moons import warnings warnings.filterwarnings('ignore') # Style global plt.rcParams.update({ 'figure.dpi': 110, 'axes.spines.top': False, 'axes.spines.right': False, 'font.size': 11, }) rng = np.random.default_rng(42) ``` ## Distribution normale multivariée ### Définition et paramètres Un vecteur aléatoire $\mathbf{X} = (X_1, \ldots, X_p)^\top$ suit une **loi normale multivariée** $\mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})$ si sa densité est : $$f(\mathbf{x}) = \frac{1}{(2\pi)^{p/2} |\boldsymbol{\Sigma}|^{1/2}} \exp\!\left(-\frac{1}{2}(\mathbf{x} - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})\right)$$ Les paramètres sont : - **$\boldsymbol{\mu} \in \mathbb{R}^p$** : vecteur des moyennes marginales - **$\boldsymbol{\Sigma} \in \mathbb{R}^{p \times p}$** : matrice de covariance, symétrique définie positive La matrice de covariance $\boldsymbol{\Sigma}$ contient sur la diagonale les variances $\text{Var}(X_i) = \sigma_i^2$, et hors-diagonale les covariances $\text{Cov}(X_i, X_j) = \rho_{ij}\sigma_i\sigma_j$ où $\rho_{ij}$ est la corrélation entre $X_i$ et $X_j$. ### Ellipses de confiance En dimension 2, les courbes de niveau de la densité sont des **ellipses**. L'ellipse à $(1 - \alpha)$% de confiance est l'ensemble des points $\mathbf{x}$ vérifiant : $$(\mathbf{x} - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) \leq \chi^2_{p,\, 1-\alpha}$$ où $\chi^2_{p, 1-\alpha}$ est le quantile $(1-\alpha)$ du chi-deux à $p$ degrés de liberté. Cette distance au centre est la **distance de Mahalanobis**, qui généralise la distance euclidienne en tenant compte des corrélations. ```{code-cell} python3 :tags: [hide-input] def ellipse_confiance(ax, mu, sigma, niveau=0.95, couleur='steelblue', label=None): """Trace une ellipse de confiance multivariée.""" chi2_val = stats.chi2.ppf(niveau, df=2) vals, vecs = np.linalg.eigh(sigma) angle = np.degrees(np.arctan2(vecs[1, 1], vecs[0, 1])) width, height = 2 * np.sqrt(chi2_val * vals[::-1]) ell = Ellipse(xy=mu, width=width, height=height, angle=angle, edgecolor=couleur, facecolor=couleur, alpha=0.15, linewidth=2) ax.add_patch(ell) ell2 = Ellipse(xy=mu, width=width, height=height, angle=angle, edgecolor=couleur, facecolor='none', linewidth=2, label=label) ax.add_patch(ell2) fig, axes = plt.subplots(1, 3, figsize=(14, 4)) configs = [ (np.array([0, 0]), np.array([[1, 0], [0, 1]]), "Indépendantes\n(ρ = 0)"), (np.array([0, 0]), np.array([[2, 1.6], [1.6, 2]]), "Corrélées\n(ρ = 0.8)"), (np.array([0, 0]), np.array([[3, -1.5], [-1.5, 1]]), "Anti-corrélées\n(ρ = -0.87)"), ] for ax, (mu, sigma, titre) in zip(axes, configs): n = 300 X = rng.multivariate_normal(mu, sigma, n) ax.scatter(X[:, 0], X[:, 1], alpha=0.35, s=15, color='steelblue') ellipse_confiance(ax, mu, sigma, 0.68, 'steelblue', '68%') ellipse_confiance(ax, mu, sigma, 0.95, 'orangered', '95%') ellipse_confiance(ax, mu, sigma, 0.99, 'seagreen', '99%') ax.set_title(titre, fontsize=11) ax.set_xlabel('$X_1$'); ax.set_ylabel('$X_2$') ax.set_aspect('equal') axes[0].legend(fontsize=9, loc='upper right') fig.suptitle('Ellipses de confiance — Loi normale bivariée', fontsize=13, y=1.01) plt.tight_layout() plt.show() ``` La forme des ellipses révèle immédiatement la structure de dépendance : des ellipses allongées selon une diagonale indiquent une corrélation forte entre les variables. --- ## Analyse en Composantes Principales (ACP) ### Principe et objectif L'ACP est une **transformation linéaire orthogonale** qui réexprime les données dans un nouveau système d'axes — les **composantes principales** — ordonnés par variance décroissante. La première composante principale (CP1) est la direction qui maximise la variance des projections ; la deuxième (CP2) est orthogonale à CP1 et maximise la variance résiduelle, etc. **Quand l'utiliser ?** - Visualiser des données en grande dimension (projection sur CP1, CP2) - Réduire la dimensionnalité avant un modèle (débruitage, décorrélation) - Détecter la structure latente dans les données - Diagnostiquer des problèmes (multicolinéarité) ### Prétraitement indispensable Avant toute ACP, il faut **centrer et réduire** les variables (soustraction de la moyenne, division par l'écart-type). Sans réduction, les variables à grande variance domineront artificiellement les premières composantes. La réduction est impérative dès que les variables sont exprimées dans des unités différentes. ### Algorithme (rappel) 1. Centrer-réduire la matrice $\mathbf{X}$ ($n \times p$) → $\mathbf{Z}$ 2. Calculer la matrice de corrélation $\mathbf{R} = \frac{1}{n-1}\mathbf{Z}^\top\mathbf{Z}$ 3. Décomposer en valeurs propres : $\mathbf{R} = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^\top$ 4. Les **loadings** sont les colonnes de $\mathbf{V}$ (vecteurs propres) 5. Les **scores** sont $\mathbf{F} = \mathbf{Z}\mathbf{V}$ (coordonnées dans le nouveau repère) ### Choix du nombre de composantes Deux critères courants : - **Scree plot** : chercher le "coude" dans la courbe des valeurs propres - **Règle de Kaiser** : conserver les composantes dont la valeur propre $> 1$ (c'est-à-dire celles qui expliquent plus de variance qu'une variable originale) - **Variance cumulée** : viser 70-80% de variance expliquée cumulée ```{code-cell} python3 :tags: [hide-input] # Chargement du jeu de données Vin (13 variables chimiques, 3 cépages) wine = load_wine() X_wine = pd.DataFrame(wine.data, columns=wine.feature_names) y_wine = wine.target labels_wine = wine.target_names # ACP avec sklearn scaler = StandardScaler() X_scaled = scaler.fit_transform(X_wine) pca = PCA() scores = pca.fit_transform(X_scaled) variance_expliquee = pca.explained_variance_ratio_ variance_cumulee = np.cumsum(variance_expliquee) fig, axes = plt.subplots(1, 3, figsize=(15, 4)) # --- Scree plot --- ax = axes[0] ax.bar(range(1, len(variance_expliquee)+1), variance_expliquee * 100, color='steelblue', alpha=0.7, label='Individuelle') ax.plot(range(1, len(variance_cumulee)+1), variance_cumulee * 100, 'o-', color='orangered', label='Cumulée') ax.axhline(80, color='gray', linestyle='--', linewidth=1, label='80%') ax.axvline(np.argmax(variance_cumulee >= 0.80) + 1, color='gray', linestyle=':') # Règle de Kaiser : val propre > 1 ↔ variance > 100/p % kaiser_seuil = 100 / X_wine.shape[1] ax.axhline(kaiser_seuil, color='seagreen', linestyle='--', linewidth=1, label=f'Kaiser ({kaiser_seuil:.1f}%)') ax.set_xlabel('Composante principale') ax.set_ylabel('Variance expliquée (%)') ax.set_title('Scree plot') ax.legend(fontsize=9) # --- Projection CP1-CP2 --- ax = axes[1] couleurs = ['#e74c3c', '#3498db', '#2ecc71'] for i, (nom, col) in enumerate(zip(labels_wine, couleurs)): mask = y_wine == i ax.scatter(scores[mask, 0], scores[mask, 1], c=col, label=nom, alpha=0.75, s=40) ax.set_xlabel(f'CP1 ({variance_expliquee[0]*100:.1f}%)') ax.set_ylabel(f'CP2 ({variance_expliquee[1]*100:.1f}%)') ax.set_title('Projection CP1–CP2') ax.legend(fontsize=9) # --- Biplot : scores + loadings --- ax = axes[2] scale = np.sqrt(pca.explained_variance_[0]) * 2 for i, (nom, col) in enumerate(zip(labels_wine, couleurs)): mask = y_wine == i ax.scatter(scores[mask, 0] / scale, scores[mask, 1] / scale, c=col, alpha=0.4, s=25) loadings = pca.components_.T # shape (p, n_components) for j, feat in enumerate(wine.feature_names): ax.annotate('', xy=(loadings[j, 0], loadings[j, 1]), xytext=(0, 0), arrowprops=dict(arrowstyle='->', color='orangered', lw=1.5)) ax.text(loadings[j, 0]*1.1, loadings[j, 1]*1.1, feat[:8], fontsize=7, ha='center', va='center', color='orangered') ax.axhline(0, color='gray', linewidth=0.5) ax.axvline(0, color='gray', linewidth=0.5) ax.set_xlabel(f'CP1 ({variance_expliquee[0]*100:.1f}%)') ax.set_ylabel(f'CP2 ({variance_expliquee[1]*100:.1f}%)') ax.set_title('Biplot (scores + loadings)') plt.suptitle('ACP sur le jeu de données Vin (13 variables, 178 échantillons)', fontsize=13, y=1.01) plt.tight_layout() plt.show() print(f"Variance expliquée par CP1 : {variance_expliquee[0]*100:.1f}%") print(f"Variance expliquée par CP1+CP2 : {variance_cumulee[1]*100:.1f}%") print(f"Nombre de composantes pour 80% de variance : " f"{np.argmax(variance_cumulee >= 0.80) + 1}") ``` ### Interprétation du biplot Le **biplot** superpose deux informations : - Les **scores** (points) : position de chaque individu dans le nouveau repère - Les **loadings** (flèches) : contribution de chaque variable originale aux composantes Des flèches longues indiquent des variables bien représentées sur le plan. Des flèches proches (angle faible) correspondent à des variables fortement corrélées entre elles. Des flèches opposées (angle ≈ 180°) correspondent à une anti-corrélation. ### Analyse des loadings ```{code-cell} python3 :tags: [hide-input] pca3 = PCA(n_components=3) pca3.fit(X_scaled) loadings_df = pd.DataFrame( pca3.components_.T, index=wine.feature_names, columns=['CP1', 'CP2', 'CP3'] ) fig, ax = plt.subplots(figsize=(10, 5)) loadings_df.sort_values('CP1').plot(kind='barh', ax=ax, width=0.7, color=['steelblue', 'orangered', 'seagreen']) ax.axvline(0, color='black', linewidth=0.8) ax.set_title('Loadings des 3 premières composantes principales — Données Vin', fontsize=12) ax.set_xlabel('Loading (corrélation variable–composante)') ax.legend(title='Composante', fontsize=9) plt.tight_layout() plt.show() ``` ```{admonition} Interprétation des loadings Les loadings peuvent être interprétés comme des corrélations entre les variables originales et les composantes. Un loading de 0.8 signifie que la variable est corrélée à 80% avec cette composante. En pratique, on nomme une composante d'après les variables qui ont les loadings les plus élevés (en valeur absolue). ``` --- ## Analyse Factorielle des Correspondances (AFC) L'AFC est l'analogue de l'ACP pour les **tableaux de contingence** (variables catégorielles). Elle décompose le $\chi^2$ d'un tableau de contingence en dimensions successives et permet de visualiser les relations entre modalités de lignes et de colonnes. ### Principe Pour un tableau de contingence $\mathbf{N}$ ($r \times c$) : 1. Calculer le tableau des **profils** (fréquences relatives par ligne et colonne) 2. Calculer la statistique du $\chi^2$ et l'**inertie totale** $= \chi^2 / n$ 3. Décomposition en valeurs singulières du tableau des résidus standardisés 4. Les **axes factoriels** maximisent l'inertie expliquée ```{code-cell} python3 :tags: [hide-input] # Exemple : AFC sur un tableau fumeurs × catégorie socio-professionnelle np.random.seed(42) tableau = pd.DataFrame({ 'Non-fumeur': [135, 85, 210, 60], 'Fumeur occ.': [ 45, 30, 70, 25], 'Fumeur rég.': [ 20, 55, 35, 40], 'Gros fumeur': [ 5, 25, 10, 30], }, index=['Cadre', 'Employé', 'Ouvrier', 'Étudiant']) # AFC manuelle (SVD) N = tableau.values.astype(float) n = N.sum() P = N / n # tableau de fréquences r = P.sum(axis=1) # masses lignes c = P.sum(axis=0) # masses colonnes Dr = np.diag(r) Dc = np.diag(c) # Matrice des résidus standardisés S = np.diag(1/np.sqrt(r)) @ (P - np.outer(r, c)) @ np.diag(1/np.sqrt(c)) U, d, Vt = np.linalg.svd(S) inertie_totale = (d**2).sum() inertie_axe = d**2 print(f"Inertie totale (= χ²/n) : {inertie_totale:.4f}") print(f"χ² associé : {inertie_totale * n:.2f}") print(f"Inertie axe 1 : {inertie_axe[0]/inertie_totale*100:.1f}%") print(f"Inertie axes 1+2 : {inertie_axe[:2].sum()/inertie_totale*100:.1f}%") # Coordonnées lignes et colonnes F_r = np.diag(1/np.sqrt(r)) @ U[:, :2] @ np.diag(d[:2]) # lignes F_c = np.diag(1/np.sqrt(c)) @ Vt[:2, :].T @ np.diag(d[:2]) # colonnes fig, ax = plt.subplots(figsize=(7, 5)) for i, nom in enumerate(tableau.index): ax.scatter(F_r[i, 0], F_r[i, 1], s=120, color='steelblue', zorder=5) ax.annotate(nom, (F_r[i, 0], F_r[i, 1]), textcoords='offset points', xytext=(6, 4), fontsize=11, color='steelblue') for j, nom in enumerate(tableau.columns): ax.scatter(F_c[j, 0], F_c[j, 1], s=120, marker='D', color='orangered', zorder=5) ax.annotate(nom, (F_c[j, 0], F_c[j, 1]), textcoords='offset points', xytext=(6, -10), fontsize=10, color='orangered') ax.axhline(0, color='gray', linewidth=0.5) ax.axvline(0, color='gray', linewidth=0.5) ax.set_xlabel(f'Axe 1 ({inertie_axe[0]/inertie_totale*100:.1f}% d\'inertie)') ax.set_ylabel(f'Axe 2 ({inertie_axe[1]/inertie_totale*100:.1f}% d\'inertie)') ax.set_title('AFC — Catégorie professionnelle × habitudes tabagiques', fontsize=12) legend_elems = [ mpatches.Patch(color='steelblue', label='Catégories (lignes)'), mpatches.Patch(color='orangered', label='Statut tabagique (colonnes)'), ] ax.legend(handles=legend_elems, fontsize=9) plt.tight_layout() plt.show() ``` ```{note} En AFC, les modalités proches sur le graphique sont associées : les cadres se retrouvent du côté "non-fumeur" tandis que les employés se rapprochent des "gros fumeurs". Les modalités proches de l'origine sont banales (elles ne dépassent pas la fréquence attendue sous indépendance). ``` --- ## Clustering non supervisé Le clustering consiste à **regrouper des individus similaires** sans connaissance préalable des groupes. Contrairement à la classification supervisée, il n'y a pas d'étiquettes d'entraînement. ### K-means #### Algorithme L'algorithme de Lloyd (communément appelé k-means) procède par itérations : 1. Initialiser $k$ centroïdes aléatoirement (ou par k-means++) 2. **Affecter** chaque point au centroïde le plus proche (distance euclidienne) 3. **Mettre à jour** chaque centroïde : nouvelle position = moyenne du cluster 4. Répéter jusqu'à convergence (les affectations ne changent plus) L'algorithme minimise l'**inertie intra-cluster** (aussi appelée WCSS, *Within-Cluster Sum of Squares*) : $$\text{WCSS} = \sum_{j=1}^{k} \sum_{\mathbf{x}_i \in C_j} \|\mathbf{x}_i - \boldsymbol{\mu}_j\|^2$$ #### Méthode du coude et silhouette ```{code-cell} python3 :tags: [hide-input] # Données simulées : 4 clusters bien séparés X_blobs, y_true = make_blobs(n_samples=400, centers=4, cluster_std=0.9, random_state=42) # Méthode du coude k_range = range(2, 11) wcss = [] sil_scores = [] for k in k_range: km = KMeans(n_clusters=k, random_state=42, n_init=10) labels = km.fit_predict(X_blobs) wcss.append(km.inertia_) sil_scores.append(silhouette_score(X_blobs, labels)) fig, axes = plt.subplots(1, 3, figsize=(15, 4)) # Coude ax = axes[0] ax.plot(list(k_range), wcss, 'o-', color='steelblue', linewidth=2) ax.axvline(4, color='orangered', linestyle='--', label='k optimal = 4') ax.set_xlabel('Nombre de clusters k') ax.set_ylabel('Inertie intra-cluster (WCSS)') ax.set_title('Méthode du coude') ax.legend() # Silhouette ax = axes[1] ax.plot(list(k_range), sil_scores, 'o-', color='seagreen', linewidth=2) ax.axvline(4, color='orangered', linestyle='--', label='k optimal = 4') ax.set_xlabel('Nombre de clusters k') ax.set_ylabel('Score silhouette moyen') ax.set_title('Score silhouette') ax.legend() # Clustering final avec k=4 km_final = KMeans(n_clusters=4, random_state=42, n_init=10) labels_km = km_final.fit_predict(X_blobs) ax = axes[2] palette = ['#e74c3c', '#3498db', '#2ecc71', '#f39c12'] for j in range(4): mask = labels_km == j ax.scatter(X_blobs[mask, 0], X_blobs[mask, 1], c=palette[j], alpha=0.65, s=30, label=f'Cluster {j+1}') ax.scatter(km_final.cluster_centers_[:, 0], km_final.cluster_centers_[:, 1], c='black', marker='X', s=150, zorder=10, label='Centroïdes') ax.set_title('K-means (k=4)') ax.legend(fontsize=9) plt.suptitle('Sélection du k optimal et résultat K-means', fontsize=13, y=1.01) plt.tight_layout() plt.show() print(f"Score silhouette pour k=4 : {silhouette_score(X_blobs, labels_km):.3f}") ``` #### Silhouette plot détaillé ```{code-cell} python3 :tags: [hide-input] fig, ax = plt.subplots(figsize=(8, 5)) km4 = KMeans(n_clusters=4, random_state=42, n_init=10) labels4 = km4.fit_predict(X_blobs) sil_vals = silhouette_samples(X_blobs, labels4) y_lower = 10 couleurs_sil = ['#e74c3c', '#3498db', '#2ecc71', '#f39c12'] for j in range(4): sil_j = np.sort(sil_vals[labels4 == j]) taille = len(sil_j) y_upper = y_lower + taille ax.fill_betweenx(np.arange(y_lower, y_upper), 0, sil_j, alpha=0.7, color=couleurs_sil[j], label=f'Cluster {j+1}') ax.text(-0.05, y_lower + taille / 2, str(j+1), fontsize=9) y_lower = y_upper + 10 ax.axvline(silhouette_score(X_blobs, labels4), color='black', linestyle='--', linewidth=1.5, label=f'Moyenne = {silhouette_score(X_blobs, labels4):.3f}') ax.set_xlabel('Coefficient silhouette') ax.set_ylabel('Individus (par cluster)') ax.set_title('Silhouette plot — K-means k=4') ax.legend(fontsize=9) plt.tight_layout() plt.show() ``` ```{admonition} Interpréter le score silhouette Le coefficient silhouette $s_i$ d'un point $i$ compare sa distance moyenne aux points de son cluster ($a_i$) à la distance moyenne au cluster le plus proche ($b_i$) : $$s_i = \frac{b_i - a_i}{\max(a_i, b_i)} \in [-1, 1]$$ - $s_i \approx 1$ : le point est bien placé dans son cluster - $s_i \approx 0$ : le point est à la frontière entre deux clusters - $s_i < 0$ : le point est probablement mal classé ``` ### Clustering hiérarchique agglomératif (CAH) La CAH construit un **dendrogramme** en fusionnant itérativement les clusters les plus proches. Contrairement au k-means, elle n'impose pas de nombre de clusters *a priori* et produit une hiérarchie complète. #### Méthodes de liaison La définition de la distance entre clusters varie selon la méthode : - **Ward** : minimise l'augmentation de l'inertie totale lors de la fusion (recommandée pour des clusters compacts et de taille similaire) - **Complete** (lien maximum) : distance entre les deux points les plus éloignés - **Average** (UPGMA) : distance moyenne entre tous les couples de points - **Single** (lien minimum) : distance entre les deux points les plus proches (sensible aux chaînes) ```{code-cell} python3 :tags: [hide-input] fig, axes = plt.subplots(1, 4, figsize=(16, 5)) methodes = ['ward', 'complete', 'average', 'single'] titres = ['Ward', 'Complete (max)', 'Average', 'Single (min)'] for ax, meth, titre in zip(axes, methodes, titres): Z = linkage(X_blobs[:80], method=meth) # 80 points pour lisibilité dendrogram(Z, ax=ax, no_labels=True, color_threshold=0.7*max(Z[:, 2])) ax.set_title(f'Liaison : {titre}', fontsize=10) ax.set_xlabel('Individus') ax.set_ylabel('Distance de fusion') plt.suptitle('Dendrogrammes — Comparaison des méthodes de liaison (CAH)', fontsize=13, y=1.01) plt.tight_layout() plt.show() ``` ```{code-cell} python3 :tags: [hide-input] # CAH complète et partition en 4 clusters Z_ward = linkage(X_blobs, method='ward') labels_cah = fcluster(Z_ward, t=4, criterion='maxclust') - 1 fig, axes = plt.subplots(1, 2, figsize=(12, 5)) # Dendrogramme tronqué ax = axes[0] dendrogram(Z_ward, ax=ax, no_labels=True, truncate_mode='lastp', p=20, color_threshold=Z_ward[-3, 2]) ax.axhline(Z_ward[-3, 2], color='orangered', linestyle='--', linewidth=1.5, label=f'Coupe → 4 clusters') ax.set_title('Dendrogramme Ward (tronqué, 20 dernières fusions)') ax.set_ylabel('Distance de fusion') ax.legend() # Résultat ax = axes[1] for j in range(4): mask = labels_cah == j ax.scatter(X_blobs[mask, 0], X_blobs[mask, 1], c=palette[j], alpha=0.65, s=30, label=f'Cluster {j+1}') ax.set_title('Résultat CAH Ward (4 clusters)') ax.legend(fontsize=9) plt.suptitle('Clustering hiérarchique agglomératif', fontsize=13, y=1.01) plt.tight_layout() plt.show() print(f"Score silhouette CAH Ward : {silhouette_score(X_blobs, labels_cah):.3f}") print(f"Score silhouette K-means : {silhouette_score(X_blobs, labels_km):.3f}") ``` ### DBSCAN DBSCAN (*Density-Based Spatial Clustering of Applications with Noise*) est un algorithme basé sur la **densité** plutôt que sur la distance aux centroïdes. Il est particulièrement adapté aux clusters de forme arbitraire et identifie automatiquement les **points de bruit**. **Paramètres clés :** - **$\epsilon$ (eps)** : rayon de voisinage — distance maximale pour qu'un point soit considéré comme voisin - **minPts (min_samples)** : nombre minimal de voisins pour qu'un point soit un **point cœur** **Classification des points :** - **Point cœur** : au moins `minPts` voisins dans le rayon $\epsilon$ - **Point frontière** : dans le voisinage d'un point cœur, mais pas assez de voisins lui-même - **Point de bruit** : ni cœur ni frontière (label = −1 dans sklearn) ```{code-cell} python3 :tags: [hide-input] # DBSCAN est particulièrement adapté aux formes non convexes X_moons, y_moons = make_moons(n_samples=400, noise=0.07, random_state=42) X_anneaux = np.vstack([ rng.standard_normal((200, 2)) * 0.5, rng.standard_normal((200, 2)) * 1.5 + 2.5, ]) fig, axes = plt.subplots(2, 3, figsize=(14, 8)) datasets = [ (X_moons, 'Deux lunes'), (X_blobs, 'Blobs gaussiens'), (X_anneaux, 'Distributions imbriquées'), ] params_dbscan = [ {'eps': 0.15, 'min_samples': 5}, {'eps': 0.8, 'min_samples': 5}, {'eps': 0.35, 'min_samples': 5}, ] for col, (X_d, titre), params in zip(range(3), datasets, params_dbscan): # K-means km_comp = KMeans(n_clusters=2, random_state=42, n_init=10).fit(X_d) lk = km_comp.labels_ axes[0, col].scatter(X_d[:, 0], X_d[:, 1], c=lk, cmap='Set1', alpha=0.6, s=20) axes[0, col].set_title(f'{titre}\nK-means (k=2)', fontsize=10) # DBSCAN db = DBSCAN(**params).fit(X_d) ld = db.labels_ n_bruit = np.sum(ld == -1) couleurs_db = np.where(ld == -1, 'gray', np.where(ld == 0, '#e74c3c', '#3498db')) axes[1, col].scatter(X_d[:, 0], X_d[:, 1], c=couleurs_db, alpha=0.6, s=20) axes[1, col].set_title(f'{titre}\nDBSCAN (ε={params["eps"]}, {n_bruit} bruits)', fontsize=10) axes[0, 0].set_ylabel('K-means', fontsize=11) axes[1, 0].set_ylabel('DBSCAN', fontsize=11) plt.suptitle('Comparaison K-means vs DBSCAN sur différentes géométries', fontsize=13, y=1.01) plt.tight_layout() plt.show() ``` ```{admonition} Limites du K-means Le K-means suppose des clusters convexes et de variance homogène. Il est inadapté aux formes en lune, en anneau ou allongées. DBSCAN gère ces cas mais requiert un bon réglage de ε et minPts. Pour choisir ε, on trace la courbe des distances au k-ième plus proche voisin (k = minPts) et on cherche le "coude". ``` --- ## Réduction de dimension : t-SNE et UMAP ### Comparaison avec l'ACP | Méthode | Type | Global/Local | Déterministe | Utilisable pour modélisation | |---------|------|-------------|--------------|------------------------------| | ACP | Linéaire | Global | Oui | Oui | | t-SNE | Non-linéaire | Local | Non (stochastique) | Non (embedding seulement) | | UMAP | Non-linéaire | Global + Local | Presque | Avec précautions | **t-SNE** (*t-distributed Stochastic Neighbor Embedding*) minimise la divergence KL entre distributions de similarités dans l'espace d'origine et dans l'espace 2D. Il préserve les **structures locales** mais distord les distances globales. Les distances entre clusters sur un t-SNE ne sont **pas** interprétables directement. **UMAP** (*Uniform Manifold Approximation and Projection*) repose sur une théorie différentielle des variétés riemanniennes. Il est plus rapide que t-SNE, préserve mieux la structure globale, et peut être utilisé pour projeter de nouveaux points. ```{code-cell} python3 :tags: [hide-input] from sklearn.manifold import TSNE # ACP vs t-SNE sur les données Vin pca2 = PCA(n_components=2) X_pca = pca2.fit_transform(X_scaled) tsne = TSNE(n_components=2, perplexity=30, random_state=42, max_iter=1000) X_tsne = tsne.fit_transform(X_scaled) fig, axes = plt.subplots(1, 2, figsize=(12, 4)) couleurs_vin = ['#e74c3c', '#3498db', '#2ecc71'] for i, (nom, col) in enumerate(zip(labels_wine, couleurs_vin)): mask = y_wine == i axes[0].scatter(X_pca[mask, 0], X_pca[mask, 1], c=col, label=nom, alpha=0.75, s=40) axes[1].scatter(X_tsne[mask, 0], X_tsne[mask, 1], c=col, label=nom, alpha=0.75, s=40) axes[0].set_title(f'ACP — CP1 ({variance_expliquee[0]*100:.1f}%) + CP2 ({variance_expliquee[1]*100:.1f}%)') axes[0].set_xlabel('CP1'); axes[0].set_ylabel('CP2') axes[0].legend() axes[1].set_title('t-SNE (perplexité=30)') axes[1].set_xlabel('Dim 1'); axes[1].set_ylabel('Dim 2') axes[1].legend() plt.suptitle('ACP vs t-SNE — Données Vin', fontsize=13, y=1.01) plt.tight_layout() plt.show() ``` --- ## MANOVA La **MANOVA** (*Multivariate Analysis of Variance*) est l'extension multivariée de l'ANOVA. Là où l'ANOVA teste si la moyenne d'**une** variable dépend d'un facteur de groupe, la MANOVA teste si le **vecteur de moyennes** de plusieurs variables dépend simultanément d'un facteur. **Hypothèses :** - $H_0$ : $\boldsymbol{\mu}_1 = \boldsymbol{\mu}_2 = \ldots = \boldsymbol{\mu}_k$ (tous les vecteurs moyens sont égaux) - $H_1$ : au moins deux vecteurs moyens diffèrent **Statistiques de test** (plusieurs variantes) : Lambda de Wilks $\Lambda$, Trace de Pillai, Trace d'Hotelling, Valeur propre maximale de Roy. ```{code-cell} python3 from statsmodels.multivariate.manova import MANOVA # MANOVA sur les données Iris (3 groupes, 4 variables) iris = load_iris() df_iris = pd.DataFrame(iris.data, columns=['sepal_length', 'sepal_width', 'petal_length', 'petal_width']) df_iris['species'] = iris.target_names[iris.target] manova = MANOVA.from_formula( 'sepal_length + sepal_width + petal_length + petal_width ~ species', data=df_iris ) result = manova.mv_test() print(result) ``` ```{note} Le Lambda de Wilks proche de 0 indique que les groupes sont bien séparés dans l'espace multivarié. Ici, $p < 0.001$ confirme que les trois espèces d'iris diffèrent significativement sur l'ensemble des quatre mesures morphologiques. ``` --- ## Comparaison complète des méthodes de clustering ```{code-cell} python3 :tags: [hide-input] from sklearn.mixture import GaussianMixture # Comparaison sur les données Vin (après ACP, 2 premières CP) X_wine_2d = PCA(n_components=2).fit_transform(X_scaled) methodes_clustering = { 'K-means (k=3)' : KMeans(n_clusters=3, random_state=42, n_init=10).fit_predict(X_wine_2d), 'DBSCAN' : DBSCAN(eps=0.5, min_samples=5).fit_predict(X_wine_2d), 'GMM (k=3)' : GaussianMixture(n_components=3, random_state=42).fit_predict(X_wine_2d), } Z_wine = linkage(X_wine_2d, method='ward') labels_cah_wine = fcluster(Z_wine, t=3, criterion='maxclust') - 1 methodes_clustering['CAH Ward (k=3)'] = labels_cah_wine fig, axes = plt.subplots(2, 2, figsize=(12, 9)) axes_flat = axes.flatten() for ax, (nom_meth, labels_meth) in zip(axes_flat, methodes_clustering.items()): n_clusters = len(set(labels_meth)) - (1 if -1 in labels_meth else 0) cmap_vals = plt.cm.Set1(np.linspace(0, 0.8, n_clusters + 1)) label_map = {l: i for i, l in enumerate(sorted(set(labels_meth)))} for l in sorted(set(labels_meth)): mask = labels_meth == l nom_cluster = 'Bruit' if l == -1 else f'Cluster {l+1}' col = 'gray' if l == -1 else cmap_vals[label_map[l]] ax.scatter(X_wine_2d[mask, 0], X_wine_2d[mask, 1], c=[col], alpha=0.65, s=35, label=nom_cluster) # Vraies étiquettes en contour for i, (nom, col) in enumerate(zip(labels_wine, couleurs_vin)): mask_true = y_wine == i ax.scatter(X_wine_2d[mask_true, 0], X_wine_2d[mask_true, 1], edgecolors=col, facecolors='none', s=60, linewidth=1, alpha=0.4) if -1 not in labels_meth: sil = silhouette_score(X_wine_2d, labels_meth) ax.set_title(f'{nom_meth}\n(silhouette = {sil:.3f})', fontsize=10) else: ax.set_title(f'{nom_meth}\n({np.sum(labels_meth == -1)} bruits)', fontsize=10) ax.set_xlabel('CP1'); ax.set_ylabel('CP2') ax.legend(fontsize=8) plt.suptitle('Comparaison des méthodes de clustering — Données Vin (projection ACP)', fontsize=13, y=1.01) plt.tight_layout() plt.show() ``` --- ## Résumé | Méthode | Type | Hypothèses | Points forts | Limites | |---------|------|-----------|-------------|---------| | ACP | Réduction dim. | Linéarité, variance maximale | Interprétable, rapide | Liaisons non-linéaires perdues | | AFC | Réduction dim. | Tableau de contingence | Variables catégorielles | Données quantitatives non adaptées | | K-means | Clustering | Clusters convexes, sphériques | Rapide, scalable | Sensible aux outliers, k fixé | | CAH | Clustering | Aucune sur la forme | Dendrogramme, k libre | Lent sur grands jeux, irréversible | | DBSCAN | Clustering | Densité uniforme | Formes arbitraires, bruit | Paramètres sensibles | | t-SNE | Visualisation | — | Structures locales | Non interprétable quantitativement | ```{admonition} Bonnes pratiques en analyse multivariée 1. **Toujours centrer-réduire** avant l'ACP si les variables ont des unités différentes. 2. **Vérifier la stabilité** du clustering : relancer k-means plusieurs fois, comparer avec CAH. 3. **Ne pas sur-interpréter les distances** sur un t-SNE : seuls les regroupements locaux sont fiables. 4. **Valider les clusters** sur des données externes quand c'est possible. 5. **L'ACP est exploratoire** : elle ne teste pas d'hypothèse. Pour des conclusions inférentielles, utiliser MANOVA ou des tests sur les scores. ```