Python 应用无监督学习（三）-优快云博客

原文：annas-archive.org/md5/6b15c463e64a9f03f0d968a77b424918

译者：飞龙

协议：CC BY-NC-SA 4.0

第八章：市场篮子分析

学习目标

到本章结束时，你将能够：

使用交易级别的数据
在合适的情境中使用市场篮子分析
运行 Apriori 算法并构建关联规则
对关联规则进行基本的可视化
解读市场篮子分析的关键指标

在本章中，我们将探讨一种基础且可靠的算法，用于分析交易数据。

介绍

在本章中，我们将完全改变方向。前一章探讨了主题模型，重点是自然语言处理、文本数据和应用相对较新开发的算法。大多数数据科学从业者会同意，自然语言处理（包括主题模型）处于数据科学的前沿，是一个活跃的研究领域。我们现在已经明白，主题模型可以并且应该在任何文本数据可能带来洞察或增长的地方应用，包括社交媒体分析、推荐引擎和新闻过滤。

本章将带我们进入零售领域，探索一种基础且可靠的算法，用于分析交易数据。虽然这个算法可能不是最前沿的，也不是最流行的机器学习算法之一，但它在零售领域无处不在，并且其影响力不可否认。它带来的洞察易于解读，立刻可以采取行动，并且对确定分析的下一步非常有帮助。如果你从事零售或交易数据分析工作，那么深入了解市场篮子分析将对你非常有益。

市场篮子分析

假设你为一家零售商工作，卖着几十种产品，你的老板走过来，问你以下问题：

哪些产品最常一起购买？
产品应该如何在商店中组织和定位？
我们如何识别最适合通过优惠券打折的产品？

你可能会感到完全困惑，因为这些问题非常广泛，而且似乎无法通过单一的算法和数据集直接回答。然而，所有这些问题以及更多问题的答案就是市场篮子分析。市场篮子分析背后的基本理念是识别并量化哪些商品或商品组经常一起购买，从而为顾客行为和产品关系提供洞察。

在深入分析之前，值得定义一下“市场篮子”这个术语。市场篮子是一个经济系统中永久存在的商品集合。在这里，永久并不一定意味着传统意义上的永久。它意味着，直到该商品被从目录中移除之前，它将始终可以购买。前述定义中的商品指的是任何商品、服务或一个群体的组成部分，包括自行车、给房子刷漆或一个网站。最后，经济系统可以是一个公司、一系列活动或一个国家。市场篮子的最简单例子是杂货店，它是由一系列食品和饮料商品组成的系统。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_01.jpg

图 8.1：一个示例市场篮子，其中经济系统是肉店，永久集合的商品是肉店提供的所有肉类产品

即使不使用任何模型或分析，某些产品之间的关系也是显而易见的。让我们以肉类和蔬菜的关系为例。通常，市场篮子分析模型会返回比“肉类和蔬菜”更具体的关系，但为了讨论的方便，我们将其概括为“肉类和蔬菜”。好吧，肉类和蔬菜之间确实存在关系。那么呢？我们知道这些是常见的食材，通常会一起购买。我们可以利用这一信息，将蔬菜和肉类分别摆放在商店的两端，你会发现这两样商品常常被摆放在商店的对立面，这迫使顾客走完整个商店，从而增加他们购买额外商品的可能性，这些商品如果顾客不需要走遍整个商店，可能就不会买。

零售公司面临的一个难题是如何有效地打折商品。让我们考虑另一个显而易见的关系：花生酱和果冻。在美国，花生酱和果冻三明治非常受欢迎，尤其是在孩子们中间。当购物篮里有花生酱时，可以假设果冻也很可能在其中。既然我们知道花生酱和果冻是一起购买的，那么同时对两者打折就没有意义。如果我们希望顾客购买这两件商品，我们只需打折其中一件商品，知道如果我们能让顾客购买打折商品，他们很可能也会购买另一件商品，即使它是原价。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_02.jpg

图 8.2：市场篮子分析的可视化

就像前一章中的主题模型一样，市场篮子分析的关键是识别频繁出现的组。在这里，我们寻找的是频繁出现的产品组，而在主题模型中，我们寻找的是频繁出现的词组。因此，正如它可以应用于主题模型一样，词汇聚类也可以应用于市场篮子分析。主要的不同之处在于，市场篮子分析中的聚类是微观的，每个聚类只有少数几个产品，并且在计算概率指标时，聚类中项目的顺序是至关重要的。我们将在本章后面深入探讨这些指标以及它们是如何计算的。

从前面两个例子中可以明显看出，在市场篮子分析中，零售商能够发现顾客购买的产品之间的关系——这些关系有时显而易见，有时又出乎意料。一旦这些关系被揭示出来，就可以用来指导和改善决策过程。市场篮子分析的一个重要特点是，尽管这种分析最初是在零售领域开发、讨论并应用的，但它同样可以应用于许多不同类型的企业。

执行这种分析的唯一要求是数据必须是一个项集合的列表。在零售案例中，这通常是一个包含多次交易的列表，每次交易中包含一组已购买的产品。另一个替代应用的例子是分析网站流量。对于网站流量，我们把网站视作产品，所以列表中的每个元素就是某个个体在特定时间段内访问的所有网站集合。不用说，市场篮子分析的应用远远超出了零售领域的主应用。

用例

在传统零售应用中，有三个主要的用例：定价优化、优惠券和折扣推荐以及商店布局。如前所述，零售商可以利用模型揭示的产品关联，策略性地在商店内摆放商品，从而促使顾客购买更多商品，并因此花费更多的钱。如果两个或更多产品之间的关系足够强大——即该产品组合在数据集中出现的频率很高，并且组合中的单个产品在其他时候很少单独出现——那么这些产品就可以放在商店的远离彼此的地方，而不会显著影响顾客购买这两种产品的几率。通过迫使顾客走遍整个商店去购买这两种产品，零售商增加了顾客注意到并购买其他商品的机会。同样，零售商也可以通过将两种相关性较弱或非基础性产品放在一起，提高顾客购买这两种商品的几率。显然，商店布局受许多因素的影响，但市场篮子分析无疑是其中一个重要因素：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_03.jpg

图 8.3：产品关联如何帮助设计高效且有利可图的商店布局

定价提升和优惠券及折扣推荐是同一个问题的两个方面。可以简单地理解为在哪里提高价格，在哪里降低价格。考虑两个强相关商品的情况。这两件商品很可能会在同一笔交易中被购买，因此提高其中一件商品的价格是增加交易利润的一种方式。如果这两件商品之间的关联足够强，价格的提高几乎不会影响客户购买这两件商品的可能性。以类似的方式，零售商可以通过折扣或优惠券促使客户购买与其他商品有弱关联的商品。

例如，零售商可以将单个客户的购买历史与所有交易的市场篮分析结果进行比较，找出某些客户购买的商品与他们未购买的商品之间的弱关联。通过这个比较，零售商可以为这些客户提供折扣，推荐模型认为与他们之前购买的商品相关的尚未购买的商品。如果你曾在交易结束时收到打印出来的优惠券，极有可能这些商品与刚刚完成的交易中的商品是相关的。

市场篮分析的一个非传统但可行的应用是增强在线广告和搜索引擎优化。假设我们可以访问个人访问的网站列表。利用市场篮分析，我们可以找到网站之间的关系，并利用这些关系来策略性地排序和分组搜索引擎查询结果中的网站。在很多方面，这与商店布局的应用场景类似。

通过对市场篮分析的总体了解和对其应用场景的清晰认识，我们现在可以深入研究这些模型中使用的数据。

重要的概率度量

市场篮分析是建立在几个概率度量的计算基础上的。这里讨论的五个主要度量是支持度、置信度、提升度、杠杆度和确信度。在深入研究交易数据和具体的市场篮分析模型（包括Apriori 算法和关联规则）之前，我们应该花些时间定义并探讨这些度量，使用一个小的、虚构的交易数据集来说明。我们从编造一些数据开始。

练习 39：创建样本交易数据

由于这是本章的第一个练习，让我们设置环境。本章将使用与第七章、主题建模中相同的环境要求。如果任何一个包没有加载，就像前一章那样，使用pip通过命令行安装它们。我们将使用的一个库是mlxtend，它可能对你来说比较陌生。它是一个机器学习扩展库，包含了许多有用的辅助工具，包括集成、堆叠和市场篮分析模型。本次练习没有实际输出，我们将简单地创建一个示例交易数据集，用于后续的练习。

打开一个使用 Python 3 的 Jupyter 笔记本。
安装以下库：matplotlib.pyplot，用于绘制模型的结果，mlxtend.frequent_patterns，用于运行模型，mlxtend.preprocessing，用于对数据进行编码和准备，以适应模型，numpy，用于处理数组，pandas，用于处理数据框：

注意
```
import matplotlib.pyplot as plt
import mlxtend.frequent_patterns
import mlxtend.preprocessing
import numpy
import pandas
```

创建 10 个虚拟交易，内容为杂货店项目。这些数据将以列表的形式出现，这种数据结构在后面讨论格式化交易数据以适应模型时会非常有用：

example = [
    ['milk', 'bread', 'apples', 'cereal', 'jelly', 
     'cookies', 'salad', 'tomatoes'],
    ['beer', 'milk', 'chips', 'salsa', 'grapes', 
     'wine', 'potatoes', 'eggs', 'carrots'],
    ['diapers', 'baby formula', 'milk', 'bread', 
     'chicken', 'asparagus', 'cookies'],
    ['milk', 'cookies', 'chicken', 'asparagus', 
     'broccoli', 'cereal', 'orange juice'],
    ['steak', 'asparagus', 'broccoli', 'chips', 
     'salsa', 'ketchup', 'potatoes', 'salad'],
    ['beer', 'salsa', 'asparagus', 'wine', 'cheese', 
     'crackers', 'strawberries', 'cookies'],
    ['chocolate cake', 'strawberries', 'wine', 'cheese', 
     'beer', 'milk', 'orange juice'],
    ['chicken', 'peas', 'broccoli', 'milk', 'bread', 
     'eggs', 'potatoes', 'ketchup', 'crackers'],
    ['eggs', 'bread', 'cheese', 'turkey', 'salad', 
     'tomatoes', 'wine', 'steak', 'carrots'],
    ['bread', 'milk', 'tomatoes', 'cereal', 'chicken', 
     'turkey', 'chips', 'salsa', 'diapers']
]

这个简单的数据集将使得解释和理解概率度量变得更加容易。

支持度

支持度简单来说就是项目集在数据中出现的概率，可以通过计算项目集出现的交易次数并将该次数除以总交易数来得到。需要注意的是，项目集可以是一个单独的项目，也可以是一组项目。尽管支持度非常简单，但它是一个重要的度量指标，因为它是用于确定项目集之间关联的可信度和强度的主要指标之一。例如，可能有两个项目只在彼此之间出现，表明它们的关联非常强，但在一个包含 100 个交易的数据集中，只有两次出现并不令人印象深刻。因为该项目集只在 2%的交易中出现，而 2%在原始出现次数中算是很小的，因此该关联不能被视为显著，因此在决策中可能无法使用。

请注意，由于支持度是一个概率值，它的范围将在[0,1]之间。如果项目集包含两个项目，X 和 Y，且 N 为总交易数，则公式如下所示。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_04.jpg

图 8.4：支持度公式

让我们暂时回到练习 39中制作的数据，创建样本交易数据，并将商品集定义为牛奶和面包。我们可以轻松地查看这 10 个交易，并统计牛奶和面包商品集出现的次数——这是 4 次。鉴于总共有 10 个交易，牛奶和面包的支持度是 4 除以 10，即 0.4。是否足够大，这取决于数据集本身，我们将在后续部分进行讨论。

置信度

置信度度量可以通过条件概率来理解，它基本上是指在购买了产品 A 的前提下，购买产品 B 的概率。置信度通常表示为 A https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_01.png B，并表达为包含 A 的交易中同时包含 B 的比例。因此，置信度是通过将交易全集筛选为包含 A 的交易，然后计算这些交易中包含 B 的比例来得出的。与支持度类似，置信度是一个概率值，因此其范围是[0,1]。使用与支持度部分相同的变量定义，以下是置信度的公式：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_05.jpg

图 8.5：置信度公式

为了演示置信度，我们将使用啤酒和葡萄酒这两个商品。具体来说，我们来计算啤酒 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_02.png 葡萄酒的置信度。首先，我们需要找出包含啤酒的交易。有 3 个这样的交易，它们是交易 2、6 和 7。现在，在这些交易中，有多少包含葡萄酒？答案是所有的交易都包含葡萄酒。因此，啤酒 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_03.png 葡萄酒的置信度是 1。每次顾客购买了啤酒，他们也购买了葡萄酒。这可能很明显，但为了识别可操作的关联，更高的置信度值是更好的：

提升度和杠杆度

我们将同时讨论接下来的两个度量，提升度和杠杆度，因为尽管它们的计算方式不同，但都试图回答相同的问题。与置信度一样，提升度和杠杆度也表示为 A https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_04.png B。我们要回答的问题是，是否可以通过一个物品，比如 A，来推断另一个物品，比如 B？换句话说，如果一个人购买了产品 A，我们能否在一定程度上确定他们是否会购买产品 B？这些问题通过将 A 和 B 的支持度在假设 A 和 B 不独立的标准情况下与假设两者独立的情况进行比较来回答。提升度计算这两种情况的比率，因此其范围是[0, 无限]。当提升度等于 1 时，两个产品是独立的，因此在购买产品 A 时，无法得出关于产品 B 的任何结论：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_06.jpg

图 8.6：提升度公式

Leverage 计算两种情况之间的差异，因此其范围是[-1, 1]。Leverage 等于零可以解释为与 lift 等于一相同的含义：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_07.jpg

图 8.7：杠杆公式

这些指标的值衡量项目之间关系的强度和方向。如果 lift 值为 0.1，我们可以说两个项目之间的关系在负方向上很强。也就是说，可以认为当购买一个产品时，购买第二个产品的机会会减少。正相关和负相关被独立性点所分隔，正如前面所说，lift 的独立性点为 1，leverage 的独立性点为 0，而值越远离这些点，关联越强。

信念

最后要讨论的指标是信念，它比其他指标稍微不直观。信念是指在 X 和 Y 独立的情况下，X 发生但 Y 不发生的预期频率与错误预测频率的比值。错误预测频率定义为 1 减去 X 的置信度 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_05.png Y。记住，置信度可以定义为 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_06.png，这意味着 .https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_07.png。分子也可以视为 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_08.png。两者的唯一区别是分子假设 X 和 Y 之间是独立的，而分母则没有。理想情况下，值大于 1，因为这意味着当 X 和 Y 之间的关联是随机偶然（换句话说，X 和 Y 是独立的）时，产品或项目集合 X 和 Y 之间的关联更常是错误的。再强调一遍，这表明 X 和 Y 之间的关联是有意义的。值为 1 表示独立性，而小于 1 的值则意味着 X 和 Y 之间的随机关系比定义为 X https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_09.png Y 的关系更常见。在这种情况下，关系可能是反向的（换句话说，Y https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_10.png X）。信念的范围是[0, ∞]，其形式如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_08.jpg

图 8.8：信念公式

让我们再次回到啤酒和葡萄酒这两个产品，但为了说明本次情况，我们将考虑葡萄酒 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_11.png 啤酒的相反关联。Support(Y)，或在本例中为 Support(啤酒)，是 3/10，而 Confidence X https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_12.png Y，或在本例中为 Confidence(葡萄酒 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_13.png 啤酒)，是 3/4。因此，Conviction(葡萄酒 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_14.png 啤酒)为(1-3/10) / (1-3/4) = (7/10) * (4/1)。我们可以得出结论，如果葡萄酒和啤酒是独立的，那么葡萄酒 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_15.png 啤酒的关联会错误出现 2.8 次。因此，之前提到的葡萄酒和啤酒的关联是有效的。

练习 40：计算指标

在本练习中，我们使用练习 39中的虚拟数据，创建样本事务数据，来计算之前描述的五个指标，我们将在讲解 Apriori 算法和关联规则时再次使用这些指标。我们将评估的关联是牛奶 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_16.png 面包。

注意

本章中的所有练习都需要在同一个 Jupyter 笔记本中完成。

定义并打印构成所有五个指标基础的频率，即 Frequency(牛奶)、Frequency(面包)和 Frequency(牛奶, 面包)。还需定义 N 为数据集中交易的总数：

N = len(example)
f_x = sum(['milk' in i for i in example]) # milk
f_y = sum(['bread' in i for i in example]) # bread
f_x_y = sum([
    all(w in i for w in ['milk', 'bread']) 
    for i in example
])
print(
    "N = {}\n".format(N) + 
    "Freq(x) = {}\n".format(f_x) + 
    "Freq(y) = {}\n".format(f_y) + 
    "Freq(x, y) = {}".format(f_x_y)
)

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_09.jpg

图 8.9：频率截图

计算并打印 Support(牛奶 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_17.png 面包)：
```
support = f_x_y / N
print("Support = {}".format(round(support, 4)))
```
x 到 y 的支持度为 0.4。根据经验，如果我们使用的是完整的交易数据集，那么这个支持值在许多情况下会被认为是非常大的。
计算并打印 Confidence(牛奶 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_17.png 面包)：
```
confidence = support / (f_x / N)
print("Confidence = {}".format(round(confidence, 4)))
```
x 到 y 的置信度为 0.5714。这意味着，给定已经购买了x，Y 被购买的概率略高于 50%。
计算并打印 Lift(牛奶 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_19.png 面包)：
```
lift = confidence / (f_y / N)
print("Lift = {}".format(round(lift, 4)))
```
x 到 y 的提升度为 1.1429。
计算并打印 Leverage(牛奶 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_19.png 面包)：
```
leverage = support - ((f_x / N) * (f_y / N))
print("Leverage = {}".format(round(leverage, 4)))
```
x 到 y 的杠杆度为 0.05。提升度和杠杆度都可以用来说明关联x到y是正向的（换句话说，x意味着y），但弱。也就是说，值分别接近 1 和 0。
计算并打印 Conviction(牛奶 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_19.png 面包)：
```
conviction = (1 - (f_y / N)) / (1 - confidence)
print("Conviction = {}".format(round(conviction, 4)))
```
1.1667的置信度值可以解释为，如果牛奶和面包是独立的，那么牛奶 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_19.png 面包的关联会错误出现1.1667次。

在深入了解 Apriori 算法和实际数据上的关联规则学习之前，我们将先探索事务数据，并加载和准备一些零售数据进行建模。

事务数据的特点

在市场篮分析中使用的数据是交易数据，或者任何类似交易数据的数据。最基本的交易数据包含某种交易标识符，如发票号或交易号，以及与该标识符相关的产品列表。恰好这两项基本要素就是进行市场篮分析所需的全部内容。然而，交易数据很少——甚至可以说从未——以这种基本形式存在。交易数据通常还包括定价信息、日期和时间、客户标识符等许多其他信息：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_10.jpg

图 8.10：每个可用的产品将映射到多个发票号码

由于交易数据的复杂性，数据清洗至关重要。在市场篮分析的背景下，数据清洗的目标是过滤掉所有不必要的信息，包括移除与分析无关的变量，并清理掉有问题的交易。完成这两步清洗所使用的技术有所不同，具体取决于交易数据文件的情况。为了避免在数据清洗中陷入困境，接下来的练习将使用来自 UCI 机器学习库的在线零售数据集的一个子集，活动将使用完整的数据集。这既限制了数据清洗的讨论，又为我们提供了一个机会，讨论当数据集大小变化时，结果如何变化。这一点很重要，因为如果你为零售商工作并进行市场篮分析，你需要理解并清楚地说明，随着数据量的增加，产品之间的关系可能会发生变化，而且很可能会发生变化。在讨论此数据集所需的具体清洗过程之前，让我们先加载在线零售数据集。

练习 41：加载数据

在本次练习中，我们将加载并查看一个示例的在线零售数据集。该数据集最初来自 UCI 机器学习库，可以在github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-Python/tree/master/Lesson08/Exercise39-Exercise45找到。下载数据集后，请保存并记下路径。现在，我们开始练习。本练习的输出结果是交易数据，未来建模练习中将使用这些数据，并通过一些探索性图形帮助我们更好地理解我们正在处理的数据。

注意

数据集来自archive.ics.uci.edu/ml/datasets/online+retail#。可以从github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-Python/tree/master/Lesson08/Exercise39-Exercise45下载。陈大庆、梁赛、郭坤，《面向在线零售行业的数据挖掘：基于 RFM 模型的客户细分案例研究》，《数据库营销与客户战略管理杂志》，第 19 卷，第 3 期，197-208 页，2012 年。

UCI 机器学习库 [http://archive.ics.uci.edu/ml]。加利福尼亚州欧文市：加利福尼亚大学信息与计算机科学学院。

使用pandas中的read_excel函数加载数据。请注意，Excel 文件的第一行包含列名：
```
online = pandas.read_excel(
    io="~/Desktop/Online Retail.xlsx", 
    sheet_name="Online Retail", 
    header=0
)
```
注意

Online Retail.xlsx的路径应根据文件在系统中的位置进行修改。
打印出数据框的前 10 行。请注意，数据中包含一些与市场篮子分析无关的列：
```
online.head(10)
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_11.jpg

图 8.11：原始在线零售数据
打印出数据框中每一列的类型。此信息在执行特定的清理任务时将非常有用：
```
online.dtypes
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_12.jpg

图 8.12：数据集中每列的数据类型

获取数据框的维度，以及唯一发票号和客户标识的数量：

print(
     "Data dimension (row count, col count): {dim}"
     .format(dim=online.shape)
)
print(
     "Count of unique invoice numbers: {cnt}"
     .format(cnt=online.InvoiceNo.nunique())
)
print(
     "Count of unique customer ids: {cnt}"
     .format(cnt=online.CustomerID.nunique())
)

输出结果如下：

Data dimension (row count, col count): (541909, 8)
Count of unique invoice numbers: 25900
Count of unique customer ids: 4372

在本次练习中，我们已经加载了数据并进行了初步探索性工作。

数据清理与格式化

现在数据集已经加载，让我们深入探讨具体的数据清理过程。由于我们将数据筛选为仅包含发票号和商品项，我们将数据清理工作集中在数据集的这两列上。请记住，市场篮分析旨在识别所有客户在一段时间内购买商品之间的关联。因此，数据清理的主要工作是删除包含非正数商品的交易。这种情况可能发生在交易作废、商品退货或行政操作等场景中。这些类型的交易将通过两种方式进行筛选。首先，已取消的交易其发票号前面带有“C”，所以我们将识别这些特定的发票号并将其从数据中删除。另一种方法是删除所有商品数量为零或负数的交易。执行这两步后，数据将仅保留发票号和商品描述这两列，并删除任何包含至少一个缺失值的行。

数据清理的下一阶段涉及将数据转换为适合建模的格式。在此及后续的练习中，我们将使用完整数据的一个子集。这个子集通过获取前 5000 个唯一的发票号来创建。将数据缩减到前 5000 个唯一发票号后，我们将数据结构更改为运行模型所需的格式。请注意，目前数据是长格式，每个商品占一行。期望的格式是一个列表的列表，类似于本章早些时候所提到的虚构数据。每个子集列表表示一个唯一的发票号，因此，在本例中，外部列表应包含 5000 个子列表。子列表的元素是所有属于该发票号的商品。按照描述的清理过程，我们接下来开始练习。

练习 42：数据清理与格式化

在这个练习中，我们将执行之前描述的数据清理步骤。在处理过程中，我们将通过打印出数据的当前状态并计算一些基本的汇总指标来监控数据的变化。确保在加载数据的同一个笔记本中执行数据清理。

创建一个指示列，标明发票号是否以 “C” 开头：

online['IsCPresent'] = (
    online['InvoiceNo']
    .astype(str)
    .apply(lambda x: 1 if x.find('C') != -1 else 0)
)

筛选出所有商品数量为零或负数的交易，使用第一步创建的列删除所有以“C”开头的发票号，将 DataFrame 子集化为 InvoiceNo 和 Description，最后删除所有包含至少一个缺失值的行。将 DataFrame 重命名为 online1：

online1 = (
    online
    # filter out non-positive quantity values
    .loc[online["Quantity"] > 0]
    # remove InvoiceNos starting with C
    .loc[online['IsCPresent'] != 1]
    # column filtering
    .loc[:, ["InvoiceNo", "Description"]]
    # dropping all rows with at least one missing value
    .dropna()
)

打印出过滤后的 DataFrame online1 的前 10 行：
```
online1.head(10)
```
https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_13.jpg

图 8.13：清理后的在线零售数据集

打印出清理后的 DataFrame 的维度以及唯一发票号的数量：

print(
    "Data dimension (row count, col count): {dim}"
    .format(dim=online1.shape)
)
print(
    "Count of unique invoice numbers: {cnt}"
    .format(cnt=online1.InvoiceNo.nunique())
)

输出结果如下：

Data dimension (row count, col count): (530693, 2)
Count of unique invoice numbers: 20136

请注意，我们已经删除了大约 10,000 行和 5,800 个发票号。

将发票号从 DataFrame 中提取为列表。删除重复元素，生成唯一发票号的列表。通过打印唯一发票号列表的长度来确认处理是否成功。与步骤 4的输出进行比较：

invoice_no_list = online1.InvoiceNo.tolist()
invoice_no_list = list(set(invoice_no_list))
print(
    "Length of list of invoice numbers: {ln}"
    .format(ln=len(invoice_no_list))
)

输出结果如下：

Length of list of invoice numbers: 20136

从第五步的列表中取出，仅保留前 5,000 个元素。打印新列表的长度以确认它确实是预期的 5,000 长度：

subset_invoice_no_list = invoice_no_list[0:5000]
print(
    "Length of subset list of invoice numbers: {ln}"
    .format(ln=len(subset_invoice_no_list))
)

输出结果如下：

Length of subset list of invoice numbers: 5000

仅保留前一步列表中的发票号，过滤online1 DataFrame：

online1 = online1.loc[online1["InvoiceNo"].isin(subset_invoice_no_list)]

打印出online1的前 10 行：
```
online1.head(10)
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_14.jpg

图 8.14：仅包含 5,000 个唯一发票号的清理后的数据集

打印出 DataFrame 的维度和唯一发票号的数量，以确认过滤和清理过程是否成功：

print(
    "Data dimension (row count, col count): {dim}"
    .format(dim=online1.shape)
)
print(
    "Count of unique invoice numbers: {cnt}"
    .format(cnt=online1.InvoiceNo.nunique())
)

输出结果如下：

Data dimension (row count, col count): (129815, 2)
Count of unique invoice numbers: 5000

将online1中的数据转换为上述所述的列表形式，称为invoice_item_list。实现这一过程的方法是遍历唯一的发票号，在每次迭代时提取项目描述作为一个列表，并将该列表附加到更大的invoice_item_list列表中。打印列表中的前四个元素：

invoice_item_list = []
for num in list(set(online1.InvoiceNo.tolist())):
    # filter dataset down to one invoice number
    tmp_df = online1.loc[online1['InvoiceNo'] == num]
    # extract item descriptions and convert to list
    tmp_items = tmp_df.Description.tolist()
    # append list invoice_item_list
    invoice_item_list.append(tmp_items)

print(invoice_item_list[1:5])

输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_15.jpg

图 8.15：列表中包含四个元素的子列表，每个子列表包含属于单个发票的所有项目

注意

此步骤可能需要几分钟才能完成。

数据编码

虽然清理数据至关重要，但数据准备过程最重要的部分是将数据调整为正确的形式。在运行模型之前，当前以列表的形式存在的数据需要被编码并重新转换为 DataFrame。为此，我们将利用mlxtend的preprocessing模块中的TransactionEncoder。编码器输出的是一个多维数组，每一行的长度等于事务数据集中唯一项目的总数，元素为布尔变量，指示该特定项目是否与该行所表示的发票号相关联。数据编码后，我们可以将其重新转换为 DataFrame，其中行是发票号，列是事务数据集中的唯一项目。

在以下练习中，数据编码将使用 mlxtend 完成，但也可以不使用任何包进行编码，方法非常简单。第一步是将嵌套列表展平，返回一个包含原始嵌套列表中每个值的单一列表。接下来，去除重复的产品，并且如果需要，可以按字母顺序对数据进行排序。在进行实际编码之前，我们通过将所有元素初始化为 false 来初始化最终的 DataFrame，行数等于数据集中发票号码的数量，列名为非重复的产品名称列表。

在这种情况下，我们有 5,000 笔交易和超过 3,100 个唯一产品。因此，DataFrame 中包含超过 15,000,000 个元素。实际的编码是通过遍历每笔交易和每笔交易中的每个商品来完成的。如果交易包含某个产品，便将初始化数据集中第 i 行和第 j 列的值从 false 改为 true。由于我们需要遍历 15,000,000 个单元格，这个双重循环并不高效。虽然有一些方法可以提高性能，包括在 mlxtend 中实现的一些方法，但为了更好地理解这个过程，通过双重循环的方法是很有帮助的。以下是一个示例函数，用于在不借助任何包（除了 pandas）的情况下从头开始进行编码：

def manual_encoding(ll):
    # unlist the list of lists input
    # result is one list with all the elements of the sublists
    list_dup_unsort_items = [element for sub in ll for element in sub]
    # two cleaning steps:
    #     1\. remove duplicate items, only want one of each item in list
    #     2\. sort items in alphabetical order
    list_nondup_sort_items = sorted(list(set(list_dup_unsort_items)))

    # initialize DataFrame with all elements having False value
    # name the columns the elements of list_dup_unsort_items
    manual_df = pandas.DataFrame(
        False, 
        index=range(len(ll)), 
        columns=list_dup_unsort_items
    )

    # change False to True if element is in individual transaction list
    # each row is represents the contains of an individual transaction
    # (sublist from the original list of lists)
    for i in range(len(ll)):
        for j in ll[i]:
            manual_df.loc[i, j] = True

    # return the True/False DataFrame
    return manual_df

练习 43：数据编码

在本练习中，我们通过对上一练习中生成的嵌套列表进行编码，继续数据准备过程，以便以特定方式运行模型。

初始化并拟合事务编码器。打印出结果数据的示例：
```
online_encoder = mlxtend.preprocessing.TransactionEncoder()
online_encoder_array = online_encoder.fit_transform(invoice_item_list)
print(online_encoder_array)
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_16.jpg

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_16.jpg)

图 8.16：包含表示每笔交易中产品存在的布尔变量的多维数组
将编码后的数组重构为名为 online_encoder_df 的 DataFrame。打印出一个预定义的 DataFrame 子集，包含 true 和 false 值：
```
online_encoder_df = pandas.DataFrame(
    online_encoder_array, 
    columns=online_encoder.columns_
)
# this is a very big table, so for more 
# easy viewing only a subset is printed
online_encoder_df.loc[
    4970:4979, 
    online_encoder_df.columns.tolist()[0:8]
]
```
输出结果将类似于以下内容：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_17.jpg

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_17.jpg)

图 8.17：将编码数据的小部分重构为 DataFrame
打印出编码后 DataFrame 的维度。它应该有 5,000 行，因为生成它的数据已经过滤为 5,000 个唯一的发票号码：
```
print(
    "Data dimension (row count, col count): {dim}"
    .format(dim=online_encoder_df.shape)
)
```
输出结果将类似于以下内容：
```
Data dimension (row count, col count): (5000, 3334)
```

数据现已准备好用于建模。在接下来的部分，我们将探索 Apriori 算法。

活动 18：加载并准备完整的在线零售数据

在本活动中，我们的任务是加载并准备一个大型事务数据集进行建模。最终输出将是一个适当编码的数据集，每个独特的事务占一行，每个独特的商品占一列。如果某个商品出现在某个事务中，那么该数据框中的该元素将标记为真。

本活动将大致重复前几次练习，但将使用完整的在线零售数据集文件。无需执行新的下载，但你需要先前下载文件的路径。请在单独的 Jupyter 笔记本中执行此活动。

以下步骤将帮助你完成此活动：

加载在线零售数据集文件：

注

此数据集来自 archive.ics.uci.edu/ml/datasets/online+retail#。它可以从 https://github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-Python/tree/master/Lesson08/Activity18-Activity20 下载。Daqing Chen、Sai Liang Sain 和 Kun Guo, 数据挖掘在在线零售行业中的应用：基于 RFM 模型的客户细分案例研究，发表于《数据库营销与客户战略管理杂志》，第 19 卷，第 3 期，197-208 页，2012 年。

UCI 机器学习库 [http://archive.ics.uci.edu/ml]。加利福尼亚州尔湾：加利福尼亚大学信息与计算机科学学院。
清理并准备建模数据，包括将清理后的数据转换为列表的列表。
对数据进行编码并将其重塑为数据框：

注

本活动的解决方案可以在第 366 页找到。

输出将类似于以下内容：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_18.jpg

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_18.jpg)

图 8.18：从完整的在线零售数据集中清理、编码和重塑后的数据框的一个子集

Apriori 算法

Apriori 算法是一种数据挖掘方法，用于识别和量化事务数据中的频繁项集，是关联规则学习的基础组件。在下一节中，将讨论如何将 Apriori 算法的结果扩展到关联规则学习。Apriori 算法中作为频繁项集的最小值是模型的输入，因此是可调节的。频率在此通过支持度来量化，因此输入模型的值是分析中接受的最小支持度。模型随后识别所有支持度大于或等于输入给定的最小支持度的项集。请注意，最小支持度参数不是通过网格搜索可以优化的参数，因为 Apriori 算法没有评估指标。相反，最小支持度参数是根据数据、使用案例和领域专业知识来设置的。

Apriori 算法背后的主要思想是 Apriori 原则：任何频繁项集的子集必须本身也是频繁的。

另一个值得提及的方面是推论：不频繁项集的超集不可能是频繁的。

让我们举一些例子。如果项集 {锤子、锯子和钉子} 是频繁的，那么根据 Apriori 原则以及显而易见的道理，任何更简单的项集，例如 {锤子、锯子}，也一定是频繁的。相反，如果同样的项集 {锤子、锯子、钉子} 是不频繁的，那么增加复杂性，比如在项集 {锤子、锯子、钉子} 中加入木材 {锤子、锯子、钉子、木材}，也不会使该项集变得频繁。

计算事务数据库中每个项集的支持度值，并仅返回支持度大于或等于预设的最小支持度阈值的项集，可能看起来很简单，但实际上并非如此，因为需要进行大量计算。例如，考虑一个包含 10 个独特项的项集。这将导致 1,023 个单独的项集，需要计算它们的支持度值。现在，试着推算一下我们的工作数据集，它包含 3,135 个独特项。我们需要为这些项集计算支持度值的数量将是巨大的。计算效率是一个重大问题。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_19.jpg

图 8.19：项集如何构建以及 Apriori 原则如何大幅减少计算需求的映射（所有灰色节点为不频繁项集）。

为了解决计算需求，Apriori 算法被定义为一个自下而上的模型，包含两个步骤。这些步骤包括通过向已存在的频繁项集中添加项目来生成候选项集，并将这些候选项集与数据集进行测试，以确定这些候选项集是否也是频繁的。对于包含不频繁项集的项集，不会计算支持度值。这个过程会一直重复，直到不再有候选项集存在：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_20.jpg

图 8.20：假设最小支持度阈值为 0.4，图示显示了一般的 Apriori 算法结构。

前述结构包括建立项集、计算支持度值、过滤掉不频繁项集、创建新项集并重复此过程。

存在一个清晰的树状结构，作为识别候选项集的路径。所使用的具体搜索技术，是为遍历树状数据结构而设计的宽度优先搜索，这意味着搜索过程的每一步都专注于完全搜索树的一层，然后再移动到下一层，而不是逐分支进行搜索。

算法的高层步骤如下：

定义频繁项集。首先，这通常是单个项目的集合。
通过将频繁项集组合在一起，推导候选项集。每次增加一个项集的大小。也就是说，从一个项集的集合开始，逐渐增加到两个项集、三个项集，依此类推。
计算每个候选项集的支持度值。
创建一个新的频繁项集，由支持度值超过指定阈值的候选项集组成。

重复步骤 1到步骤 4，直到没有更多的频繁项集；也就是说，直到我们遍历了所有的组合。

Apriori 算法的伪代码如下：

L1 = {frequent items}
For k = 1 and L1 != empty set do
    Ck+1 = candidate item sets derived from Lk
    For each transaction t in the dataset do
        Increment the count of the candidates in Ck+1 that appear in t
    Compute the support for the candidates in Ck+1 using the appearance counts
    Lk+1 = the candidates in Ck+1 meeting the minimum support requirement
        End
Return L = UkLk = all frequent item sets with corresponding support values

尽管遵循 Apriori 原则，这个算法仍然可能面临显著的计算挑战，具体取决于事务数据集的大小。目前有几种策略被接受，以进一步减少计算需求。

计算修正

事务减少是一种减少计算负担的简单方法。请注意，在生成每个候选项集之后，必须扫描整个事务数据集，以统计每个候选项集的出现次数。如果我们能缩小事务数据集的大小，数据集扫描的大小将大幅减少。事务数据集的缩小通过意识到任何在第 i次迭代中不包含频繁项集的事务，在后续迭代中也不会包含频繁项集。因此，一旦每个事务不包含频繁项集，它就可以从未来扫描中使用的事务数据集中移除。

对事务数据集进行抽样并测试每个候选项集，是减少扫描事务数据集计算每个项集支持度所需的计算量的另一种方法。在实施这一方法时，重要的是要降低最小支持度要求，以确保最终数据中没有遗漏应包含的项集。由于抽样后的事务数据集会自然导致支持度值较小，因此如果将最小支持度保持在原值，将错误地从模型输出中移除那些应该是频繁项集的项集。

一种类似的方法是分区。在这种情况下，数据集被分割成若干个独立的数据集，在每个数据集上评估每个候选项集。如果某个项集在其中一个分区中频繁出现，那么它在完整的交易数据集中也被认为是频繁的。每个分区会被连续扫描，直到确定某个项集的频率。

无论是否使用这些技术，Apriori 算法的计算需求通常都会相当庞大。正如现在应该清楚的，算法的核心，支持度的计算，并不像本文讨论的其他模型那样复杂。

练习 44：执行 Apriori 算法

mlxtend使得执行 Apriori 算法变得简单。因此，本次练习将重点讲解如何操作输出的数据集以及如何解读结果。你将回忆起清洗和编码后的交易数据被定义为online_encoder_df。请在之前所有练习运行过的相同笔记本中执行本练习，因为我们将继续使用该笔记本中已经建立的环境、数据和结果。（因此，你应该使用包含 5000 条记录的缩减数据集的笔记本，而不是活动中使用的完整数据集。）

使用mlxtend运行 Apriori 算法，不改变任何默认参数值：
```
mod = mlxtend.frequent_patterns.apriori(online_encoder_df)
mod
```
输出是一个空的 DataFrame。默认的最小支持度值为 0.5，因此，由于返回了一个空的 DataFrame，我们知道所有项集的支持度都低于 0.5。根据交易的数量和可用项的多样性，没有项集支持度超过 0.5 并不罕见。
重新运行 Apriori 算法，但将最小支持度设置为 0.01。这个最小支持度值的含义是，在分析 5000 笔交易时，项集需要出现 50 次才被认为是频繁的。如前所述，最小支持度可以设置为[0,1]范围内的任何值。没有最优的最小支持度值；该值的设置完全是主观的。许多企业有自己的特定显著性阈值，但没有行业标准或优化此值的方法：
```
mod_minsupport = mlxtend.frequent_patterns.apriori(
    online_encoder_df,
    min_support=0.01
)
mod_minsupport.loc[0:6]
```
输出将类似于以下内容：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_21.jpg

图 8.21：使用 mlxtend 运行的 Apriori 算法的基本输出

注意，输出中项集的标识是数字形式的，这使得结果难以解读。
重新运行 Apriori 算法，使用与步骤 2中相同的最小支持度，但这次将use_colnames设置为 True。这样将用实际的项名称替代数字标识：
```
mod_colnames_minsupport = mlxtend.frequent_patterns.apriori(
    online_encoder_df, 
    min_support=0.01,
    use_colnames=True
)
mod_colnames_minsupport.loc[0:6]
```
输出将类似于以下内容：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_22.jpg

图 8.22：Apriori 算法输出，使用实际商品名称代替数字表示

这个 DataFrame 包含了每个支持度值大于指定最小支持度值的商品集。也就是说，这些商品集出现的频率足够高，可能具有一定的意义，因而可以采取行动。
在步骤 3的输出中添加一列，包含商品集的大小，这有助于过滤和进一步分析：
```
mod_colnames_minsupport['length'] = (
    mod_colnames_minsupport['itemsets'].apply(lambda x: len(x))
)
mod_colnames_minsupport.loc[0:6]
```
输出结果将类似于以下内容：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_23.jpg

图 8.23：Apriori 算法输出，外加一个包含商品集长度的额外列
查找包含’10 COLOUR SPACEBOY PEN’的商品集的支持度：
```
mod_colnames_minsupport[
    mod_colnames_minsupport['itemsets'] == frozenset(
        {'10 COLOUR SPACEBOY PEN'}
    )
]
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_24.jpg

图 8.24：输出的 DataFrame 被过滤为单一商品集

这个单行 DataFrame 给出了这个特定商品集的支持度值，该商品集包含一个商品。支持度值表示该商品集出现在 1.5%的交易中。
返回所有长度为 2、支持度在[0.02, 0.021]范围内的商品集
```
mod_colnames_minsupport[
    (mod_colnames_minsupport['length'] == 2) & 
    (mod_colnames_minsupport['support'] >= 0.02) &
    (mod_colnames_minsupport['support'] < 0.021)
] 
```
输出结果将类似于以下内容：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_25.jpg

图 8.25：通过长度和支持度过滤后的 Apriori 算法输出 DataFrame

这个 DataFrame 包含了所有支持度值在步骤开始时指定范围内的商品集（成对购买的商品）。这些商品集出现在约 2.0%到 2.1%的交易中。

注意，当进行support过滤时，最好指定一个范围而非具体的值，因为很有可能选择的值没有商品集。前面的输出有 18 个商品集。请记住这一点，并记住商品集中的具体商品，因为当我们扩展到完整数据时，可能会运行相同的过滤，并且我们希望执行对比。
绘制支持度值。请注意，这个图表中不会有小于 0.01 的支持度值，因为 0.01 是最小支持度值：
```
mod_colnames_minsupport.hist("support", grid=False, bins=30)
plt.title("Support")
```
输出结果将类似于以下图表：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_26.jpg

图 8.26：Apriori 算法返回的支持度值分布

最大支持度值大约为 0.14，约为 700 笔交易。看似较小的值，考虑到可用商品的数量，实际上可能并不小。商品数量较多时，通常会导致较低的支持度值，因为商品组合的变化性增加。

希望你能想到更多的方式来利用这些数据，并且从支持零售业务的角度来考虑。我们将在下一节通过使用 Apriori 算法结果来生成关联规则，从而生成更多有用的信息。

活动 19：在完整的在线零售数据集上运行 Apriori 算法

想象你在一家在线零售商工作。你被提供了上个月的所有交易数据，并被要求找出在至少 1%的交易中出现的所有项集。确定符合条件的项集后，你接着被要求识别支持度值的分布。支持度值的分布将告诉所有相关方是否存在高概率一起购买的商品组，以及支持度值的平均值。让我们为公司领导和战略家收集所有信息。

在本次活动中，你将对完整的在线零售数据集运行 Apriori 算法。

注意

该数据集来自archive.ics.uci.edu/ml/datasets/online+retail#。你可以从github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-Python/tree/master/Lesson08/Activity18-Activity20下载。Daqing Chen, Sai Liang Sain, 和 Kun Guo, 《在线零售行业的数据挖掘：基于 RFM 模型的数据挖掘客户细分的案例研究》，《数据库营销与客户战略管理期刊》，第 19 卷，第 3 期，页 197-208，2012 年。

UCI 机器学习库[http://archive.ics.uci.edu/ml]。加利福尼亚州欧文市：加利福尼亚大学信息与计算机科学学院。

确保你在与之前活动相同的笔记本中完成此活动（换句话说，使用完整数据集的笔记本，而不是使用你为练习所用的 5000 个项目的子集的笔记本）。

这也将为你提供一个机会，将结果与仅使用 5000 个交易生成的结果进行比较。这是一个有趣的活动，因为它提供了一些关于随着更多数据的收集，数据如何变化的见解，同时也为支持度值在使用分区技术时如何变化提供了一些见解。请注意，练习中的做法并不是分区技术的完美代表，因为 5000 个交易数是一个任意的抽样数。

注意

本章的所有活动都需要在同一笔记本中完成。

以下步骤将帮助你完成该活动：

在完整数据上使用合理的参数设置运行 Apriori 算法。
将结果筛选至包含10 COLOUR SPACEBOY PEN的项集。将其支持度值与练习 44，执行 Apriori 算法的结果进行比较。
添加另一列，包含项集长度。然后，筛选出那些长度为 2 且支持度在[0.02, 0.021]范围内的项集。将其与练习 44，执行 Apriori 算法的结果进行比较。
绘制support值。

注意

本次活动的解答可以在第 367 页找到。

本次活动的输出将类似于以下内容：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_27.jpg

图 8.27：支持度值的分布

关联规则

关联规则学习是一种机器学习模型，旨在发掘交易数据中隐藏的模式（换句话说，关系），这些数据描述了任何零售商的客户的购物习惯。关联规则的定义在之前定义和解释常见的概率度量时已经有所暗示。

考虑虚拟的频繁项集 {牛奶, 面包}。可以从这个项集形成两个关联规则：牛奶 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_19.png 面包和面包 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_Formula_08_19.png 牛奶。为了简便，关联规则中的第一个项集称为前件，而第二个项集称为后件。一旦关联规则被识别出来，就可以计算之前讨论的所有度量，来评估这些关联规则的有效性，确定这些规则是否可以在决策过程中发挥作用。

关联规则的建立基于支持度和置信度。支持度，正如我们在上一节中讨论的，识别哪些项目集是频繁出现的，而置信度则衡量某个特定规则的真实性频率。置信度通常被称为有趣度的度量，因为它是决定是否应该形成关联的指标。因此，建立关联规则是一个两步过程。首先识别频繁的数据集，然后评估候选关联规则的置信度，如果该置信度值超过某个任意的阈值，则该规则就成为一个关联规则。

关联规则学习的一个主要问题是发现虚假的关联，这在潜在规则数量庞大的情况下是非常可能发生的。虚假关联被定义为那些在数据中出现的规律性令人惊讶的关联，尽管这些关联完全是偶然发生的。为了清楚地表达这个观点，假设我们处于一个拥有 100 条候选规则的情境中。如果我们在 0.05 的显著性水平上进行独立性统计检验，我们仍然会面临 5%的概率，即使没有关联，仍然会发现关联。进一步假设所有的 100 条候选规则都不是有效的关联。由于 5%的概率，我们仍然会期望发现 5 条有效的关联规则。现在，将这些假设的候选规则列表规模扩大到百万或十亿级别，那么这 5%的概率就会产生一个巨大的关联数量。这个问题与几乎所有模型面临的统计显著性和错误问题类似。值得指出的是，确实存在一些技术可以用来应对虚假关联问题，但这些技术既没有在常用的关联规则库中得到一致的应用，也不在本章的讨论范围内。

现在让我们将已掌握的关联规则学习知识应用到在线零售数据集上。

练习 45：推导关联规则

在本次练习中，我们将为在线零售数据集推导关联规则并探索相关度量。确保在与之前练习相同的笔记本中完成此练习（换句话说，使用 5,000 项子集的笔记本，而不是活动中的完整数据集）。

使用mlxtend库为在线零售数据集推导关联规则。使用置信度作为有趣性度量，将最小阈值设置为 0.6，并返回所有的度量，而不仅仅是支持度。计算返回的关联规则数量：
```
rules = mlxtend.frequent_patterns.association_rules(
    mod_colnames_minsupport, 
    metric="confidence",
    min_threshold=0.6, 
    support_only=False
)
rules.loc[0:6]
```
输出结果类似于以下内容：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_28.jpg

图 8.28: 仅使用 5,000 笔交易生成的关联规则的前 7 行
按如下方式打印关联规则的数量：
```
print("Number of Associations: {}".format(rules.shape[0]))
```
找到的关联规则数量为 5,070 条。

注：

关联规则的数量可能不同。
尝试运行模型的另一个版本。选择任何最小阈值和有趣性度量。计算并探索返回的规则：
```
rules2 = mlxtend.frequent_patterns.association_rules(
    mod_colnames_minsupport, 
    metric="lift",
    min_threshold=50, 
    support_only=False
)
rules2.loc[0:6]
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_29.jpg

图 8.29: 关联规则的前 7 行
按如下方式打印关联规则的数量：
```
print("Number of Associations: {}".format(rules2.shape[0]))
```
使用 lift 度量和最小阈值为 50 时找到的关联规则数量为 26，比步骤 2中的数量明显少。我们将看到，50 是一个相当高的阈值，因此返回的关联规则较少并不令人惊讶。
将置信度与支持度作图并识别数据中的特定趋势：
```
rules.plot.scatter("support", "confidence", alpha=0.5, marker="*")
plt.xlabel("Support")
plt.ylabel("Confidence")
plt.title("Association Rules")
plt.show()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_30.jpg

图 8.30：置信度与支持度的关系图

请注意，没有任何关联规则同时具有极高的置信度和极高的支持度。这应该是可以理解的。如果一个项集有很高的支持度，那么这些项很可能会与许多其他项一起出现，这就使得置信度很高的可能性非常低。
查看置信度的分布：
```
rules.hist("confidence", grid=False, bins=30)
plt.title("Confidence")
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_31.jpg

图 8.31：置信度值的分布
现在，查看提升值的分布：
```
rules.hist("lift", grid=False, bins=30)
plt.title("Lift")
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_32.jpg

图 8.32：提升值的分布

如前所述，此图显示 50 是一个较高的阈值，因为在该值之上的点并不多。
现在，查看杠杆值的分布：
```
rules.hist("leverage", grid=False, bins=30)
plt.title("Leverage")
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_33.jpg

图 8.33：杠杆值的分布

现在，查看定罪度的分布：

plt.hist(
    rules[numpy.isfinite(rules['conviction'])].conviction.values, 
    bins = 30
)
plt.title("Conviction")

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_34.jpg

图 8.34：定罪度值的分布

四个分布中有趣的是，在图表的上端出现了不同大小的尖峰，这意味着有一些非常强的关联规则。随着置信度值增大，置信度的分布逐渐下降，但在最高值附近，分布略微上升。提升分布有最明显的尖峰。定罪度分布图在 50 左右有一个小的尖峰，也许更准确地说是一个小的突起。最后，杠杆分布在较高的值处并没有显著的尖峰，但它确实显示了一个长尾，包含一些非常高的杠杆值。

花一些时间探索模型发现的关联规则。产品配对对你有意义吗？当你改变模型参数值时，关联规则的数量发生了什么变化？你是否理解这些规则在尝试改善任何零售业务时可能产生的影响？

活动 20：在完整的在线零售数据集上查找关联规则

让我们继续在活动 19 完整在线零售数据集上的 Apriori 中设定的场景。公司领导回来说，知道每个项集在数据集中出现的频率很好，但我们可以采取哪些项集的行动？哪些项集可以用来改变商店布局或调整定价？为了找到这些答案，我们推导出了完整的关联规则。

在本次活动中，我们将从完整的在线零售交易数据集推导关联规则。确保你在使用完整数据集的笔记本中完成此活动（换句话说，就是使用完整零售数据集的笔记本，而不是练习中使用的包含 5,000 个商品子集的笔记本）。

以下步骤将帮助我们完成此活动：

在完整数据集上拟合关联规则模型。使用指标置信度，并设定最小阈值为 0.6。
计算关联规则的数量。这个数量与练习 45 的步骤 1中找到的数量是否不同？推导关联规则？
绘制置信度与支持度的关系图。
查看置信度、提升度、杠杆度和确信度的分布情况。

注意

本次活动的解决方案可以在第 370 页找到。

到本次活动结束时，你将获得关于提升度、杠杆度和确信度的图表。

总结

市场篮子分析用于分析和提取来自交易或类似交易的数据的见解，这些见解可以帮助推动多个行业的增长，最著名的就是零售行业。这些决策可能包括如何布置零售空间、打折哪些产品以及如何定价。市场篮子分析的核心支柱之一是建立关联规则。关联规则学习是一种机器学习方法，用于发现消费者购买商品之间足够强的关联，这些关联可以在商业决策中加以利用。关联规则学习依赖于 Apriori 算法，以计算高效的方式找到频繁项集。这些模型与传统的机器学习模型不同，因为它们不进行预测，结果不能通过单一指标来评估，且参数值不是通过网格搜索选择的，而是由特定问题的领域需求来决定。尽管如此，所有机器学习模型的核心目标——模式提取，在这里依然存在。在本章结束时，你应该能舒适地评估和解读概率指标，能够运行并调整使用mlxtend的 Apriori 算法和关联规则学习模型，并了解这些模型在商业中的应用。你应该知道，附近超市中的商品陈列和定价很可能是根据你和其他顾客过去的行为做出的决策！

在下一章，我们将探讨使用核密度估计的热点分析，毫无疑问，这是所有统计学和机器学习中最常用的算法之一。

第九章：热点分析

学习目标

在本章结束时，你将能够：

了解空间建模的一些应用
在适当的背景下部署热点模型
构建核密度估计模型
执行热点分析并可视化结果

在本章中，我们将学习核密度估计，并学习如何进行热点分析。

引言

让我们考虑一个假设的情景：一种新的疾病已开始在你所在的国家的多个社区传播，政府正在努力找出如何应对这一健康紧急情况。应对这一健康紧急情况的关键是流行病学知识，包括患者的位置和疾病的传播方式。能够定位和量化问题区域（通常称为热点）可以帮助卫生专业人员、政策制定者和应急响应团队制定最有效的应对策略。这个情景突出了热点建模的众多应用之一。

热点建模是一种用于识别人口在地理区域分布的方式的方法；例如，前面提到的疾病感染者在全国范围内的分布。创建这种分布依赖于代表性样本数据的可用性。请注意，人口可以是任何在地理术语上可定义的事物，包括但不限于犯罪、感染疾病的个体、具有特定人口特征的人群或飓风：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_01.jpg

图 9.1：一个虚构的火灾位置数据示例，展示了一些潜在的热点

热点分析非常流行，主要是因为它易于可视化结果，并且易于阅读和解释这些可视化图像。报纸、网站、博客和电视节目都利用热点分析来支持其中的论点、章节和话题。虽然它可能不像最流行的机器学习模型那样知名，但主要的热点分析算法——核密度估计，无疑是最广泛使用的分析技术之一。核密度估计是一种热点分析技术，用于估计特定地理事件的真实人口分布。

空间统计

空间统计学是统计学的一个分支，专注于分析具有空间特性的数据显示，包括地理或拓扑坐标。它与时间序列分析相似，目标是分析在某一维度上发生变化的数据。在时间序列分析中，数据变化的维度是时间，而在空间统计学中，数据则在空间维度上变化。空间统计学涵盖了多种技术，但我们这里关注的技术是核密度估计。

如大多数统计分析的目标一样，在空间统计学中，我们试图通过采样地理数据并利用这些数据生成见解和做出预测。地震分析是空间统计分析常见的一个应用领域。通过收集地震位置数据，可以生成标识高低地震可能性的地图，这可以帮助科学家确定未来地震发生的可能位置及强度预期。

概率密度函数

核密度估计使用概率密度函数（PDF）的概念，这是统计学中的一个基础概念。概率密度函数是一个描述连续随机变量行为的函数。也就是说，它表示随机变量取某一范围值的可能性或概率。以美国男性的身高为例，通过使用美国男性身高的概率密度函数，我们可以确定某位美国男性身高在 1.9 米到 1.95 米之间的概率：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_02.jpg

图 9.2：标准正态分布

可能是统计学中最流行的密度函数就是标准正态分布，它是以零为中心，标准差为一的正态分布。

与密度函数不同，统计学家或数据科学家通常只能获得从一个未知的总体分布中随机采集的样本值。这正是核密度估计的应用场景，它是一种利用样本数据估计随机变量的未知概率密度函数的技术：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_03.jpg

图 9.3：三种正态分布的混合

在商业中使用热点分析

我们已经提到了一些可以利用热点建模对行业产生重要影响的方式。以下是热点建模的常见应用案例。

在报告传染性疾病时，卫生组织和媒体公司通常使用热点分析来传达疾病的地理分布及其根据地理位置的传播可能性。通过使用热点分析，这些信息可以被可靠地计算并传播。热点分析非常适合处理健康数据，因为可视化图表非常直接。这意味着数据被故意或无意地误解的可能性相对较低。

热点分析还可以用于预测某些事件在地理上的发生概率。越来越多的研究领域正在利用热点分析的预测能力，其中一个例子就是环境科学领域，包括自然灾害和极端天气事件的研究。例如，地震就以难以预测而闻名，因为重大地震之间的时间间隔可能较长，而所需的用于追踪和测量地震的机械设备相对较新。

在公共政策和资源部署方面，热点分析在分析人口统计数据时可以产生重大影响。确定应该部署哪些资源（无论是金钱还是人力）可能是具有挑战性的；然而，考虑到资源往往是特定于人口的，热点分析是一种有用的技术，因为它可以用于确定某些人口统计特征的分布。这里的人口统计特征指的是我们可以找到高中毕业生、来自特定全球区域的移民或年收入超过 10 万美元的个人的地理分布。

核密度估计

热点分析的主要方法之一是核密度估计。核密度估计通过样本数据和两个被称为核函数和带宽值的参数来构建估计密度。估计的密度与任何分布一样，本质上是对随机变量行为的一个指导。在这里，我们指的是随机变量取特定值的频率，https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_Formula_01.png。在处理通常为地理数据的热点分析时，估计的密度回答了这个问题：特定的经纬度对出现的频率是多少？。如果某个特定的经纬度对，https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_Formula_02.png，以及其他附近的经纬度对出现频率较高，那么使用样本数据构建的估计密度将预示着该经纬度对周围区域的可能性较高。

核密度估计被称为一种平滑算法，因为估计密度的过程实际上是通过忽略样本数据中的异常和离群值来估计数据的潜在形状。换句话说，核密度估计去除了数据中的噪声。该模型的唯一假设是数据确实属于某个可解释且有意义的密度，能够从中提取洞察并付诸实践。也就是说，存在一个真实的潜在分布。

可以说，比本书中的任何其他话题，核密度估计体现了统计学的基本思想，即利用有限大小的样本数据对总体做出推断。我们假设样本数据包含数据点的聚类，这些聚类意味着总体中高可能性的区域。创建高质量的真实总体密度估计的好处是，估计的密度可以用来从总体中抽取更多的数据。

在这段简要介绍之后，你可能会有以下两个问题：

什么是带宽值？
什么是核函数？

我们接下来会回答这两个问题。

带宽值

核密度估计中最关键的参数被称为带宽值，它对估计结果的影响不可过高估计。带宽值的高层次定义是，它决定了平滑的程度。如果带宽值低，那么估计的密度将具有有限的平滑度，这意味着密度会捕捉到样本数据中的所有噪声。如果带宽值高，那么估计的密度将非常平滑。过于平滑的密度会移除估计密度中那些真实的且不只是噪声的特征。

在更多的统计学或机器学习语言中，带宽参数控制着偏差-方差权衡。也就是说，低带宽值会导致高方差，因为密度对样本数据的方差非常敏感。低带宽值限制了模型适应样本数据中未在总体中出现的空缺的能力。使用低带宽值估计的密度往往会出现过拟合（这也被称为过度平滑的密度）。当使用高带宽值时，所得的密度会出现欠拟合，并且估计的密度会有较高的偏差（这也被称为过度平滑的密度）。

练习 46：带宽值的影响

在本练习中，我们将拟合九种不同的模型，每个模型使用不同的带宽值，来处理练习中创建的样本数据。这里的目标是巩固我们对带宽参数影响的理解，并明确表示，如果需要准确的估计密度，则必须小心选择带宽值。请注意，找到最优带宽值将是下一部分的主题。所有练习将在 Jupyter notebook 中使用 Python 3 完成；确保所有包的安装通过pip进行。安装mpl_toolkits中的basemap模块的最简单方法是使用Anaconda。下载和安装Anaconda的说明可以在本书的开头找到：

加载本章练习所需的所有库。在这里，matplotlib库用于创建基本图形；basemap库用于创建涉及位置数据的图形；numpy库用于处理数组和矩阵；pandas库用于处理 DataFrame；scipy库用于 Python 中的科学计算；seaborn库用于创建更具吸引力和复杂的图形；sklearn库用于访问数据、处理数据和运行模型。此外，确保图形以内联方式运行并设置为seaborn，以便所有图形以seaborn图形的形式呈现：
```
get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
import mpl_toolkits.basemap
import numpy
import pandas
import scipy.stats
import seaborn
import sklearn.datasets
import sklearn.model_selection
import sklearn.neighbors
seaborn.set()
```

创建一些样本数据（vals），通过混合三个正态分布。除了样本数据外，还定义真实的密度曲线（true_density）以及数据将被绘制的范围（x_vec）：

x_vec = numpy.linspace(-30, 30, 10000)[:, numpy.newaxis]
vals = numpy.concatenate((
    numpy.random.normal(loc=1, scale=2.5, size=500), 
    numpy.random.normal(loc=10, scale=4, size=500), 
    numpy.random.normal(loc=-12, scale=5, size=500)
))[:, numpy.newaxis]
true_density = (
    (1 / 3) * scipy.stats.norm(1, 2.5).pdf(x_vec[:, 0]) + 
    (1 / 3) * scipy.stats.norm(10, 4).pdf(x_vec[:, 0]) +
    (1 / 3) * scipy.stats.norm(-12, 5).pdf(x_vec[:, 0])
)

定义一个元组列表，用于指导创建多图形。每个元组包含特定子图的行和列索引，以及用于在该特定子图中创建估计密度的带宽值：

position_bandwidth_vec = [
    (0, 0, 0.1), (0, 1, 0.4), (0, 2, 0.7), 
    (1, 0, 1.0), (1, 1, 1.3), (1, 2, 1.6), 
    (2, 0, 1.9), (2, 1, 2.5), (2, 2, 5.0)
]

创建九个图形，每个图形使用不同的带宽值。第一个图形（索引为(0, 0)）将具有最低带宽，最后一个图形（索引为(2, 2)）将具有最高带宽。这些值不是绝对最低或绝对最高的带宽值，而仅仅是前一步中定义的列表中的最小和最大值：

fig, ax = plt.subplots(3, 3, sharex=True, sharey=True, figsize=(12, 9))
fig.suptitle('The Effect of the Bandwidth Value', fontsize=16)
for r, c, b in position_bandwidth_vec:
    kde = sklearn.neighbors.KernelDensity(bandwidth=b).fit(vals)
    log_density = kde.score_samples(x_vec)
    ax[r, c].hist(vals, bins=50, density=True, alpha=0.5)
    ax[r, c].plot(x_vec[:, 0], numpy.exp(log_density), '-', linewidth=2)
    ax[r, c].set_title('Bandwidth = {}'.format(b))

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_04.jpg

图 9.4：一个 3 x 3 的子图矩阵；每个子图都使用九个带宽值之一创建了一个估计的密度

注意，在较低带宽值时，密度曲线明显过度拟合数据。随着带宽值的增加，估计的密度变得更加平滑，直到明显低估数据。视觉上看，最优带宽可能在 1.6 左右。

下一步是设计一个算法来识别最优带宽值，从而使得估计的密度最为合理，因此也是最可靠和可操作的。

选择最优带宽

如前一个练习中提到的，我们可以通过仅仅比较几种密度来接近选择最优带宽。然而，这既不是选择参数值的最有效方法，也不是最可靠的方法。

有两种标准方法可以优化带宽值，这两种方法将在未来的练习和活动中出现。第一种方法是插件法（或公式化方法），它是确定性的，并且没有在样本数据上进行优化。插件方法通常实现更快，代码更简单，解释起来也更容易。然而，这些方法有一个重大缺点，即与在样本数据上进行优化的方法相比，其准确性往往较差。这些方法还具有分布假设。最流行的插件方法是 Silverman 规则和 Scott 规则。默认情况下，seaborn包（将在未来的练习中使用）使用 Scott 规则作为确定带宽值的方法。

第二种，也是更强健的方法，是通过搜索预定义的带宽值网格来找到最优的带宽值。网格搜索是一种经验性的方法，常用于机器学习和预测建模中，以优化模型的超参数。这个过程从定义带宽网格开始，带宽网格只是待评估的带宽值的集合。使用网格中的每个带宽值来创建估计的密度；然后，使用伪对数似然值对估计的密度进行评分。最优带宽值是具有最大伪对数似然值的那个。可以把伪对数似然值看作是获取我们确实获得数据点的概率和未获取数据点的概率。理想情况下，这两个概率应该都很大。考虑一下获取我们确实获得数据点的概率较低的情况。在这种情况下，意味着样本中的数据点是异常的，因为在真实分布下，获取我们确实获得的点的概率应该不高。

现在，让我们实现网格搜索方法来优化带宽值。

练习 47：使用网格搜索选择最优带宽

在本练习中，我们将为练习 46中创建的样本数据生成估计密度，带宽值的影响，并使用网格搜索和交叉验证方法确定最优带宽值。为了进行带有交叉验证的网格搜索，我们将使用sklearn，这是本书中一直使用的工具。这个练习是练习 1 的延续，因为我们使用的是相同的样本数据，并继续探索带宽值：

定义带宽值的网格和网格搜索交叉验证模型。理想情况下，应该使用逐一剔除交叉验证方法，但为了使模型在合理的时间内运行，我们将采用 10 倍交叉验证。按如下方式拟合模型：

bandwidths = 10 ** numpy.linspace(-1, 1, 100)
grid = sklearn.model_selection.GridSearchCV(
    estimator=sklearn.neighbors.KernelDensity(kernel="gaussian"),
    param_grid={"bandwidth": bandwidths},
    cv=10 #sklearn.model_selection.LeaveOneOut().get_n_splits(vals)
)
grid.fit(vals)

从模型中提取最优带宽值，如下所示：
```
best_bandwidth = grid.best_params_["bandwidth"]
print(
    "Best Bandwidth Value: {}"
    .format(best_bandwidth)
)
```
最优带宽值应该大约为 2。我们可以将最优带宽值解释为生成最大伪对数似然值的带宽值。请注意，根据网格中包含的值，最优带宽值可能会有所变化。

绘制样本数据的直方图，叠加真实密度和估计密度。在这种情况下，估计密度将是最优估计密度：

fig, ax = plt.subplots(figsize=(14, 10))
ax.hist(vals, bins=50, density=True, alpha=0.5, label='Sampled Values')
ax.fill(
      x_vec[:, 0], true_density,
      fc='black', alpha=0.3, label='True Distribution'
      )
log_density = numpy.exp(grid.best_estimator_.score_samples(x_vec))
ax.plot(
        x_vec[:, 0], log_density,
        '-', linewidth=2, label='Kernel = Gaussian'
        )
ax.legend(loc='upper right')

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_05.jpg

图 9.5：随机样本的直方图，叠加了真实密度和最优估计密度

估计密度没有明显的过拟合或欠拟合，且它确实捕捉到了三个聚类。可以说，它可能会更好地映射到真实密度，但这只是由有局限性的模型生成的估计密度。

现在让我们继续第二个问题：什么是核函数，它在模型中扮演什么角色？

核函数

另一个需要设置的参数是核函数。核是一个非负函数，控制密度的形状。像主题模型一样，我们在一个非负环境中工作，因为负的可能性或概率是没有意义的。

核函数通过系统地加权点来控制估计密度的形状。这种加权的系统方法相当简单；与许多其他数据点接近的数据点会被加权，而孤立或远离其他数据点的数据点则会被减权。加权的数据点将对应于最终估计密度中较高可能性的点。

可以使用多种函数作为核，但六种常见选择是高斯核、顶帽核、埃潘尼基诺夫核、指数核、线性核和余弦核。这些函数分别代表独特的分布形状。请注意，在每个公式中，参数h表示带宽值：

高斯：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_06.jpg

图 9.6：高斯核函数的公式

Tophat 核：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_07.jpg

图 9.7：Tophat 核函数的公式

Epanechnikov 核：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_08.jpg

图 9.8：Epanechnikov 核函数的公式

指数型：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_09.jpg

图 9.9：指数核函数的公式

线性型：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_10.jpg

图 9.10：线性核函数的公式

余弦型：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_11.jpg

图 9.11：余弦核函数的公式

这里是六种核函数的分布形状：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_12.jpg

图 9.12：六种核函数的整体形状

核函数的选择并非完全不重要，但它绝对没有带宽值的选择那么重要。一个合理的行动方案是，对于所有的密度估计问题使用高斯核函数，这也是我们在接下来的练习和活动中将要做的。

练习 48：核函数的影响

本练习的目标是理解核函数的选择如何影响密度估计的质量。就像我们在探索带宽值的影响时一样，我们将保持其他所有参数不变，使用在前两个练习中生成的相同数据，并使用先前指定的六种核函数运行六个不同的核密度估计模型。六个估算的密度之间应该能看到明显的差异，但这些差异应该比使用不同带宽值时的密度差异稍微小一些：

定义一个与之前定义的类似的元组列表。每个元组包括子图的行列索引，以及用于创建密度估计的核函数：

position_kernel_vec = [
    (0, 0, 'gaussian'), (0, 1, 'tophat'), 
    (1, 0, 'epanechnikov'), (1, 1, 'exponential'), 
    (2, 0, 'linear'), (2, 1, 'cosine'), 
]
Fit and plot six kernel density estimation models using a different kernel function for each. To truly understand the differences between the kernel functions, we will set the bandwidth value to the optimal bandwidth value found in Exercise 2 and not adjust it:
fig, ax = plt.subplots(3, 2, sharex=True, sharey=True, figsize=(12, 9))
fig.suptitle('The Effect of Different Kernels', fontsize=16)
for r, c, k in position_kernel_vec:
    kde = sklearn.neighbors.KernelDensity(
        kernel=k, bandwidth=best_bandwidth
        ).fit(vals)
    log_density = kde.score_samples(x_vec)
    ax[r, c].hist(vals, bins=50, density=True, alpha=0.5)
    ax[r, c].plot(x_vec[:, 0], numpy.exp(log_density), '-', linewidth=2)
    ax[r, c].set_title('Kernel = {}'.format(k.capitalize()))

输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_13.jpg

图 9.13：一个 3 x 2 的子图矩阵，每个子图展示使用六种核函数之一估算的密度

在这六种核函数中，高斯核函数产生了最合理的估算密度。更进一步，注意到不同核函数下的估算密度之间的差异小于使用不同带宽值时的差异。这证明了之前的观点，即带宽值是更重要的参数，应该在模型构建过程中重点关注。

在我们大致理解的基础上，接下来让我们以高层次的方式讨论核密度估计的推导过程。

核密度估计推导

我们跳过正式的数学推导，而采用直观的流行推导方法。核密度估计将样本中的每个数据点转化为其自身的分布，其宽度由带宽值控制。然后，将这些单独的分布相加，以创建所需的密度估计。这个概念相对容易演示；然而，在接下来的练习中之前，让我们尝试以抽象的方式思考它。对于包含许多样本数据点的地理区域，单独的密度将会重叠，并且通过加总这些密度，将在估计密度中创建高概率点。同样，对于包含少量甚至没有样本数据点的地理区域，单独的密度不会重叠，因此在估计密度中对应的将是低概率点。

练习 49：模拟核密度估计的推导

这里的目标是演示将单独分布相加，以创建随机变量的总体估计密度的概念。我们将通过从一个样本数据点开始，逐步建立这个概念，随后增加更多的样本数据点。此外，还将应用不同的带宽值，从而进一步巩固我们对带宽值对这些单独密度影响的理解：

定义一个函数来评估正态分布。输入值包括代表随机变量范围的网格 X、采样数据点 m 和带宽 b：

def eval_gaussian(x, m, b):
    numerator = numpy.exp(
        -numpy.power(x - m, 2) / (2 * numpy.power(b, 2))
    )
    denominator = b * numpy.sqrt(2 * numpy.pi)
    return numerator / denominator

将单个样本数据点绘制为直方图，并与不同带宽值下的单独密度进行比较：

m = numpy.array([5.1])
b_vec = [0.1, 0.35, 0.8]
x_vec = numpy.linspace(1, 10, 100)[:, None]
fig, ax = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(15, 10))
for i, b in enumerate(b_vec):
    ax[0, i].hist(m[:], bins=1, fc='#AAAAFF', density=True)
    ax[0, i].set_title("Histogram: Normed")
    evaluation = eval_gaussian(x_vec, m=m[0], b=b)

    ax[1, i].fill(x_vec, evaluation, '-k', fc='#AAAAFF')
    ax[1, i].set_title("Gaussian Dist: b={}".format(b))

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_14.jpg

图 9.14：展示一个数据点及其在不同带宽值下的单独密度

在这里，我们看到已经建立的结论，即较小的带宽值会产生非常狭窄的密度，容易导致数据过拟合。

重新生成步骤 2中的工作，但这次扩展到 16 个数据点：

m = numpy.random.normal(4.7, 0.88, 16)
n = len(m)
b_vec = [0.1, 0.35, 1.1]
x_vec = numpy.linspace(-1, 11, 100)[:, None]
fig, ax = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(15, 10))
for i, b in enumerate(b_vec):
    ax[0, i].hist(m[:], bins=n, fc='#AAAAFF', density=True)
    ax[0, i].set_title("Histogram: Normed")

    sum_evaluation = numpy.zeros(len(x_vec))

    for j in range(n):
        evaluation = eval_gaussian(x_vec, m=m[j], b=b) / n
        sum_evaluation += evaluation[:, 0]

        ax[1, i].plot(x_vec, evaluation, '-k', linestyle="dashed")
    ax[1, i].fill(x_vec, sum_evaluation, '-k', fc='#AAAAFF')
    ax[1, i].set_title("Gaussian Dist: b={}".format(b))

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_15.jpg

图 9.15：展示 16 个数据点、它们在不同带宽值下的单独密度，以及它们单独密度的总和

如预期的那样，使用最小带宽值的图形呈现出严重过拟合的估计密度。也就是说，估计密度捕捉到了样本数据中的所有噪声。在这三种密度中，第二种密度，即带宽值设置为 0.35 时的估计密度，最为合理。

活动 21：估计一维密度

在这第一个活动中，我们将生成一些虚拟示例数据，并使用核密度估计估计密度函数。带宽值将通过网格搜索交叉验证进行优化。目标是通过在一个简单的单维案例中运行模型，巩固我们对这一有用方法论的理解。我们将再次使用 Jupyter notebooks 来完成这项工作。

假设我们将要创建的示例数据描述的是美国某州的房价。暂时忽略以下示例数据中的数值。问题是，*房价的分布是什么样的，我们能否提取出房子价格落在某个特定范围内的概率？*这些问题以及更多的问题都可以通过核密度估计来回答。

完成该活动的步骤如下：

打开一个新的笔记本，并安装所有必要的库。
从标准正态分布中采样 1,000 个数据点。将 3.5 加到样本的最后 625 个值上（即，375 到 1,000 之间的索引）。设置随机状态为 100。为此，使用numpy.random.RandomState设置一个随机状态为 100，以保证相同的采样值，然后使用randn(1000)调用随机生成数据点。
将 1,000 个样本数据绘制成直方图，并在其下方添加散点图。
定义一个带宽值的网格。然后，定义并拟合一个网格搜索交叉验证算法。
提取最佳带宽值。
重新绘制步骤 3中的直方图，并叠加估计的密度。

输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_16.jpg

图 9.16：带有最佳估计密度叠加的随机样本直方图

注意

这个活动的解决方案可以在第 374 页找到。

热点分析

首先，热点是数据点浓度较高的区域，例如犯罪率异常高的特定街区，或遭受异常多龙卷风影响的某些地区。热点分析是通过使用采样数据，在总体中寻找这些热点的过程。这个过程通常通过利用核密度估计来完成。

热点分析可以通过四个高层次步骤来描述：

收集数据：数据应包括对象或事件的位置。如我们简要提到的，进行分析并获得可操作结果所需的数据量相对灵活。理想状态是拥有一个能够代表总体的样本数据集。
确定基础地图：下一步是确定哪种基础地图最适合项目的分析和展示需求。在这张基础地图上，将叠加模型的结果，以便更容易地以城市、街区或区域等更易理解的术语表达热点的位置。
执行模型：在此步骤中，您将选择并执行一种或多种提取空间模式以识别热点的方法。对于我们来说，这种方法将是——毫无悬念——核密度估计。
创建可视化：热点图是通过将模型结果叠加在基础地图上生成的，旨在支持解决任何悬而未决的业务问题。

从可用性角度来看，热点分析的一个主要问题是热点的统计显著性并不是特别容易确定。关于统计显著性的大部分问题围绕热点的存在展开。也就是说，发生可能性的波动是否真的构成统计学上的显著波动？需要注意的是，执行核密度估计并不要求统计显著性，且我们在后续过程中将完全不涉及显著性问题。

尽管“热点”一词传统上用来描述一组地理位置数据点，但它不限于位置数据。任何类型的数据都可能有热点，无论这些数据是否被称为热点。在接下来的练习中，我们将对一些非位置数据进行建模，以找出热点，这些热点是特征空间中发生可能性较高或较低的区域。

练习 50：加载数据并使用 Seaborn 进行建模

在本次练习中，我们将使用 seaborn 库来拟合和可视化核密度估计模型。这将应用于位置数据和非位置数据。开始建模之前，我们加载数据，这些数据是与 sklearn 一起自动加载的加利福尼亚住房数据集。该数据集来源于 1990 年美国人口普查，描述了当时加利福尼亚的住房情况。数据集中的每一行描述了一个人口普查块组。人口普查块组的定义与本次练习无关，因此我们将跳过对其的定义，专注于更多的实操编码与建模。需要提到的是，所有变量都是按人口普查块进行聚合的。例如，MedInc 是每个人口普查块的家庭收入中位数。关于此数据集的更多信息，请访问 scikit-learn.org/stable/datasets/index.html#california-housing-dataset：

使用 fetch_california_housing() 加载加利福尼亚住房数据集。使用 pandas 将数据转换为 DataFrame 并打印出 DataFrame 的前五行：
```
housing = sklearn.datasets.fetch_california_housing()
df = pandas.DataFrame(housing['data'], columns=housing['feature_names'])
print("Dataframe Dimensions: {dims}".format(dims=df.shape))
df.head()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_17.jpg

图 9.17：来自 sklearn 的加利福尼亚住房数据集的前五行
根据HouseAge特征过滤数据框，该特征表示每个普查区块的房屋中位数年龄。仅保留HouseAge小于或等于 15 的行，并将数据框命名为dfLess15。打印数据框的前五行；然后，将数据框缩减为仅包含经度和纬度特征：
```
dfLess15 = df[df['HouseAge'] <= 15.0]
dfLess15 = dfLess15[['Latitude', 'Longitude']]
print(
    "Less Than Or Equal To 15 Years Dataframe Dimensions: {dims}"
    .format(dims=dfLess15.shape)
)
dfLess15.head()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_18.jpg

图 9.18：过滤后的数据集前五行，仅包含HouseAge列值小于或等于 15 的行
使用seaborn拟合并可视化基于经度和纬度数据点构建的核密度估计模型。seaborn拟合这些模型的方法使用了 Scott 规则。该模型有四个输入，它们是求估计密度的两列的名称（即经度和纬度）、这些列所属的数据框，以及密度估计的方法（即kde或核密度估计）：
```
seaborn.jointplot("Longitude", "Latitude", dfLess15, kind="kde")
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_19.jpg

图 9.19：包含二维估计密度以及 dfLess15 数据集的边际密度的联合图

如果我们将这些结果叠加到加利福尼亚州的地图上，我们将看到热点位于南加州，包括洛杉矶和圣地亚哥，湾区，包括旧金山，以及在一定程度上被称为中央谷地的地区。这个seaborn图形的一个优点是，我们可以得到二维估计密度以及经度和纬度的边际密度。
基于HouseAge特征创建另一个过滤后的数据框；这次仅保留HouseAge大于 40 的行，并将数据框命名为dfMore40。此外，移除所有列，保留经度和纬度。然后，打印数据框的前五行：
```
dfMore40 = df[df['HouseAge'] > 40.0]
dfMore40 = dfMore40[['Latitude', 'Longitude']]
print(
    "More Than 40 Years Dataframe Dimensions: {dims}"
    .format(dims=dfMore40.shape)
)
dfMore40.head()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_20.jpg

图 9.20：过滤后的数据集顶部，仅包含HouseAge列中大于 40 的行
重复步骤 3的过程，但这次使用新的过滤后的数据框：
```
seaborn.jointplot("Longitude", "Latitude", dfMore40, kind="kde")
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_21.jpg

图 9.21：包含二维估计密度以及 dfMore40 数据集的边际密度的联合图

这个估算密度要紧凑得多，因为数据几乎完全聚集在两个区域。这些区域是洛杉矶和湾区。将其与步骤 3中的图形进行比较，我们注意到住房开发已经遍布全州。此外，新建住房开发在更多普查区块中出现的频率更高。
我们再次创建一个新的过滤后的 DataFrame。这次仅保留HouseAge小于或等于 5 的行，并将该 DataFrame 命名为dfLess5。绘制Population和MedInc的散点图，方法如下：
```
dfLess5 = df[df['HouseAge'] <= 5]
x_vals = dfLess5.Population.values
y_vals = dfLess5.MedInc.values
fig = plt.figure(figsize=(10, 10))
plt.scatter(x_vals, y_vals, c='black')
plt.xlabel('Population', fontsize=18)
plt.ylabel('Median Income', fontsize=16)
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_22.jpg

图 9.22：中位收入与人口的散点图，HouseAge 列中的值为 5 或以下

使用另一个seaborn函数来拟合核密度估计模型。同样，使用 Scott 规则找到最优带宽。重新绘制直方图并叠加估算密度，如下所示：

fig = plt.figure(figsize=(10, 10))
ax = seaborn.kdeplot(
    x_vals, 
    y_vals,
    kernel='gau',
    cmap='Blues', 
    shade=True, 
    shade_lowest=False
)
plt.scatter(x_vals, y_vals, c='black', alpha=0.05)
plt.xlabel('Population', fontsize=18)
plt.ylabel('Median Income', fontsize=18)
plt.title('Density Estimation With Scatterplot Overlay', size=18)

输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_23.jpg

图 9.23：与步骤 6 中创建的散点图相同，叠加了估算密度

在这里，估算密度显示，人口较少的普查区块比人口较多的普查区块更有可能拥有较低的中位收入，而不是较高的中位收入。本步骤的目的是展示如何在非位置数据上使用核密度估计。

在展示热点分析结果时，应当使用某种地图，因为热点分析通常是基于位置数据进行的。获取可以叠加估算密度的地图并不是一个简单的过程。由于版权问题，我们将使用非常基础的地图，称为底图，来叠加我们的估算密度。如何将本章获得的知识扩展到更复杂和详细的地图将留给你自己去做。地图环境的下载和安装可能也会很复杂且耗时。

练习 51：与底图一起工作

这个练习利用了mpl_toolkits中的basemap模块。basemap是一个地图绘制库，可以用来创建基础的地图或地理区域的轮廓。这些地图可以叠加核密度估计的结果，从而清晰地看到热点的位置。

首先，通过在 Jupyter Notebook 中运行import mpl_toolkits.basemap来检查basemap是否已安装。如果加载没有错误，那么你已经准备好了，无需进一步操作。如果调用失败，则使用pip安装basemap，方法是运行python3 -m pip install basemap。在重新启动任何已打开的 Notebook 之后，你应该就可以正常使用。请注意，pip安装只有在安装了 Anaconda 的情况下才有效。

本练习的目标是重新建模并重新绘制练习 50中位置数据的图表，使用sklearn的核密度估计功能和basemap的映射能力。从名为dfLess15的筛选后的 DataFrame 中提取经纬度值，如下所示：

形成将要叠加估算密度的位置网格。位置网格是定义随机变量范围的一维向量的二维位置等价物，在练习 1 中已涉及。
```
xgrid15 = numpy.sort(list(dfLess15['Longitude']))
ygrid15 = numpy.sort(list(dfLess15['Latitude']))
x15, y15 = numpy.meshgrid(xgrid15, ygrid15)
print("X Grid Component:\n{}\n".format(x15))
print("Y Grid Component:\n{}\n".format(y15))
xy15 = numpy.vstack([y15.ravel(), x15.ravel()]).T
print("Grid:\n{}\n".format(xy15))
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_24.jpg

图 9.24：表示 dfLess15 数据集的网格的 x 和 y 分量

定义并拟合一个核密度估计模型。设置带宽值为 0.05，以节省运行时间；然后，为位置网格中的每个点创建似然值：

kde = sklearn.neighbors.KernelDensity(
    bandwidth=0.05, 
    metric='minkowski',
    kernel='gaussian', 
    algorithm='ball_tree'
)
kde.fit(dfLess15.values)
log_density = kde.score_samples(xy15)
density = numpy.exp(log_density)
density = density.reshape(x15.shape)
print("Shape of Density Values:\n{}\n".format(density.shape))

注意，如果你打印出似然值的形状，它是 3,287 行和 3,287 列，总共有 10,804,369 个似然值。这与预设的经纬度网格（称为xy15）中的值数量相同。

创建加利福尼亚州的轮廓，并叠加在步骤 2中估算的密度值上：
```
fig = plt.figure(figsize=(10, 10))
fig.suptitle(
    """
    Density Estimation:
    Location of Housing Blocks
    Where the Median Home Age <= 15 Years
    """, 
    fontsize=16
)
the_map = mpl_toolkits.basemap.Basemap(
    projection='cyl',
    llcrnrlat=y15.min(), urcrnrlat=y15.max(),
    llcrnrlon=x15.min(),urcrnrlon=x15.max(),
    resolution='c'
)
the_map.drawcoastlines(linewidth=1)
the_map.drawcountries(linewidth=1)
the_map.drawstates(linewidth=1)
levels = numpy.linspace(0, density.max(), 25)
plt.contourf(x15, y15, density, levels=levels, cmap=plt.cm.Reds)
plt.show()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_25.jpg

图 9.25：将 dfLess15 的估算密度叠加到加利福尼亚州的轮廓上

0.05 的值是故意设置为稍微过拟合数据的。你会注意到，与在练习 50，使用 Seaborn 加载数据和建模中组成密度的大型聚类不同，这里的估算密度由许多更小的聚类组成。这个稍微过拟合的密度可能比之前的版本更有帮助，因为它能更清楚地显示出高似然值的普查区块真正的位置。前一个密度中的一个高似然区域是南加州，但南加州是一个巨大的地区，拥有庞大的人口和许多市政区。请记住，在使用结果做出商业决策时，可能需要特定的精度级别，并且如果样本数据可以支持这样的结果，应该提供该精度或粒度。

重复步骤 1，但使用dfMore40 DataFrame：

xgrid40 = numpy.sort(list(dfMore40['Longitude']))
ygrid40 = numpy.sort(list(dfMore40['Latitude']))
x40, y40 = numpy.meshgrid(xgrid40, ygrid40)
print("X Grid Component:\n{}\n".format(x40))
print("Y Grid Component:\n{}\n".format(y40))
xy40 = numpy.vstack([y40.ravel(), x40.ravel()]).T
print("Grid:\n{}\n".format(xy40))

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_26.jpg

图 9.26：表示 dfMore40 数据集的网格的 x 和 y 分量

使用步骤 4中建立的网格，重复步骤 2：

kde = sklearn.neighbors.KernelDensity(
    bandwidth=0.05, 
    metric='minkowski',
    kernel='gaussian', 
    algorithm='ball_tree'
)
kde.fit(dfMore40.values)
log_density = kde.score_samples(xy40)
density = numpy.exp(log_density)
density = density.reshape(x40.shape)
print("Shape of Density Values:\n{}\n".format(density.shape))

使用步骤 3中估算的密度值，重复操作：

fig = plt.figure(figsize=(10, 10))
fig.suptitle(
    """
    Density Estimation:
    Location of Housing Blocks
    Where the Median Home Age > 40 Years
    """, 
    fontsize=16
)
the_map = mpl_toolkits.basemap.Basemap(
    projection='cyl',
    llcrnrlat=y40.min(), urcrnrlat=y40.max(),
    llcrnrlon=x40.min(),urcrnrlon=x40.max(),
    resolution='c'
)
the_map.drawcoastlines(linewidth=1)
the_map.drawcountries(linewidth=1)
the_map.drawstates(linewidth=1)
levels = numpy.linspace(0, density.max(), 25)
plt.contourf(x40, y40, density, levels=levels, cmap=plt.cm.Reds)
plt.show()

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_27.jpg

图 9.27：将 dfMore40 的估算密度叠加到加利福尼亚州的轮廓上

这个估计密度是我们在练习 50中重新做的那个，使用 Seaborn 加载数据和建模。虽然第 3 步中的密度会为对房地产或人口普查感兴趣的人提供更多细节，但这个密度实际上与练习 50中的对应密度差别不大。聚类主要集中在洛杉矶和旧金山，几乎没有在其他地方出现数据点。

活动 22：伦敦犯罪分析

在此活动中，我们将对来自data.police.uk/data/的伦敦犯罪数据进行核密度估计的热点分析。由于处理地图数据的困难，我们将使用seaborn来可视化分析结果。不过，如果你觉得有信心并且能够运行练习 51中的所有图表，与底图一起工作，那么鼓励你尝试使用地图。

对这个犯罪数据进行热点分析的动机有两个方面。我们首先需要确定某些类型的犯罪在高概率区域发生的位置，以便能够最大化地分配警察资源。然后，作为后续分析，我们需要确定某些类型犯罪的热点是否随时间变化。这两个问题都可以通过核密度估计来回答。

注意事项

这个数据集是从data.police.uk/data/下载的。

你可以从 Packt GitHub 下载，网址是github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-Python/tree/master/Lesson09/Activity21-Activity22。

或者，为了直接从源下载数据，前往前面的警察网站，勾选大都市警察局，然后设置日期范围为2018 年 7 月至2018 年 12 月。接下来，点击生成文件，然后点击立即下载并将下载的文件命名为metro-jul18-dec18。确保你知道如何获取或能够找到下载目录的路径。

这个数据集包含根据开放政府许可证 v3.0 许可的公共部门信息。

以下是完成活动的步骤：

加载犯罪数据。使用你保存下载目录的路径，创建年份-月份标签的列表，使用read_csv命令逐个加载文件，然后将这些文件合并在一起。
打印完整（六个月）和合并数据集的诊断信息。
将数据框架（DataFrame）缩小到四个变量（经度、纬度、月份和犯罪类型）。
使用seaborn中的jointplot函数，为 2018 年 7 月、9 月和 12 月的自行车盗窃数据拟合并可视化三个核密度估计模型。
重复第 4 步；这次，使用 2018 年 8 月、10 月和 11 月的商店盗窃犯罪数据。
重复第 5 步；这次，使用 2018 年 7 月、10 月和 12 月的入室盗窃犯罪数据。

第 6 步 的最后输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_28.jpg

图 9.28：2018 年 12 月入室盗窃的联合和边缘密度估计

再次强调，本活动中找到的密度应该在地图上叠加显示，以便我们能够看到这些密度覆盖的具体区域。如果您有合适的地图平台，请尝试自行叠加结果在地图上显示。如果没有，您可以访问在线地图服务，并使用经度和纬度对来获取关于具体位置的洞察。

注意

此活动的解决方案可在第 377 页找到。

总结

核密度估计是一种经典的统计技术，与直方图技术属于同一类。它允许用户从样本数据中推断出对特定对象或事件的人群进行洞察和预测。这种推断以概率密度函数的形式呈现，这非常好，因为结果可以解读为概率或可能性。这个模型的质量取决于两个参数：带宽值和核函数。正如讨论的那样，成功利用核密度估计的关键组成部分是设置一个最佳的带宽值。最常用的方法是使用网格搜索交叉验证，以伪对数似然作为评分指标来确定最佳带宽。核密度估计的优点在于其简单性和适用性广泛。

在犯罪学、流行病学、气象学和房地产等多个领域中，经常可以找到核密度估计模型。无论您从事哪个领域的业务，核密度估计都应该适用。

在本书中，我们探讨了如何在 Python 库的支持下，使用无监督学习技术的最佳实践，并从非结构化数据中提取有意义的信息。现在，您可以自信地使用 Python 构建自己的模型。

第十一章：附录

关于

本节内容旨在帮助学生完成书中列出的活动。它包括学生需要执行的详细步骤，以完成并实现书中的目标。

第一章：聚类介绍

活动 1：实现 k-means 聚类

解决方案：

使用 pandas 加载 Iris 数据文件，pandas 是一个通过使用 DataFrame 使数据处理变得更加容易的包：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import silhouette_score
from scipy.spatial.distance import cdist
iris = pd.read_csv('iris_data.csv', header=None)
iris.columns = ['SepalLengthCm', 'SepalWidthCm', 'PetalLengthCm', 'PetalWidthCm', 'species']

将 X 特征和提供的 y 物种标签分开，因为我们希望将其视为无监督学习问题：

X = iris[['SepalLengthCm', 'SepalWidthCm', 'PetalLengthCm', 'PetalWidthCm']]
y = iris['species']

了解我们的特征是什么样的：
```
X.head()
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_01_22.jpg

图 1.22：数据的前五行

将我们之前创建的 k_means 函数拿出来参考：

def k_means(X, K):
#Keep track of history so you can see k-means in action
    centroids_history = []
    labels_history = []
    rand_index = np.random.choice(X.shape[0], K)  
    centroids = X[rand_index]
    centroids_history.append(centroids)
    while True:
# Euclidean distances are calculated for each point relative to centroids, #and then np.argmin returns
# the index location of the minimal distance - which cluster a point    is #assigned to
        labels = np.argmin(cdist(X, centroids), axis=1)
        labels_history.append(labels)
#Take mean of points within clusters to find new centroids:
        new_centroids = np.array([X[labels == i].mean(axis=0)
                                for i in range(K)])
        centroids_history.append(new_centroids)

        # If old centroids and new centroids no longer change, k-means is complete and end. Otherwise continue
        if np.all(centroids == new_centroids):
            break
        centroids = new_centroids

    return centroids, labels, centroids_history, labels_history

将我们的 Iris X 特征 DataFrame 转换为 NumPy 矩阵：
```
X_mat = X.values
```

在鸢尾矩阵上运行我们的 k_means 函数：

centroids, labels, centroids_history, labels_history = k_means(X_mat, 3)

查看我们通过查看每个样本的预测物种列表得到的标签：
```
print(labels)
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_01_23.jpg

图 1.23：预测物种列表
可视化我们在数据集上实现的 k-means 方法：
```
plt.scatter(X['SepalLengthCm'], X['SepalWidthCm'])
plt.title('Iris - Sepal Length vs Width')
plt.show()
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_01_24.jpg

图 1.24：执行的 k-means 实现的图

如下所示可视化鸢尾物种的簇：
```
plt.scatter(X['SepalLengthCm'], X['SepalWidthCm'], c=labels, cmap='tab20b')
plt.title('Iris - Sepal Length vs Width - Clustered')
plt.show()
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_01_25.jpg

图 1.25：鸢尾物种的簇
使用 scikit-learn 实现计算轮廓系数（Silhouette Score）：
```
# Calculate Silhouette Score
silhouette_score(X[['SepalLengthCm','SepalWidthCm']], labels)
```
你将得到一个大约等于 0.369 的 SSI。由于我们只使用了两个特征，这是可以接受的，结合最终图中展示的聚类成员可视化。

第二章：层次聚类

活动 2：应用连接标准

解决方案：

可视化我们在 练习 7 中创建的 x 数据集，构建层次结构：

from scipy.cluster.hierarchy import linkage, dendrogram, fcluster
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt
%matplotlib inline
# Generate a random cluster dataset to experiment on. X = coordinate points, y = cluster labels (not needed)
X, y = make_blobs(n_samples=1000, centers=8, n_features=2, random_state=800)
# Visualize the data
plt.scatter(X[:,0], X[:,1])
plt.show()

输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_02_20.jpg

图 2.20：生成的簇数据集的散点图

创建一个包含所有可能连接方法超参数的列表：

methods = ['centroid', 'single', 'complete', 'average', 'weighted']

遍历你刚才创建的列表中的每种方法，展示它们在同一数据集上的效果：

for method in methods:
    distances = linkage(X, method=method, metric="euclidean")
    clusters = fcluster(distances, 3, criterion="distance") 
    plt.title('linkage: ' + method)
    plt.scatter(X[:,0], X[:,1], c=clusters, cmap='tab20b')
    plt.show()

输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_02_21.jpg

图 2.21：所有方法的散点图

分析：

从前面的图中可以看出，通过简单地更改连接标准，可以显著改变聚类的效果。在这个数据集中，质心法和平均法最适合找到合理的离散簇。这一点从我们生成的八个簇的事实中可以看出，质心法和平均法是唯一显示出使用八种不同颜色表示的簇的算法。其他连接类型则表现不佳，尤其是单链接法。

活动 3：比较 k-means 与层次聚类

解决方案：

从 scikit-learn 导入必要的包（KMeans、AgglomerativeClustering 和 silhouette_score），如下所示：

from sklearn.cluster import KMeans
from sklearn.cluster import AgglomerativeClustering
from sklearn.metrics import silhouette_score
import pandas as pd
import matplotlib.pyplot as plt

将葡萄酒数据集读入 pandas DataFrame 并打印一个小样本：
```
wine_df = pd.read_csv("wine_data.csv")
print(wine_df.head)
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_02_22.jpg

图 2.22：葡萄酒数据集的输出
可视化葡萄酒数据集以理解数据结构：
```
plt.scatter(wine_df.values[:,0], wine_df.values[:,1])
plt.title("Wine Dataset")
plt.xlabel("OD Reading")
plt.ylabel("Proline")
plt.show()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_02_23.jpg

图 2.23：原始葡萄酒数据的聚类图
在葡萄酒数据集上使用 sklearn 实现的 k-means，已知有三种葡萄酒类型：
```
km = KMeans(3)
km_clusters = km.fit_predict(wine_df)
```

使用 sklearn 实现的层次聚类对葡萄酒数据集进行处理：

ac = AgglomerativeClustering(3, linkage='average')
ac_clusters = ac.fit_predict(wine_df)

如下所示，绘制 k-means 预测的聚类：
```
plt.scatter(wine_df.values[:,0], wine_df.values[:,1], c=km_clusters)
plt.title("Wine Clusters from Agglomerative Clustering")
plt.xlabel("OD Reading")
plt.ylabel("Proline")
plt.show()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_02_24.jpg

图 2.24：k-means 聚类的聚类图
如下所示，绘制层次聚类预测的聚类：
```
plt.scatter(wine_df.values[:,0], wine_df.values[:,1], c=ac_clusters)
plt.title("Wine Clusters from Agglomerative Clustering")
plt.xlabel("OD Reading")
plt.ylabel("Proline")
plt.show()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_02_25.jpg

图 2.25：凝聚聚类的聚类图

比较每种聚类方法的轮廓得分：

print("Silhouette Scores for Wine Dataset:\n")
print("k-means Clustering: ", silhouette_score(X[:,11:13], km_clusters))
print("Agg Clustering: ", silhouette_score(X[:,11:13], ac_clusters))

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_02_26.jpg

图 2.26：葡萄酒数据集的轮廓得分

从之前的轮廓得分可以看出，凝聚聚类在分离聚类时，平均簇内距离上略微优于 k-means 聚类。然而，并非所有版本的凝聚聚类都有这个效果。你可以尝试不同的连接类型，观察不同轮廓得分和聚类结果如何变化！

第三章：邻域方法与 DBSCAN

活动 4：从零实现 DBSCAN

解决方案：

生成一个随机的聚类数据集，如下所示：

from sklearn.cluster import DBSCAN
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
X_blob, y_blob = make_blobs(n_samples=500, centers=4, n_features=2, random_state=800)

可视化生成的数据：
```
plt.scatter(X_blob[:,0], X_blob[:,1])
plt.show()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_03_14.jpg

图 3.14：生成数据的聚类图

创建从零开始的函数，允许你在数据集上调用 DBSCAN：

def scratch_DBSCAN(x, eps, min_pts):
    """
    param x (list of vectors): your dataset to be clustered
    param eps (float): neigborhood radius threshold
    param min_pts (int): minimum number of points threshold for a nieghborhood to be a cluster
    """
     # Build a label holder that is comprised of all 0s
    labels = [0]* x.shape[0]
    # Arbitrary starting "current cluster" ID    
    C = 0

    # For each point p in x...
    # ('p' is the index of the datapoint, rather than the datapoint itself.)
    for p in range(0, x.shape[0]):

        # Only unvisited points can be evaluated as neighborhood centers
        if not (labels[p] == 0):
            continue

        # Find all of p's neighbors.
        neighbors = neighborhood_search(x, p, eps)

        # If there are not enough neighbor points, then it is classified as noise (-1).
        # Otherwise we can use this point as a neighborhood cluster
        if len(neighbors) < min_pts:
            labels[p] = -1    
        else: 
            C += 1
            neighbor_cluster(x, labels, p, neighbors, C, eps, min_pts)

    return labels
def neighbor_cluster(x, labels, p, neighbors, C, eps, min_pts):
    # Assign the cluster label to original point
    labels[p] = C

    # Look at each neighbor of p (by index, not the points themselves) and evaluate
    i = 0
    while i < len(neighbors):    

        # Get the next point from the queue.        
        potential_neighbor_ix = neighbors[i]

        # If potential_neighbor_ix is noise from previous runs, we can assign it to current cluster
        if labels[potential_neighbor_ix] == -1:
            labels[potential_neighbor_ix] = C

        # Otherwise, if potential_neighbor_ix is unvisited, we can add it to current cluster
        elif labels[potential_neighbor_ix] == 0:
            labels[potential_neighbor_ix] = C

            # Further find neighbors of potential neighbor
            potential_neighbors_cluster = neighborhood_search(x, potential_neighbor_ix, eps)

            if len(potential_neighbors_cluster) >= min_pts:
                neighbors = neighbors + potential_neighbors_cluster      

        # Evaluate next neighbor
        i += 1        
def neighborhood_search(x, p, eps):
    neighbors = []

    # For each point in the dataset...
    for potential_neighbor in range(0, x.shape[0]):

        # If a nearby point falls below the neighborhood radius threshold, add to neighbors list
        if np.linalg.norm(x[p] - x[potential_neighbor]) < eps:
            neighbors.append(potential_neighbor)

    return neighbors

使用你创建的 DBSCAN 实现，在生成的数据集中查找聚类。根据第五步中的性能，随意调整超参数：
```
labels = scratch_DBSCAN(X_blob, 0.6, 5)
```

从零开始可视化你实现的 DBSCAN 聚类性能：

plt.scatter(X_blob[:,0], X_blob[:,1], c=labels)
plt.title("DBSCAN from Scratch Performance")
plt.show()

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_03_15.jpg

图 3.15：DBSCAN 实现的聚类图

正如你可能注意到的，定制实现的运行时间相对较长。这是因为为了清晰起见，我们探索了该算法的非矢量化版本。接下来，你应该尽量使用 scikit-learn 提供的 DBSCAN 实现，因为它经过高度优化。

活动 5：比较 DBSCAN 与 k-means 和层次聚类

解决方案：

导入必要的包：

from sklearn.cluster import KMeans, AgglomerativeClustering, DBSCAN
from sklearn.metrics import silhouette_score
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

从第二章 Hierarchical Clustering加载葡萄酒数据集，并再次熟悉数据的外观：
```
# Load Wine data set
wine_df = pd.read_csv("../CH2/wine_data.csv")
# Show sample of data set
print(wine_df.head())
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_03_16.jpg

图 3.16: 葡萄酒数据集的前五行
可视化数据：
```
plt.scatter(wine_df.values[:,0], wine_df.values[:,1])
plt.title("Wine Dataset")
plt.xlabel("OD Reading")
plt.ylabel("Proline")
plt.show()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_03_17.jpg

图 3.17: 数据的图

使用 k-means、凝聚聚类和 DBSCAN 生成聚类：

# Generate clusters from K-Means
km = KMeans(3)
km_clusters = km.fit_predict(wine_df)
# Generate clusters using Agglomerative Hierarchical Clustering
ac = AgglomerativeClustering(3, linkage='average')
ac_clusters = ac.fit_predict(wine_df)

评估 DSBSCAN 超参数的几个不同选项及其对轮廓分数的影响：

db_param_options = [[20,5],[25,5],[30,5],[25,7],[35,7],[35,3]]
for ep,min_sample in db_param_options:
    # Generate clusters using DBSCAN
    db = DBSCAN(eps=ep, min_samples = min_sample)
    db_clusters = db.fit_predict(wine_df)
    print("Eps: ", ep, "Min Samples: ", min_sample)
    print("DBSCAN Clustering: ", silhouette_score(wine_df, db_clusters))

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_03_18.jpg

图 3.18: 打印聚类的轮廓分数

根据最高轮廓分数生成最终聚类 (eps: 35, min_samples: 3)：

# Generate clusters using DBSCAN
db = DBSCAN(eps=35, min_samples = 3)
db_clusters = db.fit_predict(wine_df)

可视化使用三种方法生成的聚类：

plt.title("Wine Clusters from K-Means")
plt.scatter(wine_df['OD_read'], wine_df['Proline'], c=km_clusters,s=50, cmap='tab20b')
plt.show()
plt.title("Wine Clusters from Agglomerative Clustering")
plt.scatter(wine_df['OD_read'], wine_df['Proline'], c=ac_clusters,s=50, cmap='tab20b')
plt.show()
plt.title("Wine Clusters from DBSCAN")
plt.scatter(wine_df['OD_read'], wine_df['Proline'], c=db_clusters,s=50, cmap='tab20b')
plt.show()

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_03_19.jpg

图 3.19: 使用不同算法绘制聚类图

评估每种方法的轮廓分数：

# Calculate Silhouette Scores
print("Silhouette Scores for Wine Dataset:\n")
print("K-Means Clustering: ", silhouette_score(wine_df, km_clusters))
print("Agg Clustering: ", silhouette_score(wine_df, ac_clusters))
print("DBSCAN Clustering: ", silhouette_score(wine_df, db_clusters))

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_03_20.jpg

图 3.20: Silhouette 分数

如您所见，DBSCAN 并不总是自动适合您的聚类需求。使其与其他算法不同的一个关键特征是将噪声用作潜在聚类。在某些情况下，这很好，因为它可以去除离群值，但是，可能存在调整不足的情况，会将太多点分类为噪声。您能通过调整超参数来提高轮廓分数吗？

第四章: 维度缩减和 PCA

活动 6: 手动 PCA 与 scikit-learn 对比

解决方案

导入 pandas、numpy 和 matplotlib 绘图库以及 scikit-learn 的 PCA 模型：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

加载数据集，并根据以前的练习选择仅选择萼片特征。显示数据的前五行：
```
df = pd.read_csv('iris-data.csv')
df = df[['Sepal Length', 'Sepal Width']]
df.head()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_04_43.jpg

图 4.43: 数据的前五行
计算数据的协方差矩阵：
```
cov = np.cov(df.values.T)
cov
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_04_44.jpg

图 4.44: 数据的协方差矩阵
使用 scikit-learn API 转换数据，仅使用第一个主成分。将转换后的数据存储在 sklearn_pca 变量中：
```
model = PCA(n_components=1)
sklearn_pca = model.fit_transform(df.values)
```

使用手动 PCA 仅使用第一个主成分来转换数据。将转换后的数据存储在 manual_pca 变量中。

eigenvectors, eigenvalues, _ = np.linalg.svd(cov, full_matrices=False)
P = eigenvectors[0]
manual_pca = P.dot(df.values.T)

在同一图上绘制 sklearn_pca 和 manual_pca 的值以可视化差异：

plt.figure(figsize=(10, 7));
plt.plot(sklearn_pca, label='Scikit-learn PCA');
plt.plot(manual_pca, label='Manual PCA', linestyle='--');
plt.xlabel('Sample');
plt.ylabel('Transformed Value');
plt.legend();

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_04_45.jpg

图 4.45: 数据的图

请注意，两个图几乎完全相同，唯一的区别是一个是另一个的镜像，且两者之间存在偏移。显示sklearn_pca和manual_pca模型的主成分：
```
model.components_
```
输出如下：
```
array([[ 0.99693955, -0.07817635]])
```
现在打印P：
```
P
```
输出如下：
```
array([-0.99693955,  0.07817635])
```
注意符号的差异；值是相同的，但符号不同，产生了镜像的结果。这只是约定上的差异，并无实质意义。

将manual_pca模型乘以-1并重新绘制：

manual_pca *= -1
plt.figure(figsize=(10, 7));
plt.plot(sklearn_pca, label='Scikit-learn PCA');
plt.plot(manual_pca, label='Manual PCA', linestyle='--');
plt.xlabel('Sample');
plt.ylabel('Transformed Value');
plt.legend();

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_04_46.jpg

图 4.46: 重新绘制的数据

现在，我们需要做的就是处理两个之间的偏移。scikit-learn API 在变换之前会减去数据的均值。在使用手动 PCA 进行变换之前，从数据集中减去每列的均值：
```
mean_vals = np.mean(df.values, axis=0)
offset_vals = df.values - mean_vals
manual_pca = P.dot(offset_vals.T)
```
将结果乘以-1：
```
manual_pca *= -1
```

重新绘制单独的sklearn_pca和manual_pca值：

plt.figure(figsize=(10, 7));
plt.plot(sklearn_pca, label='Scikit-learn PCA');
plt.plot(manual_pca, label='Manual PCA', linestyle='--');
plt.xlabel('Sample');
plt.ylabel('Transformed Value');
plt.legend();

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_04_47.jpg

图 4.47: 重新绘制数据

最终的图将展示通过这两种方法完成的降维实际上是相同的。不同之处在于协方差矩阵的符号差异，因为这两种方法只是使用了不同的特征作为比较的基准。最后，两数据集之间也存在偏移，这是由于在执行 scikit-learn PCA 变换之前已将均值样本减去。

活动 7: 使用扩展的鸢尾花数据集进行主成分分析（PCA）

解决方案：

导入pandas和matplotlib。为了启用 3D 绘图，您还需要导入Axes3D：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from mpl_toolkits.mplot3d import Axes3D # Required for 3D plotting

读取数据集并选择花萼长度、花萼宽度和花瓣宽度列：
```
df = pd.read_csv('iris-data.csv')[['Sepal Length', 'Sepal Width', 'Petal Width']]
df.head()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_04_48.jpg

图 4.48: 花萼长度、花萼宽度和花瓣宽度

绘制三维数据：

fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(df['Sepal Length'], df['Sepal Width'], df['Petal Width']);
ax.set_xlabel('Sepal Length (mm)');
ax.set_ylabel('Sepal Width (mm)');
ax.set_zlabel('Petal Width (mm)');
ax.set_title('Expanded Iris Dataset');

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_04_49.jpg

图 4.49: 扩展的鸢尾花数据集图

创建一个PCA模型，不指定主成分数：
```
model = PCA()
```
将模型拟合到数据集：
```
model.fit(df.values)
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_04_50.jpg

图 4.50: 拟合到数据集的模型

显示特征值或explained_variance_ratio_：

model.explained_variance_ratio_

输出如下：

array([0.8004668 , 0.14652357, 0.05300962])

我们希望减少数据集的维度，但仍保持至少 90% 的方差。为了保持 90% 的方差，所需的最小主成分数是多少？

前两个主成分需要至少 90% 的方差。前两个主成分提供了数据集中 94.7% 的方差。
创建一个新的PCA模型，这次指定所需的主成分数，以保持至少 90% 的方差：
```
model = PCA(n_components=2)
```

使用新模型转换数据：

data_transformed = model.fit_transform(df.values)

绘制转换后的数据：
```
plt.figure(figsize=(10, 7))
plt.scatter(data_transformed[:,0], data_transformed[:,1]);
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_04_51.jpg

图 4.51：转换数据的图

将转换后的数据恢复到原始数据空间：

data_restored = model.inverse_transform(data_transformed)

在一个子图中绘制恢复后的三维数据，在第二个子图中绘制原始数据，以可视化去除部分方差的效果：

fig = plt.figure(figsize=(10, 14))
# Original Data
ax = fig.add_subplot(211, projection='3d')
ax.scatter(df['Sepal Length'], df['Sepal Width'], df['Petal Width'], label='Original Data');
ax.set_xlabel('Sepal Length (mm)');
ax.set_ylabel('Sepal Width (mm)');
ax.set_zlabel('Petal Width (mm)');
ax.set_title('Expanded Iris Dataset');
# Transformed Data
ax = fig.add_subplot(212, projection='3d')
ax.scatter(data_restored[:,0], data_restored[:,1], data_restored[:,2], label='Restored Data');
ax.set_xlabel('Sepal Length (mm)');
ax.set_ylabel('Sepal Width (mm)');
ax.set_zlabel('Petal Width (mm)');
ax.set_title('Restored Iris Dataset');

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_04_52.jpg

图 4.52：扩展和恢复后的鸢尾花数据集图

看看图 4.52，我们可以看到，就像我们在二维图中做的那样，我们已经去除了数据中的大量噪声，但保留了关于数据趋势的最重要信息。可以看出，通常情况下，花萼长度与花瓣宽度成正比，并且在图中似乎有两个数据簇，一个在另一个上方。

注意

在应用 PCA 时，重要的是要考虑所建模数据的大小以及可用的系统内存。奇异值分解过程涉及将数据分解为特征值和特征向量，且可能非常占用内存。如果数据集过大，你可能无法完成该过程，或者会遭遇显著的性能损失，甚至可能导致系统崩溃。

第五章：自编码器

活动 8：使用 ReLU 激活函数建模神经元

解决方案：

导入numpy和 matplotlib：

import numpy as np
import matplotlib.pyplot as plt

允许在标签中使用 latex 符号：
```
plt.rc('text', usetex=True)
```
将 ReLU 激活函数定义为 Python 函数：
```
def relu(x):
    return np.max((0, x))
```
定义神经元的输入（x）和可调权重（theta）。在此示例中，输入（x）将是-5 到 5 之间线性间隔的 100 个数字。设置theta = 1：
```
theta = 1
x = np.linspace(-5, 5, 100)
x
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_05_35.jpg

图 5.35：打印输入
计算输出（y）：
```
y = [relu(_x * theta) for _x in x]
```

绘制神经元的输出与输入的关系图：

fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111)
ax.plot(x, y)
ax.set_xlabel('$x$', fontsize=22);
ax.set_ylabel('$h(x\Theta)$', fontsize=22);
ax.spines['left'].set_position(('data', 0));
ax.spines['top'].set_visible(False);
ax.spines['right'].set_visible(False);
ax.tick_params(axis='both', which='major', labelsize=22)

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_05_36.jpg

图 5.36：神经元与输入的关系图

现在，将theta = 5并重新计算并存储神经元的输出：
```
theta = 5
y_2 = [relu(_x * theta) for _x in x]
```
现在，将theta = 0.2并重新计算并存储神经元的输出：
```
theta = 0.2
y_3 = [relu(_x * theta) for _x in x]
```

在一张图上绘制神经元的三条不同输出曲线（theta = 1，theta = 5，theta = 0.2）：

fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111)
ax.plot(x, y, label='$\Theta=1$');
ax.plot(x, y_2, label='$\Theta=5$', linestyle=':');
ax.plot(x, y_3, label='$\Theta=0.2$', linestyle='--');
ax.set_xlabel('$x\Theta$', fontsize=22);
ax.set_ylabel('$h(x\Theta)$', fontsize=22);
ax.spines['left'].set_position(('data', 0));
ax.spines['top'].set_visible(False);
ax.spines['right'].set_visible(False);
ax.tick_params(axis='both', which='major', labelsize=22);
ax.legend(fontsize=22);

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_05_37.jpg

图 5.37：神经元的三条输出曲线

在这个活动中，我们创建了一个基于 ReLU 的人工神经网络神经元模型。我们可以看到，这个神经元的输出与 sigmoid 激活函数的输出非常不同。对于大于 0 的值，没有饱和区域，因为它仅返回函数的输入值。在负方向上，有一个饱和区域，如果输入小于 0，则返回 0。ReLU 函数是一种非常强大且常用的激活函数，在某些情况下，它比 sigmoid 函数更强大。ReLU 通常是一个很好的首选激活函数。

活动 9：MNIST 神经网络

解决方案：

在这个活动中，您将训练一个神经网络来识别 MNIST 数据集中的图像，并强化您的神经网络训练技能：

导入pickle、numpy、matplotlib，以及从 Keras 中导入Sequential和Dense类：

import pickle
import numpy as np
import matplotlib.pyplot as plt
from keras.models import Sequential
from keras.layers import Dense

加载mnist.pkl文件，其中包含来自 MNIST 数据集的前 10,000 张图像及其对应的标签，这些数据可以在随附的源代码中找到。MNIST 数据集是一系列 28 x 28 灰度图像，表示手写数字 0 到 9。提取图像和标签：
```
with open('mnist.pkl', 'rb') as f:
    data = pickle.load(f)

images = data['images']
labels = data['labels']
```

绘制前 10 个样本及其对应的标签：

plt.figure(figsize=(10, 7))
for i in range(10):
    plt.subplot(2, 5, i + 1)
    plt.imshow(images[i], cmap='gray')
    plt.title(labels[i])
    plt.axis('off')

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_05_38.jpg

图 5.38：前 10 个样本

使用独热编码对标签进行编码：
```
one_hot_labels = np.zeros((images.shape[0], 10))
for idx, label in enumerate(labels):
    one_hot_labels[idx, label] = 1

one_hot_labels
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_05_39.jpg

图 5.39：独热编码结果
为神经网络输入准备图像。提示：此过程包含两个独立的步骤：
```
images = images.reshape((-1, 28 ** 2))
images = images / 255.
```
在 Keras 中构建一个神经网络模型，该模型接受准备好的图像，包含一个具有 600 个单位的隐藏层，并使用 ReLU 激活函数，输出层的单元数与类别数相同。输出层使用softmax激活函数：
```
model = Sequential([
    Dense(600, input_shape=(784,), activation='relu'),
    Dense(10, activation='softmax'),
])
```

使用多类交叉熵、随机梯度下降以及准确度作为性能指标来编译模型：

model.compile(loss='categorical_crossentropy',
              optimizer='sgd',
              metrics=['accuracy'])

训练模型。需要多少个 epoch 才能在训练数据上达到至少 95%的分类准确率？让我们来看看：
```
model.fit(images, one_hot_labels, epochs=20)
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_05_40.jpg

图 5.40：训练模型

需要 15 个 epoch 才能在训练集上达到至少 95%的分类准确率。

在这个示例中，我们使用分类器训练时的数据来测量神经网络分类器的性能。一般来说，这种方法不应该使用，因为它通常报告比实际模型应该有的准确度更高。在监督学习问题中，有多种交叉验证技术应该被使用。由于这是一本关于无监督学习的书，交叉验证超出了本书的范围。

活动 10：简单的 MNIST 自编码器

解决方案：

导入 pickle、numpy 和 matplotlib，以及从 Keras 导入 Model、Input 和 Dense 类：

import pickle
import numpy as np
import matplotlib.pyplot as plt
from keras.models import Model
from keras.layers import Input, Dense

从提供的 MNIST 数据集样本（随附源代码 mnist.pkl）加载图像：

with open('mnist.pkl', 'rb') as f:
    images = pickle.load(f)['images']

准备图像以输入神经网络。作为提示，此过程有两个独立步骤：
```
images = images.reshape((-1, 28 ** 2))
images = images / 255.
```

构建一个简单的自编码器网络，使图像在编码阶段后缩小到 10 x 10：

input_stage = Input(shape=(784,))
encoding_stage = Dense(100, activation='relu')(input_stage)
decoding_stage = Dense(784, activation='sigmoid')(encoding_stage)
autoencoder = Model(input_stage, decoding_stage)

使用二元交叉熵损失函数和 adadelta 梯度下降法编译自编码器：

autoencoder.compile(loss='binary_crossentropy',
              optimizer='adadelta')

拟合编码器模型：
```
autoencoder.fit(images, images, epochs=100)
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_05_41.jpg

图 5.41: 训练模型

计算并存储编码阶段的前五个样本的输出：

encoder_output = Model(input_stage, encoding_stage).predict(images[:5])

将编码器输出重塑为 10 x 10（10 x 10 = 100）像素并乘以 255：
```
encoder_output = encoder_output.reshape((-1, 10, 10)) * 255
```
计算并存储解码阶段的前五个样本的输出：
```
decoder_output = autoencoder.predict(images[:5])
```

将解码器的输出重塑为 28 x 28，并乘以 255：

decoder_output = decoder_output.reshape((-1, 28, 28)) * 255

绘制原始图像、编码器输出和解码器：

images = images.reshape((-1, 28, 28))
plt.figure(figsize=(10, 7))
for i in range(5):
    plt.subplot(3, 5, i + 1)
    plt.imshow(images[i], cmap='gray')
    plt.axis('off')
    plt.subplot(3, 5, i + 6)
    plt.imshow(encoder_output[i], cmap='gray')
    plt.axis('off')   

    plt.subplot(3, 5, i + 11)
    plt.imshow(decoder_output[i], cmap='gray')
    plt.axis('off')

输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_05_42.jpg

图 5.42: 原始图像、编码器输出和解码器

到目前为止，我们已经展示了如何使用编码和解码阶段的简单单层隐藏层来将数据降到低维空间。我们还可以通过在编码和解码阶段添加额外的层来使这个模型更加复杂。

活动 11: MNIST 卷积自编码器

解决方案：

导入 pickle、numpy、matplotlib，以及从 keras.models 导入 Model 类，导入从 keras.layers 的 Input、Conv2D、MaxPooling2D 和 UpSampling2D：

import pickle
import numpy as np
import matplotlib.pyplot as plt
from keras.models import Model
from keras.layers import Input, Conv2D, MaxPooling2D, UpSampling2D

加载数据：

with open('mnist.pkl', 'rb') as f:
    images = pickle.load(f)['images']

将图像重新缩放为 0 和 1 之间的值：
```
images = images / 255.
```
我们需要重塑图像，以便为卷积阶段添加一个单独的深度通道。将图像重塑为 28 x 28 x 1 的形状：
```
images = images.reshape((-1, 28, 28, 1))
```
定义输入层。我们将使用与图像相同的输入形状：
```
input_layer = Input(shape=(28, 28, 1,))
```

添加一个具有 16 层或滤波器的卷积阶段，使用 3 x 3 的权重矩阵、ReLU 激活函数，并使用相同的填充，这意味着输出的长度与输入图像相同：

hidden_encoding = Conv2D(
    16, # Number of layers or filters in the weight matrix
    (3, 3), # Shape of the weight matrix
    activation='relu',
    padding='same', # How to apply the weights to the images
)(input_layer)

向编码器添加一个最大池化层，使用 2 x 2 的内核：
```
encoded = MaxPooling2D((2, 2))(hidden_encoding)
```

添加一个解码卷积层：

hidden_decoding = Conv2D(
    16, # Number of layers or filters in the weight matrix
    (3, 3), # Shape of the weight matrix
    activation='relu',
    padding='same', # How to apply the weights to the images
)(encoded)

添加上采样层：

upsample_decoding = UpSampling2D((2, 2))(hidden_decoding)

添加最后的卷积阶段，使用一个层，按照初始图像深度进行：

decoded = Conv2D(
    1, # Number of layers or filters in the weight matrix
    (3, 3), # Shape of the weight matrix
    activation='sigmoid',
    padding='same', # How to apply the weights to the images
)(upsample_decoding)

通过将网络的第一层和最后一层传递给 Model 类来构建模型：
```
autoencoder = Model(input_layer, decoded)
```
显示模型的结构：
```
autoencoder.summary()
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_05_43.jpg

图 5.43: 模型结构

使用二元交叉熵损失函数和 adadelta 梯度下降法编译自编码器：

autoencoder.compile(loss='binary_crossentropy',
              optimizer='adadelta')

现在，我们来拟合模型；再次传递图像作为训练数据和期望输出。由于卷积网络计算时间较长，因此训练 20 个周期：
```
autoencoder.fit(images, images, epochs=20)
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_05_44.jpg

图 5.44：训练模型

计算并存储前五个样本的编码阶段输出：

encoder_output = Model(input_layer, encoded).predict(images[:5])

为可视化重新调整编码器输出的形状，每个图像的大小为 X*Y：
```
encoder_output = encoder_output.reshape((-1, 14 * 14, 16))
```

获取前五张图像的解码器输出：

decoder_output = autoencoder.predict(images[:5])

将解码器输出调整为 28 x 28 的大小：

decoder_output = decoder_output.reshape((-1, 28, 28))

将原始图像调整回 28 x 28 的大小：
```
images = images.reshape((-1, 28, 28))
```

绘制原始图像、平均编码器输出和解码器：

plt.figure(figsize=(10, 7))
for i in range(5):
    plt.subplot(3, 5, i + 1)
    plt.imshow(images[i], cmap='gray')
    plt.axis('off')

    plt.subplot(3, 5, i + 6)
    plt.imshow(encoder_output[i], cmap='gray')
    plt.axis('off')   

    plt.subplot(3, 5, i + 11)
    plt.imshow(decoder_output[i], cmap='gray')
    plt.axis('off')

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_05_45.jpg

图 5.45：原始图像、编码器输出和解码器

在此活动结束时，你将开发一个包含卷积层的自编码器神经网络。注意解码器表示中的改进。与全连接神经网络层相比，这种架构在性能上具有显著优势，非常适用于处理基于图像的数据集和生成人工数据样本。

第六章：t-分布随机邻域嵌入（t-SNE）

活动 12：葡萄酒 t-SNE

解决方案：

导入 pandas、numpy、matplotlib 以及来自 scikit-learn 的 t-SNE 和 PCA 模型：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE

使用随附源代码中的 wine.data 文件加载葡萄酒数据集，并显示数据集的前五行：
```
df = pd.read_csv('wine.data', header=None)
df.head()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_24.jpg

图 6.24：葡萄酒数据集的前五行。
第一列包含标签；提取该列并将其从数据集中移除：
```
labels = df[0]
del df[0]
```

执行 PCA 降维，将数据集缩减到前六个组件：

model_pca = PCA(n_components=6)
wine_pca = model_pca.fit_transform(df)

确定这六个组件描述的数据中的方差量：
```
np.sum(model_pca.explained_variance_ratio_)
```
输出如下：
```
0.99999314824536
```
创建一个使用指定随机状态和 verbose 值为 1 的 t-SNE 模型：
```
tsne_model = TSNE(random_state=0, verbose=1)
tsne_model
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_25.jpg

图 6.25：创建 t-SNE 模型。
将 PCA 数据拟合到 t-SNE 模型：
```
wine_tsne = tsne_model.fit_transform(wine_pca.reshape((len(wine_pca), -1)))
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_26.jpg

图 6.26：拟合 PCA 数据 t-SNE 模型
确认 t-SNE 拟合数据的形状是二维的：
```
wine_tsne.shape
```
输出如下：
```
(172, 8)
```
创建二维数据的散点图：
```
plt.figure(figsize=(10, 7))
plt.scatter(wine_tsne[:,0], wine_tsne[:,1]);
plt.title('Low Dimensional Representation of Wine');
plt.show()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_27.jpg

图 6.27：二维数据的散点图

创建一个带有类别标签的二维数据次级散点图，以可视化可能存在的任何聚类：

MARKER = ['o', 'v', '^',]
plt.figure(figsize=(10, 7))
plt.title('Low Dimensional Representation of Wine');
for i in range(1, 4):
    selections = wine_tsne[labels == i]
    plt.scatter(selections[:,0], selections[:,1], marker=MARKER[i-1], label=f'Wine {i}', s=30);
    plt.legend();
plt.show()

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_28.jpg

图 6.28：二维数据的次级图

请注意，虽然酒类之间有重叠，但也可以看到数据中存在一些聚类。第一类酒主要位于图表的左上角，第二类酒位于右下角，第三类酒则位于前两者之间。这个表示方法当然不能用于高信心地对单个酒样进行分类，但它展示了一个总体趋势以及在之前无法看到的高维数据中包含的一系列聚类。

活动 13：t-SNE 酒类与困惑度

解决方案：

导入 pandas、numpy、matplotlib 和 t-SNE 与 PCA 模型（来自 scikit-learn）：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE

加载酒类数据集并检查前五行：
```
df = pd.read_csv('wine.data', header=None)
df.head()
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_29.jpg

图 6.29：酒类数据的前五行
第一列提供了标签；从数据框中提取它们并将其存储在一个单独的变量中。确保该列从数据框中移除：
```
labels = df[0]
del df[0]
```

对数据集执行 PCA 并提取前六个主成分：

model_pca = PCA(n_components=6)
wine_pca = model_pca.fit_transform(df)
wine_pca = wine_pca.reshape((len(wine_pca), -1))

构建一个循环，遍历困惑度值（1、5、20、30、80、160、320）。对于每个循环，生成一个对应困惑度的 t-SNE 模型，并绘制标记酒类的散点图。注意不同困惑度值的效果：

MARKER = ['o', 'v', '^',]
for perp in [1, 5, 20, 30, 80, 160, 320]:
    tsne_model = TSNE(random_state=0, verbose=1, perplexity=perp)
    wine_tsne = tsne_model.fit_transform(wine_pca)
    plt.figure(figsize=(10, 7))
    plt.title(f'Low Dimensional Representation of Wine. Perplexity {perp}');
    for i in range(1, 4):
        selections = wine_tsne[labels == i]
        plt.scatter(selections[:,0], selections[:,1], marker=MARKER[i-1], label=f'Wine {i}', s=30);
        plt.legend();

1的困惑度值无法将数据分离成任何特定的结构：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_30.jpg

图 6.30：困惑度值为 1 的图表

将困惑度增加到 5 会导致非常非线性的结构，难以分离，并且很难识别任何聚类或模式：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_31.jpg

图 6.31：困惑度为 5 的图表

困惑度为 20 最终开始显示某种马蹄形结构。虽然在视觉上明显，但实现起来仍然可能有些困难：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_32.jpg

图 6.32：困惑度为 20 的图表

困惑度为 30 显示出相当好的结果。投影结构之间存在一定的线性关系，并且酒类之间有一些分离：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_33.jpg

图 6.33：困惑度为 30 的图表

最后，活动中的最后两张图片展示了随着困惑度的增加，图表如何变得越来越复杂和非线性：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_34.jpg

图 6.34：困惑度为 80 的图表

这是困惑度为 160 的图表：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_35.jpg

图 6.35：困惑度为 160 的图表

查看每个困惑度值的单独图形，困惑度对数据可视化的影响立刻显现出来。非常小或非常大的困惑度值会产生一系列不寻常的形状，无法显示出任何持续的模式。最合理的值似乎是 30，它产生了我们在前一个活动中看到的最线性的图形。

在这次活动中，我们展示了在选择困惑度时需要小心，并且可能需要一些迭代来确定正确的值。

活动 14：t-SNE 葡萄酒与迭代

解决方案：

导入pandas、numpy、matplotlib，以及从 scikit-learn 导入的t-SNE和PCA模型：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE

加载葡萄酒数据集并检查前五行：
```
df = pd.read_csv('wine.data', header=None)
df.head()
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_36.jpg

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_36.jpg)

图 6.36：葡萄酒数据集的前五行
第一列提供标签；从 DataFrame 中提取这些标签并存储到一个单独的变量中。确保将该列从 DataFrame 中移除：
```
labels = df[0]
del df[0]
```

在数据集上执行 PCA 并提取前六个主成分：

model_pca = PCA(n_components=6)
wine_pca = model_pca.fit_transform(df)
wine_pca = wine_pca.reshape((len(wine_pca), -1))

构建一个循环，遍历迭代值（250、500、1000）。对于每次循环，生成一个具有相应迭代次数的 t-SNE 模型，并生成相同迭代次数且没有进度值的模型：

MARKER = ['o', 'v', '1', 'p' ,'*', '+', 'x', 'd', '4', '.']
for iterations in [250, 500, 1000]:
    model_tsne = TSNE(random_state=0, verbose=1, n_iter=iterations, n_iter_without_progress=iterations)
    mnist_tsne = model_tsne.fit_transform(mnist_pca)

构建一个葡萄酒类别的散点图。注意不同迭代值的效果：

    plt.figure(figsize=(10, 7))
    plt.title(f'Low Dimensional Representation of MNIST (iterations = {iterations})');
    for i in range(10):
        selections = mnist_tsne[mnist['labels'] == i]
        plt.scatter(selections[:,0], selections[:,1], alpha=0.2, marker=MARKER[i], s=5);
        x, y = selections.mean(axis=0)
        plt.text(x, y, str(i), fontdict={'weight': 'bold', 'size': 30})

输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_37.jpg

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_37.jpg)

图 6.37：250 次迭代的葡萄酒类别散点图

这是 500 次迭代的图：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_38.jpg

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_38.jpg)

图 6.38：500 次迭代的葡萄酒类别散点图

这是 1,000 次迭代的图：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_39.jpg

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_06_39.jpg)

图 6.39：1,000 次迭代的葡萄酒类别散点图

再次，我们可以看到随着迭代次数的增加，数据结构的改善。即使在像这样的相对简单的数据集中，250 次迭代也不足以将数据结构投影到低维空间中。

正如我们在相应活动中观察到的，设置迭代参数时需要找到一个平衡点。在这个例子中，250 次迭代是不足够的，至少需要 1,000 次迭代才能最终稳定数据。

第七章：主题建模

活动 15：加载和清理 Twitter 数据

解决方案：

导入必要的库：

import langdetect
import matplotlib.pyplot
import nltk
import numpy
import pandas
import pyLDAvis
import pyLDAvis.sklearn
import regex
import sklearn

从github.com/TrainingByPackt/Applied-Unsupervised-Learning-with-Python/tree/master/Lesson07/Activity15-Activity17加载 LA Times 健康 Twitter 数据(latimeshealth.txt)：

注意
```
path = '<Path>/latimeshealth.txt'
df = pandas.read_csv(path, sep="|", header=None)
df.columns = ["id", "datetime", "tweettext"]
```

运行快速探索性分析以确定数据大小和结构：

def dataframe_quick_look(df, nrows):
print("SHAPE:\n{shape}\n".format(shape=df.shape))
print("COLUMN NAMES:\n{names}\n".format(names=df.columns))
print("HEAD:\n{head}\n".format(head=df.head(nrows)))
dataframe_quick_look(df, nrows=2)

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_07_54.jpg

图 7.54：数据的形状、列名和数据头

提取推文文本并转换为列表对象：
```
raw = df['tweettext'].tolist()
print("HEADLINES:\n{lines}\n".format(lines=raw[:5]))
print("LENGTH:\n{length}\n".format(length=len(raw)))
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_07_55.jpg

图 7.55：标题及其长度

编写函数，执行语言检测、基于空格的分词、替换屏幕名和 URL 为 SCREENNAME 和 URL。该函数还应删除标点、数字以及 SCREENNAME 和 URL 的替换。将所有内容转换为小写，除了 SCREENNAME 和 URL。它应删除所有停用词，执行词形还原，并保留五个或更多字母的单词：

注意

def do_language_identifying(txt):
    	try:
           the_language = langdetect.detect(txt)
    	except:
        	the_language = 'none'
    	return the_language
def do_lemmatizing(wrd):
    	out = nltk.corpus.wordnet.morphy(wrd)
    	return (wrd if out is None else out)
def do_tweet_cleaning(txt):
# identify language of tweet
# return null if language not english
    	lg = do_language_identifying(txt)
    	if lg != 'en':
        	return None
# split the string on whitespace
    	out = txt.split(' ')
# identify screen names
# replace with SCREENNAME
    	out = ['SCREENNAME' if i.startswith('@') else i for i in out]
# identify urls
# replace with URL
    	out = ['URL' if bool(regex.search('http[s]?://', i)) else i for i in out]
      # remove all punctuation
    	out = [regex.sub('[^\\w\\s]|\n', '', i) for i in out]
      # make all non-keywords lowercase
    	keys = ['SCREENNAME', 'URL']
    	out = [i.lower() if i not in keys else i for i in out]
      # remove keywords
    	out = [i for i in out if i not in keys]
      # remove stopwords
    	list_stop_words = nltk.corpus.stopwords.words('english')
    	list_stop_words = [regex.sub('[^\\w\\s]', '', i) for i in list_stop_words]
    	out = [i for i in out if i not in list_stop_words]
      # lemmatizing
    	out = [do_lemmatizing(i) for i in out]
      # keep words 4 or more characters long
    	out = [i for i in out if len(i) >= 5]
    	return out

对每条推文应用第 5 步中定义的函数：
```
clean = list(map(do_tweet_cleaning, raw))
```
删除等于 None 的输出列表元素：
```
clean = list(filter(None.__ne__, clean))
print("HEADLINES:\n{lines}\n".format(lines=clean[:5]))
print("LENGTH:\n{length}\n".format(length=len(clean)))
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_07_56.jpg

图 7.56：删除 None 后的标题和长度
将每条推文的元素转换回字符串。使用空格连接：
```
clean_sentences = [" ".join(i) for i in clean]
print(clean_sentences[0:10])
```
输出列表的前 10 个元素应如下所示：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_07_57.jpg

图 7.57：用于建模的清洁推文
保持笔记本打开以供将来建模。

活动 16：潜在狄利克雷分配与健康推文

解决方案：

指定 number_words、number_docs 和 number_features 变量：

number_words = 10
number_docs = 10
number_features = 1000

创建词袋模型，并将特征名分配给另一个变量以供以后使用：

vectorizer1 = sklearn.feature_extraction.text.CountVectorizer(
    analyzer=»word»,
    max_df=0.95, 
    min_df=10, 
    max_features=number_features
)
clean_vec1 = vectorizer1.fit_transform(clean_sentences)
print(clean_vec1[0])
feature_names_vec1 = vectorizer1.get_feature_names()

输出如下：

(0, 320)    1

确定最佳主题数：

def perplexity_by_ntopic(data, ntopics):
    output_dict = {
        «Number Of Topics": [], 
        «Perplexity Score»: []
    }
    for t in ntopics:
        lda = sklearn.decomposition.LatentDirichletAllocation(
            n_components=t,
            learning_method="online",
            random_state=0
        )
        lda.fit(data)
        output_dict["Number Of Topics"].append(t)
        output_dict["Perplexity Score"].append(lda.perplexity(data))
    output_df = pandas.DataFrame(output_dict)
    index_min_perplexity = output_df["Perplexity Score"].idxmin()
    output_num_topics = output_df.loc[
        index_min_perplexity,  # index
        «Number Of Topics"  # column
    ]
    return (output_df, output_num_topics)
df_perplexity, optimal_num_topics = perplexity_by_ntopic(
    clean_vec1, 
    ntopics=[i for i in range(1, 21) if i % 2 == 0]
)
print(df_perplexity)

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_07_58.jpg

图 7.58：主题数与困惑分数数据框架

使用最佳主题数拟合 LDA 模型：

lda = sklearn.decomposition.LatentDirichletAllocation(
    n_components=optimal_num_topics,
    learning_method="online",
    random_state=0
)
lda.fit(clean_vec1)

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_07_59.jpg

图 7.59：LDA 模型

创建并打印词-主题表：

def get_topics(mod, vec, names, docs, ndocs, nwords):
    # word to topic matrix
    W = mod.components_
    W_norm = W / W.sum(axis=1)[:, numpy.newaxis]
    # topic to document matrix
    H = mod.transform(vec)
    W_dict = {}
    H_dict = {}
    for tpc_idx, tpc_val in enumerate(W_norm):
        topic = «Topic{}".format(tpc_idx)
        # formatting w
        W_indices = tpc_val.argsort()[::-1][:nwords]
        W_names_values = [
            (round(tpc_val[j], 4), names[j]) 
            for j in W_indices
        ]
        W_dict[topic] = W_names_values
        # formatting h
        H_indices = H[:, tpc_idx].argsort()[::-1][:ndocs]
        H_names_values = [
        (round(H[:, tpc_idx][j], 4), docs[j]) 
            for j in H_indices
        ]
        H_dict[topic] = H_names_values
    W_df = pandas.DataFrame(
        W_dict, 
        index=["Word" + str(i) for i in range(nwords)]
    )
    H_df = pandas.DataFrame(
        H_dict,
        index=["Doc" + str(i) for i in range(ndocs)]
    )
    return (W_df, H_df)
W_df, H_df = get_topics(
    mod=lda,
    vec=clean_vec1,
    names=feature_names_vec1,
    docs=raw,
    ndocs=number_docs, 
    nwords=number_words
)
print(W_df)

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_07_60.jpg

图 7.60：健康推文数据的词-主题表

打印文档-主题表：
```
print(H_df)
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_07_61.jpg

图 7.61：文档主题表
创建双图可视化：
```
lda_plot = pyLDAvis.sklearn.prepare(lda, clean_vec1, vectorizer1, R=10)
pyLDAvis.display(lda_plot)
```
https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_07_62.jpg

图 7.62：在健康推文上训练的 LDA 模型的直方图和双图
保持笔记本打开以供将来建模。

活动 17：非负矩阵分解

解决方案：

创建适当的词袋模型，并将特征名作为另一个变量输出：

vectorizer2 = sklearn.feature_extraction.text.TfidfVectorizer(
    analyzer="word",
    max_df=0.5, 
    min_df=20, 
    max_features=number_features,
    smooth_idf=False
)
clean_vec2 = vectorizer2.fit_transform(clean_sentences)
print(clean_vec2[0])
feature_names_vec2 = vectorizer2.get_feature_names()

定义并使用活动二中的主题数（n_components）值拟合 NMF 算法：

nmf = sklearn.decomposition.NMF(
    n_components=optimal_num_topics,
    init="nndsvda",
    solver="mu",
    beta_loss="frobenius",
    random_state=0, 
    alpha=0.1, 
    l1_ratio=0.5
)
nmf.fit(clean_vec2)

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_07_63.jpg

图 7.63：定义 NMF 模型

获取主题-文档和单词-主题的结果表格。花几分钟时间探索单词分组，并尝试定义抽象的主题：
```
W_df, H_df = get_topics(
    mod=nmf,
    vec=clean_vec2,
    names=feature_names_vec2,
    docs=raw,
    ndocs=number_docs, 
    nwords=number_words
)
print(W_df)
```
https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_07_64.jpg

图 7.64：带有概率的单词-主题表
调整模型参数并重新运行 步骤 3 和 步骤 4。

第八章：市场购物篮分析

活动 18：加载和准备完整在线零售数据

解决方案：

加载在线零售数据集文件：

import matplotlib.pyplot as plt
import mlxtend.frequent_patterns
import mlxtend.preprocessing
import numpy
import pandas
online = pandas.read_excel(
    io="Online Retail.xlsx", 
    sheet_name="Online Retail", 
    header=0
)

清洗并准备数据进行建模，包括将清洗后的数据转化为列表的列表：

online['IsCPresent'] = (
    online['InvoiceNo']
    .astype(str)
    .apply(lambda x: 1 if x.find('C') != -1 else 0)
)
online1 = (
    online
    .loc[online["Quantity"] > 0]
    .loc[online['IsCPresent'] != 1]
    .loc[:, ["InvoiceNo", "Description"]]
    .dropna()
)
invoice_item_list = []
for num in list(set(online1.InvoiceNo.tolist())):
    tmp_df = online1.loc[online1['InvoiceNo'] == num]
    tmp_items = tmp_df.Description.tolist()
    invoice_item_list.append(tmp_items)

对数据进行编码并将其重新构建为 DataFrame：

online_encoder = mlxtend.preprocessing.TransactionEncoder()
online_encoder_array = online_encoder.fit_transform(invoice_item_list)
online_encoder_df = pandas.DataFrame(
    online_encoder_array, 
    columns=online_encoder.columns_
)
online_encoder_df.loc[
    20125:20135, 
    online_encoder_df.columns.tolist()[100:110]
]

输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_35.jpg

图 8.35：从完整在线零售数据集中构建的清洗、编码和重新构建的 DataFrame 子集

活动 19：在完整在线零售数据集上运行 Apriori 算法

解决方案：

使用合理的参数设置在完整数据上运行 Apriori 算法：
```
mod_colnames_minsupport = mlxtend.frequent_patterns.apriori(
    online_encoder_df, 
    min_support=0.01,
    use_colnames=True
)
mod_colnames_minsupport.loc[0:6]
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_36.jpg

图 8.36：使用完整在线零售数据集的 Apriori 算法结果
将结果过滤到包含 10 COLOUR SPACEBOY PEN 的物品集。将其支持度值与 练习 44 中的支持度值进行比较，执行 Apriori 算法：
```
mod_colnames_minsupport[
    mod_colnames_minsupport['itemsets'] == frozenset(
        {'10 COLOUR SPACEBOY PEN'}
    )
]
```
输出结果如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_37.jpg

图 8.37：包含 10 支 COLOUR SPACEBOY PEN 的物品集结果

支持度值确实发生了变化。当数据集扩展到包含所有交易时，该物品集的支持度从 0.015 增加到 0.015793。也就是说，在用于练习的缩小数据集中，这个物品集出现在 1.5% 的交易中，而在完整数据集中，它出现在大约 1.6% 的交易中。
添加一个新的列，包含物品集的长度。然后，过滤出长度为 2 且支持度在 [0.02, 0.021] 范围内的物品集。是否与 练习 44 中找到的物品集相同？执行 Apriori 算法，步骤 6？
```
mod_colnames_minsupport['length'] = (
    mod_colnames_minsupport['itemsets'].apply(lambda x: len(x))
)
mod_colnames_minsupport[
    (mod_colnames_minsupport['length'] == 2) & 
    (mod_colnames_minsupport['support'] >= 0.02) &
    (mod_colnames_minsupport['support'] < 0.021)
]
```
https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_38.jpg

图 8.38：基于长度和支持度过滤结果的部分

结果确实发生了变化。在查看具体的物品集及其支持度值之前，我们看到这个经过过滤的 DataFrame 比前一个练习中的 DataFrame 物品集要少。当我们使用完整数据集时，符合过滤标准的物品集较少；也就是说，只有 14 个物品集包含 2 个物品，并且支持度值大于或等于 0.02，且小于 0.021。在前一个练习中，有 17 个物品集符合这些标准。

绘制 support 值：

mod_colnames_minsupport.hist("support", grid=False, bins=30)
plt.title("Support")

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_27.jpg

图 8.39：支持度值的分布

该图展示了完整交易数据集中的支持度值分布。正如你可能已经猜测的，分布是右偏的；也就是说，大多数项集的支持度较低，并且在高端范围内有一个长尾。考虑到存在大量独特的项集，单个项集在高比例交易中出现的情况并不令人惊讶。凭借这些信息，我们可以告诉管理层，即便是最突出的项集也仅在大约 10%的交易中出现，而绝大多数项集的出现频率不到 2%。这些结果可能不会支持改变商店布局，但很可能会对定价和折扣策略提供指导。通过公式化一些关联规则，我们可以获得更多有关如何构建这些策略的信息。

活动 20：在完整在线零售数据集上查找关联规则

解决方案：

将关联规则模型拟合到完整数据集上。使用置信度指标，最小阈值为 0.6：
```
rules = mlxtend.frequent_patterns.association_rules(
    mod_colnames_minsupport, 
    metric="confidence",
    min_threshold=0.6, 
    support_only=False
)
rules.loc[0:6]
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_40.jpg

图 8.40：基于完整在线零售数据集的关联规则
计算关联规则的数量。这个数量是否与练习 45，推导关联规则，步骤 1中找到的数量不同？
```
print("Number of Associations: {}".format(rules.shape[0]))
```
存在498条关联规则。
绘制置信度与支持度的关系图：
```
rules.plot.scatter("support", "confidence", alpha=0.5, marker="*")
plt.xlabel("Support")
plt.ylabel("Confidence")
plt.title("Association Rules")
plt.show()
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_41.jpg

图 8.41：置信度与支持度的关系图

该图揭示了该数据集中一些关联规则，它们具有相对较高的支持度和置信度值。
查看提升值、杠杆值和信念值的分布：
```
rules.hist("lift", grid=False, bins=30)
plt.title("Lift")
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_42.jpg

图 8.42：提升值的分布

rules.hist("leverage", grid=False, bins=30)
plt.title("Leverage")

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_43.jpg

图 8.43：杠杆值的分布

plt.hist(
    rules[numpy.isfinite(rules['conviction'])].conviction.values, 
    bins = 30
)
plt.title("Conviction")

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_08_44.jpg

图 8.44：信念值的分布

得出关联规则后，我们可以向管理层提供更多信息，其中最重要的一点是，大约有七个项集在支持度和置信度上都有较高的值。查看支持度与置信度的散点图，看看哪些七个项集与其他项集有所区别。这七个项集的提升值也很高，从提升直方图中可以看出。看来我们已经识别出一些可以采取行动的关联规则，这些规则可以用来推动商业决策。

第九章：热点分析

活动 21：在一维空间中估算密度

解决方案：

打开一个新笔记本并安装所有必要的库。

get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
import numpy
import pandas
import seaborn
import sklearn.datasets
import sklearn.model_selection
import sklearn.neighbors
seaborn.set()

从标准正态分布中采样 1,000 个数据点。将 3.5 加到样本的最后 625 个值上（即，索引范围从 375 到 1,000）。为此，使用numpy.random.RandomState设置随机状态为 100，以保证采样值一致，然后使用randn(1000)调用随机生成数据点：
```
rand = numpy.random.RandomState(100)
vals = rand.randn(1000)  # standard normal
vals[375:] += 3.5
```
将 1,000 个数据点的样本数据绘制为直方图，并在其下方添加散点图：
```
fig, ax = plt.subplots(figsize=(14, 10))
ax.hist(vals, bins=50, density=True, label='Sampled Values')
ax.plot(vals, -0.005 - 0.01 * numpy.random.random(len(vals)), '+k', label='Individual Points')
ax.legend(loc='upper right')
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_29.jpg

图 9.29：随机样本的直方图，下方是散点图

定义一组带宽值。然后，定义并拟合一个网格搜索交叉验证算法：

bandwidths = 10 ** numpy.linspace(-1, 1, 100)
grid = sklearn.model_selection.GridSearchCV(
    estimator=sklearn.neighbors.KernelDensity(kernel="gaussian"),
    param_grid={"bandwidth": bandwidths},
    cv=10
)
grid.fit(vals[:, None])

提取最佳带宽值：

best_bandwidth = grid.best_params_["bandwidth"]
print(
    "Best Bandwidth Value: {}"
    .format(best_bandwidth)
)

重新绘制步骤 3中的直方图，并叠加估计的密度：

fig, ax = plt.subplots(figsize=(14, 10))
ax.hist(vals, bins=50, density=True, alpha=0.75, label='Sampled Values')
x_vec = numpy.linspace(-4, 8, 10000)[:, numpy.newaxis]
log_density = numpy.exp(grid.best_estimator_.score_samples(x_vec))
ax.plot(
     x_vec[:, 0], log_density, 
     '-', linewidth=4, label='Kernel = Gaussian'
)
ax.legend(loc='upper right')

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_30.jpg

图 9.30：随机样本的直方图，叠加了最佳估计的密度

活动 22：分析伦敦的犯罪情况

解决方案：

加载犯罪数据。使用保存下载目录的路径，创建年份-月份标签的列表，使用read_csv命令逐个加载单独的文件，然后将这些文件合并在一起：

base_path = (
    "~/Documents/packt/unsupervised-learning-python/"
    "lesson-9-hotspot-models/metro-jul18-dec18/"
    "{yr_mon}/{yr_mon}-metropolitan-street.csv"
)
print(base_path)
yearmon_list = [
    "2018-0" + str(i) if i <= 9 else "2018-" + str(i) 
    for i in range(7, 13)
]
print(yearmon_list)
data_yearmon_list = []
for idx, i in enumerate(yearmon_list):
    df = pandas.read_csv(
        base_path.format(yr_mon=i), 
        header=0
    )

    data_yearmon_list.append(df)

    if idx == 0:
        print("Month: {}".format(i))
        print("Dimensions: {}".format(df.shape))
        print("Head:\n{}\n".format(df.head(2)))
london = pandas.concat(data_yearmon_list)

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_31.jpg

图 9.31：单个犯罪文件的示例

该打印信息仅针对加载的第一个文件，即 2018 年 7 月伦敦大都会警察局的犯罪信息。该文件包含近 100,000 条记录。你会注意到数据集中包含大量有趣的信息，但我们将重点关注Longitude（经度）、Latitude（纬度）、Month（月份）和Crime type（犯罪类型）。

打印完整（六个月）和合并数据集的诊断信息：

print(
    "Dimensions - Full Data:\n{}\n"
    .format(london.shape)
)
print(
    "Unique Months - Full Data:\n{}\n"
    .format(london["Month"].unique())
)
print(
    "Number of Unique Crime Types - Full Data:\n{}\n"
    .format(london["Crime type"].nunique())
)
print(
    "Unique Crime Types - Full Data:\n{}\n"
    .format(london["Crime type"].unique())
)
print(
    "Count Occurrences Of Each Unique Crime Type - Full Type:\n{}\n"
    .format(london["Crime type"].value_counts())
)

输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_32.jpg

图 9.32：完整犯罪数据集的描述信息

将数据框限制为四个变量（Longitude，Latitude，Month和Crime type）：
```
london_subset = london[["Month", "Longitude", "Latitude", "Crime type"]]
london_subset.head(5)
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_33.jpg

图 9.33：数据框中犯罪数据的子集，仅保留经度、纬度、月份和犯罪类型列
使用seaborn的jointplot函数，拟合并可视化 2018 年 7 月、9 月和 12 月自行车盗窃的三种核密度估计模型：
```
crime_bicycle_jul = london_subset[
    (london_subset["Crime type"] == "Bicycle theft") & 
    (london_subset["Month"] == "2018-07")
]
seaborn.jointplot("Longitude", "Latitude", crime_bicycle_jul, kind="kde")
```
输出如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_34.jpg

图 9.34：2018 年 7 月自行车盗窃的联合密度和边际密度估计
```
crime_bicycle_sept = london_subset[
    (london_subset["Crime type"] == "Bicycle theft") & 
    (london_subset["Month"] == "2018-09")
]
seaborn.jointplot("Longitude", "Latitude", crime_bicycle_sept, kind="kde")
```
输出如下：

![图 9.35：2018 年 9 月自行车盗窃的联合分布和边际分布估计]

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_35.jpg)

图 9.35：2018 年 9 月自行车盗窃的联合分布和边际分布估计
```
crime_bicycle_dec = london_subset[
    (london_subset["Crime type"] == "Bicycle theft") & 
    (london_subset["Month"] == "2018-12")
]
seaborn.jointplot("Longitude", "Latitude", crime_bicycle_dec, kind="kde")
```
输出结果如下：

![图 9.36：2018 年 12 月自行车盗窃的联合分布和边际分布估计]

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_36.jpg)

图 9.36：2018 年 12 月自行车盗窃的联合分布和边际分布估计

从月份到月份，自行车盗窃的密度保持相当稳定。密度之间有细微的差异，这在所预期之内，因为这些估计的密度是基于三个月的样本数据。根据这些结果，警察或犯罪学家应该对预测未来最可能发生自行车盗窃的地点充满信心。
重复步骤 4；这次，使用 2018 年 8 月、10 月和 11 月的商店盗窃数据：
```
crime_shoplift_aug = london_subset[
    (london_subset["Crime type"] == "Shoplifting") & 
    (london_subset["Month"] == "2018-08")
]
seaborn.jointplot("Longitude", "Latitude", crime_shoplift_aug, kind="kde")
```
输出结果如下：

![图 9.37：2018 年 8 月商店盗窃事件的联合分布和边际分布估计]

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_37.jpg)

图 9.37：2018 年 8 月商店盗窃事件的联合分布和边际分布估计
```
crime_shoplift_oct = london_subset[
    (london_subset["Crime type"] == "Shoplifting") & 
    (london_subset["Month"] == "2018-10")
]
seaborn.jointplot("Longitude", "Latitude", crime_shoplift_oct, kind="kde")
```
输出结果如下：

![图 9.38：2018 年 10 月商店盗窃事件的联合分布和边际分布估计]

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_38.jpg)

图 9.38：2018 年 10 月商店盗窃事件的联合分布和边际分布估计
```
crime_shoplift_nov = london_subset[
    (london_subset["Crime type"] == "Shoplifting") & 
    (london_subset["Month"] == "2018-11")
]
seaborn.jointplot("Longitude", "Latitude", crime_shoplift_nov, kind="kde")
```
输出结果如下：

![图 9.39：2018 年 11 月商店盗窃事件的联合分布和边际分布估计]

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_39.jpg)

图 9.39：2018 年 11 月商店盗窃事件的联合分布和边际分布估计

与自行车盗窃的结果类似，商店盗窃的密度在各个月份之间保持相当稳定。2018 年 8 月的密度看起来与其他两个月有所不同；然而，如果你查看经纬度值，你会发现密度非常相似，只是发生了平移和缩放。这是因为可能存在一些离群点，迫使绘图区域变得更大。

重复步骤 5；这次使用 2018 年 7 月、10 月和 12 月的入室盗窃数据：

crime_burglary_jul = london_subset[
    (london_subset["Crime type"] == "Burglary") & 
    (london_subset["Month"] == "2018-07")
]
seaborn.jointplot("Longitude", "Latitude", crime_burglary_jul, kind="kde")

输出结果如下：

![图 9.40：2018 年 7 月入室盗窃的联合分布和边际分布估计]

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_40.jpg)

图 9.40：2018 年 7 月入室盗窃的联合分布和边际分布估计

crime_burglary_oct = london_subset[
    (london_subset["Crime type"] == "Burglary") & 
    (london_subset["Month"] == "2018-10")
]
seaborn.jointplot("Longitude", "Latitude", crime_burglary_oct, kind="kde")

输出结果如下：

![图 9.41：2018 年 10 月入室盗窃的联合分布和边际分布估计]

](https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/app-unspr-lrn-py/img/C12626_09_41.jpg)

图 9.41：2018 年 10 月入室盗窃的联合分布和边际分布估计

crime_burglary_dec = london_subset[
    (london_subset["Crime type"] == "Burglary") & 
    (london_subset["Month"] == "2018-12")
]
seaborn.jointplot("Longitude", "Latitude", crime_burglary_dec, kind="kde")