Python 机器学习系统构建指南第二版（二）
原文：
annas-archive.org/md5/1799e60b9ee143f87b98f9fc0705a0c9译者：飞龙
第六章：分类 II – 情感分析
对于公司而言，密切监控重大事件的公众反应至关重要，如产品发布或新闻稿。借助 Twitter 的实时访问和用户生成内容的易获取性，现在可以对推文进行情感分类。情感分析有时也称为意见挖掘，它是一个活跃的研究领域，许多公司已经在销售此类服务。由于这表明市场显然存在，我们有动力使用在上一章中构建的分类技术，来构建我们自己的情感分类器。
绘制我们的路线图
推文的情感分析特别困难，因为 Twitter 对字符数的限制为 140 个字符。这导致了特殊的语法、创造性的缩写，并且句子通常不完整。分析句子的典型方法是汇总段落中的情感信息，然后计算文档的整体情感，这种方法在这里行不通。
很显然，我们并不会尝试构建一个最先进的情感分类器。相反，我们的目标是：
-
使用这个场景作为引入另一个分类算法 朴素贝叶斯 的载体
-
解释 词性标注 (POS) 如何工作以及它如何帮助我们
-
展示一些来自 scikit-learn 工具箱的其他小技巧，这些技巧不时会派上用场
获取 Twitter 数据
自然地，我们需要推文及其相应的标签，以判断一条推文是包含积极、消极还是中性情感。在本章中，我们将使用 Niek Sanders 提供的语料库，他手动标注了超过 5,000 条推文，并已授权我们在本章中使用这些数据。
为了遵守 Twitter 的服务条款，我们不会提供任何来自 Twitter 的数据，也不会展示任何真实的推文。相反，我们可以使用 Sanders 的手动标注数据，其中包含推文 ID 及其手动标注的情感，并使用他的脚本 install.py 获取相应的 Twitter 数据。由于该脚本与 Twitter 服务器兼容，因此下载超过 5,000 条推文的数据需要相当长的时间。所以，最好立刻开始运行它。
数据包含四个情感标签：
>>> X, Y = load_sanders_data()
>>> classes = np.unique(Y)
>>> for c in classes: print("#%s: %i" % (c, sum(Y==c)))
#irrelevant: 490
#negative: 487
#neutral: 1952
#positive: 433
在 load_sanders_data() 函数中，我们将不相关和中性标签一起处理为中性，并删除所有非英语推文，最终得到 3,362 条推文。
如果你在这里获得不同的计数，那是因为在此期间，推文可能被删除或设置为私人状态。在这种情况下，你也可能会看到与接下来的章节所展示的数字和图表略有不同。
介绍朴素贝叶斯分类器
朴素贝叶斯可能是最优雅的机器学习算法之一，并且具有实际应用价值。尽管它的名字里有“朴素”二字，但从它的分类表现来看，它并不那么“朴素”。它对无关特征表现出强大的鲁棒性，能够巧妙地忽略它们。它学习速度快，预测也同样迅速。它不需要大量存储。那么，为什么它被称为“朴素”呢？
朴素一词是为了表示朴素贝叶斯所依赖的一个假设，这个假设是特征之间相互独立。实际上，特征之间很少完全独立，这也是现实世界应用中的常见问题。然而，即便假设不成立，朴素贝叶斯在实践中仍能提供很高的准确性。
了解贝叶斯定理
朴素贝叶斯分类的核心不过是记录每个特征给出哪一类的证据。特征的设计决定了使用哪种模型进行学习。所谓的伯努利模型只关注布尔特征：一个单词是否在推文中出现一次或多次并不重要。相反，多项式模型使用单词计数作为特征。为了简化起见，我们将使用伯努利模型来解释如何使用朴素贝叶斯进行情感分析。接下来，我们将使用多项式模型来设置和调优我们的实际分类器。
让我们假设以下变量含义，用以解释朴素贝叶斯：
| 变量 | 含义 |
|---|---|
| 这是推文的类别（正面或负面） | |
| 单词“awesome”至少出现在推文中一次 | |
| 单词“crazy”至少出现在推文中一次 |
在训练过程中，我们学习了朴素贝叶斯模型，这是在已知特征的情况下，某个类别的概率！了解贝叶斯定理。这个概率可以写作！了解贝叶斯定理。
由于我们无法直接估计！了解贝叶斯定理，我们采用了一个技巧，这是贝叶斯发现的：
如果我们将！了解贝叶斯定理替换为“awesome”和“crazy”两个单词的概率，并将！了解贝叶斯定理视为我们的类别！了解贝叶斯定理，我们可以得到这个关系，它帮助我们后来推导数据实例属于指定类别的概率：
这使得我们可以通过其他概率来表达 ：
我们也可以将其描述为：
先验 和 证据 容易确定：
-
是类别
的先验概率，在不知道数据的情况下。我们可以通过简单地计算所有训练数据实例中属于该类别的比例来估计这个量。
-
是特征的证据或概率
和
。
难点在于似然度的计算 。它是描述在已知数据实例的类别是
的情况下，看到特征值
和
的可能性。要估计这个，我们需要做一些思考。
天真假设
从概率理论中，我们还知道以下关系：
然而，单靠这一点并没有太大帮助，因为我们用另一个困难问题（估计 ）来处理一个问题（估计
）。
然而，如果我们天真地假设 和
彼此独立，那么
简化为
，我们可以将其写为：
将一切合并在一起，我们得到一个相当简洁的公式：
有趣的是，尽管从理论上讲，在我们有心情时随意调整假设并不正确，但在这种情况下，实际应用中它的效果出奇地好。
使用朴素贝叶斯进行分类
给定一个新的推文，剩下的唯一任务就是简单地计算概率：
然后选择概率较高的类别 。
由于对于两个类别，分母 是相同的，我们可以忽略它而不改变最终类别。
然而，请注意，我们现在不再计算任何真实的概率。相反，我们在估算给定证据下，哪个类别更有可能。这也是朴素贝叶斯如此稳健的另一个原因：它更关心的是哪个类别更有可能，而不是实际的概率。简而言之，我们可以写：
这只是说明我们正在计算argmax之后的部分，针对所有类别（在我们的例子中是pos和neg），并返回得到最高值的类别。
但是，对于以下示例，我们将坚持使用真实的概率并进行一些计算，以便看看朴素贝叶斯是如何工作的。为了简化起见，我们将假设 Twitter 只允许使用之前提到的两个单词：“awesome”和“crazy”，并且我们已经手动分类了一些推文：
| 推文 | 类别 |
|---|---|
| awesome | 正面推文 |
| awesome | 正面推文 |
| awesome crazy | 正面推文 |
| crazy | 正面推文 |
| crazy | 负面推文 |
| crazy | 负面推文 |
在这个例子中，我们将“crazy”这条推文同时放在正面和负面推文中，以模拟现实世界中经常会遇到的一些模糊情况（例如，“热衷足球”与“疯狂的傻瓜”）。
在这个例子中，我们共有六条推文，其中四条是正面的，二条是负面的，得出的先验概率如下：
这意味着，在不了解推文本身的任何信息的情况下，假设这条推文是正面推文是明智的。
仍然缺少的部分是计算和
，这些是条件概率，分别针对两个特征
和
，并且是基于类别
计算的。
这是通过计算我们见过具体特征的推文数，再除以已被标记为类别的推文总数来得出的。假设我们想知道在已知推文类别为正面的情况下，出现“awesome”的概率，我们将得到：
因为四条正面推文中有三条包含“awesome”这个词。显然，正面推文中不包含“awesome”的概率就是它的逆：
类似地，对于其余情况（省略没有出现该词的推文）：
为了完整性，我们还将计算证据，以便在接下来的示例推文中看到实际概率。对于和
这两个具体值，我们可以按如下方式计算证据：
这导致了以下值：
现在我们有了所有的数据来分类新的推文。剩下的工作就是解析推文并提取特征：
| 推文 | 类别概率 | 分类 | ||
|---|---|---|---|---|
| "极棒" | 1 | 0 | 正面 | |
| "疯狂" | 0 | 1 | 负面 | |
| "极棒 疯狂" | 1 | 1 | 正面 |
到目前为止，一切顺利。对于简单的推文分类，似乎能够给推文分配正确的标签。然而，问题仍然是，我们应该如何处理那些在训练语料库中没有出现过的词？毕竟，使用前面的公式，新词总是会被分配零的概率。
考虑到未见过的词和其他异常情况
当我们之前计算概率时，实际上我们是在自欺欺人。我们并没有计算真正的概率，而只是通过分数来进行粗略的估算。我们假设训练语料库会告诉我们关于真实概率的全部真相，但事实并非如此。仅仅六条推文的语料库显然不能告诉我们所有关于曾经写过的推文的信息。例如，肯定有包含“文本”一词的推文，只是我们从未见过它们。显然，我们的估算非常粗略，我们应该对此加以考虑。在实践中，这通常通过所谓的加一平滑来实现。
提示
加一平滑有时也被称为加性平滑或拉普拉斯平滑。注意，拉普拉斯平滑与拉普拉斯算子平滑无关，后者是与多边形网格的平滑相关的。如果我们不是通过1来平滑，而是通过可调参数alpha<0来平滑，那就叫做 Lidstone 平滑。
这是一种非常简单的技术，它为所有特征出现次数加 1。其基本假设是，即使我们在整个语料库中没有见过某个单词，也有可能是我们的推文样本刚好没有包含这个单词。所以，通过加一平滑，我们假装每个出现的单词比实际出现的次数多见了一次。这意味着，我们现在计算的不是，而是
。
为什么我们在分母中加 2？因为我们有两个特征：“awesome”和“crazy”的出现。由于每个特征加 1，我们必须确保最终结果仍然是一个概率。事实上，我们得到的总概率为 1：
计算算术下溢
还有一个障碍。实际上，我们处理的概率要比在玩具示例中遇到的要小得多。通常，我们也有比仅仅两个特征更多的特征，而这些特征需要相互相乘。这将很快导致 NumPy 提供的精度不再足够：
>>> import numpy as np
>>> np.set_printoptions(precision=20) # tell numpy to print out more digits (default is 8)
>>> np.array([2.48E-324])
array([ 4.94065645841246544177e-324])
>>> np.array([2.47E-324])
array([ 0.])
那么，2.47E-324 这样的数字到底有多可能出现呢？为了回答这个问题，我们只需要想象一个条件概率为 0.0001，然后将其中 65 个概率相乘（意味着我们有 65 个低概率特征值），你就会遇到算术下溢：
>>> x = 0.00001
>>> x**64 # still fine
1e-320
>>> x**65 # ouch
0.0
Python 中的浮点数通常使用 C 中的 double 类型实现。要检查你的平台是否是这样，你可以通过以下方式进行确认：
>>> import sys
>>> sys.float_info
sys.float_info(max=1.7976931348623157e+308, max_exp=1024, max_10_exp=308, min=2.2250738585072014e-308, min_exp=-1021, min_10_exp=-307, dig=15, mant_dig=53, epsilon=2.220446049250313e-16, radix=2, rounds=1)
为了缓解这个问题，可以切换到诸如mpmath（code.google.com/p/mpmath/）这样的数学库，它们支持任意精度。然而，它们的速度不够快，无法替代 NumPy。
幸运的是，有一个更好的方法可以解决这个问题，这与我们可能还记得的一个学校里的好关系有关：
如果我们将其应用到我们的案例中，我们得到如下结果：
由于概率位于 0 和 1 之间，概率的对数则位于-∞和 0 之间。别为这个困扰。较大的数值仍然是正确类别的更强指示——只是现在它们变成了负数。
但有一个注意事项：公式的分子部分（即分数上方的部分）实际上并没有对数。我们只有概率的乘积。幸运的是，在我们的案例中，我们并不关心概率的实际值。我们只是想知道哪个类别具有最高的后验概率。幸运的是，如果我们发现，那么我们也总是能得到
。
快速查看前面的图表可以发现曲线是单调递增的，即从左到右时，曲线永远不会下降。所以让我们把这个带入前面提到的公式：
这将最终得到适用于两个特征的公式，它将为我们提供最佳类别，适用于我们在实践中看到的实际数据：
当然，只有两个特征的话，我们不会非常成功，所以，让我们重写代码以允许任意数量的特征：
就这样，我们准备好使用来自 scikit-learn 工具包的第一个分类器。
如前所述，我们刚刚学习了 Naïve Bayes 的伯努利模型。与布尔特征不同，我们也可以使用单词出现次数，即多项式模型。由于这个方法提供了更多的信息，通常也能带来更好的性能，因此我们会在实际数据中使用这个模型。不过，注意的是，底层的公式会有一些变化。不过，不用担心，Naïve Bayes 的基本原理依然不变。
创建我们的第一个分类器并进行调优
Naïve Bayes 分类器位于 sklearn.naive_bayes 包中。这里有不同种类的 Naïve Bayes 分类器：
-
GaussianNB：这个分类器假设特征呈正态分布（高斯分布）。它的一个应用场景可能是根据身高和体宽来判断性别。在我们的例子中，我们有推文文本，从中提取词频。这些显然不是高斯分布的。 -
MultinomialNB：这个分类器假设特征是出现次数，这正是我们接下来使用的情况，因为我们将在推文中使用词频作为特征。实际上，这个分类器也能很好地处理 TF-IDF 向量。 -
BernoulliNB：这个分类器与MultinomialNB类似，但更适合处理二进制单词出现与否，而非词频。
由于我们将主要关注词频，因此对于我们的目的来说，MultinomialNB 分类器最为合适。
先解决一个简单的问题
正如我们在查看推文数据时所看到的，推文不仅仅是正面或负面的。实际上，大多数推文不包含任何情感，而是中立或无关的，包含例如原始信息（例如，“新书：构建机器学习……http://link”）。这导致了四个类别。为了避免任务过于复杂，现在我们只关注正面和负面的推文。
>>> # first create a Boolean list having true for tweets
>>> # that are either positive or negative
>>> pos_neg_idx = np.logical_or(Y=="positive", Y=="negative")
>>> # now use that index to filter the data and the labels
>>> X = X[pos_neg_idx]
>>> Y = Y[pos_neg_idx]
>>> # finally convert the labels themselves into Boolean
>>> Y = Y=="positive"
现在，我们在X中有原始推文文本，在Y中有二元分类，0表示负面推文，1表示正面推文。
我们刚刚提到过，我们将使用单词出现次数作为特征。但我们不会直接使用它们的原始形式，而是使用我们的强力工具TfidfVectorizer，将原始推文文本转换为 TF-IDF 特征值，然后将其与标签一起用于训练我们的第一个分类器。为了方便，我们将使用Pipeline类，它允许我们将向量化器和分类器连接在一起，并提供相同的接口：
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.naive_bayes import MultinomialNB
from sklearn.pipeline import Pipeline
def create_ngram_model():
tfidf_ngrams = TfidfVectorizer(ngram_range=(1, 3), analyzer="word", binary=False)
clf = MultinomialNB()
return Pipeline([('vect', tfidf_ngrams), ('clf', clf)])
create_ngram_model()返回的Pipeline实例现在可以像普通分类器一样用于拟合和预测。
由于我们没有那么多数据，应该进行交叉验证。然而，这一次，我们不会使用KFold（它将数据划分为连续的折叠），而是使用ShuffleSplit。它会将数据打乱，但不会阻止相同的数据实例出现在多个折叠中。对于每个折叠，我们会跟踪精准率-召回率曲线下的面积和准确度。
为了保持实验的灵活性，我们将一切封装在一个train_model()函数中，该函数接受一个创建分类器的函数作为参数。
from sklearn.metrics import precision_recall_curve, auc
from sklearn.cross_validation import ShuffleSplit
def train_model(clf_factory, X, Y):
# setting random_state to get deterministic behavior
cv = ShuffleSplit(n=len(X), n_iter=10, test_size=0.3, random_state=0)
scores = []
pr_scores = []
for train, test in cv:
X_train, y_train = X[train], Y[train]
X_test, y_test = X[test], Y[test]
clf = clf_factory()
clf.fit(X_train, y_train)
train_score = clf.score(X_train, y_train)
test_score = clf.score(X_test, y_test)
scores.append(test_score)
proba = clf.predict_proba(X_test)
precision, recall, pr_thresholds = precision_recall_curve(y_test, proba[:,1])
pr_scores.append(auc(recall, precision))
summary = (np.mean(scores), np.std(scores),
np.mean(pr_scores), np.std(pr_scores))
print("%.3f\t%.3f\t%.3f\t%.3f" % summary)
将所有内容组合起来，我们可以训练我们的第一个模型：
>>> X, Y = load_sanders_data()
>>> pos_neg_idx = np.logical_or(Y=="positive", Y=="negative")
>>> X = X[pos_neg_idx]
>>> Y = Y[pos_neg_idx]
>>> Y = Y=="positive"
>>> train_model(create_ngram_model, X, Y)
0.788 0.024 0.882 0.036
使用朴素贝叶斯和向量化的 TF-IDF 三元组特征，我们的第一次尝试得到了 78.8%的准确率和 88.2%的平均 P/R AUC。当我们查看中位数的 P/R 图表（即表现最接近平均水平的训练/测试拆分）时，它表现出比我们在上一章看到的图表更加令人鼓舞的行为。
一开始，结果相当令人鼓舞。当我们意识到在情感分类任务中，100%的准确率可能永远无法实现时，结果就变得更加令人印象深刻。对于某些推文，甚至人类在分类标签上也往往达不到一致。
使用所有类别
我们再次简化了任务，因为我们只使用了正面或负面的推文。这意味着，我们假设有一个完美的分类器，首先对推文是否包含情感进行分类，然后将结果传递给我们的朴素贝叶斯分类器。
那么，如果我们还对推文是否包含任何情感进行分类，我们的表现如何呢？为了找出答案，我们首先编写一个便捷函数，返回一个修改后的类别数组，提供我们希望解释为正面的情感列表：
def tweak_labels(Y, pos_sent_list):
pos = Y==pos_sent_list[0]
for sent_label in pos_sent_list[1:]:
pos |= Y==sent_label
Y = np.zeros(Y.shape[0])
Y[pos] = 1
Y = Y.astype(int)
return Y
请注意，我们现在谈论的是两个不同的正面情感。推文的情感可以是正面的，这与训练数据的类别不同。例如，如果我们想要了解我们如何区分具有情感的推文和中立推文，我们可以这样做：
>>> Y = tweak_labels(Y, ["positive", "negative"])
在Y中，对于所有正面或负面的推文，我们现在有1（正类），对于中立和无关的推文，我们有0（负类）。
>>> train_model(create_ngram_model, X, Y, plot=True)
0.750 0.012 0.659 0.023
看一下以下的图表：
如预期的那样，P/R AUC 显著下降，目前只有 66%。准确率仍然很高，但这仅仅是因为我们有一个高度不平衡的数据集。在 3,362 条推文中，只有 920 条是正面或负面的，占约 27%。这意味着，如果我们创建一个总是将推文分类为不包含任何情感的分类器，我们就已经有 73%的准确率了。这是另一个在训练和测试数据不平衡时，总是查看精度和召回率的案例。
那么，朴素贝叶斯分类器在将正面推文与其他推文以及负面推文与其他推文进行分类时表现如何呢？一句话：差。
== Pos vs. rest ==
0.873 0.009 0.305 0.026
== Neg vs. rest ==
0.861 0.006 0.497 0.026
如果你问我，这几乎是不可用的。看下面的 P/R 曲线，我们也会发现没有可用的精度/召回率权衡，就像我们在上一章做的那样：
调整分类器的参数
当然，我们还没有充分探索当前的设置，应该进行更多的调查。大致有两个领域，我们可以调整参数：TfidfVectorizer和MultinomialNB。由于我们对哪个领域的探索没有直觉，让我们尝试分配参数值。
我们将首先查看TfidfVectorizer参数：
-
使用不同的 NGram 设置：
-
单字（1,1）
-
单字和双字（1,2）
-
单字、双字和三字（1,3）
-
-
调整
min_df：1 或 2 -
使用
use_idf和smooth_idf探索 IDF 在 TF-IDF 中的影响：False或True -
是否移除停用词，通过将
stop_words设置为english或None -
是否使用词频的对数值（
sublinear_tf） -
是否跟踪词频，还是仅跟踪词是否出现，通过将
binary设置为True或False
现在我们将查看MultinomialNB分类器：
-
通过设置
alpha来选择平滑方法：-
加一平滑或拉普拉斯平滑：1
-
Lidstone 平滑：0.01，0.05，0.1，或 0.5
-
无平滑：0
-
一种简单的方法是为所有合理的探索值训练一个分类器，同时保持其他参数不变并检查分类器的结果。由于我们不知道这些参数是否相互影响，正确的方法是为所有可能的参数值组合训练一个分类器。显然，这对我们来说太繁琐了。
因为这种参数探索在机器学习任务中经常发生，scikit-learn 为此提供了一个专门的类，叫做GridSearchCV。它接受一个估算器（具有分类器类似接口的实例），在我们这个案例中将是Pipeline实例，并且接受一个包含参数及其潜在值的字典。
GridSearchCV期望字典的键遵循特定的格式，以便能够为正确的估算器设置参数。格式如下：
<estimator>__<subestimator>__...__<param_name>
例如，如果我们想指定要探索的TfidfVectorizer（在Pipeline描述中称为vect）的min_df参数的期望值，我们需要这样说：
param_grid={"vect__ngram_range"=[(1, 1), (1, 2), (1, 3)]}
这将告诉GridSearchCV尝试将一元组到三元组作为TfidfVectorizer的ngram_range参数的候选值。
然后，它使用所有可能的参数值组合训练估算器。在这里，我们确保它在训练数据的随机样本上进行训练，使用ShuffleSplit，它生成随机的训练/测试分割迭代器。最后，它提供最佳估算器，作为成员变量best_estimator_。
由于我们希望将返回的最佳分类器与当前的最佳分类器进行比较，我们需要以相同的方式评估它。因此，我们可以通过cv参数传递ShuffleSplit实例（因此，GridSearchCV中的CV）。
最后一块缺失的部分是定义GridSearchCV应如何确定最佳估算器。这可以通过提供期望的评分函数来完成（惊讶！）传递给score_func参数。我们可以自己编写一个，或者从sklearn.metrics包中选择一个。我们肯定不应该使用metric.accuracy，因为我们的数据存在类别不平衡（包含情感的推文远少于中立推文）。相反，我们希望对两个类别都有良好的精度和召回率，即具有情感的推文和没有正面或负面观点的推文。一种结合精度和召回率的评价指标是所谓的F 值，它在metrics.f1_score中实现：
将所有内容整合后，我们得到以下代码：
from sklearn.grid_search import GridSearchCV
from sklearn.metrics import f1_score
def grid_search_model(clf_factory, X, Y):
cv = ShuffleSplit(
n=len(X), n_iter=10, test_size=0.3,random_state=0)
param_grid = dict(vect__ngram_range=[(1, 1), (1, 2), (1, 3)],
vect__min_df=[1, 2],
vect__stop_words=[None, "english"],
vect__smooth_idf=[False, True],
vect__use_idf=[False, True],
vect__sublinear_tf=[False, True],
vect__binary=[False, True],
clf__alpha=[0, 0.01, 0.05, 0.1, 0.5, 1],
)
grid_search = GridSearchCV(clf_factory(),
param_grid=param_grid,
cv=cv,
score_func=f1_score,
verbose=10)
grid_search.fit(X, Y)
return grid_search.best_estimator_
执行此操作时我们需要耐心：
clf = grid_search_model(create_ngram_model, X, Y)
print(clf)
由于我们刚刚请求了一个参数，覆盖![调整分类器的参数]的参数组合，每种组合都将在 10 折交叉验证中进行训练：
... waiting some hours ...
Pipeline(clf=MultinomialNB(
alpha=0.01, class_weight=None, fit_prior=True),
clf__alpha=0.01,
clf__class_weight=None,
clf__fit_prior=True,
vect=TfidfVectorizer(
analyzer=word, binary=False,
charset=utf-8, charset_error=strict,
dtype=<type 'long'>,input=content,
lowercase=True, max_df=1.0,
max_features=None, max_n=None,
min_df=1, min_n=None, ngram_range=(1, 2),
norm=l2, preprocessor=None, smooth_idf=False,
stop_words=None,strip_accents=None,
sublinear_tf=True,token_pattern=(?u)\b\w\w+\b,
token_processor=None, tokenizer=None,
use_idf=False, vocabulary=None),
vect__analyzer=word, vect__binary=False,
vect__charset=utf-8,
vect__charset_error=strict,
vect__dtype=<type 'long'>,
vect__input=content, vect__lowercase=True,
vect__max_df=1.0,vect__max_features=None,
vect__max_n=None, vect__min_df=1,
vect__min_n=None, vect__ngram_range=(1, 2),
vect__norm=l2, vect__preprocessor=None,
vect__smooth_idf=False, vect__stop_words=None,
vect__strip_accents=None, vect__sublinear_tf=True,
vect__token_pattern=(?u)\b\w\w+\b,
vect__token_processor=None, vect__tokenizer=None,
vect__use_idf=False, vect__vocabulary=None)
0.795 0.007 0.702 0.028
最佳估算器确实将 P/R AUC 提高了近 3.3 个百分点，现在为 70.2，设置如前面代码所示。
此外，针对正面推文与其他推文的比较，负面推文与其他推文的比较，若我们配置向量化器和分类器使用刚刚发现的参数，结果会得到显著改善：
== Pos vs. rest ==
0.889 0.010 0.509 0.041
== Neg vs. rest ==
0.886 0.007 0.615 0.035
看一下以下图表：
的确，P/R 曲线看起来好多了（注意，这些图是来自于不同折分类器的中位数，因此 AUC 值略有差异）。不过，我们可能仍然不会使用这些分类器。是时候换点新的东西了…
清洗推文
新的约束导致了新的形式。Twitter 在这方面也不例外。因为文本必须适应 140 个字符，人们自然会发展出新的语言简写，以用更少的字符表达相同的意思。到目前为止，我们忽略了各种各样的表情符号和缩写。让我们看看在考虑到这些因素后，我们能做出多大的改进。为了这个目标，我们将需要提供我们自己的preprocessor()来处理TfidfVectorizer。
首先，我们在字典中定义一系列常见的表情符号及其替换词。虽然我们可以找到更多不同的替换词，但我们选择明显的正面或负面词汇来帮助分类器：
emo_repl = {
# positive emoticons
"<3": " good ",
":d": " good ", # :D in lower case
":dd": " good ", # :DD in lower case
"8)": " good ",
":-)": " good ",
":)": " good ",
";)": " good ",
"(-:": " good ",
"(:": " good ",
# negative emoticons:
":/": " bad ",
":>": " sad ",
":')": " sad ",
":-(": " bad ",
":(": " bad ",
":S": " bad ",
":-S": " bad ",
}
# make sure that e.g. :dd is replaced before :d
emo_repl_order = [k for (k_len,k) in reversed(sorted([(len(k),k) for k in emo_repl.keys()]))]
然后，我们将缩写定义为正则表达式，并给出它们的扩展（\b表示单词边界）：
re_repl = {
r"\br\b": "are",
r"\bu\b": "you",
r"\bhaha\b": "ha",
r"\bhahaha\b": "ha",
r"\bdon't\b": "do not",
r"\bdoesn't\b": "does not",
r"\bdidn't\b": "did not",
r"\bhasn't\b": "has not",
r"\bhaven't\b": "have not",
r"\bhadn't\b": "had not",
r"\bwon't\b": "will not",
r"\bwouldn't\b": "would not",
r"\bcan't\b": "can not",
r"\bcannot\b": "can not",
}
def create_ngram_model(params=None):
def preprocessor(tweet):
tweet = tweet.lower()
for k in emo_repl_order:
tweet = tweet.replace(k, emo_repl[k])
for r, repl in re_repl.items():
tweet = re.sub(r, repl, tweet)
return tweet
tfidf_ngrams = TfidfVectorizer(preprocessor=preprocessor,
analyzer="word")
# ...
当然，这里可以使用更多的缩写。但仅凭这一有限的集合，我们已经能在情感与非情感的分类上提高了半个点，目前准确率为 70.7%：
== Pos vs. neg ==
0.808 0.024 0.885 0.029
== Pos/neg vs. irrelevant/neutral ==
0.793 0.010 0.685 0.024
== Pos vs. rest ==
0.890 0.011 0.517 0.041
== Neg vs. rest ==
0.886 0.006 0.624 0.033
考虑单词类型
到目前为止，我们的期望是仅仅使用单词本身（采用袋装词汇方法）就足够了。然而，仅凭我们的直觉，可能中立的推文包含更多名词，而积极或消极的推文则更加多彩，包含更多形容词和动词。那么，如果我们也利用推文中的语言学信息呢？如果我们能找出推文中有多少单词是名词、动词、形容词等，分类器也许能够将这一信息纳入考虑。
确定单词类型
这就是词性标注，或 POS 标注的核心。POS 标注器解析一个完整的句子，目的是将其排列成依存树，其中每个节点对应一个单词，父子关系决定了该单词依赖哪个单词。通过这棵树，它可以做出更为明智的决策，例如判断“book”是名词（“这是一本好书。”）还是动词（“你能帮我订机票吗？”）。
你可能已经猜到，NLTK 在这个领域也将发挥作用。的确，它提供了各种各样的解析器和标注器。我们将使用的 POS 标注器nltk.pos_tag()实际上是一个完整的分类器，经过使用 Penn Treebank 项目中的人工标注句子训练而成（www.cis.upenn.edu/~treebank）。它的输入是一个单词标记列表，输出是一个包含元组的列表，每个元组包含原句的部分内容及其词性标注。
>>> import nltk
>>> nltk.pos_tag(nltk.word_tokenize("This is a good book."))
[('This', 'DT'), ('is', 'VBZ'), ('a', 'DT'), ('good', 'JJ'), ('book', 'NN'), ('.', '.')]
>>> nltk.pos_tag(nltk.word_tokenize("Could you please book the flight?"))
[('Could', 'MD'), ('you', 'PRP'), ('please', 'VB'), ('book', 'NN'), ('the', 'DT'), ('flight', 'NN'), ('?', '.')]
词性标注缩写来自于 Penn Treebank（改编自www.anc.org/OANC/penn.html）：
| POS 标签 | 描述 | 示例 |
|---|---|---|
| CC | 并列连词 | 或者 |
| CD | 基数词 | 2，第二 |
| DT | 限定词 | 这 |
| EX | 存在性 there | 那里 有 |
| FW | 外来词 | 幼儿园 |
| IN | 介词/从属连词 | 在，的，像 |
| JJ | 形容词 | 酷 |
| JJR | 形容词，比较级 | 更酷 |
| JJS | 形容词，最高级 | 最酷的 |
| LS | 列表标记 | 1) |
| MD | 情态动词 | 可以，将 |
| NN | 名词，单数或不可数名词 | 书 |
| NNS | 名词复数 | 书籍 |
| NNP | 专有名词，单数 | Sean |
| NNPS | 专有名词，复数 | 维京人 |
| PDT | 预定限定词 | 两个男孩 |
| POS | 所有格结尾 | 朋友的 |
| PRP | 人称代词 | 我，他，它 |
| PRP$ | 所有格代词 | 我的，他的 |
| RB | 副词 | 然而，通常，当然，这里，好 |
| RBR | 副词，比较级 | 更好 |
| RBS | 副词，最高级 | 最好 |
| RP | 小品词 | 放弃 |
| TO | 到 | 到 去，到 他那里 |
| UH | 感叹词 | 嗯嗯 |
| VB | 动词，基本形式 | 拿 |
| VBD | 动词，过去时 | 拿 |
| VBG | 动词，动名词/现在分词 | 拿 |
| VBN | 动词，过去分词 | 已拿 |
| VBP | 动词，单数现在时，非第三人称 | 拿 |
| VBZ | 动词，第三人称单数现在时 | 拿 |
| WDT | 疑问限定词 | 哪个 |
| WP | 疑问代词 | 谁，什么 |
| WP$ | 所有格疑问代词 | 谁的 |
| WRB | 疑问副词 | 哪里，什么时候 |
通过这些标签，从pos_tag()的输出中过滤出所需的标签是相当容易的。我们只需数出所有标签以NN开头的名词，VB开头的动词，JJ开头的形容词，和RB开头的副词。
使用 SentiWordNet 成功作弊
虽然前面提到的语言学信息最有可能帮助我们，但我们可以做得更好来收获它：SentiWordNet (sentiwordnet.isti.cnr.it)。简单来说，它是一个 13MB 的文件，为大多数英语单词分配了正负值。更复杂的说，对于每个同义词集合，它记录了正面和负面的情感值。以下是一些例子：
| POS | ID | PosScore | NegScore | 同义词集合 | 描述 |
|---|---|---|---|---|---|
| a | 00311354 | 0.25 | 0.125 | 勤奋的#1 | 以细心和努力为特征；"做了一个勤奋的尝试来修理电视机" |
| a | 00311663 | 0 | 0.5 | 粗心#1 | 特征是缺乏关注、考虑、预见或彻底性；不小心… |
| n | 03563710 | 0 | 0 | 移植物#1 | 永久性地植入组织的假体 |
| v | 00362128 | 0 | 0 | kink#2 曲线#5 卷曲#1 | 形成弯曲、曲线或扭结；"雪茄烟雾在天花板上卷曲" |
利用POS列中的信息，我们能够区分名词"book"和动词"book"。PosScore和NegScore一起帮助我们确定单词的中立性，即 1-PosScore-NegScore。SynsetTerms列出该集合中所有的同义词。对于我们的任务来说，我们可以忽略ID和Description列。
同义词集术语后面会有一个数字，因为有些词会出现在多个同义词集中。例如，"fantasize"表达了两种截然不同的意思，因此会导致不同的分数：
| POS | ID | PosScore | NegScore | SynsetTerms | Description |
|---|---|---|---|---|---|
| v | 01636859 | 0.375 | 0 | fantasize#2 fantasise#2 | 在脑海中描绘；"他在幻想理想的妻子" |
| v | 01637368 | 0 | 0.125 | fantasy#1 fantasize#1 fantasise#1 | 沉溺于幻想；"他说他计划创办自己的公司时就是在幻想" |
为了找出应该使用哪些同义词集，我们需要真正理解推文的含义，这超出了本章的范围。专注于这一挑战的研究领域叫做词义消歧。在我们的任务中，我们走捷径，简单地计算所有同义词集中的分数平均值，某个词项出现在这些同义词集中。例如，"fantasize"的PosScore为 0.1875，NegScore为 0.0625。
以下函数load_sent_word_net()会为我们完成所有这些工作，并返回一个字典，其中键是形如词性/词的字符串，例如 n/implant，值为正负分数：
import csv, collections
def load_sent_word_net():
# making our life easier by using a dictionary that
# automatically creates an empty list whenever we access
# a not yet existing key
sent_scores = collections.defaultdict(list)
with open(os.path.join(DATA_DIR, SentiWordNet_3.0.0_20130122.txt"), "r") as csvfile:
reader = csv.reader(csvfile, delimiter='\t',
quotechar='"')
for line in reader:
if line[0].startswith("#"):
continue
if len(line)==1:
continue
POS, ID, PosScore, NegScore, SynsetTerms, Gloss = line
if len(POS)==0 or len(ID)==0:
continue
for term in SynsetTerms.split(" "):
# drop number at the end of every term
term = term.split("#")[0]
term = term.replace("-", " ").replace("_", " ")
key = "%s/%s"%(POS, term.split("#")[0])
sent_scores[key].append((float(PosScore),
float(NegScore)))
for key, value in sent_scores.items():
sent_scores[key] = np.mean(value, axis=0)
return sent_scores
我们的第一个估算器
现在，我们已经准备好创建我们自己的第一个特征提取器。最方便的方式是从BaseEstimator继承。它要求我们实现以下三个方法：
-
get_feature_names(): 该函数返回一个字符串列表，表示我们将在transform()中返回的特征。 -
fit(document, y=None): 由于我们没有实现分类器，因此可以忽略这一点，直接返回 self。 -
transform(documents): 该方法返回numpy.array()，包含形状为(len(documents), len(get_feature_names))的数组。也就是说，对于documents中的每个文档，它必须为get_feature_names()中的每个特征名称返回一个值。
下面是实现：
sent_word_net = load_sent_word_net()
class LinguisticVectorizer(BaseEstimator):
def get_feature_names(self):
return np.array(['sent_neut', 'sent_pos', 'sent_neg',
'nouns', 'adjectives', 'verbs', 'adverbs',
'allcaps', 'exclamation', 'question', 'hashtag', 'mentioning'])
# we don't fit here but need to return the reference
# so that it can be used like fit(d).transform(d)
def fit(self, documents, y=None):
return self
def _get_sentiments(self, d):
sent = tuple(d.split())
tagged = nltk.pos_tag(sent)
pos_vals = []
neg_vals = []
nouns = 0.
adjectives = 0.
verbs = 0.
adverbs = 0.
for w,t in tagged:
p, n = 0,0
sent_pos_type = None
if t.startswith("NN"):
sent_pos_type = "n"
nouns += 1
elif t.startswith("JJ"):
sent_pos_type = "a"
adjectives += 1
elif t.startswith("VB"):
sent_pos_type = "v"
verbs += 1
elif t.startswith("RB"):
sent_pos_type = "r"
adverbs += 1
if sent_pos_type is not None:
sent_word = "%s/%s" % (sent_pos_type, w)
if sent_word in sent_word_net:
p,n = sent_word_net[sent_word]
pos_vals.append(p)
neg_vals.append(n)
l = len(sent)
avg_pos_val = np.mean(pos_vals)
avg_neg_val = np.mean(neg_vals)
return [1-avg_pos_val-avg_neg_val, avg_pos_val, avg_neg_val,
nouns/l, adjectives/l, verbs/l, adverbs/l]
def transform(self, documents):
obj_val, pos_val, neg_val, nouns, adjectives, \
verbs, adverbs = np.array([self._get_sentiments(d) \
for d in documents]).T
allcaps = []
exclamation = []
question = []
hashtag = []
mentioning = []
for d in documents:
allcaps.append(np.sum([t.isupper() \
for t in d.split() if len(t)>2]))
exclamation.append(d.count("!"))
question.append(d.count("?"))
hashtag.append(d.count("#"))
mentioning.append(d.count("@"))
result = np.array([obj_val, pos_val, neg_val, nouns, adjectives, verbs, adverbs, allcaps, exclamation, question,
hashtag, mentioning]).T
return result
整合所有内容
然而，仅仅使用这些语言特征而不考虑单词本身，是无法带我们走得很远的。因此，我们必须将TfidfVectorizer参数与语言特征结合起来。可以使用 scikit-learn 的FeatureUnion类来实现这一点。它的初始化方式与Pipeline相同；但是，FeatureUnion不同于按顺序评估估算器并将每个估算器的输出传递给下一个估算器，它是并行处理的，并在之后将输出向量合并。
def create_union_model(params=None):
def preprocessor(tweet):
tweet = tweet.lower()
for k in emo_repl_order:
tweet = tweet.replace(k, emo_repl[k])
for r, repl in re_repl.items():
tweet = re.sub(r, repl, tweet)
return tweet.replace("-", " ").replace("_", " ")
tfidf_ngrams = TfidfVectorizer(preprocessor=preprocessor, analyzer="word")
ling_stats = LinguisticVectorizer()
all_features = FeatureUnion([('ling', ling_stats), ('tfidf', tfidf_ngrams)])
clf = MultinomialNB()
pipeline = Pipeline([('all', all_features), ('clf', clf)])
if params:
pipeline.set_params(**params)
return pipeline
在结合特征提取器上进行训练和测试，可以使正负样本的平均 P/R AUC 提高 0.4 个百分点：
== Pos vs. neg ==
0.810 0.023 0.890 0.025
== Pos/neg vs. irrelevant/neutral ==
0.791 0.007 0.691 0.022
== Pos vs. rest ==
0.890 0.011 0.529 0.035
== Neg vs. rest ==
0.883 0.007 0.617 0.033
time spent: 214.12578797340393
根据这些结果，我们可能不想使用正向对比休息和负向对比休息的分类器，而是首先使用一个分类器来确定推文是否包含情感（正向/负向与无关/中立），然后，如果包含情感，再使用正向与负向分类器来确定具体情感。
总结
恭喜你坚持到最后！我们一起学习了朴素贝叶斯如何工作，以及它为什么并不那么“天真”。尤其是在训练集数据不足以学习类别概率空间中的所有细分时，朴素贝叶斯在泛化方面表现出色。我们学习了如何将其应用于推文，并且清理粗糙的推文文本帮助很大。最后，我们意识到适当的“作弊”（前提是我们已经做了足够的工作）是可以接受的，尤其是当它带来分类器性能的提升时，正如我们在使用SentiWordNet时所经历的那样。
在下一章中，我们将讨论回归分析。
第七章：回归
你可能在高中数学课上学过回归分析。你学到的具体方法可能就是所谓的普通最小二乘（OLS）回归。这个已有 200 年历史的技术计算速度很快，并且可以应用于许多现实世界的问题。本章将首先回顾这一方法，并展示它是如何在 scikit-learn 中实现的。
然而，对于某些问题，这种方法是不够的。特别是当我们拥有许多特征时，这种方法就显得不够用了，尤其是在特征数比数据点数更多的情况下，这种方法完全无法工作。在这些情况下，我们需要更先进的方法。这些方法非常现代，过去十年取得了重大进展，名字如 Lasso、Ridge 或 ElasticNets。我们将详细介绍这些方法，它们也可以在 scikit-learn 中使用。
使用回归预测房价
我们从一个简单的问题开始——预测波士顿的房价；这是一个可以使用公开数据集的问题。我们得到了一些人口统计学和地理信息，如犯罪率或邻里的师生比例。目标是预测某个地区房屋的中位数价值。像往常一样，我们有一些训练数据，答案是已知的。
这是 scikit-learn 自带的内置数据集之一，因此非常容易将数据加载到内存中：
>>> from sklearn.datasets import load_boston
>>> boston = load_boston()
boston 对象包含多个属性；特别地，boston.data 包含输入数据，boston.target 包含房价。
我们将从一个简单的一维回归开始，尝试用一个特征来回归房价，这个特征是每个住宅区域的平均房间数，存储在位置 5（你可以查阅 boston.DESCR 和 boston.feature_names 以获取数据的详细信息）：
>>> from matplotlib import pyplot as plt
>>> plt.scatter(boston.data[:,5], boston.target, color='r')
boston.target 属性包含了平均房价（我们的目标变量）。我们可以使用你可能在高中时学到的标准最小二乘回归。我们第一次尝试的代码如下：
>>> from sklearn.linear_model import LinearRegression
>>> lr = LinearRegression()
我们从 sklearn.linear_model 模块导入 LinearRegression 并构造一个 LinearRegression 对象。这个对象的行为与我们之前使用的 scikit-learn 分类器对象类似。
>>> x = boston.data[:,5]
>>> y = boston.target
>>> x = np.transpose(np.atleast_2d(x))
>>> lr.fit(x, y)
>>> y_predicted = lr.predict(x)
这段代码中唯一不显而易见的行是对 np.atleast_2d 的调用，它将 x 从一维数组转换为二维数组。这个转换是必要的，因为 fit 方法期望其第一个参数是一个二维数组。最后，为了确保维度正确，我们需要对这个数组进行转置。
请注意，我们在 LinearRegression 对象上调用了名为 fit 和 predict 的方法，就像之前使用分类器对象时一样，尽管现在我们执行的是回归操作。这种 API 的一致性是 scikit-learn 的一个优点。
上图展示了所有的点（以点表示）和我们的拟合曲线（实线）。我们可以看到，视觉效果很好，除了少数几个离群点。
然而，从理想角度来看，我们希望定量衡量拟合的好坏。这对于能够比较不同方法非常关键。为此，我们可以测量预测值与真实值之间的接近程度。为此，我们可以使用sklearn.metrics模块中的mean_squared_error函数：
>>> from sklearn.metrics import mean_squared_error
该函数有两个参数，真实值和预测值，如下所示：
>>> mse = mean_squared_error(y, lr.predict(x))
>>> print("Mean squared error (of training data): {:.3}".format(mse))
Mean squared error (of training data): 58.4
这个值有时可能很难解读，最好取其平方根，得到均方根误差（RMSE）：
>>> rmse = np.sqrt(mse)
>>> print("RMSE (of training data): {:.3}".format(rmse))
RMSE (of training data): 6.6
使用 RMSE 的一个优势是，我们可以通过将其乘以二，快速获得误差的粗略估计。在我们的案例中，我们可以预计估算的价格与实际价格之间的差异最多为 1.3 万美元。
提示
均方根误差与预测
均方根误差大致对应于标准差的估计。由于大多数数据距离均值最多两个标准差，因此我们可以将 RMSE 翻倍来获得一个粗略的置信区间。如果误差符合正态分布，这种方法完全有效，但即使误差不是正态分布，它通常也大致正确。
像 6.6 这样的数字仍然很难直观理解。这是一个好的预测吗？回答这个问题的一种可能方法是将其与最简单的基准模型——常数模型进行比较。如果我们对输入一无所知，最好的方法就是预测输出始终是y的平均值。然后，我们可以将这个模型的均方误差与零模型的均方误差进行比较。这个思想在决定系数中得到了形式化定义，其计算公式如下：
在这个公式中，y[i]表示索引为i的元素的值，而是回归模型为同一元素提供的估计值。最后，
是y的均值，代表着零模型，即始终返回相同值的模型。这个公式大致等同于首先计算均方误差与输出方差的比率，最后计算 1 减去这个比率。这样，完美的模型得分为 1，而零模型得分为 0。请注意，可能会得到负分数，这意味着模型非常糟糕，以至于使用均值作为预测值更好。
决定系数可以使用sklearn.metrics模块中的r2_score来获得：
>>> from sklearn.metrics import r2_score
>>> r2 = r2_score(y, lr.predict(x))
>>> print("R2 (on training data): {:.2}".format(r2))
R2 (on training data): 0.31
这个指标也叫做 R²得分。如果你使用线性回归并在训练数据上评估误差，那么它确实对应于相关系数 R 的平方。然而，这个指标更为通用，正如我们讨论的那样，可能会返回一个负值。
计算决定系数的另一种方法是使用LinearRegression对象的score方法：
>>> r2 = lr.score(x,y)
多维回归
到目前为止，我们只使用了一个变量进行预测，即每个住宅的房间数。接下来，我们将使用所有数据来拟合模型，使用多维回归。我们现在尝试基于多个输入预测一个单一的输出（平均房价）。
代码看起来和之前非常相似。实际上，它现在更简单了，因为我们可以直接将boston.data的值传递给fit方法：
>>> x = boston.data
>>> y = boston.target
>>> lr.fit(x, y)
使用所有输入变量，均方根误差仅为 4.7，对应的决定系数为 0.74。这比我们之前的结果要好，说明额外的变量确实有帮助。我们不再能像之前那样轻松显示回归线，因为我们有一个 14 维的回归超平面，而不是一条单独的直线。
然而，我们可以绘制预测值与实际值的对比图。代码如下：
>>> p = lr.predict(x)
>>> plt.scatter(p, y)
>>> plt.xlabel('Predicted price')
>>> plt.ylabel('Actual price')
>>> plt.plot([y.min(), y.max()], [[y.min()], [y.max()]])
最后一行绘制了一条对角线，表示完美一致的情况。这有助于可视化。结果如下面的图所示，其中实线表示对角线（如果预测与真实值完全一致，所有点都会落在这条线上）：
回归问题的交叉验证
如果你还记得我们首次介绍分类时，强调了交叉验证在检验预测质量中的重要性。在回归中，这并不总是做的。事实上，到目前为止，我们在本章中只讨论了训练误差。如果你想自信地推断模型的泛化能力，这是一个错误。因为普通最小二乘法是一个非常简单的模型，这通常不是一个非常严重的错误。换句话说，过拟合的程度较轻。然而，我们仍然应该通过经验来测试这一点，这可以通过 scikit-learn 轻松实现。
我们将使用Kfold类来构建一个 5 折交叉验证循环，测试线性回归的泛化能力：
>>> from sklearn.cross_validation import Kfold
>>> kf = KFold(len(x), n_folds=5)
>>> p = np.zeros_like(y)
>>> for train,test in kf:
... lr.fit(x[train], y[train])
... p[test] = lr.predict(x[test])
>>> rmse_cv = np.sqrt(mean_squared_error(p, y))
>>> print('RMSE on 5-fold CV: {:.2}'.format(rmse_cv))
RMSE on 5-fold CV: 5.6
使用交叉验证，我们得到一个更为保守的估计（即，误差更大）：5.6。与分类问题一样，交叉验证的估计是我们如何将模型泛化到未见数据的更好估计。
普通最小二乘法在学习阶段很快，并且返回一个简单的模型，在预测时也非常快。因此，它通常应该是回归问题中你首先尝试的模型。然而，我们现在将看到更先进的方法，并了解为什么有时它们更为优越。
惩罚或正则化回归
本节介绍了惩罚回归，也叫做正则化回归，它是回归模型的一个重要类别。
在普通回归中，返回的拟合结果是训练数据上的最佳拟合。这可能导致过拟合。惩罚意味着我们为参数值的过度自信添加惩罚。因此，我们接受稍差的拟合，以便拥有一个更简单的模型。
另一种思考方式是，将默认设定为输入变量与输出预测之间没有关系。当我们拥有数据时，我们会改变这一观点，但添加惩罚意味着我们需要更多数据来说服我们，这确实是一个强关系。
提示
惩罚回归是关于权衡的
惩罚回归是偏差-方差权衡的另一个例子。在使用惩罚时，训练数据的拟合度会变差，因为我们引入了偏差。另一方面，我们减少了方差，倾向于避免过拟合。因此，整体结果可能更好地推广到未见过的（测试）数据。
L1 和 L2 惩罚
现在我们将详细探讨这些思想。对于不关心某些数学方面的读者，可以直接跳到下一节，了解如何在 scikit-learn 中使用正则化回归。
一般而言，问题在于我们给定了一个训练数据矩阵 X（行是观测值，每列是不同的特征），以及一个输出值向量 y。目标是得到一个权重向量，我们将其称为 b。普通最小二乘回归由以下公式给出：
也就是说，我们找到一个向量 b，使其最小化与目标 y 的平方距离。在这些方程中，我们忽略了设置截距的问题，假设训练数据已经预处理，使得 y 的均值为零。
添加惩罚项或正则化意味着我们不仅考虑训练数据的最佳拟合，还考虑向量的组成。回归中通常使用两种类型的惩罚：L1 惩罚和 L2 惩罚。L1 惩罚意味着我们通过系数的绝对值之和来惩罚回归，而 L2 惩罚则是通过系数平方和来惩罚。
当我们添加 L1 惩罚时，我们不再使用之前的方程，而是优化以下方程：
在这里，我们试图同时使误差变小，同时也使系数的值（绝对值）变小。使用 L2 惩罚意味着我们使用以下公式：
差异相当微妙：我们现在通过系数的平方来惩罚，而不是它们的绝对值。然而，结果的差异却是显著的。
提示
岭回归、套索回归和弹性网回归
这些带惩罚的模型通常都有一些非常有趣的名字。L1 惩罚模型通常被称为Lasso，而 L2 惩罚模型则被称为岭回归。当同时使用这两者时，我们称之为ElasticNet模型。
Lasso 和 Ridge 都会比无惩罚回归产生更小的系数（绝对值较小，忽略符号）。然而，Lasso 还有一个额外的特性，即它会将许多系数设为零！这意味着最终模型甚至不使用它的一些输入特征，模型是稀疏的。这一特性通常是非常受欢迎的，因为模型在单一步骤中既进行特征选择，又进行回归。
你会注意到，每当我们添加一个惩罚时，我们也会添加一个权重 α，它控制着我们希望的惩罚程度。当α接近零时，我们非常接近于无惩罚回归（事实上，如果你将α设置为零，你将仅执行普通最小二乘法 (OLS)），而当α较大时，我们得到的模型与无惩罚模型非常不同。
岭回归模型较为传统，因为 Lasso 用纸笔计算比较困难。然而，随着现代计算机的发展，我们可以像使用 Ridge 一样轻松地使用 Lasso，甚至可以将它们结合起来形成 ElasticNets。ElasticNet 有两个惩罚项，一个用于绝对值，另一个用于平方项，它解决以下方程：
这个公式是前两个公式的组合，包含了两个参数，α[1] 和 α[2]。在本章稍后，我们将讨论如何为这些参数选择合适的值。
在 scikit-learn 中使用 Lasso 或 ElasticNet
让我们调整之前的例子，使用 ElasticNets。在 scikit-learn 中，替换成 ElasticNet 回归器非常简单，和之前使用最小二乘法回归器一样：
>>> from sklearn.linear_model import ElasticNet, Lasso
>>> en = ElasticNet(alpha=0.5)
现在，我们使用的是en，而之前我们使用的是lr。这是唯一需要改变的地方。结果正是我们预期的那样。训练误差增加到了 5.0（之前是 4.6），但交叉验证误差下降到了 5.4（之前是 5.6）。我们在训练数据上牺牲了较大的误差，以获得更好的泛化能力。我们本可以使用相同的代码，尝试通过Lasso类应用 L1 惩罚，或者使用Ridge类应用 L2 惩罚。
可视化 Lasso 路径
使用 scikit-learn，我们可以轻松地可视化当正则化参数（alpha）变化时所发生的情况。我们将再次使用波士顿数据集，但这次我们将使用Lasso回归对象：
>>> las = Lasso(normalize=1)
>>> alphas = np.logspace(-5, 2, 1000)
>>> alphas, coefs, _= las.path(x, y, alphas=alphas)
对于每个 alpha 值，path 方法在 Lasso 对象上返回能够解决该参数值下 Lasso 问题的系数。由于结果随着 alpha 的变化而平滑变化，因此可以非常高效地计算。
可视化这个路径的典型方法是绘制当 alpha 减小时系数的变化值。你可以按如下方式进行：
>>> fig,ax = plt.subplots()
>>> ax.plot(alphas, coefs.T)
>>> # Set log scale
>>> ax.set_xscale('log')
>>> # Make alpha decrease from left to right
>>> ax.set_xlim(alphas.max(), alphas.min())
这将产生如下图所示的结果（我们省略了添加轴标签和标题的简单代码）：
在这个图中，x轴显示了从左到右逐渐减弱的正则化（alpha 逐渐减小）。每条线显示了不同系数随着 alpha 变化的情况。图表显示，当使用非常强的正则化时（左侧，alpha 非常高），最佳解决方案是所有值都为零。随着正则化的减弱，各个系数的值会一个接一个地首先急剧增加，然后稳定下来。到了某个点，它们都会趋于平稳，因为我们可能已经接近未惩罚的解。
P 大于 N 的情境
本节标题是一些内部术语，你现在将会学习这些。自 1990 年代起，首先在生物医学领域，然后在互联网领域，出现了 P 大于 N 的问题。这意味着特征的数量 P 大于样本的数量 N（这些字母是这些概念的常用统计缩写）。这些问题被称为P 大于 N问题。
例如，如果你的输入是一组书面文档，一种简单的方法是将字典中的每个可能单词视为一个特征，并基于这些特征进行回归（稍后我们会亲自处理类似问题）。在英语中，你有超过 2 万个单词（如果进行词干化并只考虑常见单词的话；如果跳过这个预处理步骤，单词数是它的十倍还多）。如果你只有几百或几千个样本，你将会有更多的特征而非样本。
在这种情况下，由于特征的数量大于样本的数量，因此可能会在训练数据上实现完美拟合。这是一个数学事实，与数据本身无关。实际上，你是在解决一个线性方程组，其中方程数量少于变量数量。你可以找到一组回归系数，训练误差为零（实际上，你可以找到多个完美解，无限多个）。
然而，且这是一个重大问题，零训练误差并不意味着你的解决方案会很好地泛化。事实上，它可能泛化得非常差。虽然早期的正则化可能给你一些额外的提升，但现在它是得到有意义结果的绝对必要条件。
基于文本文档的示例
我们现在转向一个来自卡内基梅隆大学诺亚·史密斯教授研究小组的研究示例。这项研究基于挖掘公司向美国证券交易委员会（SEC）提交的所谓 10-K 报告。这项申报是法律要求所有上市公司进行的。该研究的目标是基于这份公开信息预测公司股票未来的波动性。在训练数据中，我们实际上使用的是已经知道结果的历史数据。
有 16,087 个可用示例。这些特征已经为我们预处理过，表示不同的单词，总共 150,360 个。因此，我们有的特征比示例多得多，几乎是它们的十倍。在引言中提到，普通最小二乘法在这些情况下失败，我们现在通过盲目应用它来看到原因。
数据集可以通过多个来源获得 SVMLight 格式，包括本书的配套网站。这是一个 scikit-learn 可以读取的格式。正如名字所示，SVMLight 是一个支持向量机实现，也可以通过 scikit-learn 使用；目前，我们只关心文件格式：
>>> from sklearn.datasets import load_svmlight_file
>>> data,target = load_svmlight_file('E2006.train')
在前面的代码中，数据是一个稀疏矩阵（即大多数条目为零，因此只保存非零条目在内存中），而目标是一个简单的一维向量。我们可以首先查看目标的一些属性：
>>> print('Min target value: {}'.format(target.min()))
Min target value: -7.89957807347
>>> print('Max target value: {}'.format(target.max()))
Max target value: -0.51940952694
>>> print('Mean target value: {}'.format(target.mean()))
Mean target value: -3.51405313669
>>> print('Std. dev. target: {}'.format(target.std()))
Std. dev. target: 0.632278353911
所以，我们可以看到数据位于-7.9 和-0.5 之间。现在我们对数据有了大致的了解，我们可以检查使用 OLS 进行预测时会发生什么。请注意，我们可以使用与之前在波士顿示例中完全相同的类和方法：
>>> from sklearn.linear_model import LinearRegression
>>> lr = LinearRegression()
>>> lr.fit(data,target)
>>> pred = lr.predict(data)
>>> rmse_train = np.sqrt(mean_squared_error(target, pred))
>>> print('RMSE on training: {:.2}'.format(rmse_train))
RMSE on training: 0.0025
>>> print('R2 on training: {:.2}'.format(r2_score(target, pred)))
R2 on training: 1.0
均方根误差由于四舍五入误差不是完全为零，但它非常接近。决定系数为1.0。也就是说，线性模型在其训练数据上报告了完美的预测结果。
当我们使用交叉验证时（代码与我们之前在波士顿示例中使用的非常相似），我们得到的结果非常不同：RMSE 为 0.75，决定系数为负值-0.42。这意味着如果我们总是“预测”均值-3.5，我们比使用回归模型效果更好！
提示
训练和泛化误差
当特征的数量大于示例的数量时，使用 OLS 总是会得到零训练误差，除非是由于四舍五入问题。然而，这很少是模型在泛化方面表现良好的标志。事实上，你可能会得到零训练误差，却拥有一个完全无用的模型。
自然的解决方案是使用正则化来抵消过拟合。我们可以尝试使用 ElasticNet 学习器的相同交叉验证循环，并将惩罚参数设置为0.1：
>>> from sklearn.linear_model import ElasticNet
>>> met = ElasticNet(alpha=0.1)
>>> kf = KFold(len(target), n_folds=5)
>>> pred = np.zeros_like(target)
>>> for train, test in kf:
... met.fit(data[train], target[train])
... pred[test] = met.predict(data[test])
>>> # Compute RMSE
>>> rmse = np.sqrt(mean_squared_error(target, pred))
>>> print('[EN 0.1] RMSE on testing (5 fold): {:.2}'.format(rmse))
[EN 0.1] RMSE on testing (5 fold): 0.4
>>> # Compute Coefficient of determination
>>> r2 = r2_score(target, pred)
>>> print('[EN 0.1] R2 on testing (5 fold): {:.2}'.format(r2))
[EN 0.1] R2 on testing (5 fold): 0.61
现在，我们得到0.4的 RMSE 和0.61的 R2，远比仅预测均值要好。不过，这个解决方案有一个问题，那就是 alpha 的选择。当使用默认值（1.0）时，结果非常不同（而且更差）。
在这种情况下，我们作弊了，因为作者之前尝试了几个值，看看哪些值能给出好的结果。这是无效的，并且可能导致对信心的高估（我们正在查看测试数据来决定使用哪些参数值，而哪些应该永远不使用）。下一节将解释如何正确地做这件事，以及 scikit-learn 如何支持这一点。
以原则化的方式设置超参数
在前面的例子中，我们将惩罚参数设置为0.1。我们同样也可以将其设置为 0.7 或 23.9。显然，结果会因每次而异。如果我们选择一个过大的值，我们会出现欠拟合。在极端情况下，学习系统将返回所有系数为零的结果。如果我们选择一个过小的值，我们就非常接近普通最小二乘法（OLS），这会导致过拟合且泛化能力差（正如我们之前所看到的）。
我们如何选择一个合适的值呢？这是机器学习中的一个普遍问题：为我们的学习方法设置参数。一个通用的解决方案是使用交叉验证。我们选择一组可能的值，然后使用交叉验证来选择最优值。这需要更多的计算（如果我们使用五个子集，计算量是五倍），但始终适用且没有偏差。
然而，我们必须小心。为了获得泛化的估计，我们需要使用两级交叉验证：一级是估计泛化能力，二级是获得好的参数。也就是说，我们将数据分成例如五个子集。我们首先保留第一个子集，并在其他四个子集上进行学习。接下来，我们将这四个子集再次分成五个子集，用以选择参数。一旦设置好参数，我们在第一个子集上进行测试。然后，我们重复这一过程另外四次：
上图展示了如何将单一的训练子集分割成子子集。我们需要对所有其他子集重复这一过程。在这个例子中，我们有五个外部子集和五个内部子集，但并不一定非要使用相同数量的外部和内部子集，只要确保每个子集之间是分开的，你可以选择任何数量。
这会导致大量的计算，但为了正确地执行，确实是必要的。问题在于，如果你使用一部分数据来对模型做出任何决策（包括选择设置哪些参数），那么你就已经污染了这部分数据，无法再用它来测试模型的泛化能力。这是一个微妙的点，可能并不容易立刻察觉。事实上，许多机器学习用户仍然犯这个错误，过高估计他们系统的表现，因为他们没有正确地执行交叉验证！
幸运的是，scikit-learn 使得做对的事情变得非常容易；它提供了名为LassoCV、RidgeCV和ElasticNetCV的类，这些类都封装了一个内部交叉验证循环，用于优化必要的参数。代码几乎和之前的一样，唯一不同的是我们不需要为 alpha 指定任何值：
>>> from sklearn.linear_model import ElasticNetCV
>>> met = ElasticNetCV()
>>> kf = KFold(len(target), n_folds=5)
>>> p = np.zeros_like(target)
>>> for train,test in kf:
... met.fit(data[train],target[train])
... p[test] = met.predict(data[test])
>>> r2_cv = r2_score(target, p)
>>> print("R2 ElasticNetCV: {:.2}".format(r2_cv))
R2 ElasticNetCV: 0.65
这会导致大量计算，因此在等待时你可能想要喝杯咖啡（取决于你的计算机速度有多快）。通过利用多个处理器，你可能会获得更好的性能。这是 scikit-learn 的内置功能，可以通过将n_jobs参数传递给ElasticNetCV构造函数来轻松访问。要使用四个 CPU，可以使用以下代码：
>>> met = ElasticNetCV(n_jobs=4)
将n_jobs参数设置为-1以使用所有可用的 CPU：
>>> met = ElasticNetCV(n_jobs=-1)
也许你会想知道，如果 ElasticNets 有两个惩罚项，L1 和 L2 惩罚项，为什么我们只需要设置一个 alpha 值。实际上，通过分别指定 alpha 和 l1_ratio 变量（拼写为 ell-1-underscore-ratio），可以指定这两个值。然后，α1 和 α2 设置如下（其中 ρ 代表 l1_ratio）：
直观地说，alpha 设置了总体正则化的量，而l1_ratio设置了不同类型正则化（L1 和 L2）之间的权衡。
我们可以要求ElasticNetCV对象测试不同的l1_ratio值，如下所示的代码：
>>> l1_ratio=[.01, .05, .25, .5, .75, .95, .99]
>>> met = ElasticNetCV(
l1_ratio=l1_ratio,
n_jobs=-1)
这组l1_ratio值在文档中被推荐使用。它将测试几乎类似于 Ridge（当l1_ratio为 0.01 或 0.05 时）以及几乎类似于 Lasso（当l1_ratio为 0.95 或 0.99 时）的模型。因此，我们探索了各种不同选项的完整范围。
由于其灵活性和能够使用多个 CPU 的能力，当你没有任何特定原因偏好一种模型而不是其他模型时，ElasticNetCV是回归问题的一个出色的默认解决方案。
将所有这些放在一起，我们现在可以在这个大数据集上可视化预测与真实拟合：
>>> l1_ratio = [.01, .05, .25, .5, .75, .95, .99]
>>> met = ElasticNetCV(
l1_ratio=l1_ratio,
n_jobs=-1)
>>> p = np.zeros_like(target)
>>> for train,test in kf:
... met.fit(data[train],target[train])
... p[test] = met.predict(data[test])
>>> plt.scatter(p, y)
>>> # Add diagonal line for reference
>>> # (represents perfect agreement)
>>> plt.plot([p.min(), p.max()], [p.min(), p.max()])
这导致以下图表：
我们可以看到，在值范围的底端，预测结果并不很好匹配。这可能是因为在目标范围的这一端元素较少（这也意味着这只影响了少数数据点）。
最后一点：使用内部交叉验证循环来设置参数的方法在 scikit-learn 中也是可用的，可以使用网格搜索。实际上，我们在上一章中已经使用过了。
摘要
在本章中，我们从最古老的技巧开始，普通最小二乘回归。尽管有几个世纪的历史，但它仍然经常是回归问题的最佳解决方案。然而，我们也看到了更现代的方法，避免过拟合，并且在具有大量特征时可以给出更好的结果。我们使用了 Ridge、Lasso 和 ElasticNets；这些是回归问题的最先进方法。
我们再次看到了依赖训练误差来估计泛化能力的危险：它可能会给出过于乐观的估计，甚至使我们的模型在训练误差上为零，但我们知道它完全没有用。在思考这些问题时，我们引入了二级交叉验证，这是一个重要的概念，许多领域的从业者仍未完全理解。
在本章中，我们能够依赖 scikit-learn 来支持我们想要执行的所有操作，包括一种简便的方法来实现正确的交叉验证。带有内部交叉验证循环的 ElasticNet（用于参数优化，scikit-learn 中由ElasticNetCV实现）可能应该成为你回归分析的默认方法。
使用替代方法的一个原因是当你对稀疏解感兴趣时。在这种情况下，纯 Lasso 解更为合适，因为它会将许多系数设为零。它还会让你从数据中发现少数几个对输出至关重要的变量。了解这些变量的身份本身可能就很有趣，除了获得一个优秀的回归模型之外。
在下一章，我们将讨论推荐系统，这是另一个机器学习问题。我们首先的方法是使用回归来预测消费者产品评分。然后，我们将看到生成推荐的替代模型。
第八章：推荐系统
推荐系统已经成为在线服务和电商的基础之一。这种自动化系统可以为每个用户提供个性化的建议列表（无论是购买的产品、使用的功能还是新的社交连接）。在本章中，我们将看到自动化推荐生成系统的基本工作原理。基于消费者输入的推荐领域通常被称为协同过滤，因为用户通过系统进行协作，帮助彼此找到最佳产品。
在本章的第一部分，我们将看到如何利用消费者过去的产品评分来预测新的评分。我们从一些有用的想法开始，然后将它们结合在一起。在结合时，我们使用回归分析来学习它们可以如何最优地组合。这也将让我们探讨机器学习中的一个通用概念：集成学习。
在本章的第二部分，我们将探讨一种不同的推荐学习方法：购物篮分析。与我们拥有数字评分的情况不同，在购物篮分析中，我们仅拥有关于购物篮的信息，也就是说，哪些商品是一起购买的。目标是学习如何进行推荐。你可能已经在在线购物中见过类似“购买 X 的人也购买了 Y”这样的功能。我们将开发出一个类似的功能。
评分预测和推荐
如果你在过去 10 年里使用过任何在线购物系统，你可能已经见过这些推荐。有些类似于亚马逊的“购买 X 的客户也购买了 Y”。这些将在本章的购物篮分析部分中进一步探讨。其他推荐则基于预测产品的评分，比如电影的评分。
基于过去产品评分学习推荐的问题由 Netflix 大奖而闻名，Netflix 大奖是 Netflix 举办的百万美元机器学习公开挑战赛。Netflix（在美国和英国非常知名，并正在进行国际扩展）是一家电影租赁公司。传统上，你会收到邮寄来的 DVD；最近，Netflix 专注于在线电影和电视节目的流媒体播放。从一开始，Netflix 的一个独特之处在于它允许用户对看过的电影进行评分。Netflix 随后使用这些评分向用户推荐其他电影。在这个机器学习问题中，你不仅知道用户看了哪些电影，还知道用户如何评分。
2006 年，Netflix 将其数据库中大量客户对电影的评分数据公开，举行了一场公开挑战。目标是改进 Netflix 内部的评分预测算法。任何能够将其提升 10％或更多的人将赢得 100 万美元的奖金。2009 年，一个名为 BellKor's Pragmatic Chaos 的国际团队成功突破了这个标准，并获得了奖金。他们是在另一个团队 The Ensemble 的前 20 分钟成功做到的，并且同样突破了 10％的改进——这场持续了数年的竞赛最终以一个激动人心的结局落下帷幕。
提示
现实世界中的机器学习
关于 Netflix 奖已经有很多相关的讨论，你可以通过阅读相关资料了解更多。获奖的技术是先进的机器学习方法与大量数据预处理工作相结合的结果。例如，有些用户倾向于给所有电影打很高的分，而有些用户总是给出较低的评价；如果你在预处理阶段不考虑这一点，你的模型就会受到影响。为了获得好的结果，还需要进行其他的归一化处理：例如，电影的上映年份以及它收到的评价数量。好的算法很重要，但你始终需要“亲自动手”，根据你手头数据的特性调整你的方法。数据的预处理和归一化往往是机器学习过程中最耗时的部分。然而，这也是对系统最终表现产生最大影响的地方。
关于 Netflix 奖，首先要注意的是它的难度。大致来说，Netflix 使用的内部系统比没有推荐的系统（即每个电影仅给所有用户的平均值）好 10%左右。目标只是要在此基础上再提升 10%。最终，获奖系统比没有个性化的系统好大约 20%。然而，要实现这个目标，付出了巨大的时间和努力。尽管 20%的提升看起来并不多，但最终的结果是一个在实际中非常有用的系统。
不幸的是，出于法律原因，这个数据集目前已经无法获取。尽管数据是匿名的，但人们担心可能会揭示出客户的身份以及电影租赁的私人信息。不过，我们可以使用一个具有类似特征的学术数据集。这些数据来自 GroupLens，这是明尼苏达大学的一个研究实验室。
我们如何解决类似 Netflix 的评分预测问题呢？我们将看到两种不同的方法，邻域方法和回归方法。我们还会看到如何将这两种方法结合起来，得出一个统一的预测结果。
划分训练集和测试集
从高层次上看，将数据集划分为训练数据和测试数据，以便获得系统性能的原则性估计，方法与之前的章节相同：我们将取一定比例的数据点（我们将使用 10%）并将其保留用于测试；其余数据用于训练。然而，由于在此上下文中数据的结构不同，因此代码也有所不同。第一步是从磁盘加载数据，我们使用以下函数：
def load():
import numpy as np
from scipy import sparse
data = np.loadtxt('data/ml-100k/u.data')
ij = data[:, :2]
ij -= 1 # original data is in 1-based system
values = data[:, 2]
reviews = sparse.csc_matrix((values, ij.T)).astype(float)
return reviews.toarray()
请注意，这个矩阵中的零条目表示缺失的评分。
>>> reviews = load()
>>> U,M = np.where(reviews)
现在，我们使用标准的随机模块选择要测试的索引：
>>> import random
>>> test_idxs = np.array(random.sample(range(len(U)), len(U)//10))
现在，我们构建 train 矩阵，它类似于 reviews，但将测试条目的值设为零：
>>> train = reviews.copy()
>>> train[U[test_idxs], M[test_idxs]] = 0
最后，test 矩阵只包含测试值：
>>> test = np.zeros_like(reviews)
>>> test[U[test_idxs], M[test_idxs]] = reviews[U[test_idxs], M[test_idxs]]
从现在起，我们将开始处理训练数据，并尝试预测数据集中所有缺失的条目。也就是说，我们将编写代码，为每个（用户，电影）对分配一个推荐。
对训练数据进行标准化
正如我们之前讨论的那样，最好的做法是对数据进行标准化，以去除明显的电影或用户特定的效应。我们将使用一种非常简单的标准化方法，之前也使用过：转换为 z 分数。
不幸的是，我们不能简单地使用 scikit-learn 的标准化对象，因为我们必须处理数据中的缺失值（即，并非所有电影都被所有用户评分）。因此，我们希望通过实际存在的值的均值和标准差进行标准化。
我们将编写自己的类，它忽略缺失值。这个类将遵循 scikit-learn 的预处理 API：
class NormalizePositive(object):
我们希望选择标准化的轴。默认情况下，我们沿第一个轴进行标准化，但有时沿第二个轴进行标准化也会很有用。这遵循了许多其他与 NumPy 相关的函数的约定：
def __init__(self, axis=0):
self.axis = axis
最重要的方法是 fit 方法。在我们的实现中，我们计算的是非零值的均值和标准差。请记住，零值表示“缺失值”：
def fit(self, features, y=None):
如果轴为 1，我们将在转置数组上进行操作，如下所示：
if self.axis == 1:
features = features.T
# count features that are greater than zero in axis 0:
binary = (features > 0)
count0 = binary.sum(axis=0)
# to avoid division by zero, set zero counts to one:
count0[count0 == 0] = 1.
# computing the mean is easy:
self.mean = features.sum(axis=0)/count0
# only consider differences where binary is True:
diff = (features - self.mean) * binary
diff **= 2
# regularize the estimate of std by adding 0.1
self.std = np.sqrt(0.1 + diff.sum(axis=0)/count0)
return self
我们将 0.1 加到标准差的直接估计中，以避免在样本很少且所有样本可能完全相同的情况下低估标准差的值。所用的确切值对最终结果影响不大，但我们需要避免除以零的情况。
transform 方法需要维护二进制结构，如下所示：
def transform(self, features):
if self.axis == 1:
features = features.T
binary = (features > 0)
features = features - self.mean
features /= self.std
features *= binary
if self.axis == 1:
features = features.T
return features
请注意，当轴为 1 时，我们如何处理输入矩阵的转置，然后再将其转换回来，以便返回值与输入矩阵的形状相同。inverse_transform 方法执行逆操作以进行如这里所示的转换：
def inverse_transform(self, features, copy=True):
if copy:
features = features.copy()
if self.axis == 1:
features = features.T
features *= self.std
features += self.mean
if self.axis == 1:
features = features.T
return features
最后，我们添加了 fit_transform 方法，顾名思义，它结合了 fit 和 transform 两个操作：
def fit_transform(self, features):
return self.fit(features).transform(features)
我们定义的方法（fit、transform、transform_inverse 和 fit_transform）与 sklearn.preprocessing 模块中定义的对象相同。在接下来的章节中，我们将首先规范化输入，生成规范化的预测值，最后应用逆变换以获得最终预测结果。
一种基于邻域的推荐方法
邻域概念可以通过两种方式实现：用户邻域或电影邻域。用户邻域基于一个非常简单的概念：要知道一个用户如何评分电影，找到与他们最相似的用户，并查看他们的评分。我们暂时只考虑用户邻域。在本节结束时，我们将讨论如何调整代码来计算电影邻域。
我们现在要探讨的一种有趣技巧是，仅仅查看每个用户评分过的电影，即使不查看他们给出的评分。即使我们有一个二元矩阵，其中用户评分的电影用 1 表示，未评分的电影用 0 表示，我们仍然可以做出有用的预测。事后看来，这完全是合理的；我们并不是随机选择电影观看，而是选择那些我们已经有一定喜好预期的电影。我们也不是随机选择要评分的电影，而是可能只会评分那些我们感受最强烈的电影（当然，也有例外，但平均来看，这大概是对的）。
我们可以将矩阵的值可视化为一张图像，每个评分表示为一个小方块。黑色表示没有评分，灰度级表示评分值。
可视化数据的代码非常简单（你可以调整它来显示矩阵的更多部分，而不仅仅是本书中能显示的部分），如下所示：
>>> from matplotlib import pyplot as plt
>>> # Build an instance of the object we defined above
>>> norm = NormalizePositive(axis=1)
>>> binary = (train > 0)
>>> train = norm.fit_transform(train)
>>> # plot just 200x200 area for space reasons
>>> plt.imshow(binary[:200, :200], interpolation='nearest')
以下截图是该代码的输出：
我们可以看到矩阵是稀疏的——大部分格子是黑色的。我们也可以看到，有些用户评分的电影比其他用户多，而有些电影比其他电影有更多的评分。
我们现在将使用这个二元矩阵来预测电影评分。一般的算法将是（伪代码）如下：
-
对于每个用户，按相似度对其他所有用户进行排序。此步骤中，我们将使用二元矩阵，并使用相关性作为相似度的度量（将二元矩阵视为零和一使得我们可以进行这种计算）。
-
当我们需要估计一个（用户，电影）对的评分时，我们会查看所有评分过该电影的用户，并将他们分为两组：最相似的那一半和最不相似的那一半。然后我们使用最相似那一半的平均值作为预测。
我们可以使用scipy.spatial.distance.pdist函数来获取所有用户之间的距离矩阵。这个函数返回相关距离，它通过反转相关值来转化，使得更大的数值表示不那么相似。从数学上讲，相关距离是，其中
是相关值。代码如下：
>>> from scipy.spatial import distance
>>> # compute all pair-wise distances:
>>> dists = distance.pdist(binary, 'correlation')
>>> # Convert to square form, so that dists[i,j]
>>> # is distance between binary[i] and binary[j]:
>>> dists = distance.squareform(dists)
我们可以使用这个矩阵来选择每个用户的最近邻居。这些用户是最相似的用户。
>>> neighbors = dists.argsort(axis=1)
现在，我们遍历所有用户以估算所有输入的预测值：
>>> # We are going to fill this matrix with results
>>> filled = train.copy()
>>> for u in range(filled.shape[0]):
... # n_u is neighbors of user
... n_u = neighbors[u, 1:]
... # t_u is training data
... for m in range(filled.shape[1]):
... # get relevant reviews in order!
... revs = [train[neigh, m]
... for neigh in n_u
... if binary [neigh, m]]
... if len(revs):
... # n is the number of reviews for this movie
... n = len(revs)
... # consider half of the reviews plus one
... n //= 2
... n += 1
... revs = revs[:n]
... filled[u,m] = np.mean(revs )
上述代码片段中的关键部分是通过正确的值进行索引，以选择已经对电影评分的邻居。然后，我们选择距离用户最近的那一半（在rev[:n]这一行中）并对它们求平均。由于一些电影有很多评论，而其他电影的评论非常少，因此很难为所有情况找到一个统一的用户数量。选择可用数据的那一半是更通用的方法。
要得到最终结果，我们需要将预测值反归一化，具体如下：
>>> predicted = norm.inverse_transform(filled)
我们可以使用上一章学到的相同度量方法：
>>> from sklearn import metrics
>>> r2 = metrics.r2_score(test[test > 0], predicted[test > 0])
>>> print('R2 score (binary neighbors): {:.1%}'.format(r2))
R2 score (binary neighbors): 29.5%
上述代码计算了用户邻居的结果，但我们可以通过简单地转置输入矩阵来计算电影邻居。事实上，代码计算的是输入矩阵的行的邻居。
所以我们可以重新运行以下代码，只需在顶部插入以下行：
>>> reviews = reviews.T
>>> # use same code as before …
>>> r2 = metrics.r2_score(test[test > 0], predicted[test > 0])
>>> print('R2 score (binary movie neighbors): {:.1%}'.format(r2))
R2 score (binary movie neighbors): 29.8%
因此我们可以看到，结果并没有那么不同。
在本书的代码库中，邻居代码已被封装成一个简单的函数，便于重用。
一种回归方法推荐系统
邻居的替代方案是将推荐问题表述为回归问题，并应用我们在上一章学到的方法。
我们还考虑为什么这个问题不适合分类模型。我们当然可以尝试学习一个五类模型，每个类别对应一种可能的电影评分。这个方法有两个问题：
-
不同的错误类型是完全不同的。例如，把一部 5 星电影误认为 4 星电影并不是那么严重的错误，而把 5 星电影误认为 1 星电影就严重得多。
-
中间值是有意义的。即使我们的输入只有整数值，说预测值是 4.3 也是完全有意义的。我们可以看到，这与 3.5 的预测是不同的，尽管它们都四舍五入到 4。
这两个因素共同表明，分类方法不适合这个问题。回归框架更加适合。
对于基本方法，我们有两种选择：我们可以构建电影特定的模型或用户特定的模型。在我们的例子中，我们将首先构建用户特定的模型。这意味着，对于每个用户，我们将用户已评分的电影作为目标变量。输入是其他用户的评分。我们假设，这样做将为与我们用户相似的用户赋予较高的值（或者为喜欢我们用户不喜欢的电影的用户赋予负值）。
设置train和test矩阵与之前一样（包括执行标准化步骤）。因此，我们直接进入学习步骤。首先，我们按如下方式实例化一个回归器：
>>> reg = ElasticNetCV(alphas=[
0.0125, 0.025, 0.05, .125, .25, .5, 1., 2., 4.])
我们构建一个数据矩阵，矩阵中包含每一对（用户，电影）的评分。我们将其初始化为训练数据的副本：
>>> filled = train.copy()
现在，我们遍历所有用户，每次仅基于该用户提供的数据学习回归模型：
>>> for u in range(train.shape[0]):
... curtrain = np.delete(train, u, axis=0)
... # binary records whether this rating is present
... bu = binary[u]
... # fit the current user based on everybody else
... reg.fit(curtrain[:,bu].T, train[u, bu])
... # Fill in all the missing ratings
... filled[u, ~bu] = reg.predict(curtrain[:,~bu].T)
评估该方法可以像以前一样进行：
>>> predicted = norm.inverse_transform(filled)
>>> r2 = metrics.r2_score(test[test > 0], predicted[test > 0])
>>> print('R2 score (user regression): {:.1%}'.format(r2))
R2 score (user regression): 32.3%
与之前一样，我们可以调整这段代码，使用转置矩阵执行电影回归。
结合多种方法
我们现在将上述方法结合在一起，进行单一预测。直观上看，这个想法是好的，但我们如何在实践中做到这一点呢？也许，第一个想到的办法是我们可以对预测结果进行平均。这可能会给出不错的结果，但并没有理由认为所有的预测结果应该被同等对待。可能某些预测比其他预测更好。
我们可以尝试加权平均，在求和之前将每个预测值乘以一个给定的权重。然而，我们如何找到最佳的权重呢？当然，我们是通过从数据中学习它们！
提示
集成学习
我们正在使用一种机器学习中的通用技术，这种技术不仅适用于回归问题：集成学习。我们学习一个预测器的集成（即一组预测器）。然后，我们将它们组合起来得到一个单一的输出。有趣的是，我们可以将每一个预测视为一个新的特征，我们现在只是根据训练数据组合特征，这正是我们一直在做的事情。请注意，这里我们是在做回归，但相同的思路也适用于分类问题：你学习多个分类器，然后使用一个主分类器，该分类器接受所有分类器的输出并给出最终预测。不同形式的集成学习在于如何组合基础预测器。
为了结合这些方法，我们将使用一种叫做堆叠学习的技术。其思路是，你先学习一组预测器，然后将这些预测器的输出作为另一个预测器的特征。你甚至可以有多层，每一层都通过使用上一层的输出作为其预测的特征来进行学习。看看下面的图示：
为了拟合这个组合模型，我们将训练数据分成两部分。或者，我们可以使用交叉验证（原始的堆叠学习模型就是这样工作的）。然而，在这种情况下，我们有足够的数据，通过留出一部分数据来获得良好的估计。
就像在调整超参数时一样，我们需要两个层次的训练/测试划分：一个更高层次的划分，然后在训练划分内部，进行第二次划分，以便拟合堆叠学习器，如下所示：
>>> train,test = load_ml100k.get_train_test(random_state=12)
>>> # Now split the training again into two subgroups
>>> tr_train,tr_test = load_ml100k.get_train_test(train, random_state=34)
>>> # Call all the methods we previously defined:
>>> # these have been implemented as functions:
>>> tr_predicted0 = regression.predict(tr_train)
>>> tr_predicted1 = regression.predict(tr_train.T).T
>>> tr_predicted2 = corrneighbours.predict(tr_train)
>>> tr_predicted3 = corrneighbours.predict(tr_train.T).T
>>> tr_predicted4 = norm.predict(tr_train)
>>> tr_predicted5 = norm.predict(tr_train.T).T
>>> # Now assemble these predictions into a single array:
>>> stack_tr = np.array([
... tr_predicted0[tr_test > 0],
... tr_predicted1[tr_test > 0],
... tr_predicted2[tr_test > 0],
... tr_predicted3[tr_test > 0],
... tr_predicted4[tr_test > 0],
... tr_predicted5[tr_test > 0],
... ]).T
>>> # Fit a simple linear regression
>>> lr = linear_model.LinearRegression()
>>> lr.fit(stack_tr, tr_test[tr_test > 0])
现在，我们将整个过程应用于测试集并进行评估：
>>> stack_te = np.array([
... tr_predicted0.ravel(),
... tr_predicted1.ravel(),
... tr_predicted2.ravel(),
... tr_predicted3.ravel(),
... tr_predicted4.ravel(),
... tr_predicted5.ravel(),
... ]).T
>>> predicted = lr.predict(stack_te).reshape(train.shape)
评估与之前相同：
>>> r2 = metrics.r2_score(test[test > 0], predicted[test > 0])
>>> print('R2 stacked: {:.2%}'.format(r2))
R2 stacked: 33.15%
堆叠学习的结果比任何单一方法的结果都要好。将方法结合起来通常是一种简单的方式，可以获得小幅的性能提升，但结果并不会引起轰动。
通过灵活地组合多种方法，我们可以简单地通过将任何想法添加到学习器的混合中，并让系统将其折入预测中来尝试任何我们希望的想法。例如，我们可以替换最近邻代码中的邻域标准。
然而，我们必须小心不要让数据集过拟合。事实上，如果我们随意尝试太多东西，其中一些可能在这个数据集上效果很好，但无法推广到其他数据。尽管我们正在划分数据集，但我们并没有严格地进行交叉验证我们的设计决策。为了获得良好的估计，并且如果数据量充足，你应该将一部分数据留着，直到你有一个即将投入生产的最终模型。然后，在这个保留的数据上测试你的模型，可以给你一个无偏的预测，了解它在现实世界中的表现如何。
购物篮分析
到目前为止，我们讨论的方法在你拥有用户对产品喜好程度的数字评分时效果很好。但这种信息并不总是可用，因为它需要消费者的主动行为。
购物篮分析是一种推荐学习的替代模式。在这种模式下，我们的数据仅包括哪些商品是一起购买的；它不包含任何关于单个商品是否被喜好的信息。即使用户有时购买了他们后悔的商品，平均而言，知道他们的购买记录也足以帮助你构建良好的推荐系统。获取这类数据通常比评分数据更容易，因为许多用户不会提供评分，而购物篮数据则是购物的副作用。以下截图展示了亚马逊网站上托尔斯泰经典小说《战争与和平》的页面片段，演示了如何使用这些结果的常见方式：
这种学习模式不仅仅适用于实际的购物篮，自然也可以应用于任何有一组对象并且需要推荐其他对象的场景。例如，Gmail 会向正在写电子邮件的用户推荐额外的收件人，类似的技术可以应用于这种推荐（我们并不知道 Gmail 内部使用了什么技术；也许，他们像我们之前所做的那样，结合了多种技术）。或者，我们也可以使用这些方法开发一个应用程序，根据用户的浏览历史推荐访问的网页。即使我们处理的是购物，也有可能将顾客的所有购买合并成一个购物篮，而不考虑这些商品是否是一起购买的或是在不同交易中购买的。这取决于商业环境，但请记住，这些技术是灵活的，可以在许多不同的场合中发挥作用。
注意
啤酒和尿布。关于购物篮分析，常常提到的一个故事是尿布和啤酒的故事。这个故事说的是，当超市开始分析他们的数据时，他们发现尿布常常和啤酒一起购买。可以假设是父亲去超市买尿布时顺便也买了啤酒。关于这个故事是否真实，还是仅仅一个都市传说，存在很多讨论。在这种情况下，似乎是事实。1990 年代初，Osco Drug 确实发现，在傍晚时分，啤酒和尿布经常一起购买，这让当时的经理们感到非常惊讶，因为他们之前从未认为这两种产品有任何相似性。事实并非如此的是，这并没有导致商店将啤酒陈列架移到尿布区域附近。而且，我们也不知道是否真的是父亲比母亲（或祖父母）更多地购买啤酒和尿布。
获取有用的预测
这不仅仅是“购买了 X 的顾客也购买了 Y”，尽管许多在线零售商是这样表述的（参见之前给出的亚马逊截图）；一个真正的系统不能这样运作。为什么不行？因为这样的系统会被非常频繁购买的商品所欺骗，推荐的只是那些流行的商品，而没有任何个性化的推荐。
例如，在超市中，许多顾客每次购物时都会购买面包，或者购买的时间非常接近（为了说明问题，我们假设有 50%的购物包含面包）。因此，如果你关注任何特定商品，比如洗碗液，并查看与洗碗液一起常常购买的商品，你可能会发现面包常常与洗碗液一起购买。事实上，仅仅通过随机机会，每次有人购买洗碗液时，50%的几率他们也会购买面包。然而，面包与其他商品一起购买是因为每个人购买面包的频率都非常高。
我们真正寻找的是“购买了 X 的顾客，比那些没有购买 X 的普通顾客更有可能购买 Y”。如果你购买了洗碗机洗涤剂，你可能会购买面包，但不会比基准更频繁。类似地，一家书店如果只是推荐畅销书，而不考虑你已经购买的书籍，那就不算是很好的个性化推荐。
分析超市购物篮
作为例子，我们将看一个包含比利时超市匿名交易的数据集。这个数据集由哈瑟尔特大学的 Tom Brijs 提供。由于隐私问题，数据已经被匿名化，因此我们只有每个产品的编号，购物篮是由一组编号组成。该数据文件可以从多个在线来源获得（包括本书的伴随网站）。
我们首先加载数据集并查看一些统计信息（这总是个好主意）：
>>> from collections import defaultdict
>>> from itertools import chain
>>> # File is downloaded as a compressed file
>>> import gzip
>>> # file format is a line per transaction
>>> # of the form '12 34 342 5...'
>>> dataset = [[int(tok) for tok in line.strip().split()]
... for line in gzip.open('retail.dat.gz')]
>>> # It is more convenient to work with sets
>>> dataset = [set(d) for d in dataset]
>>> # count how often each product was purchased:
>>> counts = defaultdict(int)
>>> for elem in chain(*dataset):
... counts[elem] += 1
我们可以在以下表格中看到总结的结果：
| 购买次数 | 产品数量 |
|---|---|
| 仅购买一次 | 2,224 |
| 2 或 3 | 2,438 |
| 4 到 7 | 2,508 |
| 8 到 15 | 2,251 |
| 16 到 31 | 2,182 |
| 32 到 63 | 1,940 |
| 64 到 127 | 1,523 |
| 128 到 511 | 1,225 |
| 512 次或更多 | 179 |
有许多产品只被购买了少数几次。例如，33%的产品购买次数为四次或更少。然而，这仅占所有购买的 1%。这种许多产品仅被购买少数次的现象有时被称为长尾现象，随着互联网使得库存和销售小众商品变得更加便宜，这一现象变得更加显著。为了能够为这些产品提供推荐，我们需要更多的数据。
有一些开源的购物篮分析算法实现，但没有一个与 scikit-learn 或我们一直在使用的其他软件包很好地集成。因此，我们将自己实现一个经典的算法。这个算法叫做 Apriori 算法，虽然它有点老（由 Rakesh Agrawal 和 Ramakrishnan Srikant 于 1994 年发布），但它仍然有效（算法当然永远有效，只是会被更好的想法所取代）。
从形式上讲，Apriori 算法接收一个集合（即你的购物篮），并返回作为子集非常频繁的集合（即共同出现在许多购物篮中的商品）。
该算法采用自下而上的方法：从最小的候选项集（由一个单独的元素组成）开始，逐步构建，每次添加一个元素。正式来说，算法接受一个购物篮集和应考虑的最小输入（我们称之为 minsupport）。第一步是考虑所有仅包含一个元素且具有最小支持度的购物篮。然后，这些项集以所有可能的方式组合，构建出二元素购物篮。接着，筛选出仅保留那些具有最小支持度的项集。然后，考虑所有可能的三元素购物篮，并保留那些具有最小支持度的项集，如此继续。Apriori 的技巧在于，当构建一个更大的购物篮时，它只需要考虑由更小的集合构成的购物篮。
以下图示展示了该算法的示意图：
现在我们将在代码中实现这个算法。我们需要定义我们寻找的最小支持度：
>>> minsupport = 80
支持度是指一组产品一起被购买的次数。Apriori 的目标是找到具有高支持度的项集。从逻辑上讲，任何支持度高于最小支持度的项集只能由那些本身至少具有最小支持度的项组成：
>>> valid = set(k for k,v in counts.items()
... if (v >= minsupport))
我们的初始项集是单一项集（包含单个元素的集合）。特别地，所有至少具有最小支持度的单一项集都是频繁项集：
>>> itemsets = [frozenset([v]) for v in valid]
现在，我们的循环如下所示：
>>> freqsets = []
>>> for i in range(16):
... nextsets = []
... tested = set()
... for it in itemsets:
... for v in valid:
... if v not in it:
... # Create a new candidate set by adding v to it
... c = (it | frozenset([v]))
... # check If we have tested it already
... if c in tested:
... continue
... tested.add(c)
...
... # Count support by looping over dataset
... # This step is slow.
... # Check `apriori.py` for a better implementation.
... support_c = sum(1 for d in dataset if d.issuperset(c))
... if support_c > minsupport:
... nextsets.append(c)
... freqsets.extend(nextsets)
... itemsets = nextsets
... if not len(itemsets):
... break
>>> print("Finished!")
Finished!
这个方法是正确的，但比较慢。一个更好的实现有更多的基础设施来避免遍历所有数据集来获取计数（support_c）。特别地，我们跟踪哪些购物篮包含哪些频繁项集。这样可以加速循环，但也使代码更难理解。因此，我们在这里不展示它。像往常一样，您可以在本书的配套网站上找到这两种实现。网站上的代码也被封装成一个函数，可以应用于其他数据集。
Apriori 算法返回频繁项集，即那些出现频率超过某一阈值（由代码中的minsupport变量给定）的购物篮。
关联规则挖掘
频繁项集本身并不是很有用。下一步是构建关联规则。由于这个最终目标，整个购物篮分析领域有时被称为关联规则挖掘。
关联规则是一种“如果 X，则 Y”的语句，例如，“如果顾客购买了《战争与和平》，那么他们将购买《安娜·卡列尼娜》”。请注意，这条规则并不是确定性的（并不是所有购买 X 的顾客都会购买 Y），但每次都要这样表达是很繁琐的：“如果顾客购买了 X，他购买 Y 的可能性比基线高”；因此，我们说“如果 X，则 Y”，但我们是从概率的角度来理解这句话。
有趣的是，前件和结论可能都包含多个对象：购买了 X、Y 和 Z 的顾客也购买了 A、B 和 C。多个前件可能使你能够做出比单个项目更具体的预测。
你可以通过尝试所有可能的 X 蕴含 Y 的组合，从频繁项集生成规则。生成这些规则很容易。然而，你只希望得到有价值的规则。因此，我们需要衡量规则的价值。一个常用的度量叫做提升度。提升度是应用规则所得到的概率与基准概率之间的比值，公式如下：
在上述公式中，P(Y)表示包含 Y 的所有交易所占的比例，而 P(Y|X)表示在交易包含 X 的前提下，包含 Y 的交易所占的比例。使用提升度可以避免推荐畅销书的问题；对于畅销书，P(Y)和 P(Y|X)都将较大。因此，提升度会接近 1，规则将被视为无关紧要。实际上，我们希望提升度的值至少为 10，甚至可能达到 100。
请参考以下代码：
>>> minlift = 5.0
>>> nr_transactions = float(len(dataset))
>>> for itemset in freqsets:
... for item in itemset:
... consequent = frozenset([item])
... antecedent = itemset-consequent
... base = 0.0
... # acount: antecedent count
... acount = 0.0
...
... # ccount : consequent count
... ccount = 0.0
... for d in dataset:
... if item in d: base += 1
... if d.issuperset(itemset): ccount += 1
... if d.issuperset(antecedent): acount += 1
... base /= nr_transactions
... p_y_given_x = ccount/acount
... lift = p_y_given_x / base
... if lift > minlift:
... print('Rule {0} -> {1} has lift {2}'
... .format(antecedent, consequent,lift))
以下表格展示了部分结果。计数是指包含仅后件（即该商品被购买的基本比率）、前件中的所有项以及前件和后件中的所有项的交易数量。
| 前件 | 后件 | 后件计数 | 前件计数 | 前件和后件计数 | 提升度 |
|---|---|---|---|---|---|
| 1,378, 1,379, 1,380 | 1,269 | 279（0.3%） | 80 | 57 | 225 |
| 48, 41, 976 | 117 | 1026（1.1%） | 122 | 51 | 35 |
| 48, 41, 1,6011 | 16,010 | 1316（1.5%） | 165 | 159 | 64 |
例如，我们可以看到，有 80 笔交易中购买了 1,378、1,379 和 1,380 三个商品。在这些交易中，57 笔也包含了 1,269，因此估算的条件概率为 57/80 ≈ 71%。与所有交易中只有 0.3%包含 1,269 这一事实相比，这给我们带来了 255 的提升度。
必须在这些计数中拥有足够数量的交易，以便能够做出相对可靠的推断，这就是为什么我们必须首先选择频繁项集的原因。如果我们从不频繁项集生成规则，计数会非常小；因此，相关值将毫无意义（或者会受到非常大的误差范围的影响）。
请注意，从这个数据集中发现了更多的关联规则：该算法发现了 1,030 条规则（要求支持至少 80 个购物篮，并且最小提升度为 5）。与如今互联网所能处理的数据集相比，这仍然是一个较小的数据集。对于包含数百万笔交易的数据集，你可以预计会生成成千上万的规则，甚至是数百万条规则。
然而，对于每个客户或每个产品，只有少数几个规则在任何给定时刻是相关的。因此，每个客户只会收到少量推荐。
更高级的购物篮分析
现在有一些比 Apriori 更快的购物篮分析算法。我们之前看到的代码比较简单，足够满足我们的需求，因为我们只有大约 10 万个交易记录。如果我们有数百万条交易记录，可能值得使用更快的算法。不过需要注意的是，学习关联规则通常可以离线进行，在这种情况下效率不是那么大的问题。
还有一些方法可以处理时间序列信息，从而得出考虑到购买顺序的规则。例如，假设某人购买了大量派对用品后，可能会回来购买垃圾袋。因此，在第一次访问时推荐垃圾袋可能是有意义的。然而，向所有购买垃圾袋的人推荐派对用品就不太合适了。
总结
本章中，我们从使用回归进行评分预测开始。我们看到了一些不同的方式来实现这一目标，然后通过学习一组权重将它们结合成一个单一的预测。这种技术，特别是堆叠学习，作为一种集成学习方法，是一种可以在许多情况下使用的通用技术，不仅仅适用于回归。它允许你结合不同的思路，即使它们的内部机制完全不同；你也可以将它们的最终输出结合起来。
在本章的后半部分，我们转变了方向，研究了另一种产生推荐的模式：购物篮分析或关联规则挖掘。在这种模式下，我们试图发现“购买 X 的客户很可能对 Y 感兴趣”的（概率性）关联规则。这利用了仅从销售数据生成的数据，而无需用户对商品进行数值评分。目前，scikit-learn 中没有这个功能，所以我们编写了自己的代码。
关联规则挖掘需要小心，不能只是向每个用户推荐畅销书（否则，个性化推荐还有什么意义？）。为了做到这一点，我们学习了如何通过称为“规则提升度”的度量，衡量规则相对于基线的价值。
到目前为止，我们已经看到了机器学习的主要模式：分类。在接下来的两章中，我们将探索用于两种特定数据类型的技术，即音乐和图像。我们的第一个目标是构建一个音乐类型分类器。
第九章：分类 – 音乐流派分类
到目前为止，我们的训练数据实例每个都可以通过一组特征值向量轻松描述。例如，在 Iris 数据集中，花卉是通过包含花的某些部分长度和宽度值的向量表示的。在基于文本的示例中，我们可以将文本转化为词袋表示，并手动创建自己的特征来捕捉文本中的某些方面。
在本章中，当我们尝试根据流派对歌曲进行分类时，情况将有所不同。例如，我们应该如何表示一首三分钟的歌曲呢？我们是否应该取 MP3 表示中的每个单独的比特？可能不是，因为将其视为文本并创建类似于“声音片段袋”的东西肯定会太复杂。然而，我们仍然必须将一首歌曲转换为一系列足够描述它的值。
绘制我们的路线图
本章将展示如何在一个超出我们舒适区的领域中构建一个合理的分类器。首先，我们将不得不使用基于声音的特征，这比我们之前使用的基于文本的特征复杂得多。接着，我们将学习如何处理多分类问题，而目前为止我们只遇到过二分类问题。此外，我们还将了解一些新的分类性能衡量方法。
假设我们遇到一个场景，在某种原因下，我们的硬盘中出现了一堆随机命名的 MP3 文件，假设这些文件包含音乐。我们的任务是根据音乐流派将它们分类到不同的文件夹中，比如爵士、古典、乡村、流行、摇滚和金属。
获取音乐数据
我们将使用 GTZAN 数据集，该数据集常用于基准测试音乐流派分类任务。它包含 10 个不同的流派，我们为了简单起见，只使用其中的 6 个：古典音乐、爵士乐、乡村音乐、流行音乐、摇滚乐和金属乐。该数据集包含每个流派 100 首歌曲的前 30 秒。我们可以从opihi.cs.uvic.ca/sound/genres.tar.gz下载数据集。
提示
这些音轨以 22,050 赫兹（每秒 22,050 次采样）单声道的 WAV 格式录制。
转换为 WAV 格式
事实证明，如果我们以后想要在私人 MP3 收藏中测试我们的分类器，我们可能无法提取太多有意义的信息。这是因为 MP3 是一种有损的音乐压缩格式，它会去除人耳无法感知的部分。这种格式对存储非常友好，因为使用 MP3，你可以在设备上存储更多的歌曲。但对于我们的任务来说，这就不太合适了。为了分类，我们使用 WAV 文件会更简单一些，因为它们可以被scipy.io.wavfile包直接读取。因此，如果我们希望使用分类器处理 MP3 文件，我们就必须将它们转换成 WAV 格式。
提示
如果你附近没有转换工具，你可以查看 sox.sourceforge.net上的 SoX。它号称是声音处理的瑞士军刀，我们也同意这个大胆的说法。
将所有音乐文件保存为 WAV 格式的一个优点是，SciPy 工具包可以直接读取它：
>>> sample_rate, X = scipy.io.wavfile.read(wave_filename)
X现在包含了样本，而sample_rate是它们被采样的速率。让我们利用这些信息，快速查看一些音乐文件，初步了解数据的样子。
查看音乐
获取不同风格歌曲“外观”的一种非常方便的方法是为一个风格的歌曲集绘制频谱图。频谱图是歌曲中频率的可视化表示。它在y轴上显示频率的强度，在x轴上显示特定时间间隔的强度。也就是说，颜色越深，某一时间窗口内该频率的强度越大。
Matplotlib 提供了一个方便的函数specgram()，它为我们执行了大部分底层计算和绘图工作：
>>> import scipy
>>> from matplotlib.pyplot import specgram
>>> sample_rate, X = scipy.io.wavfile.read(wave_filename)
>>> print sample_rate, X.shape
22050, (661794,)
>>> specgram(X, Fs=sample_rate, xextent=(0,30))
我们刚刚读取的 WAV 文件的采样率为 22,050 Hz，包含 661,794 个样本。
如果我们现在为这些前 30 秒的不同 WAV 文件绘制频谱图，我们可以看到同一类型歌曲之间的共性，如下图所示：
仅从图像中，我们立刻能看出金属和古典歌曲在频谱上的差异。例如，金属歌曲在大部分频率范围内始终具有较高的强度（它们很有活力！），而古典歌曲则在时间上展示出更为多样化的模式。
应该可以训练一个分类器，以足够高的准确率区分至少金属与古典歌曲。然而，像乡村与摇滚这种其他音乐风格对比可能会更具挑战性。这对我们来说看起来像是一个真正的挑战，因为我们不仅需要区分两类，还要区分六类。我们需要能够合理地区分所有这些类别。
将音乐分解成正弦波分量
我们的计划是从原始样本数据（之前存储在X中的）中提取各个频率强度，并将它们输入分类器。这些频率强度可以通过应用所谓的快速傅里叶变换（FFT）来提取。由于傅里叶变换的理论超出了本章的范围，我们只看一个例子，直观地理解它的作用。稍后，我们将把它作为一个黑箱特征提取器。
例如，我们可以生成两个 WAV 文件，sine_a.wav和sine_b.wav，它们分别包含 400 Hz 和 3,000 Hz 的正弦波声音。前面提到的“瑞士军刀”SoX，就是实现这一目标的一种方法：
$ sox --null -r 22050 sine_a.wav synth 0.2 sine 400
$ sox --null -r 22050 sine_b.wav synth 0.2 sine 3000
在接下来的图表中，我们绘制了它们的前 0.008 秒。下面我们可以看到正弦波的 FFT。毫不奇怪，我们在对应的正弦波下看到 400 Hz 和 3,000 Hz 的尖峰。
现在，让我们将两者混合，将 400 Hz 的声音音量设置为 3,000 Hz 的音量的一半：
$ sox --combine mix --volume 1 sine_b.wav --volume 0.5 sine_a.wav sine_mix.wav
我们在合成声音的 FFT 图中看到两个尖峰，其中 3,000 Hz 的尖峰几乎是 400 Hz 的两倍大小。
对于真实的音乐，我们很快发现 FFT 不像前面的玩具示例那样漂亮：
使用 FFT 构建我们的第一个分类器
然而，我们现在可以使用 FFT 创建歌曲的某种音乐指纹。如果我们对几首歌曲这样做，并手动分配相应的音乐类型标签，我们就得到了可以输入到我们第一个分类器的训练数据。
提高实验灵活性
在深入分类器训练之前，让我们先思考一下实验的灵活性。尽管 FFT 中有“快”这个词，但它比我们在基于文本的章节中创建特征的速度要慢得多。而且由于我们仍处于实验阶段，我们可能想考虑如何加速整个特征创建过程。
当然，每次运行分类器时，创建每个文件的 FFT 都会是相同的。因此，我们可以缓存它，并读取缓存的 FFT 表示，而不是完整的 WAV 文件。我们通过create_fft()函数来实现这一点，后者又使用scipy.fft()来创建 FFT。为了简单起见（以及提高速度！），我们在这个示例中将 FFT 组件的数量固定为前 1,000 个。根据我们当前的知识，我们不知道这些是否是与音乐类型分类最相关的组件——只知道它们在前面的 FFT 示例中显示了最高的强度。如果以后我们想使用更多或更少的 FFT 组件，当然需要为每个音频文件重新创建 FFT 表示。
import os
import scipy
def create_fft(fn):
sample_rate, X = scipy.io.wavfile.read(fn)
fft_features = abs(scipy.fft(X)[:1000])
base_fn, ext = os.path.splitext(fn)
data_fn = base_fn + ".fft"
scipy.save(data_fn, fft_features)
我们使用 NumPy 的save()函数保存数据，该函数总是将.npy附加到文件名。每个 WAV 文件只需为训练或预测做一次此操作。
对应的 FFT 读取函数是read_fft()：
import glob
def read_fft(genre_list, base_dir=GENRE_DIR):
X = []
y = []
for label, genre in enumerate(genre_list):
genre_dir = os.path.join(base_dir, genre, "*.fft.npy")
file_list = glob.glob(genre_dir)
for fn in file_list:
fft_features = scipy.load(fn)
X.append(fft_features[:1000])
y.append(label)
return np.array(X), np.array(y)
在我们的音乐目录中，我们预期以下音乐类型：
genre_list = ["classical", "jazz", "country", "pop", "rock", "metal"]
训练分类器
让我们使用逻辑回归分类器，它已经在第六章中为我们提供了很好的效果，分类 II - 情感分析。增加的难度是，我们现在面临的是一个多类分类问题，而之前我们只需要区分两类。
需要提到的一个令人惊讶的方面是，从二分类问题切换到多分类问题时准确率的评估。在二分类问题中，我们已经了解到 50%的准确率是最差的情况，因为这个结果仅仅是通过随机猜测就能达到的。而在多分类设置下，50%的准确率可能已经很不错了。例如，在我们的六个类别中，随机猜测的结果大约只有 16.7%的准确率（假设类别大小相同）。
使用混淆矩阵来衡量多类问题中的准确性
在多类问题中，我们不仅要关注我们能多好地正确分类不同的类型。此外，我们还应该注意哪些类别之间我们存在混淆。这可以通过所谓的混淆矩阵来实现，如下所示：
>>> from sklearn.metrics import confusion_matrix
>>> cm = confusion_matrix(y_test, y_pred)
>>> print(cm)
[[26 1 2 0 0 2]
[ 4 7 5 0 5 3]
[ 1 2 14 2 8 3]
[ 5 4 7 3 7 5]
[ 0 0 10 2 10 12]
[ 1 0 4 0 13 12]]
该图展示了分类器对每个类别的测试集标签分布。对角线代表正确分类的结果。由于我们有六个类别，因此矩阵是六行六列。矩阵的第一行表示，对于 31 个古典音乐（第一行的总和），分类器预测了 26 个属于古典音乐，1 个属于爵士音乐，2 个属于乡村音乐，2 个属于金属音乐。对角线显示了正确分类的结果。在第一行中，我们看到，在（26+1+2+2）=31 首歌曲中，26 首被正确分类为古典音乐，5 首被误分类。实际上，这并不算太差。第二行则更为令人失望：在 24 首爵士歌曲中，只有 7 首被正确分类——也就是说，准确率只有 29%。
当然，我们遵循前几章中的训练/测试集拆分设置，因此我们实际上需要记录每个交叉验证折中的混淆矩阵。之后，我们需要对其进行平均和归一化处理，以便将结果的范围控制在 0（完全失败）和 1（全部分类正确）之间。
图形化的可视化通常比 NumPy 数组更易于阅读。matplotlib 的matshow()函数是我们的好朋友：
from matplotlib import pylab
def plot_confusion_matrix(cm, genre_list, name, title):
pylab.clf()
pylab.matshow(cm, fignum=False, cmap='Blues',
vmin=0, vmax=1.0)
ax = pylab.axes() ax.set_xticks(range(len(genre_list)))
ax.set_xticklabels(genre_list)
ax.xaxis.set_ticks_position("bottom")
ax.set_yticks(range(len(genre_list)))
ax.set_yticklabels(genre_list)
pylab.title(title)
pylab.colorbar()
pylab.grid(False)
pylab.xlabel('Predicted class')
pylab.ylabel('True class')
pylab.grid(False)
pylab.show()
提示
创建混淆矩阵时，请确保选择一个适当的颜色映射（matshow()的cmap参数），使得颜色的深浅变化能够立即显现其含义。尤其不推荐使用这些类型的图表的彩虹色图，比如 matplotlib 的默认jet色图，甚至Paired色图。
最终的图表看起来如下所示：
对于一个完美的分类器，我们期望从左上角到右下角呈现一条深色的对角线，其他区域则为浅色。在之前的图表中，我们可以立即看到我们的基于 FFT 的分类器离完美还有很大距离。它只正确预测了古典音乐（深色方块）。例如，对于摇滚音乐，它大多数时间都将标签预测为金属音乐。
显然，使用 FFT 点指向了正确的方向（古典音乐类别并不那么糟糕），但这还不足以得到一个不错的分类器。当然，我们可以调整 FFT 组件的数量（固定为 1,000）。但是在深入调整参数之前，我们应该先进行一些研究。结果表明，FFT 确实是一个不错的特征用于类别分类——只是它的精度还不够高。很快，我们将看到如何通过使用经过处理的 FFT 版本来提高分类性能。
然而，在我们进行这个分析之前，我们将学习另一种衡量分类性能的方法。
使用接收器操作特征（ROC）来衡量分类器性能的另一种方式
我们已经了解到，仅仅通过衡量准确率不足以真正评估一个分类器。相反，我们依赖于精准率-召回率（P/R）曲线来深入理解分类器的性能。
P/R 曲线有一个姐妹曲线，叫做接收器操作特征（ROC），它衡量分类器性能的相似方面，但提供了另一种分类性能的视角。两者的关键区别在于，P/R 曲线更适合于正类比负类更为重要，或者正类样本远少于负类样本的任务。信息检索和欺诈检测是典型的应用领域。另一方面，ROC 曲线则更好地展示了分类器整体表现。
为了更好地理解这些差异，假设我们来看一下先前训练的分类器在正确分类乡村歌曲方面的表现，如下图所示：
在左侧，我们可以看到 P/R 曲线。对于一个理想的分类器，我们希望看到的曲线是从左上角直接到右上角，再到右下角，从而形成一个面积（AUC）为 1.0 的曲线。
右侧的图展示了相应的 ROC 曲线。它绘制了真正例率与假正例率之间的关系。在这里，一个理想的分类器曲线应该从左下角走到左上角，然后到达右上角。而一个随机分类器则表现为从左下角到右上角的直线，如虚线所示，AUC 为 0.5。因此，我们不能将 P/R 曲线的 AUC 与 ROC 曲线的 AUC 进行比较。
独立于曲线，在比较同一数据集上两个不同分类器时，我们总是可以安全地假设，一个分类器的 P/R 曲线的 AUC 较高，也意味着其对应的 ROC 曲线的 AUC 较高，反之亦然。因此，我们通常不会生成两者。关于这一点的更多信息可以在 Davis 和 Goadrich（ICML，2006）撰写的非常有见地的论文《精准率-召回率与 ROC 曲线的关系》中找到。
以下表格总结了 P/R 曲线与 ROC 曲线之间的差异：
| x 轴 | y 轴 | |
|---|---|---|
| P/R | ||
| ROC |
看着这两条曲线的* x 轴和 y 轴的定义，我们可以看到，ROC 曲线的 y 轴上的真正阳性率与 P/R 图的 x *轴上的召回率是相同的。
假阳性率衡量被错误识别为阳性的真正负例的比例，在完美的情况下为 0（没有假阳性），否则为 1。与此对比，精准度则关注完全相反的内容，即我们正确分类为阳性的真正例的比例。
以后，让我们使用 ROC 曲线来衡量分类器的性能，以便更好地感知其效果。我们多类问题的唯一挑战是，ROC 和 P/R 曲线假设的是二分类问题。因此，为了我们的目的，我们将为每个音乐类型创建一张图表，展示分类器在一对其余类别分类中的表现：
from sklearn.metrics import roc_curve
y_pred = clf.predict(X_test)
for label in labels:
y_label_test = scipy.asarray(y_test==label, dtype=int)
proba = clf.predict_proba(X_test)
proba_label = proba[:,label]
# calculate false and true positive rates as well as the
# ROC thresholds
fpr, tpr, roc_thres = roc_curve(y_label_test, proba_label)
# plot tpr over fpr ...
结果是以下六个 ROC 图。如我们已经发现的，我们的第一个版本的分类器只对古典歌曲表现良好。然而，查看个别的 ROC 曲线告诉我们，大部分其他类型的表现确实不佳。只有爵士乐和乡村音乐带来了一些希望。其余的类型显然无法使用。
使用梅尔频率倒谱系数提高分类性能
我们已经了解到，FFT 指向了正确的方向，但它本身不足以最终得到一个成功的分类器，能够将我们乱序的、包含多种音乐类型的歌曲目录整理到各个单独的类型目录中。我们需要一个稍微更高级的版本。
此时，承认我们需要做更多的研究总是明智的。其他人可能曾面临类似的挑战，并且已经找到了新的方法，可能也能帮助我们。事实上，甚至每年都有一场专注于音乐类型分类的会议，由国际音乐信息检索学会（ISMIR）组织。显然，自动音乐类型分类（AMGC）已经成为音乐信息检索的一个成熟子领域。快速浏览一些 AMGC 的论文，我们看到有很多针对自动类型分类的工作，可能会帮助我们。
在许多研究中似乎成功应用的一个技术叫做梅尔频率倒谱系数（Mel Frequency Cepstral Coefficients）。梅尔频率倒谱（MFC）编码了声音的功率谱，即声音包含的每个频率的功率。它是通过对信号谱的对数进行傅里叶变换来计算的。如果这听起来太复杂，简单记住，“倒谱”这个名字源于“谱”（spectrum）一词的前四个字母倒过来。MFC 已被成功应用于语音和说话人识别。我们来看看它是否也能在我们的案例中发挥作用。
我们正处于一个幸运的情况，因为别人已经恰好需要这个并且发布了一个实现，叫做 Talkbox SciKit。我们可以从pypi.python.org/pypi/scikits.talkbox安装它。之后，我们可以调用mfcc()函数来计算 MFC 系数，方法如下：
>>> from scikits.talkbox.features import mfcc
>>> sample_rate, X = scipy.io.wavfile.read(fn)
>>> ceps, mspec, spec = mfcc(X)
>>> print(ceps.shape)
(4135, 13)
我们希望输入到分类器中的数据存储在ceps中，它包含了每个歌名为fn的歌曲的 4,135 帧的 13 个系数（mfcc()函数的nceps参数的默认值）。如果直接使用所有数据，会使分类器过载。相反，我们可以对每个系数在所有帧上进行平均。假设每首歌的开始和结束部分可能不如中间部分具有明显的音乐类型特征，我们也忽略了前后各 10%的数据。
x = np.mean(ceps[int(num_ceps*0.1):int(num_ceps*0.9)], axis=0)
不出所料，我们将使用的基准数据集只包含每首歌的前 30 秒，因此我们不需要剪掉最后 10%。不过我们还是这么做了，以确保我们的代码可以在其他可能没有截断的数据集上运行。
类似于我们使用 FFT 的工作，我们当然也希望缓存一次生成的 MFCC 特征，并在每次训练分类器时读取它们，而不是每次都重新生成。
这导致了以下代码：
def write_ceps(ceps, fn):
base_fn, ext = os.path.splitext(fn)
data_fn = base_fn + ".ceps"
np.save(data_fn, ceps)
print("Written to %s" % data_fn)
def create_ceps(fn):
sample_rate, X = scipy.io.wavfile.read(fn)
ceps, mspec, spec = mfcc(X)
write_ceps(ceps, fn)
def read_ceps(genre_list, base_dir=GENRE_DIR):
X, y = [], []
for label, genre in enumerate(genre_list):
for fn in glob.glob(os.path.join(
base_dir, genre, "*.ceps.npy")):
ceps = np.load(fn)
num_ceps = len(ceps)
X.append(np.mean(
ceps[int(num_ceps*0.1):int(num_ceps*0.9)], axis=0))
y.append(label)
return np.array(X), np.array(y)
我们通过一个每首歌只使用 13 个特征的分类器得到了以下有希望的结果：
所有音乐类型的分类表现都有所提升。古典音乐和金属音乐的 AUC 几乎达到了 1.0。实际上，下面的混淆矩阵现在看起来好多了。我们可以清楚地看到对角线，显示出分类器在大多数情况下能够正确地分类各个音乐类型。这个分类器实际上相当适用于解决我们的初始任务。
如果我们想在这方面有所改进，这个混淆矩阵会迅速告诉我们需要关注的地方：非对角线位置的非白色区域。例如，我们在一个较暗的区域中错误地将摇滚歌曲标记为爵士乐，且这种错误有相当大的概率。要解决这个问题，我们可能需要更深入地研究这些歌曲，提取诸如鼓点模式和类似的音乐风格特征。然后——在浏览 ISMIR 论文时——我们还读到了一种名为听觉滤波器带时域包络（AFTE）的特征，似乎在某些情况下优于 MFCC 特征。或许我们也该看看它们？
好消息是，只要配备了 ROC 曲线和混淆矩阵，我们可以随时借用其他专家在特征提取器方面的知识，而无需完全理解它们的内部工作原理。我们的测量工具总是会告诉我们，何时方向正确，何时需要改变。当然，作为一个渴望学习的机器学习者，我们总会有一种模糊的感觉，觉得在特征提取器的黑箱中埋藏着一个激动人心的算法，正等着我们去理解。
总结
在本章中，我们在构建音乐类型分类器时完全走出了舒适区。由于对音乐理论没有深入理解，最初我们未能训练出一个合理准确地预测歌曲音乐类型的分类器。通过快速傅里叶变换（FFT）进行尝试失败。但随后，我们使用 MFC 特征创建了一个表现出真正可用性能的分类器。
在这两种情况下，我们使用的特征我们只是了解得足够多，知道如何以及在哪里将它们放入分类器设置中。一个失败了，另一个成功了。它们之间的区别在于，在第二种情况下，我们依赖的是领域专家创建的特征。
这完全没问题。如果我们主要关心结果，有时候我们只是需要走捷径——我们只需要确保这些捷径来自于特定领域的专家。而且，因为我们学会了如何在这个新的多类别分类问题中正确地衡量性能，所以我们能够充满信心地走这些捷径。
在下一章，我们将探讨如何将你在本书其余部分学到的技术应用到这种特定类型的数据中。我们将学习如何使用 mahotas 计算机视觉包，通过传统的图像处理功能来预处理图像。
第十章：计算机视觉
图像分析和计算机视觉一直在工业和科学应用中占有重要地位。随着拥有强大相机和互联网连接的手机的普及，图像现在越来越多地由消费者生成。因此，利用计算机视觉为新情境中的用户体验提供更好的服务，成为了一个机会。
在本章中，我们将学习如何将你在本书其余部分学到的技术应用到这种特定类型的数据上。具体来说，我们将学习如何使用 mahotas 计算机视觉包从图像中提取特征。这些特征可以作为我们在其他章节中学习的分类方法的输入。我们将这些技术应用于公开可用的照片数据集。我们还将看到这些相同的特征如何用于另一个问题，即寻找相似图像的问题。
最后，在本章的结尾，我们将学习如何使用局部特征。这些方法相对较新（其中第一种取得最新成果的方法，即尺度不变特征变换（SIFT），是在 1999 年提出的），并在许多任务中取得了非常好的效果。
介绍图像处理
从计算机的角度来看，图像是一个由像素值组成的大型矩形数组。我们的目标是处理这个图像，并为我们的应用程序做出决策。
第一步是从磁盘加载图像，图像通常以 PNG 或 JPEG 等图像特定格式存储，前者是一种无损压缩格式，后者是一种有损压缩格式，优化了对照片的视觉评估。然后，我们可能希望对图像进行预处理（例如，对其进行归一化，以适应光照变化）。
本章将以分类问题为驱动。我们希望能够学习一个支持向量机（或其他）分类器，能够从图像中进行训练。因此，我们将使用一个中间表示，在应用机器学习之前，从图像中提取数值特征。
加载和显示图像
为了操作图像，我们将使用一个名为 mahotas 的软件包。你可以从pypi.python.org/pypi/mahotas获取 mahotas，并在mahotas.readthedocs.org阅读其手册。Mahotas 是一个开源软件包（MIT 许可，因此可以在任何项目中使用），由本书的作者之一开发。幸运的是，它是基于 NumPy 的。你迄今为止学到的 NumPy 知识可以用于图像处理。还有其他图像处理软件包，例如 scikit-image（skimage）、SciPy 中的 ndimage（n 维图像）模块和 OpenCV 的 Python 绑定。所有这些都可以原生支持 NumPy 数组，因此你甚至可以将来自不同软件包的功能混合使用，构建一个综合的处理管道。
我们首先导入 mahotas 库，并使用mh作为缩写，这将在本章中一直使用，如下所示：
>>> import mahotas as mh
现在，我们可以使用imread来加载图像文件，如下所示：
>>> image = mh.imread('scene00.jpg')
scene00.jpg文件（该文件包含在本书附带的代码仓库数据集中）是一个高度为h、宽度为w的彩色图像；该图像将是一个形状为(h, w, 3)的数组。第一个维度是高度，第二个是宽度，第三个是红/绿/蓝。其他系统将宽度放在第一个维度，但这是所有基于 NumPy 的包使用的惯例。数组的类型通常是np.uint8（无符号 8 位整数）。这些是您的相机拍摄的图像，或者是您的显示器能够完全显示的图像。
一些专业设备，常用于科学和技术应用中，可以拍摄具有更高位深度的图像（即对亮度变化更加敏感）。这种设备中，12 位或 16 位是常见的位深。Mahotas 可以处理所有这些类型的图像，包括浮点图像。在许多计算中，即使原始数据是无符号整数，转换为浮点数仍然有助于简化处理舍入和溢出问题。
注意
Mahotas 支持多种不同的输入/输出后端。遗憾的是，它们没有一个能够加载所有存在的图像格式（存在数百种格式，并且每种格式都有多个变体）。然而，PNG 和 JPEG 格式的图像是所有后端都支持的。我们将重点介绍这些常见格式，并参考 mahotas 文档，了解如何读取不常见的格式。
我们可以使用 matplotlib 来显示图像，这个我们已经多次使用过的绘图库，如下所示：
>>> from matplotlib import pyplot as plt
>>> plt.imshow(image)
>>> plt.show()
如下所示，这段代码展示了图像，遵循的惯例是第一个维度为高度，第二个维度为宽度。它也正确处理彩色图像。当使用 Python 进行数值计算时，我们可以受益于整个生态系统的良好协同工作：mahotas 与 NumPy 数组兼容，这些数组可以通过 matplotlib 显示；稍后我们将从图像中计算特征，以便与 scikit-learn 一起使用。
阈值处理
阈值处理是一种非常简单的操作：我们将所有高于某个阈值的像素值转换为1，将低于该阈值的像素值转换为0（或者使用布尔值，将其转换为True和False）。阈值处理中的一个重要问题是选择一个合适的阈值作为分界线。Mahotas 实现了一些从图像中选择阈值的方法。其中一种方法叫做Otsu，以其发明者命名。第一个必要的步骤是使用mahotas.colors子模块中的rgb2gray将图像转换为灰度图像。
我们也可以使用红色、绿色和蓝色通道的平均值来代替 rgb2gray，方法是调用 image.mean(2)。然而，结果会有所不同，因为 rgb2gray 为不同的颜色使用了不同的权重，以得到更为主观上更令人愉悦的结果。我们的眼睛对三种基本颜色的敏感度并不相同。
>>> image = mh.colors.rgb2grey(image, dtype=np.uint8)
>>> plt.imshow(image) # Display the image
默认情况下，matplotlib 会将这张单通道图像显示为伪彩色图像，使用红色表示高值，蓝色表示低值。对于自然图像，灰度图像更加合适。你可以通过以下方式选择灰度显示：
>>> plt.gray()
现在图像已显示为灰度图像。请注意，仅仅是像素值的解释和显示方式发生了变化，图像数据本身没有被改变。我们可以继续处理，通过计算阈值：
>>> thresh = mh.thresholding.otsu(image)
>>> print('Otsu threshold is {}.'.format(thresh))
Otsu threshold is 138.
>>> plt.imshow(image > thresh)
当应用到之前的截图时，这种方法找到了阈值为 138，它将地面与上方的天空分开，如下图所示：
高斯模糊
对图像进行模糊处理可能看起来很奇怪，但它通常能减少噪声，这有助于进一步的处理。使用 mahotas，它只需要一个函数调用：
>>> im16 = mh.gaussian_filter(image, 16)
注意，我们并没有将灰度图像转换为无符号整数：我们只是直接使用了浮动点结果。gaussian_filter 函数的第二个参数是滤波器的大小（即滤波器的标准差）。较大的值会导致更多的模糊，如下图所示：
我们可以使用左侧的截图并结合 Otsu 阈值方法（使用之前相同的代码）。现在，边界更加平滑，没有了锯齿状边缘，如下图所示：
让中心对焦
最后的示例展示了如何将 NumPy 运算符与一些简单的滤波操作结合，得到有趣的效果。我们从 Lena 图像开始，将其分割成颜色通道：
>>> im = mh.demos.load('lena')
这是一张年轻女性的图像，经常用于图像处理演示。它在下图中展示：
为了分割红色、绿色和蓝色通道，我们使用以下代码：
>>> r,g,b = im.transpose(2,0,1)
现在，我们分别对三个通道进行滤波，并通过 mh.as_rgb 将它们组合成一张图像。此函数接受三个二维数组，进行对比度拉伸使每个数组成为 8 位整数数组，然后将它们堆叠，返回一张彩色 RGB 图像：
>>> r12 = mh.gaussian_filter(r, 12.)
>>> g12 = mh.gaussian_filter(g, 12.)
>>> b12 = mh.gaussian_filter(b, 12.)
>>> im12 = mh.as_rgb(r12, g12, b12)
现在，我们将两张图像从中心到边缘进行混合。首先，我们需要构建一个权重数组 W，它在每个像素处包含一个归一化值，即该像素到中心的距离：
>>> h, w = r.shape # height and width
>>> Y, X = np.mgrid[:h,:w]
我们使用了 np.mgrid 对象，它返回大小为 (h, w) 的数组，值对应于 y 和 x 坐标。接下来的步骤如下：
>>> Y = Y - h/2\. # center at h/2
>>> Y = Y / Y.max() # normalize to -1 .. +1
>>> X = X - w/2.
>>> X = X / X.max()
现在我们使用一个高斯函数，给中心区域赋予一个高值：
>>> C = np.exp(-2.*(X**2+ Y**2))
>>> # Normalize again to 0..1
>>> C = C - C.min()
>>> C = C / C.ptp()
>>> C = C[:,:,None] # This adds a dummy third dimension to C
请注意，所有这些操作都是使用 NumPy 数组而不是 mahotas 特定的方法来完成的。最后，我们可以将两张图片结合起来，使得中心部分聚焦清晰，而边缘部分则更柔和：
>>> ringed = mh.stretch(im*C + (1-C)*im12)
基本图像分类
我们将从一个专门为本书收集的小型数据集开始。该数据集包含三个类别：建筑物、自然场景（风景）和文本图片。每个类别有 30 张图片，所有图片均使用手机摄像头拍摄，构图简单。图像与那些没有摄影训练的用户上传到现代网站上的图片类似。这个数据集可以通过本书的网站或 GitHub 代码库获得。在本章后面，我们将介绍一个更难的数据集，包含更多的图片和类别。
在进行图像分类时，我们从一个包含大量数字（像素值）的矩阵开始。现在，百万级像素已是常见情况。我们可以尝试将所有这些数字作为特征输入学习算法，但这并不是一个很好的主意。原因是，每个像素（甚至每小组像素）与最终结果之间的关系非常间接。此外，如果像素数有百万，但样本图像的数量很少，这将导致一个非常困难的统计学习问题。这是我们在第七章《回归》中讨论的 P 大于 N 类型问题的极端形式。相反，一个好的方法是从图像中计算特征，并利用这些特征进行分类。
话虽如此，我还是要指出，实际上，有一些方法是可以直接从像素值中工作并计算特征的。这些方法有自己的特征计算子模块，甚至可能尝试自动学习合适的特征。这些方法是当前研究的主题，通常在非常大的数据集（数百万张图片）上效果最佳。
我们之前使用了一个场景类别的例子。以下是文本类和建筑类的示例：
从图像中计算特征
使用 mahotas 时，从图像中计算特征非常简单。它有一个名为mahotas.features的子模块，提供了特征计算函数。
一个常用的纹理特征集是 Haralick 特征。和许多图像处理方法一样，这个名字来源于其发明者。这些特征是基于纹理的：它们能区分平滑的图像和有图案的图像，并且能区分不同的图案。使用 mahotas 时，计算这些特征非常简单，方法如下：
>>> haralick_features = mh.features.haralick(image)
>>> haralick_features_mean = np.mean(haralick_features, axis=0)
>>> haralick_features_all = np.ravel(haralick_features)
mh.features.haralick函数返回一个 4x13 的数组。第一维表示计算特征的四个可能方向（垂直、水平、对角线和反对角线）。如果我们对方向没有特别兴趣，我们可以使用所有方向的平均值（如前面的代码中的haralick_features_mean）。否则，我们可以分别使用所有特征（使用haralick_features_all）。这个决策应根据数据集的属性来决定。在我们的案例中，我们推测水平和垂直方向应该分别保留。因此，我们将使用haralick_features_all。
mahotas 中还实现了一些其他的特征集。线性二进制模式是另一种基于纹理的特征集，它对光照变化非常稳健。还有其他类型的特征，包括局部特征，我们将在本章后面讨论。
使用这些特征时，我们采用标准的分类方法，如逻辑回归，具体如下：
>>> from glob import glob
>>> images = glob('SimpleImageDataset/*.jpg')
>>> features = []
>>> labels = []
>>> for im in images:
... labels.append(im[:-len('00.jpg')])
... im = mh.imread(im)
... im = mh.colors.rgb2gray(im, dtype=np.uint8)
... features.append(mh.features.haralick(im).ravel())
>>> features = np.array(features)
>>> labels = np.array(labels)
这三个类别有非常不同的纹理。建筑物有锐利的边缘和大块区域，颜色相似（像素值很少完全相同，但变化很小）。文本由许多锐利的暗-亮过渡组成，黑色小区域在白色背景中。自然场景有更平滑的变化，带有类似分形的过渡。因此，基于纹理的分类器预计会表现良好。
作为分类器，我们将使用逻辑回归分类器，并对特征进行如下预处理：
>>> from sklearn.pipeline import Pipeline
>>> from sklearn.preprocessing import StandardScaler
>>> from sklearn.linear_model import LogisticRegression
>>> clf = Pipeline([('preproc', StandardScaler()),
('classifier', LogisticRegression())])
由于我们的数据集较小，我们可以使用如下的留一法回归：
>>> from sklearn import cross_validation
>>> cv = cross_validation.LeaveOneOut(len(images))
>>> scores = cross_validation.cross_val_score(
... clf, features, labels, cv=cv)
>>> print('Accuracy: {:.1%}'.format(scores.mean()))
Accuracy: 81.1%
对于三个类别，81%的准确率不错（随机猜测的准确率大约为 33%）。然而，通过编写我们自己的特征，我们可以做得更好。
编写你自己的特征
特征没有什么神奇之处。它只是我们从图像中计算出的一个数字。文献中已经定义了几种特征集。这些特征集通常有一个额外的优势，即它们被设计和研究为对许多不重要的因素具有不变性。例如，线性二进制模式对将所有像素值乘以一个数字或将一个常数加到所有像素值上完全不变。这使得这个特征集对图像的光照变化具有稳健性。
然而，也有可能你的特定使用案例将从一些特别设计的功能中受益。
mahotas 中没有提供的一种简单功能是颜色直方图。幸运的是，这个功能很容易实现。颜色直方图将颜色空间划分为一组区间，然后计算每个区间内有多少个像素。
图像采用 RGB 格式，即每个像素有三个值：R 表示红色，G 表示绿色，B 表示蓝色。由于这些组件每个都是 8 位值，总共有 1700 万种不同的颜色。我们将通过将颜色分组到箱子中，将这个数字减少到仅 64 种颜色。我们将编写一个函数来封装这个算法，如下所示：
def chist(im):
为了对颜色进行分箱，我们首先将图像除以 64，按如下方式向下取整像素值：
im = im // 64
这使得像素值的范围从 0 到 3，总共产生了 64 种不同的颜色。
按照以下方式分离红色、绿色和蓝色通道：
r,g,b = im.transpose((2,0,1))
pixels = 1 * r + 4 * b + 16 * g
hist = np.bincount(pixels.ravel(), minlength=64)
hist = hist.astype(float)
转换为对数尺度，如下代码片段所示。这并非严格必要，但能够产生更好的特征。我们使用 np.log1p，它计算 log(h+1)。这确保了零值保持为零值（从数学角度看，零的对数是未定义的，如果你尝试计算，NumPy 会给出警告）。
hist = np.log1p(hist)
return hist
我们可以很容易地将之前的处理代码调整为使用我们编写的函数：
>>> features = []
>>> for im in images:
... image = mh.imread(im)
... features.append(chist(im))
使用我们之前使用的相同交叉验证代码，我们获得了 90% 的准确率。然而，最佳结果来自于将所有特征结合在一起，我们可以按如下方式实现：
>>> features = []
>>> for im in images:
... imcolor = mh.imread(im)
... im = mh.colors.rgb2gray(imcolor, dtype=np.uint8)
... features.append(np.concatenate([
... mh.features.haralick(im).ravel(),
... chist(imcolor),
... ]))
通过使用所有这些特征，我们获得了 95.6% 的准确率，如以下代码片段所示：
>>> scores = cross_validation.cross_val_score(
... clf, features, labels, cv=cv)
>>> print('Accuracy: {:.1%}'.format(scores.mean()))
Accuracy: 95.6%
这是一个完美的例子，说明了好算法是容易的部分。你总是可以使用 scikit-learn 中的最先进的分类实现。真正的秘密和附加价值通常体现在特征设计和工程上。这就是你对数据集的知识变得非常有价值的地方。
使用特征来查找相似图像
通过用相对较少的特征表示图像的基本概念，不仅可以用于分类。例如，我们还可以用它来查找与给定查询图像相似的图像（就像我们之前在文本文档中做的那样）。
我们将计算与之前相同的特征，唯一不同的是：我们将忽略图像的边缘区域。原因是由于构图的业余性质，图像的边缘常常包含无关的元素。当特征在整个图像上计算时，这些元素会被考虑在内。通过简单地忽略它们，我们可以得到稍微更好的特征。在监督学习示例中，这并不那么重要，因为学习算法会学习哪些特征更具信息量并相应地加权。当以无监督的方式工作时，我们需要更加小心，确保我们的特征能够捕捉到数据中重要的元素。这在循环中实现如下：
>>> features = []
>>> for im in images:
... imcolor = mh.imread(im)
... # ignore everything in the 200 pixels closest to the borders
... imcolor = imcolor[200:-200, 200:-200]
... im = mh.colors.rgb2gray(imcolor, dtype=np.uint8)
... features.append(np.concatenate([
... mh.features.haralick(im).ravel(),
... chist(imcolor),
... ]))
我们现在标准化特征并计算距离矩阵，如下所示：
>>> sc = StandardScaler()
>>> features = sc.fit_transform(features)
>>> from scipy.spatial import distance
>>> dists = distance.squareform(distance.pdist(features))
我们将仅绘制数据的一个子集（每第十个元素），这样查询图像会显示在顶部，返回的“最近邻”图像显示在底部，如下所示：
>>> fig, axes = plt.subplots(2, 9)
>>> for ci,i in enumerate(range(0,90,10)):
... left = images[i]
... dists_left = dists[i]
... right = dists_left.argsort()
... # right[0] is same as left[i], so pick next closest
... right = right[1]
... right = images[right]
... left = mh.imread(left)
... right = mh.imread(right)
... axes[0, ci].imshow(left)
... axes[1, ci].imshow(right)
结果如以下截图所示：
很明显，系统并不完美，但至少可以找到与查询在视觉上相似的图像。在除了一个案例外，找到的图像都来自与查询相同的类别。
分类更困难的数据集
前一个数据集是一个使用纹理特征进行分类的简单数据集。事实上，从商业角度来看，许多有趣的问题相对简单。然而，有时我们可能会面临一个更困难的问题，需要更好、更现代的技术来获得良好的结果。
我们现在将测试一个公共数据集，它具有相同的结构：若干张照片分为少数几个类别。这些类别是动物、汽车、交通工具和自然场景。
与我们之前讨论的三类问题相比，这些类别更难区分。自然场景、建筑物和文本的纹理差异很大。然而，在这个数据集中，纹理和颜色不再是图像类别的明显标记。以下是来自动物类的一个例子：
这是来自汽车类别的另一个例子：
这两个对象背景都是自然背景，并且对象内部有较大的平滑区域。这是比简单数据集更具挑战性的问题，因此我们需要使用更先进的方法。第一个改进是使用一个稍微更强大的分类器。scikit-learn 提供的逻辑回归是一个带惩罚项的逻辑回归，包含一个可调参数C。默认情况下，C = 1.0，但这可能不是最佳选择。我们可以使用网格搜索来找到这个参数的最佳值，具体方法如下：
>>> from sklearn.grid_search import GridSearchCV
>>> C_range = 10.0 ** np.arange(-4, 3)
>>> grid = GridSearchCV(LogisticRegression(), param_grid={'C' : C_range})
>>> clf = Pipeline([('preproc', StandardScaler()),
... ('classifier', grid)])
数据集中的数据并不是按随机顺序组织的：相似的图像彼此靠近。因此，我们使用交叉验证策略，考虑数据已被洗牌，这样每个折叠（fold）都有一个更具代表性的训练集，如下所示：
>>> cv = cross_validation.KFold(len(features), 5,
... shuffle=True, random_state=123)
>>> scores = cross_validation.cross_val_score(
... clf, features, labels, cv=cv)
>>> print('Accuracy: {:.1%}'.format(scores.mean()))
Accuracy: 72.1%
对于四个类别来说，这还算不错，但我们现在将看看是否能通过使用不同的特征集来做得更好。事实上，我们将看到，结合这些特征与其他方法，才能获得最好的结果。
局部特征表示
计算机视觉领域的一个相对较新的发展是基于局部特征的方法。局部特征是在图像的小区域内计算的，与我们之前考虑的在整个图像上计算的特征不同。Mahotas 支持计算一种这种特征，加速稳健特征（SURF）。还有其他几种，最著名的是 SIFT 的原始提案。这些特征旨在对旋转或光照变化具有鲁棒性（即，它们在光照变化时仅略微改变其值）。
使用这些特征时，我们必须决定在哪里计算它们。常用的三种计算位置如下：
-
随机地
-
在网格中
-
检测图像中的有趣区域（这是一种被称为关键点检测或兴趣点检测的技术）
所有这些方法都是有效的，在合适的情况下，会得到很好的结果。Mahotas 支持这三种方法。如果你有理由认为你的兴趣点将对应于图像中重要区域，那么使用兴趣点检测效果最好。
我们将使用兴趣点方法。使用 mahotas 计算特征非常简单：导入正确的子模块并调用 surf.surf 函数，如下所示：
>>> from mahotas.features import surf
>>> image = mh.demos.load('lena')
>>> image = mh.colors.rgb2gray(im, dtype=np.uint8)
>>> descriptors = surf.surf(image, descriptor_only=True)
descriptors_only=True 标志意味着我们只关心描述符本身，而不是它们的像素位置、大小或方向。或者，我们也可以使用密集采样方法，使用 surf.dense 函数如下所示：
>>> from mahotas.features import surf
>>> descriptors = surf.dense(image, spacing=16)
这返回计算在相距 16 像素的点上的描述符值。由于点的位置是固定的，兴趣点的元信息并不十分有趣，默认情况下不会返回。在任何情况下，结果（描述符）是一个 n×64 的数组，其中 n 是采样点的数量。点的数量取决于图像的大小、内容以及你传递给函数的参数。在这个例子中，我们使用的是默认设置，每张图像得到几百个描述符。
我们不能直接将这些描述符输入到支持向量机、逻辑回归器或类似的分类系统中。为了使用来自图像的描述符，有几种解决方案。我们可以直接对它们求平均，但这样做的结果并不好，因为这样丢弃了所有关于位置的特定信息。在那种情况下，我们只会得到一个基于边缘度量的全局特征集。
我们将在这里使用的解决方案是词袋模型，这是一种非常新的思想。它首次于 2004 年以这种形式发布。这是一个事后看非常显而易见的想法：它非常简单实现，并且能够取得非常好的结果。
在处理图像时，说到单词可能会显得有些奇怪。如果你理解为你没有写出那些彼此容易区分的单词，而是口头发出的音频，可能更容易理解。现在，每次说出一个单词时，它的发音会略有不同，而且不同的发言者会有自己的发音方式。因此，一个单词的波形每次说出来时都不会完全相同。然而，通过对这些波形进行聚类，我们可以期望恢复大部分结构，使得给定单词的所有实例都在同一个聚类中。即使过程不完美（而且确实不会完美），我们仍然可以谈论将波形分组成单词。
我们对图像数据执行相同的操作：我们将所有图像中看起来相似的区域聚类在一起，并将这些区域称为视觉词汇。
注意
使用的词数通常对算法的最终性能影响不大。自然，如果词数非常少（例如 10 或 20，当你拥有几千张图像时），那么整体系统的表现会很差。同样，如果你有太多的词（例如，比图像的数量多得多），系统的表现也会不佳。然而，在这两者之间，通常会有一个较大的平稳区间，在这个区间内，你可以选择词数而不会对结果产生太大影响。作为经验法则，如果你有很多图像，使用 256、512 或 1024 这样的值应该能给你一个不错的结果。
我们将通过以下方式开始计算特征：
>>> alldescriptors = []
>>> for im in images:
... im = mh.imread(im, as_grey=True)
... im = im.astype(np.uint8)
... alldescriptors.append(surf.dense(image, spacing=16))
>>> # get all descriptors into a single array
>>> concatenated = np.concatenate(alldescriptors)
>>> print('Number of descriptors: {}'.format(
... len(concatenated)))
Number of descriptors: 2489031
这导致了超过 200 万个局部描述符。现在，我们使用 k-means 聚类来获取质心。我们可以使用所有描述符，但为了提高速度，我们将使用较小的样本，如下所示：
>>> # use only every 64th vector
>>> concatenated = concatenated[::64]
>>> from sklearn.cluster import KMeans
>>> k = 256
>>> km = KMeans(k)
>>> km.fit(concatenated)
完成此操作后（这需要一些时间），km对象包含有关质心的信息。我们现在回到描述符，并按以下方式构建特征向量：
>>> sfeatures = []
>>> for d in alldescriptors:
... c = km.predict(d)
... sfeatures.append(
... np.array([np.sum(c == ci) for ci in range(k)])
... )
>>> # build single array and convert to float
>>> sfeatures = np.array(sfeatures, dtype=float)
这个循环的最终结果是，sfeatures[fi, fj]表示图像fi中包含元素fj的次数。虽然使用np.histogram函数可以更快地计算这个值，但正确设置参数会有些棘手。我们将结果转换为浮点数，因为我们不想使用整数运算（及其舍入语义）。
结果是，每个图像现在都由一个大小相同的特征数组表示（在我们的例子中是 256 个聚类）。因此，我们可以按以下方式使用我们的标准分类方法：
>>> scores = cross_validation.cross_val_score(
... clf, sfeatures, labels, cv=cv)
>>> print('Accuracy: {:.1%}'.format(scores.mean()))
Accuracy: 62.6%
结果比之前更差！我们什么也没得到吗？
事实上，我们有，因为我们可以将所有特征结合起来，达到 76.1%的准确率，具体如下：
>>> combined = np.hstack([features, features])
>>> scores = cross_validation.cross_val_score(
... clf, combined, labels, cv=cv)
>>> print('Accuracy: {:.1%}'.format(scores.mean()))
Accuracy: 76.1%
这是我们获得的最佳结果，优于任何单一的特征集。之所以如此，是因为局部的 SURF 特征与我们之前的全局图像特征有足够的区别，能够提供新的信息，从而提升了综合结果。
总结
我们学习了在机器学习环境中处理图像的经典特征方法：通过将数百万个像素转换为少量数值特征，我们能够直接使用逻辑回归分类器。我们在其他章节中学到的所有技术突然变得可以直接应用于图像问题。我们在使用图像特征查找数据集中相似图像时，看到了一个例子。
我们还学习了如何在词袋模型中使用局部特征进行分类。这是一种非常现代的计算机视觉方法，在提高分类精度的同时，对图像中的许多无关因素（如光照，甚至同一图像中的不均匀光照）具有很强的鲁棒性。我们还使用了聚类作为分类中的一个有用中间步骤，而不是仅仅作为最终目标。
我们专注于 mahotas，它是 Python 中主要的计算机视觉库之一。还有其他同样维护良好的库。Skimage（scikit-image）在精神上类似，但有不同的功能集。OpenCV 是一个非常优秀的 C++库，具有 Python 接口。这些库都可以与 NumPy 数组配合使用，你可以混合搭配不同库中的函数，构建复杂的计算机视觉管道。
在下一章中，你将学习一种不同形式的机器学习：降维。如我们在前几章中所见，包括在本章中使用图像时，计算上生成许多特征是非常容易的。然而，通常我们希望通过减少特征数量来提高速度和可视化，或者改善我们的结果。在下一章中，我们将看到如何实现这一目标。
