使用时间序列分类的金融故事
让我们教一台机器来描述市场是如何运动的。
Source: Pexels
说到叙事，金融市场特别有吸引力，因为它们只会随着时间的推移沿着一个轴移动。一个人可能会感到饥饿而快乐，但在金融领域，事情会起起落落；变强或变弱；膨胀或收缩。
作为人类，当我们看市场走势图时，我们可以很容易地得出一个市场在一段时间内表现的基本故事:
但市场用噪音弥补了这种简单。在上面的例子中，下面的三幅图像描述了一个相似的场景，第四幅描述了一个不同的场景:
对我们来说，这一切都很清楚。但是你能教计算机把这些输入转化成适当的叙述吗？
一种方法是，对于给定的财务表现时间序列，决定哪种给定的叙述最能描述它。
请注意，虽然这是一个自然语言生成(NLG)问题，但机器吐出的实际文本并不完全是由机器生成的；大部分是由人类预先写好的。
一些 NLG 从业者现在认为这种方法相当过时。这部分是因为确定选定的叙述让人感觉要么文本不够丰富，要么有太多的叙述不值得努力。(如果你不得不写 1000 种不同的叙述，这真的省去了麻烦吗？)
我们自以为是。我们可能相信人类渴望无限的叙事多样性，但即使在小说中，我们实际上也被很少几种类型的情节所吸引(七种，根据布克的说法；六，据冯内古特)。
在金融学中，最简单的叙述形式只考虑起点和终点，而忽略了中间的一切。例如:“黄金价格在一个月内上涨”。让我们称之为两点叙述。
这种形式的叙述是可以的，但是当讨论市场时，很大程度上取决于时间框架和背景。比方说，对于股票，你可能想更详细地了解一年来的涨跌情况；一个下午，不太多。以下是稍微丰富一些的例子:
现在，让我们来看看由资产管理公司施罗德发布的这些写得很好的总结要点:
请注意，第一个和第三个项目符号是 3 点记叙文，第二个和第四个是 2 点记叙文。
我们不介意叙述简单，因为句子的后半部分丰富了文本。第一部分是对发生的事情的叙述，第二部分——非常松散地——谈论为什么。
回到我们自己的问题，我们将限制自己不到 20 种叙事类型，包括 2 点和 3 点叙事。编写文本组件不会花很长时间，我们可以很容易地为相同的场景提供不同的措辞来进行混合。
我们的任务
很明显，这是一个分类的问题，因为我们问:给定一个新系列，它最像我们选择的哪种叙事类型？因此，我们将构建某种形式的分类器。
但是我们用什么来训练它呢？
我们没有任何真实世界的标签数据(也就是说，一大堆过去的性能曲线，每一条都标有最合适的叙述类型)。获取这些数据既费时又费力。因此，我们将制造贴有标签的样品，看起来非常接近真实的东西。
对于本文的剩余部分，以下是计划:
1.我们将为我们的数据生成基础:每个故事的大量样本。
2.我们将在每个数据样本中添加噪声，这样看起来就像金融市场的行为方式。
3.我们将使用现在有噪声的数据来训练分类器。
4.我们将在人造数据和真实数据上测试它，看看我们会得到什么。
1。生成叙述样本[这是容易的部分]
假设在一个月内，一个指数从 100 开始，下降，然后再次上升，以高于开始的位置结束。现在我们知道这是一个三点叙事。
但是，有各种各样的事情可以变化，同时仍然与故事保持一致:指数何时达到最低点？这个低点有多低？而适当上升又是什么意思呢？
我们的第一点总是固定的:在时间 0，我们将该值归一化为 100。但其余的可以摇摆不定，上面的问题说明了参数，把我们从一个叙事的一般描述到一个具体的实例。
因此，每种叙事类型都有几个参数，这些参数可以在一个范围内变化。为了生成大量的样本，我们从这些范围中随机选取值。比方说，对于自下而上的叙述，我们可能会得到一堆类似这样的例子:
但这些看起来一点也不像现实生活中金融市场的运行方式。在现实生活中，市场是不稳定的，所以它们会来回波动。我们需要做的是把我们整洁的叙事样本，并添加一些噪音。
2。添加噪声【这是统计部分】
我的问题是我需要添加足够多的噪声。
如果我加入的噪音太少，这些线会看起来不自然的直，完全不像真实市场的行为方式。但是如果我添加太多，或者不控制随机性，那么我华丽的叙述点可能会淹没在噪音中。
就我个人而言，我不知道这个问题的答案。
所以我求助于我丈夫的瑞士军刀，他建议我使用一种叫做布朗桥的东西。
这将是对它是什么以及它为什么起作用的一个没有痛苦的解释。出于这个目的，姑且说你严重醉酒。
酒吧刚刚关门，你的家在一英里外，你需要在日出前赶到那里。无论你多么努力想走直，你都会摇摇晃晃。当你试图移动时，有时你向前迈一步，有时后退一步。喝得这么醉，你完全不记得你之前做了什么。
一个小时后，你有什么进展吗？你追踪了一条随机路径，所以我们不能肯定。但是由于你沿着随机路径走的方式，我们知道你离起点的距离像钟形曲线一样分布，其方差取决于你的步幅和你走的步数。
又一个小时过去了。因为你没有变得更清醒，你离起点的距离再次像钟形曲线一样分布，但是现在方差是原来的两倍。就这样继续下去。
我的朋友，你正在经历布朗运动。
布朗桥是**强制结束点的布朗运动:**它只看布朗运动在特定时间(日出)结束于特定点(比如你家)的轨迹。所以，这又是一个统计过程——我们无法确切知道在开始和结束之间的任何给定时间点你在哪里。但这是一个过程，是以在特定时间具有特定价值为条件的。
在开始和结束之间，你仍然像个醉鬼一样走路。但是因为你的未来已经被预言(你将在日出时结束在家里)，你不像以前那样自由自在地漫游。你在中距离左右最自由，越接近终点，你的自由就越小。
好了，别再用喝醉的比喻了，回到我们的话题。
我们希望固定的叙述点没有噪音，这样它们就不会移动。在这两者之间，我们希望它变得嘈杂，这样我们的股票市场运动就会跳跃。布朗桥正是这样做的:
3。分类【这是机器学习部分】
因为我们的数据是一个时间序列，我们完全处于时间序列分类的世界中。
但这是一个很大的世界。为了选择一种分类方法，我们将用一把大砍刀把它不相关的部分切掉。
首先，我们删减了用于预测目的的方法(例如，根据过去的数据猜测股票未来是涨还是跌)。那不是我们的任务。
然后，有一整类算法试图学习曲线的明确特征；2016 年对这一学派的评论见[1]。
显式特征对于那些对自己的数据没有太多假设的人来说非常有用，他们希望识别关键属性，以便能够生成类似的样本，或者描述哪些品质“造就”了某种具有代表性的事物。同样，也不是我们:
首先，我们已经对数据做了很多假设。另外，在样本和真实数据中，我们认为曲线的大部分只是随机噪声。
转向神经网络方法。在这里，我们分类，但不能完全描述一个类与另一个类的区别。最方便的是，今年 3 月发表的一篇文章回顾了时间序列分类的深度学习的不同方法[2]。
它有一个方便的图表，可以帮助我们进一步削减:和以前一样，我们对学习如何创建新样本的解决方案不感兴趣，我们对特征工程也不感兴趣:
这给我们留下了端到端的方法。此外，这个名字承载着让生活变得简单的承诺，我觉得这非常有吸引力。
上述论文的结果表明，对于我们的一维问题(=所有东西都在一个轴上运动，这个轴就是时间)，在大多数数据集上表现最好的端到端方法是 ResNet。
作为一个最先进的图像分类器，这并不令人震惊。(这里有 FastAi 解释它，这里有 PyTorch 源代码。我最终用 ResNet-34 获得了足够好的结果。
4。结果【这是好玩的部分】
我几乎要让它听起来好像一切都立即工作。因为在现实世界中，任何人都只会关心:这有用吗？
起初，当然没有。但是出于非常有趣的原因。
答:关于捏造数据
当我看到最初的结果时，成功率并不大。但是很快就清楚了，对于许多曲线来说，尽管分类不同于的意图，但它的解释同样有效:
一旦我去除了噪音，重新训练了一个模型，当然，成功率激增；但这些数据不再像市场实际运行的那样。这意味着，一旦真正的噪音被加入，某些叙事类型之间的区别往往会变得模糊。
我发现这相当美妙:分类器本身往往不会出错*；它获得了我所希望的显著意义，只是以一种我认为不可信的方式。但这确实提出了一个问题:我到底应该致力于改善什么？*
有一件事几乎立刻引起了注意，那就是，仅仅以时间序列作为输入，该模型很难推断出序列开始和结束之间的关系。但这意味着失去重要的信息，万一一个故事的结尾是相对于起点来定义的。
考虑到 ResNet 的工作方式，这并不奇怪，而且通过将输入更改为双倍序列也很容易解决:
进一步的测试显示，对它来说，最难学习的事情之一是什么时候忽略噪音，看什么时候东西保持“差不多一样”。因此，解决性能问题的方法不是摆弄网络，而是更仔细地观察用于定义叙事的参数。我永远也不会完全摆脱因为能够讲述一条曲线的多个故事而产生的模棱两可。因此，我的目标是让它在捏造的数据上做得很好，但不是争取接近 100%的分类。最终，它在嘈杂的测试数据中达到了 85%以上的分类，错误很可能出现在密切相关的叙述中。
b .关于真实的过往数据
有了真实数据，一些正常化是必要的，以解释日本股市对油价的不同反应。例如，一些市场比其他市场有更多的内在差异。这些调整是利用我从 Investing.com得到的历史数据完成的。
但是，即使试图做到严格，这样的调整也不是没有问题。例如，发达市场的股票一直受到永无止境的超级宽松货币政策的支撑，这导致了长期以来创纪录的低波动性。这意味着，即使分类器的噪声参数根据“真实市场数据”进行了调整，它们仍然偏向特定的经济环境。
尽管如此，结果还是不错的。
很自然，我想给这种说法加上一个数字，但这又让我们回到了起点:获取标签数据这个非常常见的问题。
因此，我将以最近几个月真实市场数据的一个小样本作为结束:
**
结论
前面，我提到了为什么人们不喜欢将解释局限于一些固定的叙述的方法。但我从没说过是什么让它如此伟大。
我们不仅仅在金融市场上有折线图。我们到处都有。
由于叙述结构不会改变，将现有系统应用到另一个领域的成本将是编写少量的自然语言组件，并可能微调参数(特别是噪声模式)。这似乎不是一个巨大的代价。
无论如何——原谅这个蹩脚的双关语——这是一个很好的学习曲线。希望你也喜欢它！
来源:
[1] A. Bagnall，J. Lines，A. Bostrom，J. Large，e .基奥，The great time series class ification bake off:最近算法进展的回顾和实验评估 (2016)，数据挖掘和知识发现。
[2] H. Ismail Fawaz，G. Forestier，J. Weber，L. Idoumghar，和 P.-A. Muller，时间序列分类的深度学习:综述 (2019)，Arxiv。
从头开始在 20 行代码中查找人脸边缘
我在读一本机器学习的书，了解到边缘是机器学习图片中是否有物体(在这种情况下是人脸)的重要特征输入。看左边只有边缘的图形，人眼很容易分辨出是一张脸，不是吗？这同样有助于机器。
最初我认为寻找边缘本身需要一些“人工智能”。但我记得 Python PIL 库有 find_edges 过滤器，这绝对不是机器学习功能。然后，我开始思考如何找到边缘，并最终找到了一个非常简单的方法，只需要从头开始构建大约 20 行 Python 代码。
逻辑
我的逻辑很简单如下。
- 从上到下，从左到右扫描所有像素
- 如果像素与右侧或底部的相邻像素非常不同，则将其标记为边缘。
履行
我用 Python 代码实现了这一点，但是算法本身是独立于语言的。
以下是我们需要的库:
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
import math
如果你是 Python 新手，先安装 Python3 ，然后用 pip 安装需要的库。
pip install -U Pillow numpy matplotlib
首先，我们将图像读取为像素的 2D 阵列。每个像素是 3 个值[红、绿、蓝]的阵列，每个颜色值是 0 到 255，例如像素值[0，0，0]是黑色。
filename = 'old.jpg'
im = Image.open(filename)
im_array = np.asarray(im)
如果图像太大，您可以先调整其大小。
n = 2 # resize times
im = im.resize( (int(im.size[0]/n), int(im.size[1]/n)) )
其次，我们定义一个函数来度量两个像素之间的差异。我使用三个颜色值之间的均方差的平方根，就像均方根偏差( RMSD )。
difference between pixel a and b
例如，黑色像素 A 为[0，0，0]，白色像素 B 为[255，255，255]，因此 A 和 B 的平方根差为:
sqrt(([255–0]**2 + [255–0]**2 + [255–0]**2)/3) = sqrt(255**2) = 255。
注:x**2 表示 x*x。
这里有一个窍门。我们需要首先将像素值转换为 int 类型，即int(a[i])，以进行减法运算，因为像素值为 ubyte[0–255]，减法运算可能会变为负值，从而导致类型溢出问题。
然后我们从上到下，从左到右迭代图像数组。如果像素与其右侧或底部相邻像素之间的平方根差大于预定义的阈值，则将其标记为边缘像素，并将其设置为黑色[0，0，0]，否则将其设置为白色[255，255，255]作为背景。
main part
您可以微调图像的阈值diff_threshold。10 是我用来创建顶部图像的值。增加较少边的阈值。例如，通过增加到 20，我得到了一个新的图像如下。
最后，使用 matplotlib 将新的边缘图像保存到本地。
plt.imshow(new_array)
plt.savefig('new.jpg')
如果您想在弹出窗口中显示新图像，请使用下面的代码。
plt.imshow(new_array)
plt.show()
你可以在这里找到组合的完整代码。
笔记
请注意，这种算法不是专门针对人脸检测的。它可以找到所有图像的边缘。你自己试试。
感谢阅读。如有任何疑问，您可以通过 Linkedin 联系我。
仅用空间查找最低公共祖先子树和最短依赖路径
使用空间作为一个解决所有问题的工具
Photo by Paula May on Unsplash
在之前的帖子:如何用 spaCy 和 StanfordNLP 寻找最短依赖路径中，我讲过如何用 spaCy 和 NetworkX 提取最短依赖路径(SDP)。
但是使用 NetworkX 有个问题。我们无法获得头部实体或尾部实体的索引。例如，我们有下面的句子。三联是Convulsions->(caused_by)->fever。但是这句话里有两个fever。
Convulsions that occur after DTaP are caused by a fever, and fever may cause headache.
NetworkX 解决方案
一个解决方案是为每个令牌添加一个索引，并指定要查找哪个fever。
import spacy
import networkx as nx
nlp = spacy.load("en_core_web_sm")
doc = nlp(u'Convulsions that occur after DTaP are caused by a fever, and fever may cause headache.')# Add pair to edges
edges = []
for token in doc:
for child in token.children:
**edges.append(('{0}-{1}'.format(token.text, token.i),
'{0}-{1}'.format(child.text, child.i)))**# Construct Graph with nextworkx
graph = nx.Graph(edges)# Get the length and path
**entity1 = 'Convulsions-0'
entity2 = 'fever-9'**print(nx.shortest_path_length(graph, source=entity1, target=entity2))
print(nx.shortest_path(graph, source=entity1, target=entity2))##### output #####
3
['Convulsions-0', 'caused-6', 'by-7', 'fever-9']
下面的edges看起来像。
In [6]: edges
Out[6]:
[('Convulsions-0', 'occur-2'),
('occur-2', 'that-1'),
('occur-2', 'after-3'),
('caused-6', 'Convulsions-0'),
('caused-6', 'DTaP-4'),
('caused-6', 'are-5'),
('caused-6', 'by-7'),
('caused-6', ',-10'),
('caused-6', 'and-11'),
('caused-6', 'cause-14'),
('by-7', 'fever-9'),
('fever-9', 'a-8'),
('cause-14', 'fever-12'),
('cause-14', 'may-13'),
('cause-14', 'headache-15'),
('cause-14', '.-16')]
这样，我们可以确保尾部实体令牌是fever-9而不是fever-12。
这个解决方案有点麻烦，因为 NetworkX 只接受字符串类型，我们必须在字符串中包含这样的信息。
只有空间的解决方案
sapCy 中的令牌类非常强大。它在每个标记中都有索引信息。但是怎么用 spaCy 找 SDP 呢？
经过一些研究，我发现我们可以利用 get_lca_matrix 函数。
In [11]: doc = nlp(u"This is a test")
...: lca_matrix = doc.get_lca_matrix()In [12]: lca_matrix
Out[12]:
array([[0, 1, 1, 1],
[1, 1, 1, 1],
[1, 1, 2, 3],
[1, 1, 3, 3]], dtype=int32)
Doc.get_lca_matrix:计算给定
*Doc*的最低公共祖先(lca)矩阵。返回包含祖先整数索引的 LCA 矩阵，如果没有找到共同祖先，则返回*-1*。
我们可以使用这个函数来查找 SDP。
import spacy
nlp = spacy.load("en_core_web_sm")doc = nlp(u'Convulsions that occur after DTaP are caused by a fever, and fever may cause headache.')**def get_sdp_path(doc, subj, obj, lca_matrix):
** lca = lca_matrix[subj, obj]
current_node = doc[subj]
subj_path = [current_node]
if lca != -1:
if lca != subj:
while current_node.head.i != lca:
current_node = current_node.head
subj_path.append(current_node)
subj_path.append(current_node.head)current_node = doc[obj]
obj_path = [current_node]
if lca != -1:
if lca != obj:
while current_node.head.i != lca:
current_node = current_node.head
obj_path.append(current_node)
obj_path.append(current_node.head)
return subj_path + obj_path[::-1][1:]
# set head entity index and tail entity index
**head = 0
tail = 9**
**sdp = get_sdp_path(doc, head, tail, doc.get_lca_matrix())** print(sdp)##### output #####
[Convulsions, caused, by, fever]
get_sdp_path()可以找到头实体和尾实体之间的 SDP。我们唯一需要的是输入头部实体索引和尾部实体索引。
在get_sdp_path()函数中，它实际上首先找到从头实体到 LCA 令牌的 SDP 路径，然后找到从尾实体到 LCA 令牌的 SDP 路径。最后，我们将两个子树组合在一起并返回结果。
查看我的其他帖子 中 同 一个分类查看 ！
GitHub:bramble Xu LinkedIn:徐亮 博客:bramble Xu
参考
使用微软 NNI 查找 scikit-learn 管道的超参数
Source: https://www.pexels.com
什么是超参数搜索？
任何机器学习算法的目标都是在样本数据中找到模式，然后使用这些模式来预测给定的未知样本(测试)数据的结果。
机器学习算法学习的模式被编码在模型的权重(也称为参数)中。然而，机器学习算法可以应用于不同类型的数据，即使对于数据“类型”相同的情况，数据的分布也可能不同。
一个具体的例子，让我们看看 scikit-learn 的支持向量分类器(SVC)的 API。
Source : Screenshot of https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html
根据上图，您可以理解为不同类型和分布的输入数据集找到 C、核、度、Gamma 和 Coef0 的适当组合并不是一项简单的任务。
这些类型的参数(机器学习算法的输入，即上述示例中的 C、内核等)被称为超参数，找到它们的任务被称为超参数搜索或调整。
什么是微软 NNI？
微软 NNI(神经网络智能)是一个开源的神经架构搜索和超参数调整的 AutoML 工具包。
我了解超参数调优，但什么是 AutoML？
实际上，在许多方面，超参数搜索/调整是一个叫做 AutoML 的大领域的一个子领域。也就是说，根据你在哪里阅读 AutoML，你可能会得到不同的印象。在这里，我提供了通常与 AutoML 相关的三个主要特性
**超参数搜索。**我的神经网络有多少层、神经元、宽度和深度？我的 SVM 的 C & gamma 值应该是多少？
网络架构搜索(NAS) 。神经网络中的哪种层组合可以产生最佳性能。例如，使用 NAS 发现了名为 EfficientNet 的最新领先 ImageNet 网络。我把这个定义过于简化了，但是希望你能明白其中的要点！
基于云的 AutoML 。云提供商公开了一项服务，并告诉我提供您的训练数据集并指定任务(分类或回归)，我将为您构建模型。回家睡觉吧！。这个定义意味着一个典型的机器学习管道的大部分(如果不是全部的话)方面(数据规范化、标准化、找到合适的算法、优化等)都是为你考虑的。这是一个自动化的机器学习。
还可以看到，上述特征很容易相互重叠。例如，有人可能认为网络架构搜索将/应该主要包括超参数搜索。
要进行超参数和架构搜索，一个简单的方法是尝试大量的组合，看看哪一个胜出！….但这不会有成效，不是吗？。这就是为什么有算法和技术的研究，为其他机器学习算法找到超参数。
这些超参数和架构搜索算法通常被称为 调谐器 ，因为大多数时候它们从一组参数开始，然后根据看到的结果选择下一组参数。简而言之，他们反复调整它们。
微软 NNI 被称为工具包，因为
- 它已经实现了许多已知的艺术调谐器，并为您编写自己的框架。
- 它为您的实验提供编写、运行、记录和可视化支持。
- 您不仅可以在本地机器上运行您的实验，而且它还可以与各种编排系统(如 Kubeflow、OpenPAI 等)很好地集成。还有足够的扩展点来集成新的或专有的机器学习平台。
- 它支持大多数机器学习框架(Tensorflow、PyTorch、MXNET、Chainer 等)以及库(scikit-learn、xgboost、lightgbm 等)。
因此，它不仅提供了运行和可视化的工具，而且还可以扩展(因此是框架的概念)进行定制。
编写和运行一个简单的 NNI 实验(试用版)
尽管 GitHub repo 列出了支持的框架(Tensorflow、scikit-learn 等)，但对于大多数用例，微软 NNI 并不知道您正在使用的框架，甚至是您的培训模型。
下图显示了撰写典型试验所涉及的 3 种类型的文件:
Source : Screenshot of a trial example from https://github.com/microsoft/nni
Source: Screenshot of search-space.json
search-space.json该文件用于指定您正在搜索的超参数的范围！
这是您最终要编写的一种示例代码。
在这个例子中，您有一个 SVC 模型。进行试验包括 4 个步骤:
1)你加载你的数据。
2)你请求 NNI 给你试验参数。
3)将参数应用到正在训练的模型
4)你回分数。
基于config.yml文件中的optimize_mode(见下文)，调谐器将通过建议下一组模型的超参数来帮助最大化或最小化分数。
这里是一个例子config.yml的文件，它被用来指定这个试验的配置选项。即在哪里运行、应该运行多长时间以及使用哪个调谐器。所以让我们看看 config.yml 文件。
Source: Screenshot of example config.yml
在配置文件中，您还可以指定到search_space.json和main.py.的路径
运行实验的时间:
nnictl create --config <path_to>/config.yml
在执行上述命令后，您将在终端上看到如下所示的响应:
INFO: Starting restful server...
INFO: Successfully started Restful server!
INFO: Setting local config...
INFO: Successfully set local config!
INFO: Starting experiment...
INFO: Successfully started experiment!
--------------------------------------------------------------------
The experiment id is egchD4qy
The Web UI urls are: http://223.255.255.1:8080 http://127.0.0.1:8080
--------------------------------------------------------------------You can use these commands to get more information about the experiment
--------------------------------------------------------------------
commands description
1\. nnictl experiment show show the information of experiments
2\. nnictl trial ls list all of trial jobs
3\. nnictl top monitor the status of running experiments
4\. nnictl log stderr show stderr log content
5\. nnictl log stdout show stdout log content
6\. nnictl stop stop an experiment
7\. nnictl trial kill kill a trial job by id
8\. nnictl --help get help information about nnictl
--------------------------------------------------------------------
你现在可以用你的浏览器来查看 NNI 进行的各种实验。每个试验本质上对应于炒作参数的一个组合。
还可以选择使用命令(如上所示)来获取终端本身的实验和试验状态。命令应该是不言自明的。
sci kit-学习管道
scikit-learn 是任何机器学习工程师最重要的库之一。它与许多优秀算法的实现捆绑在一起。scikit-learn 的一个我非常喜欢的特性是管道的概念。
管道旨在反映典型的机器学习实验，该实验由对输入特征的各种变换步骤组成，并最终以估计器(分类器或回归器)终止。转换器和估算器被称为流水线的组件。
在上图中，组件是缩放、降维、学习算法和预测模型。看看这些组件如何实现fit, transform and predict方法。
如果您与接口兼容，您还可以实现一个定制的转换器或估计器。
以下是如何使用 python 编写管道的示例:
Source: Author
上图中显示了两条管道-
a)make _ svc _ pipeline-标准化数据，然后应用 SVC
b) make_pca_svc_pipeline —规格化 pca(用于 dim 减少)和最终 svc
请注意，在这两个管道中，我是如何将所有的超参数(pca_components & svc_hparams) 作为函数的输入参数传递的。这样，我现在可以从 NNI 调谐器获得这些，并把它们传递给函数。
下面是 NNI 网络用户界面的一个标签的截图，显示了一个调谐器最后尝试的各种组合。
非常酷的可视化！我现在想编写由其他分类器组成的管道，例如，不使用 SVC，我想尝试 xgboost、lightgbm 等，但是我开始发现对于我的任务来说，这个过程是相当重复的，即对于每个实验编写config.yml, search-space.json and pipelines + entrypoint script to get the parameters from NNI Tuner, return the score and invoke the pipelines.
使用微软 NNI 进一步简化寻找 scikit-learn 管道的超级参数的任务不是很棒吗？
这就是我最终构建的东西，一个使用自省和配置文件来自动化上述步骤的小包。目前，我支持分类(很快也会支持回归)算法。
如果您对如何使用该软件包感兴趣，请继续阅读本文的其余部分。
scikit-nni:一个小 python 包，用于自动搜索 scikit-learn 管道的超参数
pip install scikit-nni
sci kit-nni(【https://github.com/ksachdeva/scikit-nni】T2)是一个辅助工具(和一个软件包)，它:
- 生成 NNI 所需的配置(config.yml & search-space.json)
- 根据您的规范自动构建 scikit-learn 管道，并成为微软 NNI 运行的实验/试用代码
scikit-nni 的特点
- 使用微软 NNI 搜索 scikit-learn 管道的超参数
- 定义管道不需要代码
- 内置数据源读取器，用于读取 npz 文件进行分类
- 支持使用自定义数据源读取器
- 定义 NNI 配置和搜索空间的单一配置文件
该软件包及其使用方法在 https://github.com/ksachdeva/scikit-nni有很好的文档记录，所以请查看完整的文档和示例。
简而言之，这是你最终会做的事-
在使用微软 NNI 时，你不用写 3 个(或更多)文件，你只需写一个配置文件。一个例子是-
看看search-space.json 和config.yaml是如何嵌入到这个配置文件中的。他们遵循与微软 NNI 相同的语法和惯例，所以我没有发明任何新的语法。最大的区别在于，您不是编写 python 代码，而是在配置文件本身中指定 scikit-learn 管道。阅读上述 YAML 文件中的注释。
你会意识到，这种自动化不仅有助于找到超参数，而且现在很容易尝试各种管道，而不必经历用代码编写它们、从 NNI 检索参数、发回分数和编写入口点脚本的麻烦。
一旦您编写了这个配置文件，然后发出以下命令:
sknni generate-experiment --spec example/pca_svc.nni.yml --output-dir experiments
像往常一样运行您的 NNI 实验，即现在发出以下命令
nnictl create --config experiments/pca-svc-classification/config.yml
然而，您需要做的是编写一个返回训练和测试数据的 callable。在存储库的文档中，我展示了如何根据项目的需求编写这样一个数据源。在不久的将来，我将为分类和回归需求捆绑一些常见的数据源。
将来会添加更多的数据源，但是您可以根据您的项目需求编写一个数据源。定制数据源的一个简单例子是
结束语
- 微软 AutoML 是一个真正设计良好的，可扩展的和文档化的解决方案
- 感谢 NNI，使用最先进的调谐器现在只是一个简单的配置问题，你不必狩猎和学习如何使用各种调谐器算法。
- scikit-nni 利用微软 nni 来自动执行超参数搜索任务，而不需要在 python 代码中定义管道和实验。它还通过使用一致的语法将文件数量减少到 1 个来简化配置文件。
scikit-nni 非常新，请随时提交您可能遇到的任何问题(https://github.com/ksachdeva/scikit-nni)以及进一步改进它的想法。贡献总是被感激的！。
找到重要的差异
确定波士顿最昂贵的社区
A view of Back Bay, Boston
你运行一个分析，你会得到这样的结果:
平均值实际上是不同的吗？当有方法将统计的严谨性放在发现的后面时，不应该对“这看起来足够不同”做出决定。结果可能在两个方向上都令人惊讶——微小的差异可能具有重大意义；巨大的差异可以缺少它。
当我在我正在从事的 MSBA 项目中了解到这些概念时，我认为做一个案例研究来解释这些概念会很有趣。
在 Github 上的项目代号是 这里是 。
波士顿的街区都一样贵吗？
我以前写的一篇文章试图预测我居住的城市波士顿的房价。提高性能的一个关键因素是为属性所在的邻域添加一个变量。
如果邻域很重要，这意味着我们应该能够首先从统计上证明平均值之间存在差异，然后进入下一个层次，了解邻域之间的比较。
数据探索
为简单起见，我将提取最昂贵且位于市中心的街区:芬威/肯莫尔、后湾、比肯山、西区、北端、市中心和南端:
我还通过按邻域移除值的顶部和底部的 10%来将值修剪到中间的 80%,以考虑异常值和极值。由于分析是基于受异常值影响的平均值进行的，我们希望尽可能创建最具代表性的数据集。我们将分析每平方英尺的价格，作为最昂贵的替代。
让我们先来看看数据的分布:
我们看到不同群体之间存在一些差异。比肯山的中间值最高，但市区的第 75 个百分点高于比肯山。西区、北端、南端和肯莫尔都有类似的分布，但在第 25 和第 75 百分位之间有不同程度的分离。哪些是真正不同的？你能自信地说吗？
当看平均值时，我们看到比箱线图多一点的变化。视觉上，我们会说比肯山最贵，北端最便宜。这是真的吗？
确定差异
ANOVA 或方差分析衡量的是数据中一系列组的平均值相同的可能性。它通过测量“组间方差”和“组内方差”来做到这一点。
组间方差
组间方差和:每组的平均值减去总数据集平均值的平方，再除以自由度。
让我们把它翻译成英语。修整数据集中所有值的平均值为 788.55 美元。按邻域划分的平均值如上图所示。对于每一个平均值，我们减去 788.55 美元，对结果求平方，然后求和。自由度是样本平均数(7)减 1。
组内方差
简单地说，对于每个组，这是组的数量(7)减去 1 倍组的方差。然后，我们将这些值相加，除以观察总数(25，028)减去样本平均数(7)。
方差分析
与 ANOVA 相关的统计值 F 值是组间方差的除以组内方差的*，并创建一个分布，我们根据该分布计算所有组具有相同均值的可能性(p 值)😗
解释上述内容的方法是，如果 p 值(Pr(>F))小于我们的显著性水平(通常为 0.05)，我们“拒绝零假设”，即所有平均值都相等，并得出组间平均值不同的结论。
这并没有告诉我们哪些群体是不同的，只是说有差异。
这个发现给我留下了疑问。主要是，哪些不一样？
成对比较和 Tukey HSD
成对比较和 Tukey HSD(“真正的显著差异”…是的，确实是这个名字)是一种查看特定值对并确定其差异显著性的方法。
再说一次，让我们把它变成现实。笔架山比市区贵吗？肯莫尔比北端贵吗？
下面左边的成对表的解读方式是对于每一对(即左上角的北端和西端)，均值相等的概率是多少？较低的值意味着平均值不太可能相等。
对于 Tukey HSD，它计算平均值的差异并围绕该差异创建一个置信区间，查看 95%置信区间(“lwr”和“upr”)是否包含范围内的 0。如果有，那么平均值没有显著差异，因为 0-表示没有差异-是一个看似合理的值。否则，差异是显著的。如果“diff”为负，则第一个邻域比第二个更便宜，反之则为正。
Pairwise Comparison (left) and Tukey HSD (right)
我更喜欢视觉，所以我们来看一下图中的 Tukey HSD:
但这意味着什么呢？
比肯山看起来比市中心更贵。现在，我们可以说它是。95%置信区间不包含 0，差异具有统计学意义。我们现在有把握地知道，在 2018 年，每平方英尺的评估值在灯塔山最高，其次是市中心，然后是后湾。
事实上，在我们的样本社区中，唯一没有统计学意义的一对是肯莫尔和北端；我们不能断定这两个街区的平均值在每平方英尺的水平上是不同的。
业务含义和结论
方差分析和相关统计在商业中的应用非常广泛。例如，如果我们关注市场中不同公司的表现，我们想问的第一个问题是:
- 是否有机会在市场中脱颖而出？
如果所有平均值在统计上是相同的，那么公司不得不超越业绩来真正影响客户赢得率(服务、产品范围等)。).换句话说，如果产品性能是赢得合同的基准标准，那么下一级投资资本应该投向哪里才能真正实现差异化？
一旦我们确定了差异化的机会所在，我们就想知道:
- 公司差异化在哪里？
这时，我们可以进行成对比较或 Tukey HSD 来确定公司之间的差异。文章开头图表中 A 组的 5.95 分和 C 组的 5.83 分实际上有区别吗？
结果是在给定的市场中理解:
- 哪里有差异化的机会？
- 如今的公司差异化在哪里？
- 公司应该如何分配资本(M&A、R&D 等)。)是在差异化标准上加倍努力，还是加快对表现不佳的改进？
我坚信学习的实际应用，虽然对我每天都知道和经历的社区进行分析很有趣，但对商业的影响以及这些可以做什么的探索更重要也更有趣。
递归程序设计
如何通过假装你已经解决了问题来解决问题
尽管递归的概念经常在大多数项目的早期就被引入到编程中，但在第一次遇到它时，它可能看起来很奇怪，并且可能会令人不快。这看起来几乎是自相矛盾的:我们如何能够用一个问题的解决方案找到另一个问题的解决方案？
Recursion can be a bit of a headache
对于那些试图掌握递归概念的人来说，我经常觉得首先认识到递归不仅仅是一种编程实践是有益的——它是一种解决问题的哲学，适用于可以处理和部分解决的问题，而问题的其余部分保持同样的形式，但在某些方面更容易或更小。这不仅仅适用于编程中的函数；我们可以用递归构造简单的日常问题。举个例子，就拿我来说，写这篇文章:假设我想写 1000 字左右，如果我的目标是每次打开都写 100 字，那么第一次写 100 字，给自己留 900 字。下一次，我写 100 字，只剩 800 了。我可以继续这样写，直到我剩下 0 个字。每次，我都部分地解决了问题，剩下的问题正在减少。
撰写我的帖子的代码可能是这样的:
write_words(words_left):
if words left > 0:
write_100_words()
words_left = words_left - 100
write_words(words_left)
我也可以迭代地实现这个算法:
write_words(words_left):
while words_left > 0:
write_100_words()
words_left = words_left - 100
如果你用任一实现完成函数调用write_words(1000),你会发现它们有完全相同的行为。事实上，每一个我们可以用递归解决的问题，我们也可以用迭代来解决(for和while循环)。那么我们为什么会选择使用递归呢？
Me writing this post recursively
为什么是递归？
信不信由你，一旦我们掌握了它，有些问题用递归比用迭代更容易解决。有时候递归效率更高，有时候可读性更强；有时递归既不更快也不更可读，但实现起来更快。有些数据结构，如树，非常适合递归算法。甚至有一些编程语言没有循环的概念——纯函数式语言如 Haskell 完全依赖递归来迭代解决问题。要点很简单:你不必理解递归就能成为一名程序员，但是你必须理解递归才能成为一名优秀的 T4 程序员。事实上，我甚至可以说理解递归是成为一个好的问题解决者的一部分，除了编程！
递归的本质
一般来说，使用递归，我们试图将一个更复杂的问题分解成一个简单的解决步骤和一个更简单的相同问题的剩余部分。然后，我们可以重复这个过程，每次都朝着解决方案迈出相同的一步，直到我们用一个非常简单的解决方案解决我们的问题(称为基础案例)。我们基础案例的简单解决方案与我们实现目标的步骤聚合在一起，然后形成我们原始问题的解决方案。
We can solve P by breaking the problem into a step towards the solution and a remaining smaller problem of the same form as the original, until we reach a simple solution to a small problem (a base case).
假设我们得到一些数据类型的实际数据，称之为 dₒ.递归的思想是假装我们已经解决了问题，或者根据我们需要定义的某种难度度计算了比 dₒ简单的这种数据类型的所有形式的期望函数f。然后，如果我们能找到一种方法，用一个或多个 f(d) s 来表达 f(dₒ) ，其中所有这些 d s 都比 dₒ 更容易(具有更小的度),那么我们就找到了一种方法来简化和求解 f(dₒ) 。我们重复这个过程，希望在某个时候，剩下的 f(d) s 会变得如此简单，以至于我们可以轻松地对它们实现一个固定的、封闭的解决方案。然后，当我们对越来越简单的问题的解决方案聚集并级联到顶部时，我们对原始问题的解决方案就会显现出来。
在上面写这个帖子的例子中，数据是这个文档中包含的等待写的文本，难度的度是文档的长度。这是一个有点做作的例子，但假设我已经解决了如何写 900 字的问题 f(900) ，那么解决 f(1000) 我需要做的就是写 100 字，然后执行我的 900 字的解决方案， f(900) 。
一个更好的例子是考虑斐波那契数列，其中第一个斐波那契数列是 0，第二个是 1，第个 ᵗʰ斐波那契数列等于前两个的和。假设我们有一个斐波纳契函数，它告诉我们 n ᵗʰ斐波纳契数:
fib(n):
if n == 0:
return 0
if n == 1:
return 1
else:
return fib(n-1) + fib(n-2)
这个函数的执行是什么样子的？让我们试试fib(4):
Visualisation of a recursion tree showing the recursive computation that leads to fib(4) = 3, notice that computation is performed depth-first.
递归解决问题时采用的一个有用的口头禅是’假装它’直到你成功’，也就是说，假装你已经解决了一个简单情况的问题，然后尝试减少更大的问题，以使用这个更简单情况的解决方案。如果一个问题适合递归，实际上应该只有少量的简单情况需要你显式地解决，也就是说，这种化简为简单问题的方法可以用来解决所有其他情况。这在斐波纳契例子fib中举例说明，其中定义fib(n)我们只是像已经计算过fib(n-1)和fib(n-2)一样行动，正如我们所希望的，这级联并减少问题到逐渐更简单的情况，直到我们到达fib(0)和fib(1)，它们有固定和简单的解决方案。
Fake it ’til you make it!
递归策略
递归有点微妙，实际上取决于你试图解决什么问题。然而，我们可以想出一些大致的步骤，这些步骤或多或少可以把我们引向正确的方向。这一战略包含三个步骤:
- 订购您的数据
- 解决小案子
- 破大案
正如我以前说过的，我认为在我们学习的时候举一个例子是有用的，但是记住递归是依赖于问题的，所以在这里试着把重点放在一般的原则上。我们将使用反转字符串的简单示例，即我们希望编写函数reverse使得reverse('Hello world') = 'dlrow olleH'。我建议回去看看这些步骤是如何应用于斐波那契函数的，然后在其他一些例子中尝试一下(网上有很多练习)。
订购您的数据
这一步绝对是开始递归解决问题的关键，然而它经常被忽视或被隐含地执行。无论我们在操作什么数据，无论是数字、字符串、列表、二叉树还是人，都有必要明确地找到一个合适的排序，给我们一个方向，让问题变得更小。这种排序完全取决于问题，但一个好的开始是考虑明显的排序:数字有自己的排序，字符串和列表可以按长度排序，二叉树可以按深度排序，人可以按无限多种合理的方式排序，例如身高、体重或在组织中的排名。如前所述，这种排序应该与我们试图解决的问题的难度相对应。
Order that data, yee-haw!
一旦我们对数据进行了排序，我们就可以认为它是我们可以减少的东西。事实上，我们可以把我们的排序写成一个序列:
0 ， 1 ， 2 ，…， n 为整数(即对于整数数据 d ，度(d) = d )
[]，[■]，[■，■]，…，[■，… , ■]用于列表
(注意 len = 0，len = 1，…，len = n，即用于列表数据 d ，度(d) = len(d) )
从右到左，我们通过一般(“大”)情况，到基本(“小”)情况。对于我们的reverse例子，我们正在操作一个字符串，我们可以选择字符串的长度作为问题的排序或度。
解决小案子
这通常是容易的部分。一旦我们有了正确的排序，我们需要查看排序中的最小元素，并决定如何处理它们。通常有一个显而易见的解决方案:在reverse(s)的情况下，一旦我们到了len(s) == 0并且我们有了reverse('')，那么我们就知道我们的答案，因为反转空字符串不会做任何事情，也就是说，我们只是返回空字符串，因为我们没有要移动的字符。一旦我们解决了我们的基本情况，并且我们知道我们的排序，那么解决一般情况就像以这样一种方式减少问题一样简单，即我们正在操作的数据的度向基本情况移动。我们需要小心，不要遗漏任何小的情况:它们被称为基本情况的原因是因为它们覆盖了排序的基础——在更复杂的递归问题中，遗漏基本情况是常见的，因此归约步骤越过了我们排序的合理末端，并开始对无意义的数据进行操作，或者导致错误。
破大案
这里，我们在排序中向右处理数据，也就是高度数据。通常，我们考虑任意次的数据，并旨在通过将其简化为包含较小次的相同问题的表达式来找到解决问题的方法，例如，在我们的 Fibonacci 示例中，我们从任意 n 开始，并将fib(n) 简化为fib(n-1) + fib(n-2)，这是一个包含我们开始时的问题的两个实例的表达式，具有较小的次(n-1 和 n-2
Big case?
当谈到reverse时，我们可以考虑长度为 n 的任意字符串，我们可以假设我们的reverse函数对所有长度小于 *n 的字符串都有效。*我们如何利用这一点来解决长度为 n 的字符串的问题呢？嗯，我们可以把字符串反过来，除了最后一个字符，然后把最后一个字符放在前面。在代码中:
reverse(string) = reverse(string[-1]) + reverse(string[:-1])
其中string[-1]对应最后一个字符，string[:-1] 对应没有最后一个字符的字符串(这些都是 pythonisms)。最后一个reverse(string[:-1])项是我们的原始问题，但是是在一个长度为 n-1 的字符串上操作的，也就是说，我们已经将我们的原始问题表述为向解决方案迈进了一步，并结合了相同的降阶问题。
将我们的reverse函数的解放在一起，我们得到如下结果:
reverse(string):
if len(string) == 0:
return ''
else:
return string[-1] + reverse(string[:-1])
Visualisation of recursive function reverse operating on some sample data.
通常需要考虑不止一种递归情况，因为给定数据类型的数据可以采用稍微不同的形式，但这完全取决于问题。举个例子，如果我们想展平一个项目列表，其中一些项目本身就是列表，我们需要区分从列表中取出的项目是单个项目还是子列表，这至少会导致两种递归情况。
最终提示
提高递归的唯一真正方法是练习。在网上看看成千上万的递归问题中的一些，或者挑战自己提出你认为可能适合递归的问题。一旦你不可避免地掌握了递归的窍门，记住如果你发现自己很难递归地解决问题，那么尝试迭代。除了学习成为更好的程序员，递归还是一种解决问题的方法，让你的生活变得更容易。如果一个问题不适合递归，它就是不适合递归；当你花更多的时间处理那些适合递归或迭代方法的问题时，你会对此有一种感觉。
有时在更困难的递归问题中，我们上面看到的策略中的第 2 步和第 3 步采取了更循环的反馈循环过程的形式。如果你不能快速找到问题的整体解决方案，最好的过程是解决你能想到的递归/'大’案例，解决你能想到的基本/'小’案例，然后看看你的方法如何在不同的数据片段上突破。这将揭示任何缺失的基础和递归案例，或者任何相互之间交互不良并需要重新思考的案例。
最后，回想一下，知道您的排序是解决递归问题的最重要的步骤，您的目标始终是涵盖这种排序的向右(递归)和最左(基本)情况，以解决给定类型的所有数据的问题。
这是一个总结——感谢阅读！
如果你喜欢这篇递归介绍，请随时联系我( 汤姆·格里格 )，对未来的博文有任何想法、疑问或建议！
现在，回到我的数据科学帖子，敬请关注！
使用 Trie(前缀树)查找字符串中的所有单词
如何使用 Trie 来存储前缀，以便在字符列表中查找句子。
A Trie example, Wikipedia Commons. https://en.wikipedia.org/wiki/Trie
找到字符列表中的所有单词是一个有趣的挑战，可能有很多方法可以做到。例如，我们以字符串“thedogatethecat”为例，它可以通过两种方式进行解析。第一个是“那只狗吃了猫”&第二个是“那只做门猫”，两个版本在语法分析上都是正确的，但其中只有一个对我们人类是正确的。
下面的解决方案将找到两个版本，并不会试图找出哪一个在语法等方面更正确。
这个想法是使用英语中的所有单词创建一个前缀树或 Trie 树，然后逐个字母地找出 Trie 树中是否有单词。很明显，有一些关于复杂性的问题，但是使用 Trie，我们将我们的词汇表中的搜索复杂性降低到 O(M ),其中 M 是您要搜索的文本的长度，而不是搜索 O(N ),其中 N 是词汇表的长度。值得注意的是，Trie 的空间复杂度更高，如果你对它感兴趣，可以在 Google 上搜索这个话题。
然而，最终的时间复杂度略高于 O(M ),因为我们是逐字母递归搜索，所以我把计算留给读者作为练习。
在运行下面的代码之前，请通过“pip install pytrie”来安装 PyTrie 。我们用字符串中所有可能的单词填充 Trie，显然，我们可以添加额外的单词，但这不会增加这个例子的价值。对于真实世界的例子，请使用适当的英语词典。我们一个字母一个字母地迭代，直到找到一个有效的前缀，递归地发送字符串的其余部分，继续搜索，在返回的路上收集所有结果。
输出应该很清楚，第一个例子是测试 Trie 中的’ t ‘，’ th ‘，然后是’ The ‘，一旦找到，字符串的其余部分将从单词’ the '的末尾开始递归处理。
我们可以看到，我们要找的两个句子被算法找到了。在收集了所有有效的句子之后，我们大概可以使用一个英语单词嵌入模型来看看我们的句子在全球范围内是否有意义。
Ori Cohen 博士拥有计算机科学博士学位，专注于机器学习。他领导着 Zencity.io 的研究团队，试图积极影响市民的生活。
寻找贝叶斯乐高——第二部分
一次拯救一只脚
Photo by Markus Spiske on Unsplash
上一集《寻找贝氏乐高积木》，我和我亲爱的朋友乔试图找到一种方法来帮助乔树立信心，让他相信他已经捡起了他的孩子们留下的所有乐高积木。我的第一个建议是，我们通过对乔进行实验来确定他在任何给定的房间范围内捡一个乐高积木的概率，从而确定乔在捡乐高积木方面有多好。乔说，“不，继续这个计划！我真的不喜欢拿起乐高玩具。你就不能找个更好的办法吗？”
“告诉我更多关于你如何收拾孩子的事.”
乔向我解释说，他通常会在孩子们上床睡觉后再去拿乐高玩具。乔喜欢每天晚上读点书或看看电视来放松一下。在这本书的两章之间，或者在电视广告中，他会拿起乐高玩具快速扫一扫。第一次扫描似乎获得了大部分的乐高积木，每次扫描都获得了前一次扫描中遗漏的部分。他一直这样做，直到他收工睡觉。
“乔，我想我们可以这样做。”
下次乔拿起乐高玩具时，我会让他做一点小工作。我让他数一数每次他捡了多少个并记录下来。他很不情愿地同意了，但是任何可以帮助他获得信心的事情，当他第二天早上起来走向咖啡机的时候，他都不会踩到乐高积木。
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')np.random.seed(42) # It's nice to replicate even virtual experimentsimport pymc3 as pm
import scipy.stats as stats
print('Running on PyMC3 v{}'.format(pm.__version__))
> Running on PyMC3 v3.6
乔汇报道
几天后，我收到乔的电子邮件:他已经按照要求做了。Joe 报告说，他扫了六次，每次都捡了以下数量的乐高积木:
Removed[1:]
> array([135, 48, 11, 4, 3, 0])
那天晚上，我邀请乔过来吃晚饭，做一些贝叶斯分析。
与 PyMC 共进晚餐
在吃完我们的羽衣甘蓝、鳄梨和黑豆玉米煎饼碗的晚餐后，乔和我舒适地围着大屏幕进行分析。
我向乔解释说，每次他打扫房间时，就好像每个乐高都有可能被捡起来。我们不知道到底有多少乐高积木，这就是我们试图用概率来确定的。然后，Joe 可以将这个数字与他获得的总数进行比较，看看他是否想再进行一次传递。因为我花了几天时间来解决这个问题，所以我画了一个模型，然后回去解释它:
在这种情况下， p 表示乔捡起一个乐高积木的概率。概率 p 被建模为贝塔分布(在第 1 部分中我们介绍了为什么使用贝塔分布)。剩余的乐高数是 N_i ，在 i 扫完被乔捡起的乐高后；乔总共扫了 n 次。每扫一次，我就移除一个乐高积木。我们试图确定的是 N₀ ，乔的后代最初留下的乐高积木的数量。 R₀ 被设置为零处理边界条件，其中乔在第零次扫描时拾取零个乐高积木。
为了完成这个模型，我们需要一些关于 N₀ 、 α 、 β 的先验分布。与第一集一样，我们选择参数，使建模前的分布信息较弱:
先验概率分布的参数与第一集相同。我问乔他的孩子可能会留下的绝对最大乐高玩具。Joe 估计大约有 2000 个，所以为了安全起见，我选择了最大值，N _ max3000 个，并使用 0 和 N_max 之间的先验概率的均匀分布来模拟留下的乐高积木的初始数量。
现在剩下的就是在 PyMC 中制作模型，并检查结果，看看它们是否能帮助 Joe。
模型设置和拟合
从模型到 PyMC 代码的转换实际上是一种直接的转换，并且适合默认的不掉头采样器。经过一番折腾后，调优样本的数量超过了默认值，目标接受率也比正常情况下设置得稍高一些；这样做是为了处理这个问题中的一些有问题的后验。很好，PyMC 库通过模型拟合检查和突出可能的问题，很好地帮助用户调整这些参数。
我想让乔在每次清扫后大概知道还剩多少乐高积木。这将有助于 Joe 更好地决定何时停止，因为我知道 Joe 非常想减少他在乐高玩具收集上投入的精力。这样，我们将从最少的信息开始，仅是在第一次扫描中删除的数字。此外，还包括了剩余乐高玩具估计数量的确定性计算，主要是为了让 Joe 更容易确定剩余的乐高玩具数量。即使 n_remaining_est 是确定性的，它也是随机变量，因为 n_remaining_est 的两个输入变量是随机变量。
sweep_model = pm.Model()with sweep_model:
n_0_est = pm.DiscreteUniform('n_0_est', lower=0, upper=3000)
alpha_sweep = pm.HalfCauchy('alpha_sweep', beta = 2.)
beta_sweep = pm.HalfCauchy('beta_sweep', beta = 2.)
p_sweep = pm.Beta('p_sweep', alpha=alpha_sweep, beta=beta_sweep)
removed_sweep_1 = pm.Binomial('removed_sweep_1', p=p_sweep,
n=n_0_est, observed=Removed[1] )
n_remaining_est = pm.Deterministic('n_remaining_est', n_0_est - removed_sweep_1)with sweep_model:
trace = pm.sample(50000, tune=7500, random_seed=123, progressbar=True, nuts_kwargs={"target_accept":0.95})
模型评估
完成模型拟合后，Joe 和我快速浏览了一下轨迹图。我们关注的是乐高玩具的初始数量和第一次扫描后剩余的估计数量。
不幸的是，仅仅一次清扫后，剩余的乐高积木数量与清扫前相比并没有减少多少。从数学直觉来看，这是有意义的:有两个参数需要估计，捡起一个乐高的概率和乐高的初始数量，但只有一个已知变量，第一遍捡起的乐高数量。仅仅用一个已知方程来求解两个未知变量是很困难的。
pm.traceplot(trace, varnames=['n_0_est', 'n_remaining_est'], combined=True);
pm.summary(trace, varnames=['n_0_est', 'n_remaining_est']).round(2)
第二次扫描
在这一点上，Joe 不确定这项工作是否值得。我鼓励乔，我们应该继续，因为他的脚如果没有别的。我们坚持下去，根据 Joe 提供的第一个和第二个观察到的扫描数据创建一个模型。模型创建代码开始看起来有点重复，但是为了快速实验，我做了一些代码的“复制-粘贴”并运行模型，想着如果它工作了，我以后会重构它。
在此模型中，我们从估计的初始乐高数中减去第一次扫描后剩余的乐高数，并将其用作第二次扫描剩余乐高数的输入。我们正在模拟重复的、连续的二项式试验。
sweep_model = pm.Model()with sweep_model:
n_0_est = pm.DiscreteUniform('n_0_est', lower=0, upper=3000)
alpha_sweep = pm.HalfCauchy('alpha_sweep', beta = 2.)
beta_sweep = pm.HalfCauchy('beta_sweep', beta = 2.)
p_sweep = pm.Beta('p_sweep', alpha=alpha_sweep, beta=beta_sweep)
removed_sweep_1 = pm.Binomial('removed_sweep_1', p=p_sweep,
n=n_0_est, observed=Removed[1] )
n_1_est = n_0_est - removed_sweep_1
removed_sweep_2 = pm.Binomial('removed_sweep_2', p=p_sweep,
n=n_1_est, observed=Removed[2] )
n_remaining_est = pm.Deterministic('n_remaining_est', n_0_est - removed_sweep_1 - removed_sweep_2)with sweep_model:
trace = pm.sample(50000, tuning=7500, random_seed=123, progressbar=True,
nuts_kwargs={"target_accept":0.95})
“乔，我们现在有所进展了。”
快速浏览一下模型结果，我们看到乐高剩余部分的估计范围已经缩小到 9 到 62 之间，置信度为 95%。在这一点上，Joe 受到了鼓舞，认为清点所有这些乐高积木可能是值得的。他想继续，但我有其他计划。
有点晚了，我决定是时候重构这段代码了。乔回家过夜，我们答应明天聚一聚。有了一点干净的代码，我希望能够移动得更快。
pm.traceplot(trace, varnames=['n_0_est', 'n_remaining_est'], combined=True);
pm.summary(trace, varnames=['n_0_est', 'n_remaining_est']).round(2)
通过数学清理代码
第二天晚上，乔手里拿着外卖来了。我抓起我的笔记本电脑，与乔分享了一夜之间发生的事情。
“这确实是一个多项问题.”
多项式分布模拟 k 个可能事件中的某个事件发生若干次 n 的次数。它是二项分布的推广。想象重复投掷一个六面骰子(6 个可能的事件，1 到 6，骰子的每一面一个)，并计算它每一面朝上落地的次数。多项式模型允许每个事件有不同的发生概率。继续 6 面骰子的类比，如果骰子以某种方式具有不同的每一面向上降落的概率。这就是多项式分布模型的过程。乔有点迷惑地看着我，所以我继续说。
在我们的场景中，每个落在 n 清扫或中的乐高积木都还没有被拾起。每扫一次，加上没被捡的事件，就像是骰子的一面。乐高在第一次扫描中被捡起的概率是未知的，被建模为 p 。乐高在第二次扫描中被捡起的概率是乐高在第一次扫描中没有被捡起的概率和在第二次扫描中被捡起的概率。剩余乐高积木的概率就是在所有的扫描中，将所有的乐高积木被选中的概率相加后剩下的概率。
这里 R 是每个给定概率的计数向量(也就是说，它是每次扫描中捡出的乐高积木的估计数量)， N₀ 是试验次数(或者在我们的例子中，是最初留下的乐高积木的估计数量)， p_i 是在扫描中捡出一个乐高积木的概率， i ， p_remains 是一个乐高积木没有被捡出的概率
虽然 PyMC3 具有多项分布，但它不具有这种“不完全多项分布”，多项总体的大小是未知的，并且是待估计的。已经有几篇关于从观察到的图中估计多项人口规模的论文；我最喜欢的是 Sanathanan 的。我致力于采用贝叶斯方法并使用 PyMC3。经过一点试验，我确定可行的解决方案是推广之前尝试的重复二项式方法。我写了一个封装了移除模型的函数，称这个新模型为RemovalMultinomial。虽然我肯定有更好的方法来做到这一点，乔和我已经准备好了。
# Done quickly to make code a bit easier to deal with
# Doesn't do much / any error checking of the inputs, assumes observed information
def RemovalMultinomial(basename, n, p, observed):
remaining = n
results = []
for idx, p_idx in enumerate(p):
i = idx+1
name = basename + '_{}'.format(i)
removed = pm.Binomial(name, n=remaining, p=p_idx, observed=observed[idx])
results.append(removed)
remaining = remaining - removed
remaining = pm.Deterministic( basename+'_remaining', remaining)
results.append(remaining)
return results
有了这个新函数，模型的设置就短得多，只需通过扫描次数就可以参数化。我们将用 3 次扫描来做这个模型，比上次多一次。这有助于乔和我了解我们一次又一次地收集了多少信息。此外，这可能有助于 Joe 根据他自己的风险模型确定他何时可以在传球中“提前”停止，这些风险模型是指地板上残留的一定密度的乐高玩具，如果他踩到它们，会给他带来疼痛。现在，乔一直在寻找乐高积木，直到他经过时，没有看到任何要捡起的乐高积木。
n_sweeps = 3sweep_model = pm.Model()with sweep_model:
n_0_est = pm.DiscreteUniform('n_0_est', lower=0, upper=3000)
alpha_sweep = pm.HalfCauchy('alpha_sweep', beta = 2.) # Weakly regularizing prior
beta_sweep = pm.HalfCauchy('beta_sweep', beta = 2.) # Weakly regularizing prior
p_sweep = pm.Beta('p_sweep', alpha=alpha_sweep, beta=beta_sweep)
p_multinomial = [ pm.Deterministic('p_{}'.format(str(i+1)), p_sweep) for i in range(0,n_sweeps) ]
removed = RemovalMultinomial('rr', n=n_0_est, p=p_multinomial, observed=Removed[1:] )with sweep_model:
trace = pm.sample(50000, tuning=7500, random_seed=123, progressbar=True,
nuts_kwargs={"target_accept":0.95})
“这看起来相当不错！”乔说
乔和我都很兴奋。在这最后一轮，经过三次捡乐高积木的扫描后，现在有 95%的置信区间是零(！)和剩余的十三块乐高积木。乔对这些结果非常满意，并认为这可以帮助他决定是继续传球还是收工。
pm.plot_posterior(trace, varnames=['n_0_est', 'rr_remaining']);
pm.summary(trace, varnames=['n_0_est', 'rr_remaining']).round(2)
一个最终模型
在最后一个模型中，我们想测试 Joe 目前的启发式方法对停止有多好。乔通常会在经过时停止寻找乐高积木，因为他发现没有更多的乐高积木可捡。现在我们有了数学，我们可以量化，至少基于这些结果和模型的假设，Joe 采用的停止策略的质量。我们创建了一个所有六个回合的模型，并使用它来估计剩余的乐高积木的数量。这里没有显示代码，但是可以在文章末尾链接的笔记本中找到。
plt.figure(figsize=(4,3), dpi=300)
plt.hist(trace['rr_remaining'], density=True, bins=range(0,5), align='left')
plt.xlabel('Estimated Legos Remaining')
plt.xticks(range(0,5))
plt.ylabel('Probability')
plt.title('Probability Density Function of Estimated Legos Remaining')
plt.grid(True)
plt.show();
有些模型是有用的
此时，Joe 和我对模型结果非常满意。看起来，至少在这个例子中，有大约五分之一(20%)的机会至少有一个乐高剩余。Joe 说，如果我把它清理干净并粘贴到移动应用程序中，他就可以拿走并使用它。乔是一个非常有趣的人！我同意将此代码签入 GitHub，并向 Joe 展示如何使用它。
一切都好，乔回家了。
然后几天后，电话来了，“汉克，第二天我还是发现乐高玩具，而且比我认为模型预测的要多得多。”
在简单模型中发现商业价值
现实世界中的数据科学
复杂模型失败时，一般线性模型如何成功
直到今天，我在野外看到的最有用的预测模型是一般的线性模型。这不是最强大的模型，它没有做出最准确或最精确的预测，房间里的工程师们也没有对谈论它感到兴奋。然而，业务用户信任并理解它的输出，并且该模型被吸收到他们的决策过程中。
我相信有一个很好的解释，为什么这个模型被非技术团队所接受，并在许多复杂模型失败的地方取得了成功。该模型完成了比重复推荐或标签更有价值的事情——它促进了需要做出重要产品决策的团队之间的对话。尽管这个模型很简单，但它的影响波及了整个企业。
我从这一结果中看到了一些有价值的教训，在整篇文章中，我的目标是阐明这些教训。
第一课:给非技术团队一个关心的理由
在技术团队和业务团队之间往往有一堵墙。聪明人生活在这堵墙的两边，但总的来说，他们说着不同的语言，被不同的动机所驱使。对工程师来说似乎重要的事情对产品经理来说可能显得微不足道，反之亦然。
如果我们对什么是真正重要的看法不一致，很自然的，电线会交叉，沟通会发生错误。对于我参与的许多项目来说，这的确是事实。这种现象并不仅限于技术人员和非技术人员；当顾问为他们的客户工作时，这种情况经常发生。如果我们不能持续地在期望值上保持一致，那么客户想要的和实际交付的路径自然会出现分歧。
我们可以从中学到的是:你讨论模型的方式很重要。如果你面对的是技术观众而不是非技术观众，你可以用不同的方式组织讨论，注意强调对你的观众来说重要的观点。您可能对业务开发演示中讨论的细节不感兴趣，所以不要错误地假设任何人都和您一样对您的模型的技术细节感兴趣。
我相信这有助于解释为什么一般的线性模型能够吸引决策者的注意力，即使“更好的”模型被忽略。线性模型易于理解，它给了企业一个关注它的理由，因为它清楚地说明了过去的决策、客户行为和业务绩效之间的关系。在这种情况下，该模型能够解决业务用户的问题，并且能够更好地为他们的决策过程提供信息。
第二课:强调结果，而不是技术
随着我们的技术团队所拥有的工具不断进步，我们考虑技术团队和非技术团队之间日益扩大的鸿沟变得越来越重要。例如，想想那些喜欢你不喜欢的爱好的人，想想当他们开始谈论你听不懂的行话时，你有多快就听不进他们的谈话。现在考虑一下技术人员的词汇和他们的非技术同行的词汇有什么不同。这与智力没有任何关系，但这确实意味着根据你属于哪个团队，会有不同的“商业方言”。
尽管工作场所的不同学科之间存在差异，但在那里工作的团队都存在，因为他们服务于一些有益于业务的功能。就像你在国外旅行时可能会和不讲你的语言的人找到共同点一样，你也可以和你的同事找到共同点，不管你是做什么工作的。关键是不要失去你们共同的目标。
因此，当开发一个旨在推进业务目标的预测模型时，您应该后退一步，考虑该模型对您的业务对手意味着什么。你可能爱上了一个新的框架，或者熬夜思考一个新的算法，但是这对业务真的重要吗？
我认为这无关紧要。是的，你应该和其他同样对你感兴趣的领域感兴趣的同行一起庆祝你的成功，但是预测模型只有在被使用时才是有用的。如果企业确信你的模型能产生结果，他们才会使用你的模型，所以关注结果似乎是一个成功的策略。
给我的客户的业务用户留下深刻印象的一般线性模型并没有给他们留下深刻印象，因为它的实现或编写它的编程语言。这给他们留下了深刻的印象，因为它产生的结果可以融入他们的决策过程。结果以他们理解的行话呈现给业务部门，这是数据科学团队在门口检查自我的一个例子，以便在共同点上满足业务需求，从而推进整个业务的目标。
第三课:把结果解释为指导，而不是答案
我倾向于假设我的模型中存在大量的偏差。这些年来与数据的密切合作对我产生了清醒的影响，现在我知道香肠是如何制作的了。建立一个模型需要很多假设；在生成数据、收集数据、清理数据、处理数据时有一些假设，最后还有一个假设，即分析历史数据可以准确预测未来的行为。就像“电话”游戏一样，这种假设系统几乎可以确保输出是有缺陷的。
出于这个原因，我明白了把模型的结果当作一盏指路明灯，而不是某个真理，通常会更好。我们在雾中飞行，我们模型的结果是指引我们的灯，但它们不一定是最终的目的地。
我发现这种方法引起了企业决策者的共鸣。如果数据是有缺陷的，并且数据通过的所有中间系统都是有缺陷的，那么将模型的结果解释为纯粹的真理是不明智的。此外，许多人不喜欢被告知该做什么，如果你告诉非技术人员听你的模型，他们就会忽视你。
然而，如果你把你的模型当作对话的向导，并且你能证明你的模型所暗示的关系是真实的，那么你就处于有利地位。对于非技术团队来说，这是一个味道更好的配方，他们可能不一定相信数据是所有事情的答案。
这是一般线性模型如此成功的原因之一。开发这一模型的分析团队并没有声称拥有所有答案，而是小心翼翼地专注于他们的模型所确定的关系。商业用户对该模型的结果反应良好，我认为他们的积极反应可能是由于该模型量化了他们在某种程度上已经意识到的某些关系。
第四课:好的模型并不总是保证好的决策
我留给你们的最后一则轶事表明，即使是设计良好的模型也可能导致糟糕的结果，这取决于对它的解释。就我提到的一般线性模型而言，在它取得初步成功几个月后，关于如何利用该模型帮助量化的关系，不同的观点开始出现。
例如，该模型指出的一个观点是，购买产品“X”的客户与倾向于导致其他有利可图的购买的行为相关联。总之，谈话最终变得过度集中在如何销售更多的产品“X”上。最终的商业决策是为了不惜一切代价卖出更多的产品“X ”,这违背了分析团队的建议，并且结果远不如业务预期的那样有利可图。他们选择的策略导致了他们的产品市场中的非预期行为，这点击了模型所揭示的关系的重置按钮。
在我看来，这不是模型的失败。该模型正确地识别了存在的关系。问题是，该企业尝到了最初见解的成功滋味，却草率行事了。他们暂时接受了数据驱动的流程，但一旦他们觉得事情朝着正确的方向发展，他们很快就会转身将这些数据驱动的流程从他们的决策中删除。对话很快变成了政治话题，业务团队内部形成了不同的阵营。人们加强了他们的防御阵地，随之而来的沟通失误降低了预测模型最初的成功。
我敢肯定，从这种情况下，我们还可以学到更多的东西。我在这里的目的是说明技术团队，特别是分析和面向数据的团队，可以做些什么来提高他们的工作被业务利用的几率。虽然这里描述的结果并不完全是一个童话故事的结局，但该模型在渗透到战略性商业对话的能力方面是独一无二的成功。
最终，我认为如果他们的业务部门继续与构建预测模型的技术团队合作，该公司会受益更多。也许他们已经意识到他们所做的一些糟糕的假设，即把所有的鸡蛋放在增加产品 X 销售的篮子里。总的来说，我认为分析团队应该在决策过程中拥有更大的发言权。不过，那是改天再聊。
我希望你喜欢阅读！谢谢收听。
寻找循环模式:如何用 Python 实现 STFT 的教程
如果你听说过 Python 和傅立叶名词，你可能会发现这篇文章很有用:在这里，我将探索一种简单的方法来实现 Python 中的短时傅立叶变换，以便运行频率分析来检测给定信号中的循环模式。我们开始工作吧！
去年，我和我研究组的同事应邀参加了一个 Kaggle 的比赛。这个想法对我来说听起来非常有趣，因为它是关于地震数据的，即时间序列，这在某种程度上与我习惯处理的信号类别有关——出于好奇，我将小波应用于从大脑信号中检测睡意。顺便说一下，当 Stéphane Mallat 将小波分析应用于地震和地质数据时，它就变得很有名了。
但是在我们开始处理小波之前(我打算在不久的将来回到小波)，这里的主要思想是熟悉周期信号。我的意思是，一个时间序列代表一些周期性事件，不时重复。因此，这通常可以被称为一种模式。
你可能会发现以最不同的形式表现的模式，但有一个是特别的圆:白天和晚上的时间，一年中的季节，水力发电产生的电力，所有这些都与一个圆形源有关(或者至少是一个近似的，因为地球是一个扁球体，它以椭圆路径围绕太阳运行。
从数学的角度来看，任何循环模式都可以用正弦和余弦三角函数来表示，正如法国数学家和物理学家约瑟夫·傅立叶所证明的那样，他在十九世纪初建立了今天所知的傅立叶级数和傅立叶变换。
The sine and cosine functions are the projections of a circular movement over a vertical and horizontal axis, respectively.
我假设你已经熟悉三角信号及其与频率的关系。在任何情况下，重要的是要注意，周期性或循环信号将在给定的时间内重复自身，称为周期 ( T )。而频率 ( f )则由每秒的重复次数组成，以赫兹(Hz)为计量单位，以此向德国物理学家海因里希·赫兹致敬。
在接下来的章节中，我们将以一种非常实用的方式探索这些概念(我还假设您已经习惯了 Python)。我们将首先弄清楚如何生成和可视化一些周期波(或信号)，然后移动到短时傅立叶变换实现。最后，我们将看到一个例子，说明如何将这种分析应用于一个非三角函数，而是一个节奏函数，如大脑信号。
1.产生周期性信号
在数据科学——这里我考虑的是与之相关的所有学科，比如模式识别、信号处理、机器学习等等——深入了解我们的数据(也就是信号)在插入分析或分类工具时的表现总是有用的。问题是，起初我们通常对我们的数据知之甚少——在现实世界中，这些数据可能是复杂和嘈杂的——事实上，寻找分析工具就是试图理解它们。
在这种情况下，最好的替代方法是使用简单且众所周知的数据来探索分析工具。诺贝尔奖获得者理查德·费曼在 1988 年去世后，在他的黑板上留下了下面这句名言:“我不能创造的东西，我不理解”——我们将首先产生自己的周期信号。
时间序列通常是抽样过程的结果。当这个信号产生时，正如这里的情况，我们还必须确定信号应该出现的时间。下面的函数模拟采样过程，并将初始和最终时间样本作为输入。换句话说，它返回我们的信号将被放置的 x 轴基准。
让我们记住正弦信号是由函数产生的:
其中 f 是信号 y(t) 的频率，单位为赫兹， t 是上述代码生成的时间序列。然后将该等式转录为 Python:
几乎与生成信号同样重要的是提供一种可视化信号的方法。使用 Matplotlib 库，可以修改以下函数来查看任何时间序列:
编写一个包装函数来创建和可视化一个正弦函数，并将频率和幅度作为输入参数，这听起来也很实用:
最后，我们现在将一起执行这些函数，以便获得一些正弦波:
然后给出以下结果:
它可能看起来很简单，但重要的是要有一个我们确切知道它是如何产生的信号——我们持有它的模型，即它的公式。在上面的第一张图中观察到，频率为 0.1Hz 意味着信号将每 10 秒重复一次。在第二张图中，我们复制了信号的振幅，现在信号每秒钟重复一次。
将我们用 numpy.concatenate 方法定义的函数放在一起，如下所示，我们可以生成更复杂的波形:
从我们的第一张图表可以清楚地看出，当我们处理低频和单一频率时，我们仍然只能通过时域信号中的视觉检查来识别循环模式(例如，在一秒钟内计数顶部)。然而，对于更复杂的波形(见上图中 1s 和 2s 之间的间隔),显然需要一种工具来识别任何循环模式。
2.短时傅立叶变换
要完全理解频域中的时间序列表示，需要深化线性代数和微积分，以便将信号展开成一系列函数。
不过，这里不会采用这种方法:那些已经理解这种数学抽象的人会发现用 Python 实现频率分析的实用代码。同理，第一次面对周期性信号的人，可能会从这些简单的例子开始加深自己的理论理解。
在我们开始编码之前，我将提供一个简短的定义*(请不要害怕这些符号……如果你愿意，你可以直接向下翻页跳过这一部分):傅立叶已经证明，可以将大量函数展开为纯调和函数级数*(即正弦和余弦之和)。顺便说一下，从欧拉公式我们知道我们可以将调和函数表示为一个复指数:
现在考虑我们可以将给定信号 y(t) 展开成傅立叶级数:
其中 c_ω 为信号经过傅里叶变换得到的傅里叶系数；而 e^iωt 是代表频率ω的谐波基函数。
重要的是要记住，每个系数——我们的计算机将为我们计算的系数——都与一个谐波基函数相关联，而谐波基函数又将该系数与特定频率相关联。然后，如果我们知道每个系数的值，我们可以知道信号 y(t) 中存在的每个频率分量的强度。
应该考虑的另一个重要方面是，这些基函数具有无限的持续时间，因此不适合表示具有不连续性的信号，或者也需要其时间或空间位置的信号(即，当我们试图知道给定频率或节奏在信号的哪个周期中出现得最多时)。
短时傅立叶变换 (STFT) 就是克服这一点的方法。其策略包括将变换的每个基函数乘以*窗口函数 w(t)。*因此，后者限于给定的时间段，即仅在有限的时间间隔内具有非空值。
考虑到傅立叶变换中的函数空间可以由复数形式的三角函数定义，对于 STFT，我们创建了由频率(ω)和位置(τ)定义的新基:
必须根据分析兴趣选择 w(t) 的形状和长度。形状取决于用来产生它的函数，并决定了它的频率分辨能力。另一方面，长度( N )定义了窗口间隔，并因此定义了其时间分辨率。
可以使用和尝试许多窗口功能(可以在这里找到其主要功能和差异的概述)。就我们的目的而言，或者在犹豫该选择哪一个时，汉恩的窗口可能是一个令人满意的选择:
现在让我们回到我们的 Python 代码。一旦我们理解了 STFT 所依赖的基本原理，我们就可以利用 SciPy 库中的 信号模块来实现频谱图——它由绘制 STFT 系数的平方幅度组成。
因此，使用下面的代码，我们将计算第一个信号的 STFT(向上翻页，寻找信号 1 )。此外，我们必须考虑到 SciPy 实现将始终将频率轴设置为采样频率的一半，这也是我设置 plt.ylim() 命令的原因。
好吧，我同意你的观点，这个结果和我预期的相差甚远。好消息是(几乎)一切都是调整的问题。为了使事情变得简单，让我们编写下面的函数:
为了测试它，我们将创建一个新信号，其值易于理解和操作，一个振幅为 1000 的 100Hz 正弦波。由于我们将保持它持续 10s，以相同的采样率 1kHz ，我们可以使用之前定义的数组时间。
从一个简单的信号开始的目的是，它允许我们检查我们的代码中是否有重大错误。还好不是这样！正如预期的那样，绘图频率范围上升到 500Hz，是采样频率的一半。此外，我们可以精确地看到 100Hz 分量，其值被归一化(你会在这里找到归一化策略的细节)。
更进一步，我们可以为我们的 STFT 指定其他参数。为了创建下面的图表，我为参数 nperseg 设置了不同的值，这些值对应于窗口函数的大小。对于左图， nperseg = 64 ，而对于右图， nperseg = 2048 。请注意，窗口越大，频率分辨率越精确。反之亦然。
没有特定的规则来确定每个片段的样本数。然而，考虑采样率的 1/4 是有用的。例如，对于 1kHz 采样速率， nperseg 参数默认设置为 256 。如果我们分析我们的复合信号，这种权衡会变得更加清晰，复合信号是由 0.25 至 123Hz 的分量产生的:
在上面的图中，信号的所有已知频率成分都变得明显，尽管时间分辨率较低。对于某些应用来说， 250ms 并不重要。对其他人来说，这可能是限制性的。我必须预见到，没有一种分析工具可以像海森堡原理所证明的那样，在时间和频率上同时提供准确性——但这将是下一篇文章的主题。
3.大脑信号分析
为了举例说明 STFT 分析如何在“现实生活”中发挥作用，我们现在将应用我们在从脑电图(EEG)中获得的大脑信号片段上构建的函数。如果你想自己尝试一下，数据在这里有。
由于信号之前是在 Matlab 上处理的，我们首先必须使用 loadmat 方法将数据结构转换成 Python 字典。知道信号被标记为 sD_W ，我们可以继续加载它:
与处理正弦信号时不同，现在我们不能仅通过查看其时域图来推断 EEG 信号中存在的节律(即频率分量)。我们甚至可以说在 20 秒后会发生一些事情，但是我们不能给出任何细节。
另一方面，上面的声谱图向我们显示，在大约 20s 的事件之后，从 20s 到 36s，然后从 44s 到 56s，出现了一个良好建立的低频节奏。与任何数据科学活动一样，解释总是依赖于上下文。在这种情况下，该 EEG 信号来自接受多导睡眠描记术的受试者。所识别的频率范围对应于阿尔法节律，指示患者可能是清醒的和放松的，但是闭着眼睛。最后，50Hz 附近的黑暗区域对应于信号采集期间应用的带通滤波器，以便过滤来自交流电源的 50Hz 干扰。
4.结论
这篇文章的目标不仅是展示如何用 Python 实现 STFT，还简要介绍了这个强大的分析工具背后的理论——它支持更复杂的数学思想。
此外，通过能够理解在这种分析中应该观察的主要参数，并借鉴大脑信号的例子，我希望读者可以定制代码，并将这样的工具应用于他/她可能找到的任何循环信号:从生物医学或地震信号，甚至到股票市场价格或恒星的运动。
原载于 2019 年 2 月 23 日http://tbnsilveira . info。
寻找捐赠者:PySpark 分类项目
了解如何使用 Apache PySpark 来增强您的分类预测
Picture from Unsplash
介绍
本文的目的是对机器学习中的分类问题做一个温和的介绍，并通过一个全面的指南来成功地使用 PySpark 开发一个类预测。
所以没有进一步的行动，让我们投入进去吧！
分类
如果你想深入了解分类问题、它的主要算法以及如何使用机器学习技术来处理它们，我强烈建议你去看看下面这篇文章，在这篇文章中，我详细地解释了这个概念。
了解机器如何分类
towardsdatascience.com](/supervised-learning-basics-of-classification-and-main-algorithms-c16b06806cd3)
什么是分类？
分类是监督学习的一个子类，其目标是基于过去的观察来预测新实例的分类类别标签(离散的、无序的值、组成员)。
分类问题主要有两种类型:
- 二元分类:典型的例子是邮件垃圾邮件检测，其中每封邮件都是垃圾邮件→ 1 封垃圾邮件；或者不是→ 0。
- 多类分类:像手写字符识别(类从 0 到 9)。
下面的例子很有代表性地解释了二元分类法:
有两个类别，圆和十字，以及两个特征，X1 和 X2。该模型能够找到每个数据点的特征与其类别之间的关系，并在它们之间设置边界线，因此当提供新数据时，它可以估计其所属的类别，给定其特征。
Figure by Author
在这种情况下，新数据点落入圆形子空间，因此，模型将预测其类别为圆形。
分类主要算法
为了预测某些样本的类别，可以使用几种分类算法。事实上，在开发我们的机器学习模型时，我们会对其中的一定数量进行训练和评估，我们会保留那些预测性能更好的模型。
一些最常用的算法的非穷尽列表如下:
- 逻辑回归
- 决策树
- 随机森林
- 支持向量机
- k-最近邻(KNN)
分类评估指标
当对事件进行预测时，我们可以得到四种结果:
- 真阳性:TP
- 真阴性:TN
- 误报:FP
- 假阴性:FN
所有这些都在下面的分类矩阵中表示:
准确性衡量分类器做出正确预测的频率。它是正确预测的数量与预测总数(测试数据点的数量)的比率。
Precision 告诉我们我们归类为某一类的事件的比例，实际上是那个类。它是真阳性与所有阳性的比率。
回忆(敏感度)告诉我们实际上属于某一类别的事件中有多少被我们归类为该类别。它是真阳性与所有阳性的比率。
特异性是被正确识别为阴性的类别占阴性类别总数的比例。
对于分类分布有偏差的分类问题，准确性本身并不是一个合适的度量。反而精度和召回率更有代表性。
这两个指标可以结合起来得到 F1 分数，它是精确度和召回分数的加权平均值(调和平均值)。这个分数的范围从 0 到 1，1 是可能的最佳 F1 分数(我们在处理比率时采用调和平均值)。
皇家对空观察队
最后，我们将在项目中使用的指标是接收机操作特性或 ROC 。
ROC 曲线告诉我们该模型区分两个类别的能力有多强。它可以得到从 0 到 1 的值(€[0，1])。模型越好，越接近 1 值。
Figure by Author
从上图中可以看出，我们的分类模型将在类别和之间绘制一个分离边界:
- 落在阈值左边的每个样本将被归类为负类。
- 落在阈值右边的每个样本将被分类为阳性类别，
预测的分布如下:
Figure by Author
灵敏度&特异性之间的权衡
当我们降低阈值时，我们最终会预测更多的正值，并提高灵敏度。因此，特异性降低。
当我们提高阈值时，我们最终会预测更多的负值并提高特异性。因此降低了灵敏度。
如敏感性 ⬇️ 特异性 ⬆️
如特异性 ⬇️ 敏感性 ⬆️
为了优化分类性能，我们考虑(1-特异性)而不是特异性。所以，当敏感性增加时，(1-特异性)也会增加。这就是我们计算 ROC 的方法。
Figure by Author
性能示例
如前所述，赋值器越接近 1，模型的预测性能就越好，类之间的重叠区域就越小。
Figure by Author
寻找捐助者项目
可以在以下文章中找到该项目的完整介绍:
使用此分类模型查找并预测谁将向慈善机构捐款！
towardsdatascience.com](/classification-project-finding-donors-853db66fbb8c)
在本文中，我们将关注该项目的 PySpark 实现。
总之，在整个项目中，我们将使用许多不同的监督算法，使用从 1994 年美国人口普查中收集的数据来精确预测个人收入。
然后，我们将从初步结果中选择最佳候选算法，并进一步优化该算法以最佳地模拟数据。
我们实现的目标是构建一个模型，准确预测个人收入是否超过 50，000 美元。根据我们之前的研究，我们发现最有可能向慈善机构捐款的人是那些年收入超过 5 万美元的人。
因此，我们面临着一个二元分类问题，我们想确定一个人的年收入是否超过 5 万美元(第 1 类)或没有(第 0 类)。
这个项目的数据集来源于 UCI 机器学习知识库。
数据
人口普查数据集由大约 45222 个数据点组成，每个数据点有 13 个特征。
特性
age:年龄workclass:工人阶级(私营、自营企业、自营企业、联邦政府、地方政府、州政府、无薪、从未工作)education_level:教育水平(学士、专科学校、11 年级、高级研究生、专业学校、acdm 协会、voc 协会、9 年级、7-8 年级、12 年级、硕士、1-4 年级、10 年级、博士、5-6 年级、学前班)education-num:完成的教育年数marital-status:婚姻状况(已婚-未婚-配偶、离婚、未婚、分居、丧偶、已婚-配偶不在、已婚-配偶)occupation:工作职业(技术支持、工艺维修、其他服务、销售、行政管理、专业教授、搬运工人、清洁工、机器操作员、行政文员、农业渔业、运输搬运、私人服务、保安服务、武装部队)relationship:关系状态(妻子、亲生子女、丈夫、非家庭成员、其他亲属、未婚)- 种族(白人、亚洲太平洋岛民、美洲印第安爱斯基摩人、其他人、黑人)
sex:性别(女，男)capital-gain:货币资本收益capital-loss:货币资本损失hours-per-week:每周平均工作时间native-country:本土国家(美国、柬埔寨、英国、波多黎各、加拿大、德国、美国外围地区(关岛-USVI 等)、印度、日本、希腊、中国、古巴、伊朗、洪都拉斯、菲律宾、意大利、波兰、牙买加、越南、墨西哥、葡萄牙、爱尔兰、法国、多米尼加共和国、老挝、厄瓜多尔、台湾、海地、哥伦比亚、匈牙利、危地马拉、尼加拉瓜、苏格兰、泰国、南斯拉夫、萨尔瓦多、特立尼达和多巴哥&多巴哥、秘鲁、香港、荷兰)
目标变量
income:收入阶层(< =50K，> 50K)
导入数据和探索性数据分析(EDA)
我们将从导入数据集并显示数据的前几行开始，以对探索性数据分析进行初步近似。
**# File location and type**
file_location = "/FileStore/tables/census.csv"
file_type = "csv"**# CSV options**
infer_schema = "true"
first_row_is_header = "true"
delimiter = ","**# The applied options are for CSV files. For other file types, these will be ignored.**
df = spark.read.format(file_type) \
.option("inferSchema", infer_schema) \
.option("header", first_row_is_header) \
.option("sep", delimiter) \
.load(file_location)display(df)
我们现在将通过使用。describe()方法。
**# Display Dataset's Summary**
display(df.describe())
让我们也找出数据集的模式。
**# Display Dataset's Schema**
display(df.describe())
准备数据
因为我们想预测个人年收入是否超过 50K 美元，我们将把标签“收入”替换为“> 50K”。
为此，我们将创建一个新列，根据个人年收入是否超过 5 万美元，该列的值将为 1 或 0。然后我们将删除收入栏。
**# Import pyspark functions**
from pyspark.sql import functions as F**# Create add new column to the dataset**
df = df.withColumn('>50K', F.when(df.income == '<=50K', 0).otherwise(1))**# Drop the Income label**
df = df.drop('income')**# Show dataset's columns** df.columns
向量化数字特征和一次性编码分类特征
为了对模型的训练进行处理，Apache Spark 中的特征必须被转换成向量。这个过程将使用我们现在将要探索的某些类来完成。
首先，我们将导入相关的库和方法。
from pyspark.ml import Pipeline
from pyspark.ml.feature import StringIndexer, OneHotEncoder, VectorAssembler
from pyspark.ml.classification import (DecisionTreeClassifier, GBTClassifier, RandomForestClassifier, LogisticRegression)
from pyspark.ml.evaluation import BinaryClassificationEvaluator
现在，我们将选择分类特征。
**# Selecting categorical features**
categorical_columns = [
'workclass',
'education_level',
'marital-status',
'occupation',
'relationship',
'race',
'sex',
'hours-per-week',
'native-country',
]
为了对这种分类特征进行一次性编码，我们将首先将它们通过索引器，然后传递给编码器。
**# The index of string values multiple columns**
indexers = [
StringIndexer(inputCol=c, outputCol="{0}_indexed".format(c))
for c in categorical_columns]**# The encode of indexed values multiple columns**
encoders = [OneHotEncoder(dropLast=False,inputCol=indexer.getOutputCol(),
outputCol="{0}_encoded".format(indexer.getOutputCol()))
for indexer in indexers]
现在，我们将把分类编码的特征与数字特征结合起来，并用它们两者构成一个向量。
**# Vectorizing encoded values**
categorical_encoded = [encoder.getOutputCol() for encoder in encoders]
numerical_columns = ['age', 'education-num', 'capital-gain', 'capital-loss']
inputcols = categorical_encoded + numerical_columns
assembler = VectorAssembler(inputCols=inputcols, outputCol="features")
现在，我们将建立一个管道来自动化这个阶段。
pipeline = Pipeline(stages=indexers + encoders+[assembler])
model = pipeline.fit(df)**# Transform data**
transformed = model.transform(df)
display(transformed)
最后，我们将只选择具有相关特征的数据集。
**# Transform data**
final_data = transformed.select('features', '>50K')
初始化模型
对于这个项目，我们将研究三种不同分类算法的预测性能:
- 决策树
- 随机森林
- 梯度增强树
**# Initialize the classification models**
dtc = DecisionTreeClassifier(labelCol='>50K', featuresCol='features')
rfc = RandomForestClassifier(numTrees=150, labelCol='>50K', featuresCol='features')
gbt = GBTClassifier(labelCol='>50K', featuresCol='features', maxIter=10)
拆分数据
我们将在训练和测试数据之间执行经典的 80/20 拆分。
train_data, test_data = final_data.randomSplit([0.8,0.2])
训练模型
dtc_model = dtc.fit(train_data)
rfc_model = rfc.fit(train_data)
gbt_model = gbt.fit(train_data)
获得预测
dtc_preds = dtc_model.transform(test_data)
rfc_preds = rfc_model.transform(test_data)
gbt_preds = gbt_model.transform(test_data)
评估模型的性能
如前所述，我们的评估者将是中华民国。我们将初始化它的类，并将预测值传递给它，以便获得值。
my_eval = BinaryClassificationEvaluator(labelCol='>50K')**# Display Decision Tree evaluation metric**
print('DTC')
print(my_eval.evaluate(dtc_preds))
**# Display Random Forest evaluation metric**
print('RFC')
print(my_eval.evaluate(rfc_preds))
**# Display Gradien Boosting Tree evaluation metric**
print('GBT')
print(my_eval.evaluate(gbt_preds))
最佳预测器是梯度提升树。实际上，0.911 是一个非常好的值，当显示其预测值时，我们将看到以下内容:
提高模型性能
我们将尝试通过执行网格搜索交叉验证技术来做到这一点。有了它，我们将使用先前超参数值集的不同组合来评估模型的性能。
我们将调整的超参数有:
- 最大深度
- 最大箱数
- 最大迭代次数
**# Import libraries**
from pyspark.ml.tuning import ParamGridBuilder, CrossValidator**# Set the Parameters grid**
paramGrid = (ParamGridBuilder()
.addGrid(gbt.maxDepth, [2, 4, 6])
.addGrid(gbt.maxBins, [20, 60])
.addGrid(gbt.maxIter, [10, 20])
.build())**# Iinitializing the cross validator class**
cv = CrossValidator(estimator=gbt, estimatorParamMaps=paramGrid, evaluator=my_eval, numFolds=5)**# Run cross validations. This can take about 6 minutes since it is training over 20 trees**
cvModel = cv.fit(train_data)
gbt_predictions_2 = cvModel.transform(test_data)
my_eval.evaluate(gbt_predictions_2)
我们在预测性能上获得了微小的改进。而计算时间，几乎到了 20 分钟。因此，在这些情况下，我们应该分析这种改进是否值得努力。
结论
在整篇文章中，我们做了一个端到端的机器学习分类项目。我们还学习并获得了一些关于分类模型的见解，以及使用 PySpark、其方法和实现来开发一个具有良好性能的分类模型的关键。
我们还学会了一旦识别出一个性能良好的模型，如何调整我们的算法。
在接下来的文章中，我们将学习如何在 PySpark 中开发回归模型。所以，如果你对这个话题感兴趣，我强烈建议你继续关注！
如果你喜欢这篇文章，那么你可以看看我关于数据科学和机器学习的其他文章 这里 。
如果你想了解更多关于机器学习、数据科学和人工智能的知识 请在 Medium 上关注我，敬请关注我的下一篇帖子！
寻找多莉，隐马尔可夫模型和简化生活！
前几天和我的一个同事聊天，我刚刚意识到我可能是健忘的。或者至少比她“更多”的记忆——“更少”？在你开始把我的失忆和迪士尼海底总动员/多莉系列中可怜的多莉的失忆联系起来之前，T5 被描绘成(有点不公平和不准确地说一条实际上拥有良好记忆力的鱼！)，让我澄清一下，我以为我是无记忆的，尽管像少数概率分布。如果你面临概率挑战，不要担心——那是因为概率不在你这边(双关语！)我马上就要讲到了——但基本上她是在告诉我她不会忘记别人说过的事情或发生过的事情！？这引发了一系列关于生活和生活窍门的讨论，以及后来更多的反思和一连串的思考——我刚刚意识到，随着生活的进程变得势不可挡，在我无数的愚蠢和失败中，我经常拯救自己的东西——是我大脑的“无记忆”性质(或者更确切地说是思维过程！).像马尔可夫模型一样没有记忆！
那么这些马尔可夫模型是什么？马尔可夫模型是简单的系统(尽管他们的名字！)保持其状态有规律地变化。只是它每次都进入一个新的状态，只是基于现在，而不是看过去。当它改变状态时，它向世界表达(输出)一些信息。所以它去哪里和说什么纯粹是基于现在，使它没有记忆。事实证明，这些马尔可夫系统中的一些倾向于对它们的状态保密——虽然它在改变状态时会说些什么，但它从不透露它的状态。这些相当神秘的无记忆系统被称为隐马尔可夫模型(HMM)。下图显示了一个简单的例子。
A Simple Example of Hidden Markov Models to predict weather from clothing patterns
有趣的是，这些看似简单的健忘 hmm 其实是我们身边一些令人难以置信的技术的核心——当你对着亚马逊 Echo 设备喊“Alexa”或要求谷歌将句子翻译成一种新语言时。很长一段时间以来，这些语音识别或翻译系统都是由 HMMs 驱动的(我作为一名数据科学家也曾在此基础上开发过无数可能的应用程序，但我们不会以简化事物为目标！).
如果人工智能(AI)系统可以在没有太多过去包袱的情况下做得很好(顺便说一下，国家本身可以有很多关于过去的信息编码在里面！)，有没有我们人类可以借鉴的地方？我知道这听起来很奇怪，因为人工智能系统通常是在人类知识的基础上训练的，事实上人工智能系统现在正试图让自己“长期和短期记忆”就像人类更好地完成工作一样，我在说人类希望向无记忆人工智能系统学习。
嗯，我认为毫无疑问，世界正变得越来越复杂，我开始相信，当我们在这个每时每刻都被信息冲击的复杂世界中航行时，对人类的救赎在于简化。简化我们看待周围世界的方式。正如一句流行的名言所说，一个人无法控制外部世界，但你可以控制如何做出反应。
《T4 人》中的尤瓦尔谈到了 21 世纪人类的三大追求——神性、不朽和幸福。尤瓦尔还相当怪异地补充道，对幸福的追求可能是通过生化手段！虽然我认为这些预测似乎有一点可怕的必然性，但就目前而言，我敢说简化我们的思维系统是生存的关键，也是通往幸福的希望之路。
所以，忘掉记忆，忘记过去，在你的人生追求中，做好你目前的状态。事实上，对我们来说，要表现得善良，多一点宽容和不报复的本性，或者表现出毅力和耐心，就需要一种无记忆的元素。让我用一句哲学的话来结束我的演讲——当我们说时间治愈时，这仅仅意味着我们的记忆会随着时间消逝，痛苦也会消失。所以，你越快养成无记忆，你就能越快前进:)当你做决定时，不考虑过去的包袱，这确实简化了生活！
做马尔可夫！
PS:我已经决定在 14 年后再次写作，这几乎是一个自我放逐的时期，因为我曾经用泰米尔语写短篇小说和博客。现在，这是一个担心和后悔的原因吗？错过了记录无数的观察，推论和在此期间生命旅程的细节的机会吗？嗯，我是马尔科夫。我没有记忆；)
使用 Tensorflow 对象检测器、Pytorch 特征提取器和 Spotify 的 airy 查找熟悉的面孔
在一些不同的帖子中，我已经使用对象检测器和类似的暹罗网络将面部识别管道或类似的东西放在一起进行特征提取，以找到类似的图像，但我还没有真正深入研究如何以更实际的大规模方式使用这些特征向量。我在这里的意思是，你不希望在整个数据库中进行成对的比较，这是不实际的。所以在这篇文章中，我将演示如何使用 Spotify 的近似最近邻居库(骚扰)来根据一些初始图像找到相似的游戏角色。
Annoy 是 Spotify 开发的一个库，用于帮助大规模推荐音乐。它还有一些其他有用的属性，因为您可以预先计算索引，稍后可以调用这些索引来查找类似的项目，这在您进行大规模工作时很有帮助。它通过传递你所关心的任何东西(音乐、文本、动画人物的图片)的某种向量表示来工作，并基于此建立模型和生成索引。
这个管道的最简单的版本是使用一些神经网络进行特征提取，然后将这些特征向量传递给 airy，然后使用 airy 来查找相似的图像。这取决于您的用例。在这个具体的例子中，我想看看我是否会因为在不同的图像中返回相似或相同的字符而感到恼火。我对此进行了测试，并最终通过建立面部识别网络和传递面部生成的特征向量来尝试和集中管道以找到相似的字符，从而完成了两个阶段的方法。
对于这个项目的数据集，我将使用相同的命运大订单数据集，我已经在一些职位。主要是因为它是我最近建立的约 410 张图片的小数据集中最大的，所以它是这类事情的一个很好的测试案例。
进入这个项目时，我发现关键(像许多数据科学项目一样)是找到一种好的方法来表示我想要解决的问题的数据集。对于这种类型的任务，将 3D 像素阵列压缩成 1D 特征向量在机械上是正确的，但是确保最终的 1D 特征向量描述了您想要的也是重要的。
请随意查看我为这个在这里使用的笔记本。这些笔记本不是超级可读的，因为我很快就把东西处理完了。我也不像往常一样在 repo 中提供模型/数据集文件。
版本 1:完整图像特征提取
因为这是最容易做到的事情，并且是为后面的测试设置我的通用管道的好方法，所以我使用了 pytorch 预训练网络(一个 ResNet101，它为每个矢量生成 1000 个特征以输入到 aroy ),并向它传递完整的图像，然后由 aroy 使用这些图像来查找类似的图像。这个初始过程产生了如下结果。最左边的图像是基础图像，接下来的 4 个是“数据库”中最相似的 4 个。
Left most image is the seed image followed by 4 of the top ones. (there was a duplicate in the dataset for this one)
看起来特征提取器从原始图像中获得了许多类似的深色细节，其中红色和黑色是焦点。然而，在这种情况下，我试图看看我是否可以得到类似的字符，所以这不是最好的输出。
这是另一个图像示例。
Left most image is the seed image followed by 4 of the top ones.
这一个有点棘手，因为在主图像中有两个角色，但是当整个图像被传递到烦人的模型中时，相似的图像有相似的颜色味觉(可能…)。
因此，从机械上来说，这一过程是可行的，但我要说的是输出很差。下一步是想办法让输出更连贯。在这种情况下，因为我想建立的是一些显示类似的字符给定一些基础图像，为什么不尝试，使它的图像是基于只是字符，而不是完整的图像？我可以想到很多方法来做到这一点，但对我来说最简单的是建立一个单一的类对象检测器来识别图像中的人脸，将该人脸传递到一个预训练的网络进行特征提取，然后将该特征向量传递到 are。当骚扰拉相似的图像，我可以让它返回基本图像，所以基本上它会找到数据库中的哪些面部与新提取的面部最相似，并返回这些面部出现的图像。
基于人脸识别的相似度
对于这个项目，我建立了一个简单的 1 类对象检测器，只识别图像中的动画人物的脸。我使用 Labelimg 来快速注释图像，我非常确定我只花了大约 20 分钟来标记我的测试和训练分割的 400 张图像。在这一点上，我已经做了相当多的工作，只有一个单独的类可以大大加快这个过程。在训练探测器时，我使用了一个更快的 RCNN-Inceptionv2 模型，该模型具有来自 Tensorflow 模型 zoo 的 coco 权重，并训练该模型大约 3 个小时。我从周五午夜开始训练，一直到周六凌晨 3 点左右，这打乱了我的睡眠时间表，因为那时我还在忙其他的事情。
对象检测器训练得相当快，输出看起来相当干净。这是令人振奋的，因为这将是在这个管道中找到更多的角色特定的相似图像的关键部分。
Example of some cropped image outputs
当使用对象检测器从原始数据集中裁剪头部时，我保存了头部到其原始图像的 csv 映射。我的想法是，我可以在头像上运行一个特征提取器，并将其存储在骚扰模型中，到时候我可以将骚扰输出与原始图像进行匹配。
Pytorch 和 Annoy 的特征提取
现在我可以从图像中提取头部，我所要做的就是将这些头部通过特征提取器(再次是 ResNet101)，然后将这些特征向量传递给 aroy。
作为一个演示，这里有一张之前的图片，原始图片模型有一些问题。这是对象检测器检测图像中的两张脸的示例输出。因此，每张图像都有从中提取的特征，然后与更大的数据库进行匹配。
第一个输出来自主图像左边的角色(他出现在数据集中), 4 个相似图像的前两个和最后一个图像是该角色。这是对该图像的原始图像输入的改进，其中有 0 个匹配。
第二个角色(右边的一个)实际上并没有出现在数据库中…但是她的面部特征与左边的一个基本相同，所以四个图像 2 匹配(第一相似和第四相似)
Left is the base image, followed by 4 most similar
与只使用基本图像相比，这似乎是一个很好的改进，因为目标是返回相似的字符。现在让我们看看我之前为基本模型使用的另一个示例图像。
Left is the base image, followed by 4 most similar
因此，在这张照片中，他们看起来更有针对性，而不是仅仅得到一堆红色和黑色的图像。虽然第一、第二和第四个相似图像具有不同的特征，但是第三个图像具有相同的特征。这一次，所有人都至少与基本图像的性别相同。另一个只是让所有的男性角色和一个女性角色的基本图像配对。虽然这不是一个很好的结果，但似乎比以前的版本有所改进。
结束语
在查看了这两条管道的输出后，我觉得这些结果是可以接受的，但不是很好。使用基本图像返回具有相似感觉但不一定相似特征的图像。虽然面部检测器有助于将输出集中到相似的字符，但它通常不会返回总体相似风格的图像。
original full image model
虽然 4 个返回图像中的 2 个不是相同的字符，但我确实喜欢中间的第二个结果，因为它与基础图像有相似的“感觉”。
headshot based model
正如我之前提到的，基于头像的模型很好地聚焦于所讨论的角色。在这种情况下，所有字符都是相同的。然而，它不符合原始图像的感觉。我真正想要的是两者的某种结合，我可以得到相似的角色和相似的整体形象(基本上我自私地想要鱼和熊掌兼得)。
经过一些试验后，我发现我能够非常接近那个。
New model!
正如我在本文开始时提到的，从这条管道获得更好的输出基本上可以归结为修改哪些数据被压缩到最终的特征向量中，这些向量被传递到 annoy 中。虽然在大多数情况下鱼与熊掌兼得是不可行的，但在这种情况下却是可行的！我认为，这个“新模型”比其他两个做得更好，因为它获得了所有正确的字符(击败了基本图像模型)，并且显示的图像比头像模型更接近基本图像。
通常情况下，这只是一种我必须从新的角度解决问题的情况。
Still really enjoy big hero 6 and code to Immortals as a theme song for my life
我只是不得不重新思考我在最终的特征向量中编码了什么信息。我最终做的是将来自检测到的头像和基础图像的信息传递给一个组合的特征向量，该向量捕获了关于角色面部(以获得相似的角色)和基础图像(以获得整体“感觉”)的信息。
然而，这个最终模型并不那么简单，我花了一点时间才弄明白，所以我会再发一篇后续文章来保持这篇文章的合理长度。
因此，下次请继续收听，我将介绍我是如何将面部和基本图像信息组合成密集的表示形式，以让 Spotify 的 annoy 在性格和感觉方面找到相似的图像。
跟进博文 此处
再一次，在这里你可以随意查看我为这个使用的笔记本。这些笔记本不是超级可读的，因为我很快就把东西处理完了。我也不像往常一样在 repo 中提供模型/数据集文件。
寻找车道线-用于车道检测的简单管道。
识别道路的车道是人类驾驶员执行的非常常见的任务。这对于保持车辆在车道的限制内是很重要的。这对于自动驾驶汽车来说也是非常关键的任务。并且利用简单的计算机视觉技术，非常简单的车道检测流水线是可能的。本文将使用 Python 和 OpenCV 描述可用于简单车道检测的简单管道。这个练习是作为“uda city 无人驾驶汽车 nano degree”的一部分完成的。
请注意，该管道有其自身的局限性(如下所述),并且可以改进。改进将在后续文章中描述。
车道检测管道:
- 将原始图像转换为灰度。
- 使灰度图像变暗(这有助于降低道路变色区域的对比度)
- 将原始图像转换到 HLS 颜色空间。
- 将黄色从 HLS 中分离出来，得到黄色遮罩。(适用于黄色车道标线)
- 将白色从 HLS 中分离出来，得到白色面膜。(适用于白色车道标线)
- 逐位或黄色和白色掩码以获得公共掩码。
- 逐位并使用暗化图像进行屏蔽。
- 应用轻微高斯模糊。
- 应用 canny 边缘检测器(调整阈值—试凑法)获得边缘。
- 定义感兴趣的区域。这有助于剔除 canny 边缘检测器检测到的不需要的边缘。
- 检索霍夫线。
- 巩固和外推霍夫线，并绘制在原始图像上。
原始测试图像
Original Test Images
转换为灰度
将原始图像转换成灰度有它的好处。我们必须找到黄色和白色的车道，将原始图像转换为灰度可以增加车道相对于道路的对比度。
Original Images to Grayscale
使灰度图像变暗
这样做的目的是降低道路变色部分的对比度。
Darkened grayscale images
将原始图像转换到 HLS 颜色空间
原始图像是 RGB 的，但我们也应该探索其他颜色空间，如 HSV 和 HLS 。当并排看时，可以很容易地看到，我们可以从道路的 HLS 颜色空间中获得更好的颜色对比度。这可能有助于更好的颜色选择，进而有助于车道检测。
RGB vs HSV vs HLS
颜色选择
这里我们使用 OpenCV 的 inRange 来获取 thresh hold 值之间的掩码。经过反复试验，我们可以找出阈值的范围。
对于黄色面具:
- 色调值介于 10 和 40 之间。
- 我们使用较高的饱和度值(100–255)来避免山丘的黄色。
对于白色面具:
- 我们使用更高的亮度值(200-255)作为白色蒙版。
我们按位或两个掩码来得到组合掩码。
下图显示了组合蒙版是逐位的，并带有暗化图像。
高斯模糊
高斯模糊(高斯平滑)是预处理步骤，用于减少图像中的噪声(或平滑图像)。我们使用这个预处理步骤来移除许多检测到的边缘，并且只保留图像中最突出的边缘。
OpenCV 的 GaussianBlur 期望模糊的内核大小(奇数值)。在尝试了一些值之后，我们使用了 7。
Gaussian Blur
应用 Canny 边缘检测
现在我们将 Canny 边缘检测应用于这些高斯模糊图像。Canny 边缘检测是基于梯度变化检测边缘的算法。并不是说 Canny 边缘检测的第一步是使用默认的核大小 5 进行图像平滑，我们仍然在上一步中应用显式高斯模糊。Canny 边缘检测的其他步骤包括:
- 寻找图像的强度梯度
- 非最大抑制
- 滞后阈值
这个链接很好的解释了 Canny 边缘检测。
选择感兴趣区域
即使在应用 Canny 边缘检测之后，仍然有许多边缘被检测到，它们不是车道。感兴趣区域是一个多边形，它定义了图像中我们感兴趣的边缘所在的区域。
注意，图像中的坐标原点是图像的左上角。行坐标自上而下增加，列坐标自左而右增加。
这里的假设是摄像机保持在恒定的位置，并且车道是平坦的，因此我们可以“猜测”感兴趣的区域。
Region of Interest applied to Canny Edge Detected Images.
霍夫变换线检测
霍夫变换是一种通过识别直线上所有点来找出直线的技术。这是通过将线表示为点来实现的。并且点被表示为线/正弦曲线(取决于笛卡尔/极坐标系统)。如果多条直线/正弦曲线通过该点，我们可以推断这些点位于同一条直线上。
Line in polar co-ordinates and sinusoidal that represent point , intersection point represent line.
更多信息可在这里找到。
从感兴趣区域图像中找出 Hough 线，然后在原始图像上绘制 Hough 线。
Hough Lines
霍夫线的合并和外推
我们需要追踪完整的车道标志。为此，我们首先需要区分左车道和右车道。有一个简单的方法来识别左车道和右车道。
- **左车道:**随着列坐标值的增加，行坐标值减少。所以梯度一定是负的。
- **右车道:**随着列坐标值的增加，行坐标值也增加。所以梯度一定是正的。
- 我们将忽略垂直线。
在识别左车道和右车道霍夫线之后，我们将外推这些线。我们做两件事:
- 车道检测到多条线，我们将对这些线进行平均
- 有一些部分检测到的线，我们将推断它们。
将管道应用于视频
现在让我们将这个管道应用于视频。
管道在直道上非常好用。
但是这对于弯曲的车道来说效果不好。
缺点
- 基于直线的霍夫线不适用于弯曲的道路/车道。
- 有许多试错法来获得正确的超参数。此外，感兴趣的区域假设摄像机停留在相同的位置，并且车道是平坦的。因此，在决定多边形顶点时涉及到“猜测”工作或硬编码。
- 一般来说，有许多道路可能没有车道标志，这是行不通的。
未来的改进
- 代替直线，使用在弯曲道路上有用的更高程度的曲线将是有益的。
- 即使当我们使用来自先前帧的信息时，仅仅平均车道可能不是很好的策略。可以是加权平均值，或者某种形式的优先级值也可以。
该代码可在 GitHub 上获得
如果你喜欢这篇文章，请一定给我一个掌声，并关注我以获得我未来文章的更新。
还有，随时在LinkedIn上联系我或者在Twitter上关注我。
如果你喜欢我的作品，请考虑赞助我，它会帮助我推出更多这样的作品。
在芝加哥寻找低风险餐厅享受美食
Photo by Robin Stickel on Unsplash
作为一个美食爱好者，我总是对探索餐馆很感兴趣。因此，在这个项目中，我决定满足自己的渴望，处理食物关节数据。我建立了这个项目的目标，根据给定的位置探索芝加哥的餐馆。我去了 Kaggle，找到了完美的数据集。
来自芝加哥市公开数据
www.kaggle.com](https://www.kaggle.com/chicago/chicago-food-inspections)
我开发了一个 Kaggle 内核，在那里我使用各种可视化工具探索了数据集，然后创建了一个方法来获取距离当前位置(使用纬度和经度定义)一定数量公里附近的餐馆。
使用来自多个数据源的数据
www.kaggle.com](https://www.kaggle.com/bhanotkaran22/facilities-in-chicago-risk-and-result-analysis)
导入库和数据集
第一步是导入将要使用的库。我导入numpy和pandas进行数据操作。为了绘制图表和地图，我使用了plotly。注意，为了使用地图，我们需要一个访问令牌。我把令牌保存在input目录下的一个方法中，在执行过程中把它复制到working目录，使用它，最后会删除它。
最后，我导入数据集并将其保存到dataset变量中。有很多列。要开始了解每一列的值，请参考 Kaggle 数据集的Data部分并向下滚动。
仪表板可见。在右上方，可以看到Views。选择第三个选项，你会看到每一列的名称、类型、发行版等等。
数据分析
第一步是识别所有设施，并对每个设施进行最近的检查。我还将删除所有Risk、Facility Type、DBA Name、Latitude或Longitude为空的行。一些企业已经停止运营或不存在，因此也可以被移除。我将创建一个新列Name，它从AKA Name和DBA Name中提取姓名，并优先选择AKA Name。Name还将包括每个工厂的地址。
latest_data是我将在余下的探索中使用的数据集。
我还将为绘图中使用的颜色创建列。All、Risk 1 (High)、Risk 2 (Medium)和Risk 3 (Low)的风险颜色分别为黑色、红色、黄色和绿色。检查时，Pass和Pass w/ Conditions设置为绿色，而Fail、No Entry和Not Ready设置为红色。
所有的设置现在都完成了。让我们继续对我们的数据进行可视化。
数据可视化
设施类型
首先，我将为这次食品检查中包括的各种类型的设施绘制一个饼图。有一些设施不到 1%，因此，我会将它们全部组合在一起，并将其视为Others。
Plotly饼图提供了许多样式选项，如高度、宽度、内圆大小等。
大约 55%的设施是餐馆。
风险分析
接下来，我将基于上面创建的Risk Color列在芝加哥地图上绘制设施列表。
看来红点的数量在图中确实很高。这意味着大多数设施处于高风险状态。我们也可以从数据中证实这一点。
高风险贷款的数量多于中低风险贷款的总数。
成功和失败
接下来，我将探究最近的检查中有多少是Pass和多少是Fail。我认为Pass w/ Conditions是好的，而No Entry和Not Ready是坏的。
正如我们从上面的柱状图中所看到的，尽管许多设施处于高风险中，但它们中的大多数也通过了检查。
最大检查结果为
Pass。
地图上的检查结果如下图所示。
大多数设施通过了检查，只有极少数没有通过。
逐区分析
每个设施都有一个与之相关联的Ward。我根据病房对所有设施进行了分组，然后在柱状图上绘制了合格与不合格数。
似乎 36 号病房通过检查的次数最多。
大多数检查是针对 36 号病房的，它具有最大的通过率。
多数违反
对于许多餐馆来说，在一定数量的检查中存在违规行为。因此，我决定在所有违规行为中找出关键的违规行为。
每个违规都作为句子添加到Violations栏的一个段落中。我必须从每个列值中提取所有的数字，创建一个数字列表，并计算数值计数。完成后，我选择了超过 100 的违规，并绘制了一个条形图。
最大违规发生在 41 号，它没有正确使用和存放抹布。
所有的可视化现在都完成了，让我们创建一个方法来找到我们的餐馆。
探索附近的餐馆
我将创建一个方法，该方法将接受一个位置列表，用户的当前位置(以纬度和经度表示)，您想要搜索的风险级别，以及需要搜索的餐馆的距离(以公里为单位)。
让我们首先创建一个方法，根据两个位置的纬度和经度值来获取它们之间的距离。
完成后，我们可以使用上面使用的代码来绘制所有设施点，但有基于数据值的限制。
我已经添加了默认值，以防某些值没有提供。在该方法中，我遵循以下步骤:
- 过滤位置以仅包含餐馆
- 根据
risk_level选择要搜索的风险级别。这可能是Low、Medium、High或All - 过滤搜索距离内的餐馆并重置索引
- 在地图上标出所有点
我使用以下命令来运行该方法:
get_plot(latest_data, 41.8781, -87.6298, 'Medium', 5)
我取latest_data中所有最近的数据，设定搜索中心为芝加哥中部，风险等级为Medium，搜索距离为 5 公里。
我们可以看到，这些地点被标为绿色，表示低风险，黄色表示中等风险。它看起来像一个半径为 5 公里的半圆。您可以调整搜索参数以获得想要的结果。
结论
在本文中，我探索了数据科学的数据可视化方面，以及如何用它来以一种简单的方式理解数据。
您可能还喜欢:
[## 使用 Python 中的多个绘图库在地图上绘制商业位置
比较地图打印库
towardsdatascience.com](/plotting-business-locations-on-maps-using-multiple-plotting-libraries-in-python-45a00ea770af) [## 使用 React 和 Flask 创建一个完整的机器学习 web 应用程序
用 UI 和后端构建完整的 ML 应用程序
towardsdatascience.com](/create-a-complete-machine-learning-web-application-using-react-and-flask-859340bddb33) [## 使用机器学习预测心脏病的存在
机器学习在医疗保健中的应用
towardsdatascience.com](/predicting-presence-of-heart-diseases-using-machine-learning-36f00f3edb2c) [## 在 Kaggle 上使用 CNN 进行深度学习以识别疟疾细胞
医学领域的深度学习
towardsdatascience.com](/deep-learning-to-identify-malaria-cells-using-cnn-on-kaggle-b9a987f55ea5)
我希望你喜欢这篇文章。请和我分享你的想法和建议。我很乐意收到你的来信。
发现意义:小区域分析的情况和数据可用性的进展
作为纽约市的居民和教师，我的经历促使我着手一个研究项目，对纽约市的街区进行小区域分析比较。当时我并不知道，我最初的尝试发生在早几个月的时候，在之前的一篇中有详细描述。
接近并不意味着相似
你不必相信我的话。首先，看看这两张在纽约大通银行同时发现的收据就知道了。
I remembered seeing the photo circulating social media while I was living there, and had to dig it up for this article. I can’t access the original post through Facebook anymore, but the photo was extremely easy to find even after all these years — I still remembered the caption it had: “a tale of two New Yorks.”
一位匿名银行家是百万富翁，另一位的账户里连 100 美元都没有。他们在同一天使用了同一个自动取款机！考虑到纽约市的经济不平等，这并不奇怪，但这是这种不平等的明显表现。
小区域分析的必要性和好处
Economic, education, and health needs in New York City neighborhoods from Case Studies — Scientific Figure on ResearchGate
在小范围内发现有意义的差异有明显的好处。小区域分析通常涉及分析大型数据库，以比较通常被视为整体的群体的子集内的资源可用性和使用情况。在这种情况下，整个就是组成纽约市的五个区。这个想法是发现差异有助于确定社区的需求。这将最终导致资源的重新分配或调查为什么分配的资源没有被一个子群体使用，并找到一种方法使它们更容易获得。
小区域分析在历史上一直是一件很难实现的事情，尤其是对于我所追求的细粒度。可靠的、有统计学意义的数据很难获得。这种情况正在改变，纽约市是带头冲锋的城市之一，因为该市承认，其巨大的多样性给平等带来了巨大的挑战，而且它有资源应对其中的一些挑战。
开放数据来到纽约
2012 年，纽约市通过了“开放数据法”，要求到 2018 年底，所有市政府的公共数据都可以通过一个单一的在线门户免费获取。
世界上最大的城市之一的公共记录的大规模发布是迈向数据民主化的一大步。纽约人已经将这种资源很好地用于 T2 的 T3 航站楼。
http://storymaps.esri.com/stories/2016/wealth-divides/index.html
一年前，认输并比较各区要容易得多，但我坚持捕捉较小社区之间的差异——这在当时导致了一些相当不起眼的分析。这一次，由于开放数据倡议，这将容易得多。虽然不是所有的数据对我分析小于全市或区范围的区域都有用，但大部分是有用的。
纽约市的怪癖
如果你从未逛过纽约，可能很难理解为什么我这么想“放大”纽约的街区。我是说，至少布鲁克林是布鲁克林，布朗克斯是布朗克斯。有哈莱姆区，某处有唐人街，中央公园周围有价值百万的公寓。但仅此而已，对吧？
嗯……数据科学家徐仅在布鲁克林的威廉斯堡社区就能识别出 75 个不同的社区。她通过分析纽约市出租车和豪华轿车委员会的数据做到了这一点，正如她在这篇文章中概述的那样。诚然，这些“社区”中的大多数并不居住在布鲁克林，小区域分析通常集中在居住地，75 仍然很多。这是她的地图，显示了她所识别的五个社区的出租车数据。
Five Communities of Williamsburg
在纽约住了十年之前，我不知道街区变化有多快(跨越时空)。我想我从电影和电视节目《老友记》中获得了关于这座城市的大部分信息。虽然我知道这些是不切实际的描绘，但当我第一次冒险进入这座城市的不同地方时，它仍然让我感到惊讶。
不到半英里外的世界
2015 年，我住在离纽约长老会医院两个街区远的地方，那里是哥伦比亚大学医学中心的所在地，当时我正在和一名住院医生约会。这套公寓是新装修的，保养得很好，地下室有洗衣房。这个社区安静、高档、专业，有相当多的种族多样性。地铁站就在我住的街区的尽头，我步行一分钟就能到达附近的一家墨西哥小吃店和一家星巴克，附近还有一些不错的餐馆，虽然价格有点高。
Economic, education, and health needs in New York City neighborhoods from Case Studies — Scientific Figure on ResearchGate.
当住院医生和我分手时，我是那个必须搬家的人，因为我们住在哥伦比亚大学的宿舍里。这很容易，我最好的朋友独自住在八个街区外，我的新公寓每月便宜 700 多美元，尽管哥伦比亚的住房有部分补贴。赢。
这个社区也很有趣。除了街对面的塔可钟之外，在几分钟的步行距离内没有任何大型连锁餐馆。有一些非常好的波多黎各和多明尼加小餐馆，不过，大多数都很实惠。无论白天黑夜，任何时候都会有一群人坐在外面的小凳子或折叠椅上演奏加勒比音乐。哦，这个社区几乎完全是拉丁裔的 T4 人。街上听到的西班牙语比英语多。我的朋友是波多黎各人，所以她感觉非常自在。另一方面，我绝对是周围唯一的亚裔美国人和白人。我确实会说一些西班牙语，但是是轻微的美国口音加上墨西哥口音/方言。加勒比西班牙语几乎是另一种语言…**
无论如何，更便宜的房租也意味着，总的来说，这个社区更穷。我现在正好在两个地铁站的中间，从大门到城市的其他地方大约需要走 7-8 分钟。相信我，在炎热潮湿的夏天和寒风凛冽的冬天，这几分钟非常重要。我注意到我每天见到的警察单位数量急剧增加。街道更加肮脏，建筑物更加破败。我在通勤开始和结束时看到的人似乎都很疲惫。这似乎是一件奇怪的事情，但我也累了——因为我每月节省的 350 美元(700 美元中我的份额)是有代价的。
This is the block I actually lived on. On the right a few buildings down is my old apartment building. The restaurant in the center, was one of my go-to spots. This photo was among the top twenty photos when I googled “Washington Heights.” Unfortunately, it belongs to a news article about a stabbing.
我的公寓里有坏掉的炉子和烤箱，一个我们必须定期维修以便冲洗的厕所，墙上有大量的蟑螂和老鼠，一扇通往消防通道的窗户没有锁好，一个散热器将蒸汽泄漏到整个公寓，使它成为一个桑拿浴室，顺便说一下，也为黑色霉菌提供了理想的条件。如果问题很容易解决，大楼的负责人就很好，否则就不可避免地会很忙。此外，洗衣服从去地下室的几次旅行变成了去四个街区外最近的自助洗衣店的两个小时的旅行。
请记住，从高档社区到普通建筑的巨大变化都发生在八个街区之内。就规模而言，通常认为大约 20 个纽约市街区等于 1 英里，平均一个街区大约需要步行一分钟。我还没有核实，但听起来是对的。
*
**
A definitely-not-completely-accurate-but-sometimes-funny-sometimes-just-weird takes on Brooklyn and Manhattan neighborhoods by Urbane demonstrates how neighborhoods have identities unto themselves.*
在布鲁克林教书
我在纽约教了四年高中数学、科学和计算机科学，我在一所公立特许学校工作，它是多元化特许学校联盟的成员，所以我的课堂是社会经济、种族和民族多元化的。很自然，我的不同学生从纽约市的许多不同社区来到学校。虽然大多数学生来自布鲁克林，但包括斯塔滕岛在内的每个区都有代表。一些学生的通勤单程时间超过两个小时。
面对全国重新种族隔离的趋势，我仍然钦佩学校对多样性的承诺。然而，在那里教书让我面对了一些令人不安的事实。我让一些学生带着所有最新的小玩意来学校，去街对面昂贵的咖啡店闲逛，而其他学生则报名参加任何周末聚会的俱乐部(比如我们的野营俱乐部)，因为这意味着他们可以在周末得到食物。纽约市各社区的经济差异非常显著，学生在学校的成绩往往与这些差异相关。
虽然整合学校可能是所必需的*，以消除《不让一个孩子掉队法案》( 2001 年由老布什签署)和《力争上游》( 2009 年由奥巴马签署)试图消除的成绩差距，但这还不够。*
当在学校学习、参与和集中注意力的能力依赖于一系列其他因素时，高质量的教师和学校资源只能做这么多，其中许多因素由学生的家庭和社区决定或与之相关。
不同社区之间的不平等可能是巨大的，甚至在小区域内也是如此。小区域分析是一种工具，有助于全面了解不同社区的需求，并在制定公共政策和资源分配决策时提供指导。虽然纽约市在这方面并不独特，但它是一个这些不平等经常被放大的城市，而且似乎是一个现在转向公开解决这些问题的战略的城市。这个开放数据项目将如何发展将会很有趣…
寻找意义:社会科学中小区域预测分析的数据问题(下)
如果你还没有看过我这个项目的第一部分，请随意看看。虽然没有必要理解我在这里提出的问题，但我解释了什么是小区域分析，我对这个项目的灵感，我如何收集和准备我的数据，以及我在做这件事时面临的困难。这是一项正在进行的工作，在接下来的几天里将会有这个系列的第 3 部分。
本文原载于 2018 年 8 月。这是原物稍加编辑后的拷贝，被取了下来。
在这里，我将谈谈我在对纽约市五个区的 59 个指定社区的公共卫生、环境和教育数据建模时遇到的困难。
Community Districts of New York City
我将在这篇文章中讨论这三个主题。
- 我考虑和尝试的不同模型
- 有限数据的过拟合问题
- 使用机器学习进行解释和分析，而不是预测
甚至在开始使用模型之前，我就发现了每个特性与教育程度列之间的相关性。有些很有趣。看看不同社区的健康相关指标——尤其是“含糖饮料”的相关性。
在某些情况下，空气质量和“社区健康”指标(一个社区是否繁荣的指标)也有很强的相关性。
最后，种族与国家比率有着同等的相关性…
数据建模
开始时最明显的模型是简单的线性回归。虽然可以将 scikit-learn 的LinearRegression用于多个输出(在这种情况下，没有高中文凭的成年人的比率、有高中文凭(可能有大学文凭)的成年人的比率以及有大学或更高学位的成年人的比率)，但该模型假设输出是独立的，因此基本上创建了三个完全独立的线性回归模型，每个模型对应一个输出。
我继续做了，完成了对所有数据的线性回归，得到了 0.95 的 R 值(三个输出的平均值)。这一切都很好，但有两个突出的问题…
第一个问题是，三个产出显然不是独立的；事实上，它们是如此紧密地联系在一起，以至于在每种情况下，这三个比率的总和应该是 100(这是百分比的趋势)。因此，创建三个完全独立的线性模型实际上并不能准确描述情况。
seaborn pairplot of the educational attainment data
第二个问题是，虽然这个模型似乎告诉我们，模型中包含的特征似乎确实与教育成果数据相关，但当我将我的数据分为训练集和测试集时，无论大小如何，我总是得到训练数据大于 0.95 的 R 值，而测试数据通常是负值——我得到的最高值是 0.49，这完全是基于数据的随机分割的运气。所以这里肯定存在过度拟合的问题。只有 59 个样本，因为只有 59 个社区区，这似乎是一个不可避免的问题，无论模型。
以几种不同的方式在第一个相关教育成就特征的问题上取得进展。首先，我研究了可以处理多个输出的模型，同时考虑到它们是相关的。唯一有希望的途径似乎是尖端的神经网络。然而，这些只是在研究论文中提到过，还没有生产模型，所以在这一点上它真的只是在理论阶段…所以对我来说是一个死胡同。
唯一的另一种方法是将我的三个输出合并成一个单一的分数。最简单的方法是给每个教育程度分配一个值:0 代表没有完成高中学业，1 代表高中学历和大学学历，2 代表大学学历或更高。对于每个地区，我将每个百分比乘以其相应的值，然后将它们加在一起。这为每个地区创建了一个单一的，虽然有点抽象的“教育分数”，然后我可以用它来运行一个单一的线性回归。
educational_attainment['edu_score'] = educational_attainment.apply(
lambda row: row['Eduhsdegreeorsomecollege_rate'] + 2 *
row['Educollegedegreeandhigher_rate'], axis=1)
在整个数据集上，这种线性回归产生了 0.96 的 R 值；然而，当分成训练集和测试集时，它陷入了同样的过度拟合陷阱。
有序回归
我并不期望过拟合问题会有更好的结果，所以决定使用mord模块尝试一些有序回归。所以我用我的教育成绩，用平均值和标准差把数据分成 4 类。
# Get the mean and standard deviation (indices 1 and 2 respectively)
statistics = educational_attainment['edu_score'].describe()# Define my categories
def get_ordinal(row):
if row.loc['edu_score'] >= (statistics[1] + statistics[2]):
ordinal = 1
elif row.loc['edu_score'] >= (statistics[1]):
ordinal = 2
elif row.loc['edu_score'] >= (statistics[1] - statistics[2]):
ordinal = 3
else:
ordinal = 4
return ordinal
educational_attainment['ordinal'] = educational_attainment.apply(
lambda row: get_ordinal(row), axis=1)
令我惊讶的是，当我使用mord.LogisticIT()(使用即时阈值变量的有序逻辑回归)时，这在过度拟合方面稍微好一些。当作为一个整体用于数据集时，R 值为 0.93。当分成训练集和测试集时，测试数据集的 R 值至少是总是正值，有时高达 0.5——仍然没有什么值得夸耀的，但却是一个相当大的改进。
对抗这种过度拟合的另一种方法可能是找到不同年份的相同数据，或者找到更多组邻域的相同数据。然而，看到找到这个数据有多么困难，我将把它留到以后再说。因此，虽然我们似乎可以找到相关性，但该模型的任何预测能力本质上都是无意义的，因此没有因果关系可以合理建立。
更有希望的途径似乎是使用 L2 正则化线性回归(也称为岭回归)来计算特征重要性，以伴随我们的原始相关性，以追求理解和解释数据，而不是试图建立一个纯粹的预测模型。对这种分析有用的其他模型是随机 L1 (Lasso)回归和 scikit-learn 的额外树回归。以下是我在实施这些策略时发现的结果。
岭回归、随机套索和额外树
这些模型是解释数据的绝佳工具，因为这些模型是稳定的并且有用的特征往往具有非零系数/分数——这是算法中的冗余和平均试验和/或惩罚高系数以确保特征不会被过度呈现的结果。
作为另一层，我使用了包含对sklearn的支持的eli5包来实现PermutationImportance，其中，一列接一列的值被打乱，以查看模型仍然可以预测目标的程度。
# Using robustly scaled columns...
def get_scores(X, Y, scale_method):
index_tuples = []
model_data = []
for level in Y.columns:
ridge_model = RidgeCV()
et_model = ExtraTreesRegressor(n_estimators=50,
bootstrap=True)
randL_model = RandomizedLasso()
models = [ridge_model, et_model]
y = Y[level] for model in models:
model.fit(X, y)
score = model.score(X, y)
model_name = f'{model}'.split('(')[0]
try:
coefs = model.coef_
except:
try:
importances = model.feature_importances_
except:
importances = np.array(None)
else:
importances = np.absolute(coefs)
finally:
perm = PermutationImportance(model).fit(X, y)
perm_importances = perm.feature_importances_
index_tuple1 = (level, 'importances', score,
model_name)
index_tuple2 = (level, 'perm_importances', score,
model_name)
if importances.any():
index_tuples.append(index_tuple1)
model_data.append(importances)
index_tuples.append(index_tuple2)
model_data.append(perm_importances)
当我对这些不同模型的得分进行平均，并筛选出作为社区和环境健康指标的特征时，我得到了下面的结果。
The most important factors of Environmental and Community Health in predicting Educational Attainment of adults in the community. “Extremes” is an aggregation of rates for “Did Not Complete HS” and “College Degree or Higher”
在继续这项研究的过程中，我计划通过查看相关性和 p 值来选择适当的特征，从而减少特征的数量，希望将过度拟合的影响降至最低。我还计划继续寻找其他数据来源——尤其是可以显示随时间变化的数据。这是一项正在进行的工作。生活也是如此。
用 Apache Spark 大海捞针
“HearU Music” by Yi Mao is licensed under CC BY-NC-ND 4.0
TL；DR: 客户流失对企业来说是一个真实的问题，在不断增长的(大)数据中预测哪个用户可能会流失可能很困难。Apache Spark 允许数据科学家在大数据中轻松地进行大规模数据清理/建模/预测。
大多数数据科学家都知道处理数据并不总是一帆风顺。像清理、输入缺失值、特征工程、建模和预测这样的过程，即使数据小到可以放入笔记本电脑的内存中，其本身也是一个巨大的问题。如果数据比这大得多，事情很容易变得更复杂。解决这个问题的一个非常常见的方法是将数据放在 SQL 或非 SQL 数据库中，在汇总数据并将汇总的数据移动到本地工作站进行建模之前，在那里进行大部分的争论/清理。
然而，有时数据科学家需要大量数据输入到他们的模型中，并根据大数据进行训练和预测。对于传统的库，如 Python Pandas、scikit-learn 或 R dplyr，这并不容易，因为我们只有有限的内存来容纳。
Apache Spark 是大数据生态系统中最大的明星之一。它允许数据科学家使用熟悉的工具工作，但允许 Spark 做所有繁重的工作，如并行化和任务扩展。它提供了类似 Spark 数据帧的工具，类似于 R 数据帧或 Pandas 数据帧。如果您更喜欢传统的 SQL，您可以使用 SQL 来处理数据，而不是使用数据框。Spark 通过 MLlib 库支持许多开箱即用的机器学习算法。它还通过 Spark Streaming 为数据工程师提供流支持，最后，它通过 GraphX 支持原生图形处理。Spark 是数据科学家/工程师在处理大数据时的一把瑞士军刀。
在本帖中，我们将研究一个商业案例，这是许多数据科学家的共同任务，因为它对商业营销或战略工作有非常直接的影响。我们将尝试预测可能会为音乐流媒体平台流失的用户。基于良好的预测，通过促销或折扣来保持用户的参与，以流失用户，这一点至关重要。
在本例中，我们有完整数据集(12GB)的一个小子集(128MB)。由于我在本地集群上工作(仅指我的笔记本电脑，而不是一组服务器)，我将对小数据集进行分析，但我们将探索的所有方法都适用于较大的数据，没有什么根本不同，只是 Spark 将在大多数情况下处理并行性。
我们使用的数据来自 Udacity 数据科学家纳米学位计划 Apache Spark 顶点计划
Import necessary libraries. We will import some of them again for clarity later
Spark 需要一个叫做 Spark Session 的东西，这是你的代码和主节点之间的驱动。让我们创建一个，如果还没有的话，或者如果已经有了就返回它。
Creating/getting Spark session. Note that our cluster is a local one.
加载和清理 Sparkify 数据集
我们将用户活动数据集mini_sparkify_event_data.json导入 Spark，然后清除一些缺失的值。
Read data into Spark and do few simple checks on data
让我们浏览一下前几行，感受一下我们的数据
data.take(5)[Row(artist='Martha Tilston', auth='Logged In', firstName='Colin', gender='M', itemInSession=50, lastName='Freeman', length=277.89016, level='paid', location='Bakersfield, CA', method='PUT', page='NextSong', registration=1538173362000, sessionId=29, song='Rockpools', status=200, ts=1538352117000, userAgent='Mozilla/5.0 (Windows NT 6.1; WOW64; rv:31.0) Gecko/20100101 Firefox/31.0', userId='30'), Row(artist='Five Iron Frenzy', auth='Logged In', firstName='Micah', gender='M', itemInSession=79, lastName='Long', length=236.09424, level='free', location='Boston-Cambridge-Newton, MA-NH', method='PUT', page='NextSong', registration=1538331630000, sessionId=8, song='Canada', status=200, ts=1538352180000, userAgent='"Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/37.0.2062.103 Safari/537.36"', userId='9'), Row(artist='Adam Lambert', auth='Logged In', firstName='Colin', gender='M', itemInSession=51, lastName='Freeman', length=282.8273, level='paid', location='Bakersfield, CA', method='PUT', page='NextSong', registration=1538173362000, sessionId=29, song='Time For Miracles', status=200, ts=1538352394000, userAgent='Mozilla/5.0 (Windows NT 6.1; WOW64; rv:31.0) Gecko/20100101 Firefox/31.0', userId='30'), Row(artist='Enigma', auth='Logged In', firstName='Micah', gender='M', itemInSession=80, lastName='Long', length=262.71302, level='free', location='Boston-Cambridge-Newton, MA-NH', method='PUT', page='NextSong', registration=1538331630000, sessionId=8, song='Knocking On Forbidden Doors', status=200, ts=1538352416000, userAgent='"Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/37.0.2062.103 Safari/537.36"', userId='9'), Row(artist='Daft Punk', auth='Logged In', firstName='Colin', gender='M', itemInSession=52, lastName='Freeman', length=223.60771, level='paid', location='Bakersfield, CA', method='PUT', page='NextSong', registration=1538173362000, sessionId=29, song='Harder Better Faster Stronger', status=200, ts=1538352676000, userAgent='Mozilla/5.0 (Windows NT 6.1; WOW64; rv:31.0) Gecko/20100101 Firefox/31.0', userId='30')]
每行标识一个用户活动，如当前艺术家、歌曲、会话 Id、用户 Id，以及用户正在从哪个设备收听、unix 时间戳、性别、用户是否登录和付费用户、用户详细信息…等等
数据清理
我们应该检查数据集中的空值。根据结构，我们可以省略一些列或估算(如果我们的分析需要的话)。让我们从统计和视觉两方面来看看
Checking Null and NaN values in our dataset
我们有大量的空值，特别是在某些列中，比如artist、length或song。
我们还观察到，在列组之间的空值计数方面可能存在相关性，例如artist、length或song都具有相同数量的空值。对于firstName、gender、lastName、location、registration、userAgent都是如此。
其他列都没有任何缺失值。
同样，让我们直观地检查缺失值，看看我们之前的相关性声明是否得到支持。如果这是真的，我们可以说不仅空计数是相同的，而且空值出现的行数也应该是相同的。
由于 pyspark 没有可视化库，我们采样 spark 数据帧并将其转换为 pandas 数据帧，然后使用 Seaborn 进行可视化。
plt.figure(figsize=(12, 12))
sns.heatmap(data.sample(False, 0.1, 42).toPandas().isnull())
Using Python Seaborn heatmap we can visualise NULL values in the dataset
上面的热图支持我们的说法，列在缺失值、数字和索引方面是相关的。
我想检查的最后一件事是，如果firstName和其他类似字段为空，那么artist和类似字段是否也为空？
sns.heatmap(data.filter(data.firstName.isNull()).toPandas().isnull())
Filelds like artist, location… etc all are null when fistName is null
是的，它们是相关的，所以不仅空计数相等的组内部相关，而且组间空模式也相关。用简单的英语来说，如果firstName为空，那么gender也为空，因为它们的空计数相等，但是上面的图表也表明artist col 也将为空
根据下面的两个表，当用户注销或作为来宾收听时，firstName、gender、lastName、location、registration、userAgent和auth列具有空值。因为我们真正感兴趣的是那些已经注册并可能流失的用户，所以我们可以删除这些行，因为它们不是登录用户。
data.filter(data.firstName.isNull())\
.groupBy('firstName','gender', 'lastName',
'location', 'registration','userAgent', 'auth')\
.count()\.show()
data.groupBy('auth').count().show()
Rows where firstName … etc are null they all are either Guest or Logged out.
# Drop the rows with null values for the columns below
data = data.na.drop(subset=['firstName','gender',
'lastName','location','registration',
'userAgent', 'auth'])sns.heatmap(data.toPandas().isnull())
Removing logged out users gives us a cleaner dataset
现在，我们有了一个更清晰的用户登录数据集。
ts 列是 unix epoch 格式的时间戳，我们最好转换成时间戳格式，并添加到timestamp列中
data = data\
.withColumn('timestamp', from_unixtime(data.ts/1000))data.take(2)[Row(artist='Martha Tilston', auth='Logged In', firstName='Colin', gender='M', itemInSession=50, lastName='Freeman', length=277.89016, level='paid', location='Bakersfield, CA', method='PUT', page='NextSong', registration=1538173362000, sessionId=29, song='Rockpools', status=200, **ts=1538352117000**, userAgent='Mozilla/5.0 (Windows NT 6.1; WOW64; rv:31.0) Gecko/20100101 Firefox/31.0', userId='30', **timestamp='2018-10-01 00:01:57'**), Row(artist='Five Iron Frenzy', auth='Logged In', firstName='Micah', gender='M', itemInSession=79, lastName='Long', length=236.09424, level='free', location='Boston-Cambridge-Newton, MA-NH', method='PUT', page='NextSong', registration=1538331630000, sessionId=8, song='Canada', status=200, **ts=1538352180000**, userAgent='"Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/37.0.2062.103 Safari/537.36"', userId='9', **timestamp='2018-10-01 00:03:00'**)]
对于我们的模型，我们将删除一些我们不需要的字段，可能它们在其他列中有重复的信息或者只是不相关。在这一点上，我决定不在我的模型中使用它们，但是由于模型构建是一个迭代过程，如果我相信添加这些列会有所不同，我可能会在以后决定不删除这些列。
cols_to_drop = ['userAgent', 'artist', 'firstName',
'lastName', 'location', 'song',
'ts', 'registration', 'length',
'method', 'status']data = data.drop(*cols_to_drop)data.show()
User music after adding **timestamp** and removing irrelevant columns for prediction
下表告诉我们，免费和付费用户都可以取消他们的服务，但只有付费用户可以降级服务。嗯，有道理。我相信我们可以进一步过滤我们的数据框架来清理几行，因为像设置或帮助这样的页面不会告诉我们太多关于客户流失的信息。
Pages vs User Level. Next Song is the dominant as expected. There are 52 Cancellations (churns)
pages_to_keep = ['NextSong', 'Downgrade',
'Cancellation Confirmation', 'Upgrade',
'Add to Playlist', 'Cancel', 'Submit Upgrade', ]data = data.filter(col('page').isin(pages_to_keep))data.groupBy('auth', 'page').count().orderBy('auth').show()
After filtering out pages like Help, we are left with these above. Cancellation Confirmation is the churn KPI
探索性数据分析
Cancellation Confirmation是对流失最明确的 KPI。Downgrade页面被访问了数百次，但是数据并没有清楚地告诉我们用户在访问该页面后是否被降级，比如我们没有Downgrade Confirmed后面跟有Downgrade页面。此外，数据集中的Downgrade行数远远大于唯一用户数。因此，我坚持使用Cancellation Confirmation作为客户流失指示器。
# In our final dataset we do not have any null values
data.select([count(when(col(c).isNull(), c)).alias(c)
for c in data.columns])\
.head()\
.asDict(){'auth': 0,
'gender': 0,
'itemInSession': 0,
'level': 0,
'page': 0,
'sessionId': 0,
'userId': 0,
'timestamp': 0}
我们的数据集中有多少不同的用户？
data.select('userId').distinct().count()>> 225
让我们来看看一个随机用户的听力活动。 *itemInSession* 连续收听时单调递增。我们可以提取每个会话的最大值，以了解用户参与度(稍后在功能工程部分)
data.filter(col('userId') == 30).show(20)
User activity for user with userId 30. Note that itemInSession increases when the page changes. sessionID stays the same for a continuous music streaming session.
让我们简要地看一下被搅动的用户/订户。看起来free和paid的客户都会流失。我们将这些数据保存到churn_data中，因为我们将在以后加入特征工程部分时使用这些数据
churn_data = data.filter(col('page')=='Cancellation Confirmation')
churn_data.show(10)
List of churned users, there are both free and paid customers. We store the userId for joining later with other summarised dataset
我们稍后将需要一个用户列表，当在特征工程部分加入数据集时，创建一个新的列名churned
churn_data_summary = churn_data.select('userId')\
.withColumn('churned', F.lit(True))churn_data_summary.show(10)
Churned userIds
特征工程
在这一部分中，我们将构建一些特征/预测器来预测不稳定的用户。这通常意味着使用现有的列创建新的列，并进行转换，以便我们可以输入到 ML 管道中。
平均会话间隔时间
衡量用户参与度的一个有用功能是查看用户收听 Sparkify 平台的频率。我们会查看用户有多少个会话以及会话之间的平均时间。如果用户经常访问我们的应用程序/平台，他们可能会更加投入，因此，对于投入的用户来说，平均会话间隔时间应该更短。
userId 10 has 6 streaming sessions and mean time between sessions is 7 days.
在我们的模型中，我们只需要将sessionsFreqDay作为我们的特征，而userId是连接下面其他数据集的关键，剩余的列是不必要的，所以我们要过滤掉
listen_freq = listen_freq.select('userId', 'sessionsFreqDay')listen_freq.show(10)
we use userId as key when joining with other data. sessionsFreqDay is a new feature engineered predictor
会话计数，平均歌曲计数/会话
另外两个好的预测指标是给定用户的总会话数和每次会话收听的平均歌曲数。如果用户有更多的会话，并且每个会话听很多首歌，他们可能会更投入。
下面是另一个聚合数据集，显示了每个用户的收听活动摘要、他们收听了多少个会话以及所有会话的平均歌曲数。
user with id 10, has 6 sessions and on average he/she listened to 133.7 songs/session, whereas user with id 100006 has only 1 session and listened to 45 songs. User 100006 is engaged less.
下面是上面用于user_engagement的 pyspark 代码的 SQL 等价物，我只是在这里添加作为参考，它给出了完全相同的数据，但是使用了更容易阅读的参考。对于那些更喜欢 SQL 而不是 Dataframe 方法的人来说，他们会很高兴地看到他们几乎可以只用 SQL 来进行数据操作。
Same result is achieved with SQL
用户功能
像 *gender* 和 *level* 这样的用户特征在预测客户流失时会很有用，因为免费用户和付费用户或者女性和男性用户的行为可能会有所不同。让我们将这些特征保存到 *user_feature* 数据帧中
user_features = data.select('gender', 'userid', 'level')\
.distinct().orderBy('userId')user_features.show(10)
User Features
最后，让我们将上面创建的所有功能加入到一个joined数据框架中
churned column is the target variable and all other columns (except userId) are predictors
现在我们有 5 个流失预测列:
gender: M 或 Flevel: 付费或免费sessionCount:用户听了多少次音乐？每个会话有多首歌曲meanSongCount:对于给定的会话，每个用户平均听了多少首歌？(一次会议持续多长时间)sessionFreqDay:会话之间的平均天数是多少？更低意味着更多的参与。
在继续之前，让我们做一个理智检查，看看我们的预测是否有意义。下面是每个churned和gender组的sessionCount、meanSongCount和SessionFreqDay的平均值
joined.groupBy('churned','gender')\
.agg({'sessionCount': 'mean',
'meanSongCount': 'mean',
'sessionsFreqDay': 'mean'})\
.orderBy('churned','gender')\
.show()
根据上表，易激动用户倾向于拥有较低的sessionCount，他们在会话中听较少的音乐，但是更频繁地使用平台*(易激动用户的会话间隔为 2 - 2.5 天，非易激动用户为 4 - 4.8 天)。最后一个观察结果与我们的预期有点相反，但这仍然是一个有用的预测。一个有趣的发现是，女性用户比男性用户更喜欢这个平台，这对我们的模型也很有用*
joined.groupBy('churned','level')\
.agg({'sessionCount': 'mean',
'meanSongCount': 'mean',
'sessionsFreqDay': 'mean'})\
.orderBy('churned','level')\
.show()
当我们从免费或付费用户的角度看数据时，也有类似的观察。sessionCount和meanSongCount对于非搅动用户来说都高于预期，但是搅动用户更频繁地使用该平台。这是否意味着他们不是因为缺乏兴趣，而是因为找不到想要的音乐而烦恼，并因此离开了平台？这将是一个有趣的想法，供以后探讨。但是我们不会在这篇博文中看到它。问问题总是一个好习惯，如果有矛盾的发现，问更多的问题来发现潜在的事实。
我们有一个矛盾发现的另一个原因可能只是统计学上的怪癖。辛普森悖论在某些情况下可能会困扰我们。为了确保我们没有陷入辛普森悖论的统计谬误，让我们一起来看看 T10 和 T11。
joined.groupBy('churned','level', 'gender')\
.agg({'sessionCount': 'mean',
'meanSongCount': 'mean',
'sessionsFreqDay': 'mean'})\
.orderBy('churned','level', 'gender')\
.show()
通过level和gender的划分仍然显示，被搅动的用户有更频繁的参与。从这个角度来看，不存在辛普森悖论
Spark 中的机器学习模型只能处理数字，因此我们需要将我们的分类特征转换成数字。
我们将使用StringIndexer，将数字分配给level和gender，并使用OneHotEncoderEstimator将它们编码成SparseVector。
Converting Strings into indexes then One Hot Encoding into Sparse Vector
# Dropping gender and level (temporarily), just to fit the cell into # screen for readabilityjoined.drop('gender', 'level').show()
String Indexers created level_idx and gender_idx and OneHotEncoderEstimator created dummy columns
如果数字列中的数据是偏斜的，预测值强度会降低，因此我们还应该检查数字列sessionCount、meanSongCount和sessionsFreqDay的分布。让我们目测一下
joined_pandas = joined.select('sessionCount',
'meanSongCount',
'sessionsFreqDay')\
.toPandas()f, axes = plt.subplots(2, 3, figsize=(14, 7), sharex=False)sns.distplot( joined_pandas["sessionCount"] , color="skyblue", ax=axes[0, 0])sns.distplot( joined_pandas["meanSongCount"] , color="olive", ax=axes[0, 1])sns.distplot( joined_pandas["sessionsFreqDay"] , color="gold", ax=axes[0, 2])
All numeric features seems left skewed. We better transform them before modelling.
所有看起来都向左倾斜，让我们应用log和sqrt变换来处理倾斜
Original distribution (top row) and features after transformation. Transformation helped to fix skew.
这看起来更好，我们回到我们的 spark 数据框架joined并在那里应用转换。
joined = joined\
.withColumn('logSessionCount', F.log('sessionCount'))\
.withColumn('sqrtMeanSongCount', F.sqrt('meanSongCount'))\
.withColumn('sqrtSessionsFreqDay', F.sqrt('sessionsFreqDay'))joined.cache()
joined.take(1)
New Columns logSessionCount, sqrtMeanSongCount and sqrtSessionsFreqDay are created
最后，我们需要创建两列，label和features。Spark ML 需要label列，它必须是数字，这意味着我们需要通过转换成整数将churned列转换成label
我们还需要所有的预测器都在一个向量中，Spark ML 处理向量，所以我们需要连接所有这 5 个预测器。
Spark needs a features vector, therefore we add all our predictors into a Dense Vector using VectorAssembler
nonScaledFeatures is a DenseVector created from other features
仅关注nonScaledFeatures和label:
joined_vector.select('nonScaledFeatures', 'label')\
.show(10, truncate=False)
机器学习模型更好地处理标准化/规模化数据，因此我们需要在进入建模阶段之前标准化/规模化我们的特征
scaled = StandardScaler(inputCol="nonScaledFeatures",
outputCol="features",
withStd=True,
withMean=True)\
.fit(joined_vector)\
.transform(joined_vector)scaled.select('features', 'label').show(10, truncate=False)
建模
现在我们在向量中有了缩放/标准化的特征，并标记为预测因子，我们可以进入建模阶段。
这里，我们需要将完整的数据集分成训练集、测试集和测试几种机器学习方法。我们将评估各种模型的准确性，必要时调整参数。我们将根据测试准确性确定获胜模型，并在验证集上报告结果。
我们评估三个模型，DecisionTreeClassifier和RandomForestClassifier和GBTClassifier。RandomForestClassifer和GBTClassifier将许多DecisionTreeClassifiers结合起来，通过使分类器不容易过度拟合来获得更好的精度。DecisionTreeClassifer一般倾向于过度拟合。
我们的第一个模型是带有默认输入的DecisionTreeClassifier
train, test = scaled.randomSplit([0.8, 0.2], seed=42)tree = DecisionTreeClassifier()
prediction_tree = tree.fit(train).transform(test)
prediction_tree.groupBy('label', 'prediction').count().show()
DecisionTreeClassifier predicted 9 churns correctly, but missed 8 churns
为了更深入的了解，我们可以看看每个预测的概率值
prediction_tree.select('label', 'prediction', 'probability')\
.show(20, False)
labels, predictions and probability values for each user
流失事件相当稀疏，基本决策树分类器能够成功预测 9 个流失事件，错过 8 个，错误预测 2 个客户将流失，但没有。
现在我们来看一下带有默认超参数的RandomForestClassifier
prediction_forest = RandomForestClassifier()\
.fit(train)\
.transform(test)prediction_forest.groupBy('label', 'prediction').count().show()
RandomForestClassifier predictions
与DecisionTreeClassifier相比，RandomForestClassifier在这种情况下并没有表现得特别好
最后让我们看看梯度增强树分类器
from pyspark.ml.classification import GBTClassifiergbt = GBTClassifier(featuresCol='features', labelCol='label')
pred_gbt = gbt.fit(train).transform(test)pred_gbt.groupBy('label', 'prediction').count().show()
GBTClassifier has better prediction capability in this test set
GBTClassifier能够超越简单的DecisionTreeClassifier和RandomForestClasssifier。它能够成功地预测 11 次搅拌，只错过了 6 次搅拌。
模型评估、管道、交叉验证和网格搜索
尽管我们通过创建几个表来检查我们的预测准确性，但是 Apache Spark ML 库已经内置了评估器来评估我们的模型是否成功。在我们的例子中，我们的标签是我们可以使用的唯一一个二进制类BinaryClassificationEvaluator。
另一方面，管道对于创建健壮的 ML 管道和避免测试数据泄漏到模型训练中非常有用。管道也让我们的代码更加整洁。
交叉验证允许我们将训练数据集划分为折叠，并在每个折叠上进行训练，以找到更好的模型。
许多模型、转换器和索引器都有超参数。GridSearch 允许我们尝试许多不同的超参数，并使用这些超参数的笛卡尔乘积自动构建许多模型。找到理想的参数被称为模型调整，GridSeach 使这变得轻而易举。
让我们再次获取连接的数据，并在输入管道之前应用log或sqrt转换
joined = user_features\
.join(churn_data_summary,
on=['userId'],
how='left')\
.join(user_engagement,
on=['userId'],
how='left')\
.join(listen_freq,
on=['userId'],
how='left')\
.orderBy('userId')\
.fillna(False, subset=['churned'])\
.withColumn('churned',
col('churned').cast('integer'))\
.withColumn('logSessionCount', F.log('sessionCount'))\
.withColumn('sqrtMeanSongCount', F.sqrt('meanSongCount'))\
.withColumn('sqrtSessionsFreqDay', F.sqrt('sessionsFreqDay'))joined.cache()
joined.show(5)
cleaned dataset, ready for Pipeline
现在定义 ML 管道的每个阶段
Stages define the order of transformations to be applied
[StringIndexer_68305429c317,
StringIndexer_c98ace58d733,
OneHotEncoderEstimator_0216ca1e6769,
VectorAssembler_45b1ace90533,
MinMaxScaler_bfcb63301c7e,
StandardScaler_dd543c1f0b07,
RandomForestClassifier_dce35260ff2f]
创建一个网格搜索对象，我们定义不同的RandomForestClassifer参数，如maxDepth或featureSubsetStrategy
还要定义评估器，该评估器将评估由网格搜索创建的模型的性能，并且将选择具有基于评估器的 bets 度量的模型。
"Number of models to be tested with 5 crossfolds = {}".format(5*len(params))
> Number of models to be tested with 5 crossfolds = 60
由于我们正在安装 60 个模型，这可能需要一些时间。Spark 可以并行拟合模型以加速训练，但在这种情况下，我们串行进行，因为我们只有本地集群。
train, test = joined.randomSplit([0.8, 0.2], seed=42)
model = cv.fit(train)model.avgMetrics
Evaluator results for each model 3 x 4 hyperparameters
我们确实有很多模型，但是cv.transform() 会根据评估器自动选择最佳模型。
让我们用最好的模型来预测
predictions = model.transform(test)
predictions.groupBy('churned', 'prediction').count().show()
Prediction results with RandomForestClassifier with optimum params
在这种情况下，与普通的RandomForestClassifier相比，我们没有得到不同的结果，但是我们看到了如何使用管道、交叉验证和评估器使我们的生活变得更容易。
结论
作为数据科学家，我们总是试图用数据回答问题，并遵循数学和统计的严谨性。在这篇文章中，问题是预测用户流失，使用一种系统的方法，我们取得了一些成功。这是第一次迭代，通常模型构建需要多次迭代，通过添加更好的预测器或更多数据来进一步优化预测。
我们还看到了如何使用 Apache Spark 来处理数据、创建新功能、构建模型并对其进行评估。尽管这一转换基于一小部分数据，但同样的转换可以应用于更大的数据集，而无需更改代码。
这是一个很长的阅读，但希望你喜欢它
扫一扫
722
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
