Résumé: Illustration des algorithmes de classification non supervisée de Scikit-learn sur des données "jouet": réallocation dynamique ($k$-means), DBSCAN, mélange gaussien.
Tutoriel adapté de: https://github.com/ageron/handson-ml/blob/master/08_dimensionality_reduction.ipynb
Machine Learning avec Scikit-Learn, Mise en oeuvre et cas concrets, par Aurélien Géron.
Ce tutoriel permet d'illustrer sur des données simulées ou "jouet"les principales fonctions de classification non supervisées présentent dans la librairie Scikit-learn. L'algorithme $k$-means y est particulièrement développé ainsi que différentes versions de DBSCAN et les modèles de mélange gaussien. En revanche, la classification ascendante hiérarchique est à rechercher dans la librairie scipy L'usage de cet algorrithme est illustré dans un autre tutoriel.
In [ ]:
import numpy as np
# initialisation du générateur de nombres aléatoires
np.random.seed(42)
# Graphiques
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
# Ignorer les warnings inutiles (cf. SciPy issue #5998)
import warnings
warnings.filterwarnings(action="ignore", module="scipy", message="^internal gelsd")
In [ ]:
from sklearn.datasets import make_blobs
In [ ]:
# Définition des centres
blob_centers = np.array(
[[ 0.2, 2.3],
[-1.5 , 2.3],
[-2.8, 1.8],
[-2.8, 2.8],
[-2.8, 1.3]])
# Définition des dispersions
blob_std = np.array([0.4, 0.3, 0.1, 0.1, 0.1])
# Générations des *blobs*
X, y = make_blobs(n_samples=2000, centers=blob_centers, cluster_std=blob_std, random_state=7)
Représentation graphique.
In [ ]:
def plot_clusters(X, y=None):
plt.scatter(X[:, 0], X[:, 1], c=y, s=1)
plt.xlabel("$x_1$", fontsize=14)
plt.ylabel("$x_2$", fontsize=14, rotation=0)
In [ ]:
plt.figure(figsize=(8, 4))
plot_clusters(X)
plt.show()
In [ ]:
from sklearn.cluster import KMeans
In [ ]:
k = 5
# Choix des options
kmeans = KMeans(n_clusters=k, random_state=42)
# Exécution de l'algorithme
y_pred = kmeans.fit_predict(X)
In [ ]:
# affichage des classes
kmeans.labels_
In [ ]:
# qui sont aussi dans la variables `y_pred`
y_pred
La fonction fit_predict calcule aussi l'affectation à la classe dont le barycentre est le plus proche pour de nouvelles observations.
In [ ]:
X_new = np.array([[0, 2], [3, 2], [-3, 3], [-3, 2.5]])
kmeans.predict(X_new)
La fonction transform calcule elle les distances aux barycentres plutôt que l'affectation.
In [ ]:
kmeans.transform(X_new)
Affichage des barycentres des classes
In [ ]:
kmeans.cluster_centers_
In [ ]:
def plot_data(X):
plt.plot(X[:, 0], X[:, 1], 'k.', markersize=2)
def plot_centroids(centroids, weights=None, circle_color='w', cross_color='k'):
if weights is not None:
centroids = centroids[weights > weights.max() / 10]
plt.scatter(centroids[:, 0], centroids[:, 1],
marker='o', s=30, linewidths=8,
color=circle_color, zorder=10, alpha=0.9)
plt.scatter(centroids[:, 0], centroids[:, 1],
marker='x', s=50, linewidths=50,
color=cross_color, zorder=11, alpha=1)
def plot_decision_boundaries(clusterer, X, resolution=1000, show_centroids=True,
show_xlabels=True, show_ylabels=True):
mins = X.min(axis=0) - 0.1
maxs = X.max(axis=0) + 0.1
xx, yy = np.meshgrid(np.linspace(mins[0], maxs[0], resolution),
np.linspace(mins[1], maxs[1], resolution))
Z = clusterer.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
plt.contourf(Z, extent=(mins[0], maxs[0], mins[1], maxs[1]),
cmap="Pastel2")
plt.contour(Z, extent=(mins[0], maxs[0], mins[1], maxs[1]),
linewidths=1, colors='k')
plot_data(X)
if show_centroids:
plot_centroids(clusterer.cluster_centers_)
if show_xlabels:
plt.xlabel("$x_1$", fontsize=14)
else:
plt.tick_params(labelbottom=False)
if show_ylabels:
plt.ylabel("$x_2$", fontsize=14, rotation=0)
else:
plt.tick_params(labelleft=False)
In [ ]:
plt.figure(figsize=(8, 4))
plot_decision_boundaries(kmeans, X)
plt.show()
Ce graphique montre bien que l'algorithme est typiquement basé sur des distances euclidiennes et ne prend pas en compte la dispersion plus ou moins grande des classes. Certains points à la frontière ne semblent pas affectés à la bonne classe.
In [ ]:
def plot_clusterer_comparison(clusterer1, clusterer2, X, title1=None, title2=None):
clusterer1.fit(X)
clusterer2.fit(X)
plt.figure(figsize=(10, 3.2))
plt.subplot(121)
plot_decision_boundaries(clusterer1, X)
if title1:
plt.title(title1, fontsize=14)
plt.subplot(122)
plot_decision_boundaries(clusterer2, X, show_ylabels=False)
if title2:
plt.title(title2, fontsize=14)
Chaque exécution est limitée à une seule initilisation.
In [ ]:
kmeans_rnd_init1 = KMeans(n_clusters=5, init="random", n_init=1, algorithm="full", random_state=11)
kmeans_rnd_init2 = KMeans(n_clusters=5, init="random", n_init=1, algorithm="full", random_state=19)
plot_clusterer_comparison(kmeans_rnd_init1, kmeans_rnd_init2, X,
"Solution 1", "Solution 2 avec différentes initialisations")
plt.show()
L'évaluation de l'inertie du nuage permet de comparer les qualités de ces classifications.
In [ ]:
kmeans.inertia_
Inertie qui est la somme des carrés des distances entre chaque observation et le barycentre de sa classe.
In [ ]:
X_dist = kmeans.transform(X)
np.sum(X_dist[np.arange(len(X_dist)), kmeans.labels_]**2)
In [ ]:
kmeans.score(X)
Alors que le score() est moins l'inertie.
Une meilleure solution serait celle qui minimise l'inertie. La première exécution est donc meilleure car par défaut l'algorithme lance 10 exécutions et conserve celle qui minimise l'inertie.
In [ ]:
kmeans_rnd_init1.inertia_
In [ ]:
kmeans_rnd_init2.inertia_
Exécution avec 10 initialisations
In [ ]:
kmeans_rnd_10_inits = KMeans(n_clusters=5, init="random", n_init=10, algorithm="full", random_state=11)
kmeans_rnd_10_inits.fit(X)
qui fournit sans doute la meilleure classification connaissant le bon nombre de classes.
In [ ]:
plt.figure(figsize=(8, 4))
plot_decision_boundaries(kmeans_rnd_10_inits, X)
plt.show()
Arthur and Vassilvitskii (2006) ont proposé une autre initialisation des barycentres de l'algorithme:
En principe ce procédé permet de réduire le nombre n_init d'initialisations et de donc de compenser le temps pris par cette procédure d'initialisation.
Remarque: cette initialisation est obtenue par init="k-means++" qui est prise par défaut. Une "initialisation optimale"conduit à des résultats identiques.
In [ ]:
KMeans()
In [ ]:
good_init = np.array([[-3, 3], [-3, 2], [-3, 1], [-1, 2], [0, 2]])
kmeans = KMeans(n_clusters=5, init=good_init, n_init=1, random_state=42)
kmeans.fit(X)
kmeans.inertia_
Elkan (2003) a proposé une accélération de l'algorithme en exploitant l'inégalité triangulaire pour économiser des calculs de distance: pour trois points A, B et C, la distance AC ≤ AB + BC. Cela permet de garder en mémoire des bornes inféreures et supérieures des distances entre observations et barycentres. L'option algorithm="elkan" exécute ce procédé par défaut mais seulement pour des données denses. Pour des données creuses (sparse) l'option full est utilisée par défaut.
In [ ]:
%timeit -n 50 KMeans(algorithm="elkan").fit(X)
In [ ]:
%timeit -n 50 KMeans(algorithm="full").fit(X)
In [ ]:
kmeans_k3 = KMeans(n_clusters=3, random_state=42)
kmeans_k8 = KMeans(n_clusters=8, random_state=42)
plot_clusterer_comparison(kmeans_k3, kmeans_k8, X, "$K=3$", "$K=8$")
plt.show()
In [ ]:
kmeans_k3.inertia_
In [ ]:
kmeans_k8.inertia_
Mais le critère d'inertie ne permet pas de choisir car celle-ci décroît mécaniquement avec le nombre de classes. Nénamoins la recherche d'un coude dans cette décroissance peut être une indication d'un nombre minimal de classes.
In [ ]:
kmeans_per_k = [KMeans(n_clusters=k, random_state=42).fit(X)
for k in range(1, 10)]
inertias = [model.inertia_ for model in kmeans_per_k]
In [ ]:
plt.figure(figsize=(8, 3.5))
plt.plot(range(1, 10), inertias, "bo-")
plt.xlabel("$k$", fontsize=14)
plt.ylabel("Inertia", fontsize=14)
plt.axis([1, 8.5, 0, 1300])
plt.show()
Attention, ces données "jouets" sont relativement "faciles" à segmenter en classes. des données réelles ne conduisent pas à des graphes aussi limpides. Le choix de 4 classes n'est pas mauvais en soit; la principale difficulté pour segmenter ces données réside dans les dispersions différentes des groupes.
Tracé des silouhettes des classes en fonction de $K$.
In [ ]:
from sklearn.metrics import silhouette_score
In [ ]:
silhouette_score(X, kmeans.labels_)
In [ ]:
silhouette_scores = [silhouette_score(X, model.labels_)
for model in kmeans_per_k[1:]]
In [ ]:
plot_decision_boundaries(kmeans_per_k[4-1], X)
plt.show()
In [ ]:
plt.figure(figsize=(8, 3))
plt.plot(range(2, 10), silhouette_scores, "bo-")
plt.xlabel("$k$", fontsize=14)
plt.ylabel("Silhouette score", fontsize=14)
plt.axis([1.8, 8.5, 0.55, 0.7])
plt.show()
In [ ]:
from sklearn.metrics import silhouette_samples
from matplotlib.ticker import FixedLocator, FixedFormatter
plt.figure(figsize=(11, 9))
for k in (3, 4, 5, 6):
plt.subplot(2, 2, k - 2)
y_pred = kmeans_per_k[k - 1].labels_
silhouette_coefficients = silhouette_samples(X, y_pred)
padding = len(X) // 30
pos = padding
ticks = []
for i in range(k):
coeffs = silhouette_coefficients[y_pred == i]
coeffs.sort()
cmap = matplotlib.cm.get_cmap("Spectral")
color = cmap(i / k)
plt.fill_betweenx(np.arange(pos, pos + len(coeffs)), 0, coeffs,
facecolor=color, edgecolor=color, alpha=0.7)
ticks.append(pos + len(coeffs) // 2)
pos += len(coeffs) + padding
plt.gca().yaxis.set_major_locator(FixedLocator(ticks))
plt.gca().yaxis.set_major_formatter(FixedFormatter(range(k)))
if k in (3, 5):
plt.ylabel("Cluster")
if k in (5, 6):
plt.gca().set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1])
plt.xlabel("Silhouette Coefficient")
else:
plt.tick_params(labelbottom=False)
plt.axvline(x=silhouette_scores[k - 2], color="red", linestyle="--")
plt.title("$k={}$".format(k), fontsize=16)
plt.show()
In [ ]:
X1, y1 = make_blobs(n_samples=1000, centers=((4, -4), (0, 0)), random_state=42)
X1 = X1.dot(np.array([[0.374, 0.95], [0.732, 0.598]]))
X2, y2 = make_blobs(n_samples=250, centers=1, random_state=42)
X2 = X2 + [6, -8]
X = np.r_[X1, X2]
y = np.r_[y1, y2]
In [ ]:
plot_clusters(X)
In [ ]:
kmeans_good = KMeans(n_clusters=3, init=np.array([[-1.5, 2.5], [0.5, 0], [4, 0]]), n_init=1, random_state=42)
kmeans_bad = KMeans(n_clusters=3, random_state=42)
kmeans_good.fit(X)
kmeans_bad.fit(X)
In [ ]:
plt.figure(figsize=(10, 3.2))
plt.subplot(121)
plot_decision_boundaries(kmeans_good, X)
plt.title("Inertia = {:.1f}".format(kmeans_good.inertia_), fontsize=14)
plt.subplot(122)
plot_decision_boundaries(kmeans_bad, X, show_ylabels=False)
plt.title("Inertia = {:.1f}".format(kmeans_bad.inertia_), fontsize=14)
plt.show()
In [ ]:
from sklearn.datasets import make_moons
In [ ]:
X, y = make_moons(n_samples=1000, noise=0.05, random_state=42)
In [ ]:
from sklearn.cluster import DBSCAN
In [ ]:
dbscan = DBSCAN(eps=0.05, min_samples=5)
dbscan.fit(X)
In [ ]:
dbscan.labels_[:10]
In [ ]:
len(dbscan.core_sample_indices_)
In [ ]:
dbscan.core_sample_indices_[:10]
In [ ]:
dbscan.components_[:3]
In [ ]:
np.unique(dbscan.labels_)
In [ ]:
dbscan2 = DBSCAN(eps=0.2)
dbscan2.fit(X)
In [ ]:
def plot_dbscan(dbscan, X, size, show_xlabels=True, show_ylabels=True):
core_mask = np.zeros_like(dbscan.labels_, dtype=bool)
core_mask[dbscan.core_sample_indices_] = True
anomalies_mask = dbscan.labels_ == -1
non_core_mask = ~(core_mask | anomalies_mask)
cores = dbscan.components_
anomalies = X[anomalies_mask]
non_cores = X[non_core_mask]
plt.scatter(cores[:, 0], cores[:, 1],
c=dbscan.labels_[core_mask], marker='o', s=size, cmap="Paired")
plt.scatter(cores[:, 0], cores[:, 1], marker='*', s=20, c=dbscan.labels_[core_mask])
plt.scatter(anomalies[:, 0], anomalies[:, 1],
c="r", marker="x", s=100)
plt.scatter(non_cores[:, 0], non_cores[:, 1], c=dbscan.labels_[non_core_mask], marker=".")
if show_xlabels:
plt.xlabel("$x_1$", fontsize=14)
else:
plt.tick_params(labelbottom=False)
if show_ylabels:
plt.ylabel("$x_2$", fontsize=14, rotation=0)
else:
plt.tick_params(labelleft=False)
plt.title("eps={:.2f}, min_samples={}".format(dbscan.eps, dbscan.min_samples), fontsize=14)
In [ ]:
plot_dbscan(dbscan, X, size=100)
In [ ]:
plt.figure(figsize=(9, 3.2))
plt.subplot(121)
plot_dbscan(dbscan, X, size=100)
plt.subplot(122)
plot_dbscan(dbscan2, X, size=600, show_ylabels=False)
plt.show()
In [ ]:
from sklearn.datasets import load_iris
# from IPython.display import IFrame
# IFrame('http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', width=300, height=200)
In [ ]:
data = load_iris()
X = data.data
y = data.target
data.target_names
In [ ]:
plt.figure(figsize=(9, 3.5))
plt.subplot(121)
plt.plot(X[y==0, 0], X[y==0, 1], "yo", label="Iris-Setosa")
plt.plot(X[y==1, 0], X[y==1, 1], "bs", label="Iris-Versicolor")
plt.plot(X[y==2, 0], X[y==2, 1], "g^", label="Iris-Virginica")
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.legend(fontsize=12)
plt.subplot(122)
plt.scatter(X[:, 0], X[:, 1], c="k", marker=".")
plt.xlabel("Petal length", fontsize=14)
plt.tick_params(labelleft=False)
plt.show()
Un mélange gaussien est bien adapté à ces données.
In [ ]:
from sklearn.mixture import GaussianMixture
In [ ]:
y_pred = GaussianMixture(n_components=3, random_state=42).fit(X).predict(X)
mapping = np.array([2, 0, 1])
y_pred = np.array([mapping[cluster_id] for cluster_id in y_pred])
In [ ]:
plt.plot(X[y_pred==0, 0], X[y_pred==0, 1], "yo", label="Cluster 1")
plt.plot(X[y_pred==1, 0], X[y_pred==1, 1], "bs", label="Cluster 2")
plt.plot(X[y_pred==2, 0], X[y_pred==2, 1], "g^", label="Cluster 3")
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.legend(loc="upper right", fontsize=12)
plt.show()
In [ ]:
np.sum(y_pred==y)
In [ ]:
np.sum(y_pred==y) / len(y_pred)
In [ ]:
X1, y1 = make_blobs(n_samples=1000, centers=((4, -4), (0, 0)), random_state=42)
X1 = X1.dot(np.array([[0.374, 0.95], [0.732, 0.598]]))
X2, y2 = make_blobs(n_samples=250, centers=1, random_state=42)
X2 = X2 + [6, -8]
X = np.r_[X1, X2]
y = np.r_[y1, y2]
In [ ]:
gm = GaussianMixture(n_components=3, n_init=10, random_state=42)
gm.fit(X)
Valeurs des paramètres estimés par l'algorithme.
In [ ]:
gm.weights_
In [ ]:
gm.means_
In [ ]:
gm.covariances_
La convergence est elle atteinte?
In [ ]:
gm.converged_
En combien d'itérations?
In [ ]:
gm.n_iter_
La fonction predict() permet de prédire l'appartenance d'une observation à une classe ou encore avec predict_proba(), la probabilité d'appartenance à chacune des classes.
In [ ]:
gm.predict(X)
In [ ]:
gm.predict_proba(X)
De même que pour de nouvelles observations.
In [ ]:
X_new, y_new = gm.sample(6)
X_new
In [ ]:
y_new
Estimation du logarithme de la fonction de densité entout point par score_samples()
In [ ]:
gm.score_samples(X)
Graphique des frontières des classes et des contours des densités.
In [ ]:
from matplotlib.colors import LogNorm
def plot_gaussian_mixture(clusterer, X, resolution=1000, show_ylabels=True):
mins = X.min(axis=0) - 0.1
maxs = X.max(axis=0) + 0.1
xx, yy = np.meshgrid(np.linspace(mins[0], maxs[0], resolution),
np.linspace(mins[1], maxs[1], resolution))
Z = -clusterer.score_samples(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
plt.contourf(xx, yy, Z,
norm=LogNorm(vmin=1.0, vmax=30.0),
levels=np.logspace(0, 2, 12))
plt.contour(xx, yy, Z,
norm=LogNorm(vmin=1.0, vmax=30.0),
levels=np.logspace(0, 2, 12),
linewidths=1, colors='k')
Z = clusterer.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
plt.contour(xx, yy, Z,
linewidths=2, colors='r', linestyles='dashed')
plt.plot(X[:, 0], X[:, 1], 'k.', markersize=2)
plot_centroids(clusterer.means_, clusterer.weights_)
plt.xlabel("$x_1$", fontsize=14)
if show_ylabels:
plt.ylabel("$x_2$", fontsize=14, rotation=0)
else:
plt.tick_params(labelleft=False)
In [ ]:
plt.figure(figsize=(8, 4))
plot_gaussian_mixture(gm, X)
plt.show()
En vue d'estimer des modèles plus parcimonieux, c'est-à-dire avec moins de paramètres, des contraintes peuvent être imposées aux matrices de covariance en spécifiant le paramètre covariance_type:
"full" (defaut): pas de contrainte sur les matrices de covariance, chaque ellipsoïde peut avoir des formes spécifiques, "tied": même matrices de covariance ou même forme des nauges,"spherical": nuages sphériques mais de diamètres différents,"diag": nuages ellipsoïdaux parallèles mais de dimensions différentes.
In [ ]:
gm_full = GaussianMixture(n_components=3, n_init=10, covariance_type="full", random_state=42)
gm_tied = GaussianMixture(n_components=3, n_init=10, covariance_type="tied", random_state=42)
gm_spherical = GaussianMixture(n_components=3, n_init=10, covariance_type="spherical", random_state=42)
gm_diag = GaussianMixture(n_components=3, n_init=10, covariance_type="diag", random_state=42)
gm_full.fit(X)
gm_tied.fit(X)
gm_spherical.fit(X)
gm_diag.fit(X)
In [ ]:
def compare_gaussian_mixtures(gm1, gm2, X):
plt.figure(figsize=(9, 4))
plt.subplot(121)
plot_gaussian_mixture(gm1, X)
plt.title('covariance_type="{}"'.format(gm1.covariance_type), fontsize=14)
plt.subplot(122)
plot_gaussian_mixture(gm2, X, show_ylabels=False)
plt.title('covariance_type="{}"'.format(gm2.covariance_type), fontsize=14)
In [ ]:
compare_gaussian_mixtures(gm_tied, gm_spherical, X)
plt.show()
In [ ]:
compare_gaussian_mixtures(gm_full, gm_diag, X)
plt.tight_layout()
plt.show()
Comparer ces graphiques, lequel semble le plus adapté à la forme des nuages?
Comme la clasisficaiton ascendante hiérarchique, les mélanges gaussiens peuvent être utilisés pour détecter des anomalies ou observations atypiques c'est-à-dire présentes dans des zones de faible densité. Un paramètre doit être réglé, il définit par exemple le taux d'anomalies à détecter.
In [ ]:
densities = gm.score_samples(X)
density_threshold = np.percentile(densities, 4)
anomalies = X[densities < density_threshold]
In [ ]:
plt.figure(figsize=(8, 4))
plot_gaussian_mixture(gm, X)
plt.scatter(anomalies[:, 0], anomalies[:, 1], color='r', marker='*')
plt.ylim(ymax=5.1)
plt.show()
L'inertie, qui décroît avec le nombre de paramètres du modèles donc de classes, de même que le score de silouhette adapté à des nuages sphériques ne sont pas adaptés à l'optimisation dex choix de modèles: nombre de classes et contraintes sur les matrices de covariance. Le modèle de mélange est un modèle statistique usuel estimé par maximum de la log vraisemblance. La stratégie de sélection de modèle par pénalisation de la vraisemblance du modèle linéaire général peut être utilisée pour le modèle de mélange gaussien. Cette sytratégie conduit à la minimisation des critères BIC (Bayesian Information Criterion) et AIC (Akaike Information Criterion) qui introduisent une pénalisation fonction du nombre de paramètres.
${BIC} = {\log(m)p - 2\log({\hat L})}$
${AIC} = 2p - 2\log(\hat L)$
Minimiser le BIC ou l'AIC conduit trouver un "meilleur" équilibre entre qualité d'ajustement du modèle (Biais) et variance des paramètres.
In [ ]:
gm.bic(X)
In [ ]:
gm.aic(X)
Let's train Gaussian Mixture models with various values of $k$ and measure their BIC:
In [ ]:
gms_per_k = [GaussianMixture(n_components=k, n_init=10, random_state=42).fit(X)
for k in range(1, 11)]
In [ ]:
bics = [model.bic(X) for model in gms_per_k]
aics = [model.aic(X) for model in gms_per_k]
In [ ]:
plt.figure(figsize=(8, 3))
plt.plot(range(1, 11), bics, "bo-", label="BIC")
plt.plot(range(1, 11), aics, "go--", label="AIC")
plt.xlabel("$k$", fontsize=14)
plt.ylabel("Information Criterion", fontsize=14)
plt.axis([1, 9.5, np.min(aics) - 50, np.max(aics) + 50])
plt.legend()
plt.show()
Chercher le meilleur choix pour le nombre de classes et le covariance_type:
In [ ]:
min_bic = np.infty
for k in range(1, 11):
for covariance_type in ("full", "tied", "spherical", "diag"):
bic = GaussianMixture(n_components=k, n_init=10,
covariance_type=covariance_type,
random_state=42).fit(X).bic(X)
if bic < min_bic:
min_bic = bic
best_k = k
best_covariance_type = covariance_type
In [ ]:
best_k
In [ ]:
best_covariance_type
In [ ]: