Python 无监督学习实用指南：1~5（3）

Python 无监督学习实用指南：1~5（2）https://developer.aliyun.com/article/1426881

权变矩阵

一个非常简单而强大的工具，可以在已知真实情况时显示聚类算法的表现，它是权变矩阵C[m]。 如果存在m类，则C[m] ∈ ℜ^(m×m)和每个元素C[m](i, j)代表已分配给群集j的Y[true] = i的样本数。 因此，一个完美的权变矩阵是对角线的，而所有其他单元格中元素的存在则表明了聚类误差。

在我们的案例中，我们获得以下信息：

from sklearn.metrics.cluster import contingency_matrix
cm = contingency_matrix(kmdff['diagnosis'].apply(lambda x: 0 if x == 'B' else 1), kmdff['prediction'])

上一个代码片段的输出可以显示为热图（变量cm是2×2矩阵）：

权变矩阵的图形表示

该结果表明，几乎所有良性样本均已正确聚类，而适度百分比的恶性样本已被错误地分配给第一个聚类。 我们已经使用其他度量进行了确认，但是类似于分类任务中的混淆矩阵，列联矩阵可以立即可视化最难分离的类别，从而帮助数据科学家寻找更有效的解决方案。

K 最近邻

K 最近邻（KNN）是属于称为基于实例的学习类别的方法。 在这种情况下，没有参数化模型，而是样本的重新排列以加快特定查询的速度。 在最简单的情况下（也称为暴力搜索），假设我们有一个X数据集，其中包含M个样本x[i] ∈ ℜ^N。 给定距离函数d(x[i], x[j])，则可以定义测试样本的半径邻域x[i]如下：

集合ν(x[i])是一个以x[i]为中心的球，包括所有距离小于或等于的样本R。 另外，也可以只计算最接近的k个邻居，即更接近x[i]的k个样本（通常， 该集合是ν(x[i])的子集，但当k非常大）。 该过程很简单，但不幸的是，从计算的角度来看太昂贵了。 实际上，对于每个查询，有必要计算M^2个N维距离（即，假设每距离N个运算） ，复杂度为O(NM^2)，这是使暴力破解方法遭受维度诅咒的条件。 例如，在N = 2和M = 1,000 的情况下，复杂度为O(2 * 10^6)，但是当N = 1,000和M = 10,000时，其变为O(10^11)。 例如，如果每个操作需要 1 纳秒，那么查询将需要 100 秒，这在许多实际情况下超出了可容忍的限制。 此外，对于 64 位浮点值，成对距离矩阵每次计算将需要约 764MB，再次考虑到任务的性质，这可能是一个过多的请求。

由于这些原因，仅当M非常小时并且在所有其他情况下都依赖于稍微复杂的结构时，KNN 具体实现才使用暴力搜索。 第一种替代方法基于 kd 树，这是将二叉树自然扩展到多维数据集的方法。

在下图中，表示了由3维向量组成的部分 kd 树：

具有 3 维向量的 kd 树示例

kd 树的构建非常简单。 给定一个根样本(a[1], a[2], ..., a[n])，考虑第一个特征，因此左分支包含b[1] < a[1]，以此类推，和右分支c[1] > a[1]，以此类推。 该过程将继续执行第二个特征，第三个特征，依此类推，直到到达叶节点为止（分配给叶的样本数量是需要调整的超参数。 该参数称为leaf_size，默认值为 30 个样本）。

当维度N不太大时，计算复杂度变为O(N log M)，这比暴力搜索要好得多。 例如，在N = 1,000和M = 10,000的情况下，计算复杂度变为O(4,000) << O(10^11)。 不幸的是，当N大时，kd 树查询变为O(NM)，因此，考虑前面的示例，O(10^7)，它比蛮横搜索更好，但有时对于实时查询而言仍然太昂贵。

KNN 中常用的第二个数据结构是球树。 在这种情况下，根节点由R[0]球表示，精确定义为样本的邻域：

选择第一个球以便捕获所有样本。 此时，将其他较小的球嵌套到β[R0]中，以确保每个样本始终属于一个球。 在下图中，有一个简单的球树的示意图：

一个简单的球树的例子

由于每个球都由其中心c[j]完全确定，因此对测试样本x[i]的查询要求计算距离d(x[i], c[j])。 因此，从底部（最小的球所在的位置）开始，执行完整的扫描。 如果没有一个球包含样本，则级别会增加，直到达到根节点为止（记住一个样本可以属于一个球）。 由于球的特性，计算复杂度现在始终为O(N log M)（也就是说，给定中心和半径，可以通过一次距离计算来检查样本的隶属度） 。 确定正确的球后，样本x[i]的邻居需要计算有限数量的成对距离（该值小于叶大小，因此与数据集的维数相比通常可以忽略不计）。

当然，这些结构是在训练阶段构建的，在生产阶段不会对其进行修改。 这意味着要仔细选择最小半径或分配给叶节点的样本数。 实际上，由于查询通常需要多个邻居k，因此仅当k < |ν(x[i])|达到最佳值时，才能实现最优。* 。 换句话说，我们想在同一子结构中找到所有包含x[i]的邻居。 每当k > |ν(x[i])|，该算法还必须检查相邻结构并合并结果。 当然，当叶子大小太大（与样本总数M相比）时，这些树的优势就消失了，因为有必要计算太多的成对距离才能回答查询。 必须根据软件的生产使用情况来正确选择叶子大小。

例如，如果推荐系统需要具有 100 个邻居的初始查询和具有 10 个邻居的几个（例如 5 个）后续查询，则等于 10 的叶子大小将优化优化阶段，但在第一个查询上会产生负面影响。 相反，选择等于 100 的叶子大小将减慢所有 10 个邻居查询的速度。 权衡可能是 25，这减少了第一个查询的负担，但对细化查询的成对距离的计算产生了中等程度的负面影响。

现在，我们可以基于 Olivetti 人脸数据集（由 scikit-learn 直接提供）分析一个简短示例。 它由代表不同人物肖像的 400 张64×64灰度图像组成。 让我们从如下加载数据集开始：

from sklearn.datasets import fetch_olivetti_faces
faces = fetch_olivetti_faces()
X = faces['data']

变量X包含数据集的展开版本（400 个 4,096 维实例已经在 0 和 1 之间标准化）。 在这一点上，我们可以训练一个NearestNeighbor模型，假设使用 10 个样本（参数n_neighbors）和半径等于 20（参数radius）的默认查询。 我们保留默认的leaf_size (30)，并使用p=2（欧几里得距离）明确设置 Minkowski 度量。 该算法基于一棵球树，但我邀请读者同时测试不同的指标和 kd 树。 现在，我们可以创建NearestNeighbors实例并继续训练模型：

from sklearn.neighbors import NearestNeighbors
knn = NearestNeighbors(n_neighbors=10, metric='minkowski', p=2, radius=20.0, algorithm='ball_tree')
knn.fit(X)

训练好模型后，请使用嘈杂的测试脸来查找最近的 10 个邻居，如下所示：

import numpy as np
i = 20
test_face = X[i] + np.random.normal(0.0, 0.1, size=(X[0].shape[0]))

测试面绘制在以下屏幕截图中：

嘈杂的测试面

可以使用仅提供测试样本的方法kneighbors()来执行具有默认邻居数的查询（在邻居数不同的情况下，必须调用该函数，并同时提供参数n_neighbors）。 如果参数为return_distance=True，则该函数返回包含distances, neighbors的元组，如下所示：

distances, neighbors = knn.kneighbors(test_face.reshape(1, -1))

查询结果显示在以下屏幕截图中：

测试样本的最近邻居及其相对距离

第一个样本始终是测试样本（在这种情况下，它被去噪，因此其距离不为零）。 可以看到，即使距离是一个累积函数，第二个和第四个样本是指同一个人，而其他样本共享不同的解剖元素。 当然，欧几里德距离并不是衡量图像之间差异的最合适方法，但是该示例在一定程度上证实了当图像相当相似时，全局距离也可以为我们提供用于查找相似的样本有价值的工具。

现在，使用radius_neighbors()设置radius=100的方法执行半径查询，如下所示：

import numpy as np
distances, neighbors = knn.radius_neighbors(test_face.reshape(1, -1), radius=100.0)
sd, sd_arg = np.sort(distances[0]), np.argsort(distances[0])

以下屏幕快照显示了包含前 20 个邻居的结果：

使用半径查询的前 50 个邻居

有趣的是，距离并没有很快变化（第二个样本具有d=8.91，第五个d=10.26）。 这主要是由于两个因素：第一个是样本之间的全局相似性（就几何元素和色调而言），第二个可能与欧氏距离对 4,096 维向量的影响有关。 正如在谈到聚类基本原理时所解释的那样，高维样本可能缺乏可区分性（尤其是当p >> 1时）。 在这种情况下，图片不同部分的平均效果可能会产生与分类系统完全不兼容的结果。 特别是，深度学习模型倾向于通过使用可以学会检测不同级别特定特征的卷积网络来避免此陷阱。 我建议以不同的指标重复此示例，并观察p对半径查询样本所显示的实际差异的影响。

向量量化

向量量化（VQ）是一种利用无监督学习对样本x[i] ∈ ℜ^N或整个数据集X进行有损压缩的方法。（为简单起见，我们假设多维样本被展开）。 主要思想是找到带有许多条目C << N的密码本Q，并将每个元素与条目q[i] ∈ Q相关联。 在单个样本的情况下，每个条目将代表一个或多个特征组（例如，可以是均值），因此，该过程可以描述为一种变换T，其一般表示为 ：

码本被定义为Q = (q[1], q[2], ..., q[C])。 因此，给定一个由一组特征集合（例如，一组两个连续元素）组成的综合数据集，VQ 将关联一个码本条目：

由于使用汇总整个组的固定值的组合表示输入样本，因此将过程定义为量化。 类似地，如果输入是数据集X，则转换将按样本组进行操作，就像任何标准聚类过程一样。 主要区别在于目的：使用 VQ 代表每个聚类及其质心，从而减少数据集的方差。 这个过程是不可逆的。 一旦执行了转换，就不可能重建原始聚类（唯一可行的过程是基于具有相同原始均值和协方差的分布的采样，但是重建显然是近似的）。

让我们从显示一个非常简单的高斯数据集的示例开始，如下所示：

import numpy as np
nb_samples = 1000
data = np.random.normal(0.0, 1.5, size=(nb_samples, 2))
n_vectors = 16
qv = np.random.normal(0.0, 1.5, size=(n_vectors, 2))

我们的目标是用 16 个向量表示数据集。 在以下屏幕截图中，该图显示了初始配置的示例：

VQ 示例的向量的初始配置

当我们使用随机数时，相同代码的后续执行会产生不同的初始配置。 该过程遍历所有样本，选择最接近的量化向量，并将其距离减小固定量delta=0.05，如下所示：

import numpy as np
from scipy.spatial.distance import cdist
delta = 0.05
n_iterations = 1000
for i in range(n_iterations):
    for p in data:
        distances = cdist(qv, np.expand_dims(p, axis=0))
        qvi = np.argmin(distances)
        alpha = p - qv[qvi]
        qv[qvi] += (delta * alpha)
distances = cdist(data, qv)
Y_qv = np.argmin(distances, axis=1)

除了固定的for循环外，还可以使用while循环来检查量化向量是否已达到稳态（比较t和t + 1）。 以下屏幕快照显示了该过程结束时的结果：

量化向量的最终配置（左）。 每个量化向量的影响范围（右）

正如预期的那样，量化向量已达到最终配置，其中每个量化向量都代表数据集的一小部分（如右图所示）。 在这一点上，给定一个点，最接近的向量将代表它。 有趣的是，全局方差并未受到影响，但是，选择任何子集后，内部方差会大大降低。 向量的相对位置反映了数据集的密度，因为区域中的更多样本吸引了更多向量。 这样，通过构建距离矩阵，可以获得粗略的密度估计（例如，当向量与向量的近邻的平均距离较高时，意味着底层区域的密度较小）。 我们将在第 6 章，“异常检测”中更详细地讨论此主题。

现在让我们考虑一个示例，该示例具有一个代表浣熊图片的单个样本。 由于过程可能很长，因此第一步是加载示例 RGB 图像（由 SciPy 提供）并将其大小调整为 192×256，如下所示：

from scipy.misc import face
from skimage.transform import resize
picture = resize(face(gray=False), output_shape=(192, 256), mode='reflect')

以下屏幕快照显示了原始图片（已在[0, 1]范围内标准化）：

VQ 示例的 RGB 图片样本

我们想用 24 个使用2×2正方形区域计算的向量执行 VQ（由包含2×2×3特征的展开向量表示）。 但是，我们将使用 K 均值算法来查找质心，而不是从头开始执行该过程。 第一步是收集所有正方形区域，如下所示：

import numpy as np
square_fragment_size = 2
n_fragments = int(picture.shape[0] * picture.shape[1] / (square_fragment_size**2))
fragments = np.zeros(shape=(n_fragments, square_fragment_size**2 * picture.shape[2]))
idx = 0
for i in range(0, picture.shape[0], square_fragment_size):
    for j in range(0, picture.shape[1], square_fragment_size):
        fragments[idx] = picture[i:i + square_fragment_size, 
                                 j:j + square_fragment_size, :].flatten()
        idx += 1

此时，可以使用 24 个量化向量进行 K 均值聚类，如下所示：

from sklearn.cluster import KMeans
n_qvectors = 24
km = KMeans(n_clusters=n_qvectors, random_state=1000)
km.fit(fragments)
qvs = km.predict(fragments)

在训练结束时，变量qvs将包含与每个正方形区域关联的质心的索引（可通过实例变量cluster_centers_获得）。

现在可以使用质心构建量化的图像，如下所示：

import numpy as np
qv_picture = np.zeros(shape=(192, 256, 3))
idx = 0   for i in range(0, 192, square_fragment_size):
    for j in range(0, 256, square_fragment_size):
        qv_picture[i:i + square_fragment_size,
                   j:j + square_fragment_size, :] = \
            km.cluster_centers_[qvs[idx]].\
                reshape((square_fragment_size, square_fragment_size, 3))   idx += 1

量化的图像显示在以下屏幕截图中：

用 24 个向量量化的图片

结果显然是原始图像的有损压缩版本。 每个组都可以用一个索引来表示（例如，在我们的示例中，它可以是 8 位整数），该索引指向码本中的条目（km.cluster_centers_）。 因此，如果最初有192×256×3 = 1,474,560个 8 位值，则在量化之后，我们有 12,288 个 8 位索引（2×2×3块的数量），再加上 24 个 12 维量化向量。 为了了解 VQ 对图像的影响，绘制原始图像和处理后图像的 RGB 直方图非常有用，如以下直方图所示：

原始图像的 RGB 直方图（顶部）和量化的版本（底部）

对于不熟悉直方图的读者，我们可以简单地将其描述为具有 X 数据集和固定数量的桶。 每个单元分配一个范围（从min(X)开始，以max(X)结束），并且每个范围(a, b)与样本数量相关，从而a ≤ x < b 。 结果图与生成X的实际概率分布的近似成比例。 在我们的情况下，在 x 轴上，每个通道（8 位）的每个像素都有所有可能的值，而 y 轴表示估计的频率（Nx/像素总数）。

可以看到，量化减少了信息量，但是直方图往往会重现原始信息。 增加量化向量的数量具有减少近似值的效果，从而产生差异较小的直方图。 对该主题的完整分析超出了本书的范围。 但是，我邀请读者使用其他图像和不同数量的量化向量来测试该过程。 也可以将原始图像的（协）方差（或熵）与量化版本进行比较，并找到保留 80% 的方差的阈值。 例如，仅考虑红色通道，并使用频率计数来近似每个值（0÷255）的概率，我们可以获得以下信息：

import numpy as np
hist_original, _ = np.histogram(picture[:, :, 0].flatten() * 255.0, bins=256)
hist_q, _ = np.histogram(qv_picture[:, :, 0].flatten() * 255.0, bins=256)
p_original = hist_original / np.sum(hist_original)
H_original = -np.sum(p_original * np.log2(p_original + 1e-8))
p_q = hist_q / np.sum(hist_q)
H_q = -np.sum(p_q * np.log2(p_q + 1e-8))
print('Original entropy: {0:.3f} bits - Quantized entropy: {1:.3f} bits'.format(H_original, H_q))

上一个代码段的输出如下：

Original entropy: 7.726 bits - Quantized entropy: 5.752 bits

由于信息量与熵成正比，因此我们现在已经确认，24 个量化向量（具有2×2正方形块）能够解释红色通道原始熵的大约 74% （即使三个通道都不是）。 独立地，可以通过对三个熵求和来获得总熵的粗略近似。 该方法可以有效地用于在压缩强度和最终结果质量之间进行权衡。

总结

在本章中，我们从相似性的概念及其度量方法入手，解释了聚类分析的基本概念。 我们讨论了 K 均值算法及其优化的变体 KMeans++ ，并分析了乳腺癌威斯康星州数据集。 然后，我们讨论了最重要的评估指标（无论是否了解基本事实），并且了解了哪些因素会影响绩效。 接下来的两个主题是 KNN（一种非常著名的算法，可用于在给定查询向量的情况下查找最相似的样本），以及 VQ（一种利用聚类算法以查找样本的有损表示形式的技术）（例如， 图片）或数据集。

在下一章中，我们将介绍一些最重要的高级聚类算法，展示它们如何轻松解决非凸问题。

问题

1. 如果两个样本的 Minkowski 距离（p = 5）等于 10，那么您能说出它们的曼哈顿距离吗？
   
2. 对 K 均值的收敛速度产生负面影响的主要因素是数据集的维数。 它是否正确？
   
3. 可以积极影响 K 均值表现的最重要因素之一是聚类的凸度。 它是否正确？
   
4. 聚类应用的同质性得分等于 0.99。 这是什么意思？
   
5. 调整后的兰德得分等于 -0.5 是什么意思？
   
6. 考虑到前面的问题，不同数量的聚类能否产生更好的分数？
   
7. 基于 KNN 的应用平均每分钟需要 100 个 5-NN 基本查询。 每分钟执行 2 个 50-NN 查询（每个查询需要 4 秒，叶子大小为 25），并在紧接其后执行 2 秒的阻塞任务。 假设没有其他延迟，则每分钟叶子大小= 50 可以执行多少个基本查询？
   
8. 球形树结构不适合管理高维数据，因为它遭受了维数的诅咒。 它是否正确？
   
9. 获得了一个数据集，该数据集从 3 个二维高斯分布中采样了 1,000 个样本：N([-1.0, 0.0], diag[0.8, 0.2])，N([0.0, 5.0], diag[0.1, 0.1])和N([-0.8, 0.0], diag[0.6, 0.3])。 集群中最可能的数量是？
   
10. 可以使用 VQ 压缩文本文件吗（例如，构建具有 10,000 个单词的字典，该单词在[0.0, 1.0]范围内均匀映射，将文本拆分为标记，然后将其转换为浮点序列）？
    

进一步阅读

• On the Surprising Behavior of Distance Metrics in High Dimensional Space, Aggarwal C. C., Hinneburg A., Keim D. A., ICDT, 2001
• K-means++: The Advantages of Careful Seeding, Arthur D., Vassilvitskii S., Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, 2007
• Visualizing Data using t-SNE, van der Maaten L., Hinton G., Journal of Machine Learning Research 9, 2008
• Robust Linear Programming Discrimination of Two Linearly Inseparable Sets, Bennett K. P., Mangasarian O. L., Optimization Methods and Software 1, 1992
• Breast cancer diagnosis and prognosis via linear programming, Mangasarian O. L., Street W.N, Wolberg W. H., Operations Research, 43(4), pages 570-577, July-August 1995
• V-Measure: A conditional entropy-based external cluster evaluation measure, Rosenberg A., Hirschberg J., Proceedings of the 2007 Joint Conference on Empirical Methods in Natural Language Processing and Computational Natural Language Learning, 2007

三、高级聚类

在本章中，我们将继续探索可用于非凸任务的更复杂的聚类算法（即，例如，K 均值无法同时获得内聚和分离力。 几何形状）。 我们还将展示如何将基于密度的算法应用于复杂的数据集，以及如何正确选择超参数并根据所需结果评估表现。 这样，数据科学家就可以准备面对各种问题，排除价值较低的解决方案，而只关注最有前途的解决方案。

特别是，我们将讨论以下主题：

• 谱聚类
• 均值漂移
• 具有噪声的基于密度的应用空间聚类（DBSCAN）
• 其他评估指标：Calinski-Harabasz 指数和群集不稳定性
• K-Medoids
• 在线聚类（小批量 K 均值和使用层次结构的平衡迭代归约和聚类（BIRCH））

技术要求

本章中提供的代码要求：

• Python3.5+（强烈建议您使用 Anaconda 发行版）
• 库：

• SciPy 0.19+
• NumPy 1.10+
• Scikit-Learn 0.20+
• Pandas 0.22+
• Matplotlib 2.0+
• Seaborn 0.9+

数据集可以通过 UCI 获得。 可以从这里下载 CSV 文件，并且不需要任何预处理，除了添加将在加载阶段出现的列名。

Github 存储库中提供了示例。

谱聚类

可以管理非凸类的最常见算法家族之一是谱聚类。 主要思想是将数据集X投影到可以通过超球体捕获群集的空间（例如，使用 K-means）。 该结果可以通过不同的方式实现，但是由于该算法的目标是消除通用形状区域的凹面，因此第一步始终是X表示为图{V，E}，其中顶点V ≡ X，加权边通过参数w[ij] ≥ 0表示每对样本x[i], x[j] ∈ X的邻近性。 生成的图可以是完整的（完全连接的），也可以仅在某些样本对之间具有边（即，不存在的权重的权重设置为零）。 在下图中，有一个局部图的示例：

图的示例：点x[0]是连接到x[1]的唯一点

可以使用两种主要策略来确定权重w[ij]：KNN 和径向基函数（RBF）。 第一个基于上一章中讨论的相同算法。 考虑到邻居数k，数据集表示为球树或 kd 树，对于每个样本x[i]，计算集合kNN(x[i])。 此时，给定另一个样本x[j]，权重计算如下：

在这种情况下，该图不包含有关实际距离的任何信息，因此，考虑到 KNN 中使用的相同距离函数d(·)，最好表示w[ij]如下：

此方法简单且相当可靠，但是生成的图未完全连接。 通过采用如下定义的 RBF 可以轻松实现这种条件：

通过这种方式，将根据其距离自动对所有夫妇加权。 由于 RBF 是高斯曲线，因此当x[i] = x[j]时，它等于1，并且与平方距离成比例地减小d(x[i], x[j])（表示为差异的范数）。 参数γ确定半铃曲线的振幅（通常，默认值为γ = 1）。 当γ < 1时，振幅增加，反之亦然。 因此，γ < 1表示对距离的灵敏度较低，而γ > 1的情况下，RBF 下降得更快，如以下屏幕截图所示：

二维 RBF，作为x和 0 之间的距离的函数，γ = 0.1、1.0 和 5.0

当γ = 0.1时，x = 1（相对于 0.0）的权重约为 0.9。 对于γ = 1.0，该值约为 0.5；对于γ = 5.0，该值几乎为零。 因此，在调整频谱聚类模型时，考虑γ的不同值并选择产生最佳表现的值（例如，使用第 2 章，“聚类基础知识”）。 一旦创建了图，就可以使用对称亲和矩阵W = {w[ij]} 来表示。 对于 KNN，W通常比较稀疏，可以使用专门的库进行有效地存储和操作。 相反，对于 RBF，它始终是密集的，并且如果X ∈ R^(N×M)，则它需要存储N^2个值 。

不难证明到目前为止我们分析过的程序等同于将X分割为多个内聚区域。 实际上，让我们考虑例如具有通过 KNN 获得的亲和度矩阵的图G。 连接的组件C[i]是一个子图，其中每对顶点x[a]， x[b] ∈ C[i]通过属于C[i]的顶点的路径连接，并且没有连接任何顶点的边C[i]的顶点具有不属于C[i]的顶点。 换句话说，连接的组件是一个内聚子集C[i] G，它代表群集选择的最佳候选者。 在下图中，有一个从图中提取的连接组件的示例：

从图中提取的连通组件的示例

在原始空间中，点x[0]，x[2]和x[3]通过x[1]连接到x[n]，x[m]和x[q]。 这可以表示非常简单的非凸几何，例如半月形。 实际上，在这种情况下，凸度假设对于最佳分离不再是必需的，因为，正如我们将要看到的那样，这些分量被提取并投影到具有平坦几何形状的子空间中（可通过诸如此类的算法轻松管理） 以 K 均值表示）。

当使用 KNN 时，此过程更加明显，但是，通常，我们可以说，当区域间距离（例如，两个最近点之间的距离）与平均内部区域相当时，可以合并两个区域。 Shi 和 Malik（在《归一化剪切和图像分割》中）提出了解决此问题的最常用方法之一。这称为归一化剪切。 整个证明超出了本书的范围，但是我们可以讨论主要概念。 给定一个图，可以构建归一化的图拉普拉斯算子，定义为：

对角矩阵D被称为度矩阵，每个元素d[i][i]是相应行的权重之和，可以证明以下陈述：

• 在对L进行特征分解之后（考虑非正规图拉普拉斯L[u] = D-W并求解方程L[u]·v =λDv，很容易计算特征值和特征向量 ），则空特征值始终以多重性出现p。
• 如果G是无向图（因此wij] ≥ 0, ∀i, j），则连接的组件数等于p（空特征值的多重性）。
• 如果A包含于ℜ^N和Θ是A的可数子集（即X是可计数的子集，因为样本数始终是有限的），向量v∈ℜ^N对于Θ被称为指标向量，给定θ[i] ∈ Θ，如果θ[i] ∈ A，v^(i) = 1，否则v^(i) = 0。 例如，如果我们有两个向量a = (1, 0)和b = (0, 0)（因此，Θ = {a, b}），我们认为A = {(1, n)}其中n ∈ [1, 10]，向量v = (1, 0)是一个指标向量，因为a ∈ A和b ∉ A。
• L的第一个p特征向量（对应于空特征值）是每个连接的分量C[1], C[2], ..., C[p]。

因此，如果数据集由M个样本x[i] ∈ ℜ^N以及图G与亲和力矩阵W^(M×M)相关联，Shi 和 Malik 建议建立矩阵B ∈ ℜ^(M×p)包含第一个p特征向量作为列，并使用诸如 K 均值的更简单方法对行进行聚类。 实际上，每一行代表样本在p维子空间上的投影，其中非凸性由可以封装在规则球中的子区域表示。

现在，我们应用频谱聚类以分离由以下代码段生成的二维正弦数据集：

import numpy as np
nb_samples = 2000
X0 = np.expand_dims(np.linspace(-2 * np.pi, 2 * np.pi, nb_samples), axis=1)
Y0 = -2.0 - np.cos(2.0 * X0) + np.random.uniform(0.0, 2.0, size=(nb_samples, 1))
X1 = np.expand_dims(np.linspace(-2 * np.pi, 2 * np.pi, nb_samples), axis=1)
Y1 = 2.0 - np.cos(2.0 * X0) + np.random.uniform(0.0, 2.0, size=(nb_samples, 1))
data_0 = np.concatenate([X0, Y0], axis=1)
data_1 = np.concatenate([X1, Y1], axis=1)
data = np.concatenate([data_0, data_1], axis=0)

数据集显示在以下屏幕截图中：

谱聚类示例的正弦数据集

我们尚未指定任何基本事实； 但是，目标是将两个正弦曲线分开（非正弦曲线）。 很容易检查捕获正弦曲线的球是否还会包含许多属于另一个正弦曲线子集的样本。 为了显示纯 K 均值和频谱聚类之间的差异（scikit-learn 实现 Shi-Malik 算法，然后进行 K 均值聚类），我们将训练两个模型，后者使用 RBF（affinity参数，其中γ = 2.0（gamma参数）。 当然，我邀请读者也测试其他值和 KNN 相似性。 以下片段显示了基于 RBF 的解决方案：

from sklearn.cluster import SpectralClustering, KMeans
km = KMeans(n_clusters=2, random_state=1000)
sc = SpectralClustering(n_clusters=2, affinity='rbf', gamma=2.0, random_state=1000)
Y_pred_km = km.fit_predict(data)
Y_pred_sc = sc.fit_predict(data)

结果显示在以下屏幕截图中：

原始数据集（左）。 谱聚类结果（中心）。 K 均值结果（右）

如您所见，K 均值将数据集沿x轴划分为两个球，而谱聚类成功地正确分离了两个正弦曲线。 只要群集的数量和X的维数都不太大（在这种情况下，拉普拉斯算子的本征分解会变得非常昂贵），该算法就非常强大。 此外，由于该算法基于图切割程序，因此，当群集数为偶数时，它非常适合。

MeanShift

让我们考虑一个数据集X∈^(M×N)（MN维样本）是从多元数据生成过程p_data。 均值漂移算法应用于聚类问题的目的是找到p_data最大的区域，并将周围子区域中包含的样本关联到同一集群。 由于p_data是概率密度函数（PDF），因此将其表示为以一小部分参数（例如均值和方差）为特征的常规 PDF（例如，高斯）的总和。 以这种方式，可以认为 PDF 由样本生成的概率最高。 我们还将在第 5 章，“软聚类和高斯混合模型”，和第 6 章，“异常检测”中讨论该过程。 出于我们的目的，将问题重构为一个迭代过程，可以更新均值向量（质心）的位置，直到达到最大值为止。 当质心到达其最终位置时，将使用标准邻域函数将样本分配给每个聚类。

该算法的第一步是确定近似p_data的合适方法。 一种经典方法基于 Parzen 窗口（将在第六章，“异常检测”中进行讨论）。 就目前而言，可以说 Parzen 窗口是一个非负内核函数f(·)，其特征是称为带宽的参数（有关更多详细信息，请参阅原始论文《关于概率密度函数和模式的估计》）。 顾名思义，此参数的作用是加宽或限制 Parzen 窗口接近其最大值的区域。 考虑到与高斯分布的类比，带宽具有与方差相同的作用。 因此，较小的带宽将产生在均值附近非常峰值的函数，而较大的值与较平坦的函数关联。 不难理解，在这种特定情况下，群集的数量由带宽和相反的方式隐式确定。 因此，大多数实现（例如 scikit-learn）仅采用一个参数，然后计算另一个参数。 考虑到该算法已设计为适用于概率分布，自然的选择是指定所需带宽或让实现检测最佳带宽。 这个过程看起来比施加特定数量的群集更为复杂，但是，在许多实际情况下，尤其是当至少部分地了解了基本事实时，测试不同带宽的结果会更容易。

均值平移的最常见选择是用n个扁平核的总和来近似数据生成过程（n是形心数）：

因此，收敛之后，每个样本都由最接近的质心表示。 不幸的是，这种近似导致了分段函数，该分段函数不太可能代表实际过程。 因此，最好基于相同的基础内核采用平滑的 Parzen 窗口K(·)：

K(·)是平方距离（例如对于标准球）和带宽h的函数。 可以使用许多可能的候选函数，但是，当然最明显的是高斯核（RBF），其中h^2发挥方差的作用。 现在得到的近似值非常平滑，n峰对应于质心（即均值）。 定义函数后，就可以计算质心的最佳位置x[1], x[2], ..., x[n]。

给定质心和邻域函数（为简单起见，我们假定使用半径为h和K(x) ≠ 0, ∀x ∈ B[r]的标准球B[h]），相应的均值漂移向量定义为：

可以看到，m(·)是加权为K(·)的所有邻域样本的平均值。 显然，由于K(·)是对称的并且具有一定的距离，所以x[i]达到实际平均值。 带宽的作用是限制x[i]周围的区域。 现在应该更加清楚，较小的值会强制算法引入更多的质心，以便将所有样本分配给一个群集。 相反，较大的带宽可能导致单个群集具有最终配置。 迭代过程从初始质心猜测开始x[1]^(0), x[2]^(0), ..., x[n]^(0)等，并使用以下规则校正向量：

以前的公式很简单； 在每一步中，质心都会移动（移动）到m(·)附近。 这样，由于m(·)与相对于x[i]计算的邻域密度成比例。 ，当x[i]到达概率最大的位置m(·) → m[final]时，不需要更多更新 。 当然，收敛速度受样本数量的强烈影响。 对于非常大的数据集，该过程可能会变得很慢，因为每个均值漂移向量的计算都需要对邻域进行预先计算。 另一方面，当聚类标准由数据密度定义时，该算法非常有用。

作为示例，现在让我们考虑具有500二维样本的合成数据集，该样本由三个带有对角协方差矩阵的多元高斯生成，如下所示：

import numpy as np
nb_samples = 500
data_1 = np.random.multivariate_normal([-2.0, 0.0], np.diag([1.0, 0.5]), size=(nb_samples,))
data_2 = np.random.multivariate_normal([0.0, 2.0], np.diag([1.5, 1.5]), size=(nb_samples,))
data_3 = np.random.multivariate_normal([2.0, 0.0], np.diag([0.5, 1.0]), size=(nb_samples,))
data = np.concatenate([data_1, data_2, data_3], axis=0)

数据集显示在以下屏幕截图中：

MeanShift 算法示例的样本数据集

在这种情况下，我们知道基本事实，但是我们想测试不同的带宽并比较结果。 由于生成高斯粒子彼此非常接近，因此可以将某些外部区域识别为群集。 为了将研究重点放在最佳参数上，我们可以观察到平均方差（考虑不对称性）为 1，因此可以考虑值h = 0.9，1.0，1.2和 1.5。 此时，我们可以实例化 scikit-learn 类MeanShift，将h值通过参数bandwidth 传递，如下所示：

from sklearn.cluster import MeanShift
mss = []
Y_preds = []
bandwidths = [0.9, 1.0, 1.2, 1.5]
for b in bandwidths:
   ms = MeanShift(bandwidth=b)
    Y_preds.append(ms.fit_predict(data))
    mss.append(ms)

密度分析后，训练过程会自动选择质心的数量和初始位置。 不幸的是，这个数字通常大于最后一个数字（由于局部密度差异）； 因此，该算法将优化所有质心，但在完成操作之前，将执行合并过程以消除所有与其他质心太近（即重复的质心）的质心。 Scikit-learn 提供了参数bin_seeding，可以通过根据带宽对样本空间进行离散化（合并）来加快这项研究。 这样，有可能以合理的精度损失来减少候选者的数量。

下图显示了这四个训练过程结束时的结果：

MeanShift 在不同带宽下的聚类结果

如您所见，带宽的微小差异会导致群集数量不同。 在我们的情况下，最佳值为h=1.2，它产生的结果是确定了三个不同的区域（以及一个包含潜在异常值的额外聚类）。 最大聚类的质心大致对应于实际均值，但是聚类的形状与任何高斯分布都不相似。 这是可以通过采用其他方法解决的缺陷（在第 5 章 ， “软聚类和高斯混合模型”中进行了讨论）。 实际上，均值偏移适用于局部邻域，并且p_data不被认为属于特定分布。 因此，最终结果是将数据集非常准确地分割为高度密集的区域（注意不再需要最大分隔），这也可以从多个标准分布的叠加中得出。 没有任何先前的假设，我们不能期望结果非常规则，但是，将该算法与 VQ 进行比较，很容易注意到分配是基于找到每个密集 Blob 的最佳代表的想法。 因此，由高斯N(μ, Σ)产生的一些点以低概率分配给质心比μ更具代表性的另一个聚类。

DBSCAN

DBSCAN 是基于数据集密度估计的另一种聚类算法。 但是，与均值移位相反，没有直接参考数据生成过程。 在这种情况下，实际上，过程通过自下而上的分析建立了样本之间的关系，从X由高密度区域（气泡 ）由低密度的分隔。 因此，DBSCAN 不仅需要最大的分隔约束，而且为了确定群集的边界，它强制执行这种条件。 而且，此算法不允许指定所需的群集数，这是X结构的结果，但是，类似于均值移位，可以控制进程的粒度。

特别是，DBSCAN 基于两个基本参数：ε，它表示球B[ε](x[i])的半径，以样本x[i]为中心，和n[min]，这是B[ε](x[i])中必须包含的最小样本数，以便考虑x[i]作为核心点（即可以被视为集群的实际成员的点）。 形式上，给定一个函数N(·)，该函数对集合中包含的样本数进行计数，则在以下情况下将样本x[i] ∈ X称为核心点 ：

所有点x[j] ∈ B[ε](x[i])定义为从x[i]直接密度可达的。 这样的条件是点之间最强的关系，因为它们都属于以x[i]为中心的同一球。B[ε](x[i])包含的样本总数足够大，可以将邻域视为密集的子区域。 此外，如果存在x[i] → x[i + 1] → ... → x[j]，其中x[i + 1]是从x[i]直接密度可达的（适用于所有顺序对），x[j]定义为x[i]中密度可达的。 该概念非常直观，可以通过考虑以下图表立即理解：

如果n[min] = 4，则点x[2]从x[0]密度可达

如果我们将样本的最小数量设置为等于 4，则x[0]，x[1]和x[2]是核心点，x[1]从x[0]直接密度可达。 x[2]从x[1]直接密度可达。 因此，x[2]从x[0]密度可达。 换句话说，这意味着可以从x[0]开始并以x[2]结尾，定义一系列重叠的密集球N(·) ≥ n[min]。 通过添加更多定义，可以将此概念扩展到属于球的所有其他点：给定点x[k]，点x[i]和x[j]是密度连通的，如果x[i]和x[j]从x[k]密度可达。

重要的是要了解这种情况比密度可达性弱，因为为了保证密集链，有必要考虑第三点，该点代表两个点之间的连接器 密集的次区域。 实际上，可能有两个密度连接点a和b，而a不能从b达到密度。 （反之亦然）。 只要仅在一个方向上移动就满足最小数量的样本条件（即，属于一个球的样本不是均匀分布，而是倾向于以较小的超量累积），就会发生这种情况。

因此，例如，如果N(a) >> n[min]和N(a[1]) << N(a)，转换a → a[1]可以允许构建一个球B[ε](a)，其中还包含a[1]（以及许多其他要点）。 但是，在逆转换中a[1] →a，B[ε](a[1])的密度不足以建立直接的密度可达性条件。

因此，当在两个方向之一上移动时，较长的序列可能会被破坏，从而导密集度可达性丧失。 现在应该更清楚地知道，在x[i]和x[j]两点之间的密度连接可以使我们避免此问题，假设还有一个点可以同时达到x[i]和x[j]。

所有具有x[i]和x[j] ∈ X将分配给同一群集C[t]。 此外，如果x[k] ∈ C[t]，则所有密度可达的点x[p] ∈ X * x[k]也将属于同一群集。 从任何其他点x[i] ∈ X不可达到密度的点x[n]被定义为噪声点。 因此，与其他算法相反，DBSCAN 输出n群集以及一个包含所有噪声点的附加集（不必将其视为离群值，而应视为不属于任何密集子区域的点）。 当然，由于噪声点没有标签，因此其数目应相当低； 因此，重要的是要调整参数ε和n[min]，以达到双重目的：最大化内聚和分离度，避免过多的点被标记为噪点 。 没有实现此目标的标准规则，因此，我建议在做出最终决定之前测试不同的值。

最后，重要的是要记住，DBSCAN 可以处理非凸几何形状，并且与均值移位相反，它假设存在由低密度区域包围的高密度区域。 而且，它的复杂性与所采用的 KNN 方法（强力，球树或 kd 树）严格相关。 通常，当数据集不太大时，平均性能大约为O(N log N)，但可能趋于O(N^2)当N非常大时。 要记住的另一个重要元素是样本的大小。 正如我们已经讨论的那样，高维度量可以减少两点的可分辨性，从而对 KNN 方法的表现产生负面影响。 因此，当维数很高时，应避免（或至少仔细分析）DBSCAN，因为生成的群集不能有效地表示实际的密集区域。

在显示具体示例之前，最好先介绍一种在未知的真实情况下可以采用的进一步评估方法。

Calinski-Harabasz 分数

假设将聚类算法应用于包含M个样本的数据集X，以便将其分割为n[c]个群集，C[i]由重心μ[i]表示，i = 1..n。 我们可以将群集内分散度（WCD）定义如下：

如果x[i]是N维列向量，则X[k] ∈ ℜ^(N×N)。 不难理解，WCD(k)编码有关群集的伪方差的全局信息。 如果满足最大内聚条件，我们预计质心周围的分散性有限。 另一方面，即使WCD(k)也可能受到包含异常值的单个群集的负面影响。 因此，我们的目标是在每种情况下都将WCD(k)最小化。 以类似的方式，我们可以将集群间分散度（BCD）定义为：

在上一个公式中， N(C[i])是分配给群集C[i]的元素数和μ是整个数据集的全局质心。 考虑到最大分离的原理，我们希望有一个远离全局质心的密集区域。 BCD(k)精确地表达了这一原理，因此我们需要对其进行最大化以实现更好的表现。

Calinski-Harabasz 分数定义为：

由于不考虑质心计算是聚类算法的一部分，因此引入了对预测标签的显式依赖性。 分数没有绝对含义，但是有必要比较不同的值以了解哪种解决方案更好。 显然， CH[k](·)越高，聚类表现越好，因为这种条件意味着更大的分离度和更大的内部凝聚力。

使用 DBSCAN 分析旷工数据集

旷工数据集（按照本章开头的说明进行下载）由 740 条记录组成，其中包含有关请假几天的员工的信息。 共有 20 个属性，分别代表年龄，服务时间，教育程度，习惯，疾病，纪律衰竭，交通费用，从家到办公室的距离等（这些字段的完整说明，请参见这里）。 我们的目标是预处理数据并应用 DBSCAN，以发现具有特定语义内容的密集区域。

第一步是按以下方式加载 CSV 文件（必须更改占位符，以指向文件的实际位置）：

import pandas as pd
data_path = '<data_path>\Absenteeism_at_work.csv'
df = pd.read_csv(data_path, sep=';', header=0, index_col=0).fillna(0.0)
print(df.count())

上一条命令的输出如下：

Reason for absence                 740
Month of absence                   740
Day of the week                    740
Seasons                            740
Transportation expense             740
Distance from Residence to Work    740
Service time                       740
Age                                740
Work load Average/day              740
Hit target                         740
Disciplinary failure               740
Education                          740
Son                                740
Social drinker                     740
Social smoker                      740
Pet                                740
Weight                             740
Height                             740
Body mass index                    740
Absenteeism time in hours          740
dtype: int64

其中一些特征是分类的，并使用连续的整数进行编码（例如Reason for absence，Month of absence等）。 由于这些值可能会在没有确切语义原因的情况下影响距离（例如Month=12大于Month=10，但两个月的距离相等），因此在继续下一步之前，我们需要对所有这些特征进行单热编码 （新特征将添加到列表末尾）。 在以下代码段中，我们使用get_dummies() pandas 函数来执行编码； 然后删除原始列：

import pandas as pd
cdf = pd.get_dummies(df, columns=['Reason for absence', 'Month of absence', 'Day of the week', 'Seasons', 'Disciplinary failure', 'Education', 'Social drinker', 'Social smoker'])
cdf = cdf.drop(labels=['Reason for absence', 'Month of absence', 'Day of the week', 'Seasons', 'Disciplinary failure', 'Education', 'Social drinker', 'Social smoker']).astype(np.float64)

一键式编码的结果通常会在方法之间产生差异，因为许多特征将被限制为 0 或 1，而其他特征（例如，年龄）的范围可能会更大。 因此，最好对方法进行标准化（在不影响标准差的情况下，由于它们与现有信息内容成正比，因此保持不变是有帮助的）。 可以使用StandardScaler类设置参数with_std=False来完成此步骤，如下所示：

from sklearn.preprocessing import StandardScaler
ss = StandardScaler(with_std=False)
sdf = ss.fit_transform(cdf)

在这一点上，像往常一样，我们可以使用 t-SNE 算法来减少数据集的维数（使用n_components=2）并可视化结构。 数据框dff将包含原始数据集和 t-SNE 坐标，如下所示：

from sklearn.manifold import TSNE
tsne = TSNE(n_components=2, perplexity=15, random_state=1000)
data_tsne = tsne.fit_transform(sdf)
df_tsne = pd.DataFrame(data_tsne, columns=['x', 'y'], index=cdf.index) 
dff = pd.concat([cdf, df_tsne], axis=1)

生成的绘图显示在以下屏幕截图中：

旷工数据集的 t-SNE 二维表示

在进行任何考虑之前，重复一下 t-SNE 产生最佳的低维表示很重要，但是始终必须在原始数据集上测试算法，以检查 t-SNE 标识的邻居是否对应于实际的聚集体。 特别是，考虑到 DBSCAN 的结构，考虑到 t-SNE 表示形式，ε值可能是合理的，但是当移至更高维度的空间时，这些球不再能够捕获相同的样本。 但是，先前的图显示了被空白空间包围的密集区域的存在。 不幸的是，密度极不可能是均匀的（这是 DBSCAN 的建议要求之一，因为 ε和n[min]不能改变，但是在这种情况下，我们假设所有斑点的密度都是恒定的。

为了找到适合我们目的的最佳配置，我们绘制了群集数，噪声点数，轮廓分数和 Calinski-Harabasz 分数作为ε的函数，采用了p = 2，p = 4，p = 8和p = 12，如下图所示：

评估指标作为ε的函数

Silhouette 和 Calinski-Harabasz 均基于凸群集的假设（例如，色散显然是一种假设样本围绕质心呈放射状分布的度量），因此在非凸情况下其期望值通常较小 。 但是，我们要最大化两个分数（轮廓→ 1和 Calinski-Harabasz → ∞），同时避免大量聚类。 考虑到我们的最初目标（寻找以一组特定特征为特征的凝聚聚类），我们选择了ε = 25和 Minkowski 度量，其中p = 12， 这会产生合理数量的群集（13）和 22 个噪声点。 在第 2 章，“聚类基本原理”中，我们证明了，当p → ∞时（但效果对于p > 2已经可见），距离趋向于最大的特征差异。

因此，应该始终通过上下文分析来证明这种选择是合理的。 在这种情况下，我们可以假设每个（非）凸斑点代表一个由特定特征（具有所有其他特征的次要贡献）主导的类别，因此p = 12（导致 17 个群集） 对于中等粗粒度的分析（考虑有 20 个属性），这可能是一个很好的权衡。 此外，ε = 22.5与最高的 Calinski-Harabasz 得分之一 129.3 和轮廓得分约等于 0.2 关联。 特别是，后者的值表示总体聚类是合理正确的，但可能存在重叠。 由于基础几何很可能是非凸的，因此考虑到具有相应峰值的 Calinski-Harabasz 分数，这样的结果是可以接受的（通常在凸形场景中不是这样）。 较大的ε值会产生略高的轮廓分数（小于 0.23），但是所得的群集数和 Calinski-Harbasz 分数均不受所得构型的影响。 必须清楚这一选择尚未得到任何外部证据的证实，必须通过对结果的语义分析加以验证。 如果需要进行细粒度的分析，则可以使用具有更多群集和更多噪声点的配置（因此，读者可以玩转这些值并提供一个结果的解释）。 但是，此示例的最终目标仍然是相同的：分割数据集，以便每个群集包含特定的（可能是唯一的）属性。

现在，我们可以实例化DBSCAN模型，并使用包含规范化特征的数组sdf对其进行训练。 配置为ε = 25（参数eps）和n[min] = 3 （参数min_samples），以及 Minkowski 度量（metric='minkowski'）和p=12。

现在，我们可以执行以下集群：

from sklearn.cluster import DBSCAN
from sklearn.metrics import silhouette_score, calinski_harabaz_score
ds = DBSCAN(eps=25, min_samples=3, metric='minkowski', p=12)
Y_pred = ds.fit_predict(sdf)
print('Number of clusters: {}'.format(np.max(Y_pred) + 1))
print('Number of noise points: {}'.format(np.sum(Y_pred==-1)))
print('Silhouette score: {:.3f}'.format(silhouette_score(dff, Y_pred, metric='minkowski', p=12)))
print('Calinski-Harabaz score: {:.3f}'.format(calinski_harabaz_score(dff, Y_pred)))

由于DBSCAN用标签-1标记了噪声点，因此上一个代码段的输出如下：

Number of clusters: 13
Number of noise points: 22
Silhouette score: 0.2
Calinski-Harabaz score: 129.860

生成的绘图显示在以下屏幕截图中：

旷工数据集的聚类结果

如您所见（我建议运行代码以便获得更好的视觉确认），已成功检测出大多数孤立区域（即使在 t-SNE 图中没有内聚），并且已将样本分配给了相同群集。 我们还可以观察到两个基本结果：在 t-SNE 表示中，噪声点（带有叉号的标记）不是孤立的，并且某些群集被部分拆分。 这不是算法的失败，而是降维的直接结果。 在原始空间中，所有噪声点实际上都没有与任何其他样本紧密连接，但在 t-SNE 图中它们可能看起来重叠或接近某些斑点。 但是，我们对高密度和准粘结的非凸区域感兴趣，幸运的是，它们在二维图中也显示为连通。

现在让我们考虑两个不同的区域（为简单起见，将分析限制为单次热编码后的前 10 个属性）。 第一个是二维区域x < -45，如下所示：

sdff = dff[(dff.x < -45.0)]
print(sdff[sdff.columns[0:10]].describe())

以下屏幕截图显示了输出的打印精美版本：

对应于子数据集x < -45的统计度量

有两个因素可以立即引起我们的注意：运输费用（这似乎标准化为 179 的值）和子孙的数量（考虑到平均值和标准差，对于大多数样本而言，其为 0）。 我们还考虑服务时间和从住所到工作的距离，这可以帮助我们找到群集的语义标签。 所有其他参数的判别力都较小，因此在此简要分析中将它们排除在外。 因此，我们可以假设这样的子集群包含大约 40 岁的没有孩子的人，服务时间较长，居住在离办公室很远的地方（我请读者检查总体统计数据以确认这一点），并且交通费用标准化（ 例如一次性支出汽车）。

现在让我们将该结果与-20 < x < 20和y < 20的区域进行比较，如下：

sdff = dff[(dff.x > 20.0) & (dff.y > -20.0) & (dff.y < 20.0)]
print(sdff[sdff.columns[0:10]].describe())

相应的输出如下：

对应于子数据集-20 < x < -20和y < 20的统计度量

在这种情况下，运输费用会更大，而从住所到工作的距离大约是前一个示例的一半（还要考虑标准差）。 此外，平均儿子数为 1，雇员中有两个孩子的雇员比例适中，服务时间约为 12，标准差为 3.6。 我们可以推断出，该集群包含所有年龄在（28-58）之间有家庭的（已婚）人的所有样本，这些人有家庭，办公室相对较近，但旅行费用较高（例如，由于使用出租车服务）。 这样的员工倾向于避免加班，但是他们的平均工作量几乎与前面示例中观察到的相同。 即使没有正式的确认，我们也可以假设这样的员工通常效率更高，而第一批员工包含生产型员工，但是他们需要更多的时间来实现他们的目标（例如，由于旅行时间更长）。

这显然不是详尽的分析，也不是一组客观的陈述。 目的是展示如何通过观察样本的统计特征来找到聚类的语义内容。 在现实生活中，所有观察都必须由专家（例如，HR 经理）进行验证，以便了解分析的最后部分（特别是语义上下文的定义）是否正确或是否正确。 需要使用更多的群集，不同的指标或其他算法。 作为练习，我邀请读者分析包含单个群集的所有区域，以完成大图并测试与不同类别相对应的人工样本的预测（例如，非常小的年轻人，有三个孩子的雇员， 等等）。

作为表现指标的群集不稳定性

群集不稳定性是 Von Luxburg 提出的一种方法（在《群集稳定性：概述》中）可以用以下方法衡量算法的优缺点： 关于特定数据集。 它可以用于不同的目的（例如，调整超参数或找到最佳数目的群集），并且它相对容易计算。 该方法基于这样的想法，即满足最大内聚和分离要求的聚类结果也应该对数据集的噪声扰动具有鲁棒性。 换句话说，如果将数据集X分割为群集集C，则派生数据集X[n]（基于特征的细微扰动）应映射到同一群集集。 如果不满足此条件，则通常有两种可能性：噪声扰动太强或算法对小变化过于敏感，因此不稳定。 因此，我们定义了原始数据集X的一组k扰动（或二次采样）版本：

如果我们应用产生相同数目群集n[c]的算法A，我们可以定义A(X[i])和A(X[j])之间的距离度量d(·)，该值测量不一致分配的数量（即A(X[i])），并且可以表示为返回对应于每个点的赋值的向量函数，因此d(·)可以简单地计算不同标签的数量，假设算法（如果需要）以相同的方式播种，并且数据集显然没有被改组，则算法的不稳定性（对于k*X的噪声变化）定义为：

因此，不稳定性是几对噪声变化的聚类结果之间的平均距离。 当然，该值不是绝对的，因此可以得出的规则是：选择产生最小不稳定的配置。 同样重要的是，这种方法不能与之前讨论的其他方法相提并论，因为它是基于其他超参数（噪声变化的数量，噪声均值和方差，二次采样率等），因此可以产生不同的结果。 当A和X固定时。 特别是噪声的大小会极大地改变不稳定性，因此在确定高斯噪声的μ和Σ之前，有必要评估X的均值和协方差矩阵。在我们的示例中（基于“旷工”数据集中的 DBSCAN 聚类），我们从加性噪声项n[i] ~ N(E[X], Cov(X) / 4)开始创建了 20 个扰动版本。然后应用从均匀分布U(0, 1)中采样的乘法掩码。 这样，一些噪声项将被随机抵消或减少，如以下代码所示：

import numpy as np
data = sdf.copy()
n_perturbed = 20
n_data = []
data_mean = np.mean(data, axis=0)
data_cov = np.cov(data.T) / 4.0
for i in range(n_perturbed):
    gaussian_noise = np.random.multivariate_normal(data_mean, data_cov, size=(data.shape[0], ))
    noise = gaussian_noise * np.random.uniform(0.0, 1.0, size=(data.shape[0], data.shape[1]))
    n_data.append(data.copy() + noise)

在这种情况下，我们想将不稳定性计算为ε的函数，但是可以使用任何其他算法和超参数重复该示例。 此外，我们采用归一化的汉明距离，该距离与两个聚类结果之间不一致分配的数量成正比，如下所示：

from sklearn.cluster import DBSCAN
from sklearn.metrics.pairwise import pairwise_distances
instabilities = []
for eps in np.arange(5.0, 31.0, 1.5):
    Yn = []
    for nd in n_data:
        ds = DBSCAN(eps=eps, min_samples=3, metric='minkowski', p=12)
        Yn.append(ds.fit_predict(nd))
    distances = []
    for i in range(len(Yn)-1):
        for j in range(i, len(Yn)):
            d = pairwise_distances(Yn[i].reshape(-1, 1), Yn[j].reshape(-1, 1), 'hamming')
            distances.append(d[0, 0])
    instability = (2.0 * np.sum(distances)) / float(n_perturbed ** 2)
    instabilities.append(instability)

结果如下图所示：

应用于旷工数据集的 DBSCAN 的集群不稳定性，作为ε的函数

Python 无监督学习实用指南：1~5（4）https://developer.aliyun.com/article/1426884