构建-Python-机器学习系统-三-

布客飞龙
105 阅读1小时+

构建 Python 机器学习系统（三）

原文：Building Machine Learning Systems with Python

协议：CC BY-NC-SA 4.0

九、分类二——情感分析

对于公司来说，密切监控公众对关键事件的接受是至关重要的，例如产品发布或新闻发布。随着用户生成的内容在推特上的实时访问和容易访问，现在可以对推文进行情感分类。有时也被称为意见挖掘，这是一个活跃的研究领域，几个公司已经在销售这种服务。由于这表明显然存在市场，我们有动力使用上一章构建的分类肌肉来构建我们自己的本土情感分类器。

绘制我们的路线图

由于推特的每条消息的大小限制，对推文的情绪分析尤其困难。这导致了一种特殊的句法，创造性的缩写，以及很少格式良好的句子。典型的方法是分析句子，汇总每个段落的情感信息，然后计算文档的整体情感，这种方法在这里不起作用。

显然，我们不会试图构建最先进的情感分类器。相反，我们希望执行以下操作:

  • 以此场景为载体，引入另一种分类算法朴素贝叶斯
  • 解释词性 ( 词性)标注是如何工作的，以及它如何帮助我们
  • 展示 scikit-learn 工具箱中更多有用的技巧

获取推特数据

自然，我们需要描述情感的推文和相应的标签。在这一章中，我们将使用来自 Niek Sanders 的语料库，他出色地完成了手动将 5000 多条推文标记为积极、消极或中立的工作，并授予我们在这一章中使用它的权限。

为了遵守推特的服务条款，我们不会提供任何来自推特的数据，也不会在本章中显示任何真实的推文。相反，我们可以使用桑德的手标数据，其中包含推文 id 和他们的手标情绪。我们将使用推特的应用编程接口逐个获取相应的推文。为了不让你太厌烦，只需执行相应 Jupyter 笔记本的第一部分，这将启动下载过程。为了很好地使用推特的服务器，下载 5000 多条推特的所有数据需要相当长的时间，这意味着立即开始是个好主意。

数据自带四个情感标签，由load_sanders_data()返回:

>>> X_orig, Y_orig = load_sanders_data()
>>> classes = np.unique(Y_orig)
>>> for c in classes: print("#%s: %i" % (c, sum(Y_orig == c)))
#irrelevant: 437
#negative: 448
#neutral: 1801
#positive: 391

在load_sanders_data()里面，我们是把无关的、中性的标签一起当作中性的对待，把所有非英文的推文都去掉，结果推文 3077 条。

如果你在这里得到不同的计数，这是因为，在此期间，推文可能已经被删除或设置为私人的。在这种情况下，您可能还会得到与接下来几节中显示的数字和图表略有不同的数字和图表。

引入朴素贝叶斯分类器

天真的贝叶斯可能是最优雅的机器学习算法之一，具有实用价值。尽管它的名字，但当你看它的分类表现时，它并不那么幼稚。事实证明，它对不相关的特性相当稳健，但它会善意地忽略这些特性。它学得很快，预测也一样。它不需要大量存储空间。那么，为什么叫幼稚呢？

添加天真是为了解释天真贝叶斯优化工作所需的一个假设。假设特征是不相关的。然而，现实世界的应用程序很少会出现这种情况。然而，它在实践中仍然返回非常好的准确性，即使当独立性假设不成立时。

了解贝叶斯定理

本质上，朴素贝叶斯分类只不过是跟踪哪个特征为哪个类别提供证据。特征的设计方式决定了用来学习的模型。所谓伯努利模型只关心布尔特征；一个词在一条推文中出现一次还是多次并不重要。相比之下，多项式模型使用字数作为特征。为了简单起见，我们将使用伯努利模型来解释如何使用朴素贝叶斯进行情感分析。然后我们将使用多项式模型来建立和调整我们的真实世界分类器。

让我们假设以下变量的含义，我们将用来解释朴素贝叶斯:

| 变量 | 表示 | | c | 这是推文的类别(正面或负面——对于这个解释，我们忽略中性标签) | | F<sub>1</sub> | “棒极了”这个词在推文中至少出现过一次 | | F<sub>2</sub> | 疯狂这个词在推文中至少出现过一次 |

在训练中，我们学习了朴素贝叶斯模型，这是当我们已经知道特征和时类的概率。这个概率写成。

由于我们不能直接估计，我们应用了一个技巧，这个技巧被贝叶斯发现:

如果我们将替换为“可怕”和“疯狂”的概率，并将视为我们的类，我们会得到一个关系，帮助我们稍后检索属于指定类的数据实例的概率:

这允许我们通过其他概率来表达:

我们也可以这样描述:

在先和证据很容易确定:

  • 是不知道数据的类的先验概率。我们可以通过简单地计算属于该特定类的所有训练数据实例的分数来估计这个数量。
  • 是特征和的证据或概率

棘手的部分是可能性的计算。如果我们知道数据实例的类是，那么这个值描述了看到和特征值的可能性。要估计这一点，我们需要做一些思考。

天真

从概率论中，我们还知道以下关系:

然而，这本身并没有多大帮助，因为我们用另一个难题(估计)来处理一个难题(估计)。

然而，如果我们天真地假设和相互独立，简化为，我们可以这样写:

把所有的东西放在一起，我们得到了以下非常容易管理的公式:

有趣的是，尽管在我们有心情的时候简单地调整我们的假设在理论上是不正确的，但在这种情况下，它在现实世界的应用中被证明是非常有效的。

用朴素贝叶斯分类

给定一条新的推文，剩下的唯一部分就是计算概率:

然后选择类作为概率较高的。

至于两个类，分母是一样的，我们可以简单的忽略它，不改变胜者类。

然而，请注意，我们不再计算任何真实的概率。相反，我们正在估计哪一类更有可能给出证据。这是朴素贝叶斯如此稳健的另一个原因:它对真实概率不感兴趣，只对信息感兴趣，哪个类更有可能。简而言之，我们可以这样写:

简单来说就是我们在计算 ( pos 和 neg 的所有类在 argmax 之后的部分，并返回得到最高值的类。

但是，对于下面的例子，让我们坚持真实概率，做一些计算，看看朴素贝叶斯是如何工作的。为了简单起见，我们将假设推特只允许前面提到的两个词，可怕和疯狂，并且我们已经手动分类了一些推文:

推文级
可怕的正面推文
可怕的正面推文
可怕的疯狂正面推文
疯狂的正面推文
疯狂的负面推文
疯狂的负面推文

在这个例子中，我们在正面和负面的推文中都有疯狂的推文，以模仿你在现实世界中经常会发现的一些歧义(例如，对足球疯狂对疯狂的白痴)。

在这种情况下，我们总共有六条推文，其中四条是积极的，两条是消极的，这导致了以下优先事项:

这意味着，在对推文本身一无所知的情况下，假设推文是正面的是明智的。

我们仍然缺少和的计算，这是两个特征和的概率，在类中有条件。

这是通过我们已经看到具体特征的推文数量除以已经用类标记的推文数量来计算的。假设我们想知道在推文中看到 awesome 发生的概率，知道它的类是正的，我们将有:

这是因为在四条积极的推文中，有三条包含了“棒极了”这个词。显然，在正面推文中没有令人敬畏的可能性是相反的:

同样，对于其余部分(省略一个词没有出现在推文中的情况):

为了完整起见，我们还将计算证据，以便我们可以在下面的示例推文中看到真实概率。对于和两个具体数值，我们可以计算出如下证据:

这导致以下值:

现在我们有了对新推文进行分类的所有数据。剩下的唯一工作就是解析这条推文，分析它的特点:

| 推文 | | | 类概率 | 分类 | | 可怕的 | one | Zero | | 积极的 | | 疯狂的 | Zero | one | | 否定的；消极的；负面的；负的 | | 可怕的疯狂 | one | one | | 积极的 |

目前为止，一切顺利。琐碎推文的分类似乎给推文分配了正确的标签。然而，问题仍然是，我们应该如何对待训练语料库中没有出现的单词。毕竟，有了前面的公式，新单词总会被赋予零的概率。

解释看不见的单词和其他奇怪的东西

当我们早先计算概率时，我们实际上欺骗了自己。我们不是在计算真实的概率，而只是通过分数进行粗略的近似。我们假设训练语料库会告诉我们真实概率的全部真相。它没有。一个只有六条推文的语料库显然不能给我们所有关于每一条被写过的推文的信息。例如，肯定有包含文字的推文。只是我们从未见过他们。显然，我们的近似值非常粗略，我们应该考虑到这一点。这在实践中经常用所谓的加一平滑来完成。

Add-one smoothing is sometimes also referred to as additive smoothing or Laplace smoothing. Note that Laplace smoothing has nothing to do with Laplacian smoothing, which is related to the smoothing of polygon meshes. If we do not smooth by 1 but by an adjustable parameter, alpha>0, it is called Lidstone smoothing.

这是一种非常简单的技术，可以将一个要素添加到所有要素实例中。它有一个潜在的假设，即使我们没有在整个语料库中看到一个给定的单词，我们的推文样本仍然有可能不包括这个单词。因此，通过加一平滑，我们假装我们比实际看到的更多地看到了每一个事件。这意味着我们现在不用计算，而是做。

为什么分母要加 2？因为我们有两个特点:牛逼和疯狂的发生。因为我们为每个特征加 1，我们必须确保最终结果还是一个概率。事实上，我们得到 1 作为总概率:

算术欠流核算

还有另一个路障。实际上，我们处理的概率比我们在玩具例子中处理的概率小得多。通常情况下，我们也有许多不仅仅是两个特征，它们是相互相乘的。这将很快导致 NumPy 提供的浮点精度不再满足要求:

>>> import numpy as np
>>> # tell numpy to print out more digits (default is 8)
>>> np.set_printoptions(precision=20) 
>>> 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(mpmath.org)。然而，它们的速度不足以作为 NumPy 的替代品。

幸运的是，有一个更好的方法来解决这个问题，这与我们可能还记得的学校里的一段美好关系有关(也称为 logsum-trick):

如果我们将此公式应用于我们的案例，我们会得到以下结果:

由于概率在 0 和 1 之间的区间内，概率的对数在-和 0 之间的区间内。不要为此烦恼。更高的数字仍然是正确类别的一个更强有力的指标——只是它们现在是负数:

但是有一个警告:我们实际上没有公式命名者的日志(分数的顶部)。我们只有概率的乘积。幸运的是，在我们的例子中，我们对概率的实际值不感兴趣。我们只是想知道哪个类的后验概率最高。我们很幸运，因为如果我们发现以下情况:

那么我们也将永远拥有:

。

快速看一下前面的图表，可以发现曲线是严格单调递增的，也就是说，当我们从左向右时，它总是向上的。让我们把这个放在前面提到的公式中:

这将最终检索出两个特性的公式，这两个特性将为我们提供最好的类，也是我们将在实践中看到的真实数据:

当然，如果只有两个特性，我们将不会非常成功，所以让我们重写它，以允许任意数量的特性:

我们在这里，准备使用 scikit-learn 工具包中的第一个分类器。

如前所述，我们刚刚学习了朴素贝叶斯的伯努利模型。除了布尔特征，我们还可以使用单词出现的次数，也称为多项式模型。由于这提供了更多的信息，并且通常还会带来更好的性能，因此我们将把它用于我们的真实数据。但是，请注意，基础公式会有一些变化。然而，不用担心，因为朴素贝叶斯的工作原理仍然是一样的。

创建我们的第一个分类器并对其进行调整

天真的贝叶斯分类器存在于sklearn.naive_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_orig=="positive", Y_orig =="negative")

>>> # now use that index to filter the data and the labels
>>> X = X_orig [pos_neg_idx]
>>> Y = Y_orig [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(params=None):
    tfidf_ngrams = TfidfVectorizer(ngram_range=(1, 3), 
                                   analyzer="word", binary=False)
    clf = MultinomialNB()
    pipeline = Pipeline([('tfidf', tfidf_ngrams), ('clf', clf)])
    if params:
        pipeline.set_params(**params)
    return pipeline

由create_ngram_model()返回的Pipeline实例现在可以用于拟合和预测，就像我们有一个正常的分类器一样。稍后，我们将传递一个参数字典作为params，这将帮助我们创建自定义管道。

由于我们没有那么多数据，我们应该进行交叉验证。但是，这次我们不会使用KFold，它将数据划分为连续的折叠；相反，我们将使用ShuffleSplit。它为我们打乱了数据，但并不阻止同一数据实例被多次折叠。然后，对于每个折叠，我们跟踪精度-召回曲线下的区域，以确保准确性。

为了让我们的实验保持敏捷，让我们把所有的东西都包装在一个train_model()函数中，该函数将一个函数作为创建分类器的参数:

from sklearn.metrics import precision_recall_curve, auc
from sklearn. model_selection import ShuffleSplit

def train_model(clf_factory, X, Y):
    # setting random_state to get deterministic behavior
    cv = ShuffleSplit(n_splits=10, test_size=0.3, 
                       random_state=0)

    scores = []
    pr_scores = []

    for train, test in cv.split(X, Y):
        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.mean(pr_scores))
        print("Mean acc=%.3ftMean P/R AUC=%.3f" % summary)

把所有东西放在一起，我们可以训练我们的第一个模型:

 >>> X_orig, Y_orig = load_sanders_data()
 >>> pos_neg_idx = np.logical_or(Y_orig =="positive", Y_orig =="negative")
 >>> X = X_orig[pos_neg_idx]
 >>> Y = Y_orig [pos_neg_idx]
 >>> Y = Y_orig =="positive"
 >>> train_model(create_ngram_model, X, Y)
 Mean acc=0.777    Mean P/R AUC=0.885

我们首次尝试在向量化的 TF-IDF 三元组特征上使用朴素贝叶斯，得到了 77.7%的准确率和 88.5%的平均市盈率 AUC。查看中位数的市盈率图表(表现与平均值最相似的列车/测试分割)，它显示了比我们在上一章中看到的图表更令人鼓舞的行为。请注意，0.90 的曲线的 AUC 与 0.885 的平均 P/R 略有不同，因为该曲线取自训练运行的中位数，而平均 P/R AUC 是所有 AUC 分数的平均值。同样的原理也适用于后续图像:

首先，结果相当令人鼓舞。当我们意识到 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
    return Y.astype(int)

请注意，我们现在谈论的是两种不同的积极因素。推文的情绪可以是积极的，这要与训练数据的类别区分开来。例如，如果我们想知道如何区分有情绪的推文和中性的推文，我们可以做:

>>> X = X_orig
>>> Y = tweak_labels(Y_orig, ["positive", "negative"])

在Y中，我们现在有1(积极类)用于所有积极或消极的推文，还有0(消极类)用于中性和不相关的推文:

>>> train_model(create_ngram_model, X, Y)
Mean acc=0.734    Mean P/R AUC=0.661

看看下面的情节:

不出所料，市盈率 AUC 大幅下降，目前仅为 66%。准确性仍然很高，但这只是因为我们有一个高度不平衡的数据集。在总共 3077 条推文中，只有 839 条是正面或负面的，约占 27%。这意味着，如果我们创建一个分类器，总是将推文分类为不包含任何情绪，我们将已经有 73%的准确率。如果训练和测试数据不平衡，这是另一个总是看精度和回忆的原因。

那么，天真的贝叶斯分类器将如何对积极的推文和其他推文进行分类，以及对消极的推文和其他推文进行分类？一个字:差:

== Pos vs. rest == 
Mean acc=0.87  Mean P/R AUC=0.305 
== Neg vs. rest == 
Mean acc=0.852 Mean P/R AUC=0.49

如果你问我的话，我觉得很难用。查看以下图表中的市盈率曲线，我们也将发现没有可用的精度/召回权衡，正如我们在上一章中所做的那样:

调整分类器的参数

当然，我们还没有充分探索当前的设置，应该进行更多的调查。大致有两个区域可以玩旋钮:TfidfVectorizer和MultinomialNB。由于我们没有真正的直觉，我们应该探索哪个领域，让我们尝试扫描超参数。

我们将首先看到TfidfVectorizer参数:

  • 使用不同的 ngrams 设置:
    • unigrams (1，1)
    • 单图和双图(1，2)
    • 单项式、二元式和三元式(1，3)
  • 玩min_df : 1或2
  • 使用use_idf和smooth_idf : False或True探索以色列国防军在 TF-以色列国防军中的影响
  • 通过将stop_words设置为english 或None 是否删除停止词
  • 是否使用字数的对数(sublinear_tf)
  • 通过将binary设置为True或False，是跟踪字数还是简单地跟踪单词是否出现

现在我们将看到MultinomialNB分类器:

  • 通过设置alpha使用哪种平滑方法
  • 加一或拉普拉斯平滑:1
  • 丽德斯通平滑:0.01, 0.05, 0.1, or 0.5
  • 无平滑:0

一个简单的方法是为所有这些合理的探索值训练一个分类器，同时保持其他参数不变并检查分类器的结果。由于我们不知道这些参数是否会相互影响，要做好这一点，我们需要为所有参数值的每个可能组合训练一个分类器。显然，这对我们来说太繁琐了。

因为这种参数探索在机器学习任务中经常出现，scikit-learn 有一个专门的类，叫做GridSearchCV。它需要一个估计器(带有类分类器接口的实例)，在我们的例子中是Pipeline实例，以及一个带有潜在值的参数字典。

GridSearchCV期望字典的关键字遵循某种格式，以便能够设置正确估计器的参数。格式如下:

<estimator>__<subestimator>__...__<param_name>

例如，如果我们想要为TfidfVectorizer(在Pipeline描述中被命名为tfidf的ngram_range参数指定想要探索的值，我们将不得不说:

param_grid={"tfidf__ngram_range"=[(1, 1), (1, 2), (1, 3)]}

这将告诉GridSearchCV尝试将三元组的单值作为TfidfVectorizer的ngram_range参数的参数值。

然后，它用所有可能的参数值组合训练估计量。我们确保它使用ShuffleSplit对训练数据的随机样本进行训练，这将生成随机训练/测试拆分的迭代器。最后，它以成员变量best_estimator_的形式提供了最佳估计量。

当我们想要将返回的最佳分类器与我们当前的最佳分类器进行比较时，我们需要以同样的方式对其进行评估。因此，我们可以使用cv参数传递ShuffleSplit实例(因此，GridSearchCV中的CV)。

最后缺少的部分是定义GridSearchCV应该如何确定最佳估计量。这可以通过使用make_scorer辅助函数向scoring参数提供所需的评分函数来实现。我们可以自己写一个，也可以从sklearn.metrics包里挑一个。我们当然不应该拿metric.accuracy来说事，因为我们的阶级不平衡(我们包含情感的推文比中立的少得多)。相反，我们希望对这两个类、带有情绪的推文和没有正面或负面意见的推文都有很好的精度和召回率。结合了精确度和召回率的一个度量是 F-measure ，它被实现为metrics.f1_score:

将所有内容放在一起后，我们得到以下代码:

from sklearn. model_selection import GridSearchCV
from sklearn.metrics import make_scorer, f1_score

def grid_search_model(clf_factory, X, Y):
    cv = ShuffleSplit(n_splits=10, test_size=0.3, random_state=0)

    param_grid = dict(tfidf__ngram_range=[(1, 1), (1, 2), (1, 3)],
            tfidf__min_df=[1, 2],
            tfidf__stop_words=[None, "english"],
            tfidf__smooth_idf=[False, True],
            tfidf__use_idf=[False, True],
            tfidf__sublinear_tf=[False, True],
            tfidf__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,
            scoring=make_scorer(f1_score), verbose=10)
    grid_search.fit(X, Y) 

    return grid_search.best_estimator_

执行此操作时，我们必须保持耐心:

print("== Pos/neg vs. irrelevant/neutral ==")
X = X_orig
Y = tweak_labels(Y_orig, ["positive", "negative"])
clf = grid_search_model(create_ngram_model, X, Y)
print(clf)

由于我们刚刚请求了一个参数，扫描参数组合，每个参数组合被训练 10 次:

... waiting some 20 minutes  ...
Pipeline(memory=None,
    steps=[('tfidf', TfidfVectorizer(analyzer='word', binary=True, 
       decode_error='strict',
       dtype=<class 'numpy.int64'>, encoding='utf-8', input='content',
       lowercase=True, max_df=1.0, max_features=None, min_df=1,
    ngram_range=(1, 2), norm='l2', preprocessor=None, 
smooth_idf=False, vocabulary=None)), 
('clf', MultinomialNB(alpha=0.01, class_prior=None, fit_prior=True))])

为了能够将这些数字与我们之前的方法进行比较，我们将创建一个best_params字典，然后将其传递给分类器工厂，然后运行与之前相同的代码，对 10 倍 CV 拆分进行训练并输出平均分数:

best_params = dict(all__tfidf__ngram_range=(1, 2),
                   all__tfidf__min_df=1,
                   all__tfidf__stop_words=None,
                   all__tfidf__smooth_idf=False,
                   all__tfidf__use_idf=False,
                   all__tfidf__sublinear_tf=True,
                   all__tfidf__binary=False,
                   clf__alpha=0.01,
                )
print("== Pos/neg vs. irrelevant/neutral ==")
X = X_orig
Y = tweak_labels(Y_orig, ["positive", "negative"])
train_model(lambda: create_ngram_model(best_params), X, Y)

结果如下:

== Pos/neg vs. irrelevant/neutral ==
Mean acc=0.791    Mean P/R AUC=0.681

最佳估计器确实将市盈率 AUC 从 65.8%提高到 68.1%，其设置如前面的代码所示。

此外，如果我们用我们刚刚发现的参数配置矢量器和分类器，那么正面推文对其余推文和负面推文对其余推文的破坏性结果会得到改善。只有阳性和阴性分类显示出稍差的性能:

看看下面的情节:

事实上，市盈率曲线看起来要好得多(注意，曲线来自折叠分类器的中间，因此，AUC 值略有差异)。然而，我们可能仍然不会使用这些分类器。是时候做些完全不同的事情了...

清理推文

新的约束导致新的形式。在这方面，推特也不例外。因为文本必须适合 280 个字符，人们自然会开发新的语言快捷方式，用更少的字符说出相同的内容。到目前为止，我们忽略了所有不同的表情符号和缩写。让我们看看考虑到这一点，我们能提高多少。为了这个努力，我们必须提供我们自己的preprocessor()到TfidfVectorizer。

首先，我们在字典中定义一系列常见的表情符号及其替换。虽然我们可以找到更明显的替代词，但我们会用明显的肯定或否定的词来帮助分类器:

    emo_repl = {
     # positive emoticons
     "&lt;3": " good ",
     ":d": " good ", # :D in lower case
     ":dd": " good ", # :DD in lower case
     "8)": " good ",
     ":-)": " good ",
     ":)": " good ",
     ";)": " good ",
     "(-:": " good ",
     "(:": " good ", 
     # negative emoticons:
     ":/": " bad ",
     ":&gt;": " 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"brb": "are",
     r"bub": "you",
     r"bhahab": "ha",
     r"bhahahab": "ha",
     r"bdon'tb": "do not",
     r"bdoesn'tb": "does not",
     r"bdidn'tb": "did not",
     r"bhasn'tb": "has not",
     r"bhaven'tb": "have not",
     r"bhadn'tb": "had not",
     r"bwon'tb": "will not",
     r"bwouldn'tb": "would not",
     r"bcan'tb": "can not",
     r"bcannotb": "can not",
 }

def create_ngram_model_emoji(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")
    clf = MultinomialNB()
    pipeline = Pipeline([('tfidf', tfidf_ngrams), ('clf', clf)])
    if params:
        pipeline.set_params(**params)

    return pipeline

请记住，我们创建了一个参数字典，我们发现最好使用GridSearchCV。我们将把它传递给create_ngram_model_emoji，这样我们的新模型将在我们已经发现的基础上进行改进。由于train_model需要一个分类器工厂，因为它会被一次又一次地实例化，我们使用 Python 的lambda创建工厂:

print("== Pos/neg vs. irrelevant/neutral ==")
X = X_orig
Y = tweak_labels(Y_orig, ["positive", "negative"])
train_model(lambda: create_ngram_model_emoji(best_params), X, Y)

当然，这里可以使用更多的缩写。但是已经有了这个有限的集合，我们得到了大约半个点的积极与消极的改善，以及情绪与无情绪的改善。复制前面的表格并为新方法填写数字，我们看到它为每个分类器都批准了一点:

考虑到单词类型

到目前为止，我们的希望是，简单地使用相互独立的单词和单词包的方法就足够了。然而，仅仅从我们的直觉来看，我们知道中性推文可能包含更高比例的名词，而正面或负面推文更丰富多彩，需要更多的形容词和动词。如果我们也使用推文的语言信息呢？如果我们能找出一条推文中有多少单词是名词、动词、形容词等等，分类器可能也会考虑到这一点。

确定单词类型

这就是词性标注，或者说词性标注的意义所在。POS tagger 分析一个句子，用词类标记每个单词，例如，book 这个词是名词(这是一本好书)还是动词(能不能请你订航班？).

你可能已经猜到 NLTK 也会在这个领域发挥作用。事实上，它很容易打包成各种解析器和标记器。我们将使用的词性标注器nltk.pos_tag()实际上是一个成熟的分类器，它是使用佩恩树库项目中手动标注的句子训练的。它将单词标记列表作为输入，并输出元组列表，其中每个元素包含原始句子的部分及其词性标记:

>>> 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'), ('?', '.')]

POS 标签缩写取自佩恩树库(改编自www.anc.org/OANC/penn.h…):

POS 标签描述例
CC并列连词或者
CD基数第二
DT限定词这
EX那里存在有有
FW外来词幼儿园
IN介词/从属连词在，在，像
JJ形容词凉爽的
JJR形容词，比较的冷却器
JJS形容词，最高级最凉快的
LS列表制作人1)
MD情态的可能，威尔
NN名词，单数或复数书
NNS名词复数书
NNP专有名词，单数肖恩
NNPS专有名词，复数海盗
PDT前限定词两个男孩
POS所有格结尾朋友的
PRP人称代名词我，他，它
PRP$所有格代名词我的，他的
RB副词然而，通常，自然地，在这里，很好
RBR副词，比较的较好的
RBS副词，最高级最好的
RP颗粒放弃放弃
TO到去去，去他
UH感叹词嗯嗯嗯
VB动词，基本形式拿
VBD动词，过去式拿
VBG动词、动名词/现在分词拿
VBN动词，过去分词拿
VBP动词，唱。存在，非 3d拿
VBZ动词，第三人称歌唱。礼物接受
WDTwh-确定哪个
WPwh 代词谁，什么
WP$所有格 wh 代词谁的
WRBwh-abverb何时何地

有了这些标签，从pos_tag()的输出中过滤出想要的标签是非常容易的。我们只需计算所有标签以NN开头的名词、VB开头的动词、 开头的形容词和RB开头的副词。

使用 SentiWordNet 成功作弊

虽然语言信息，如前一节所述，很可能会对我们有所帮助，但我们可以用一些更好的东西来收获它:SentiWordNet(http://SentiWordNet . isti . CNR . it)。简单来说就是一个 13 MB 的文件，我们要从网站下载，解压，放入 Jupyter 笔记本的数据目录；它给大多数英语单词赋予了一个正值和负值。这意味着对于每个同义词集，它都记录了积极和消极的情感值。一些例子如下:

POSID姿势得分否定得分同步组术语描述
a00311354Zero point two fiveZero point one two five勤奋好学#1以小心和努力为特点的；努力修理电视机
a00311663ZeroZero point five粗心的#1以缺乏关注、考虑、深谋远虑或彻底为特点的；不小心...
n03563710ZeroZero植入物#1永久放置在组织中的假体
v00362128ZeroZero扭结#2 曲线#5 卷曲#1形成卷曲、曲线或扭结；雪茄烟雾在天花板上袅袅上升

有了 POS 一栏的信息，我们就能区分名词书和动词书了。PosScore和NegScore一起会帮助我们确定这个词的中性，就是 1-poss core-neggle。SynsetTerms列出集合中所有同义词。我们可以放心地忽略ID和Description列。

synset 术语附加了一个数字，因为有些术语在不同的 synset 中出现多次。例如，innovate 传达了两种不同的意思，这也导致了不同的分数:

| POS | ID | 姿势得分 | 否定得分 | 同步组术语 | 描述 | | v | 01636859 | Zero point three seven five | Zero | 幻想#2 幻想 | 在头脑中描绘；他幻想着理想的妻子 | | v | 01637368 | Zero | Zero point one two five | 幻想#1 幻想#1 幻想 | 沉溺于幻想；当他说他计划创办自己的公司时，他是在幻想 |

为了找出应该采用哪种 synsets，我们需要真正理解推文的含义，这超出了本章的范围。专注于这一挑战的研究领域被称为词义消歧。 对于我们的任务，我们采取简单的路线，简单地对找到一个术语的所有合成集进行平均。对于幻想来说，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

然而，在组合特征器上的训练和测试有点令人失望。相对于其余部分，我们在积极方面提高了 1%，但在其他方面却有所下降:

有了这些结果，如果我们不能在市盈率 AUC 中获得一个重要的职位，我们可能不想为昂贵得多的 SentiWord 方法付出代价。相反，我们可能会选择 GridCV + Cleaning 方法，首先使用确定推文是否包含情感的分类器(pos/neg 对无关/中性)，然后在确实包含情感的情况下，使用肯定对否定的分类器来确定实际的情感。

摘要

恭喜你坚持到最后！我们一起学习了朴素贝叶斯是如何工作的，以及为什么它一点也不幼稚。特别是对于训练集，我们没有足够的数据来学习类概率空间中的所有小生境，朴素贝叶斯在推广方面做得很好。我们学会了如何将它应用到推文中，清理粗糙的推文文本有很大帮助。最后，我们意识到一点点欺骗是可以的(只有在我们完成了应尽的工作之后)。然而，由于我们意识到昂贵得多的分类器并没有奖励我们一个改进得多的分类器，我们回到了便宜的分类器。

在第 10 章、主题建模中，我们将学习如何使用潜在狄利克雷分配(也称为主题建模)从文档中提取主题。这将有助于我们通过分析所涵盖主题的相似程度来比较文档。

十、主题建模

在第 6 章、聚类-查找相关帖子中，我们使用聚类对文本文档进行了分组。这是一个非常有用的工具，但并不总是最好的。聚类结果是每个文本恰好属于一个聚类。这本书是关于机器学习和 Python 的。它应该与其他 Python 相关的作品还是与机器相关的作品分组？在实体书店，我们需要选择一个单一的地方来存放这本书。然而，在一家互联网商店，答案是这本书是关于机器学习和 Python 的，这本书应该列在这两个部分。当然，这并不意味着这本书会在所有章节中列出。我们不会把这本书和其他烘焙书籍放在一起。

在本章中，我们将学习不将文档聚集成完全独立的组，而是允许每个文档引用几个主题的方法。这些主题将从文本文档集合中自动识别。这些文档可以是整本书或较短的文本，如博客文章、新闻故事或电子邮件。

我们还希望能够推断出这样一个事实，即这些文件可能有对它们至关重要的主题，而只是顺便提及其他主题。这本书经常提到绘图，但它不像机器学习那样是一个中心话题。这意味着文档中的主题是它们的核心，而其他主题则是次要的。处理这些问题的机器学习子领域称为主题建模，是本章的主题。特别是，我们将了解以下内容:

  • 有哪些主题模型，特别是关于潜在狄利克雷分配 ( LDA )
  • 如何使用gensim包构建主题模型
  • 主题模型如何作为不同应用程序的中间表示有用
  • 我们如何建立整个英语维基百科的主题模型

潜在狄利克雷分配

不幸的是，机器学习中有两种方法，首字母是 LDA:潜在狄利克雷分配，这是一种主题建模方法，以及线性判别分析，这是一种分类方法。它们完全不相关，除了首字母 LDA 可以指任何一个。在某些情况下，这可能会令人困惑。scikit-learn 工具有一个子模块sklearn.lda，实现线性判别分析。目前，scikit-learn 没有实现潜在的 Dirichlet 分配。

第一个话题model我们要看的是潜狄利克雷分配。LDA 背后的数学思想相当复杂，这里就不赘述了。

对于那些感兴趣并且足够冒险的人，维基百科提供了这些算法背后的所有等式:en.wikipedia.org/wiki/Latent…。

但是，我们可以从高层次直观地理解 LDA 背后的思想。LDA 属于一类被称为生成模型的模型，因为它们有一种寓言解释了数据是如何生成的。这个生成性的故事是对现实的简化，当然是为了让机器学习更容易。在 LDA 寓言中，我们首先通过给单词分配概率权重来创建主题。每个主题将为不同的单词分配不同的权重。例如，一个 Python 主题会给单词变量分配高概率，给喝醉的单词分配低概率。当我们希望生成一个新文档时，我们首先选择它将使用的主题，然后从这些主题中混合单词。

例如，假设我们只有三个书籍讨论的主题:

  • 机器学习
  • 计算机编程语言
  • 烘烤

对于每个主题，我们都有一个与之相关的单词列表。这本书将前两个主题混合在一起，也许各占 50%。混合物不需要相等；也可以是 70/30 分割。当我们生成实际的文本时，我们一个字一个字地生成；首先我们决定这个词来自哪个话题。这是基于主题权重的随机决定。一旦选择了一个主题，我们就从该主题的单词列表中生成一个单词。准确地说，我们用题目给出的概率选择一个英语单词。同一个单词可以由多个主题生成。例如，重量是机器学习和烘焙中的常用词(尽管含义不同)。

在这个模型中，单词的顺序并不重要。这是一个包字模型，我们在上一章已经看到了。这是对语言的粗略简化，但它通常足够好，因为仅仅知道文档中使用了哪些单词及其频率就足以做出机器学习决策。

在现实世界中，我们不知道主题是什么。我们的任务是收集一些文本，并对这个寓言进行逆向工程，以便发现有哪些主题，同时找出每个文档使用哪些主题。

构建主题模型

不幸的是，scikit-learn 没有实现潜在的 Dirichlet 分配。因此，我们将使用 Python 中的gensim包。Gensim 是由英国的机器学习研究员和顾问 radim řehůřek 开发的。

作为输入数据，我们将使用来自美联社 ( 美联社)的新闻报道集合。这是一个用于文本建模研究的标准数据集，在一些关于主题模型的初始工作中使用过。下载数据后，我们可以通过运行以下代码来加载它:

import gensim 
from gensim import corpora, models 
corpus = corpora.BleiCorpus('./data/ap/ap.dat', './data/ap/vocab.txt')

corpus变量保存所有文本文档，并以易于处理的格式加载它们。我们现在可以建立一个主题model，使用这个对象作为输入:

model = models.ldamodel.LdaModel( 
              corpus, 
              num_topics=100, 
              id2word=corpus.id2word)

该单个构造函数调用将统计推断corpus中存在哪些主题。我们可以通过多种方式探索最终的模型。我们可以看到使用model[doc]语法的topics文档列表，如下例所示:

doc = corpus.docbyoffset(0) 
topics = model[doc] 
print(topics) 
[(3, 0.023607255776894751), 
 (13, 0.11679936618551275), 
 (19, 0.075935855202707139), 
.... 
 (92, 0.10781541687001292)]

结果几乎肯定会在我们的电脑上看起来不同！学习算法使用了一些随机数，每次在同一个输入数据上学习一个新的题目model，结果都不一样。重要的是，如果您的数据表现良好，模型的一些定性属性将在不同的运行中保持稳定。例如，如果您正在使用主题来比较文档，就像我们在这里所做的那样，那么相似性应该是健壮的，并且只有轻微的变化。另一方面，不同主题的顺序也会完全不同。

结果的格式是一个配对列表:(topic_index, topic_weight)。我们可以看到每个文档只使用了几个主题(在前面的例子中，主题0、1、2没有权重；这些话题的权重是0。主题model是一个稀疏模型，因为虽然有很多可能的主题，但是对于每个文档，只使用了其中的几个。严格来说，这不是真的，因为在 LDA 模型中，所有主题都有非零概率，但其中一些主题的概率非常小，我们可以将其四舍五入为零，作为一个很好的近似值。

我们可以通过绘制每个文档涉及的主题数量的直方图来进一步探讨这一点:

num_topics_used = [len(model[doc]) for doc in corpus] 
fig,ax = plt.subplots() 
ax.hist(num_topics_used)

你会得到如下图:

Sparsity means that while you may have large matrices and vectors, in principle, most of the values are zero (or so small that we can round them to zero as a good approximation). Therefore, only a few things are relevant at any given time. Often problems that seem too big to solve are actually feasible because the data is sparse. For example, even though any web page can link to any other web page, the graph of links is actually very sparse as each web page will link to a very tiny fraction of all other web pages.

在上图中，我们可以或多或少地看到大多数文档涉及大约 10 个主题。

这在很大程度上是由于所用参数的值，即alpha参数。alpha 的确切含义有点抽象，但是 alpha 值越大，每个文档的主题就越多。

Alpha 需要是大于零的值，但通常设置为较小的值，通常小于 1。alpha的值越小，每个文档期望讨论的主题就越少。默认情况下，gensim 会将alpha设置为1/num_topics，但是您可以通过在LdaModel构造函数中将它作为参数进行显式设置，如下所示:

model = models.ldamodel.LdaModel( 
              corpus, 
              num_topics=100, 
              id2word=corpus.id2word, 
              alpha=1)

在这种情况下，这是一个比默认值更大的 alpha，这应该会导致每个文档有更多的主题。正如我们在下面给出的组合直方图中所看到的，gensim 的行为符合我们的预期，并为每个文档分配了更多的主题:

现在，我们可以在前面的柱状图中看到，许多文档涉及到 20 到 25 不同的主题。如果您将该值设置得更低，您将观察到相反的情况(从在线存储库中下载代码将允许您使用这些值)。

这些话题是什么？从技术上讲，正如我们前面讨论的，它们是单词的多项式分布，这意味着它们为词汇表中的每个单词分配了一个概率。概率较高的单词比概率较低的单词更容易与该主题相关联:

我们的大脑不太擅长用概率分布进行推理，但我们可以很容易地理解一系列单词。因此，通常用权重最高的单词列表来总结主题。

在下表中，我们显示了前十个主题:

| 题号 | 话题 | | one | 着装军事苏联总统新州卡卢奇上尉州领导人立场政府 | | Two | 科赫赞比亚卢萨卡一党橙色科赫一党政府市长新政治 | | three | 人权土耳其侵犯皇家汤普森威胁新州写道花园总统 | | four | 比尔雇员实验莱文税收联邦措施立法参议院主席告密者赞助 | | five | 俄亥俄州 7 月干旱耶稣灾难百分之哈特福德密西西比州北部河谷作物弗吉尼亚州 | | six | 美国百分之十亿年总统世界各国人民i布什新闻 | | seven | b休斯宣誓书美国盎司平方英尺护理延迟指控不切实际的布什 | | eight | 尤特·杜卡基斯布什公约农业补贴乌拉圭百分比秘书长i告诉 | | nine | 克什米尔政府人民斯利那加印度甩市两查谟克什米尔集团莫斯利巴基斯坦 | | Ten | 工人越南爱尔兰工资移民百分比议价最后一个岛警察赫顿i |

虽然乍一看令人望而生畏，但当通读单词列表时，我们可以清楚地看到，主题不仅仅是随机的单词，相反，这些是逻辑组。我们也可以看到，这些话题指的是比较老的新闻条目，从苏联还存在，戈尔巴乔夫还是总书记的时候开始。我们还可以将主题表示为单词云，使更有可能的单词变大。例如，这是一个关于警察的主题的可视化:

我们还可以看到，有些词可能应该删除，因为它们信息不多；这些都是空话。在构建主题建模时，过滤掉停止词可能会很有用，否则你可能会得到一个完全由停止词组成的主题。我们可能还希望将文本预处理为词干，以便规范化复数和动词形式。这个过程在第 6 章、聚类-查找相关帖子中有介绍，具体可以参考。如果你感兴趣，可以从本书的配套网站下载代码，尝试所有这些变体来绘制不同的图片。

按主题比较文档

主题本身对于构建上一个截图中显示的那种带有单词的小短文非常有用。这些可视化可以用于浏览大量文档。例如，网站可以将不同的主题显示为不同的单词云，允许用户点击查看文档。事实上，它们只是以这种方式被用来分析大量的文档集合。

然而，话题往往只是达到另一个目的的中间工具。现在，我们已经估计了每个文档中有多少来自每个主题，我们可以在主题空间中比较这些文档。这仅仅意味着，如果两个文档谈论相同的主题，我们就说它们是相似的，而不是逐字比较。

这可能非常强大，因为两个共享很少单词的文本文档实际上可能指的是同一个主题！他们可能只是使用不同的结构来引用它(例如，一个文档可能会提到英国，而另一个文档会使用缩写 UK)。

Topic models are good on their own to build visualizations and explore data. They are also very useful as an intermediate step in many other tasks.

在这一点上，我们可以通过使用主题来定义相似性，来重复寻找与输入查询最相似的帖子的任务。之前我们通过直接比较两个文档的单词向量来比较它们，现在我们可以通过比较它们的主题向量来比较两个文档。

为此，我们将把文档投影到主题空间。也就是说，我们希望有一个总结文档的主题向量。这是 第五章降维中讨论的类型的降维的另一个例子。在这里，我们向您展示主题模型是如何准确地用于这一目的的；一旦为每个文档计算了主题，我们就可以对主题向量执行操作，而忘记原始单词。如果主题有意义，它们可能比生字更有信息。此外，这可能带来计算上的优势，因为比较主题权重向量比输入词汇表(包含数千个术语)大的向量要快得多。

使用gensim，我们研究了如何计算语料库中所有文档对应的主题。我们现在将为所有文档计算这些，并将其存储在 NumPy 数组中，并计算所有成对距离:

from gensim import matutils 
topics = matutils.corpus2dense(model[corpus], num_terms=model.num_topics)

现在，topics是一个主题矩阵。我们可以使用 SciPy 中的pdist函数来计算所有成对距离。有几个可用的距离。默认情况下，它使用欧几里德距离。也就是说，通过一次函数调用，我们计算出sum((topics[ti] - topics[tj])**2)的所有值:

from scipy.spatial import distance 
distances = distance.squareform(distance.pdist(topics))

现在，我们将使用最后一个小技巧；我们将把distance矩阵的对角线元素设置为无穷大，以确保它看起来比任何其他元素都大:

 for ti in range(len(topics)): 
     distances[ti,ti] = np.inf

我们完了！对于每个文档，我们可以很容易地查找最近的元素(这是一种最近邻分类器):

def closest_to(doc_id): 
    return distances[doc_id].argmin()

This will not work if we have not set the diagonal elements to a large value: the function will always return the same element as it is the one most similar to itself (except in the weird case where two elements have exactly the same topic distribution, which is very rare unless they are exactly the same).

例如，这里有一个可能的查询文档(它是我们集合中的第二个文档):

From: geb@cs.pitt.edu (Gordon Banks) 
Subject: Re: request for information on "essential tremor" and 
 Indrol? 

In article <1q1tbnINNnfn@life.ai.mit.edu> sundar@ai.mit.edu 
 writes: 

Essential tremor is a progressive hereditary tremor that gets 
 worse 
when the patient tries to use the effected member.  All limbs, 
 vocal 
cords, and head can be involved.  Inderal is a beta-blocker and 
is usually effective in diminishing the tremor.  Alcohol and 
 mysoline 
are also effective, but alcohol is too toxic to use as a 
 treatment. 
-- 
------------------------------------------------------------------
 ---------- 
Gordon Banks  N3JXP      | "Skepticism is the chastity of the 
 intellect, and 
geb@cadre.dsl.pitt.edu   |  it is shameful to surrender it too 
 soon." 
  ----------------------------------------------------------------
 ------------

如果我们要求与closest_to(1)最相似的文档，结果我们会收到以下文档:

From: geb@cs.pitt.edu (Gordon Banks) 
Subject: Re: High Prolactin 

In article <93088.112203JER4@psuvm.psu.edu> JER4@psuvm.psu.edu 
 (John E. Rodway) writes: 
>Any comments on the use of the drug Parlodel for high prolactin 
 in the blood? 
> 

It can suppress secretion of prolactin. Is useful in cases of galactorrhea. 
Some adenomas of the pituitary secret too much. 

-- 
------------------------------------------------------------------
 ---------- 
Gordon Banks  N3JXP      | "Skepticism is the chastity of the 
 intellect, and 
geb@cadre.dsl.pitt.edu   |  it is shameful to surrender it too 
 soon."

系统返回同一作者讨论药物的帖子。

建模整个维基百科

虽然最初的 LDA 实现可能很慢，这限制了它们在小文档集合中的使用，但是现代算法可以很好地处理非常大的数据集合。根据gensim的文档，我们将为整个英语维基百科建立一个主题模型。这需要几个小时，但只需一台笔记本电脑就可以完成！有了一组机器，我们可以让它运行得更快，但我们将在后面的章节中研究这种处理环境。

首先，我们从dumps.wikimedia.org下载整个维基百科转储。这是一个大文件(目前超过 14 GB)，因此可能需要一段时间，除非您的互联网连接非常快。然后，我们将使用gensim工具对其进行索引:

python -m gensim.scripts.make_wiki enwiki-latest-pages-articles.xml.bz2 wiki_en_output

在命令行上运行前一行，不要在 Python shell 上运行！几个小时后，索引将保存在同一个目录中。此时，我们可以构建最终的主题模型。这个过程看起来和它对小 AP 数据集做的完全一样。我们首先导入几个包:

import logging, gensim

现在，我们使用标准的 Python 日志模块(使用gensim打印状态消息)设置日志。严格来说，这一步并不是必须的，但是多做一点输出来了解发生了什么是很好的:

logging.basicConfig(format='%(asctime)s : %(levelname)s : %(message)s', level=logging.INFO)

现在我们加载预处理数据:

id2word = gensim.corpora.Dictionary.load_from_text( 
              'wiki_en_output_wordids.txt') 
mm = gensim.corpora.MmCorpus('wiki_en_output_tfidf.mm')

最后，我们像前面一样构建 LDA 模型:

model = gensim.models.ldamodel.LdaModel( 
          corpus=mm, 
          id2word=id2word, 
          num_topics=100)

这又需要几个小时。您将在控制台上看到进度，这可以指示您还需要等待多长时间。

一旦完成，我们可以将主题model保存到一个文件中，这样我们就不必重做它了:

model.save('wiki_lda.pkl')

如果您退出会话并稍后返回，您可以使用以下命令再次加载模型(当然是在适当的导入之后):

model = gensim.models.ldamodel.LdaModel.load('wiki_lda.pkl')

model对象可以用来探索文档集合并构建topics矩阵，就像我们之前所做的那样。

我们可以看到，这仍然是一个稀疏的模型，即使我们比以前有更多的文档(在我们写这篇文章时超过 400 万):

lens = (topics > 0).sum(axis=1) 
print('Mean number of topics mentioned: {0:.4}'.format(np.mean(lens))) 
print('Percentage of articles mentioning <10 topics: {0:.1%}'.format( 
               np.mean(lens <= 10))) 
Mean number of topics mentioned: 6.244 
Percentage of articles mentioning <10 topics: 95.1%

因此，平均文件提到6.244话题和95.1的百分比提到10或更少的话题。

我们可以问维基百科上谈论最多的话题是什么。我们将首先计算每个主题的总权重(通过汇总所有文档的权重)，然后检索与权重最高的主题相对应的单词。这是使用以下代码执行的:

weights = topics.sum(axis=0) 
words = model.show_topic(weights.argmax(), 64)

使用与我们之前构建可视化相同的工具，我们可以看到谈论最多的主题与音乐相关，并且是一个非常连贯的主题。整整 18%的维基百科页面与这个主题部分相关(维基百科中所有单词的 5.5%被分配到这个主题)。请看下面的截图:

These plots and numbers were obtained when the book was being written. As Wikipedia keeps changing, your results will be different. We expect that the trends will be similar, but the details may vary.

或者，我们可以看看最少被谈论的话题:

words = model.show_topic(weights.argmin(), 64)

参考以下截图:

最少被提及的话题很难解释，但它的许多关键词指的是非洲的一些地方。只有 2.1%的文档涉及到它，它只代表了 0.1%的单词。

选择主题数量

到目前为止，在本章中，我们使用了固定数量的主题进行分析，即 100 个。这是一个纯粹任意的数字；我们可以使用 20 或 200 个主题。幸运的是，对于许多用途来说，这个数字并不重要。如果你打算只使用主题作为中间步骤，就像我们之前在寻找类似帖子时所做的那样，系统的最终行为很少对模型中使用的主题的确切数量非常敏感。这意味着，只要你使用足够多的主题，无论你使用 100 个主题还是 200 个主题，这个过程产生的推荐都不会有很大的不同。然而，100 通常是一个足够好的数字(而 20 对于一般的文本文档集合来说太少了)，但是如果我们有更多的文档，我们本可以使用更多的。

设置alpha值也是如此。虽然玩弄它可以改变话题，但就这种变化而言，最终结果还是很稳健的。自然，这取决于你的数据的确切性质，应该根据经验进行测试，以确保结果确实是稳定的。

Topic modeling is often an end toward a goal. In that case, it is not always very important which parameter values are used exactly. A different number of topics or values for parameters such as alpha will result in systems whose end results are almost identical in their final results.

另一方面，如果你打算直接探索这些主题，或者构建一个可视化工具来展示它们，你可能应该尝试一些值，看看哪个给你最有用或最有吸引力的结果。还有一些统计概念，如困惑，可以用来确定一系列模型中哪一个最适合数据，从而做出更明智的决定。

或者，有几种方法可以根据数据集自动确定主题的数量。一个流行的模型叫做分层狄利克雷过程 ( HDP )。同样，其背后的完整数学模型很复杂，超出了本书的范围。然而，我们可以告诉你的是，不是像在 LDA 生成方法中那样首先固定主题，而是主题本身与数据一起生成，一次一个。每当作者创建一个新文档时，他们可以选择使用已经存在的主题或者创建一个全新的主题。当创建了更多的主题时，创建新主题而不是重用现有主题的可能性就会降低，但可能性总是存在的。

这意味着我们拥有的文档越多，我们最终得到的主题也就越多。这是那些一开始不直观，但经过思考后完全有意义的陈述之一。我们正在对文档进行分组，我们拥有的示例越多，就越能分解它们。如果我们只有几篇新闻文章的例子，那么体育将是一个话题。然而，随着我们拥有的越来越多，我们开始把它分解成单独的形式:曲棍球、足球等等。因为我们有了更多的数据，我们可以开始区分细微差别，关于单个团队甚至单个球员的文章。人也是如此。在一个由许多不同背景的人组成的小组里，有几个电脑人，你可能会把他们放在一起；在稍微大一点的小组中，你将为程序员和系统管理员举行单独的聚会；而在现实世界中，我们甚至为 Python 和 Ruby 程序员举办了不同的聚会。

HDP 在gensim有售。使用它是微不足道的。为了修改我们为 LDA 编写的代码，我们只需要将对gensim.models.ldamodel.LdaModel的调用替换为对HdpModel构造函数的调用，如下所示:

hdp = gensim.models.hdpmodel.HdpModel(mm, id2word)

就是这样(只是计算时间长一点——没有免费午餐)。现在，我们可以像使用 LDA 模型一样使用这个模型，只是我们不需要指定主题的数量。

摘要

在本章中，我们讨论了主题建模。主题建模比聚类更灵活，因为这些方法允许每个文档在多个组中部分呈现。为了探索这些方法，我们使用了一个新的包，gensim。

主题建模最初是为文本开发的，在文本的情况下更容易理解，但是在第 12 章、计算机视觉中，我们将看到这些技术如何应用于图像。主题模型在现代计算机视觉研究中非常重要。事实上，与前几章不同，这一章非常接近机器学习算法研究的前沿。最初的 LDA 算法发表在 2003 年的一份科学杂志上，但是gensim用来处理维基百科的方法直到 2010 年才被开发出来，而 HDP 算法是从 2011 年开始的。研究还在继续，你可以找到很多名字很棒的变体和模式，比如印度自助餐流程(不要和中餐厅流程混淆，中餐厅流程是不同的模式)或者弹球分配(弹球是日本游戏的一种，是吃角子老虎和弹球的交叉)。

我们现在已经经历了一些主要的机器学习模式:分类、聚类和主题建模。

在第 11 章、分类三-音乐流派分类中，我们回到分类，但这次我们将探索先进的算法和方法。

十一、分类三——音乐类型分类

到目前为止，我们很幸运，每个训练数据实例都可以很容易地用特征值向量来描述。例如，在 Iris 数据集中，花由包含花的某些方面的长度和宽度值的向量表示。在基于文本的示例中，我们可以将文本转换为单词包表示，并手动创建自己的特征来捕捉文本的某些方面。

这一章会有所不同，当我们试图根据歌曲的流派来分类时。例如，我们如何表现一首三分钟长的歌曲？我们应该把它的 MP3 表现的个别位？可能不会，因为像文本一样对待它并创建一个类似于一包声音片段的东西肯定会非常复杂。不知何故，我们将不得不把一首歌转换成能充分描述它的价值矢量。

绘制我们的路线图

本章将向您展示我们如何在超出我们舒适范围的领域中提出一个合适的分类器。首先，我们将不得不使用基于声音的功能，这比我们目前使用的基于文本的功能要复杂得多。然后我们将学习如何处理比以前更多的课程。此外，我们将了解衡量分类绩效的新方法。

让我们假设一个场景，出于某种原因，我们在硬盘上发现了一堆随机命名的 MP3 文件，这些文件被假设包含音乐。我们的任务是根据音乐类型将它们分类到不同的文件夹中，如爵士乐、古典音乐、乡村音乐、流行音乐、摇滚音乐和金属音乐。

获取音乐数据

我们将使用 GTZAN 数据集，该数据集经常用于音乐流派分类任务的基准测试。它被组织成 10 个不同的流派，为了简单起见，我们将只使用其中的 6 个:古典、爵士、乡村、流行、摇滚和金属。数据集包含每个流派 100 首歌曲的前 30 秒。我们可以从opihi.cs.uvic.ca/sound/genre…下载数据集。

我们可以直接用 Python 下载和提取它，这很好，尤其是如果你使用的是 Windows，它没有附带 tarball 解压程序。

在整个 Jupyter 笔记本中，我们将利用优秀的pathlib库，该库自 3.4 版本以来就是 Python 的一部分。它允许简单的路径和文件操作:

from pathlib import Path
DATA_DIR = "data"
if not Path(DATA_DIR).exists():
    os.mkdir(DATA_DIR)
import urllib.request
genre_fn = 'http://opihi.cs.uvic.ca/sound/genres.tar.gz'
# The division operator of Path instances is overloaded to behave 
# like os.path.join(), which makes it very convenient to use.
urllib.request.urlretrieve(genre_fn, Path(DATA_DIR) / 'gen-res.tar.gz')

现在我们已经下载了它，我们使用tarfile模块提取它:

import tarfile
cwd = os.getcwd()
os.chdir(DATA_DIR)

try:
    f = tarfile.open('genres.tar.gz', 'r:gz')
    try: 
        f.extractall()
    finally: 
        f.close()
finally:
    os.chdir(cwd)

转换为 WAV 格式

果然，如果我们想在我们的私人 MP3 收藏上测试我们的分类器，我们将无法提取太多的意义。这是因为 MP3 是一种有损音乐压缩格式，它会剪切掉人耳无法感知的部分。这很适合存储，因为有了 MP3，你可以在你的设备上存储 10 倍多的歌曲。然而，对于我们的努力来说，情况并不那么好。对于分类，我们使用 WAV 文件会更容易，因为它们可以被scipy.io.wavfile包直接读取。因此，如果我们想使用我们的分类器，我们必须转换我们的 MP3 文件。

如果附近没有转换工具，您可能想查看 SoX:sox.sourceforge.net。它声称是声音处理的瑞士军刀，我们同意这个大胆的说法。

然而，GTZAN 数据集附带的音乐文件不是 MP3 格式，而是 AU 格式，这意味着我们必须一个文件一个文件地转换它。以下片段是 Jupyter 笔记本中可能出现的一个巧妙技巧:它方便地允许我们运行系统命令，例如 Python 环境中的sox声音转换器。我们只需在命令行前面加上感叹号(!)并使用花括号传递 Python 表达式:

GENRE_DIR = Path(DATA_DIR) / 'genres' 
# You need to adapt the SOX_PATH accordingly on your system
SOX_PATH = r'C:\Program Files (x86)\sox-14-4-2'
for au_fn in Path(GENRE_DIR).glob('**/*.au'):
   print(au_fn)
    !"{SOX_PATH}/sox.exe" {au_fn} {au_fn.with_suffix('.wav')}

当然，所有这些都可以在普通的 Linux 或 Windows 外壳中完成，但是需要更多的外壳专业知识。

我们所有的音乐文件都是 WAV 格式的一个优点是它可以被 SciPy 工具包直接读取:

>>> sample_rate, X = scipy.io.wavfile.read(wave_filename)

X现在包含样本，sample_rate是获取样本的速率。让我们利用这些信息来浏览一些音乐文件，以了解数据是什么样子的。

看音乐

快速了解不同流派的歌曲“外观”的一个非常方便的方法是为一个流派中的一组歌曲绘制一个声谱图。声谱图是歌曲中出现的频率的视觉表示。它在 x 轴上显示指定时间间隔内 y 轴上频率的强度。在下面的声谱图中，这意味着在歌曲的特定时间窗口中，颜色越亮，频率越强。

Matplotlib 提供了方便的specgram()功能，为我们执行大部分的幕后计算和绘图:

>>> import scipy.io.wavfile
>>> 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), cmap='hot')

我们刚刚读入的 WAV 文件是以22050 Hz 的速率采样的，并且包含661794个样本。

如果我们现在为不同的 WAV 文件绘制这前 30 秒的声谱图，我们可以看到相同流派的歌曲之间存在共性，如下图所示:

只要看一眼图像，我们就能立即看到例如金属和古典歌曲之间的光谱差异。虽然金属歌曲在大部分频谱上一直具有高强度(它们充满活力！)，古典歌曲呈现出更加多样化的格局。

应该可以训练出一个分类器，它至少能以足够高的准确度区分金属和古典歌曲。不过，乡村和摇滚等其他类型的组合可能会构成更大的挑战。这对我们来说似乎是一个真正的挑战，因为我们不仅需要辨别两个类别，还需要辨别六个类别。我们需要能够合理地区分它们。

将音乐分解成正弦波成分

我们的计划是从原始样本读数(之前存储在X中)中提取单个频率强度，并将其输入分类器。这些频率强度可以通过应用快速傅立叶变换 ( 快速傅立叶变换)来提取，快速傅立叶变换将波信号转换成其频率分量的系数。由于快速傅立叶变换背后的理论超出了本章的范围，让我们只看一个例子来了解它实现了什么。稍后，我们将把它当作一个黑盒特征提取器。

比如我们生成两个 WAV 文件，sine_a.wav和sine_b.wav，分别包含 400 Hz 和 3000Hz 正弦波的声音。前面提到的瑞士军刀sox，是在命令行上实现这一点的一种方法(或者通过在前面加上感叹号直接从 Jupyter 获得):

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 秒。在下面的图像中，我们可以看到正弦波的快速傅立叶变换。毫不奇怪，我们在相应的正弦波下面看到400 Hz 和3000 Hz 的尖峰。

现在，让我们将两者混合，使 400 赫兹的声音只有 3000 赫兹声音的一半:

sox --combine mix --volume 1 sine_b.wav --volume 0.5 sine_a.wav 
sine_mix.wav

我们在组合声音的快速傅立叶变换图中看到两个尖峰，其中 3000 赫兹的尖峰几乎是 400 赫兹的两倍:

对于真正的音乐，我们很快会发现 FFT 看起来并不像前面的例子中那么漂亮:

使用快速傅立叶变换构建我们的第一个分类器

我们现在可以使用 FFT 创建歌曲的音乐指纹。如果我们对几首歌曲这样做，并手动分配它们相应的流派作为标签，我们就有了可以输入到第一个分类器中的训练数据。

增加实验灵活性

在我们深入分类器训练之前，让我们考虑一下实验敏捷性。虽然我们在快速傅立叶变换中有“快速”这个词，但它比我们基于文本的章节中的功能创建要慢得多。因为我们还处于实验阶段，我们可能会考虑如何加快整个特征创建过程。

当然，每次运行分类器时，每个文件的快速傅立叶变换的创建都是相同的。因此，我们可以缓存它并读取缓存的 FFT 表示，而不是完整的 WAV 文件。我们用create_fft()函数来做，该函数反过来使用scipy.fft()来创建快速傅立叶变换。为了简单(和速度！)，让我们将本例中 FFT 组件的数量固定为前 1000 个。以我们目前的知识，我们不知道这些对于音乐类型分类是否是最重要的，只知道它们在前面的快速傅立叶变换例子中显示了最高的强度。如果我们以后想使用更多或更少的快速傅立叶变换组件，我们必须为每个声音文件重新创建快速傅立叶变换表示:

import numpy as np
import scipy

def create_fft(fn):
    sample_rate, X = scipy.io.wavfile.read(fn)

    fft_features = abs(scipy.fft(X)[:1000])
    np.save(Path(fn).with_suffix('.fft'), fft_features)

for wav_fn in Path(GENRE_DIR).glob('**/*.wav'):
    create_fft(wav_fn)

我们使用 NumPy 的save()函数保存数据，该函数总是将.npy附加到文件名中。我们只需要为训练或预测所需的每个 WAV 文件做一次。

对应的 FFT 读取功能为read_fft():

def read_fft(genre_list, base_dir=GENRE_DIR):
    X = []
    y = []
    for label, genre in enumerate(genre_list):
        genre_dir = Path(base_dir) / genre
        for fn in genre_dir.glob("*.fft.npy"):
            fft_features = np.load(fn)

            X.append(fft_features[:1000])
            y.append(label)

    return np.array(X), np.array(y)

在我们的加扰music目录中，我们期待以下音乐流派:

GENRES = ["classical", "jazz", "country", "pop", "rock", "metal"]

训练分类器

让我们使用逻辑回归分类器，它已经在第 9 章、分类二–情感分析中为我们提供了很好的服务:

from sklearn.linear_model.logistic import LogisticRegression
def create_model():
    return LogisticRegression()

仅提一个令人惊讶的方面:当第一次从二进制转换到多类分类时的准确率评估。在二进制分类问题中，我们知道 50%的准确率是最坏的情况，因为这可以通过随机猜测来实现。在多类设置中，50%已经很好了。以我们的六个流派为例，随机猜测的结果只有 16.7%(假设班级规模相等)。

完整的培训过程如下所示:

from collections import defaultdict
from sklearn.metrics import precision_recall_curve, roc_curve, \
                            confusion_matrix
from sklearn.metrics import auc
from sklearn.model_selection import ShuffleSplit

def train_model(clf_factory, X, Y):
    labels = np.unique(Y)

    cv = ShuffleSplit(n_splits=1, test_size=0.3, random_state=0)

    train_errors = []
    test_errors = []

    scores = []
    pr_scores = defaultdict(list)
    precisions = defaultdict(list)
    recalls = defaultdict(list)
    thresholds = defaultdict(list)

    roc_scores = defaultdict(list)
    tprs = defaultdict(list)
    fprs = defaultdict(list)

    clfs = [] # used to later get the median

    cms = []

    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)
        clfs.append(clf)

        train_score = clf.score(X_train, y_train)
        test_score = clf.score(X_test, y_test)
        scores.append(test_score)

        train_errors.append(1 - train_score)
        test_errors.append(1 - test_score)

        y_pred = clf.predict(X_test)
        cm = confusion_matrix(y_test, y_pred) # will be explained soon
        cms.append(cm)

        for label in labels:
            y_label_test = np.asarray(y_test == label, dtype=int)
            proba = clf.predict_proba(X_test)
            proba_label = proba[:, label]

            precision, recall, pr_thresholds = preci-sion_recall_curve(
                y_label_test, proba_label)
            pr_scores[label].append(auc(recall, precision))
            precisions[label].append(precision)
            recalls[label].append(recall)
            thresholds[label].append(pr_thresholds)

            fpr, tpr, roc_thresholds = roc_curve(y_label_test, 
                                                        pro-ba_label)
            roc_scores[label].append(auc(fpr, tpr))
            tprs[label].append(tpr)
            fprs[label].append(fpr)

    all_pr_scores = np.asarray(pr_scores.values()).flatten()
    summary = (np.mean(scores), np.std(scores),
                 np.mean(all_pr_scores), np.std(all_pr_scores))
    print("%.3f\t%.3f\t%.3f\t%.3f\t" % summary)

    return np.mean(train_errors), np.mean(test_errors), np.asarray(cms)

整个训练调用如下:

X, Y = read_fft(GENRES)
train_avg, test_avg, cms = train_model(create_model, X, Y)

在多类问题中使用混淆矩阵来测量精确度

对于多类问题，我们不应该只对如何正确地对体裁进行分类感兴趣。我们还应该调查哪些体裁我们彼此混淆。这可以通过适当命名的混淆矩阵来完成，您可能已经注意到这是培训过程的一部分:

>>> cm = confusion_matrix(y_test, y_pred)

如果我们打印出混淆矩阵，我们会看到如下内容:

 [[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 as plt

def plot_confusion_matrix(cm, genre_list, name, title):
 plt.clf()
 plt.matshow(cm, fignum=False, cmap='Blues', vmin=0, vmax=1.0)
 ax = plt.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)
 ax.tick_params(axis='both', which='both', bottom='off', left='off')
 plt.title(title)
 plt.colorbar()
 plt.grid(False)
 plt.show()
 plt.xlabel('Predicted class')
 plt.ylabel('True class')
 plt.grid(False)

创建混淆矩阵时，一定要选择一个颜色顺序合适的颜色映射图(matshow()的cmap参数)，这样就可以立即看到较浅或较深的颜色意味着什么。特别不推荐这类图形是彩虹色地图，例如，matplotlib实例默认 jet 甚至成对的彩色地图。

最终的图表如下所示:

对于一个完美的分类器，我们希望左上角到右下角有一个对角的深色方块，其余区域有浅色。在上图中，我们立即看到我们基于快速傅立叶变换的分类器远非完美。它只能正确预测古典歌曲(暗方)。例如，对于岩石，大多数时候它更喜欢标签金属。

显然，使用快速傅立叶变换为我们指出了正确的方向(古典流派没有那么糟糕)，但不足以获得一个像样的分类器。当然，我们可以使用快速傅立叶变换组件的数量(固定为 1000)。但是在我们深入研究参数调整之前，我们应该先做研究。在那里，我们发现快速傅立叶变换确实是流派分类的一个不错的特征——它只是不够精炼。很快，我们将看到如何通过使用它的处理版本来提高我们的分类性能。

然而，在此之前，我们将学习另一种测量分类性能的方法。

使用接收器-操作器特性测量分类器性能的另一种方法

我们已经了解到，测量精度不足以真正评估分类器。相反，我们依靠精确-回忆(P/R)曲线来更深入地理解我们的分类器是如何工作的。

有一个 P/R 曲线的姊妹曲线，称为接收器-操作者-特征 ( ROC )，它测量分类器性能的相似方面，但提供了分类性能的另一个视图。关键的区别在于，P/R 曲线更适合于正类比负类有趣得多的任务，或者正例数比负例数少得多的任务。信息检索和欺诈检测是典型的应用领域。另一方面，ROC 曲线更好地描述了分类器的总体表现。

为了更好地理解差异，让我们考虑一下先前训练的分类器在正确分类乡村歌曲方面的性能，如下图所示:

在左边，我们看到了市盈率曲线。对于一个理想的分类器，我们会让曲线从左上角直接到右上角，然后到右下角，从而产生 1.0 的曲线下面积 ( AUC )。

右图描绘了相应的 ROC 曲线。它绘制了真阳性率 ( TPR )与假阳性率 ( FPR )的关系图。这里，理想的分类器会有一条从左下角到左上角，然后到右上角的曲线。随机分类器将是从左下角到右上角的直线，如虚线所示，其 AUC 为 0.5。因此，我们不能将市盈率曲线的 AUC 与 ROC 曲线的 AUC 进行比较。

与曲线无关，当在同一数据集上比较两个不同的分类器时，我们总是可以安全地假设一个分类器的 P/R 曲线的较高 AUC 也意味着相应 ROC 曲线的较高 AUC，反之亦然。因此，我们从不费心去产生两者。关于这一点的更多信息可以在 Davis 和 Goadrich 的非常有见地的论文精确-回忆和 ROC 曲线之间的关系中找到(ICML，2006)。

下表总结了市盈率和 ROC 曲线之间的差异:

| | x 轴 | y 轴 | | 损益 | | | | 皇家对空观察队 | | |

查看“ x 和 y 轴的定义，我们看到 ROC 曲线 y 轴的 TPR 与 P/R 图 x 轴的 Recall 相同。

FPR 测量了被错误归类为阳性的真实阴性样本的比例，范围从完美情况下的0(无false阳性)到1(均为false阳性)。

接下来，让我们使用 ROC 曲线来衡量我们的分类器的性能，以获得更好的感觉。我们的多类问题的唯一挑战是 ROC 和 P/R 曲线都假设一个二元分类问题。因此，出于我们的目的，让我们为每个流派创建一个图表，显示分类器是如何执行一对一分类的:

from sklearn.metrics import roc_curve
y_pred = clf.predict(X_test)for label in labels:
    y_label_test = np.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 图(同样，完整的代码，请遵循随附的 Jupyter 笔记本)。正如我们已经发现的，我们的第一个版本的分类器只在古典歌曲上表现良好。然而，从单个 ROC 曲线来看，我们确实在其他大部分流派中表现不佳。只有爵士乐和乡村音乐提供了一些希望。其余流派显然不可用:

利用 mel 频率倒谱系数提高分类性能

我们已经了解到，快速傅立叶变换为我们指明了正确的方向，但其本身不足以最终得出一个分类器，成功地将我们的歌曲加扰目录组织成单个流派目录。我们需要一个更先进的版本。

在这一点上，我们必须做更多的研究。其他人过去可能也遇到过类似的挑战，并且已经找到了可能对我们有帮助的方法。事实上，甚至还有一个由国际音乐信息检索协会组织的年度音乐流派分类会议。显然，自动音乐流派分类 ( AMGC )是音乐信息检索的一个既定子领域。浏览一些 AMGC 论文，我们可以看到有一堆针对自动体裁分类的工作可能会帮助我们。

一种似乎在许多情况下成功应用的技术叫做 mel 频率倒频谱 ( MFC )系数。MFC 对声音的功率谱进行编码，功率谱是声音包含的每个频率的功率。它被计算为信号频谱对数的傅里叶变换。如果听起来太复杂，只需记住倒谱这个名字来源于谱，前四个字符颠倒即可。MFC 已成功应用于语音和说话人识别。让我们看看它是否也适用于我们。我们很幸运，因为其他人已经确切地需要这个，并且发布了它的实现作为python_speech_features模块的一部分。我们可以用pip轻松安装。之后，我们可以调用mfcc()函数，计算 MFC 系数，如下所示:

>>> from python_speech_features import mfcc

>>> fn = Path(GENRE_DIR) / 'jazz' / 'jazz.00000.wav'
>>> sample_rate, X = scipy.io.wavfile.read(fn)
>>> ceps = mfcc(X)
>>> print(ceps.shape)
 (4135, 13)

ceps包含歌曲的每个4135帧的13系数(默认为mfcc()的num_ceps参数)。获取所有的数据会让我们的分类器不堪重负。相反，我们可以做的是对所有帧的每个系数取平均值。假设每首歌的开头和结尾可能没有中间部分那么特定于流派，我们也忽略了第一个和最后 10%:

>>> num_ceps = ceps.shape[0]
>>> np.mean(ceps[int(num_ceps*0.1):int(num_ceps*0.9)], axis=0)
array([ 16.43787597, 7.44767565, -13.48062285, -7.49451887,
 -8.14466849, -4.79407047, -5.53101133, -5.42776074,
 -8.69278344, -6.41223865, -3.01527269, -2.75974429, -3.61836327])

果然，我们将使用的基准数据集只包含每首歌的前 30 秒，所以我们不需要剪掉最后 10%。我们无论如何都会这样做，这样我们的代码也可以在其他数据集上工作，这些数据集很可能不会被截断。

与我们使用 FFT 的工作类似，我们也希望缓存曾经生成的 MFCC 特征并读取它们，而不是每次训练分类器时都重新创建它们。

这将导致以下代码:

def create_ceps(fn):
    sample_rate, X = scipy.io.wavfile.read(fn)
    np.save(Path(fn).with_suffix('.ceps'), mfcc(X))

for wav_fn in Path(GENRE_DIR).glob('**/*.wav'):
    create_fft(wav_fn)
def read_ceps(genre_list, base_dir=GENRE_DIR):
    X = []
    y = []
    for label, genre in enumerate(genre_list):
        genre_dir = Path(base_dir) / genre
        for fn in genre_dir.glob("*.ceps.npy"):
            ceps = np.load(fn)
            num_ceps = len(ceps)
            X.append(np.mean(ceps[int(num_ceps / 10):int(num_ceps * 9 / 10)], axis=0))
            y.append(label)

    return np.array(X), np.array(y)

使用每首歌曲仅使用 13 个特征的分类器，我们得到了以下有希望的结果:

所有流派的分类性能都有所提高。古典和金属的 AUC 几乎都在 1.0。事实上，以下情节中的混乱矩阵现在看起来好多了。我们可以清楚地看到对角线，表明分类器在大多数情况下能够正确地对流派进行分类。这个分类器实际上非常适用于解决我们的初始任务:

如果我们想改进这一点，这个混淆矩阵会很快告诉我们应该关注什么:非对角线位置上的非白点。例如，我们有一个更黑暗的地方，我们错误地将摇滚歌曲贴上了爵士的标签，这种可能性很大。为了解决这个问题，我们可能需要更深入地研究歌曲，提取一些东西，比如鼓的模式和类似的特定流派的特征。然后，在浏览 ISMIR 论文的同时，我们还阅读了关于听觉滤波器组时间包络 ( AFTE )特征的文章，这些特征在某些情况下似乎优于 MFCC 特征。也许我们也应该看看他们？

好的一点是，只有配备了 ROC 曲线和混淆矩阵，我们才可以自由地在特征提取器方面引入其他专家的知识，而不必完全了解他们的内部工作方式。我们的测量工具总是会告诉我们什么时候方向是对的，什么时候该改变。当然，作为一个渴望学习的机器学习者，我们总会有一种感觉，在我们的特征提取器的黑匣子里，某个地方埋藏着一个令人兴奋的算法，就等着我们被理解。

使用张量流的音乐分类

我们可以用我们的特征输入张量流吗？当然可以！但是让我们试着利用这个机会实现另外两个目标:

  • 我们将使 TensforFlow 分类器的行为类似于sklearn分类器，以便在所有兼容函数中重用。
  • 即使神经网络可以提取任何特征，它们仍然需要被设计和训练来提取它们。在这个例子中，从原始声音文件开始，我们将向您展示，获得比倒谱系数更好的结果是不够的。

但是让我们切入正题，设置我们的超参数:

import tensorflow as tf
import numpy as np

n_epochs = 50
learning_rate = 0.01
batch_size = 128
step = 32
dropout_rate = 0.2

signal_size = 1000
signal_shape = [signal_size,1]

我们从 600 个样本开始，但为了向训练中添加更多数据，我们将把文件分成几个块:

def read_wav(genre_list, multiplicity=1, base_dir=GENRE_DIR):
    X = []
    y = []
    for label, genre in enumerate(genre_list):
        genre_dir = Path(base_dir) / genre
        for fn in genre_dir.glob("*.wav"):
            sample_rate, new_X = scipy.io.wavfile.read(fn)
            for i in range(multiplicity):
                X.append(new_X[i*signal_size:(i+1)*signal_size])
                y.append(label)

    return np.array(X).reshape((-1, signal_size, 1)), np.array(y)

从每个文件中，我们将获得 20 个较短的样本:

from sklearn.model_selection import train_test_split

X, Y = read_wav(GENRES, 20)
classes = len(GENRES)
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, 
                                     test_size=(1\. / 6.))

我们的网络将非常类似于图像分类网络。我们将使用 1D 层代替 2D 卷积层。我们还将添加一个获取池大小的参数。以前的网络使用小的 2D 池，我们可能不得不使用更大的池。

我们还将使用脱落层。由于样本不多，我们不得不避免网络泛化不良。这将帮助我们实现这一目标:

class CNN():
    def __init__(
            self,
            signal_shape=[1000,1],
            dim_W1=64,
            dim_W2=32,
            dim_W3=16,
            classes=6,
            kernel_size=5,
            pool_size=16
            ):

        self.signal_shape = signal_shape

        self.dim_W1 = dim_W1
        self.dim_W2 = dim_W2
        self.dim_W3 = dim_W3
        self.classes = classes
        self.kernel_size = kernel_size
        self.pool_size = pool_size

    def build_model(self):
        image = tf.placeholder(tf.float32, [None]+self.signal_shape, name="signal")
        Y = tf.placeholder(tf.int64, [None], name="label")

        probabilities = self.discriminate(image, training)
        cost = tf.reduce_mean(tf.nn.sparse_softmax_cross_entropy_with_logits(labels=Y, 
                logits=probabilities))
        accuracy = tf.reduce_mean(tf.cast(tf.equal(tf.argmax(probabilities, ax-is=1), 
                        Y), tf.float32), name="accuracy")

        return image, Y, cost, accuracy, probabilities

在这里，我们重用我们在第 8 章、人工神经网络和深度学习中看到的sparse_softmax_cross_entropy_with_logits cost-helper函数。提醒一下，它将整数目标与图层进行比较，以便具有最大值的节点与该目标匹配:

    def create_conv1d(self, input, filters, kernel_size, name):
        layer = tf.layers.conv1d(
                    inputs=input,
                    filters=filters,
                    kernel_size=kernel_size,
                    activation=tf.nn.leaky_relu,
                    name="Conv1d_" + name,
                    padding="same")
        return layer

    def create_maxpool(self, input, name):
        layer = tf.layers.max_pooling1d(
                    inputs=input,
                    pool_size=[self.pool_size],
                    strides=self.pool_size,
                    name="MaxPool_" + name)
        return layer

    def create_dropout(self, input, name, is_training):
        layer = tf.layers.dropout(
                    inputs=input,
                    rate=dropout_rate,
                    name="DropOut_" + name,
                    training=is_training)
        return layer

    def create_dense(self, input, units, name):
        layer = tf.layers.dense(
                inputs=input,
                units=units,
                name="Dense" + name,
                )
        layer = tf.layers.batch_normalization(
                inputs=layer,
                momentum=0,
                epsilon=1e-8,
                training=True,
                name="BatchNorm_" + name,
        )
        layer = tf.nn.leaky_relu(layer, name="LeakyRELU_" + name)
        return layer

    def discriminate(self, signal, training):
        h1 = self.create_conv1d(signal, self.dim_W3, self.kernel_size, "Layer1")
        h1 = self.create_maxpool(h1, "Layer1")

        h2 = self.create_conv1d(h1, self.dim_W2, self.kernel_size, "Layer2")
        h2 = self.create_maxpool(h2, "Layer2")
        h2 = tf.reshape(h2, (-1, self.dim_W2 * h2.shape[1]))

        h3 = self.create_dense(h2, self.dim_W1, "Layer3")
        h3 = self.create_dropout(h3, "Layer3", training)

        h4 = self.create_dense(h3, self.classes, "Layer4")
        return h4

如前所述，可以将我们的网络封装在一个sklearn对象BaseEstimator中。这些估计器有一组从构造函数中提取的参数，正如我们在这里看到的:

from sklearn.base import BaseEstimator

class Classifier(BaseEstimator):
    def __init__(self,
            signal_shape=[1000,1],
            dim_W1=64,
            dim_W2=32,
            dim_W3=16,
            classes=6,
            kernel_size=5,
            pool_size=16):
        self.signal_shape=signal_shape
        self.dim_W1=dim_W1
        self.dim_W2=dim_W2
        self.dim_W3=dim_W3
        self.classes=classes
        self.kernel_size=kernel_size
        self.pool_size=pool_size

fit方法是创建和训练模型的方法。这里我们还保存了网络和变量的状态:

    def fit(self, X, y):
        tf.reset_default_graph()

        print("Fitting (W1=%i) (W2=%i) (W3=%i) (kernel=%i) (pool=%i)"
              % (self.dim_W1, self.dim_W2, self.dim_W3, self.kernel_size, self.pool_size))

        cnn_model = CNN(
                signal_shape=self.signal_shape,
                dim_W1=self.dim_W1,
                dim_W2=self.dim_W2,
                dim_W3=self.dim_W3,
                classes=self.classes,
                kernel_size=self.kernel_size,
                pool_size=self.pool_size
                )

        signal_tf, Y_tf, cost_tf, accuracy_tf, output_tf = cnn_model.build_model()
        train_step = tf.train.AdamOptimizer(learning_rate, be-ta1=0.5).minimize(cost_tf)

        saver = tf.train.Saver()

        with tf.Session() as sess:
            sess.run(tf.global_variables_initializer())
            for epoch in range(n_epochs):
                permut = np.random.permutation(len(X_train))
                for j in range(0, len(X_train), batch_size):
                    batch = permut[j:j+batch_size]
                    Xs = X_train[batch]
                    Ys = Y_train[batch]

                    sess.run(train_step,
                            feed_dict={
                                Y_tf: Ys,
                                signal_tf: Xs
                                })
            saver.save(sess, './classifier')
        return self

然后在predict方法中恢复它们，我们将使用训练好的网络对新数据进行分类:

    def predict(self, X):
        tf.reset_default_graph()
        new_saver = tf.train.import_meta_graph("classifier.meta") 
        with tf.Session() as sess: 
            new_saver.restore(sess, tf.train.latest_checkpoint('./'))

            graph = tf.get_default_graph()
            training_tf = graph.get_tensor_by_name('is_training:0')
            signal_tf = graph.get_tensor_by_name('signal:0')
            output_tf = graph.get_tensor_by_name('LeakyRELU_Layer4/Maximum:0')

            predict = sess.run(output_tf,
                            feed_dict={
                                training_tf: False,
                                signal_tf: X
                                })
            return np.argmax(predict, axis=1)

我们还没有在这个估算器中创建一个score方法，但是我们可以使用sklearn API 从predict方法中创建一个。

现在我们将使用这个估计量，并进行网格搜索，以找到一组合适的超参数。由于我们没有很多样本，我们将在每个卷积层以及密集层上仅使用少数几个单元。我们还想提取更好的过滤器。毕竟，声音通常需要更大的过滤器来提取有意义的特征:

from sklearn.model_selection import GridSearchCV
from sklearn.metrics import accuracy_score, make_scorer

param_grid = {
    "dim_W1": [4, 8, 16],
    "dim_W2": [4, 8, 16],
    "dim_W3": [4, 8, 16],
    "kernel_size":[7, 11, 15],
    "pool_size":[8, 12, 16],
}

cv = GridSearchCV(Classifier(), param_grid, scor-ing=make_scorer(accuracy_score), cv=6)


cv.fit(X, Y)
print(cv.best_params_)

{'dim_W1': 4, 'dim_W2': 4, 'dim_W3': 16, 'kernel_size': 15, 'pool_size': 12}

现在我们已经花了几个小时找到了一组足够的超参数，我们可以用它来检查混淆矩阵:

clf = Classifier(**cv.best_params_)
clf.fit(X_train, Y_train)

Y_train_predict = clf.predict(X_train)
Y_test_predict = clf.predict(X_test)

from sklearn.metrics import confusion_matrix
cm = confusion_matrix(Y_train, Y_train_predict)
plot_confusion_matrix(cm / np.sum(cm, axis=0), GENRES, "CNN",
    "Confusion matrix of a CNN based classifier (train)")
cm = confusion_matrix(Y_test, Y_test_predict)
plot_confusion_matrix(cm / np.sum(cm, axis=0), GENRES, "CNN",
    "Confusion matrix of a CNN based classifier (test)")

请参考以下图表:

即使有最好的参数，这也表明了我们之前所说的:你也需要有足够的特性。你可以用电脑来训练一个更好的网络，但是你也需要更多的数据。

在倒频谱特征上尝试此网络(或具有不同层配置的网络)。能达到更好的分类吗？它能匹配我们之前创建的最佳LogisticRegression分类器吗？来看看吧！

摘要

在这一章中，当我们构建音乐类型分类器时，我们走出了舒适区。由于对音乐理论了解不深，一开始我们没有训练出一个分类器，用 FFT 以合理的精度预测歌曲的音乐流派。但是，然后，我们创建了一个使用 MFC 特性显示真正可用性能的分类器。

在这两种情况下，我们使用的特征，我们只了解知道如何和在哪里把它们放在我们的分类器设置。第一次失败，第二次成功。它们之间的区别在于，在第二种情况下，我们依赖于该领域专家创建的功能。

这是完全可以的。如果我们主要对结果感兴趣，我们有时只需要走捷径——我们只需要确保这些捷径来自特定领域的专家。因为我们已经学会了如何在这个新的多类分类问题中正确衡量性能，所以我们自信地走了这些捷径。

在第 12 章、计算机视觉中，将看一下图像处理、特征表示、CNN 和 GAN。*

十二、计算机视觉

图像分析和计算机视觉在工业和科学应用中一直很重要。随着具有强大摄像头和互联网连接的手机的普及，现在图像越来越多地由消费者生成。因此，有机会利用计算机视觉在新的环境中提供更好的用户体验。

在本章中，我们将了解如何将您在本书其余部分中学习的几种技术应用于这种特定类型的数据。特别是，我们将学习如何使用mahotas计算机视觉包从图像中提取特征。然后，这些特征可以用作我们在其他章节中研究的相同分类方法的输入。我们将把这些技术应用于公开的照片数据集。我们还将看到如何使用相同的特征来寻找相似的图像。我们还将学习如何使用本地功能。这些是相对通用的，并且在许多任务中获得非常好的结果(尽管它们具有较高的计算成本)。

最后，在本章的最后，我们将使用 Tensorflow 基于现有数据集生成新图像。特别是，在本章中，我们将执行以下操作:

  • 了解如何将图像表示为 NumPy 数组并对其进行操作
  • 了解如何将图像表示为一组小特征，以便在这种数据类型上使用标准分类和聚类方法
  • 学习如何使用视觉单词来生成另一种类型的特征
  • 了解如何生成与现有图像相似的新图像

介绍图像处理

从计算机的角度来看，图像是像素值的大矩形阵列。我们的目标是处理一个或多个图像，并为我们的应用程序做出决定。

根据设置，它可能是一个分类问题，一个聚类问题，或者我们在书中看到的任何其他问题类别。

第一步是从磁盘加载图像，图像通常以特定于图像的格式存储，如巴布亚新几内亚或 JPEG，前者是无损压缩格式，后者是有损压缩，是针对照片的视觉评估而优化的。然后，我们可能希望对图像进行预处理(例如，针对光照变化对图像进行归一化)。

加载和显示图像

为了操纵图像，我们将使用一个名为mahotas的包。您可以通过蟒蛇获取mahotas，并在httpT4T6】上阅读其手册://maho tas . read docs .T8】io。Mahotas 是一个开源包(它拥有麻省理工学院的许可，因此可以在任何项目中使用)，由本书的作者之一开发。它基于 NumPy。因此，您迄今为止获得的 NumPy 知识可以用于图像处理。还有其他的图像包，比如scikit-image(skipage)n 维图像(ndi image)模块在 SciPy 中，以及 OpenCV 的 Python 绑定。所有这些都与 NumPy 数组一起工作，因此您甚至可以混合和匹配来自不同包的功能来构建一个组合管道。

我们首先用mh缩写导入mahotas，我们将在本章中使用该缩写，如下所示:

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 can use a variety of different input/output backends. Unfortunately, none of them can load all image formats that exist (there are hundreds, with several variations of each). However, the loading of PNG and JPEG images is supported by all of them. We will focus on these common formats and refer you to the mahotas documentation on how to read uncommon formats.

我们可以使用matplotlib在屏幕上显示图像，我们已经使用过几次的plotting库如下:

from matplotlib import pyplot as plt 
fig,ax = plt.subplots() 
ax.imshow(image)

如下面的截图所示，这段代码使用第一个维度是高度，第二个维度是宽度的约定来显示图像。它也能正确处理彩色图像。当使用 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) 
fig,ax = plt.subplots() 
ax.imshow(image) # Display the image

默认情况下，matplotlib 会将此单通道图像显示为假彩色图像，使用红色表示高值，蓝色表示低值。对于自然图像，灰度更合适。您可以通过以下方式选择它:

plt.gray()

现在图像以灰度显示。请注意，只有解释和显示像素值的方式发生了变化，图像数据未被触及。我们可以通过计算阈值来继续我们的处理。阈值化是两类聚类的一种形式，对于这一任务有几种方法。Otsu就是这样一种thresholding方法，它试图找到两组紧凑的像素组，高于阈值和低于阈值的像素组:

thresh = mh.thresholding.otsu(image)
print('Otsu threshold is {}.'.format(thresh))
Otsu threshold is 138. 
fig,ax = plt.subplots()
ax.imshow(image > thresh)

该方法应用于上一张图像时，发现threshold为138，将地面与上方天空分开，如下图截图所示:

高斯模糊

模糊你的图像可能看起来很奇怪，但它通常有助于减少噪音，这有助于进一步的处理。有了mahotas，就只是一个函数调用:

im16 = mh.gaussian_filter(image, 16)

请注意，我们没有将灰度图像转换为无符号整数；我们只是按原样使用了浮点结果。gaussian_filter函数的第二个参数是过滤器的大小(过滤器的标准偏差)。较大的值会导致更多的模糊，如下图所示:

我们可以使用前面的截图和带有Otsu的threshold(使用前面的代码)。现在，边界更加平滑，没有锯齿边缘，如下图所示:

聚焦中心

最后一个例子向您展示了如何将 NumPy 运算符与一点点过滤混合在一起，以获得一个有趣的结果。我们从森林中一条小路的照片开始:

im = mh.imread('forest')

要分割红色、绿色和蓝色通道，我们使用以下代码。NumPy transpose方法改变多维数组中轴的顺序:

r,g,b = im.transpose(2,0,1)

现在，我们分别过滤这三个通道，并用mh.as_rgb从其中构建一个合成图像。该函数采用三个二维数组，执行对比度拉伸，使每个数组成为一个 8 位整数数组，然后堆叠它们，返回一个彩色 RGB 图像:

r24 = mh.gaussian_filter(r, 24.) 
g24 = mh.gaussian_filter(g, 24.) 
b24 = mh.gaussian_filter(b, 24.) 
im24 = mh.as_rgb(r24, g24, b24)

现在，我们将两幅图像从中心向边缘混合。首先，我们需要构建一个权重数组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)*im24)

最后，我们可以将两幅图像合并，使中心清晰聚焦，边缘更加柔和:

基本图像分类

我们将从专为本书收集的一个小数据集开始。它有三类:建筑物、自然场景(风景)和文字图片。每个类别有 30 张图片，它们都是用手机摄像头拍摄的，构图很少。这些图片类似于那些没有经过摄影训练的用户上传到现代网站的图片。该数据集可在伴随代码库中找到。在本章的后面，我们将看到一个更大的数据集，其中有更多的图像和更难分类的类别。

对图像进行分类时，我们从一个大的矩形数字数组(像素值)开始。如今，数百万像素很常见。我们可以尝试将所有这些数字作为特征输入到学习算法中。这不是一个很好的主意，除非你有很多数据。这是因为每个像素(甚至每个小像素组)与最终结果的关系非常间接。此外，拥有数百万像素，但仅作为少量示例图像，会导致非常困难的统计学习问题。这是我们在第三章、回归中讨论的 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函数返回一个 4×13 的数组。第一维指的是计算要素的四个可能方向(垂直、水平、对角线和反对角线)。如果我们对方向不特别感兴趣，我们可以使用所有方向的平均值(在前面的代码中显示为haralick_features_mean)。否则，我们可以单独使用所有功能(使用haralick_features_all)。这个决定应该由数据集的属性来决定。在我们的例子中，我们认为水平和垂直方向应该分开。因此，我们将使用haralick_features_all。

在mahotas中还实现了一些其他的特性集。线性二进制模式是另一种基于纹理的特征集，它对光照变化非常鲁棒。还有其他类型的功能，包括本地功能，我们将在本章后面讨论。

有了这些特性，我们使用标准的分类方法，如logistic regression，如下所示:

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 model_selection 
cv = model_selection.LeaveOneOut() 
scores = model_selection.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

这使得像素值的范围从零到三，这给了我们总共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，它计算对数(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 = model_selection.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))

我们将只绘制数据的子集(每 10 个元素)，这样查询将在顶部，返回的最近邻居在底部，如下面的代码所示:

fig, axes = plt.subplots(2, 9) 
for ci,i in enumerate(range(0,90,10)): 
    query = images[i] 
    dists_query = dists[i] 
    closest = dists_query.argsort() 
    # closest[0] is same as the query image, so pick next closest 
    closest = closest[1] 
    result = images[closest] 
    query = mh.imread(query) 
    result = mh.imread(result) 
    axes[0, ci].imshow(query) 
    axes[1, ci].imshow(result)

结果显示在下面的截图中(顶部图像是查询图像，底部图像是返回的结果):

很明显，该系统并不完美，但可以找到至少在视觉上与查询相似的图像。除了一种情况之外，找到的图像都来自与查询相同的类。

分类较难的数据集

之前的数据集是一个易于使用纹理特征进行分类的数据集。事实上，从商业角度来看，许多有趣的问题都相对容易。然而，有时我们可能会面临一个更棘手的问题，需要更好、更现代的技术来获得好的结果。

我们现在将测试一个公共数据集，它具有相同的结构:几张照片被分成少量的类。课程有动物、汽车、交通和自然景观。

与我们之前讨论的三类问题相比，这些类更难区分。自然场景、建筑和文本具有非常不同的纹理。然而，在这个数据集中，纹理和颜色不是图像类的清晰标记。以下是动物类的一个例子:

这是汽车课的另一个例子:

两个对象都是在自然背景下，对象内部有大而平滑的区域。这是一个比以前的数据集更难的问题，因此我们需要使用更高级的方法。第一个改进是使用稍微强大一点的分类器。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)])

数据在数据集中不是以随机的顺序组织的:相似的图像靠得很近。因此，我们使用一个交叉验证计划来考虑被打乱的数据，这样每个文件夹都有一个更具代表性的训练集，如下面的代码所示:

cv = model_selection.KFold(n_splits=5, 
                     shuffle=True, random_state=123) 
scores = model_selection.cross_val_score( 
   clf, ifeatures, labels, cv=cv) 
print('Accuracy: {:.1%}'.format(scores.mean())) 
Accuracy: 73.4%

这对于四个类来说还不错，但是我们现在将看看我们是否可以通过使用一组不同的特性做得更好。事实上，我们将看到，我们需要将这些特性与其他方法结合起来，以获得最佳的可能结果。

局部特征表示

与我们之前使用的特征不同，局部特征是在图像的一个小区域上计算的。Mahotas 支持名为加速鲁棒特征 ( SURF )的计算类型特征。这些特征被设计为对旋转或照明变化具有鲁棒性(也就是说，它们仅在照明变化时轻微改变它们的值)。

使用这些功能时，我们必须决定在哪里计算它们。通常使用三种可能性:

  • 随便地
  • 在网格中
  • 检测图像的感兴趣区域(一种称为关键点检测或兴趣点检测的技术)

所有这些都是有效的，并将在适当的情况下产生良好的结果。Mahotas 支持这三种方法。如果您有理由预期您的兴趣点将对应于图像中的重要区域，则使用兴趣点检测效果最佳。

我们将使用interest point方法。用mahotas计算特征很简单:导入正确的子模块，调用surf.surf函数如下:

descriptors = surf.surf(im, descriptor_only=True)

descriptors_only=True标志意味着我们只对局部特征本身感兴趣，而对它们的像素位置、大小或方向不感兴趣(单词descriptor经常用来指代这些局部特征)。或者，我们可以使用dense sampling方法，使用如下的surf.dense功能:

from mahotas.features import surf 
descriptors = surf.dense(im, spacing=16)

这将返回在相距 24 像素的点上计算的描述符的值。由于点的位置是固定的，所以关于兴趣点的元信息不是很有趣，并且默认不返回。在任一情况下，结果(描述符)都是 nx 64 数组，其中 n 是采样的点数。点数取决于图像的大小、内容以及传递给函数的参数。在这个例子中，我们使用默认设置，每个图像获得几百个描述符。

我们不能将这些描述符直接输入到支持向量机、逻辑回归机或类似的分类系统中。为了使用来自图像的描述符，有几种解决方案。我们可以对它们进行平均，但是这样做的结果不是很好，因为它们会丢弃所有特定于位置的信息。在这种情况下，我们将只有另一个基于边缘测量的全局特征集。

我们这里要用到的解决方案就是包字模型。它于 2004 年首次出版，但这是一个明显的后知后觉的想法；实现起来非常简单，取得了很好的效果。

在处理图像的时候说字可能会显得很奇怪。如果你认为你没有书面文字，很容易区分，而是口头音频，这可能更容易理解。现在，每次说一个单词，听起来都会略有不同，不同的说话人会有自己的发音。因此，一个单词的波形不会在每次说的时候都一样。然而，通过在这些波形上使用聚类，我们可以希望恢复大部分结构，以便给定单词的所有实例都在同一聚类中。即使过程不完美(也不会完美)，我们仍然可以谈论将波形分组为单词。

我们对图像数据执行相同的操作:我们将所有图像中看起来相似的区域聚集在一起，并将这些视觉单词称为。

The number of words used does not usually have a big impact on the final performance of the algorithm. Naturally, if the number is extremely small (10 or 20, when you have a few thousand images), then the overall system will not perform well. Similarly, if you have too many words (many more than the number of images, for example), the system will also not perform well. However, in between these two extremes, there is often a very large plateau, where you can choose the number of words without a big impact on the result. As a rule of thumb, using a value such as 256, 512, or 1,024 if you have many images should give you a good result.

我们将从计算以下特征开始:

alldescriptors = [] 
for im in images: 
    im = mh.imread(im, as_grey=True) 
    im = im.astype(np.uint8) 
    alldescriptors.append((surf.surf(im, descriptor_only=True)) 
# get all descriptors into a single array 
concatenated = np.concatenate(alldescriptors)

现在，我们使用 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.bincount(c, minlength=256)) 
# build single array and convert to float 
sfeatures = np.array(sfeatures, dtype=float)

这个循环的最终结果是sfeatures[fi, fj]是图像fi包含元素fj的次数。用np.histogram函数可以更快地计算出同样的结果，但是要让参数恰到好处有点棘手。我们将结果转换为浮点，因为我们不需要整数算术(及其舍入语义)。

结果是，每个图像现在由相同大小的单个特征阵列表示(在我们的例子中，簇的数量是 256)。因此，我们可以使用如下标准分类方法:

scores = model_selection.cross_val_score( 
   clf, sfeatures, labels, cv=cv) 
print('Accuracy: {:.1%}'.format(scores.mean())) 
Accuracy: 62.4%

这比以前更糟糕了！我们一无所获吗？

事实上，我们有，因为我们可以将所有特征组合在一起以获得百分之76.7的准确度，如下所示:

allfeatures = np.hstack([ifeatures, sfeatures]) scores = model_selection.cross_val_score( clf, allfeatures, labels, cv=cv) print('Accuracy: {:.1%}'.format(scores.mean())) 
Accuracy: 76.7%

这是我们拥有的最好的结果，比任何单一的特征集都好。这是由于局部 SURF 特征足够不同，以向我们之前拥有的全局图像特征添加新信息，并改善组合结果。

敌对网络下的图像生成

生成性对抗网络 ( GANs )是一种新的、潮流的网络类型。他们的主要吸引力是生殖方面。这意味着我们可以训练一个网络来生成一个类似于引用的新数据样本。

几年前，研究人员使用深度信念网络 ( DBN )来完成这项任务，它由一个可见层和一组内部层组成，最终会反复出现。训练这样的网络相当困难，所以人们考虑新的架构。

进入我们的 GAN。我们如何训练网络来生成类似于参考的样本？首先，我们需要设计一个发电机网络。通常，我们需要一组随机变量，这些变量将被输入到一组密集的conv2d_transpose层中。后者与conv2d层相反，从看起来像卷积输出的输入到看起来像卷积输入的输出。

现在，为了训练这个网络，我们使用对抗性部分。诀窍是训练另一个网络，即鉴别器，以检测样本是真实样本还是生成样本。一次迭代将训练鉴别器以增强其鉴别能力，之后的迭代将训练生成器以获得更接近真实图像的图像。

让我们尝试生成逼真的手写数字；我们将重用我们以前的 CNN 分类器的部分内容。我们需要在那里添加一个生成器，并改变图层，以考虑到将生成我们的图像的额外随机输入。

让我们从一些辅助函数开始:我们的cost函数的匹配辅助函数，以及将新生成的样本写入磁盘并在训练期间显示的函数。我们还为批处理规范化创建了自己的层，以简化底层计算:

import tensorflow as tf
import numpy as np

def match(logits, labels):
    logits = tf.clip_by_value(logits, 1e-7, 1\. - 1e-7)
    return tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits_v2(
        logits=logits, labels=labels))

We use the non-sparse version of the cost function helper we used before, softmax_cross_entropy_with_logits_v2. We could have used the sparse version as well, but at the cost of additional code. As there is only one output value, this is simple to handle.

让我们解释批处理规范化层。我们计算输入张量的平均值和标准差。然后，我们对这个矩阵进行归一化(这样，矩阵的平均值为0，标准差为1)。我们明确处理了 2D 矩阵和 4D 矩阵，因为我们需要明确处理坐标轴的差异:

def batchnormalize(X, eps=1e-8, g=None, b=None):
    if X.get_shape().ndims == 4:
        mean = tf.reduce_mean(X, [0,1,2])
        std = tf.reduce_mean( tf.square(X-mean), [0,1,2] )
        X = (X-mean) / tf.sqrt(std+eps)

        if g is not None and b is not None:
            g = tf.reshape(g, [1,1,1,-1])
            b = tf.reshape(b, [1,1,1,-1])
            X = X*g + b

    elif X.get_shape().ndims == 2:
        mean = tf.reduce_mean(X, 0)
        std = tf.reduce_mean(tf.square(X-mean), 0)
        X = (X-mean) / tf.sqrt(std+eps)

        if g is not None and b is not None:
            g = tf.reshape(g, [1,-1])
            b = tf.reshape(b, [1,-1])
            X = X*g + b

    else:
        raise NotImplementedError

    return X

def save_visualization(X, nh_nw, save_path='./sample.jpg'):
    from imageio import imwrite
    from matplotlib import pyplot as plt
    h,w = X.shape[1], X.shape[2]
    img = np.zeros((h * nh_nw[0], w * nh_nw[1], 3))

    for n,x in enumerate(X):
        j = n // nh_nw[1]
        i = n % nh_nw[1]
        img[j*h:j*h+h, i*w:i*w+w, :] = x

    img = img.astype(np.uint8)
    imwrite(save_path, img)
    plt.imshow(img)
    plt.show()

As we are using the sigmoid mapping function for probabilities, we need to remove values 0 and 1 from the mapping (as they map from infinity). We do that by adding 1e-7 to 0 and subtracting it from 1.

我们现在可以创建我们的类，创建我们的模型。构造器将有额外的新参数，类的数量Y，以及随机状态的大小Z:

class DCGAN():
    def __init__(
            self,
            image_shape=[28,28,1],
            dim_z=100,
            dim_y=10,
            dim_W1=1024,
            dim_W2=128,
            dim_W3=64,
            dim_channel=1,
            ):

        self.image_shape = image_shape
        self.dim_z = dim_z
        self.dim_y = dim_y

        self.dim_W1 = dim_W1
        self.dim_W2 = dim_W2
        self.dim_W3 = dim_W3
        self.dim_channel = dim_channel

现在，这是我们为生成器创建的新的特殊层，它将创建良好的图像——我们之前谈到的卷积转置层:

def create_conv2d_transpose(self, input, filters, kernel_size, name, with_batch_norm):
        layer = tf.layers.conv2d_transpose(
                    inputs=input,
                    filters=filters,
                    kernel_size=kernel_size,
                    strides=[2,2],
                    name="Conv2d_transpose_" + name,
                    padding="SAME")
        if with_batch_norm:
            layer = batchnormalize(layer)
            layer = tf.nn.relu(layer, name="RELU_" + name)
        return layer

我们的鉴别器应该返回真实图像的0和1之间的概率。为了实现这一点并允许生成器创建各种类型的图像，我们将图像以及它们在每一层上的类驱动到鉴别器中:

def discriminate(self, image, Y, reuse=False):
        with tf.variable_scope('discriminate', reuse=reuse):
            Y = tf.one_hot(Y, dim_y)
            yb = tf.reshape(Y, tf.stack([-1, 1, 1, self.dim_y]))
            image = tf.concat(axis=3, values=
                [image, yb*tf.ones([1, 28, 28, self.dim_y])])

            h1 = self.create_conv2d(image, self.dim_W3, 5, "Lay-er1", True)
            h1 = tf.concat(axis=3, values=
                [h1, yb*tf.ones([1, 14, 14, self.dim_y])])

            h2 = self.create_conv2d(h1, self.dim_W2, 5, "Layer2", True)
            h2 = tf.reshape(h2, tf.stack([-1, 7*7*128]))
            h2 = tf.concat(axis=1, values=[h2, Y])

            h3 = self.create_dense(h2, self.dim_W1, "Layer3", True)
            h3 = tf.concat(axis=1, values=[h3, Y])

            h4 = self.create_dense(h3, 1, "Layer4", True)
            return h4

正如我们之前所说的，生成器做相反的事情，从我们的类变量和随机状态到最终生成的值在0和1之间的图像:

def generate(self, Z, Y, reuse=False):
        with tf.variable_scope('generate', reuse=reuse):

            Y = tf.one_hot(Y, dim_y)
            yb = tf.reshape(Y, tf.stack([-1, 1, 1, self.dim_y]))
            Z = tf.concat(axis=1, values=[Z,Y])
            h1 = self.create_dense(Z, self.dim_W1, "Layer1", False)
            h1 = tf.concat(axis=1, values=[h1, Y])
            h2 = self.create_dense(h1, self.dim_W2*7*7, "Layer2", False)
            h2 = tf.reshape(h2, tf.stack([-1,7,7,self.dim_W2]))
            h2 = tf.concat(axis=3, values=
                [h2, yb*tf.ones([1, 7, 7, self.dim_y])])

            h3 = self.create_conv2d_transpose(h2, self.dim_W3, 5, "Layer3", True)
            h3 = tf.concat(axis=3, values=
                [h3, yb*tf.ones([1, 14,14,self.dim_y])] )

            h4 = self.create_conv2d_transpose(
                h3, self.dim_channel, 7, "Layer4", False)
            x = tf.nn.sigmoid(h4)
            return x

现在是组装零件的时候了。我们为生成器以及真实图像输入创建占位符。然后我们创建我们的图像生成器(我们将使用它来显示我们生成的图像)，然后创建我们的鉴别器。诀窍就在这里。我们同时创造了两个。一个用于实像，应返回1；另一个将被输入生成的图像并返回0。由于两者共享相同的权重，我们为第二个鉴别器传递重用标志。

我们试图为鉴别器步骤优化的成本是两个鉴别器和生成器的差异之和。如前所述，我们优化了生成器，以在鉴别器上获得一个1:

def build_model(self):
    Z = tf.placeholder(tf.float32, [None, self.dim_z])
    Y = tf.placeholder(tf.int64, [None])

    image_real = tf.placeholder(tf.float32, [None]+self.image_shape)
    image_gen = self.generate(Z, Y)

    raw_real = self.discriminate(image_real, Y, False)
    raw_gen = self.discriminate(image_gen, Y, True)

    discrim_cost_real = match(raw_real, tf.ones_like(raw_real))
    discrim_cost_gen = match(raw_gen, tf.zeros_like(raw_gen))
    discrim_cost = discrim_cost_real + discrim_cost_gen

    gen_cost = match( raw_gen, tf.ones_like(raw_gen) )

    return Z, Y, is_training, image_real, image_gen, dis-crim_cost, gen_cost

我们现在可以开始用超参数建立图表:

n_epochs = 10
learning_rate = 0.0002
batch_size = 1024
image_shape = [28,28,1]
dim_z = 10
dim_y = 10
dim_W1 = 1024
dim_W2 = 128
dim_W3 = 64
dim_channel = 1

visualize_dim=196
step = 200

和以前一样，我们阅读了 MNIST 数据集:

from sklearn.datasets import fetch_mldata
mnist = fetch_mldata('MNIST original')
mnist.data.shape = (-1, 28, 28)
mnist.data = mnist.data.astype(np.float32).reshape( [-1, 28, 28, 1]) / 255.
mnist.num_examples = len(mnist.data)

我们创建我们的图表。我们将变量分成两个列表(因为它们有一个前缀)，每个列表都有一个 Adam 优化器。我们还为样本创建变量，以检查我们的生成器是否已经开始生成可识别的图像:

from sklearn.datasets import fetch_mldata
mnist = fetch_mldata('MNIST original')
mnist.data.shape = (-1, 28, 28)
mnist.data = mnist.data.astype(np.float32).reshape( [-1, 28, 28, 1]) / 255.
mnist.num_examples = len(mnist.data)
tf.reset_default_graph()
dcgan_model = DCGAN(
        image_shape=image_shape,
        dim_z=dim_z,
        dim_W1=dim_W1,
        dim_W2=dim_W2,
        dim_W3=dim_W3,
        )
Z_tf, Y_tf, iimage_tf, image_tf_sample, d_cost_tf, g_cost_tf, = dcgan_model.build_model()
discrim_vars = list(filter(lambda x: x.name.startswith('discriminate'),
    tf.trainable_variables()))
gen_vars = list(filter(lambda x: x.name.startswith('generate'), tf.trainable_variables()))

train_op_discrim = tf.train.AdamOptimizer(learning_rate, beta1=0.5).minimize(
    d_cost_tf, var_list=discrim_vars)
train_op_gen = tf.train.AdamOptimizer(learning_rate, beta1=0.5).minimize(
    g_cost_tf, var_list=gen_vars)

Z_np_sample = np.random.uniform(-1, 1, size=(visualize_dim,dim_z))
Y_np_sample = np.random.randint(10, size=[visualize_dim])

在我们的会话中，我们首先使用 Tensorflow 将标签转换为单向编码，然后使用与之前网络相同的模式。我们为每批生成的图像生成随机数，然后按照计划优化鉴别器和生成器:

with tf.Session() as sess:
    mnist.target = tf.one_hot(mnist.target.astype(np.int8), dim_y).eval()
    Y_np_sample = tf.one_hot(Y_np_sample, dim_y).eval()

    sess.run(tf.global_variables_initializer())
    for epoch in range(n_epochs):
        permut = np.random.permutation(mnist.num_examples)
        trX = mnist.data[permut]
        trY = mnist.target[permut]
        Z = np.random.uniform(-1, 1, size=[mnist.num_examples, dim_z]).astype(np.float32)

        print("epoch: %i" % epoch)
        for j in range(0, mnist.num_examples, batch_size):
            if j % step == 0:
                print(" batch: %i" % j)

            batch = permut[j:j+batch_size]

            Xs = trX[batch]
            Ys = trY[batch]
            Zs = Z[batch]

            sess.run(train_op_discrim,
                    feed_dict={
                        Z_tf:Zs,
                        Y_tf:Ys,
                        image_tf:Xs,
                        })

            sess.run(train_op_gen,
                    feed_dict={
                        Z_tf:Zs,
                        Y_tf:Ys,
                        })

            if j % step == 0:
                generated_samples = sess.run(
                        image_tf_sample,
                        feed_dict={
                            Z_tf:Z_np_sample,
                            Y_tf:Y_np_sample,
                            })
                generated_samples = generated_samples * 255
                save_visualization(generated_samples, (7,28),
                    save_path='./sample_%03d_%04d.jpg' %
                        (epoch, j / step))
epoch: 0
 batch: 0


…
epoch: 3
 batch: 0


epoch: 9
 batch: 64000

我们可以看到非常早期的类似手指的形状。它们的进化方式非常有趣。它们从光滑到松脆再到嘈杂(对于大量的时代)。这是可以理解的，因为这些网络没有融合。因为它们是对抗性的，每次一个人学会了对另一个人有效的技巧，另一个人就会反击。例如，如果图像不够平滑，鉴别器可以对这个差异进行鉴别，如果生成器生成平滑图像，则鉴别器将继续处理其他差异。问题是发电机会忘记过去已知的把戏，所以没有办法停在有意义的点上！

摘要

我们学习了在机器学习环境中处理图像的经典的基于特征的方法；通过将一百万个像素转换成几个数字特征，我们能够直接使用逻辑回归分类器。我们在其他章节中学到的所有技术突然变得直接适用于图像问题。我们看到了一个使用图像特征在数据集中查找相似图像的例子。

我们还学习了如何使用单词包模型中的局部特征进行分类。这是一种非常现代的计算机视觉方法，并取得了良好的效果，同时对于图像的许多不相关方面(例如照明，甚至同一图像中的不均匀照明)足够鲁棒。我们还将聚类作为分类中一个有用的中间步骤，而不是目的本身。

我们关注的是mahotas，这是 Python 中主要的计算机视觉库之一。还有一些同样得到很好的维护。Skimage 在精神上是相似的，但是有一套不同的特征。OpenCV 是一个非常好的带有 Python 接口的 C++库。所有这些都可以与 NumPy 数组一起工作，您可以混合和匹配不同库中的函数来构建复杂的计算机视觉管道。

我们还尝试了一种用 Tensorflow(可用于非图像域)生成类似图像的新方法，目前流行的网络类型名为 GAN。

在第 13 章、强化学习中，我们将探讨强化学习这一深度学习的热门话题。我们将看到如何让神经网络在不告诉它任何事情的情况下学习一套规则。

avatar
文章
阅读
粉丝