Travaux pratiques - DBSCAN

(correspond à 1 séance de TP)

Références externes utiles :

L’objectif de cette séance de TP est d’illustrer la mise en pratique de la classification automatique avec DBSCAN. Dans un premier temps, nous observerons le comportement de l’algorithme dans un cas synthétique mettant en échec k-means. Nous verrons comment utiliser les heuristiques classiques pour déterminer les paramètres de l’algorithme et comment interpréter le partitionnement obtenu. Enfin, nous appliquerons DBSCAN sur un jeu de données réelles et nous illustrerons son usage pour la détection de données aberrantes.

La classification par densité dans scikit-learn

Comme nous l’avons vu dans le TP précédent, scikit-learn propose de nombreuses méthodes de classification automatique. Nous nous intéressons pour cette séance à DBSCAN (Density-Based Spatial Clustering of Applications with Noise).

Comme k-means, DBSCAN permet d’obtenir pour une matrice d’observation n_samples x n_features un vecteur correspondant aux n_samples numéros des groupes (cluster) pour ces observations. En revanche, DBSCAN est une méthode transductive : l’algorithme ne permet pas de produire des classifications pour de nouvelles données (il est nécessaire de réentraîner le modèle en incluant les nouvelles observations).

La description de l’implémentation de DBSCAN se trouve dans https://scikit-learn.org/stable/modules/generated/sklearn.cluster.DBSCAN.html

Classification avec DBSCAN de données générées

Commençons par générer un jeu de données de 500 observations dans un espace 2D. Nous allons utiliser une fonction intégrée à Scikit-learn permettant de produire des nuages de points circulaires (make_circles).

import matplotlib.pyplot as plt
import numpy as np
import sklearn.datasets
from sklearn.utils import shuffle

# Générons un nuage de points composé de deux cercles
# Le nuage contient 500 observations (`n_samples`) bruitées par
# un bruit gaussien d'écart-type 0.1 (`noise`).
# Le rapport entre le rayon du petit cercle et du grand cercle
# est de 0.3 (`factor`).
data, labels = sklearn.datasets.make_circles(n_samples=500, noise=0.1, factor=0.3, random_state=0)

print(data.shape)
# Permutation aléatatoire des lignes de la matrice (on mélange les observations)
data, labels = shuffle(data, labels)

# Affichage du nuage de points
plt.scatter(data[:,0], data[:,1], c=labels)
plt.show()

Question :

Combien de groupes comporte ce jeu de données ?

Correction

Ce jeu de données semble comporter deux groupes (un pour chaque anneau).

Question :

Réaliser un partitionnement de ce jeu de données à l’aide de \(k\)-means. À quoi peut-on s’attendre ? Que constatez-vous ?

Correction :

from sklearn.cluster import KMeans
km = KMeans(n_clusters=2)
labels = km.fit_predict(data)

plt.scatter(data[:,0], data[:,1], c=labels) and plt.show()

Les deux cercles étant séparés par une zone sans données, une méthode basée sur la densité semble appropriée. Nous pouvons créer un modèle de clustering utilisant DBSCAN en l’important depuis Scikit-learn :

from sklearn.cluster import DBSCAN

db = DBSCAN()

Les arguments du constructeur de DBSCAN sont les suivants :

  • eps: la dimension du voisinage, c’est-à-dire la distance maximale entre deux observations permettant de les considérer comme voisines l’une de l’autre,

  • min_samples: le nombre minimal de voisins que doit avoir un point central,

  • metric: la distance à considérer (par défaut, la distance Euclidienne est utilisée).

Il est possible d’appeler les méthodes suivantes :

  • .fit(X): réalise une classification automatique en utilisant la méthode DBSCAN sur la matrice d’observations X. Les résultats sont stockés dans l’attribut .labels_.

  • .fit_predict(X): identique à .fit(X) mais renvoie directement les étiquettes des groupes.

Les attributs suivants sont disponibles après l’appel à la méthode .fit():

  • core_sample_indices_: les indices des points centraux.

  • labels_: les numéros de groupe des points de la matrice d’observations.

Question :

Quelles sont les valeurs par défaut des paramètres importants de DBSCAN dans scikit-learn (\(\varepsilon\) et \(m\)) ?

Correction

D’après la documentation de DBSCAN, eps=0.5 et min_samples=5 par défaut.

Appliquons une classification automatique par DBSCAN sur notre jeu de données. Comme pour k-means, il est possible de réaliser cette étape en deux temps en appelant fit() puis en accédant à l’attribut labels_, ou de tout faire en une seule manipulation à l’aide de la méthode fit_predict():

predictions = db.fit_predict(data)
# équivalent à
# db.fit(data)
# predictions = db.labels_

# Affichage du nuage de points coloré par les prédictions
plt.scatter(data[:,0], data[:,1], c=predictions)
plt.show()

Question :

Que constatez-vous ? Sur quel paramètre est-il vraisemblablement nécessaire de jouer pour améliorer ce résultat ?

Correction

DBSCAN ne produit qu’un seul groupe qui recouvre l’intégralité du jeu de données. Il est probable que la valeur choisie pour eps soit trop élevée (tous les points sont considérés comme étant dans le voisinage de tous les autres).

Pour raffiner notre analyse, nous allons appliquer l’heuristique de Schubert, qui exploite le graphe des \(k\)-distances dans le nuage des observations. On considère pour l’instant que min_samples est fixé à sa valeur par défaut, c’est-à-dire 5. Nous devons donc tracer les graphe des 4-distances pour notre matrice d’observations.

from sklearn.neighbors import NearestNeighbors

nn = NearestNeighbors(n_neighbors=4).fit(data)
distances, _ = nn.kneighbors(data)

Question :

En vous aidant de la documentation de NearestNeighbours dans scikit-learn, expliquez ce que fait le code ci-dessus.

Correction

NearestNeighbors permet d’obtenir pour chaque point du jeu de données la liste de ses k-voisins les plus proches (en l’occurrence, k=4). La méthode .kneighbors() permet ensuite d’obtenir d’une part la liste des distances (c’est-à-dire les k-distances !) et d’autre part la liste des indices des k-voisins pour chaque point de la liste passée en paramètre. distances contient donc la liste des 4-distances pour chaque point de la matrice d’observations data.

Nous pouvons maintenant tracer le graphe des 4-distances. Pour ce faire, on ne conserve que la distance de chaque point à son quatrième voisin, puis on trie cette liste par ordre décroissant.

distances_triees = np.sort(distances[:,-1])[::-1]
plt.plot(distances_triees)
plt.xlabel("Nombre de points")
plt.ylabel("$4$-distance")
plt.show()

Question :

À partir du graphe des 4-distances, déterminez la valeur de eps adaptée pour ce jeu de données en utilisant l’heuristique vue en cours. Appliquez de nouveau DBSCAN avec ces paramètres. Affichez le nuage de points obtenu.

Correction :

db = DBSCAN(eps=0.18, min_samples=5)
predictions = db.fit_predict(data)

plt.scatter(data[:,0], data[:,1], c=predictions)
plt.title("Partitionnement obtenu par DBSCAN ($\\varepsilon=0.15, m=4$)")
plt.plot()

Question :

Combien obtenez vous de groupes ? À quoi correspondent les observations d’étiquette -1 ?

Correction

On obtient désormais deux groupes, comme attendu. Les observations d’étiquette -1 sont les points aberrants identifiés par DBSCAN (les outliers) qui n’appartiennent à aucun groupe.

Classification avec DBSCAN des données « Iris »

Le jeu de données Iris est est un grand classique de l’apprentissage statistique. Il comporte 150 observations de plantes selon quatre attributs :

  • longueur de la sépale,

  • largeur de la sépale,

  • longueur de la pétale,

  • largeur de la pétale.

Les observations appartiennent à l’une des trois classes, correspondant aux trois espèces d’iris: Setosa, Versicolour ou Virginica.

Plus de détails peuvent se trouver dans la documentation.

Iris est intégré à scikit-learn et est disponible depuis le sous-module sklearn.datasets:

from sklearn.datasets import load_iris

X, y = load_iris(return_X_y=True)

Question

Combien ce jeu de données contient-il d’observations?

Correction

print(X.shape)

Afin de simuler la présence de données aberrantes, nous allons générer aléatoirement 20 points bruitées, tirées selon une loi uniforme entre les valeurs minimales et maximales de chaque colonne de la matrice d’observation :

min_, max_ = X.min(axis=0), X.max(axis=0)
noise = np.random.rand(20, 4) * (max_ - min_) + min_
X = np.concatenate((X, noise))
y = np.concatenate((y, -1 * np.ones(20, dtype=int)))

Question :

Réalisez une analyse en composantes principales et visualisez le jeu de données Iris projeté selon ses deux premiers axes principaux.

Correction

from sklearn.decomposition import PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)
plt.scatter(X_pca[:,0], X_pca[:,1], c=y) and plt.show()

Question :

Appliquez une classification automatique à l’aide de DBSCAN (travaillez sur les données en dimension 4, pas les données projetées !). Visualisez les groupes obtenus dans le plan principal. Comparez ce résultat au partitionnement obtenu à l’aide d’un k-means.

Correction :

from sklearn.cluster import DBSCAN

neighbors = NearestNeighbors(n_neighbors=7)
neighbors_fit = neighbors.fit(X)
distances, indices = neighbors_fit.kneighbors(X)

db = DBSCAN(eps=0.8, min_samples=8)
labels = db.fit_predict(X)

print(np.unique(labels))

fig = plt.figure(figsize=(16, 6))
fig.add_subplot(121)
plt.plot(-np.sort(-distances[:,-1]))
plt.xlabel("Nombre de points")
plt.ylabel("$4$-distance")
plt.title("Graphe des $4$-distances")
fig.add_subplot(122)
plt.scatter(x=X_pca[:,0], y=X_pca[:,1], c=labels)
plt.plot()

km = KMeans(n_clusters=3)
labels = km.fit_predict(X)

plt.scatter(X_pca[:,0], X_pca[:,1], c=labels) and plt.show()

Les points correspondant à l’étiquette -1 sont ceux ayant été identifiés comme données aberrantes (outliers). Si la valeur de \(\varepsilon\) a bien été choisie, les observations bruitées injectées dans le jeu de données devraient ainsi avoir été isolées par DBSCAN.

Question :

À l’aide des fonctions de sklearn.metrics, calculez le taux de bonne détection des outliers. Jusqu’à quelle proportion de données bruitées le partitionnement obtenu par DBSCAN est-il robuste ?

Correction :

from sklearn.metrics import accuracy_score

predicted_outliers = db.labels_
predicted_outliers[db.labels_ > 0] = 0
true_outliers = y
true_outliers[y > 0] = 0

print(accuracy_score(true_outliers, predicted_outliers))