K-means一种将数据分成K个中心组的聚类算法
K-means算法假设输入的数据点围绕K个不同的中心旋转。每个中心坐标就像一个隐藏的靶心,被散乱的数据点包围着。该算法的目的是为了发掘这些隐藏的中心坐标。
我们首先通过选择K来初始化K-means,这是我们要搜索的中心坐标的数量。在我们的镖靶分析中,K被设定为2,尽管通常K可以等于任何整数。该算法随机选择了K个数据点。这些数据点被当作是真正的中心。然后,该算法通过更新选择的中心位置进行迭代,数据科学家称之为中心点。在一次迭代中,每个数据点都被分配到其最近的中心,从而形成K组。接下来,每个组的中心被更新。新的中心等于该组坐标的平均值。如果我们重复这个过程的时间足够长,组的平均值将收敛到K个代表中心(图6)。这种收敛在数学上是有保证的。然而,我们不能事先知道收敛所需的迭代次数。一个常见的技巧是,当新计算的中心都没有明显偏离其前身时,就停止迭代。
图6.K-means算法从两个随机选择的中心点迭代收敛到实际的牛眼中心点
K-means并非没有局限性。该算法的前提是我们对K的了解:要寻找的聚类数量。通常,这种知识是不可用的。另外,虽然K-means通常能找到合理的中心,但它在数学上不能保证能找到数据中可能的最佳中心。偶尔,由于在算法的初始化步骤中对随机中心点的选择不当,K-means会返回非直观的或次优的组。最后,K-means的前提是,数据中的群组实际上是围绕着K中心位置旋转的。但正如我们在本节后面所了解的,这个假设并不总是成立的。
使用scikit-learn进行K-means聚类
如果K-means算法被有效地实现,它就能在合理的时间内运行。该算法的快速实现可以通过外部的scikit-learn库来实现。Scikit-learn是一个非常流行的机器学习工具包,建立在NumPy和SciPy之上。它具有各种核心分类、回归和聚类算法--当然也包括K-means。让我们来安装这个库。然后我们导入scikit-learn的KMeans 聚类算法。
从命令行终端调用pip install scikit-learn 来安装scikit-learn库。
清单8.从scikit-learn导入KMeans
from sklearn.cluster import Kmeans
将KMeans 应用于我们的darts 数据很容易。首先,我们需要运行KMeans(n_clusters=2) ,这将创建一个能够找到两个牛眼中心的cluster_model 对象。然后我们可以通过运行cluster_model.fit_predict(darts) 来执行K-means。该方法调用将返回一个assigned_bulls_eyes 数组,该数组存储每个飞镖的牛眼索引。
清单9.使用scikit-learn进行K-means聚类
cluster_model = KMeans(n_clusters=2)
❶ 创建一个 cluster_model 对象,其中的中心数被设置为 2
❷ 使用K-means算法优化两个中心,并返回为每个飞镖分配的簇。
让我们根据它们的聚类分配给我们的飞镖上色,以验证结果(图7)。
清单10.绘制K-means聚类分配图
for bs_index in range(len(bulls_eyes)):
selected_darts = [darts[i] for i in range(len(darts))
if bs_index == assigned_bulls_eyes[i]]
x_coordinates, y_coordinates = np.array(selected_darts).T
plt.scatter(x_coordinates, y_coordinates,
color=['g', 'k'][bs_index])
plt.show()
图7.scikit-learn返回的K-means聚类结果与我们的预期一致。
我们的聚类模型已经在数据中找到了中心点。现在我们可以重新使用这些中心点来分析模型以前没有见过的新数据点。执行cluster_model.predict([x, y]) ,将一个中心点分配给由x 和y 定义的数据点。我们使用predict 方法来聚类两个新的数据点。
清单11.使用cluster_model 对新数据进行聚类
new_darts = [[500, 500], [-500, -500]]
new_bulls_eye_assignments = cluster_model.predict(new_darts)
for i, dart in enumerate(new_darts):
bulls_eye_index = new_bulls_eye_assignments[i]
print(f"Dart at {dart} is closest to bull's-eye {bulls_eye_index}")
Dart at [500, 500] is closest to bull's-eye 0
Dart at [-500, -500] is closest to bull's-eye 1
使用肘部方法选择最佳K
K-means依赖于一个输入的K。当数据中的真实聚类的数量事先不知道时,这可能是一个严重的障碍。然而,我们可以使用一种被称为肘部方法的技术来估计一个合适的K值。
肘部方法依赖于一个被称为惯性的计算值,它是每个点与其最接近的K-均值中心之间的平方距离之和。如果K是1,那么惯性就等于所有与数据集平均值的平方距离之和。正如第5节所讨论的,这个值与方差成正比。而方差又是对分散性的一种衡量。因此,如果K是1,惯性是对分散性的估计。基本上,惯性估计我们的K计算的平均值周围的总离散度。
通过估计分散性,我们可以确定我们的K值是太高还是太低。例如,想象一下,我们把K值设置为1。有可能,我们的许多数据点将被定位在离一个中心太远的地方。我们的分散性会很大,我们的惯性也会很大。当我们把K增加到一个更合理的数字时,额外的中心将导致惯性的减少。最终,如果我们过分地将K设置为等于总点数,每个数据点都将落入自己的私人集群。分散性将被消除,而惯性将下降到零(图8)。
图8.六个编号为1到6的点被绘制在二维空间。用星星标记的中心是在不同的K值下计算出来的。从每一个点到它最近的中心都画一条线。惯性是通过六条线的平方长度之和计算出来的。(A)K=1。所有六条线都从一个中心延伸出去。惯性是相当大的。(B)K= 2。5点和6点非常接近第二个中心。惯性减小。©K= 3.点1和点3在很大程度上更接近一个新形成的中心。点2和点4也大大接近了一个新形成的中心。惯性从根本上减少了。(D)K= 4。点1和点3现在与它们的中心重合。它们对惯性的贡献已经从一个非常低的值转变为零。剩下的四个点和它们相关的中心之间的距离保持不变。因此,将K从3增加到4引起了非常小的惯性减少。
有些惯性值太大。另一些则太低了。在这两者之间可能有一个恰到好处的值。我们怎样才能找到它呢?
让我们想出一个解决方案。我们首先绘制出我们的镖靶数据集在很大的K值范围内的惯性(图9)。惯性是为每个scikit-learnKMeans 对象自动计算的。我们可以通过模型的_inertia 属性访问这个存储值。
清单12.绘制K-means的惯性图
k_values = range(1, 10)
inertia_values = [KMeans(k).fit(darts).inertia_ for k in k_values]
plt.plot(k_values, inertia_values)
plt.xlabel('K')
plt.ylabel('Inertia')
plt.show()
图9.一个包含两个牛眼目标的镖靶模拟的惯性图。该图类似于一个胳膊在肘部弯曲的情况。肘部直接指向K为2。
我们已经知道,这个K值准确地抓住了我们预先编入数据集的两个中心。
如果增加现在的中心数量,这个方法还能成立吗?我们可以通过在我们的投掷飞镖模拟中增加一个额外的牛眼来了解一下。在我们将群集数量增加到三个之后,我们重新生成了我们的惯性图(图10)。
清单13.绘制3个镖靶模拟的惯性图
new_bulls_eye = [12, 0]
for _ in range(5000):
x = np.random.normal(new_bulls_eye[0], variance ** 0.5)
y = np.random.normal(new_bulls_eye[1], variance ** 0.5)
darts.append([x, y])
inertia_values = [KMeans(k).fit(darts).inertia_
for k in k_values]
plt.plot(k_values, inertia_values)
plt.xlabel('K')
plt.ylabel('Inertia')
plt.show()
图10.一个包含三个牛眼目标的镖靶模拟的惯性图。该图类似于一个胳膊弯曲的肘部。肘部的最下面部分指向3的K。
从本质上讲,我们的肘部图追踪了每个递增的K值所捕获的分散性。惯性的减少随着惯性曲线的变平而逐渐失去其影响。这种从垂直下降到较平缓角度的过渡导致了我们的图中出现了一个肘部形状。我们可以利用肘部的位置在K-means算法中选择一个合适的K。
肘部方法选择标准是一个有用的启发式方法,但它不能保证在每一种情况下都有效。在某些条件下,肘部在多个K值上慢慢趋于平缓,使得它很难选择一个有效的聚类计数。
存在着更强大的K选择方法,比如剪影分数,它可以捕捉到每个点与相邻聚类的距离。对剪影分数的彻底讨论超出了本书的范围。然而,我们鼓励你使用sklearn.metrics.silhouette_score 方法,自己去探索这个分数。
K-means聚类方法
k_means_model = KMeans(n_clusters=K)- 创建一个K-means模型来搜索K个不同的中心点。我们需要将这些中心点与输入的数据相适应。clusters = k_means_model.fit_predict(data)- 使用一个初始化的KMeans对象对输入的数据执行K-means。返回的clusters数组包含从0到K的簇ID。data[i]的簇ID等于clusters[i]。clusters = KMeans(n_clusters=K).fit_predict(data)- 在一行代码中执行K-means,并返回产生的聚类。new_clusters = k_means_model.predict(new_data)- 使用数据优化的KMeans对象中的现有中心点,找到与以前未见过的数据最近的中心点。inertia = k_means_model.inertia_- 返回与数据优化的KMeans对象相关的惯性。inertia = KMeans(n_clusters=K).fit(data).inertia_- 在一行代码中执行K-means,并返回产生的惯性。
肘部方法并不完美,但是如果数据是以K个不同的手段为中心,它的表现还算不错。当然,这是假设我们的数据集群由于中心性而不同。然而,在许多情况下,数据集群的不同是由于数据点在空间中的密度。让我们来探讨一下密度驱动的集群的概念,它不依赖于中心性。
使用密度来发现集群
假设一位天文学家在太阳系的遥远边缘发现了一颗新的行星。这颗行星,很像土星,有多个环围绕其中心不断旋转。每个环都是由成千上万的岩石形成的。我们将把这些岩石建模为由X和Y坐标定义的单个点。让我们使用scikit-learn的makes_circles (图11),生成由许多岩石组成的三个岩环。
清单14.模拟一个行星周围的环
from sklearn.datasets import make_circles
x_coordinates = []
y_coordinates = []
for factor in [.3, .6, 0.99]:
rock_ring, _ = make_circles(n_samples=800, factor=factor,
❶ make_circles函数在二维空间创建两个同心圆。小圆的半径相对于大圆的比例由因子参数决定。
图11.围绕一个中心点定位的三个岩环的模拟
图中清楚地出现了三个岩环群。让我们用K-means搜索这三个群组,将K设置为3(图12)。
清单15。使用K-means对岩环进行聚类
rocks = [[x_coordinates[i], y_coordinates[i]]
for i in range(len(x_coordinates))]
rock_clusters = KMeans(3).fit_predict(rocks)
colors = [['g', 'y', 'k'][cluster] for cluster in rock_clusters]
plt.scatter(x_coordinates, y_coordinates, color=colors)
plt.show()
图12.K-means聚类未能正确识别三个不同的岩石环。
其输出结果是完全失败的!K-means将数据剖析成三个对称的片段,而每个片段都跨越了多个环。这个解决方案与我们的直觉预期不一致,即每个环都应该属于自己的不同组。什么地方出错了?好吧,K-means假设这三个集群是由三个独特的中心定义的,但实际的环是围绕着一个中心点旋转的。集群之间的差异不是由中心性驱动的,而是由密度驱动的。每个环都是由密集的点组成的,环与环之间的空白区域是由稀疏的空间作为边界的。
我们需要设计一种算法,将数据聚集在空间的密集区域。这样做需要我们定义一个给定的区域是密集的还是稀疏的。一个简单的密度定义如下:只有当一个点位于其他Y个点的X距离内时,它才是在一个密集区域。我们将X和Y分别称为epsilon 和min_points 。下面的代码将epsilon 设为0.1,min_points 设为10。因此,如果我们的岩石位于至少10个其他岩石的0.1半径范围内,那么它们就存在于一个密集的空间区域。
清单16.指定密度参数
epsilon = 0.1
min_points = 10
让我们分析一下我们的rocks 列表中第一块岩石的密度。我们首先搜索epsilon 单位内的所有其他岩石,rocks[0] 。我们将这些相邻岩石的索引存储在一个neighbor_indices 列表中。
清单17.找出 "A "的邻居rocks[0]
neighbor_indices = [i for i, rock in enumerate(rocks[1:])
if euclidean(rocks[0], rock) <= epsilon]
现在我们比较min_points 的邻居数量,以确定rocks[0] 是否位于空间的密集区域。
清单18.检查rocks[0]
num_neighbors = len(neighbor_indices)
print(f"The rock at index 0 has {num_neighbors} neighbors.")
if num_neighbors >= min_points:
print("It lies in a dense region.")
else:
print("It does not lie in a dense region.")
The rock at index 0 has 40 neighbors.
It lies in a dense region.
索引为0的岩石位于空间的密集区。rocks[0] 的邻居是否也共享这个空间的密集区域?这是个很棘手的问题。毕竟,有可能每个邻居都有少于min_points 的邻居。根据我们严格的密度定义,我们不会认为这些邻居是密集点。然而,这将导致一种可笑的情况,即密集区只由一个点组成:rocks[0] 。我们可以通过更新我们的密度定义来避免这种荒谬的结果。让我们正式定义密度如下。
min_point如果一个点位于epsilon邻居的距离之内,那么这个点就在空间的密集区。- 一个点在空间密集区中的每一个邻居也都在该空间中聚类。
基于我们更新的定义,我们可以将rocks[0] 和它的邻居合并为一个单一的密集集群。
清单19.创建一个密集簇
dense_region_indices = [0] + neighbor_indices
dense_region_cluster = [rocks[i] for i in dense_region_indices]
dense_cluster_size = len(dense_region_cluster)
print(f"We found a dense cluster containing {dense_cluster_size} rocks")
We found a dense cluster containing 41 rocks
索引为0的岩石和它的邻居形成了一个单一的41元素的密集簇。邻居的任何一个邻居是否属于空间的密集区域?如果是,那么根据我们更新的定义,这些岩石也属于密集簇。因此,通过分析额外的邻接点,我们可以扩大dense_region_cluster 的规模。
清单20:扩展一个密集簇
dense_region_indices = set(dense_region_indices)
❶ 将dense_region_indices转换为一个集合。这使得我们可以用额外的指数来更新这个集合,而不用担心重复的问题。
我们已经迭代了邻居的邻居,并将我们的密集簇扩大了近20倍。为什么停在这里呢?我们可以通过分析新遇到的邻居的密度来进一步扩大我们的集群。迭代重复我们的分析将增加我们集群边界的广度。最终,这个边界将扩展到完全包括我们的一个岩环。然后,由于没有新的邻居可以吸收,我们可以在一个迄今为止尚未分析的rocks 元素上重复迭代分析。这种重复将导致更多致密环的聚类。
刚刚描述的程序被称为DBSCAN。DBSCAN算法根据数据的空间分布来组织数据。