Python 数据科学手册第二版（七）
原文：
zh.annas-archive.org/md5/051facaf2908ae8198253e3a14b09ec1译者：飞龙
第四十一章：深入：朴素贝叶斯分类
前四章概述了机器学习的概念。在第五部分的其余部分，我们将首先更详细地查看四种监督学习算法，然后是四种无监督学习算法。我们从第一个监督方法朴素贝叶斯分类开始。
朴素贝叶斯模型是一组极快速且简单的分类算法，通常适用于非常高维度的数据集。因为它们速度快、可调参数少，所以它们通常用作分类问题的快速基准线。本章将提供朴素贝叶斯分类器工作原理的直观解释，并在一些数据集上展示它们的几个例子。
贝叶斯分类
朴素贝叶斯分类器是建立在贝叶斯分类方法之上的。这些方法依赖于贝叶斯定理，该定理描述了统计量的条件概率关系。在贝叶斯分类中，我们感兴趣的是找到给定一些观察特征的标签L的概率，可以写作P ( L | features )。贝叶斯定理告诉我们如何用我们可以更直接计算的量来表达这一点：
P ( L | features ) = P( features |L)P(L) P( features )
如果我们试图在两个标签之间做出决策——让我们称它们为L 1和L 2——那么做出这个决定的一种方法是计算每个标签的后验概率的比率：
P(L 1 | features ) P(L 2 | features ) = P( features |L 1 ) P( features |L 2 ) P(L 1 ) P(L 2 )
现在我们所需的只是一些模型，通过这些模型我们可以计算每个标签P ( features | L i )。这样的模型被称为生成模型，因为它指定了生成数据的假设随机过程。为每个标签指定这种生成模型是这样一个贝叶斯分类器训练的主要部分。对于这样一个训练步骤的一般版本来说，这是一个非常困难的任务，但是我们可以通过对这个模型形式做一些简化的假设来简化它。
这就是“朴素贝叶斯”中的“朴素”所在：如果我们对每个标签的生成模型作出非常朴素的假设，我们可以找到每个类别的生成模型的粗略近似，然后继续贝叶斯分类。不同类型的朴素贝叶斯分类器基于关于数据的不同朴素假设，我们将在接下来的几节中讨论其中一些。
我们从标准导入开始：
In [1]: %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('seaborn-whitegrid')
高斯朴素贝叶斯
或许最容易理解的朴素贝叶斯分类器是高斯朴素贝叶斯。使用这个分类器，假设每个标签的数据都来自简单的高斯分布。想象一下我们有以下数据，显示在图 41-1 中：
In [2]: from sklearn.datasets import make_blobs
X, y = make_blobs(100, 2, centers=2, random_state=2, cluster_std=1.5)
plt.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='RdBu');
图 41-1. 高斯朴素贝叶斯分类数据^(1)
最简单的高斯模型假设数据由没有各维度间协方差的高斯分布描述。这个模型可以通过计算每个标签内点的均值和标准差来拟合，这是我们定义这种分布所需的所有内容。这种朴素高斯假设的结果显示在图 41-2 中。
图 41-2. 高斯朴素贝叶斯模型可视化^(2)
这里的椭圆代表每个标签的高斯生成模型，中心区域的概率更高。有了每个类别的生成模型，我们可以简单地计算任何数据点的似然P ( features | L 1 )，因此我们可以快速计算后验比率，并确定给定点最有可能的标签。
这个过程在 Scikit-Learn 的sklearn.naive_bayes.GaussianNB估计器中实现：
In [3]: from sklearn.naive_bayes import GaussianNB
model = GaussianNB()
model.fit(X, y);
让我们生成一些新数据并预测标签：
In [4]: rng = np.random.RandomState(0)
Xnew = [-6, -14] + [14, 18] * rng.rand(2000, 2)
ynew = model.predict(Xnew)
现在我们可以绘制这些新数据，以了解决策边界的位置（见图 41-3）。
In [5]: plt.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='RdBu')
lim = plt.axis()
plt.scatter(Xnew[:, 0], Xnew[:, 1], c=ynew, s=20, cmap='RdBu', alpha=0.1)
plt.axis(lim);
图 41-3. 高斯朴素贝叶斯分类可视化
我们可以看到分类中有一个略微弯曲的边界——一般来说，高斯朴素贝叶斯模型产生的边界将是二次的。
这种贝叶斯形式主义的一个优点是它自然地允许概率分类，我们可以使用predict_proba方法来计算：
In [6]: yprob = model.predict_proba(Xnew)
yprob[-8:].round(2)
Out[6]: array([[0.89, 0.11],
[1. , 0. ],
[1. , 0. ],
[1. , 0. ],
[1. , 0. ],
[1. , 0. ],
[0. , 1. ],
[0.15, 0.85]])
这些列分别给出了第一个和第二个标签的后验概率。如果您正在寻找分类中的不确定性估计，像这样的贝叶斯方法可能是一个很好的起点。
当然，最终的分类结果将仅仅与导致它的模型假设一样好，这就是为什么高斯朴素贝叶斯通常不会产生非常好的结果的原因。尽管如此，在许多情况下——特别是特征数量变得很大时——这种假设并不足以阻止高斯朴素贝叶斯成为一种可靠的方法。
多项式朴素贝叶斯
刚刚描述的高斯假设远非是可以用于指定每个标签生成分布的唯一简单假设。另一个有用的例子是多项式朴素贝叶斯，其中假设特征是从简单的多项分布生成的。多项分布描述了在多个类别中观察计数的概率，因此多项式朴素贝叶斯最适合表示计数或计数率的特征。
这个想法与之前完全相同，只是不再用最佳拟合的高斯来建模数据分布，而是用最佳拟合的多项分布来建模。
示例：文本分类
多项式朴素贝叶斯经常用于文本分类的一个场合，其中特征与要分类的文档中的单词计数或频率相关。我们在第四十章中讨论了从文本中提取这些特征；在这里，我们将使用通过 Scikit-Learn 提供的 20 个新闻组语料库中的稀疏单词计数特征来展示如何将这些简短文档分类到不同的类别中。
让我们下载数据并查看目标名称：
In [7]: from sklearn.datasets import fetch_20newsgroups
data = fetch_20newsgroups()
data.target_names
Out[7]: ['alt.atheism',
'comp.graphics',
'comp.os.ms-windows.misc',
'comp.sys.ibm.pc.hardware',
'comp.sys.mac.hardware',
'comp.windows.x',
'misc.forsale',
'rec.autos',
'rec.motorcycles',
'rec.sport.baseball',
'rec.sport.hockey',
'sci.crypt',
'sci.electronics',
'sci.med',
'sci.space',
'soc.religion.christian',
'talk.politics.guns',
'talk.politics.mideast',
'talk.politics.misc',
'talk.religion.misc']
为了简单起见，在这里我们将只选择其中几个类别并下载训练和测试集：
In [8]: categories = ['talk.religion.misc', 'soc.religion.christian',
'sci.space', 'comp.graphics']
train = fetch_20newsgroups(subset='train', categories=categories)
test = fetch_20newsgroups(subset='test', categories=categories)
这里是数据的一个代表性条目：
In [9]: print(train.data[5][48:])
Out[9]: Subject: Federal Hearing
Originator: dmcgee@uluhe
Organization: School of Ocean and Earth Science and Technology
Distribution: usa
Lines: 10
Fact or rumor....? Madalyn Murray O'Hare an atheist who eliminated the
use of the bible reading and prayer in public schools 15 years ago is now
going to appear before the FCC with a petition to stop the reading of the
Gospel on the airways of America. And she is also campaigning to remove
Christmas programs, songs, etc from the public schools. If it is true
then mail to Federal Communications Commission 1919 H Street Washington DC
20054 expressing your opposition to her request. Reference Petition number
2493.
为了将这些数据用于机器学习，我们需要将每个字符串的内容转换为一个数字向量。为此，我们将使用 TF-IDF 向量化器（在第四十章介绍），并创建一个管道，将其附加到多项式朴素贝叶斯分类器：
In [10]: from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.naive_bayes import MultinomialNB
from sklearn.pipeline import make_pipeline
model = make_pipeline(TfidfVectorizer(), MultinomialNB())
有了这个管道，我们可以将模型应用于训练数据，并预测测试数据的标签：
In [11]: model.fit(train.data, train.target)
labels = model.predict(test.data)
现在我们已经预测了测试数据的标签，我们可以评估它们以了解估计器的性能。例如，让我们看一下测试数据的真实标签和预测标签之间的混淆矩阵（参见图 41-4）。
In [12]: from sklearn.metrics import confusion_matrix
mat = confusion_matrix(test.target, labels)
sns.heatmap(mat.T, square=True, annot=True, fmt='d', cbar=False,
xticklabels=train.target_names, yticklabels=train.target_names,
cmap='Blues')
plt.xlabel('true label')
plt.ylabel('predicted label');
显然，即使这个非常简单的分类器可以成功地将关于空间讨论与计算机讨论分开，但它会在宗教讨论和基督教讨论之间感到困惑。这或许是可以预料的！
这里很酷的一点是，我们现在有工具来确定任何字符串的类别，只需使用此管道的predict方法。下面是一个实用函数，用于返回单个字符串的预测结果：
In [13]: def predict_category(s, train=train, model=model):
pred = model.predict([s])
return train.target_names[pred[0]]
让我们试试它：
In [14]: predict_category('sending a payload to the ISS')
Out[14]: 'sci.space'
In [15]: predict_category('discussing the existence of God')
Out[15]: 'soc.religion.christian'
In [16]: predict_category('determining the screen resolution')
Out[16]: 'comp.graphics'
图 41-4. 多项式朴素贝叶斯文本分类器的混淆矩阵
请记住，这只是一个对字符串中每个单词（加权）频率的简单概率模型；尽管如此，结果令人印象深刻。即使是非常朴素的算法，在小心使用并在大量高维数据上训练时，也可以出奇地有效。
何时使用朴素贝叶斯
由于朴素贝叶斯分类器对数据做出如此严格的假设，它们通常不如更复杂的模型表现好。尽管如此，它们有几个优点：
-
它们在训练和预测时都非常快速。
-
它们提供直观的概率预测。
-
它们通常易于解释。
-
它们具有少量（如果有的话）可调参数。
这些优势意味着朴素贝叶斯分类器通常是作为初始基线分类的不错选择。如果它表现得合适，那么恭喜你：你已经拥有了一个非常快速、易于解释的分类器来解决你的问题。如果它表现不佳，那么你可以开始探索更复杂的模型，同时具备一些关于它们应该如何表现的基础知识。
朴素贝叶斯分类器在以下情况下表现特别好：
-
当朴素假设实际上与数据匹配时（在实践中非常罕见）
-
对于非常分离的类别，当模型复杂度不那么重要时
-
对于非常高维数据，当模型复杂度不那么重要时
最后两点看似不同，但实际上是相关的：随着数据集维度的增长，任何两个点在一起的可能性大大降低（毕竟，它们必须在每个维度上都很接近才能在总体上接近）。这意味着在高维空间中，簇通常比低维空间中更为分离，平均而言。基于这个原因，像这里讨论的简单分类器往往在维度增加时表现得同样或更好：一旦你有足够的数据，即使是简单模型也可以非常强大。
^(1) 此图的全彩版本可在GitHub上找到。
^(2) 生成此图的代码可在在线附录中找到。
第四十二章：深入解析：线性回归
就像朴素贝叶斯（讨论见第四十一章）对于分类任务是一个很好的起点一样，线性回归模型对于回归任务也是一个很好的起点。这样的模型很受欢迎，因为它们可以快速拟合并且易于解释。你已经熟悉了最简单形式的线性回归模型（即将直线拟合到二维数据），但是这样的模型可以扩展到对更复杂的数据行为进行建模。
在本章中，我们将首先快速了解这个众所周知问题背后的数学知识，然后再看看线性模型如何被泛化以解决数据中更复杂的模式。
我们从标准导入开始：
In [1]: %matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
import numpy as np
简单线性回归
我们将从最熟悉的线性回归开始，即对数据进行直线拟合。直线拟合是一个形式为的模型：
y = a x + b
其中a通常被称为斜率，而b通常被称为截距。
考虑以下数据，这些数据分布在一条斜率为 2，截距为-5 的直线周围（见图 42-1）。
In [2]: rng = np.random.RandomState(1)
x = 10 * rng.rand(50)
y = 2 * x - 5 + rng.randn(50)
plt.scatter(x, y);
图 42-1. 线性回归的数据
我们可以使用 Scikit-Learn 的LinearRegression估计器来拟合这些数据并构建最佳拟合线，如图 42-2 所示。
In [3]: from sklearn.linear_model import LinearRegression
model = LinearRegression(fit_intercept=True)
model.fit(x[:, np.newaxis], y)
xfit = np.linspace(0, 10, 1000)
yfit = model.predict(xfit[:, np.newaxis])
plt.scatter(x, y)
plt.plot(xfit, yfit);
图 42-2. 一个简单的线性回归模型
数据的斜率和截距包含在模型的拟合参数中，在 Scikit-Learn 中始终以下划线结尾标记。这里的相关参数是coef_和intercept_：
In [4]: print("Model slope: ", model.coef_[0])
print("Model intercept:", model.intercept_)
Out[4]: Model slope: 2.0272088103606953
Model intercept: -4.998577085553204
我们看到结果与用于生成数据的值非常接近，这是我们所希望的。
然而，LinearRegression估计器要比这更强大——除了简单的直线拟合外，它还可以处理形式为的多维线性模型：
y = a 0 + a 1 x 1 + a 2 x 2 + ⋯
其中有多个x值。从几何上讲，这类似于在三维空间中对点拟合平面，或者在更高维度中对点拟合超平面。
这种回归的多维性使其更难以可视化，但我们可以通过构建一些示例数据，使用 NumPy 的矩阵乘法运算符来看到其中一个拟合的过程：
In [5]: rng = np.random.RandomState(1)
X = 10 * rng.rand(100, 3)
y = 0.5 + np.dot(X, [1.5, -2., 1.])
model.fit(X, y)
print(model.intercept_)
print(model.coef_)
Out[5]: 0.50000000000001
[ 1.5 -2. 1. ]
这里的y数据是从三个随机x值的线性组合构成的，线性回归恢复了用于构建数据的系数。
通过这种方式，我们可以使用单个 LinearRegression 评估器来拟合线条、平面或超平面到我们的数据。看起来这种方法仍然限制在变量之间严格的线性关系，但事实证明我们也可以放宽这一点。
基函数回归
你可以用的一个技巧是将线性回归适应变量之间的非线性关系，方法是根据 基函数 转换数据。我们之前已经见过这样的一个版本，在第三十九章和第四十章中使用的 PolynomialRegression 流水线中。这个想法是将我们的多维线性模型：
y = a 0 + a 1 x 1 + a 2 x 2 + a 3 x 3 + ⋯
并从我们的单维输入 x 中构建 x 1 , x 2 , x 3 , 等等。也就是说，我们让 x n = f n ( x ) ，其中 f n ( ) 是将我们的数据转换的某个函数。
例如，如果 f n ( x ) = x n ，我们的模型就会变成多项式回归：
y = a 0 + a 1 x + a 2 x 2 + a 3 x 3 + ⋯
注意，这仍然是 线性模型 —— 线性指的是系数 a n 从不相乘或相除。我们所做的实质上是将我们的一维 x 值投影到更高的维度，这样线性拟合可以拟合 x 和 y 之间更复杂的关系。
多项式基函数
这种多项式投影非常有用，以至于它被内置到 Scikit-Learn 中，使用 PolynomialFeatures 变换器：
In [6]: from sklearn.preprocessing import PolynomialFeatures
x = np.array([2, 3, 4])
poly = PolynomialFeatures(3, include_bias=False)
poly.fit_transform(x[:, None])
Out[6]: array([[ 2., 4., 8.],
[ 3., 9., 27.],
[ 4., 16., 64.]])
我们在这里看到，变换器已经将我们的一维数组转换为一个三维数组，其中每列包含了指数化的值。这种新的、更高维的数据表示可以被插入到线性回归中。
正如我们在第四十章中看到的，实现这一点的最干净的方法是使用一个流水线。让我们用这种方式制作一个 7 次多项式模型：
In [7]: from sklearn.pipeline import make_pipeline
poly_model = make_pipeline(PolynomialFeatures(7),
LinearRegression())
通过这种转换，我们可以使用线性模型更好地拟合x和y之间更复杂的关系。例如，这里是带噪声的正弦波（参见图 42-3）。
In [8]: rng = np.random.RandomState(1)
x = 10 * rng.rand(50)
y = np.sin(x) + 0.1 * rng.randn(50)
poly_model.fit(x[:, np.newaxis], y)
yfit = poly_model.predict(xfit[:, np.newaxis])
plt.scatter(x, y)
plt.plot(xfit, yfit);
图 42-3. 对非线性训练数据进行线性多项式拟合
我们的线性模型，通过使用七阶多项式基函数，可以很好地拟合这些非线性数据！
高斯基函数
当然，也可以使用其他基函数。例如，一个有用的模式是拟合一个不是多项式基函数的模型，而是高斯基函数的总和。结果可能看起来像是图 42-4。
图 42-4. 高斯基函数拟合非线性数据^(1)
图中阴影区域是经过缩放的基函数，将它们相加后可以重现数据中的平滑曲线。这些高斯基函数没有内置到 Scikit-Learn 中，但我们可以编写一个自定义转换器来创建它们，如此处所示，并在图 42-5 中进行了说明（Scikit-Learn 的转换器是以 Python 类的形式实现的；查看 Scikit-Learn 的源代码是了解它们如何创建的好方法）：
In [9]: from sklearn.base import BaseEstimator, TransformerMixin
class GaussianFeatures(BaseEstimator, TransformerMixin):
"""Uniformly spaced Gaussian features for one-dimensional input"""
def __init__(self, N, width_factor=2.0):
self.N = N
self.width_factor = width_factor
@staticmethod
def _gauss_basis(x, y, width, axis=None):
arg = (x - y) / width
return np.exp(-0.5 * np.sum(arg ** 2, axis))
def fit(self, X, y=None):
# create N centers spread along the data range
self.centers_ = np.linspace(X.min(), X.max(), self.N)
self.width_ = self.width_factor*(self.centers_[1]-self.centers_[0])
return self
def transform(self, X):
return self._gauss_basis(X[:, :, np.newaxis], self.centers_,
self.width_, axis=1)
gauss_model = make_pipeline(GaussianFeatures(20),
LinearRegression())
gauss_model.fit(x[:, np.newaxis], y)
yfit = gauss_model.predict(xfit[:, np.newaxis])
plt.scatter(x, y)
plt.plot(xfit, yfit)
plt.xlim(0, 10);
图 42-5. 使用自定义转换器计算的高斯基函数拟合
我包含了这个例子，只是为了明确指出多项式基函数并非魔法：如果你对数据生成过程有某种直觉，认为某种基函数可能更合适，你可以使用它。
正则化
将基函数引入线性回归中使模型更加灵活，但也很快会导致过拟合（参见第三十九章中的讨论）。例如，如果使用大量高斯基函数，图 42-6 展示了会发生什么：
In [10]: model = make_pipeline(GaussianFeatures(30),
LinearRegression())
model.fit(x[:, np.newaxis], y)
plt.scatter(x, y)
plt.plot(xfit, model.predict(xfit[:, np.newaxis]))
plt.xlim(0, 10)
plt.ylim(-1.5, 1.5);
图 42-6. 过度复杂的基函数模型会过拟合数据
将数据投影到 30 维基函数后，模型的灵活性大大增加，并在数据约束位置之间出现极值。如果我们绘制高斯基函数系数相对于它们位置的图表，我们可以看到这一现象，如图 42-7 所示。
In [11]: def basis_plot(model, title=None):
fig, ax = plt.subplots(2, sharex=True)
model.fit(x[:, np.newaxis], y)
ax[0].scatter(x, y)
ax[0].plot(xfit, model.predict(xfit[:, np.newaxis]))
ax[0].set(xlabel='x', ylabel='y', ylim=(-1.5, 1.5))
if title:
ax[0].set_title(title)
ax[1].plot(model.steps[0][1].centers_,
model.steps[1][1].coef_)
ax[1].set(xlabel='basis location',
ylabel='coefficient',
xlim=(0, 10))
model = make_pipeline(GaussianFeatures(30), LinearRegression())
basis_plot(model)
图 42-7. 过度复杂模型中高斯基函数的系数
此图的下部面板显示了每个位置的基函数的幅度。当基函数重叠时，这是典型的过拟合行为：相邻基函数的系数会急剧增加并相互抵消。我们知道这种行为是有问题的，如果我们可以通过惩罚模型参数的大值来明确限制这样的峰值，那就太好了。这样的惩罚被称为正则化，有几种形式。
岭回归（L[2]正则化）
或许最常见的正则化形式被称为岭回归或L 2 正则化（有时也称为Tikhonov 正则化）。这通过对模型系数θ n的平方和（2-范数）进行惩罚来实现。在这种情况下，模型拟合的惩罚将是：
P = α ∑ n=1 N θ n 2
其中α是一个自由参数，用于控制惩罚的强度。这种类型的惩罚模型已经内置到 Scikit-Learn 中的Ridge估计器中（参见 Figure 42-8）。
In [12]: from sklearn.linear_model import Ridge
model = make_pipeline(GaussianFeatures(30), Ridge(alpha=0.1))
basis_plot(model, title='Ridge Regression')
图 42-8. 岭（L[2]）正则化应用于过度复杂的模型（与 Figure 42-7 进行比较）
参数α本质上是一个控制生成模型复杂性的旋钮。在极限α → 0中，我们恢复了标准线性回归结果；在极限α → ∞中，所有模型响应都将被抑制。岭回归的一个优点是它特别高效地计算—几乎没有比原始线性回归模型更多的计算成本。
套索回归（L[1]正则化）
另一种常见的正则化方法被称为套索回归或L[1]正则化，它涉及对回归系数的绝对值（1-范数）的惩罚：
P = α ∑ n=1 N | θ n |
尽管这在概念上与岭回归非常相似，但结果可能出奇地不同。例如，由于其构造，套索回归倾向于偏爱可能的稀疏模型：也就是说，它更倾向于将许多模型系数设为零。
如果我们使用L 1 -归一化系数复制前面的示例，我们可以看到这种行为（参见 Figure 42-9）。
In [13]: from sklearn.linear_model import Lasso
model = make_pipeline(GaussianFeatures(30),
Lasso(alpha=0.001, max_iter=2000))
basis_plot(model, title='Lasso Regression')
图 42-9. 套索（L[1]）正则化应用于过度复杂的模型（与 Figure 42-8 进行比较）
使用套索回归惩罚，大多数系数确实为零，功能行为由可用基函数的一小部分建模。与岭回归一样，α参数调节惩罚的强度，应通过例如交叉验证确定（请参阅第三十九章讨论此问题）。
示例：预测自行车流量
举个例子，让我们看看是否能够预测西雅图弗里蒙特桥上的自行车出行次数，基于天气、季节和其他因素。我们在第二十三章已经看过这些数据，但在这里我们将自行车数据与另一个数据集结合，并尝试确定天气和季节因素——温度、降水和日照小时——对该走廊自行车流量的影响程度。幸运的是，美国国家海洋和大气管理局（NOAA）提供其日常气象站数据——我使用的是站点 ID USW00024233——我们可以轻松使用 Pandas 将这两个数据源连接起来。我们将执行简单的线性回归，将天气和其他信息与自行车计数关联起来，以估算这些参数中的任何变化如何影响给定日的骑行者数量。
特别是，这是 Scikit-Learn 工具如何在统计建模框架中使用的示例，其中假定模型的参数具有可解释的含义。正如前面讨论的那样，这不是机器学习中的标准方法，但对于某些模型是可能的。
让我们从加载两个数据集开始，以日期为索引：
In [14]: # url = 'https://raw.githubusercontent.com/jakevdp/bicycle-data/main'
# !curl -O {url}/FremontBridge.csv
# !curl -O {url}/SeattleWeather.csv
In [15]: import pandas as pd
counts = pd.read_csv('FremontBridge.csv',
index_col='Date', parse_dates=True)
weather = pd.read_csv('SeattleWeather.csv',
index_col='DATE', parse_dates=True)
为简单起见，让我们查看 2020 年之前的数据，以避免新冠肺炎大流行的影响，这显著影响了西雅图的通勤模式：
In [16]: counts = counts[counts.index < "2020-01-01"]
weather = weather[weather.index < "2020-01-01"]
接下来我们将计算每日自行车总流量，并将其放入独立的DataFrame中：
In [17]: daily = counts.resample('d').sum()
daily['Total'] = daily.sum(axis=1)
daily = daily[['Total']] # remove other columns
我们之前看到使用模式通常从一天到另一天有所不同。让我们在数据中考虑这一点，通过添加指示星期几的二进制列：
In [18]: days = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']
for i in range(7):
daily[days[i]] = (daily.index.dayofweek == i).astype(float)
同样，我们可能期待骑行者在假日有不同的行为表现；让我们也加上一个指示器：
In [19]: from pandas.tseries.holiday import USFederalHolidayCalendar
cal = USFederalHolidayCalendar()
holidays = cal.holidays('2012', '2020')
daily = daily.join(pd.Series(1, index=holidays, name='holiday'))
daily['holiday'].fillna(0, inplace=True)
我们还可能怀疑白天的小时数会影响骑行的人数。让我们使用标准的天文计算来添加这些信息（参见图 42-10）。
In [20]: def hours_of_daylight(date, axis=23.44, latitude=47.61):
"""Compute the hours of daylight for the given date"""
days = (date - pd.datetime(2000, 12, 21)).days
m = (1. - np.tan(np.radians(latitude))
* np.tan(np.radians(axis) * np.cos(days * 2 * np.pi / 365.25)))
return 24. * np.degrees(np.arccos(1 - np.clip(m, 0, 2))) / 180.
daily['daylight_hrs'] = list(map(hours_of_daylight, daily.index))
daily[['daylight_hrs']].plot()
plt.ylim(8, 17)
Out[20]: (8.0, 17.0)
图 42-10. 西雅图的日照小时可视化
我们还可以将平均温度和总降水量添加到数据中。除了降水英寸外，让我们添加一个指示某一天是否干燥（降水量为零）的标志：
In [21]: weather['Temp (F)'] = 0.5 * (weather['TMIN'] + weather['TMAX'])
weather['Rainfall (in)'] = weather['PRCP']
weather['dry day'] = (weather['PRCP'] == 0).astype(int)
daily = daily.join(weather[['Rainfall (in)', 'Temp (F)', 'dry day']])
最后，让我们添加一个从第 1 天开始递增的计数器，并测量经过了多少年。这将让我们测量每日过境量的观察到的年增长或年减少：
In [22]: daily['annual'] = (daily.index - daily.index[0]).days / 365.
现在我们的数据已经整理好了，我们可以看一下：
In [23]: daily.head()
Out[23]: Total Mon Tue Wed Thu Fri Sat Sun holiday \
Date
2012-10-03 14084.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
2012-10-04 13900.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
2012-10-05 12592.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
2012-10-06 8024.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
2012-10-07 8568.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
daylight_hrs Rainfall (in) Temp (F) dry day annual
Date
2012-10-03 11.277359 0.0 56.0 1 0.000000
2012-10-04 11.219142 0.0 56.5 1 0.002740
2012-10-05 11.161038 0.0 59.5 1 0.005479
2012-10-06 11.103056 0.0 60.5 1 0.008219
2012-10-07 11.045208 0.0 60.5 1 0.010959
有了这些东西，我们可以选择要使用的列，并对我们的数据拟合一个线性回归模型。我们将设置fit_intercept=False，因为每日标志基本上充当它们自己的特定于每天的截距：
In [24]: # Drop any rows with null values
daily.dropna(axis=0, how='any', inplace=True)
column_names = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun',
'holiday', 'daylight_hrs', 'Rainfall (in)',
'dry day', 'Temp (F)', 'annual']
X = daily[column_names]
y = daily['Total']
model = LinearRegression(fit_intercept=False)
model.fit(X, y)
daily['predicted'] = model.predict(X)
最后，我们可以通过视觉比较总的和预测的自行车流量（见图 42-11）。
In [25]: daily[['Total', 'predicted']].plot(alpha=0.5);
图 42-11. 我们模型对自行车流量的预测
从数据和模型预测不完全一致的事实来看，很明显我们已经错过了一些关键特征。我们的特征要么不完整（即，人们决定是否骑自行车去工作基于不止这些特征），要么有一些非线性关系我们没有考虑到（例如，也许在高温和低温下人们骑行较少）。然而，我们的粗略近似足以给我们一些见解，我们可以查看线性模型的系数以估算每个特征对每日自行车数量的贡献程度：
In [26]: params = pd.Series(model.coef_, index=X.columns)
params
Out[26]: Mon -3309.953439
Tue -2860.625060
Wed -2962.889892
Thu -3480.656444
Fri -4836.064503
Sat -10436.802843
Sun -10795.195718
holiday -5006.995232
daylight_hrs 409.146368
Rainfall (in) -2789.860745
dry day 2111.069565
Temp (F) 179.026296
annual 324.437749
dtype: float64
这些数字如果没有一些不确定性的度量，就很难解释。我们可以通过对数据进行自举重采样来快速计算这些不确定性：
In [27]: from sklearn.utils import resample
np.random.seed(1)
err = np.std([model.fit(*resample(X, y)).coef_
for i in range(1000)], 0)
有了这些误差估计，让我们再次查看结果：
In [28]: print(pd.DataFrame({'effect': params.round(0),
'uncertainty': err.round(0)}))
Out[28]: effect uncertainty
Mon -3310.0 265.0
Tue -2861.0 274.0
Wed -2963.0 268.0
Thu -3481.0 268.0
Fri -4836.0 261.0
Sat -10437.0 259.0
Sun -10795.0 267.0
holiday -5007.0 401.0
daylight_hrs 409.0 26.0
Rainfall (in) -2790.0 186.0
dry day 2111.0 101.0
Temp (F) 179.0 7.0
annual 324.0 22.0
这里的effect列，粗略地说，显示了骑手数量如何受到所讨论特征变化的影响。例如，一周中的某一天就有明显的区别：周末的骑手要比工作日少几千人。我们还看到，每增加一个小时的阳光，会有 409 ± 26 人选择骑行；华氏度每增加一度，就会有 179 ± 7 人选择骑自行车；晴天意味着平均增加 2,111 ± 101 名骑手，而每英寸的降雨则导致 2,790 ± 186 名骑手选择另一种交通方式。一旦考虑了所有这些影响，我们就会看到每年新的日常骑手数量有一个适度的增加，为 324 ± 22 人。
我们的简单模型几乎肯定会缺少一些相关信息。例如，正如前面提到的，非线性效应（如降水和寒冷温度的影响）以及每个变量内的非线性趋势（例如对极冷和极热温度不愿意骑车的影响）无法在简单的线性模型中考虑。此外，我们还丢弃了一些更精细的信息（例如雨天早晨和雨天下午之间的差异），并且忽略了天之间的相关性（例如一个雨天星期二对星期三的可能影响，或者在连续多日雨天后出现意外晴天的影响）。这些都是潜在的有趣效应，如果你愿意，现在你已经有了开始探索它们的工具！
^(1) 生成这个图表的代码可以在在线附录中找到。
第四十三章：深入探讨支持向量机
支持向量机（SVMs）是一种特别强大和灵活的监督算法类，适用于分类和回归。在本章中，我们将探讨 SVM 背后的直觉及其在分类问题中的应用。
我们从标准导入开始：
In [1]: %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
from scipy import stats
注
全尺寸、全彩色图像可在GitHub 的补充材料中找到。
激励支持向量机
作为我们讨论贝叶斯分类的一部分（参见第四十一章），我们了解到了描述每个潜在类分布的简单模型，并尝试使用它来概率地确定新点的标签。那是一个生成分类的例子；在这里，我们将考虑判别分类。也就是说，我们不再模拟每个类，而是简单地找到一个（在二维中为线或曲线，在多维中为流形），将类彼此分开。
作为一个例子，考虑一个分类任务的简单情况，其中两类点是完全分开的（见图 43-1）。
In [2]: from sklearn.datasets import make_blobs
X, y = make_blobs(n_samples=50, centers=2,
random_state=0, cluster_std=0.60)
plt.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='autumn');
图 43-1. 分类简单数据
线性判别分类器将尝试画一条直线分隔两组数据，并因此创建一个分类模型。对于像这样的二维数据，我们可以手动完成这项任务。但我们立即看到了一个问题：存在不止一条可能完全区分这两类的分界线！
我们可以如下绘制其中一些；图 43-2 展示了结果：
In [3]: xfit = np.linspace(-1, 3.5)
plt.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='autumn')
plt.plot([0.6], [2.1], 'x', color='red', markeredgewidth=2, markersize=10)
for m, b in [(1, 0.65), (0.5, 1.6), (-0.2, 2.9)]:
plt.plot(xfit, m * xfit + b, '-k')
plt.xlim(-1, 3.5);
图 43-2. 我们数据的三个完美线性判别分类器
这是三个非常不同的分隔符，尽管如此，它们完全可以区分这些样本。根据你选择的分界线，新的数据点（例如，在这个图中用“X”标记的点）将被分配不同的标签！显然，我们简单的“在类之间画一条线”的直觉不够好，我们需要更深入地思考。
支持向量机：最大化边缘
支持向量机提供了一种改进方法。其直觉是：与其简单地在类之间画一条零宽度线，我们可以在每条线周围绘制一定宽度的边缘，直到最近的点。这是一个展示的例子（见图 43-3）。
In [4]: xfit = np.linspace(-1, 3.5)
plt.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='autumn')
for m, b, d in [(1, 0.65, 0.33), (0.5, 1.6, 0.55), (-0.2, 2.9, 0.2)]:
yfit = m * xfit + b
plt.plot(xfit, yfit, '-k')
plt.fill_between(xfit, yfit - d, yfit + d, edgecolor='none',
color='lightgray', alpha=0.5)
plt.xlim(-1, 3.5);
最大化这个边缘的线就是我们将选择作为最优模型的线。
图 43-3. 判别分类器内“边缘”的可视化
适配支持向量机
现在让我们看看对这些数据进行实际拟合的结果：我们将使用 Scikit-Learn 的支持向量分类器（SVC）来训练一个 SVM 模型。暂时地，我们将使用线性核并将参数C设置为一个非常大的数（稍后我们将深入讨论它们的含义）：
In [5]: from sklearn.svm import SVC # "Support vector classifier"
model = SVC(kernel='linear', C=1E10)
model.fit(X, y)
Out[5]: SVC(C=10000000000.0, kernel='linear')
为了更好地可视化这里发生的情况，让我们创建一个快速便利函数，它将为我们绘制 SVM 决策边界（图 43-4）。
In [6]: def plot_svc_decision_function(model, ax=None, plot_support=True):
"""Plot the decision function for a 2D SVC"""
if ax is None:
ax = plt.gca()
xlim = ax.get_xlim()
ylim = ax.get_ylim()
# create grid to evaluate model
x = np.linspace(xlim[0], xlim[1], 30)
y = np.linspace(ylim[0], ylim[1], 30)
Y, X = np.meshgrid(y, x)
xy = np.vstack([X.ravel(), Y.ravel()]).T
P = model.decision_function(xy).reshape(X.shape)
# plot decision boundary and margins
ax.contour(X, Y, P, colors='k',
levels=[-1, 0, 1], alpha=0.5,
linestyles=['--', '-', '--'])
# plot support vectors
if plot_support:
ax.scatter(model.support_vectors_[:, 0],
model.support_vectors_[:, 1],
s=300, linewidth=1, edgecolors='black',
facecolors='none');
ax.set_xlim(xlim)
ax.set_ylim(ylim)
In [7]: plt.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='autumn')
plot_svc_decision_function(model);
图 43-4. 一个拟合到数据的支持向量机分类器，显示了边界（虚线）和支持向量（圆点）
这是最大化两组点之间间隔的分隔线。请注意，一些训练点恰好接触边界：它们在图 43-5 中被圈出来。这些点是此拟合的关键元素；它们被称为支持向量，并赋予了算法其名称。在 Scikit-Learn 中，这些点的标识存储在分类器的support_vectors_属性中：
In [8]: model.support_vectors_
Out[8]: array([[0.44359863, 3.11530945],
[2.33812285, 3.43116792],
[2.06156753, 1.96918596]])
此分类器成功的关键在于对拟合来说，只有支持向量的位置是重要的；远离边界但在正确一侧的点不会修改拟合。从技术上讲，这是因为这些点不会对用于拟合模型的损失函数产生贡献，因此它们的位置和数量并不重要，只要它们不跨越边界。
例如，如果我们绘制从这个数据集的前 60 个点和前 120 个点学习到的模型（图 43-5），我们可以看到这一点。
In [9]: def plot_svm(N=10, ax=None):
X, y = make_blobs(n_samples=200, centers=2,
random_state=0, cluster_std=0.60)
X = X[:N]
y = y[:N]
model = SVC(kernel='linear', C=1E10)
model.fit(X, y)
ax = ax or plt.gca()
ax.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='autumn')
ax.set_xlim(-1, 4)
ax.set_ylim(-1, 6)
plot_svc_decision_function(model, ax)
fig, ax = plt.subplots(1, 2, figsize=(16, 6))
fig.subplots_adjust(left=0.0625, right=0.95, wspace=0.1)
for axi, N in zip(ax, [60, 120]):
plot_svm(N, axi)
axi.set_title('N = {0}'.format(N))
图 43-5. 新训练点对 SVM 模型的影响
在左侧面板中，我们看到了 60 个训练点的模型和支持向量。在右侧面板中，我们增加了训练点的数量，但模型没有改变：左侧面板中的三个支持向量与右侧面板中的支持向量相同。这种对远点行为的确切不敏感性是 SVM 模型的一种优势之一。
如果您正在实时运行此笔记本，您可以使用 IPython 的交互式小部件来交互地查看 SVM 模型的此功能：
In [10]: from ipywidgets import interact, fixed
interact(plot_svm, N=(10, 200), ax=fixed(None));
Out[10]: interactive(children=(IntSlider(value=10, description='N', max=200, min=10),
> Output()), _dom_classes=('widget-...
超越线性边界：核支持向量机
当 SVM 与核结合时，它可以变得非常强大。我们之前在第四十二章中已经看到了核的一个版本，即基函数回归。在那里，我们将数据投影到由多项式和高斯基函数定义的更高维空间中，从而能够使用线性分类器拟合非线性关系。
在 SVM 模型中，我们可以使用相同思想的一个版本。为了证明核函数的必要性，让我们看一些不是线性可分的数据（参见图 43-6）。
In [11]: from sklearn.datasets import make_circles
X, y = make_circles(100, factor=.1, noise=.1)
clf = SVC(kernel='linear').fit(X, y)
plt.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='autumn')
plot_svc_decision_function(clf, plot_support=False);
图 43-6. 线性分类器对非线性边界的性能较差。
很明显，没有任何线性判别能永远分离这些数据。但我们可以从第四十二章的基函数回归中吸取教训，并思考如何将数据投影到更高的维度，以便线性分隔器足够。例如，我们可以使用的一个简单投影是在中间聚集上计算一个径向基函数（RBF）：
In [12]: r = np.exp(-(X ** 2).sum(1))
我们可以使用三维图来可视化这个额外的数据维度，如图 43-7 所示。
In [13]: from mpl_toolkits import mplot3d
ax = plt.subplot(projection='3d')
ax.scatter3D(X[:, 0], X[:, 1], r, c=y, s=50, cmap='autumn')
ax.view_init(elev=20, azim=30)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('r');
图 43-7. 为数据添加的第三个维度允许线性分离
我们可以看到，通过这个额外的维度，数据变得简单线性可分，通过在 r=0.7 处绘制一个分离平面。
在这种情况下，我们不得不选择并仔细调整我们的投影：如果我们没有将我们的径向基函数放在正确的位置，我们就不会看到如此清晰、线性可分的结果。一般来说，需要做出这样的选择是一个问题：我们希望以某种方式自动找到最佳的基函数来使用。
实现这一目标的一种策略是计算数据集中每个点处的基函数，并让 SVM 算法筛选结果。这种类型的基函数转换被称为核变换，因为它是基于每对点之间的相似关系（或核）。
这种策略的一个潜在问题是——将 N 点投影到 N 维空间中可能会变得非常计算密集，当 N 变大时。然而，由于一个称为 核技巧 的巧妙小程序，对核变换数据的拟合可以隐式完成——也就是说，根本不需要构建核投影的完整 N -维表示。这个核技巧内置在 SVM 中，是该方法如此强大的原因之一。
在 Scikit-Learn 中，我们可以通过将我们的线性核改为 RBF 核，使用 kernel 模型超参数来应用核化的 SVM：
In [14]: clf = SVC(kernel='rbf', C=1E6)
clf.fit(X, y)
Out[14]: SVC(C=1000000.0)
让我们使用之前定义的函数来可视化拟合并标识支持向量（参见图 43-8）。
In [15]: plt.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='autumn')
plot_svc_decision_function(clf)
plt.scatter(clf.support_vectors_[:, 0], clf.support_vectors_[:, 1],
s=300, lw=1, facecolors='none');
图 43-8. 核 SVM 对数据的拟合
使用这种核化支持向量机，我们学习到了一个适合的非线性决策边界。这种核变换策略在机器学习中经常被使用，将快速的线性方法转换为快速的非线性方法，特别适用于可以使用核技巧的模型。
调整 SVM：软化间隔
到目前为止，我们的讨论集中在非常干净的数据集上，其中存在完美的决策边界。但是如果您的数据有一定的重叠呢？例如，您可能有这样的数据（见图 43-9）。
In [16]: X, y = make_blobs(n_samples=100, centers=2,
random_state=0, cluster_std=1.2)
plt.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='autumn');
图 43-9. 具有一定重叠级别的数据
为了处理这种情况，SVM 实现中有一个“软化”间隔的修正因子：即，如果允许更好的拟合，某些点可以进入间隔。间隔的硬度由调整参数控制，通常称为C。对于很大的C，间隔是硬的，点不能位于其中。对于较小的C，间隔较软，并且可以包含一些点。
图 43-10 中显示的图表展示了通过软化间隔来改变C如何影响最终拟合的视觉效果：
In [17]: X, y = make_blobs(n_samples=100, centers=2,
random_state=0, cluster_std=0.8)
fig, ax = plt.subplots(1, 2, figsize=(16, 6))
fig.subplots_adjust(left=0.0625, right=0.95, wspace=0.1)
for axi, C in zip(ax, [10.0, 0.1]):
model = SVC(kernel='linear', C=C).fit(X, y)
axi.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='autumn')
plot_svc_decision_function(model, axi)
axi.scatter(model.support_vectors_[:, 0],
model.support_vectors_[:, 1],
s=300, lw=1, facecolors='none');
axi.set_title('C = {0:.1f}'.format(C), size=14)
图 43-10. C 参数对支持向量拟合的影响
C 的最佳值将取决于您的数据集，您应该使用交叉验证或类似的程序来调整此参数（参考第三十九章）。
示例：人脸识别
作为支持向量机在实际中的应用示例，让我们来看一下人脸识别问题。我们将使用“野外标记人脸”数据集，该数据集包含数千张各种公众人物的合并照片。Scikit-Learn 内置了该数据集的获取器：
In [18]: from sklearn.datasets import fetch_lfw_people
faces = fetch_lfw_people(min_faces_per_person=60)
print(faces.target_names)
print(faces.images.shape)
Out[18]: ['Ariel Sharon' 'Colin Powell' 'Donald Rumsfeld' 'George W Bush'
'Gerhard Schroeder' 'Hugo Chavez' 'Junichiro Koizumi' 'Tony Blair']
(1348, 62, 47)
让我们绘制几张这些人脸，看看我们正在处理的内容（见图 43-11）。
In [19]: fig, ax = plt.subplots(3, 5, figsize=(8, 6))
for i, axi in enumerate(ax.flat):
axi.imshow(faces.images[i], cmap='bone')
axi.set(xticks=[], yticks=[],
xlabel=faces.target_names[faces.target[i]])
图 43-11. 来自野外标记人脸数据集的示例
每个图像包含 62 × 47，约 3,000 个像素。我们可以简单地使用每个像素值作为特征，但通常使用某种预处理器来提取更有意义的特征更为有效；在这里，我们将使用主成分分析（见第四十五章）提取 150 个基本组分，以供支持向量机分类器使用。我们可以通过将预处理器和分类器打包到单个管道中来实现这一点：
In [20]: from sklearn.svm import SVC
from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
pca = PCA(n_components=150, whiten=True,
svd_solver='randomized', random_state=42)
svc = SVC(kernel='rbf', class_weight='balanced')
model = make_pipeline(pca, svc)
为了测试我们分类器的输出，我们将数据分割为训练集和测试集：
In [21]: from sklearn.model_selection import train_test_split
Xtrain, Xtest, ytrain, ytest = train_test_split(faces.data, faces.target,
random_state=42)
最后，我们可以使用网格搜索交叉验证来探索参数的组合。在这里，我们将调整C（控制边界硬度）和gamma（控制径向基函数核的大小），并确定最佳模型：
In [22]: from sklearn.model_selection import GridSearchCV
param_grid = {'svc__C': [1, 5, 10, 50],
'svc__gamma': [0.0001, 0.0005, 0.001, 0.005]}
grid = GridSearchCV(model, param_grid)
%time grid.fit(Xtrain, ytrain)
print(grid.best_params_)
Out[22]: CPU times: user 1min 19s, sys: 8.56 s, total: 1min 27s
Wall time: 36.2 s
{'svc__C': 10, 'svc__gamma': 0.001}
最优值集中在我们网格的中间；如果它们在边缘，我们将扩展网格以确保找到真正的最优值。
现在有了这个经过交叉验证的模型，我们可以预测测试数据的标签，这些数据模型尚未见过：
In [23]: model = grid.best_estimator_
yfit = model.predict(Xtest)
让我们看一些测试图像及其预测值（见图 43-12）。
In [24]: fig, ax = plt.subplots(4, 6)
for i, axi in enumerate(ax.flat):
axi.imshow(Xtest[i].reshape(62, 47), cmap='bone')
axi.set(xticks=[], yticks=[])
axi.set_ylabel(faces.target_names[yfit[i]].split()[-1],
color='black' if yfit[i] == ytest[i] else 'red')
fig.suptitle('Predicted Names; Incorrect Labels in Red', size=14);
图 43-12. 我们模型预测的标签
在这个小样本中，我们的最优估计器只误标了一个面孔（底部行的布什面孔被误标为布莱尔）。我们可以通过分类报告更好地了解我们估计器的性能，报告会逐标签列出恢复统计信息：
In [25]: from sklearn.metrics import classification_report
print(classification_report(ytest, yfit,
target_names=faces.target_names))
Out[25]: precision recall f1-score support
Ariel Sharon 0.65 0.73 0.69 15
Colin Powell 0.80 0.87 0.83 68
Donald Rumsfeld 0.74 0.84 0.79 31
George W Bush 0.92 0.83 0.88 126
Gerhard Schroeder 0.86 0.83 0.84 23
Hugo Chavez 0.93 0.70 0.80 20
Junichiro Koizumi 0.92 1.00 0.96 12
Tony Blair 0.85 0.95 0.90 42
accuracy 0.85 337
macro avg 0.83 0.84 0.84 337
weighted avg 0.86 0.85 0.85 337
我们还可以显示这些类别之间的混淆矩阵（见图 43-13）。
In [26]: from sklearn.metrics import confusion_matrix
import seaborn as sns
mat = confusion_matrix(ytest, yfit)
sns.heatmap(mat.T, square=True, annot=True, fmt='d',
cbar=False, cmap='Blues',
xticklabels=faces.target_names,
yticklabels=faces.target_names)
plt.xlabel('true label')
plt.ylabel('predicted label');
图 43-13. 面部数据的混淆矩阵
这帮助我们了解哪些标签可能会被估计器混淆。
对于一个现实世界的人脸识别任务，在这种任务中照片并未预先裁剪成漂亮的网格，面部分类方案唯一的区别在于特征选择：您需要使用更复杂的算法来找到面部，并提取与像素化无关的特征。对于这种应用，一个好的选择是利用OpenCV，它包括对一般图像和特别是人脸的先前实现的最新特征提取工具。
摘要
这是支持向量机背后原理的简明直观介绍。这些模型由于以下几个原因而是一种强大的分类方法：
-
它们依赖于相对较少的支持向量，因此紧凑且占用极少的内存空间。
-
一旦模型训练完成，预测阶段非常快速。
-
因为它们只受到靠近边界的点的影响，所以它们在处理高维数据时表现良好——即使是比样本更多维度的数据，这对其他算法来说是个挑战。
-
它们与核方法的集成使它们非常灵活，能够适应许多类型的数据。
然而，支持向量机（SVMs）也有几个缺点：
-
样本数量N的缩放为最坏情况下是𝒪 [ N 3 ]，或者对于高效实现是𝒪 [ N 2 ]。对于大量的训练样本，这种计算成本可能是限制性的。
-
结果强烈依赖于合适的软化参数
C的选择。必须通过交叉验证仔细选择，随着数据集增大，这可能是昂贵的。 -
结果没有直接的概率解释。可以通过内部交叉验证来估计（参见
SVC的probability参数），但这额外的估计是昂贵的。
考虑到这些特性，我通常只有在其他更简单、更快速、不需要过多调整的方法被证明不足以满足我的需求时，才会转向支持向量机（SVM）。尽管如此，如果你有足够的 CPU 周期来进行数据训练和交叉验证 SVM，这种方法可以带来出色的结果。
第四十四章：深入探讨：决策树和随机森林
之前，我们深入研究了一个简单的生成分类器（朴素贝叶斯；见 第四十一章）和一个强大的判别分类器（支持向量机；见 第四十三章）。在这里，我们将看看另一种强大的算法：一个称为 随机森林 的非参数算法。随机森林是一种 集成 方法的示例，意味着它依赖于聚合一组更简单的估算器的结果。这样的集成方法的一个令人惊讶的结果是，总和可以大于各部分之和：也就是说，多个估算器之间的多数投票的预测准确度最终可能会比任何进行投票的单个估算器的准确度更高！我们将在以下部分看到这方面的例子。
我们从标准导入开始：
In [1]: %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
推动随机森林的动机：决策树
随机森林是建立在决策树上的集成学习器的一个例子。因此，我们将首先讨论决策树本身。
决策树是极其直观的分类或标记对象的方式：你只需提出一系列旨在对分类进行精准定位的问题。例如，如果你想构建一个用于对徒步时遇到的动物进行分类的决策树，你可以构建如 图 44-1 所示的决策树。
图 44-1. 二叉决策树的示例^(1)
二元分割使其极其高效：在构造良好的树时，每个问题将使选项数量减少约一半，非常快速地将选项缩小，即使在大量类别中也是如此。当然，关键在于决定每一步要问什么问题。在决策树的机器学习实现中，问题通常采用数据中的轴对齐分割形式：即，树中的每个节点都使用一个特征内的截止值将数据分为两组。现在让我们看一个示例。
创建决策树
考虑以下二维数据，它具有四个类标签之一（参见 图 44-2）。
In [2]: from sklearn.datasets import make_blobs
X, y = make_blobs(n_samples=300, centers=4,
random_state=0, cluster_std=1.0)
plt.scatter(X[:, 0], X[:, 1], c=y, s=50, cmap='rainbow');
图 44-2. 决策树分类器的数据
基于这些数据构建的简单决策树将根据某些定量标准迭代地沿着一个或另一个轴将数据分割，并在每个级别根据其中的点的多数投票确定新区域的标签。图 44-3 展示了此数据的决策树分类器的前四个级别的可视化。
图 44-3. 决策树如何分割数据的可视化^(2)
注意，第一个分割后，上层每个点保持不变，因此无需进一步细分此分支。除了包含同一颜色的节点外，在每个级别每个区域再次沿着两个特征之一进行分割。
在 Scikit-Learn 中，可以使用DecisionTreeClassifier估计器来拟合决策树到我们的数据：
In [3]: from sklearn.tree import DecisionTreeClassifier
tree = DecisionTreeClassifier().fit(X, y)
让我们编写一个实用函数来帮助我们可视化分类器的输出：
In [4]: def visualize_classifier(model, X, y, ax=None, cmap='rainbow'):
ax = ax or plt.gca()
# Plot the training points
ax.scatter(X[:, 0], X[:, 1], c=y, s=30, cmap=cmap,
clim=(y.min(), y.max()), zorder=3)
ax.axis('tight')
ax.axis('off')
xlim = ax.get_xlim()
ylim = ax.get_ylim()
# fit the estimator
model.fit(X, y)
xx, yy = np.meshgrid(np.linspace(*xlim, num=200),
np.linspace(*ylim, num=200))
Z = model.predict(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)
# Create a color plot with the results
n_classes = len(np.unique(y))
contours = ax.contourf(xx, yy, Z, alpha=0.3,
levels=np.arange(n_classes + 1) - 0.5,
cmap=cmap, zorder=1)
ax.set(xlim=xlim, ylim=ylim)
现在我们可以看一下决策树分类的样子（参见图 44-4）。
In [5]: visualize_classifier(DecisionTreeClassifier(), X, y)
图 44-4. 决策树分类的可视化
如果您正在实时运行此笔记本，您可以使用在线附录中包含的辅助脚本来打开决策树构建过程的交互式可视化：
In [6]: # helpers_05_08 is found in the online appendix
import helpers_05_08
helpers_05_08.plot_tree_interactive(X, y);
Out[6]: interactive(children=(Dropdown(description='depth', index=1, options=(1, 5),
> value=5), Output()), _dom_classes...
注意，随着深度的增加，我们往往会得到非常奇怪形状的分类区域；例如，在深度为五时，在黄色和蓝色区域之间有一个高而瘦的紫色区域。显然，这不是真实的内在数据分布的结果，而更多地是数据的特定采样或噪声特性的结果。也就是说，即使在仅深度为五的情况下，这棵决策树明显地过拟合了我们的数据。
决策树和过拟合
这种过拟合事实上是决策树的一个普遍特性：很容易使树的深度过深，从而适应特定数据的细节，而不是它们抽取自的分布的总体特性。另一种看待这种过拟合的方法是查看在不同数据子集上训练的模型——例如，在图 44-5 中，我们训练了两棵不同的树，每棵树都使用了原始数据的一半。
图 44-5. 两棵随机决策树的示例^(3)
显然，在某些地方，两棵树产生一致的结果（例如，在四个角落），而在其他地方，两棵树给出非常不同的分类结果（例如，在任意两个簇之间的区域）。关键观察是，这种不一致往往发生在分类不确定的地方，因此通过使用这两棵树的信息，我们可能会得到更好的结果！
如果您正在实时运行此笔记本，以下函数将允许您交互地显示在数据的随机子集上训练的树的拟合情况：
In [7]: # helpers_05_08 is found in the online appendix
import helpers_05_08
helpers_05_08.randomized_tree_interactive(X, y)
Out[7]: interactive(children=(Dropdown(description='random_state', options=(0, 100),
> value=0), Output()), _dom_classes...
正如利用两棵树的信息可以改进我们的结果一样，我们可能期望利用许多树的信息进一步改进我们的结果。
集成估计器：随机森林
这种多个过拟合估计器组合以减少过拟合效应的概念，是支持一种称为bagging的集成方法的基础。Bagging 利用一个并行估计器的集合（可能是一个抓袋），每个估计器都会对数据过拟合，并对结果进行平均以找到更好的分类。随机化决策树的集成称为随机森林。
这种袋装分类可以通过 Scikit-Learn 的BaggingClassifier元估计器手动完成，如下所示（见图 44-6）。
In [8]: from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import BaggingClassifier
tree = DecisionTreeClassifier()
bag = BaggingClassifier(tree, n_estimators=100, max_samples=0.8,
random_state=1)
bag.fit(X, y)
visualize_classifier(bag, X, y)
在本例中，我们通过对训练点的随机 80%子集拟合每个估计器来随机化数据。在实践中，通过在如何选择分割时注入一些随机性来更有效地随机化决策树：这样每次都会使所有数据对拟合有贡献，但拟合结果仍具有所需的随机性。例如，在确定要分割哪个特征时，随机树可能从顶部几个特征中选择。您可以在Scikit-Learn 文档和其中的参考文献中阅读有关这些随机化策略的更多技术细节。
图 44-6. 随机决策树集成的决策边界
在 Scikit-Learn 中，这样一个优化的随机决策树集成是通过RandomForestClassifier估计器实现的，它自动处理所有随机化。你只需选择一些估计器，它将非常快速地（如果需要的话是并行的）拟合树的集成（见图 44-7）。
In [9]: from sklearn.ensemble import RandomForestClassifier
model = RandomForestClassifier(n_estimators=100, random_state=0)
visualize_classifier(model, X, y);
图 44-7. 随机森林的决策边界，这是一组优化的决策树集成
我们看到通过对一百个随机扰动模型进行平均，最终得到一个与我们关于参数空间如何分割的直觉更接近的整体模型。
随机森林回归
在前一节中，我们考虑了随机森林在分类的上下文中。随机森林也可以在回归的情况下工作（即使用连续变量而不是分类变量）。用于此目的的估计器是RandomForestRegressor，其语法与我们之前看到的非常相似。
考虑以下数据，这些数据来自快速和慢速振荡的组合（见图 44-8）。
In [10]: rng = np.random.RandomState(42)
x = 10 * rng.rand(200)
def model(x, sigma=0.3):
fast_oscillation = np.sin(5 * x)
slow_oscillation = np.sin(0.5 * x)
noise = sigma * rng.randn(len(x))
return slow_oscillation + fast_oscillation + noise
y = model(x)
plt.errorbar(x, y, 0.3, fmt='o');
图 44-8. 随机森林回归的数据
使用随机森林回归器，我们可以找到最佳拟合曲线（见图 44-9）。
In [11]: from sklearn.ensemble import RandomForestRegressor
forest = RandomForestRegressor(200)
forest.fit(x[:, None], y)
xfit = np.linspace(0, 10, 1000)
yfit = forest.predict(xfit[:, None])
ytrue = model(xfit, sigma=0)
plt.errorbar(x, y, 0.3, fmt='o', alpha=0.5)
plt.plot(xfit, yfit, '-r');
plt.plot(xfit, ytrue, '-k', alpha=0.5);
图 44-9. 随机森林模型拟合数据
这里显示了真实模型的平滑灰色曲线，而随机森林模型则通过锯齿状红色曲线展示。非参数随机森林模型足够灵活，能够拟合多期数据，而无需指定多期模型！
示例：用于分类数字的随机森林
在第三十八章中，我们通过一个使用 Scikit-Learn 提供的数字数据集的示例来工作。让我们再次使用它来看看随机森林分类器在这种情况下的应用：
In [12]: from sklearn.datasets import load_digits
digits = load_digits()
digits.keys()
Out[12]: dict_keys(['data', 'target', 'frame', 'feature_names', 'target_names',
> 'images', 'DESCR'])
为了提醒我们正在查看的内容，我们将可视化前几个数据点（参见图 44-10）。
In [13]: # set up the figure
fig = plt.figure(figsize=(6, 6)) # figure size in inches
fig.subplots_adjust(left=0, right=1, bottom=0, top=1,
hspace=0.05, wspace=0.05)
# plot the digits: each image is 8x8 pixels
for i in range(64):
ax = fig.add_subplot(8, 8, i + 1, xticks=[], yticks=[])
ax.imshow(digits.images[i], cmap=plt.cm.binary, interpolation='nearest')
# label the image with the target value
ax.text(0, 7, str(digits.target[i]))
图 44-10. 数字数据的表示
我们可以使用随机森林对数字进行分类，如下所示：
In [14]: from sklearn.model_selection import train_test_split
Xtrain, Xtest, ytrain, ytest = train_test_split(digits.data, digits.target,
random_state=0)
model = RandomForestClassifier(n_estimators=1000)
model.fit(Xtrain, ytrain)
ypred = model.predict(Xtest)
让我们看看这个分类器的分类报告：
In [15]: from sklearn import metrics
print(metrics.classification_report(ypred, ytest))
Out[15]: precision recall f1-score support
0 1.00 0.97 0.99 38
1 0.98 0.98 0.98 43
2 0.95 1.00 0.98 42
3 0.98 0.96 0.97 46
4 0.97 1.00 0.99 37
5 0.98 0.96 0.97 49
6 1.00 1.00 1.00 52
7 1.00 0.96 0.98 50
8 0.94 0.98 0.96 46
9 0.98 0.98 0.98 47
accuracy 0.98 450
macro avg 0.98 0.98 0.98 450
weighted avg 0.98 0.98 0.98 450
并且为了更直观，绘制混淆矩阵（参见图 44-11）。
In [16]: from sklearn.metrics import confusion_matrix
import seaborn as sns
mat = confusion_matrix(ytest, ypred)
sns.heatmap(mat.T, square=True, annot=True, fmt='d',
cbar=False, cmap='Blues')
plt.xlabel('true label')
plt.ylabel('predicted label');
我们发现，一个简单的未调整的随机森林能够对数字数据进行相当准确的分类。
图 44-11. 使用随机森林进行数字分类的混淆矩阵
概要
本章简要介绍了集成估计器的概念，特别是随机森林，它是随机化决策树的集成。随机森林是一种功能强大的方法，具有多个优点：
-
由于底层决策树的简单性，训练和预测都非常快。此外，由于每棵树都是独立实体，因此这两个任务可以直接并行化。
-
多棵树允许进行概率分类：估算器的多数投票给出了概率的估计（在 Scikit-Learn 中通过
predict_proba方法访问）。 -
非参数模型非常灵活，因此在其他估算器欠拟合的任务上表现良好。
随机森林的一个主要缺点是结果不易解释：也就是说，如果你想对分类模型的含义得出结论，随机森林可能不是最佳选择。
^(1) 生成此图的代码可在在线附录中找到。
^(2) 生成此图的代码可在在线附录中找到。
^(3) 生成此图的代码可在在线附录中找到。
第四十五章：深入解析：主成分分析
到目前为止，我们一直在深入研究监督学习估计器：那些根据标记的训练数据预测标签的估计器。在这里，我们开始研究几个无监督估计器，这些估计器可以突出数据的一些有趣方面，而不需要参考任何已知的标签。
在本章中，我们将探讨或许是最广泛使用的无监督算法之一，即主成分分析（PCA）。PCA 本质上是一种降维算法，但它也可以作为可视化工具、噪声过滤器、特征提取和工程等工具使用。在简短讨论了 PCA 算法的概念后，我们将探索其进一步应用的几个例子。
我们从标准导入开始：
In [1]: %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
引入主成分分析
主成分分析是一种快速灵活的无监督数据降维方法，我们在第 38 章中简要介绍过。通过查看一个二维数据集，最容易理解其行为。考虑这 200 个点（见图 45-1）。
In [2]: rng = np.random.RandomState(1)
X = np.dot(rng.rand(2, 2), rng.randn(2, 200)).T
plt.scatter(X[:, 0], X[:, 1])
plt.axis('equal');
图 45-1. 主成分分析演示数据
从眼睛的角度来看，很明显x和y变量之间存在近乎线性的关系。这让人想起我们在第 42 章中探索的线性回归数据，但这里的问题设置略有不同：与尝试从x值预测y值不同，无监督学习问题试图学习x和y值之间的关系。
在主成分分析中，通过找到数据中的一组主轴来量化这种关系，并使用这些轴来描述数据集。使用 Scikit-Learn 的PCA估计器，我们可以按如下方式计算：
In [3]: from sklearn.decomposition import PCA
pca = PCA(n_components=2)
pca.fit(X)
Out[3]: PCA(n_components=2)
拟合从数据中学习到一些量，最重要的是分量和解释方差：
In [4]: print(pca.components_)
Out[4]: [[-0.94446029 -0.32862557]
[-0.32862557 0.94446029]]
In [5]: print(pca.explained_variance_)
Out[5]: [0.7625315 0.0184779]
要理解这些数字的含义，让我们将它们视为向量在输入数据上的可视化，使用分量定义向量的方向和解释方差定义向量的平方长度（见图 45-2）。
In [6]: def draw_vector(v0, v1, ax=None):
ax = ax or plt.gca()
arrowprops=dict(arrowstyle='->', linewidth=2,
shrinkA=0, shrinkB=0)
ax.annotate('', v1, v0, arrowprops=arrowprops)
# plot data
plt.scatter(X[:, 0], X[:, 1], alpha=0.2)
for length, vector in zip(pca.explained_variance_, pca.components_):
v = vector * 3 * np.sqrt(length)
draw_vector(pca.mean_, pca.mean_ + v)
plt.axis('equal');
图 45-2. 数据中主轴的可视化
这些向量代表了数据的主轴，每个向量的长度表示描述数据分布中该轴的“重要性”的指标——更精确地说，它是数据在投影到该轴上时的方差的度量。将每个数据点投影到主轴上得到了数据的主成分。
如果我们将这些主成分与原始数据一起绘制，我们会看到在图 45-3 中显示的图形。
图 45-3. 数据中转换后的主轴^(1)
将数据轴转换为主轴的这种变换是仿射变换，这意味着它由平移、旋转和均匀缩放组成。
虽然这种寻找主成分的算法可能看起来只是一种数学上的好奇，但它事实证明在机器学习和数据探索领域有着非常广泛的应用。
PCA 作为降维
使用 PCA 进行降维涉及将一个或多个最小主成分置零，结果是数据的低维投影，保留了最大的数据方差。
这里是使用 PCA 作为降维转换的示例：
In [7]: pca = PCA(n_components=1)
pca.fit(X)
X_pca = pca.transform(X)
print("original shape: ", X.shape)
print("transformed shape:", X_pca.shape)
Out[7]: original shape: (200, 2)
transformed shape: (200, 1)
转换后的数据已经降维到了单一维度。为了理解这种降维的效果，我们可以对这个降维后的数据执行逆变换，并将其与原始数据一起绘制（参见图 45-4）。
In [8]: X_new = pca.inverse_transform(X_pca)
plt.scatter(X[:, 0], X[:, 1], alpha=0.2)
plt.scatter(X_new[:, 0], X_new[:, 1], alpha=0.8)
plt.axis('equal');
图 45-4. PCA 作为降维的可视化
浅色点是原始数据，而深色点是投影版本。这清楚地展示了 PCA 降维的含义：移除沿着最不重要的主轴或轴的信息，只留下数据方差最大的分量。被削减的方差比例（与在前述图中形成的线周围的点的散布成比例）大致是这种降维中丢失的“信息”量的度量。
这个降维后的数据集在某些意义上是“足够好”的，能够编码点之间最重要的关系：尽管将数据特征的数量减少了 50%，但数据点之间的整体关系大部分得到了保留。
PCA 用于可视化：手写数字
在只有两个维度时，降维的实用性可能并不完全明显，但当我们查看高维数据时，它变得更加清晰。为了看到这一点，让我们快速看一下我们在第四十四章中处理的数字数据集应用 PCA 的情况。
我们将从加载数据开始：
In [9]: from sklearn.datasets import load_digits
digits = load_digits()
digits.data.shape
Out[9]: (1797, 64)
请回忆，数字数据集由 8×8 像素的图像组成，这意味着它们是 64 维的。为了对这些点之间的关系有一些直观的理解，我们可以使用 PCA 将它们投影到一个更易处理的维度，比如说两个维度：
In [10]: pca = PCA(2) # project from 64 to 2 dimensions
projected = pca.fit_transform(digits.data)
print(digits.data.shape)
print(projected.shape)
Out[10]: (1797, 64)
(1797, 2)
现在我们可以绘制每个点的前两个主成分，以了解数据，如在图 45-5 中所见。
In [11]: plt.scatter(projected[:, 0], projected[:, 1],
c=digits.target, edgecolor='none', alpha=0.5,
cmap=plt.cm.get_cmap('rainbow', 10))
plt.xlabel('component 1')
plt.ylabel('component 2')
plt.colorbar();
图 45-5. 应用于手写数字数据的 PCA
回想一下这些组件的含义：完整的数据是一个 64 维点云，而这些点是每个数据点沿着最大方差方向的投影。基本上，我们在 64 维空间中找到了允许我们在两个维度上看到数据布局的最佳拉伸和旋转，而且这是无监督的方式完成的，即没有参考标签。
组件的含义是什么？
在这里，我们可以再进一步，开始问减少的维度意味着什么。这种意义可以用基向量的组合来理解。例如，训练集中的每个图像由一组 64 个像素值定义，我们将其称为向量x：
x = [ x 1 , x 2 , x 3 ⋯ x 64 ]
我们可以用像素基础来考虑这一点。也就是说，为了构建图像，我们将向量的每个元素乘以它描述的像素，然后将结果相加以构建图像：
image ( x ) = x 1 · ( pixel 1 ) + x 2 · ( pixel 2 ) + x 3 · ( pixel 3 ) ⋯ x 64 · ( pixel 64 )
我们可以想象一种降低数据维度的方法是将除了少数基向量之外的所有值归零。例如，如果我们仅使用前八个像素，我们得到数据的八维投影（见图 45-6）。然而，这并不太反映整个图像：我们几乎丢弃了 90%的像素！
图 45-6. 通过丢弃像素实现的天真降维^(2)
上排面板显示单独的像素，下排显示这些像素对图像构建的累积贡献。仅使用八个像素基础组件，我们只能构建 64 像素图像的一小部分。如果我们继续这个序列并使用所有 64 个像素，我们将恢复原始图像。
但像素级表示并不是唯一的基础选择。我们还可以使用其他基函数，每个基函数都包含来自每个像素的一些预定义贡献，并编写如下内容：
i m a g e ( x ) = mean + x 1 · ( basis 1 ) + x 2 · ( basis 2 ) + x 3 · ( basis 3 ) ⋯
PCA 可以被看作是选择最优基函数的过程，使得仅添加前几个基函数就足以适当地重构数据集中的大部分元素。主成分作为我们数据的低维表示，实际上只是乘以这一系列中每个元素的系数。图 45-7 展示了使用平均值加上前八个 PCA 基函数重建相同数字的类似描述。
图 45-7. 通过丢弃最不重要的主成分实现的更复杂的降维（与图 45-6 比较)^(3)
与像素基础不同，PCA 基础允许我们仅通过平均值加上八个组件来恢复输入图像的显著特征！每个像素在每个组件中的量是我们二维示例中向量方向的必然结果。这就是 PCA 提供数据低维表示的方式：它发现一组比输入数据的原生像素基础更有效的基础函数。
选择组件的数量
在实际应用中使用 PCA 的一个重要部分是估计需要多少个组件来描述数据。这可以通过查看组件数量作为累积解释方差比的函数来确定（参见图 45-8）。
In [12]: pca = PCA().fit(digits.data)
plt.plot(np.cumsum(pca.explained_variance_ratio_))
plt.xlabel('number of components')
plt.ylabel('cumulative explained variance');
此曲线量化了在前N个组件中包含的总 64 维方差的比例。例如，我们看到对于数字数据，前 10 个组件包含大约 75%的方差，而您需要约 50 个组件来描述接近 100%的方差。
图 45-8. 累积解释方差，用于衡量 PCA 保留数据内容的效果
这告诉我们，我们的二维投影丢失了大量信息（由解释方差度量），我们需要大约 20 个组件来保留 90%的方差。查看高维数据集的此图可以帮助您了解其特征中存在的冗余水平。
PCA 作为噪声过滤器
PCA 还可以用作噪声数据的过滤方法。其思想是：任何方差远大于噪声影响的成分应该相对不受噪声影响。因此，如果您仅使用主成分的最大子集重建数据，则应优先保留信号并丢弃噪声。
让我们看看数字数据的情况。首先，我们将绘制几个无噪声输入样本（图 45-9）。
In [13]: def plot_digits(data):
fig, axes = plt.subplots(4, 10, figsize=(10, 4),
subplot_kw={'xticks':[], 'yticks':[]},
gridspec_kw=dict(hspace=0.1, wspace=0.1))
for i, ax in enumerate(axes.flat):
ax.imshow(data[i].reshape(8, 8),
cmap='binary', interpolation='nearest',
clim=(0, 16))
plot_digits(digits.data)
图 45-9. 无噪声的数字
现在让我们添加一些随机噪声以创建一个带噪声的数据集，并重新绘制它（图 45-10）。
In [14]: rng = np.random.default_rng(42)
rng.normal(10, 2)
Out[14]: 10.609434159508863
In [15]: rng = np.random.default_rng(42)
noisy = rng.normal(digits.data, 4)
plot_digits(noisy)
图 45-10. 添加了高斯随机噪声的数字
可视化使得随机噪声的存在变得明显。让我们在嘈杂数据上训练一个 PCA 模型，并要求投影保留 50%的方差：
In [16]: pca = PCA(0.50).fit(noisy)
pca.n_components_
Out[16]: 12
这里 50%的方差相当于 12 个主成分，而原始的 64 个特征。现在我们计算这些成分，然后使用变换的逆来重构经过滤波的数字；图 45-11 展示了结果。
In [17]: components = pca.transform(noisy)
filtered = pca.inverse_transform(components)
plot_digits(filtered)
图 45-11. 使用 PCA 进行“去噪”处理的数字
这种信号保留/噪声过滤特性使得 PCA 成为非常有用的特征选择例程——例如，不是在非常高维度的数据上训练分类器，而是在较低维度的主成分表示上训练分类器，这将自动过滤输入中的随机噪声。
示例：特征脸
我们之前探讨了使用 PCA 技术作为支持向量机的特征选择器进行人脸识别的示例（见第四十三章）。现在让我们回顾一下，并深入探讨这背后的更多内容。回想一下，我们使用的是由 Scikit-Learn 提供的 Labeled Faces in the Wild（LFW）数据集：
In [18]: from sklearn.datasets import fetch_lfw_people
faces = fetch_lfw_people(min_faces_per_person=60)
print(faces.target_names)
print(faces.images.shape)
Out[18]: ['Ariel Sharon' 'Colin Powell' 'Donald Rumsfeld' 'George W Bush'
'Gerhard Schroeder' 'Hugo Chavez' 'Junichiro Koizumi' 'Tony Blair']
(1348, 62, 47)
让我们看看涵盖此数据集的主轴。由于这是一个大数据集，我们将在PCA估计器中使用"random"特征求解器：它使用随机方法来更快地近似前N个主成分，而不是标准方法，以牺牲一些准确性。这种权衡在高维数据（这里接近 3,000 维）中非常有用。我们将看一看前 150 个成分：
In [19]: pca = PCA(150, svd_solver='randomized', random_state=42)
pca.fit(faces.data)
Out[19]: PCA(n_components=150, random_state=42, svd_solver='randomized')
在这种情况下，可以通过可视化与前几个主成分相关联的图像来进行探索（这些成分在技术上称为特征向量，因此这些类型的图像通常被称为特征脸；正如您在图 45-12 中所见，它们听起来就像它们看起来那样可怕）：
In [20]: fig, axes = plt.subplots(3, 8, figsize=(9, 4),
subplot_kw={'xticks':[], 'yticks':[]},
gridspec_kw=dict(hspace=0.1, wspace=0.1))
for i, ax in enumerate(axes.flat):
ax.imshow(pca.components_[i].reshape(62, 47), cmap='bone')
图 45-12. 从 LFW 数据集学习的特征脸的可视化
结果非常有趣，并且为我们提供了关于图像变化的见解：例如，前几个特征脸（左上角）似乎与脸部的光照角度有关，而后来的主成分似乎在挑选出特定的特征，如眼睛、鼻子和嘴唇。让我们看一看这些成分的累计方差，以查看投影保留了多少数据信息（见图 45-13）。
In [21]: plt.plot(np.cumsum(pca.explained_variance_ratio_))
plt.xlabel('number of components')
plt.ylabel('cumulative explained variance');
图 45-13. LFW 数据的累计解释方差
我们选择的 150 个组件占了超过 90%的方差。这使我们相信，使用这 150 个组件，我们将恢复数据的大部分基本特征。为了更具体化，我们可以将输入图像与从这些 150 个组件重建的图像进行比较（参见图 45-14）。
In [22]: # Compute the components and projected faces
pca = pca.fit(faces.data)
components = pca.transform(faces.data)
projected = pca.inverse_transform(components)
In [23]: # Plot the results
fig, ax = plt.subplots(2, 10, figsize=(10, 2.5),
subplot_kw={'xticks':[], 'yticks':[]},
gridspec_kw=dict(hspace=0.1, wspace=0.1))
for i in range(10):
ax[0, i].imshow(faces.data[i].reshape(62, 47), cmap='binary_r')
ax[1, i].imshow(projected[i].reshape(62, 47), cmap='binary_r')
ax[0, 0].set_ylabel('full-dim\ninput')
ax[1, 0].set_ylabel('150-dim\nreconstruction');
图 45-14. LFW 数据的 150 维 PCA 重建
这里的顶部行显示了输入图像，而底部行显示了仅从约 3000 个初始特征中的 150 个进行图像重建。这种可视化清楚地说明了 PCA 特征选择在第四十三章中为何如此成功：虽然它将数据的维度减少了近 20 倍，但投影图像包含足够的信息，使我们可以通过肉眼识别每个图像中的个体。这意味着我们的分类算法只需要在 150 维数据上进行训练，而不是 3000 维数据，根据我们选择的特定算法，这可能会导致更高效的分类。
摘要
在本章中，我们探讨了主成分分析在降维、高维数据可视化、噪声过滤和高维数据特征选择中的应用。由于其多功能性和可解释性，PCA 已被证明在各种背景和学科中都非常有效。对于任何高维数据集，我倾向于从 PCA 开始，以便可视化数据点之间的关系（就像我们在数字数据中所做的那样），理解数据中的主要方差（就像我们在特征脸中所做的那样），并理解内在的维度（通过绘制解释方差比）。当然，PCA 并非对每个高维数据集都有用，但它提供了一条直观和高效的路径，以洞察高维数据。
PCA 的主要弱点是它往往受到数据中异常值的影响。因此，已经开发了几种鲁棒性较强的 PCA 变体，其中许多变体通过迭代地丢弃初始组件描述不佳的数据点来作用。Scikit-Learn 在sklearn.decomposition子模块中包括了许多有趣的 PCA 变体；一个例子是SparsePCA，它引入了一个正则化项（参见第四十二章），用于强制组件的稀疏性。
在接下来的章节中，我们将研究其他建立在 PCA 思想基础上的无监督学习方法。
^(1) 可在在线附录中找到生成此图的代码。
^(2) 可在在线附录中找到生成此图的代码。
^(3) 生成此图的代码可以在在线附录中找到。
第四十六章：深入探讨：流形学习
在上一章中，我们看到 PCA 可以用于降维，减少数据集的特征数，同时保持点之间的基本关系。虽然 PCA 灵活、快速且易于解释，但当数据中存在非线性关系时，它的表现并不理想，我们很快将看到一些例子。
为了解决这一不足，我们可以转向流形学习算法——一类无监督估计器，旨在将数据集描述为嵌入在高维空间中的低维流形。当你想到流形时，我建议想象一张纸：这是一个二维物体，存在于我们熟悉的三维世界中。
在流形学的术语中，你可以将这张纸片看作是嵌入在三维空间中的二维流形。在三维空间中旋转、重新定向或拉伸这张纸片并不会改变其平面几何特性：这些操作类似于线性嵌入。如果你将纸张弯曲、卷曲或揉皱，它仍然是一个二维流形，但是嵌入到三维空间的方式不再是线性的。流形学算法旨在学习纸张的基本二维特性，即使它被扭曲以填充三维空间。
在这里，我们将研究多种流形方法，深入探讨其中的一些技术：多维尺度法（MDS）、局部线性嵌入（LLE）和等距映射（Isomap）。
我们从标准导入开始：
In [1]: %matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
import numpy as np
流形学习：“HELLO”
为了更清晰地阐明这些概念，让我们首先生成一些二维数据，以便用来定义一个流形。这里是一个能创建“HELLO”形状数据的函数：
In [2]: def make_hello(N=1000, rseed=42):
# Make a plot with "HELLO" text; save as PNG
fig, ax = plt.subplots(figsize=(4, 1))
fig.subplots_adjust(left=0, right=1, bottom=0, top=1)
ax.axis('off')
ax.text(0.5, 0.4, 'HELLO', va='center', ha='center',
weight='bold', size=85)
fig.savefig('hello.png')
plt.close(fig)
# Open this PNG and draw random points from it
from matplotlib.image import imread
data = imread('hello.png')[::-1, :, 0].T
rng = np.random.RandomState(rseed)
X = rng.rand(4 * N, 2)
i, j = (X * data.shape).astype(int).T
mask = (data[i, j] < 1)
X = X[mask]
X[:, 0] *= (data.shape[0] / data.shape[1])
X = X[:N]
return X[np.argsort(X[:, 0])]
让我们调用这个函数，并可视化生成的数据（参见图 46-1）。
In [3]: X = make_hello(1000)
colorize = dict(c=X[:, 0], cmap=plt.cm.get_cmap('rainbow', 5))
plt.scatter(X[:, 0], X[:, 1], **colorize)
plt.axis('equal');
输出是二维的，由“HELLO”形状的点组成。这种数据形式将帮助我们直观地了解这些算法的作用。
图 46-1. 用于流形学习的数据
多维尺度法
查看这样的数据，我们可以看到数据集中特定的x和y值并不是数据最基本的描述：我们可以缩放、收缩或旋转数据，而“HELLO”仍然是显而易见的。例如，如果我们使用旋转矩阵旋转数据，x和y值会改变，但数据基本上仍然是相同的（参见图 46-2）。
In [4]: def rotate(X, angle):
theta = np.deg2rad(angle)
R = [[np.cos(theta), np.sin(theta)],
[-np.sin(theta), np.cos(theta)]]
return np.dot(X, R)
X2 = rotate(X, 20) + 5
plt.scatter(X2[:, 0], X2[:, 1], **colorize)
plt.axis('equal');
图 46-2. 旋转的数据集
这证实了x和y值不一定是数据关系中的基本要素。在这种情况下，基本的是数据集中每个点之间的距离。表示这一点的常见方法是使用距离矩阵：对于N个点，我们构造一个N × N数组，使得条目( i , j )包含点i和点j之间的距离。让我们使用 Scikit-Learn 的高效pairwise_distances函数来为我们的原始数据做到这一点：
In [5]: from sklearn.metrics import pairwise_distances
D = pairwise_distances(X)
D.shape
Out[5]: (1000, 1000)
正如承诺的那样，对于我们的N=1,000 个点，我们获得了一个 1,000 × 1,000 的矩阵，可以像这样进行可视化（参见图 46-3）。
In [6]: plt.imshow(D, zorder=2, cmap='viridis', interpolation='nearest')
plt.colorbar();
图 46-3. 点之间的成对距离可视化
如果我们类似地为我们旋转和平移的数据构建一个距离矩阵，我们会看到它是相同的：
In [7]: D2 = pairwise_distances(X2)
np.allclose(D, D2)
Out[7]: True
这个距离矩阵给出了一个与旋转和平移无关的数据表示，但在图 46-3 中的矩阵可视化并不完全直观。在那里显示的表示中，我们失去了数据中有趣结构的任何可见迹象：“HELLO”。
此外，从(x, y)坐标计算距离矩阵很简单，但将距离转换回x和y坐标却相当困难。这正是多维缩放算法的目标：给定点之间的距离矩阵，它恢复数据的D维坐标表示。让我们看看它如何处理我们的距离矩阵，使用precomputed表示我们正在传递一个距离矩阵（见图 46-4）。
In [8]: from sklearn.manifold import MDS
model = MDS(n_components=2, dissimilarity='precomputed', random_state=1701)
out = model.fit_transform(D)
plt.scatter(out[:, 0], out[:, 1], **colorize)
plt.axis('equal');
图 46-4. 从成对距离计算得到的 MDS 嵌入
MDS 算法使用描述数据点之间关系的N × N距离矩阵之一，仅恢复了我们数据的可能的二维坐标表示。
MDS 作为流形学习
当我们考虑到距离矩阵可以从任何维度的数据中计算时，这种用处变得更加明显。例如，我们不仅可以简单地将数据在二维平面上旋转，还可以使用以下函数将其投影到三维（本质上是早期使用的旋转矩阵的三维推广）。
In [9]: def random_projection(X, dimension=3, rseed=42):
assert dimension >= X.shape[1]
rng = np.random.RandomState(rseed)
C = rng.randn(dimension, dimension)
e, V = np.linalg.eigh(np.dot(C, C.T))
return np.dot(X, V[:X.shape[1]])
X3 = random_projection(X, 3)
X3.shape
Out[9]: (1000, 3)
让我们可视化这些点，看看我们在处理什么（见图 46-5）。
In [10]: from mpl_toolkits import mplot3d
ax = plt.axes(projection='3d')
ax.scatter3D(X3[:, 0], X3[:, 1], X3[:, 2],
**colorize);
图 46-5. 数据线性嵌入到三维空间中
现在，我们可以要求MDS估计器输入这三维数据，计算距离矩阵，然后确定这个距离矩阵的最优二维嵌入。结果恢复了原始数据的表示，如图 46-6 所示。
In [11]: model = MDS(n_components=2, random_state=1701)
out3 = model.fit_transform(X3)
plt.scatter(out3[:, 0], out3[:, 1], **colorize)
plt.axis('equal');
这基本上是流形学习估计器的目标：在给定高维嵌入数据的情况下，寻找保持数据内部某些关系的低维表示。在 MDS 的情况下，保持的量是每对点之间的距离。
图 46-6. 三维数据的 MDS 嵌入恢复了输入，经过了旋转和反射。
非线性嵌入：MDS 失败的地方
到目前为止，我们讨论的是线性嵌入，这本质上是将数据旋转、平移和缩放到更高维空间中。当嵌入是非线性的时候，即超出这简单操作集合时，MDS 失效了。考虑下面的嵌入，它将输入扭曲成了三维空间中的“S”形状：
In [12]: def make_hello_s_curve(X):
t = (X[:, 0] - 2) * 0.75 * np.pi
x = np.sin(t)
y = X[:, 1]
z = np.sign(t) * (np.cos(t) - 1)
return np.vstack((x, y, z)).T
XS = make_hello_s_curve(X)
这同样是三维数据，但正如我们在图 46-7 中所见，嵌入更加复杂。
In [13]: from mpl_toolkits import mplot3d
ax = plt.axes(projection='3d')
ax.scatter3D(XS[:, 0], XS[:, 1], XS[:, 2],
**colorize);
图 46-7. 数据非线性嵌入到三维空间中
数据点之间的基本关系仍然存在，但这次数据以非线性的方式进行了转换：它被包裹成了一个“S”形状。
如果我们尝试在这些数据上使用简单的 MDS 算法，它无法“展开”这个非线性嵌入，我们会丢失嵌入流形中基本的关系（见图 46-8）。
In [14]: from sklearn.manifold import MDS
model = MDS(n_components=2, random_state=2)
outS = model.fit_transform(XS)
plt.scatter(outS[:, 0], outS[:, 1], **colorize)
plt.axis('equal');
图 46-8. 应用于非线性数据的 MDS 算法；它未能恢复底层结构
最佳的二维线性嵌入并没有展开 S 形曲线，而是丢弃了原始的 y 轴。
非线性流形：局部线性嵌入
在这里我们该怎么继续？退一步来看，我们可以看到问题的根源在于，MDS 试图在构建嵌入时保持远距离点之间的距离。但如果我们改变算法，使其仅保持附近点之间的距离会怎样？结果的嵌入会更接近我们想要的。
从视觉上看，我们可以将其想象成如图 46-9 所示。
这里每条淡淡的线表示一个应该在嵌入中保留的距离。左侧是 MDS 使用的模型的表示：它试图保持数据集中每对点之间的距离。右侧是一种名为局部线性嵌入的流形学习算法使用的模型的表示：它不是保留所有距离，而是试图仅保留相邻点之间的距离（在这种情况下，每个点的最近 100 个邻居）。
考虑左侧面板，我们可以看到为什么 MDS 失败了：没有办法展开这些数据同时充分保持两点之间绘制的每条线的长度。另一方面，对于右侧面板，情况看起来更加乐观。我们可以想象以一种方式展开数据，以保持线的长度大致相同。这正是 LLE 通过全局优化反映这种逻辑的成本函数所做的。
图 46-9. MDS 和 LLE 之间点之间联系的表示^(1)
LLE 有许多变体；在这里，我们将使用改进的 LLE算法来恢复嵌入的二维流形。总的来说，改进的 LLE 比算法的其他变体更能够恢复具有很少扭曲的明确定义流形（参见图 46-10）。
In [15]: from sklearn.manifold import LocallyLinearEmbedding
model = LocallyLinearEmbedding(
n_neighbors=100, n_components=2,
method='modified', eigen_solver='dense')
out = model.fit_transform(XS)
fig, ax = plt.subplots()
ax.scatter(out[:, 0], out[:, 1], **colorize)
ax.set_ylim(0.15, -0.15);
与我们原始流形相比，结果仍然有些扭曲，但捕捉到了数据中的基本关系！
图 46-10. 局部线性嵌入可以从非线性嵌入的输入中恢复底层数据
对流形方法的一些思考
尽管这些示例可能很引人注目，但在实践中，流形学习技术往往很难处理，因此很少被用于除了高维数据的简单定性可视化之外的任何其他用途。
以下是流形学习的一些特定挑战，这些挑战都与 PCA 相比非常不利：
-
流形学习中，没有处理缺失数据的良好框架。相比之下，在 PCA 中有直接的迭代方法来处理缺失数据。
-
在流形学习中，数据中的噪声存在可以“短路”流形并显著改变嵌入的情况。相比之下，PCA 自然地从最重要的组件中过滤噪声。
-
流形嵌入结果通常高度依赖于选择的邻居数，并且通常没有一种可靠的定量方法来选择最佳邻居数。相比之下，PCA 不涉及这样的选择。
-
在流形学习中，确定全局最佳输出维度的难度很大。相比之下，PCA 可以根据解释的方差来确定输出维度的数量。
-
在流形学习中，嵌入维度的含义并不总是清晰的。在 PCA 中，主成分有一个非常明确的含义。
-
在流形学习中，流形方法的计算开销按O [ N 2 ] 或 O [ N 3 ] 的方式扩展。对于 PCA，存在一些随机化方法通常要快得多（尽管请参考megaman package以获取更多可扩展的流形学习实现）。
总结一下，流形学习方法相对于 PCA 唯一明显的优势是它们能够保留数据中的非线性关系；因此，我倾向于仅在首先使用 PCA 探索数据后才使用流形方法探索数据。
Scikit-Learn 实现了除了 LLE 和 Isomap 之外的几种常见流形学习变体（我们在前几章中使用过 Isomap，将在下一节中查看）：Scikit-Learn 文档对它们进行了很好的讨论和比较。根据我的经验，我会提出以下建议：
-
对于像我们之前看到的 S 曲线这样的玩具问题，LLE 及其变体（特别是修改的 LLE）表现非常好。这在
sklearn.manifold.LocallyLinearEmbedding中实现。 -
对于来自现实世界来源的高维数据，LLE 通常会产生较差的结果，而 Isomap 似乎通常会导致更有意义的嵌入。这在
sklearn.manifold.Isomap中实现。 -
对于高度聚类的数据，t-分布随机近邻嵌入（t-SNE）似乎效果非常好，尽管与其他方法相比速度较慢。这在
sklearn.manifold.TSNE中实现。
如果你对这些方法的运行感兴趣，我建议你在这一节的数据上分别运行每种方法。
例如：Faces 上的 Isomap
流形学习经常用于理解高维数据点之间的关系。高维数据的一个常见案例是图像：例如，每个包含 1,000 个像素的图像集可以被视为一个 1,000 维空间中的点集，其中每个像素的亮度定义了该维度上的坐标。
为了说明，让我们在一些来自 Labeled Faces in the Wild 数据集的数据上应用 Isomap，我们之前在第四十三章和第四十五章中看到过。运行此命令将下载数据集并将其缓存到您的主目录供以后使用：
In [16]: from sklearn.datasets import fetch_lfw_people
faces = fetch_lfw_people(min_faces_per_person=30)
faces.data.shape
Out[16]: (2370, 2914)
我们有 2,370 张图像，每张图像有 2,914 个像素。换句话说，这些图像可以被视为 2,914 维空间中的数据点！
让我们展示几张这些图像，以便提醒我们正在处理的内容（参见图 46-11）。
In [17]: fig, ax = plt.subplots(4, 8, subplot_kw=dict(xticks=[], yticks=[]))
for i, axi in enumerate(ax.flat):
axi.imshow(faces.images[i], cmap='gray')
图 46-11. 输入人脸示例
当我们在第四十五章遇到这些数据时，我们的目标基本上是压缩：使用组件从较低维度表示重建输入。
PCA 足够灵活，我们也可以在此上下文中使用它，我们想要绘制 2914 维数据的低维嵌入，以学习图像之间的基本关系。让我们再次查看解释的方差比率，这将为我们提供所需的线性特征数量的概念（参见图 46-12）。
In [18]: from sklearn.decomposition import PCA
model = PCA(100, svd_solver='randomized').fit(faces.data)
plt.plot(np.cumsum(model.explained_variance_ratio_))
plt.xlabel('n components')
plt.ylabel('cumulative variance');
图 46-12. 来自 PCA 投影的累积方差
对于这些数据，我们看到需要近 100 个组件才能保留 90%的方差。这告诉我们，数据在本质上是非常高维的——不能仅用少量组件线性描述。
当情况如此时，非线性流形嵌入如 LLE 和 Isomap 可能会有所帮助。我们可以使用与之前相同的模式在这些人脸上计算 Isomap 嵌入：
In [19]: from sklearn.manifold import Isomap
model = Isomap(n_components=2)
proj = model.fit_transform(faces.data)
proj.shape
Out[19]: (2370, 2)
输出是所有输入图像的二维投影。为了更好地了解投影告诉我们的内容，让我们定义一个函数，该函数将在投影位置输出图像缩略图：
In [20]: from matplotlib import offsetbox
def plot_components(data, model, images=None, ax=None,
thumb_frac=0.05, cmap='gray'):
ax = ax or plt.gca()
proj = model.fit_transform(data)
ax.plot(proj[:, 0], proj[:, 1], '.k')
if images is not None:
min_dist_2 = (thumb_frac * max(proj.max(0) - proj.min(0))) ** 2
shown_images = np.array([2 * proj.max(0)])
for i in range(data.shape[0]):
dist = np.sum((proj[i] - shown_images) ** 2, 1)
if np.min(dist) < min_dist_2:
# don't show points that are too close
continue
shown_images = np.vstack([shown_images, proj[i]])
imagebox = offsetbox.AnnotationBbox(
offsetbox.OffsetImage(images[i], cmap=cmap),
proj[i])
ax.add_artist(imagebox)
现在调用此函数，我们将在图 46-13 中看到结果。
In [21]: fig, ax = plt.subplots(figsize=(10, 10))
plot_components(faces.data,
model=Isomap(n_components=2),
images=faces.images[:, ::2, ::2])
图 46-13. LFW 数据的 Isomap 嵌入
结果很有趣。前两个 Isomap 维度似乎描述了全局图像特征：图像的整体亮度从左到右，以及脸部的一般方向从底部到顶部。这为我们提供了一些数据中一些基本特征的良好视觉指示。
从这里，我们可以继续对这些数据进行分类（也许使用流形特征作为分类算法的输入），就像我们在第四十三章中所做的那样。
示例：可视化数字中的结构
作为使用流形学习进行可视化的另一个示例，让我们看一下 MNIST 手写数字数据集。这类似于我们在第四十四章中看到的数字数据集，但每个图像的像素更多。可以使用 Scikit-Learn 工具从http://openml.org下载它：
In [22]: from sklearn.datasets import fetch_openml
mnist = fetch_openml('mnist_784')
mnist.data.shape
Out[22]: (70000, 784)
数据集由 70,000 个图像组成，每个图像有 784 个像素（即，图像为 28 × 28）。和以前一样，我们可以查看前几个图像（参见图 46-14）。
In [23]: mnist_data = np.asarray(mnist.data)
mnist_target = np.asarray(mnist.target, dtype=int)
fig, ax = plt.subplots(6, 8, subplot_kw=dict(xticks=[], yticks=[]))
for i, axi in enumerate(ax.flat):
axi.imshow(mnist_data[1250 * i].reshape(28, 28), cmap='gray_r')
图 46-14. MNIST 数字示例
这让我们了解到数据集中手写风格的多样性。
让我们在数据集中进行流形学习投影。为了加快速度，我们只使用了 1/30 的数据，大约是 ~2,000 个点（由于流形学习的比例尺度相对较差，我发现几千个样本是一个相对快速探索的良好起点，然后再转向完整的计算）。图 46-15 展示了结果。
In [24]: # Use only 1/30 of the data: full dataset takes a long time!
data = mnist_data[::30]
target = mnist_target[::30]
model = Isomap(n_components=2)
proj = model.fit_transform(data)
plt.scatter(proj[:, 0], proj[:, 1], c=target,
cmap=plt.cm.get_cmap('jet', 10))
plt.colorbar(ticks=range(10))
plt.clim(-0.5, 9.5);
图 46-15. MNIST 数字数据的 Isomap 嵌入
结果散点图展示了数据点之间的一些关系，但有点拥挤。我们可以通过逐个查看单个数字来获得更多见解（参见图 46-16）。
In [25]: # Choose 1/4 of the "1" digits to project
data = mnist_data[mnist_target == 1][::4]
fig, ax = plt.subplots(figsize=(10, 10))
model = Isomap(n_neighbors=5, n_components=2, eigen_solver='dense')
plot_components(data, model, images=data.reshape((-1, 28, 28)),
ax=ax, thumb_frac=0.05, cmap='gray_r')
图 46-16. 仅包含 MNIST 数据集中数字 1 的 Isomap 嵌入
结果让你了解到数据集中数字 1 可以呈现的各种形式。数据沿着投影空间中的一条宽曲线分布，这条曲线似乎跟数字的方向有关。通过观察图上方的部分，你会发现一些带帽子和/或底座的数字 1，尽管它们在数据集中非常稀疏。投影让我们能够识别出存在数据问题的离群值：例如，一些邻近数字的部分可能已经混入提取的图像中。
现在，这本身对于分类数字可能并不有用，但它确实帮助我们了解数据，并可能给我们关于如何继续进行的想法——比如我们可能希望在构建分类管道之前对数据进行预处理。
^(1) 在在线附录中可以找到产生此图的代码。
