KMeans 聚类算法
目录
在介绍 KMeans 文本聚类 后，我们此篇内容对 KMeans 算法做进一步详细介绍。因关于 KMeans 等经典算法的介绍，无论是原理篇还是实战篇都数不胜数，本篇内容节选自 《Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow, 2nd Edition》的第九章内容，该书无论是内容的编排上还是知识点的覆盖程度上都算的是佳作，是初学者的不二之选。
0、环境准备
加载一些必要的包、绘图设置及图片保存函数等内容
# Python ≥3.5 is required
import sys
assert sys.version_info >= (3, 5)
# Scikit-Learn ≥0.20 is required
import sklearn
assert sklearn.__version__ >= "0.20"
# Common imports
import numpy as np
import os
# to make this notebook's output stable across runs
np.random.seed(42)
# To plot pretty figures
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rc('axes', labelsize=14)
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)
# Where to save the figures
PROJECT_ROOT_DIR = "."
CHAPTER_ID = "clustering"
IMAGES_PATH = os.path.join(PROJECT_ROOT_DIR, "images", CHAPTER_ID)
os.makedirs(IMAGES_PATH, exist_ok=True)
def save_fig(fig_id, tight_layout=True, fig_extension="png", resolution=300):
path = os.path.join(IMAGES_PATH, fig_id + "." + fig_extension)
print("Saving figure", fig_id)
if tight_layout:
plt.tight_layout()
plt.savefig(path, format=fig_extension, dpi=resolution)
# Ignore useless warnings (see SciPy issue #5998)
import warnings
warnings.filterwarnings(action="ignore", message="^internal gelsd")
1、聚类 vs 分类
下面使用的数据集为经典的 iris 鸢尾花数据集，Iris数据集是常用的分类实验数据集，由 Fisher 在 1936 年收集整理。Iris 也称鸢尾花卉数据集，是一类多重变量分析的数据集。数据集包含 150 个数据样本，分为 3 类，每类 50 个数据，每个数据包含 4 个属性。可通过花萼长度，花萼宽度，花瓣长度，花瓣宽度 4 个属性预测鸢尾花卉属于（Setosa，Versicolour，Virginica）三个种类中的哪一类。
该数据集在 scikit-learn 下有提供，我们可以直接加载，并将数据分布可视化呈现出来。
from sklearn.datasets import load_iris
data = load_iris()
X = data.data
y = data.target
data.target_names
array(['setosa', 'versicolor', 'virginica'], dtype='<U10')
plt.figure(figsize=(9, 3.5))
plt.subplot(121)
plt.plot(X[y==0, 2], X[y==0, 3], "yo", label="Iris setosa")
plt.plot(X[y==1, 2], X[y==1, 3], "bs", label="Iris versicolor")
plt.plot(X[y==2, 2], X[y==2, 3], "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[:, 2], X[:, 3], c="k", marker=".")
plt.xlabel("Petal length", fontsize=14)
plt.tick_params(labelleft=False)
save_fig("classification_vs_clustering_plot")
plt.show()
上图中，左侧是通过真实标签将不同 label 的数据点用不同的颜色标记出来了，而右侧则是去掉真实标签信息的数据分布情况。一般的算法对左小角这一小堆都可以很好地聚类出来，但是上面两个簇却未必可以很好地分开。图中为了方便可视化只选择了其中的两个属性，还有另外 2 个属性（萼片的长度和宽度）。聚类算法可以很好地利用数据的所有特征信息，例如高斯混合模型使用所有特征信息可以近乎完美的将三种鸢尾花区分开。（150 个数据只有 5 个聚类错误）
from sklearn.mixture import GaussianMixture
y_pred = GaussianMixture(n_components=3, random_state=42).fit(X).predict(X)
mapping = np.array([1, 2, 0])
y_pred = np.array([mapping[cluster_id] for cluster_id in y_pred])
plt.plot(X[y_pred==0, 2], X[y_pred==0, 3], "yo", label="Cluster 1")
plt.plot(X[y_pred==1, 2], X[y_pred==1, 3], "bs", label="Cluster 2")
plt.plot(X[y_pred==2, 2], X[y_pred==2, 3], "g^", label="Cluster 3")
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.legend(loc="upper left", fontsize=12)
plt.show()
np.sum(y_pred==y)
# 145
np.sum(y_pred==y) / len(y_pred)
# 0.9666666666666667
2、K-Means
K-Means算法是一种简单的算法，能够非常快速和高效地对这类数据集进行聚类，通常只需几次迭代。1957年，Stuart Lloyd 在贝尔实验室提出了一种脉冲编码调制技术，但1982年才在公司外发表。1965年，Edward W.Forgy 发布了几乎相同的算法，因此 K-均值有时被称为 Lloyd-Forgy。
为了进行演示，下面我们随机生成 2000 个数据点，并有目的的设置 5 个中心点，具体见下图：
from sklearn.datasets import make_blobs
blob_centers = np.array(
[[ 0.2, 2.3],
[-1.5 , 2.3],
[-2.8, 1.8],
[-2.8, 2.8],
[-2.8, 1.3]])
blob_std = np.array([0.4, 0.3, 0.1, 0.1, 0.1])
X, y = make_blobs(n_samples=2000, centers=blob_centers,
cluster_std=blob_std, random_state=7)
将其可视化：
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)
plt.figure(figsize=(8, 4))
plot_clusters(X)
save_fig("blobs_plot")
plt.show()
2.1、数据拟合和预测
下面，我们使用 K-Means 算法对上述数据进行聚类分析，算法将尝试寻找每个簇的中心点，并将每个实例指派到距离最近的一个簇中
from sklearn.cluster import KMeans
k = 5
kmeans = KMeans(n_clusters=k, random_state=42)
y_pred = kmeans.fit_predict(X)
每个实例都被指派到 5 个簇中的其中一个
y_pred
# array([4, 1, 0, ..., 3, 0, 1], dtype=int32)
y_pred is kmeans.labels_
# True
下面我们可以打印出每个簇的中点点坐标信息：
kmeans.cluster_centers_
array([[ 0.20876306, 2.25551336],
[-2.80389616, 1.80117999],
[-1.46679593, 2.28585348],
[-2.79290307, 2.79641063],
[-2.80037642, 1.30082566]])
这样一来，我们也可以对新数据进行类似分类算法的预测，即将新实例指派到 5 个簇中的其中一个。
X_new = np.array([[0, 2], [3, 2], [-3, 3], [-3, 2.5]])
kmeans.predict(X_new)
# array([0, 0, 3, 3], dtype=int32)
2.2、画出决策边界
下面我们可以使用 Voronoi 划分算法给出每个簇之间的边界：
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)
plt.figure(figsize=(8, 4))
plot_decision_boundaries(kmeans, X)
save_fig("voronoi_plot")
plt.show()
# 该值下文中进行解释
kmeans.inertia_
# 211.5985372581684
Voronoi 划分后，绝大多数实例都被清楚地分配到适当的簇中，但也有少数实例可能被错误地标记（尤其是在左上角的簇和中间的簇之间的边界附近）。事实上，K-Means 算法在 blob 的直径非常不同的情况下表现得并不好，因为它在为簇分配实例时只关心到质心的距离。
2.3、硬聚类 vs 软聚类
与其将每个实例分配给一个单独的簇（称为硬聚类），不如为每个实例在每个簇下分配一个分数，这称为软聚类。分数可以是实例与质心之间的距离；反过来，它也可以是相似性分数（或亲和力），例如高斯径向基函数。在 KMeans 中，transform() 方法测量从每个实例到每个质心的距离：
kmeans.transform(X_new)
array([[0.32995317, 2.81093633, 1.49439034, 2.9042344 , 2.88633901],
[2.80290755, 5.80730058, 4.4759332 , 5.84739223, 5.84236351],
[3.29399768, 1.21475352, 1.69136631, 0.29040966, 1.71086031],
[3.21806371, 0.72581411, 1.54808703, 0.36159148, 1.21567622]])
可以很容易验证，这其实就是每个数据点与每个簇中心点之间的欧氏距离，而且通过这种技术实际可以对高维数据集做数据降维。例如原始数据集是几百维的甚至更高，经过降维后，数据仅 5 维。
np.linalg.norm(np.tile(X_new, (1, k)).reshape(-1, k, 2) - kmeans.cluster_centers_, axis=2)
array([[0.32995317, 2.81093633, 1.49439034, 2.9042344 , 2.88633901],
[2.80290755, 5.80730058, 4.4759332 , 5.84739223, 5.84236351],
[3.29399768, 1.21475352, 1.69136631, 0.29040966, 1.71086031],
[3.21806371, 0.72581411, 1.54808703, 0.36159148, 1.21567622]])
2.4、K-Means 算法详解
下面我们详细看一下 K-Means 算法的工作原理，K-Means 算法是速度最快的聚类算法中的一种，也是最为简单的一种，该算法的计算复杂度通常与实例数 m m m、簇数 k k k 和维数 n n n 成线性关系，但只有当数据具有聚类结构时才是如此。如果没有，那么在最坏的情况下，复杂度会随着实例的数量呈指数级增加。实际上，这种情况很少发生，K-Means 通常是最快的聚类算法之一。具体原理：
- 首先随机初始化 k k k 个聚类中心点，表明我们最终需要将所有数据指派到这 k k k 个簇下
- 计算每个簇下所有数据点的平均值，重新更新簇中心点的位置
在 scikit-learn 下，默认的 KMeans 算法已经是被优化过的，如果出于教学或者其他演示目的，你需要将其中的部分参数设置为 init="random", n_init=1以及 algorithm="full"，这些超参数我们下面会进行解释
我们下面给出上述数据分别在 1、2 和 3 次迭代更新下的簇中心点位移变化，实际上仅 3 次迭代就已经接近比较理想的结果了
kmeans_iter1 = KMeans(n_clusters=5, init="random", n_init=1,
algorithm="full", max_iter=1, random_state=1)
kmeans_iter2 = KMeans(n_clusters=5, init="random", n_init=1,
algorithm="full", max_iter=2, random_state=1)
kmeans_iter3 = KMeans(n_clusters=5, init="random", n_init=1,
algorithm="full", max_iter=3, random_state=1)
kmeans_iter1.fit(X)
kmeans_iter2.fit(X)
kmeans_iter3.fit(X)
KMeans(algorithm='full', copy_x=True, init='random', max_iter=3, n_clusters=5,
n_init=1, n_jobs=None, precompute_distances='auto', random_state=1,
tol=0.0001, verbose=0)
plt.figure(figsize=(10, 8))
plt.subplot(321)
plot_data(X)
plot_centroids(kmeans_iter1.cluster_centers_, circle_color='r', cross_color='w')
plt.ylabel("$x_2$", fontsize=14, rotation=0)
plt.tick_params(labelbottom=False)
plt.title("Update the centroids (initially randomly)", fontsize=14)
plt.subplot(322)
plot_decision_boundaries(kmeans_iter1, X, show_xlabels=False, show_ylabels=False)
plt.title("Label the instances", fontsize=14)
plt.subplot(323)
plot_decision_boundaries(kmeans_iter1, X, show_centroids=False, show_xlabels=False)
plot_centroids(kmeans_iter2.cluster_centers_)
plt.subplot(324)
plot_decision_boundaries(kmeans_iter2, X, show_xlabels=False, show_ylabels=False)
plt.subplot(325)
plot_decision_boundaries(kmeans_iter2, X, show_centroids=False)
plot_centroids(kmeans_iter3.cluster_centers_)
plt.subplot(326)
plot_decision_boundaries(kmeans_iter3, X, show_ylabels=False)
save_fig("kmeans_algorithm_plot")
plt.show()
2.5、K-Means 易变性
虽然该算法可以保证收敛，但它可能无法收敛到正确的解（即，它可能会收敛到局部最优解）：是否收敛取决于质心初始化。下图显示了两个次优的解决方案。
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)
kmeans_rnd_init1 = KMeans(n_clusters=5, init="random", n_init=1, algorithm="full")
kmeans_rnd_init2 = KMeans(n_clusters=5, init="random", n_init=1, algorithm="full")
plot_clusterer_comparison(kmeans_rnd_init1, kmeans_rnd_init2, X,
"Solution 1", "Solution 2 (with a different random init)")
save_fig("kmeans_variability_plot")
plt.show()
2.6、Inertia
那怎么解决这个初始化的问题呢？我们先介绍一个简单的，例如，我们通过运行其他聚类算法或者通过一些可视化工作，可大致判定各个簇的中心位置，可我们可以直接指定：
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)
print(kmeans_rnd_init1.inertia_, kmeans_rnd_init2.inertia_)
# 236.9724856977868 219.85933533452254
另一个解决方案是使用不同的随机初始化多次运行该算法，并保持最佳解。随机初始化的数量由 n_init 超参数控制：默认情况下等于10，这意味着前面描述的整个算法在调用 fit() 时运行 10 次，scikit-learn 保持最优解。但它究竟如何知道哪种解决方案是最好的呢？它使用性能指标！这个度量称为模型的惯性，它是每个实例与其最近质心之间的均方根距离。上图中左边的模型大致等于236.97，右边的模型大约等于 219.86，而上面一开始 kmeans 模型为 211.60 左右。
很容易验证，惯性是每个训练实例与其最近质心之间的平方距离之和：
X_dist = kmeans.transform(X)
np.sum(X_dist[np.arange(len(X_dist)), kmeans.labels_]**2)
# 211.59853725816856
score() 方法返回负惯性。为什么是负的？好吧，这是因为预测器的 score() 方法必须始终遵守 “great is better” 的规则。
kmeans.score(X)
# -211.59853725816856
所以我们可以通过调整超参 n_init 的值，通过多次初始化进行模型选择。下面我们介绍一个对 K-Means 算法具有重要提高作用的算法 K-Means++ 算法。
2.7、K-Means++
K-Means++ 由 David Arthur 和 Sergei Vassilvitskii 在 2006 年发表的论文中提到，它不在采用完全的随机初始化方法，而是：
- 从数据集中随机的均匀分布中选择一个中心点 c 1 c_1 c1
- 确定一个新的中心点 c i c_i ci，具体做法是以某种概率选择一个实例 x i \mathbf{x}_i xi，具体的概率计算方式为： D ( x i ) 2 D(\mathbf{x}_i)^2 D(xi)2 / ∑ j = 1 m D ( x j ) 2 \sum\limits_{j=1}^{m}{D(\mathbf{x}_j)}^2 j=1∑mD(xj)2 这里的 D ( x i ) D(\mathbf{x}_i) D(xi) 是实例 x i \mathbf{x}_i xi 与之前已选择的中心点中最近的距离。这种概率分布确保了距离已经选择的质心较远的实例更有可能被选为质心。
- 重复上述步骤，直到 k k k 个中心点被选择
K-Means++ 算法的其余部分只是普通的 K-Means。通过这种初始化，K-Means 算法不太可能收敛到一个次优解，因此可以大大减少 n_init。大多数时候，这在很大程度上弥补了初始化过程的额外复杂性。
可以通过设置 init="k-means++" 来使用 k-means++ 算法进行初始化，实际上这也是默认的初始化方法。
KMeans()
KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
n_clusters=8, n_init=10, n_jobs=None, precompute_distances='auto',
random_state=None, tol=0.0001, verbose=0)
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_
# 211.5985372581684
2.8、K-Means 加速
通过避免许多不必要的距离计算，K-Means 算法可以大大加快速度：这是通过利用三角形不等式（给定三个点 A、B 和 C，距离 AC 始终是 AC≤AB+BC）和跟踪实例和质心之间距离的上下界（见下文）来实现的，更多细节详见 Charles Elkan 的2003年论文。
为了使用 Elkan 的变体 K-Means，仅需要设置 algorithm="elkan"，注意，默认不支持稀疏数据，因此对于稠密数据，scikit-learn 默认使用 "elkan" 而稀疏数据默认使用 "full" (正则化 K-Means 算法)
%timeit -n 50 KMeans(algorithm="elkan").fit(X)
# 50 loops, best of 3: 101 ms per loop
%timeit -n 50 KMeans(algorithm="full").fit(X)
# 50 loops, best of 3: 139 ms per loop
2.9、Mini-Batch K-Means
Scikit-Learn 也实现了一个支持 mini-batch 的 K-Means 加速算法（详见这篇论文）:
from sklearn.cluster import MiniBatchKMeans
minibatch_kmeans = MiniBatchKMeans(n_clusters=5, random_state=42)
minibatch_kmeans.fit(X)
MiniBatchKMeans(batch_size=100, compute_labels=True, init='k-means++',
init_size=None, max_iter=100, max_no_improvement=10,
n_clusters=5, n_init=3, random_state=42,
reassignment_ratio=0.01, tol=0.0, verbose=0)
minibatch_kmeans.inertia_
# 211.93186531476775
如果数据集比较大，无法直接在内存中拟合，最简单的方式是使用 memmap，如同 PCA 算法一样，例如我们下面以 MNIST 数据集为例：
import urllib
from sklearn.datasets import fetch_openml
mnist = fetch_openml('mnist_784', version=1)
mnist.target = mnist.target.astype(np.int64)
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(
mnist["data"], mnist["target"], random_state=42)
接下来把它写入到 memmap:
filename = "my_mnist.data"
X_mm = np.memmap(filename, dtype='float32', mode='write', shape=X_train.shape)
X_mm[:] = X_train
minibatch_kmeans = MiniBatchKMeans(n_clusters=10, batch_size=10, random_state=42)
minibatch_kmeans.fit(X_mm)
MiniBatchKMeans(batch_size=10, compute_labels=True, init='k-means++',
init_size=None, max_iter=100, max_no_improvement=10,
n_clusters=10, n_init=3, random_state=42,
reassignment_ratio=0.01, tol=0.0, verbose=0)
如果你的数据太大，那么你就不能使用 memmap了。让我们先写一个函数来加载下一个批量的数据（在现实生活中，你将从磁盘中加载数据）：
def load_next_batch(batch_size):
return X[np.random.choice(len(X), batch_size, replace=False)]
现在我们可以训练模型了，一次只喂一个批量的数据。我们还需要实现多次初始化，并使模型保持最小的惯性：
np.random.seed(42)
k = 5
n_init = 10
n_iterations = 100
batch_size = 100
init_size = 500 # more data for K-Means++ initialization
evaluate_on_last_n_iters = 10
best_kmeans = None
for init in range(n_init):
minibatch_kmeans = MiniBatchKMeans(n_clusters=k, init_size=init_size)
X_init = load_next_batch(init_size)
minibatch_kmeans.partial_fit(X_init)
minibatch_kmeans.sum_inertia_ = 0
for iteration in range(n_iterations):
X_batch = load_next_batch(batch_size)
minibatch_kmeans.partial_fit(X_batch)
if iteration >= n_iterations - evaluate_on_last_n_iters:
minibatch_kmeans.sum_inertia_ += minibatch_kmeans.inertia_
if (best_kmeans is None or
minibatch_kmeans.sum_inertia_ < best_kmeans.sum_inertia_):
best_kmeans = minibatch_kmeans
best_kmeans.score(X)
# -211.70999744411483
Mini-batch K-Means 比 regular K-Means 快的多:
%timeit KMeans(n_clusters=5).fit(X)
# 10 loops, best of 3: 55.9 ms per loop
%timeit MiniBatchKMeans(n_clusters=5).fit(X)
# 10 loops, best of 3: 34.7 ms per loop
虽然 MiniBatchKMeans 算法比常规 K-Means 算法快得多，但其惯性通常稍差，特别是随着簇数的增加。你可以在下中看到这一点：左边的图比较了 MiniBatchKMeans 模型和常规 K-Means 模型的惯量，这些模型使用不同数量的 K K K 对前面数据集进行训练。在右边的图中，你可以看到 MiniBatchKMeans 比常规 K-Means 快得多，并且这种差异随着 K K K 的增加而增大。
from timeit import timeit
times = np.empty((100, 2))
inertias = np.empty((100, 2))
for k in range(1, 101):
kmeans_ = KMeans(n_clusters=k, random_state=42)
minibatch_kmeans = MiniBatchKMeans(n_clusters=k, random_state=42)
print("\r{}/{}".format(k, 100), end="")
times[k-1, 0] = timeit("kmeans_.fit(X)", number=10, globals=globals())
times[k-1, 1] = timeit("minibatch_kmeans.fit(X)", number=10, globals=globals())
inertias[k-1, 0] = kmeans_.inertia_
inertias[k-1, 1] = minibatch_kmeans.inertia_
plt.figure(figsize=(10,4))
plt.subplot(121)
plt.plot(range(1, 101), inertias[:, 0], "r--", label="K-Means")
plt.plot(range(1, 101), inertias[:, 1], "b.-", label="Mini-batch K-Means")
plt.xlabel("$k$", fontsize=16)
plt.title("Inertia", fontsize=14)
plt.legend(fontsize=14)
plt.axis([1, 100, 0, 100])
plt.subplot(122)
plt.plot(range(1, 101), times[:, 0], "r--", label="K-Means")
plt.plot(range(1, 101), times[:, 1], "b.-", label="Mini-batch K-Means")
plt.xlabel("$k$", fontsize=16)
plt.title("Training time (seconds)", fontsize=14)
plt.axis([1, 100, 0, 6])
save_fig("minibatch_kmeans_vs_kmeans")
plt.show()
2.10、如何确定最优的聚类数
如果我们的聚类数比 5 大或者小呢？
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$")
save_fig("bad_n_clusters_plot")
plt.show()
这两个模型貌似不太好，我们来看一下他们的惯性值
kmeans_k3.inertia_
# 653.2167190021553
kmeans_k8.inertia_
118.41983763508077
哈哈哈，k8 的惯性值才 118，所以，我们不能简单地把惯性最小化的 k k k 值取出来，因为随着我们增加 k k k，它会不断降低。实际上，簇越多，每个实例就越接近其最近的质心，因此惯性也就越低。但是，我们可以将惯性绘制为 k k k 的函数，并分析结果曲线：
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]
plt.figure(figsize=(8, 3.5))
plt.plot(range(1, 10), inertias, "bo-")
plt.xlabel("$k$", fontsize=14)
plt.ylabel("Inertia", fontsize=14)
plt.annotate('Elbow',
xy=(4, inertias[3]),
xytext=(0.55, 0.55),
textcoords='figure fraction',
fontsize=16,
arrowprops=dict(facecolor='black', shrink=0.1)
)
plt.axis([1, 8.5, 0, 1300])
save_fig("inertia_vs_k_plot")
plt.show()
如图所示，在 k = 4 k=4 k=4 处有一个拐点，这意味着少于这个值的簇将是糟糕的，而更多的簇将不会有太大帮助，并且可能会将簇削减一半。所以 k = 4 k=4 k=4 是个不错的选择。当然，在这个例子中，它并不完美，因为它意味着左下角的两个 blob 将被视为一个单独的簇，但是它仍然是一个非常好的簇。
plot_decision_boundaries(kmeans_per_k[4-1], X)
plt.show()
另一种方法是查看 “silhouette score”，即所有实例中的平均“轮廓系数”。一个实例的轮廓系数等于 （ b − a ） / max （ a ， b ） （b-a）/\max（a，b） （b−a）/max（a，b），其中 a a a 是到同一个簇中其他实例的平均距离（它是“平均簇内距离”，而 b b b 是“平均最近簇距离”，即到下一个最近的簇实例的平均距离（定义为最小化 b b b，不包括自己的簇）。轮廓系数可以在 -1 和 +1 之间变化：接近 +1 的系数表示实例位于自己的簇内，远离其他簇；系数接近 0 表示它靠近簇边界，最后接近 -1 表示实例可能被分配到错误的簇。
from sklearn.metrics import silhouette_score
silhouette_score(X, kmeans.labels_)
# 0.655517642572828
silhouette_scores = [silhouette_score(X, model.labels_)
for model in kmeans_per_k[1:]]
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])
save_fig("silhouette_score_vs_k_plot")
plt.show()
如此所示，这种可视化比上一个要丰富得多：特别是，虽然它确认了 k = 4 k=4 k=4 是一个非常好的选择，但它也强调了这样一个事实，即 k = 5 k=5 k=5 也相当不错。
当绘制每个实例的轮廓系数时，将获得更丰富的可视化效果，并按它们所分配的集群和系数的值进行排序。这称为轮廓图。每个图包含每个簇的一个刀形状。形状的高度表示簇中包含的实例数，其宽度表示簇中实例的已排序轮廓系数（越宽越好）。虚线表示平均轮廓系数。
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()
color = mpl.cm.Spectral(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)
save_fig("silhouette_analysis_plot")
plt.show()
2.11、K-Means 的局限性
尽管 K-Means 有许多优点，最显著的是它的快速性和可扩展性，但它并不完美。为了避免算法的次优，我们可以指定若干次优解。此外，当簇具有不同大小、不同密度或非球形形状时，K-Means 表现得并不好。例如，下图显示了 K-Means 如何对包含三个不同尺寸、密度和方向的椭球簇的数据集进行聚类。
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]
plot_clusters(X)
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)
KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
n_clusters=3, n_init=10, n_jobs=None, precompute_distances='auto',
random_state=42, tol=0.0001, verbose=0)
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)
save_fig("bad_kmeans_plot")
plt.show()
如图所示，这两种解决方案都不好。左侧的解决方案更好，但它仍然会切掉中间簇的 25% 并将其分配给右侧的簇。右边的解决方案很糟糕，尽管它的惯性较低。因此，根据数据的不同，不同的聚类算法可能表现得更好。对于这类椭圆团簇，高斯混合模型非常有效。
上一篇： 索爱手机IMSI序列号获取
下一篇： 有关Cookie的一些东西。