Python 数据科学手册第二版（八）
原文：
zh.annas-archive.org/md5/051facaf2908ae8198253e3a14b09ec1译者：飞龙
第四十七章：深入理解 k-Means 聚类
在之前的章节中，我们探索了用于降维的无监督机器学习模型。现在我们将转向另一类无监督机器学习模型：聚类算法。聚类算法试图从数据的属性中学习出最优的分割或离散标记的群组点。
Scikit-Learn 和其他地方提供了许多聚类算法，但可能最容易理解的算法是称为 k-means 聚类 的算法，它在 sklearn.cluster.KMeans 中实现。
我们从标准导入开始：
In [1]: %matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
import numpy as np
介绍 k-Means
k-means 算法在未标记的多维数据集中搜索预定数量的簇。它使用简单的概念来定义最优的聚类：
-
簇中心 是属于该簇的所有点的算术平均值。
-
每个点更接近它自己的簇中心，而不是其他簇中心。
这两个假设是 k-means 模型的基础。我们很快将深入了解算法如何达到这个解决方案，但现在让我们看一看一个简单数据集，并查看 k-means 的结果。
首先，让我们生成一个包含四个不同斑点的二维数据集。为了突出这是一个无监督算法，我们将在可视化中省略标签（见图 47-1）。
In [2]: from sklearn.datasets import make_blobs
X, y_true = make_blobs(n_samples=300, centers=4,
cluster_std=0.60, random_state=0)
plt.scatter(X[:, 0], X[:, 1], s=50);
图 47-1：用于演示聚类的数据
肉眼看来，相对容易选出这四个簇。k-means 算法自动执行此操作，并在 Scikit-Learn 中使用典型的估计器 API：
In [3]: from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=4)
kmeans.fit(X)
y_kmeans = kmeans.predict(X)
让我们通过按照这些标签对数据进行着色来可视化结果（图 47-2）。我们还将绘制由 k-means 估计器确定的簇中心：
In [4]: plt.scatter(X[:, 0], X[:, 1], c=y_kmeans, s=50, cmap='viridis')
centers = kmeans.cluster_centers_
plt.scatter(centers[:, 0], centers[:, 1], c='black', s=200);
令人高兴的是，k-means 算法（至少在这个简单案例中）将点分配到簇中的方式与我们通过肉眼观察的方式非常相似。但你可能会想知道这个算法是如何如此快速地找到这些簇的：毕竟，簇分配的可能组合数随数据点数量呈指数增长——全面搜索将非常、非常昂贵。对我们来说幸运的是，这样的全面搜索并不是必需的：相反，k-means 的典型方法涉及一种直观的迭代方法，称为期望最大化。
图 47-2：带有颜色指示簇的 k-means 簇中心
期望最大化
期望最大化（E–M）是数据科学中多种情境中的一个强大算法。k-means 是该算法的一个特别简单且易于理解的应用；我们将在这里简要介绍它。简而言之，在这里期望最大化方法包括以下步骤：
-
猜测一些聚类中心。
-
直到收敛重复：
-
E 步：将点分配给最近的聚类中心。
-
M 步：将聚类中心设置为其分配点的平均值。
-
这里的 E 步 或 期望步骤 之所以这样命名，是因为它涉及更新我们对每个点属于哪个聚类的期望。M 步 或 最大化步骤 之所以这样命名，是因为它涉及最大化某些定义聚类中心位置的适应函数——在本例中，通过简单地取每个聚类中数据的平均值来实现该最大化。
关于这一算法的文献非常丰富，但可以总结如下：在典型情况下，每次 E 步和 M 步的重复都会导致对聚类特征的更好估计。
我们可以将算法可视化如图 Figure 47-3 所示。对于此处显示的特定初始化，聚类在仅三次迭代中收敛。（有关此图的交互版本，请参阅在线 附录 中的代码。）
图 47-3. k-means 的 E-M 算法可视化^(1)
k-means 算法简单到我们可以用几行代码来编写它。以下是一个非常基本的实现（参见图 Figure 47-4）。
In [5]: from sklearn.metrics import pairwise_distances_argmin
def find_clusters(X, n_clusters, rseed=2):
# 1\. Randomly choose clusters
rng = np.random.RandomState(rseed)
i = rng.permutation(X.shape[0])[:n_clusters]
centers = X[i]
while True:
# 2a. Assign labels based on closest center
labels = pairwise_distances_argmin(X, centers)
# 2b. Find new centers from means of points
new_centers = np.array([X[labels == i].mean(0)
for i in range(n_clusters)])
# 2c. Check for convergence
if np.all(centers == new_centers):
break
centers = new_centers
return centers, labels
centers, labels = find_clusters(X, 4)
plt.scatter(X[:, 0], X[:, 1], c=labels,
s=50, cmap='viridis');
图 47-4. 使用 k-means 标记的数据
大多数经过良好测试的实现在底层会做更多事情，但上述函数传达了期望-最大化方法的主旨。在使用期望-最大化算法时，有几个需要注意的事项：
可能无法达到全局最优结果
首先，尽管 E-M 过程保证在每个步骤中改善结果，但不能保证它会导致 全局 最佳解。例如，如果在我们的简单过程中使用不同的随机种子，特定的起始猜测会导致糟糕的结果（参见图 Figure 47-5）。
In [6]: centers, labels = find_clusters(X, 4, rseed=0)
plt.scatter(X[:, 0], X[:, 1], c=labels,
s=50, cmap='viridis');
图 47-5. k-means 算法收敛不良的示例
在这里，E-M 方法已经收敛，但未收敛到全局最优配置。因此，通常会对算法使用多个起始猜测进行多次运行，默认情况下 Scikit-Learn 就是如此（该数字由 n_init 参数设置，默认为 10）。
必须事先选择聚类数
k-means 的另一个常见挑战是您必须告诉它您期望的聚类数：它无法从数据中学习到聚类数。例如，如果我们要求算法识别六个聚类，它将愉快地继续并找到最佳的六个聚类，如图 Figure 40-1 所示：
In [7]: labels = KMeans(6, random_state=0).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels,
s=50, cmap='viridis');
图 47-6. 簇数量选择不当的示例
结果是否有意义是一个很难明确回答的问题；一个相当直观的方法是使用轮廓分析，但我们这里不再进一步讨论。
或者，您可以使用更复杂的聚类算法，该算法对于每个簇的适应度有更好的定量衡量（例如，高斯混合模型；参见第四十八章），或者可以选择合适的簇数量（例如，DBSCAN、均值漂移或亲和力传播，这些都在sklearn.cluster子模块中提供）。
k-means 仅限于线性簇边界
k-means 的基本模型假设（点会更靠近自己的簇中心而不是其他簇）意味着如果簇具有复杂的几何结构，则该算法通常会失效。
特别地，k-means 簇之间的边界始终是线性的，这意味着对于更复杂的边界，它将失败。 考虑以下数据，以及典型k-means 方法找到的簇标签（见图 47-7）。
In [8]: from sklearn.datasets import make_moons
X, y = make_moons(200, noise=.05, random_state=0)
In [9]: labels = KMeans(2, random_state=0).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels,
s=50, cmap='viridis');
图 47-7. k-means 在非线性边界下的失败
这种情况让人想起了第四十三章中的讨论，在那里我们使用核变换将数据投影到更高的维度，从而可能实现线性分离。 我们可以想象使用相同的技巧来允许k-means 发现非线性边界。
这种基于核的k-means 的一个版本在 Scikit-Learn 中通过SpectralClustering估计器实现。 它使用最近邻图来计算数据的更高维表示，然后使用k-means 算法分配标签（参见图 47-8）。
In [10]: from sklearn.cluster import SpectralClustering
model = SpectralClustering(n_clusters=2,
affinity='nearest_neighbors',
assign_labels='kmeans')
labels = model.fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels,
s=50, cmap='viridis');
图 47-8. SpectralClustering 学习的非线性边界
我们看到，通过这种核变换方法，基于核的k-means 能够找到更复杂的簇之间的非线性边界。
对于大量样本，k-means 可能会运行缓慢。
因为k-means 的每次迭代都必须访问数据集中的每个点，所以随着样本数量的增长，该算法可能相对缓慢。 您可能会想知道是否可以放宽每次迭代使用所有数据的要求；例如，您可能只使用数据的子集来更新每个步骤的簇中心。 这就是批量式k-means 算法背后的思想，其中一种形式在sklearn.cluster.MiniBatchKMeans中实现。 其接口与标准的KMeans相同；我们将在继续讨论时看到其使用示例。
例子
虽然我们要注意算法的这些限制，但我们可以在各种情况下利用k-均值来获益。现在我们来看几个例子。
示例 1：数字上的k-均值
首先，让我们看看在我们在第四十四章 和第四十五章 中看到的相同简单数字数据上应用k-均值。在这里，我们将尝试使用k-均值来尝试识别类似的数字，而不使用原始标签信息；这可能类似于从一个没有任何先验标签信息的新数据集中提取含义的第一步。
我们将从加载数据集开始，然后找到聚类。回想一下，数字数据集包含 1,797 个样本，每个样本有 64 个特征，其中每个特征是 8 × 8 图像中一个像素的亮度。
In [11]: from sklearn.datasets import load_digits
digits = load_digits()
digits.data.shape
Out[11]: (1797, 64)
我们可以像之前一样执行聚类：
In [12]: kmeans = KMeans(n_clusters=10, random_state=0)
clusters = kmeans.fit_predict(digits.data)
kmeans.cluster_centers_.shape
Out[12]: (10, 64)
结果是 64 维空间中的 10 个聚类。请注意，聚类中心本身是 64 维点，可以解释为聚类内“典型”的数字。让我们看看这些聚类中心是什么样子的（见 图 47-9）。
In [13]: fig, ax = plt.subplots(2, 5, figsize=(8, 3))
centers = kmeans.cluster_centers_.reshape(10, 8, 8)
for axi, center in zip(ax.flat, centers):
axi.set(xticks=[], yticks=[])
axi.imshow(center, interpolation='nearest', cmap=plt.cm.binary)
图 47-9. k-均值学习到的聚类中心
我们看到，即使没有标签的情况下，KMeans 也能够找到其聚类中心可识别的数字，也许除了“1”和“8”之外。
因为k-均值对聚类的身份一无所知，0–9 标签可能会被排列。我们可以通过将每个学习到的聚类标签与聚类中找到的真实标签匹配来解决这个问题：
In [14]: from scipy.stats import mode
labels = np.zeros_like(clusters)
for i in range(10):
mask = (clusters == i)
labels[mask] = mode(digits.target[mask])[0]
现在我们可以检查我们的无监督聚类在找到数据中相似数字方面的准确性：
In [15]: from sklearn.metrics import accuracy_score
accuracy_score(digits.target, labels)
Out[15]: 0.7935447968836951
仅仅使用简单的k-均值算法，我们就为 80%的输入数字找到了正确的分组！让我们来查看这个混淆矩阵，它在 图 47-10 中可视化。
In [16]: from sklearn.metrics import confusion_matrix
import seaborn as sns
mat = confusion_matrix(digits.target, labels)
sns.heatmap(mat.T, square=True, annot=True, fmt='d',
cbar=False, cmap='Blues',
xticklabels=digits.target_names,
yticklabels=digits.target_names)
plt.xlabel('true label')
plt.ylabel('predicted label');
图 47-10. k-均值分类器的混淆矩阵
正如我们之前可视化的聚类中心所示，混淆的主要点在于数字“8”和“1”。但这仍然表明，使用k-均值，我们基本上可以建立一个数字分类器，无需参考任何已知标签！
仅仅是为了好玩，让我们尝试推动这个进展更远。我们可以在执行k-均值之前使用 t-分布随机邻居嵌入算法（在第四十六章中提到）对数据进行预处理。t-SNE 是一种非线性嵌入算法，特别擅长保留簇内的点。我们来看看它的表现：
In [17]: from sklearn.manifold import TSNE
# Project the data: this step will take several seconds
tsne = TSNE(n_components=2, init='random',
learning_rate='auto',random_state=0)
digits_proj = tsne.fit_transform(digits.data)
# Compute the clusters
kmeans = KMeans(n_clusters=10, random_state=0)
clusters = kmeans.fit_predict(digits_proj)
# Permute the labels
labels = np.zeros_like(clusters)
for i in range(10):
mask = (clusters == i)
labels[mask] = mode(digits.target[mask])[0]
# Compute the accuracy
accuracy_score(digits.target, labels)
Out[17]: 0.9415692821368948
这是一种不使用标签的 94% 分类准确率。这展示了无监督学习在谨慎使用时的强大能力：它可以从数据集中提取信息，这可能难以手工或肉眼提取。
示例 2：颜色压缩的k-均值
聚类的一个有趣应用是图像内的颜色压缩（此示例改编自 Scikit-Learn 的“使用 K-Means 进行颜色量化”）。例如，想象一下你有一幅包含数百万种颜色的图像。在大多数图像中，许多颜色将未被使用，并且图像中的许多像素将具有相似或甚至相同的颜色。
例如，考虑图像显示在图 47-11 中，这是来自 Scikit-Learn datasets模块的（为了使其工作，您必须安装PIL Python 包）：
In [18]: # Note: this requires the PIL package to be installed
from sklearn.datasets import load_sample_image
china = load_sample_image("china.jpg")
ax = plt.axes(xticks=[], yticks=[])
ax.imshow(china);
图 47-11. 输入图像
图像本身存储在一个大小为(height, width, RGB)的三维数组中，包含从 0 到 255 的整数表示的红/蓝/绿分量：
In [19]: china.shape
Out[19]: (427, 640, 3)
我们可以将这组像素视为三维色彩空间中的一组点云。我们将重新调整数据为[n_samples, n_features]的形状，并重新缩放颜色，使其介于 0 到 1 之间：
In [20]: data = china / 255.0 # use 0...1 scale
data = data.reshape(-1, 3)
data.shape
Out[20]: (273280, 3)
我们可以使用 10000 个像素的子集在此色彩空间中可视化这些像素（见图 47-12）。
In [21]: def plot_pixels(data, title, colors=None, N=10000):
if colors is None:
colors = data
# choose a random subset
rng = np.random.default_rng(0)
i = rng.permutation(data.shape[0])[:N]
colors = colors[i]
R, G, B = data[i].T
fig, ax = plt.subplots(1, 2, figsize=(16, 6))
ax[0].scatter(R, G, color=colors, marker='.')
ax[0].set(xlabel='Red', ylabel='Green', xlim=(0, 1), ylim=(0, 1))
ax[1].scatter(R, B, color=colors, marker='.')
ax[1].set(xlabel='Red', ylabel='Blue', xlim=(0, 1), ylim=(0, 1))
fig.suptitle(title, size=20);
In [22]: plot_pixels(data, title='Input color space: 16 million possible colors')
图 47-12. 在 RGB 色彩空间中的像素分布^(3)
现在让我们将这 1600 万种颜色减少到只有 16 种颜色，使用像素空间的* k -means 聚类。由于我们正在处理一个非常大的数据集，我们将使用小批量 k -means，它在数据子集上计算结果（显示在图 47-13 中）比标准 k *-means 算法更快：
In [23]: from sklearn.cluster import MiniBatchKMeans
kmeans = MiniBatchKMeans(16)
kmeans.fit(data)
new_colors = kmeans.cluster_centers_[kmeans.predict(data)]
plot_pixels(data, colors=new_colors,
title="Reduced color space: 16 colors")
图 47-13. RGB 色彩空间中的 16 个聚类^(4)
结果是原始像素的重新着色，其中每个像素分配到最接近的聚类中心的颜色。将这些新颜色在图像空间而不是像素空间中绘制，显示了这种效果（见图 47-14）。
In [24]: china_recolored = new_colors.reshape(china.shape)
fig, ax = plt.subplots(1, 2, figsize=(16, 6),
subplot_kw=dict(xticks=[], yticks=[]))
fig.subplots_adjust(wspace=0.05)
ax[0].imshow(china)
ax[0].set_title('Original Image', size=16)
ax[1].imshow(china_recolored)
ax[1].set_title('16-color Image', size=16);
图 47-14. 全彩图像（左）和 16 色图像（右）的比较
右侧面板中确实丢失了一些细节，但整体图像仍然很容易识别。在存储原始数据所需的字节方面，右侧的图像实现了约 100 万的压缩比！现在，这种方法不会与像 JPEG 这样的专用图像压缩方案匹敌，但这个例子展示了通过* k *-means 等无监督方法进行创新思维的威力。
^(1) 生成此图的代码可以在在线附录中找到。
^(2) 欲查看彩色版本以及后续图像，请参阅本书的在线版本。
^(3) 这幅图的全尺寸版本可以在GitHub上找到。
^(4) 这幅图的全尺寸版本可以在GitHub上找到。
第四十八章：深入：高斯混合模型
在前一章中探讨的 k-means 聚类模型简单且相对易于理解，但其简单性导致在实际应用中存在实际挑战。特别是，k-means 的非概率性质以及其使用简单的距离从聚类中心分配聚类成员导致在许多实际情况下性能不佳。在本章中，我们将介绍高斯混合模型，它可以被视为对 k-means 背后思想的扩展，同时也可以是一种超越简单聚类的强大工具。
我们从标准导入开始：
In [1]: %matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
import numpy as np
激励高斯混合模型：k-means 的弱点
让我们来看一些 k-means 的弱点，并思考如何改进聚类模型。正如我们在前一章中看到的，对于简单而明确分离的数据，k-means 能够找到合适的聚类结果。
例如，如果我们有简单的数据块，k-means 算法可以快速地标记这些聚类，其结果与我们可能通过眼睛观察到的相似（参见 图 48-1）。
In [2]: # Generate some data
from sklearn.datasets import make_blobs
X, y_true = make_blobs(n_samples=400, centers=4,
cluster_std=0.60, random_state=0)
X = X[:, ::-1] # flip axes for better plotting
In [3]: # Plot the data with k-means labels
from sklearn.cluster import KMeans
kmeans = KMeans(4, random_state=0)
labels = kmeans.fit(X).predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis');
图 48-1：简单数据的 k-means 标签
从直觉上讲，我们可能期望某些点的聚类分配比其他点更加确定：例如，在两个中间聚类之间似乎存在非常轻微的重叠，因此我们可能对它们之间的点的聚类分配没有完全的信心。不幸的是，k-means 模型没有内在的概率或聚类分配不确定性的衡量方法（尽管可能可以使用自举方法来估计此不确定性）。为此，我们必须考虑模型的泛化。
关于 k-means 模型的一种思考方式是，它在每个聚类的中心放置一个圆（或在更高维度中，一个超球体），其半径由聚类中最远的点定义。这个半径作为训练集内聚类分配的硬截止：任何在这个圆外的点都不被视为聚类的成员。我们可以用以下函数可视化这个聚类模型（参见 图 48-2）。
In [4]: from sklearn.cluster import KMeans
from scipy.spatial.distance import cdist
def plot_kmeans(kmeans, X, n_clusters=4, rseed=0, ax=None):
labels = kmeans.fit_predict(X)
# plot the input data
ax = ax or plt.gca()
ax.axis('equal')
ax.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis', zorder=2)
# plot the representation of the KMeans model
centers = kmeans.cluster_centers_
radii = [cdist(X[labels == i], [center]).max()
for i, center in enumerate(centers)]
for c, r in zip(centers, radii):
ax.add_patch(plt.Circle(c, r, ec='black', fc='lightgray',
lw=3, alpha=0.5, zorder=1))
In [5]: kmeans = KMeans(n_clusters=4, random_state=0)
plot_kmeans(kmeans, X)
图 48-2：k-means 模型暗示的圆形聚类
对于 k-means 的一个重要观察是，这些聚类模型必须是圆形的：k-means 没有内建的方法来处理椭圆形或椭圆形聚类。因此，例如，如果我们取同样的数据并对其进行转换，聚类分配最终变得混乱，正如你可以在 图 48-3 中看到的。
In [6]: rng = np.random.RandomState(13)
X_stretched = np.dot(X, rng.randn(2, 2))
kmeans = KMeans(n_clusters=4, random_state=0)
plot_kmeans(kmeans, X_stretched)
图 48-3：k-means 对于非圆形聚类的性能不佳
凭眼观察，我们认识到这些转换后的聚类不是圆形的，因此圆形聚类会拟合效果差。然而，k-means 不足以解决这个问题，并试图强行将数据拟合为四个圆形聚类。这导致了聚类分配的混合，其中结果的圆形重叠：尤其是在图的右下角可见。可以想象通过使用 PCA 预处理数据来处理这种情况（参见第四十五章），但实际上不能保证这样的全局操作会使各个群体圆形化。
k-means 的这两个缺点——在聚类形状上的灵活性不足和缺乏概率聚类分配——意味着对于许多数据集（特别是低维数据集），其性能可能不如人们所期望的那样好。
您可以通过泛化k-means 模型来解决这些弱点：例如，可以通过比较每个点到所有聚类中心的距离来测量聚类分配的不确定性，而不是仅关注最近的距离。您还可以想象允许聚类边界为椭圆而不是圆形，以适应非圆形聚类。事实证明，这些是不同类型聚类模型——高斯混合模型的两个基本组成部分。
泛化 E-M：高斯混合模型
高斯混合模型（GMM）试图找到最适合模拟任何输入数据集的多维高斯概率分布混合物。在最简单的情况下，GMM 可以像k-means 一样用于查找聚类（参见图 48-4）。
In [7]: from sklearn.mixture import GaussianMixture
gmm = GaussianMixture(n_components=4).fit(X)
labels = gmm.predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis');
图 48-4。数据的高斯混合模型标签
但是因为 GMM 在幕后包含一个概率模型，因此可以找到概率聚类分配——在 Scikit-Learn 中，这是通过predict_proba方法完成的。这将返回一个大小为[n_samples, n_clusters]的矩阵，用于测量任何点属于给定聚类的概率：
In [8]: probs = gmm.predict_proba(X)
print(probs[:5].round(3))
Out[8]: [[0. 0.531 0.469 0. ]
[0. 0. 0. 1. ]
[0. 0. 0. 1. ]
[0. 1. 0. 0. ]
[0. 0. 0. 1. ]]
我们可以通过将每个点的大小与其预测的确定性成比例来可视化这种不确定性；查看图 48-5，我们可以看到恰好在聚类边界上的点反映了聚类分配的不确定性：
In [9]: size = 50 * probs.max(1) ** 2 # square emphasizes differences
plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', s=size);
图 48-5。GMM 概率标签：点的大小显示概率大小
在幕后，高斯混合模型与k-means 非常相似：它使用期望最大化方法，大致如下：
-
选择位置和形状的初始猜测。
-
直到收敛为止重复：
-
E 步骤：对于每个点，找到编码每个聚类成员概率的权重。
-
M 步：对于每个聚类，根据所有数据点更新其位置、归一化和形状，利用权重。
-
其结果是，每个聚类不再关联于硬边界的球体，而是与平滑的高斯模型相关联。就像k-means 期望最大化方法一样，这种算法有时会错过全局最优解，因此在实践中使用多个随机初始化。
让我们创建一个函数，通过根据 GMM 输出绘制椭圆来帮助我们可视化 GMM 聚类的位置和形状：
In [10]: from matplotlib.patches import Ellipse
def draw_ellipse(position, covariance, ax=None, **kwargs):
"""Draw an ellipse with a given position and covariance"""
ax = ax or plt.gca()
# Convert covariance to principal axes
if covariance.shape == (2, 2):
U, s, Vt = np.linalg.svd(covariance)
angle = np.degrees(np.arctan2(U[1, 0], U[0, 0]))
width, height = 2 * np.sqrt(s)
else:
angle = 0
width, height = 2 * np.sqrt(covariance)
# Draw the ellipse
for nsig in range(1, 4):
ax.add_patch(Ellipse(position, nsig * width, nsig * height,
angle, **kwargs))
def plot_gmm(gmm, X, label=True, ax=None):
ax = ax or plt.gca()
labels = gmm.fit(X).predict(X)
if label:
ax.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis',
zorder=2)
else:
ax.scatter(X[:, 0], X[:, 1], s=40, zorder=2)
ax.axis('equal')
w_factor = 0.2 / gmm.weights_.max()
for pos, covar, w in zip(gmm.means_, gmm.covariances_, gmm.weights_):
draw_ellipse(pos, covar, alpha=w * w_factor)
有了这些基础，我们可以看看四分量 GMM 对我们的初始数据给出了什么结果（参见图 48-6）。
In [11]: gmm = GaussianMixture(n_components=4, random_state=42)
plot_gmm(gmm, X)
图 48-6. 存在圆形聚类的四分量 GMM
同样地，我们可以使用 GMM 方法拟合我们的伸展数据集；允许完全协方差模型将适合甚至是非常椭圆形、拉伸的聚类，正如我们在图 48-7 中所看到的。
In [12]: gmm = GaussianMixture(n_components=4, covariance_type='full',
random_state=42)
plot_gmm(gmm, X_stretched)
图 48-7. 存在非圆形聚类的四分量 GMM
这清楚地表明，GMM 解决了之前在k-means 中遇到的两个主要实际问题。
选择协方差类型
如果您查看前面拟合的细节，您会发现在每个拟合中设置了covariance_type选项。该超参数控制每个聚类形状的自由度；对于任何给定的问题，仔细设置这一点至关重要。默认值是covariance_type="diag"，这意味着可以独立设置每个维度上的聚类大小，生成的椭圆受限于与轴对齐。covariance_type="spherical"是一个稍微简单且更快的模型，它限制了聚类形状，使得所有维度相等。结果聚类将具有与k-means 类似的特征，尽管它并非完全等价。一个更复杂和计算开销更大的模型（特别是在维度增长时）是使用covariance_type="full"，它允许将每个聚类建模为带有任意方向的椭圆。图 48-8 表示了这三种选择对单个聚类的影响。
图 48-8. GMM 协方差类型可视化^(1)
高斯混合模型作为密度估计
尽管 GMM 通常被归类为聚类算法，但从根本上讲，它是一种用于密度估计的算法。也就是说，对某些数据进行 GMM 拟合的结果在技术上不是聚类模型，而是描述数据分布的生成概率模型。
以 Scikit-Learn 的make_moons函数生成的数据为例，介绍在第四十七章中。
In [13]: from sklearn.datasets import make_moons
Xmoon, ymoon = make_moons(200, noise=.05, random_state=0)
plt.scatter(Xmoon[:, 0], Xmoon[:, 1]);
图 48-9. GMM 应用于具有非线性边界的聚类
如果我们尝试用一个两组分的 GMM 作为聚类模型来拟合它，结果并不特别有用（见图 48-10）。
In [14]: gmm2 = GaussianMixture(n_components=2, covariance_type='full',
random_state=0)
plot_gmm(gmm2, Xmoon)
图 48-10. 对非线性聚类拟合的两组分 GMM
但是，如果我们使用更多组分并忽略聚类标签，我们会发现拟合结果更接近输入数据（见图 48-11）。
In [15]: gmm16 = GaussianMixture(n_components=16, covariance_type='full',
random_state=0)
plot_gmm(gmm16, Xmoon, label=False)
图 48-11. 使用多个 GMM 组件来建模点分布
这里的 16 个高斯分量的混合并不是为了找到数据的分离聚类，而是为了对输入数据的整体分布进行建模。这是一个生成模型，意味着 GMM 给了我们一个生成新随机数据的方法，其分布类似于我们的原始输入数据。例如，这里有 400 个新点从这个 16 组分的 GMM 拟合到我们的原始数据中绘制出来（见图 48-12）。
In [16]: Xnew, ynew = gmm16.sample(400)
plt.scatter(Xnew[:, 0], Xnew[:, 1]);
图 48-12. 从 16 组分 GMM 中绘制的新数据
GMM 作为一种灵活的方法，方便地对数据的任意多维分布进行建模。
GMM 作为生成模型的事实给了我们一种自然的方法来确定给定数据集的最优组件数。生成模型本质上是数据集的概率分布，因此我们可以简单地在模型下评估数据的似然性，使用交叉验证来避免过度拟合。另一种校正过度拟合的方法是使用一些分析标准来调整模型的似然性，例如阿卡奇信息准则（AIC）或贝叶斯信息准则（BIC）。Scikit-Learn 的GaussianMixture估计器实际上包含内置方法来计算这两者，因此使用这种方法非常容易。
让我们看看我们的 moons 数据集的 GMM 组件数对应的 AIC 和 BIC（见图 48-13）。
In [17]: n_components = np.arange(1, 21)
models = [GaussianMixture(n, covariance_type='full',
random_state=0).fit(Xmoon)
for n in n_components]
plt.plot(n_components, [m.bic(Xmoon) for m in models], label='BIC')
plt.plot(n_components, [m.aic(Xmoon) for m in models], label='AIC')
plt.legend(loc='best')
plt.xlabel('n_components');
图 48-13. AIC 和 BIC 的可视化，用于选择 GMM 组件数
最优的聚类数是能够最小化 AIC 或 BIC 的值，具体取决于我们希望使用哪种近似方法。AIC 告诉我们，我们之前选择的 16 组分可能太多了：选择大约 8-12 组分可能更合适。对于这类问题，BIC 通常推荐一个更简单的模型。
注意重要的一点：组件数量的选择衡量的是 GMM 作为密度估计器的工作效果，而不是作为聚类算法的工作效果。我鼓励您主要将 GMM 视为密度估计器，并仅在简单数据集内合适时用它进行聚类。
示例：使用 GMM 生成新数据
我们刚刚看到了使用 GMM 作为生成模型的简单示例，以便从定义为输入数据分布的模型中创建新的样本。在这里，我们将继续这个想法，并从之前使用过的标准数字语料库中生成新的手写数字。
首先，让我们使用 Scikit-Learn 的数据工具加载数字数据：
In [18]: from sklearn.datasets import load_digits
digits = load_digits()
digits.data.shape
Out[18]: (1797, 64)
接下来，让我们绘制前 50 个样本，以确切回顾我们正在查看的内容（参见图 48-14）。
In [19]: def plot_digits(data):
fig, ax = plt.subplots(5, 10, figsize=(8, 4),
subplot_kw=dict(xticks=[], yticks=[]))
fig.subplots_adjust(hspace=0.05, wspace=0.05)
for i, axi in enumerate(ax.flat):
im = axi.imshow(data[i].reshape(8, 8), cmap='binary')
im.set_clim(0, 16)
plot_digits(digits.data)
图 48-14. 手写数字输入
我们有将近 1,800 个 64 维度的数字样本，我们可以在其上构建一个混合高斯模型（GMM）以生成更多数字。在这么高维度的空间中，GMM 可能会有收敛困难，因此我们将从数据中开始使用一个可逆的降维算法。这里我们将使用简单的 PCA，要求它在投影数据中保留 99%的方差：
In [20]: from sklearn.decomposition import PCA
pca = PCA(0.99, whiten=True)
data = pca.fit_transform(digits.data)
data.shape
Out[20]: (1797, 41)
结果是 41 个维度，几乎没有信息损失的减少了近 1/3。鉴于这个投影数据，让我们使用 AIC 来确定我们应该使用多少个 GMM 组件（参见图 48-15）。
In [21]: n_components = np.arange(50, 210, 10)
models = [GaussianMixture(n, covariance_type='full', random_state=0)
for n in n_components]
aics = [model.fit(data).aic(data) for model in models]
plt.plot(n_components, aics);
图 48-15. 选择适当的 GMM 组件数量的 AIC 曲线
看起来大约使用 140 个组件可以最小化 AIC；我们将使用这个模型。让我们快速将其拟合到数据上并确认它已经收敛：
In [22]: gmm = GaussianMixture(140, covariance_type='full', random_state=0)
gmm.fit(data)
print(gmm.converged_)
Out[22]: True
现在我们可以在这个 41 维度的投影空间内绘制 100 个新点的样本，使用 GMM 作为生成模型：
In [23]: data_new, label_new = gmm.sample(100)
data_new.shape
Out[23]: (100, 41)
最后，我们可以使用 PCA 对象的逆变换来构造新的数字（参见图 48-16）。
In [24]: digits_new = pca.inverse_transform(data_new)
plot_digits(digits_new)
图 48-16. 从 GMM 估计器的基础模型中随机绘制的“新”数字
大多数结果看起来像数据集中合理的数字！
考虑我们在这里所做的：鉴于手写数字的抽样，我们已经模拟了该数据的分布，以便我们可以从数据中生成全新的样本：这些是“手写数字”，它们不会单独出现在原始数据集中，而是捕捉了混合模型建模的输入数据的一般特征。这样的手写数字的生成模型在贝叶斯生成分类器的组成部分中可以非常有用，这一点我们将在下一章看到。
^(1) 生成此图的代码可以在在线附录中找到。
第四十九章：深入探讨：核密度估计
在第四十八章中，我们讨论了高斯混合模型，这是一种聚类估计器和密度估计器之间的混合类型。回想一下，密度估计器是一种算法，它接受一个D维数据集，并生成数据抽取自其中的D维概率分布的估计。GMM 算法通过将密度表示为高斯分布的加权和来实现这一点。核密度估计（KDE）在某种意义上是将高斯混合思想推向其逻辑极限的算法：它使用每个点一个高斯分量的混合，从而得到一个基本上是非参数的密度估计器。在本章中，我们将探讨 KDE 的动机和用途。
我们从标准导入开始：
In [1]: %matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
import numpy as np
激发核密度估计：直方图
如前所述，密度估计器是一种算法，旨在模拟生成数据集的概率分布。对于一维数据，你可能已经熟悉一个简单的密度估计器：直方图。直方图将数据分成离散的箱子，计算落入每个箱子的点的数量，然后直观地可视化结果。
例如，让我们创建一些从两个正态分布中绘制的数据：
In [2]: def make_data(N, f=0.3, rseed=1):
rand = np.random.RandomState(rseed)
x = rand.randn(N)
x[int(f * N):] += 5
return x
x = make_data(1000)
我们之前看到，可以通过指定直方图的density参数来创建标准的基于计数的直方图。通过这种方式得到的是归一化的直方图，其中箱子的高度不反映计数，而是反映概率密度（参见图 49-1）。
In [3]: hist = plt.hist(x, bins=30, density=True)
图 49-1. 从正态分布组合中绘制的数据
对于等距分箱，这种归一化仅仅改变了 y 轴的比例，使得相对高度基本上与计数直方图中的相同。选择这种归一化是为了使直方图下面积等于 1，我们可以通过直方图函数的输出来确认：
In [4]: density, bins, patches = hist
widths = bins[1:] - bins[:-1]
(density * widths).sum()
Out[4]: 1.0
使用直方图作为密度估计器的一个问题是，箱子的大小和位置的选择可能导致具有不同特征的表现。例如，如果我们查看只有 20 个点的此数据的一个版本，如何绘制箱子的选择可以导致对数据的完全不同解释！考虑到这个例子，在图 49-2 中可视化。
In [5]: x = make_data(20)
bins = np.linspace(-5, 10, 10)
In [6]: fig, ax = plt.subplots(1, 2, figsize=(12, 4),
sharex=True, sharey=True,
subplot_kw={'xlim':(-4, 9),
'ylim':(-0.02, 0.3)})
fig.subplots_adjust(wspace=0.05)
for i, offset in enumerate([0.0, 0.6]):
ax[i].hist(x, bins=bins + offset, density=True)
ax[i].plot(x, np.full_like(x, -0.01), '|k',
markeredgewidth=1)
图 49-2. 直方图的问题：箱子的位置会影响解释
左侧，直方图清楚地显示这是一个双峰分布。右侧，我们看到一个长尾单峰分布。如果不看前面的代码，你可能不会猜到这两个直方图是由同一组数据构建的。有了这个认识，我们如何相信直方图所传达的直觉呢？我们如何改进这一点呢？
总结一下，我们可以把直方图看作是一堆块，我们在数据集中的每个点上都堆叠一个块。让我们直接查看这一点（见图 49-3）。
In [7]: fig, ax = plt.subplots()
bins = np.arange(-3, 8)
ax.plot(x, np.full_like(x, -0.1), '|k',
markeredgewidth=1)
for count, edge in zip(*np.histogram(x, bins)):
for i in range(count):
ax.add_patch(plt.Rectangle(
(edge, i), 1, 1, ec='black', alpha=0.5))
ax.set_xlim(-4, 8)
ax.set_ylim(-0.2, 8)
Out[7]: (-0.2, 8.0)
图 49-3. 堆叠的块直方图
我们两次分箱的问题源于这样一个事实：块堆叠的高度经常反映不出附近点的实际密度，而是由于分箱与数据点对齐的巧合。这种点与它们块之间的不对齐可能是导致这里糟糕直方图结果的一个潜在原因。但是，如果我们不是将块与分箱对齐，而是将块与它们所代表的点对齐会怎样呢？如果我们这样做，块就不会对齐，但我们可以在每个 x 轴位置上添加它们的贡献来找到结果。让我们试试这个方法（见图 49-4）。
In [8]: x_d = np.linspace(-4, 8, 2000)
density = sum((abs(xi - x_d) < 0.5) for xi in x)
plt.fill_between(x_d, density, alpha=0.5)
plt.plot(x, np.full_like(x, -0.1), '|k', markeredgewidth=1)
plt.axis([-4, 8, -0.2, 8]);
图 49-4. 一个“直方图”，其中每个块都以各个个体点为中心；这是一个核密度估计的示例
结果看起来有些杂乱，但它比标准直方图更能鲜明地反映实际数据特征。尽管如此，其粗糙的边缘既不美观，也不能反映数据的任何真实特性。为了平滑它们，我们可以决定在每个位置用一个平滑函数来取代这些块，比如一个高斯函数。让我们在每个点上使用一个标准正态曲线代替一个块（见图 49-5）。
In [9]: from scipy.stats import norm
x_d = np.linspace(-4, 8, 1000)
density = sum(norm(xi).pdf(x_d) for xi in x)
plt.fill_between(x_d, density, alpha=0.5)
plt.plot(x, np.full_like(x, -0.1), '|k', markeredgewidth=1)
plt.axis([-4, 8, -0.2, 5]);
图 49-5. 一个使用高斯核的核密度估计
这个平滑处理后的图，每个输入点处贡献一个高斯分布，更准确地反映了数据分布的形状，并且具有更低的方差（即对不同采样差异的响应更小）。
在过去的两个图中，我们所得到的就是一维核密度估计：我们在每个点的位置放置一个“核”—在前者中是一个方形或者顶帽形的核，在后者中是一个高斯核，并用它们的总和作为密度的估计。有了这个直觉，我们现在将更详细地探讨核密度估计。
实践中的核密度估计
核密度估计的自由参数包括 核函数，它指定放置在每个点处的分布的形状，以及 核带宽，它控制每个点处核的大小。实际上，可以使用许多核函数进行核密度估计：特别是，Scikit-Learn 的 KDE 实现支持六种核函数，你可以在“密度估计”部分的文档中了解更多信息。
虽然 Python 中有几个实现 KDE 的版本（特别是在 SciPy 和 statsmodels 包中），但我更倾向于使用 Scikit-Learn 的版本，因为它高效且灵活。它是在 sklearn.neighbors.KernelDensity 估计器中实现的，可以使用六种核函数和几十种距离度量来处理多维 KDE。由于 KDE 可能计算量较大，Scikit-Learn 的估计器在底层使用基于树的算法，并可以通过 atol（绝对容差）和 rtol（相对容差）参数在计算时间和准确性之间进行权衡。核带宽可以使用 Scikit-Learn 的标准交叉验证工具来确定，这很快我们会看到。
让我们首先展示一个简单的示例，使用 Scikit-Learn 的 KernelDensity 估计器复制先前的图（参见 Figure 49-6）。
In [10]: from sklearn.neighbors import KernelDensity
# instantiate and fit the KDE model
kde = KernelDensity(bandwidth=1.0, kernel='gaussian')
kde.fit(x[:, None])
# score_samples returns the log of the probability density
logprob = kde.score_samples(x_d[:, None])
plt.fill_between(x_d, np.exp(logprob), alpha=0.5)
plt.plot(x, np.full_like(x, -0.01), '|k', markeredgewidth=1)
plt.ylim(-0.02, 0.22);
图 49-6. 使用 Scikit-Learn 计算的核密度估计
此处的结果已归一化，使得曲线下面积等于 1。
通过交叉验证选择带宽
KDE 过程产生的最终估计对带宽的选择非常敏感，带宽是控制密度估计中偏差-方差权衡的旋钮。带宽太窄会导致高方差估计（即过拟合），其中单个点的存在或缺失会产生较大差异。带宽太宽会导致高偏差估计（即欠拟合），数据结构被广核模糊化。
在统计学中，有很长的历史可以快速估计基于数据的最佳带宽，基于对数据的严格假设：例如，如果你查看 SciPy 和 statsmodels 包中的 KDE 实现，你会看到基于这些规则的实现。
在机器学习环境中，我们看到这种超参数调整通常通过经验交叉验证方法完成。考虑到这一点，Scikit-Learn 的 KernelDensity 估计器设计成可以直接在包的标准网格搜索工具中使用。在这里，我们将使用 GridSearchCV 来优化前述数据集的带宽。由于我们正在查看一个如此小的数据集，我们将使用留一法交叉验证，以最小化每个交叉验证试验的训练集大小减少：
In [11]: from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import LeaveOneOut
bandwidths = 10 ** np.linspace(-1, 1, 100)
grid = GridSearchCV(KernelDensity(kernel='gaussian'),
{'bandwidth': bandwidths},
cv=LeaveOneOut())
grid.fit(x[:, None]);
现在我们可以找到使得分数最大化的带宽选择（在这种情况下，默认为对数似然）：
In [12]: grid.best_params_
Out[12]: {'bandwidth': 1.1233240329780276}
最优带宽与我们在之前示例图中使用的非常接近，那里的带宽是 1.0（即scipy.stats.norm的默认宽度）。
示例：不那么朴素的贝叶斯
此示例探讨了带 KDE 的贝叶斯生成分类，并演示了如何使用 Scikit-Learn 架构创建自定义估计器。
在第四十一章中，我们探讨了朴素贝叶斯分类，其中我们为每个类别创建了一个简单的生成模型，并使用这些模型构建了一个快速分类器。对于高斯朴素贝叶斯，生成模型是一个简单的轴对齐高斯分布。使用 KDE 等密度估计算法，我们可以去除“朴素”元素，并使用更复杂的生成模型为每个类别执行相同的分类。它仍然是贝叶斯分类，但不再是朴素的。
生成分类的一般方法如下：
-
根据标签将训练数据进行拆分。
-
对每个集合，拟合一个 KDE 以获得数据的生成模型。这允许你对于任意观测值x和标签y，计算出一个似然概率P ( x | y )。
-
根据训练集中每个类别的示例数量，计算类先验P ( y )。
-
对于未知点x，每个类别的后验概率为P ( y | x ) ∝ P ( x | y ) P ( y )。最大化这个后验概率的类别是分配给该点的标签。
算法很简单直观易懂；更难的部分是将其嵌入 Scikit-Learn 框架中，以便利用网格搜索和交叉验证架构。
这是在 Scikit-Learn 框架中实现算法的代码；我们将在代码块后面逐步分析它：
In [13]: from sklearn.base import BaseEstimator, ClassifierMixin
class KDEClassifier(BaseEstimator, ClassifierMixin):
"""Bayesian generative classification based on KDE
Parameters
----------
bandwidth : float
the kernel bandwidth within each class
kernel : str
the kernel name, passed to KernelDensity
"""
def __init__(self, bandwidth=1.0, kernel='gaussian'):
self.bandwidth = bandwidth
self.kernel = kernel
def fit(self, X, y):
self.classes_ = np.sort(np.unique(y))
training_sets = [X[y == yi] for yi in self.classes_]
self.models_ = [KernelDensity(bandwidth=self.bandwidth,
kernel=self.kernel).fit(Xi)
for Xi in training_sets]
self.logpriors_ = [np.log(Xi.shape[0] / X.shape[0])
for Xi in training_sets]
return self
def predict_proba(self, X):
logprobs = np.array([model.score_samples(X)
for model in self.models_]).T
result = np.exp(logprobs + self.logpriors_)
return result / result.sum(axis=1, keepdims=True)
def predict(self, X):
return self.classes_[np.argmax(self.predict_proba(X), 1)]
自定义估计器的解剖学。
让我们逐步分析这段代码，并讨论其关键特性：
from sklearn.base import BaseEstimator, ClassifierMixin
class KDEClassifier(BaseEstimator, ClassifierMixin):
"""Bayesian generative classification based on KDE
Parameters
----------
bandwidth : float
the kernel bandwidth within each class
kernel : str
the kernel name, passed to KernelDensity
"""
Scikit-Learn 中的每个评估器都是一个类，最方便的是这个类也应该从 BaseEstimator 类以及适当的 mixin 继承，提供标准功能。例如，在这里，BaseEstimator 包含了克隆/复制评估器以供交叉验证过程使用的必要逻辑，而 ClassifierMixin 定义了这些例程使用的默认 score 方法。我们还提供了一个文档字符串，这将被 IPython 的帮助功能捕获（参见 第一章）。
接下来是类的初始化方法：
def __init__(self, bandwidth=1.0, kernel='gaussian'):
self.bandwidth = bandwidth
self.kernel = kernel
当使用 KDEClassifier 实例化对象时，执行的实际代码是这样的。在 Scikit-Learn 中，初始化不包含任何操作，只是将传递的值按名称分配给 self。这是由于 BaseEstimator 中包含的逻辑，用于克隆和修改评估器，以供交叉验证、网格搜索和其他功能使用。同样，所有传递给 __init__ 的参数都应该是明确的：即应避免使用 *args 或 **kwargs，因为它们在交叉验证过程中无法正确处理。
下面是 fit 方法，我们在这里处理训练数据：
def fit(self, X, y):
self.classes_ = np.sort(np.unique(y))
training_sets = [X[y == yi] for yi in self.classes_]
self.models_ = [KernelDensity(bandwidth=self.bandwidth,
kernel=self.kernel).fit(Xi)
for Xi in training_sets]
self.logpriors_ = [np.log(Xi.shape[0] / X.shape[0])
for Xi in training_sets]
return self
在这里，我们找到训练数据中的唯一类别，为每个类别训练一个 KernelDensity 模型，并基于输入样本的数量计算类别先验概率。最后，fit 应该始终返回 self，以便我们可以链接命令。例如：
label = model.fit(X, y).predict(X)
注意，每次 fit 持久结果都带有下划线结尾（例如 self.logpriors_）。这是 Scikit-Learn 中的一种约定，因此您可以快速扫描评估器的成员（使用 IPython 的 tab 自动完成），并查看哪些成员是适合训练数据的。
最后，我们有预测新数据标签的逻辑：
def predict_proba(self, X):
logprobs = np.vstack([model.score_samples(X)
for model in self.models_]).T
result = np.exp(logprobs + self.logpriors_)
return result / result.sum(axis=1, keepdims=True)
def predict(self, X):
return self.classes_[np.argmax(self.predict_proba(X), 1)]
因为这是一个概率分类器，我们首先实现 predict_proba，它返回一个形状为 [n_samples, n_classes] 的类别概率数组。数组中的条目 [i, j] 是计算得到的样本 i 是类别 j 的后验概率，通过将似然性乘以类先验并进行归一化计算得到。
predict 方法使用这些概率，简单地返回具有最大概率的类别。
使用我们的自定义评估器
让我们尝试将这个自定义评估器应用于我们之前见过的问题：手写数字的分类。在这里，我们将加载数字并使用 GridSearchCV 元评估器计算一系列候选带宽的交叉验证分数（参考 第三十九章）：
In [14]: from sklearn.datasets import load_digits
from sklearn.model_selection import GridSearchCV
digits = load_digits()
grid = GridSearchCV(KDEClassifier(),
{'bandwidth': np.logspace(0, 2, 100)})
grid.fit(digits.data, digits.target);
接下来，我们可以绘制交叉验证分数作为带宽的函数（参见 图 49-7）。
In [15]: fig, ax = plt.subplots()
ax.semilogx(np.array(grid.cv_results_['param_bandwidth']),
grid.cv_results_['mean_test_score'])
ax.set(title='KDE Model Performance', ylim=(0, 1),
xlabel='bandwidth', ylabel='accuracy')
print(f'best param: {grid.best_params_}')
print(f'accuracy = {grid.best_score_}')
Out[15]: best param: {'bandwidth': 6.135907273413174}
accuracy = 0.9677298050139276
图 49-7. 基于 KDE 的贝叶斯分类器的验证曲线
这表明我们的 KDE 分类器达到了超过 96%的交叉验证准确率，而朴素贝叶斯分类器的准确率约为 80%。
In [16]: from sklearn.naive_bayes import GaussianNB
from sklearn.model_selection import cross_val_score
cross_val_score(GaussianNB(), digits.data, digits.target).mean()
Out[16]: 0.8069281956050759
这样一个生成分类器的一个好处是结果的可解释性：对于每个未知样本，我们不仅获得一个概率分类，而且获得了与其比较的点分布的完整模型！如果需要，这为特定分类的原因提供了一个直观的窗口，而像 SVM 和随机森林这样的算法往往会掩盖这些原因。
如果您希望进一步进行，这里有一些可以改进我们的 KDE 分类器模型的想法：
-
您可以允许每个类别中的带宽独立变化。
-
不应该基于它们的预测分数来优化这些带宽，而是应该基于每个类别中生成模型下训练数据的可能性（即使用
KernelDensity本身的分数而不是全局预测准确度）。
最后，如果你想要一些练习来构建自己的估计器，可以尝试使用高斯混合模型而不是 KDE 来构建类似的贝叶斯分类器。
第五十章：应用：一个人脸检测流水线
本书的这一部分探讨了许多机器学习的中心概念和算法。但是从这些概念到真实世界的应用可能是一个挑战。真实世界的数据集通常是嘈杂和异构的；它们可能具有缺失的特征，并且数据可能以难以映射到干净的[n_samples, n_features]矩阵的形式存在。在应用这里讨论的任何方法之前，您必须首先从您的数据中提取这些特征：没有适用于所有领域的公式，因此这是您作为数据科学家必须运用自己的直觉和专业知识的地方。
机器学习的一个有趣而引人注目的应用是图像，我们已经看到了一些例子，其中像素级特征用于分类。再次强调，现实世界的数据很少是如此统一的，简单的像素将不合适：这导致了大量关于图像数据特征提取方法的文献（参见第四十章）。
在本章中，我们将介绍一种特征提取技术：方向梯度直方图（HOG），它将图像像素转换为对广泛信息的敏感向量表示，而不受照明等混淆因素的影响。我们将使用这些特征来开发一个简单的人脸检测流水线，利用本书这部分中已经介绍过的机器学习算法和概念。
我们从标准导入开始：
In [1]: %matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
import numpy as np
HOG 特征
HOG 是一个简单直接的特征提取过程，最初用于图像中行人的识别。它包括以下步骤：
-
可选择地对图像进行预归一化。这导致特征对照明变化的依赖性较小。
-
将图像与两个对水平和垂直亮度梯度敏感的滤波器卷积。这些捕捉边缘、轮廓和纹理信息。
-
将图像细分为预定大小的单元格，并计算每个单元格内梯度方向的直方图。
-
通过与相邻单元格块比较来归一化每个单元格中的直方图。这进一步抑制了整个图像中照明效果的影响。
-
从每个单元格中的信息构建一个一维特征向量。
Scikit-Image 项目中内置了一个快速的 HOG 特征提取器，我们可以相对快速地尝试并可视化每个单元格内的定向梯度（见图 50-1）。
In [2]: from skimage import data, color, feature
import skimage.data
image = color.rgb2gray(data.chelsea())
hog_vec, hog_vis = feature.hog(image, visualize=True)
fig, ax = plt.subplots(1, 2, figsize=(12, 6),
subplot_kw=dict(xticks=[], yticks=[]))
ax[0].imshow(image, cmap='gray')
ax[0].set_title('input image')
ax[1].imshow(hog_vis)
ax[1].set_title('visualization of HOG features');
图 50-1。从图像计算的 HOG 特征的可视化
HOG 在行动：一个简单的人脸检测器
利用这些 HOG 特征，我们可以使用任何 Scikit-Learn 评估器构建一个简单的面部检测算法；在这里，我们将使用线性支持向量机（如果需要恢复记忆，请参阅 第四十三章）。具体步骤如下：
-
获得一组人脸缩略图，作为“正”训练样本。
-
获得一组非面部图像缩略图，作为“负”训练样本。
-
从这些训练样本中提取 HOG 特征。
-
在这些样本上训练线性 SVM 分类器。
-
为“未知”图像，通过图像上的滑动窗口，使用模型评估该窗口是否包含面部。
-
如果检测重叠，将它们合并成一个单一窗口。
让我们按照这些步骤进行并尝试一下。
1. 获得一组正训练样本
我们将从中找出一些显示各种面部的正训练样本。我们有一个易于使用的数据集——Wild 中的带标签面部数据集，可以通过 Scikit-Learn 下载：
In [3]: from sklearn.datasets import fetch_lfw_people
faces = fetch_lfw_people()
positive_patches = faces.images
positive_patches.shape
Out[3]: (13233, 62, 47)
这为我们提供了一些用于训练的 13,000 张面部图像样本。
2. 获得一组负训练样本
接下来，我们需要一组大小相似的缩略图，其中没有面部。获得这些的一种方法是从任何输入图像语料库中提取它们的缩略图，以多种比例。在这里，我们将使用一些随着 Scikit-Image 提供的图像以及 Scikit-Learn 的 PatchExtractor：
In [4]: data.camera().shape
Out[4]: (512, 512)
In [5]: from skimage import data, transform
imgs_to_use = ['camera', 'text', 'coins', 'moon',
'page', 'clock', 'immunohistochemistry',
'chelsea', 'coffee', 'hubble_deep_field']
raw_images = (getattr(data, name)() for name in imgs_to_use)
images = [color.rgb2gray(image) if image.ndim == 3 else image
for image in raw_images]
In [6]: from sklearn.feature_extraction.image import PatchExtractor
def extract_patches(img, N, scale=1.0, patch_size=positive_patches[0].shape):
extracted_patch_size = tuple((scale * np.array(patch_size)).astype(int))
extractor = PatchExtractor(patch_size=extracted_patch_size,
max_patches=N, random_state=0)
patches = extractor.transform(img[np.newaxis])
if scale != 1:
patches = np.array([transform.resize(patch, patch_size)
for patch in patches])
return patches
negative_patches = np.vstack([extract_patches(im, 1000, scale)
for im in images for scale in [0.5, 1.0, 2.0]])
negative_patches.shape
Out[6]: (30000, 62, 47)
现在我们有 30,000 个合适的图像补丁，不包含面部。让我们可视化其中的一些，以了解它们的外观（见 图 50-2）。
In [7]: fig, ax = plt.subplots(6, 10)
for i, axi in enumerate(ax.flat):
axi.imshow(negative_patches[500 * i], cmap='gray')
axi.axis('off')
我们希望这些样本足以覆盖算法可能见到的“非面部”空间。
图 50-2. 不包含面部的负图像补丁
3. 合并集合并提取 HOG 特征
现在我们有了这些正样本和负样本，我们可以将它们合并并计算 HOG 特征。这一步骤需要一些时间，因为它涉及到每个图像的非平凡计算：
In [8]: from itertools import chain
X_train = np.array([feature.hog(im)
for im in chain(positive_patches,
negative_patches)])
y_train = np.zeros(X_train.shape[0])
y_train[:positive_patches.shape[0]] = 1
In [9]: X_train.shape
Out[9]: (43233, 1215)
我们剩下 43,000 个训练样本，具有 1,215 个维度，现在我们已经将数据处理成了可以输入到 Scikit-Learn 的形式！
4. 训练支持向量机
接下来，我们将利用这里探索的工具创建缩略图补丁的分类器。对于这种高维二元分类任务，线性支持向量机是一个很好的选择。我们将使用 Scikit-Learn 的 LinearSVC，因为与 SVC 相比，它通常对大量样本具有更好的扩展性。
不过首先，让我们使用简单的高斯朴素贝叶斯估算器得到一个快速的基准：
In [10]: from sklearn.naive_bayes import GaussianNB
from sklearn.model_selection import cross_val_score
cross_val_score(GaussianNB(), X_train, y_train)
Out[10]: array([0.94795883, 0.97143518, 0.97224471, 0.97501735, 0.97374508])
我们看到在我们的训练数据上，即使是简单的朴素贝叶斯算法也可以达到 95% 以上的准确率。让我们尝试支持向量机，并对几种 C 参数进行网格搜索：
In [11]: from sklearn.svm import LinearSVC
from sklearn.model_selection import GridSearchCV
grid = GridSearchCV(LinearSVC(), {'C': [1.0, 2.0, 4.0, 8.0]})
grid.fit(X_train, y_train)
grid.best_score_
Out[11]: 0.9885272620319941
In [12]: grid.best_params_
Out[12]: {'C': 1.0}
这将使我们的准确率提高到接近 99%。让我们选择最佳的估计器，并在完整数据集上重新训练它：
In [13]: model = grid.best_estimator_
model.fit(X_train, y_train)
Out[13]: LinearSVC()
5. 在新图像中查找面部
现在我们已经建立了这个模型，让我们拿一幅新的图像来看看模型的表现。我们将简单地使用宇航员图像中的一部分，如 图 50-3 所示，运行一个滑动窗口，并评估每个补丁：
In [14]: test_image = skimage.data.astronaut()
test_image = skimage.color.rgb2gray(test_image)
test_image = skimage.transform.rescale(test_image, 0.5)
test_image = test_image[:160, 40:180]
plt.imshow(test_image, cmap='gray')
plt.axis('off');
图 50-3. 一幅我们将尝试定位面部的图像。
接下来，让我们创建一个窗口，迭代这幅图像的补丁，并为每个补丁计算 HOG 特征：
In [15]: def sliding_window(img, patch_size=positive_patches[0].shape,
istep=2, jstep=2, scale=1.0):
Ni, Nj = (int(scale * s) for s in patch_size)
for i in range(0, img.shape[0] - Ni, istep):
for j in range(0, img.shape[1] - Ni, jstep):
patch = img[i:i + Ni, j:j + Nj]
if scale != 1:
patch = transform.resize(patch, patch_size)
yield (i, j), patch
indices, patches = zip(*sliding_window(test_image))
patches_hog = np.array([feature.hog(patch) for patch in patches])
patches_hog.shape
Out[15]: (1911, 1215)
最后，我们可以取这些具有 HOG 特征的补丁，并使用我们的模型评估每个补丁是否包含面部：
In [16]: labels = model.predict(patches_hog)
labels.sum()
Out[16]: 48.0
我们看到在将近 2,000 个补丁中，我们发现了 48 个检测结果。让我们利用这些补丁的信息来显示它们在我们的测试图像中的位置，将它们绘制成矩形（参见 图 50-4）。
In [17]: fig, ax = plt.subplots()
ax.imshow(test_image, cmap='gray')
ax.axis('off')
Ni, Nj = positive_patches[0].shape
indices = np.array(indices)
for i, j in indices[labels == 1]:
ax.add_patch(plt.Rectangle((j, i), Nj, Ni, edgecolor='red',
alpha=0.3, lw=2, facecolor='none'))
图 50-4. 被确定包含面部的窗口。
所有检测到的补丁都重叠并找到了图像中的面部！对于几行 Python 代码来说效果不错。
注意事项和改进
如果你深入研究前面的代码和示例，你会发现在我们宣称拥有一个可以投入生产的人脸检测器之前，我们还有一些工作要做。我们做的工作存在几个问题，也可以做出几个改进。特别是：
我们的训练集，特别是负面特征，不太完整。
中心问题在于训练集中有许多像面部的纹理，而我们当前的模型非常容易产生假阳性。如果你尝试在完整的宇航员图像上运行算法，就会发现这一点：当前的模型在图像的其他区域导致了许多误检测。
我们可以通过向负训练集添加更多类型的图像来解决这个问题，这可能会带来一些改善。另一种选择是使用更有针对性的方法，如硬负样本挖掘，我们采用一组分类器尚未见过的新图像，找到所有表示假阳性的补丁，并明确将它们作为负实例添加到训练集中，然后重新训练分类器。
我们当前的流水线仅在一个尺度上进行搜索。
按照当前的写法，我们的算法会漏掉不是大约 62 × 47 像素的面部。可以通过使用各种大小的滑动窗口，并在将每个补丁输入模型之前使用skimage.transform.resize来直接解决这个问题。事实上，这里使用的sliding_window工具已经考虑到了这一点。
我们应该结合重叠的检测补丁。
对于一个可投入生产的流水线，我们更希望不要有 30 个相同人脸的检测结果，而是将重叠的检测组减少到一个单独的检测结果。这可以通过无监督的聚类方法（均值漂移聚类是一个很好的选择），或者通过类似于非最大抑制这样的程序化方法来实现，这是机器视觉中常见的一种算法。
流水线应该被简化
一旦解决了上述问题，创建一个更简化的流水线来输入训练图像并预测滑动窗口输出也将是一个不错的选择。这就是 Python 作为数据科学工具的优势所在：通过一点工作，我们可以将我们的原型代码打包成一个设计良好的面向对象的 API，让用户能够轻松地使用它。我将把这留给读者作为一个“练习题”。
更近期的进展：深度学习
最后，我应该补充说，在机器学习环境中，HOG 和其他程序化特征提取方法并不总是被使用。相反，许多现代目标检测流水线使用深度神经网络的变体（通常称为深度学习）：一个思考神经网络的方式是将其视为从数据中确定最佳特征提取策略的估计器，而不是依赖用户的直觉。
虽然近年来该领域取得了巨大的成果，但深度学习在概念上与前几章中探讨的机器学习模型并没有太大的不同。主要进步在于利用现代计算硬件（通常是大型强大机器集群）在更大的训练数据集上训练更加灵活的模型。但尽管规模不同，最终目标仍然是非常相似的：从数据中构建模型。
如果你对更深入的了解感兴趣，以下部分的参考文献清单应该是一个很好的起点！
更多机器学习资源
本书的这一部分快速介绍了 Python 中的机器学习，主要使用了 Scikit-Learn 库中的工具。尽管这些章节很长，但仍然太短，无法涵盖许多有趣和重要的算法、方法和讨论。在这里，我想为那些有兴趣的人提供一些关于在 Python 中学习更多关于机器学习的资源：
Scikit-Learn 网站拥有令人印象深刻的文档和示例，涵盖了这里讨论的一些模型，以及更多内容。如果你想要简要了解最重要和经常使用的机器学习算法，这是一个很好的开始。
SciPy、PyCon 和 PyData 教程视频
Scikit-Learn和其他机器学习主题是许多以 Python 为重点的会议系列（特别是 PyCon、SciPy 和 PyData 会议）教程轨道中的常青之选。这些会议大多免费在线发布其主题演讲、讨论和教程的视频，你可以通过合适的网络搜索轻松找到（例如，“PyCon 2022 视频”）。
《Python 机器学习导论》，作者 Andreas C. Müller 和 Sarah Guido（O’Reilly）
本书涵盖了这些章节讨论的许多机器学习基础知识，但特别是在涵盖 Scikit-Learn 更高级特性方面更具相关性，包括额外的估算器、模型验证方法和管道化。
《使用 PyTorch 和 Scikit-Learn 的机器学习》，作者 Sebastian Raschka（Packt）
Sebastian Raschka 的最新书籍从这些章节覆盖的基础主题开始，但深入探讨，并展示了这些概念如何应用于更复杂和计算密集的深度学习和强化学习模型，使用知名的PyTorch 库。
