TowardsDataScience 博客中文翻译 2019（三百四十九）

原文：TowardsDataScience Blog

协议：CC BY-NC-SA 4.0

流体力学遇上数据科学

原文：https://towardsdatascience.com/model-selection-yacht-hydrodynamics-data-set-ec0f8591e8e8?source=collection_archive---------27-----------------------

流体力学、数据科学、统计学

UCI 游艇水动力学知识库的统计分析

Yacht Hydrodynamics

一切都始于一个浪漫的故事，当流体力学遇到统计学……他们一见钟情……(呃)……

好吧，别走开，让我们进入我们的科学。本研究的主要目的是实践关于模型选择的统计分析，并将其应用于实际问题。我正在使用的数据集来自 UCI 机器学习库 **。**这个数据集在机器学习研究中非常有名，它已经存在很长时间了。您可以很容易地从这个存储库中找到许多关于数据集的研究论文。我将介绍我的近似模型，以及我用来评估不同模型的方法。我在这篇文章中只使用了统计方法，所以没有花哨的机器学习或深度学习术语。此外，感谢队友们的帮助:拉宾杜兰庞斯、辛西娅兴、韦德梁、陈子鸿(校对用)。( Github Repo 针对本项目)

在建模之前，让我们检查一下数据集:

**Left-**Dataset Information **Right-**Prismatic Coefficient of A Ship(Variable 2)

Dataset Examples

因此，我们需要使用前 6 个变量来预测 V7 的值，V7 代表每单位重量位移的剩余阻力。

Left — Import Library and Implement Forward Selection Right — V6 against V7 (exponential growth)

首先，将库导入 R studio 并使用正向选择对各个变量的影响力进行排名，你很容易发现 V6 排名第一，然后是 V2。

然后下一步，让我们绘制 V7 和 V6 的对比图。输出图显示 V6 和 V7 之间的正相关**(指数增长)**。

在我们进行回归之前，我先陈述一下我的假设:零假设是我的回归方程中每个系数都等于 0；换句话说，V1 到 V6 对我们的响应变量 V7 没有影响。

Null Hypothesis — Beta 0 to 7 all equal to 0

第一款测试车型(1 号车型)，V7 ~ V6 + V2 :

Model #1: V7 = b0 + b1V6 + b2V2

从汇总输出中，你可以发现 V2 系数的 Pr( > |t|)大约比 0.025 或 0.05 大 0.395，这意味着不显著**，所以我们不能拒绝我们的零假设。但是 V6 的系数非常重要，所以我们应该保留 V6 和从我们的预测模型中移除 V2 。此外，正常的 Q-Q 图显示数据集不完全符合我们的 Q-Q 线，残差图显示一种模式，这是不好的。也许我们需要对 y 应用对数变换，在 y 大的时候降低 y 的值。因此，我们使用 log(V7)作为新的响应变量。**

第二个模型的结果—从模型#1 推导而来:

Model #2: log(V7) = b0 + b1*V6

哇！我们在正确的轨道上！

• 残差图并不完美，但看起来比模型#1 更好
• 所有变量都很重要
• Adj R 从 0.65 到 0.9698 提升了不少

从残差图中，我有一种感觉，多项式的阶数越高，图就越好看。所以让我试试我的 3 号模型。

第三个模型的结果——从模型#2 推导而来:

Model #3: log(V7) = b0 + b1V6 + b2(V6)²

Emmmm……残差图仍然有一个模式

• 这个模型比 2 号模型好
• 添加更多术语肯定会增加 R ，考虑改用 AIC 或 Adj R

从汇总输出来看，V⁶之前的系数也很重要，所以让我尝试更高的多项式阶。

第四个模型的结果—从模型#3 推导而来:

Model #4: log(V7) = b0 + b1V6 + b2(V6)² + b3*(V6)³ + b4*(V6)⁴

爽！这是迄今为止最好的型号，波浪外观的残留剧情已经被成功取消。型号#4 的 Adj R 约为 0.9826。

结论:

Conclusion: from model #1 to model #4

从上到下我们可以发现，AIC 从最初的 2224 降到了 30。

Adj R² and AIC Improvement Calculation

关于模型#4 的一个有趣的事情，没有提供 AIC 值，Adj R 只增加了 0.9% 。有些人可能认为从模型#3 到模型#4 是不必要的，因为改进是如此之小，但增加一项会增加复杂性。然而，重要的是 4 号模型的 AIC 比 3 号模型的 AIC 小得多。

选型提示:

1. 如果每个型号的 Adj R 相似，请考虑使用 AIC 或 Cp，它们可能会有显著差异。
2. 使用正向选择或反向选择可以帮助计算出重要变量。

感谢您停留在此页面！

如果你也对**计算机架构、**数据科学、**文学、**感兴趣，欢迎你来看看！

使用 XGBoost 进行模型调整和特征工程

原文：https://towardsdatascience.com/model-tuning-feature-engineering-using-xgboost-ef819bccc82e?source=collection_archive---------18-----------------------

在这篇文章中，我将讨论数据预处理技术，这是我之前为初学者发布的系列文章 Tensorflow 的延伸。

构建人工智能或机器学习模型的第一步是预处理数据，以确保模型的正确表示。这个阶段是最关键的部分，最终决定了模型的健壮性和生产力。一个机器学习模型作为一个行业级应用的产品发布时，必须灌输一些核心功能:

1. 概括:模型应该能够跨应用程序对有噪声的数据执行优化和准确的操作，即输入数据的变化仍然应该产生有效的输出。
2. 正则化:在给定有偏数据的情况下，模型应该产生准确的输出，即行业级应用中的数据仓库在简单训练的基础上产生不稳定的模型，因为数据分布混乱，要么欠拟合，要么过拟合。
3. 可再现性:模型在两个不同的时间阶段或应用中对相同的输入产生相同的输出。

在本文中，我们将介绍两个预处理任务，即超参数优化和特征工程，这是两个最重要的模型调整方法，可以提高模型的性能，使其符合行业标准。

本文中使用的 python 库称为 XGBoost，是一个跨 Kaggle、Topcoder 和 Tunedit 等平台的常用渐变增强库。使用 XGBoost 的优势来自于这样一个事实，即该库通过创建弱分类器树来实现调整模型的正则化方法。(它与 scikit-learn 和 Tensorflow 兼容，这是目前使用最广泛的两个机器学习库。)

本文中使用的数据集是由 10 个参数组成的 pima-indian-diabetes 数据集。我们从导入必要的库开始

下一步是导入数据集并准备我们的训练和测试样本。此外，我们还创建了两个新变量，称为变量 _ 参数和静态 _ 参数，，它们将分别存储要调整的参数和保持不变的参数。

下一步是将这两个变量传递到我们的网格搜索函数中，并使模型符合我们的数据。

这种简单的实现可以扩展到更大的应用程序。fit()方法。

我们将使用的下一种方法叫做特征工程或分析。本文不包括自动化特征工程工具，如 FeatureTools 等。而是简单地使用 XGBoost 深入分析数据集特性的意义和重要性。

我们使用相同的数据集，但这次我们创建了一个名为 dtrain 的变量，它的类型为“xgboost.core.DMatrix ”,这是 xgboost 培训中使用的多维矩阵表示。

XGBoost 有一个非常好的属性叫做。get_dump()，它提供了关于创建的每棵树的详细信息，以及树在某个参数上的划分方式。我们使用这个属性来标出所有参数的重要性，并对它们进行比较。记得在你的工作目录中创建一个名为 featmap.txt 的空文本文件。

我们获得的图显示了参数-葡萄糖浓度的显著高值，这使其成为最重要的特征，因此可以进一步研究。

本文末尾提供了完整的代码。点击横幅！

要更深入地了解超参数调优是如何工作的，以及如何利用 sklearn 之外的库，请查看 neptune.ai 博客提供的这篇文章。

本文到此为止。在下一篇文章中，我将讨论数据可视化以及如何创建专业报告。干杯！

GPT-2 生成的 NFL 子编辑评论

原文：https://towardsdatascience.com/modeling-and-generating-nfl-subreddit-comments-99c8cfc35bf2?source=collection_archive---------24-----------------------

几周前，一个朋友告诉我关于 pushshift reddit 评论数据库的事情，这是一个包罗万象的 reddit 评论数据集，可以追溯到 2005 年。作为一个狂热的体育迷，我对研究不同美国运动队的子评论活动的差异有着浓厚的兴趣。

项目范围

随着 NFL 赛季的快速临近，我专注于 NFL 球队的 subreddits，因为我认为在这个激动人心的时刻，球迷们会特别感兴趣。特别是，我想制作一个交互式地图，让球迷可以随时查看情绪和话题，并查看在此期间发布的与他们球队相关的新闻故事。我还开始为每支 NFL 球队独立训练 32 个 GPT-2 模型，这些模型会产生 subreddit 风格的评论。该项目的 GPT-2 部分的目标只是看看该模型是否可以复制代表给定团队粉丝群的刻板态度和文化的幽默镜头。

数据采集

pushshift 评论数据库是一个不可思议的资源，但是每个月解压缩的 reddit 评论可能高达 100GB JSON 文件，所以我写了一个小的脚本来帮助解析每个解压缩的文件。该脚本一次下载一个月的评论，使用“grep”只保留所需子编辑的评论，将评论写入 mongoDB，删除 100GB 的 JSON 文件，然后下载下个月的 reddit 评论。如果您决定使用该数据集，希望它可以成为访问该数据集的良好起点。

有了数据集，我就可以开始主题建模、情感分析和 GPT 2 文本生成了。

情感分析

我使用了一个预先训练好的手套情绪分析器，因为它非常强大，简化了项目的情绪分析部分。GloVe 是给定语料库中所有单词的共现矩阵，它通过“最小化单词向量 I 和 j 的点积之差与 I 和 j 共现次数的对数之差的平方”(Trailhead)来训练其单词向量。参见下面的目标函数:

函数 f(Pij)对非常常见和不常见的单词进行加权，以确保非常常见的单词不会对单词向量产生不适当的影响。在将我的评论数据按团队和月份分类后，多亏了 GloVe，很容易就能获得情绪。

主题建模

在尝试了几种不同的主题建模技术之后，包括 LDA、LSA 和 NMF。最后，我选择了 LDA，因为它始终如一地提供最好的话题，通过了测试。在将一长串停用词传递给模型后，该模型持续产生与给定时间段内 NFL 中值得注意的事件一致的主题。

情感和主题地图

对于给定的一个月，球迷可以悬停在他们球队体育场的位置上，并查看该月的话题和情绪以及该月发表的 ESPN NFL 文章(见上文)。在 flask 应用程序中生成如上图所示的 plotly 地图的代码:

import plotly.graph_objs as go
import plotly
import pandas as pd
import jsondef create_map(df):
    data = [
        go.Scattergeo(
            lon = df['long'],
            lat = df['lat'],           
            mode="markers" ,           
            locationmode="USA-states",           
            hovertext = df['text']
        )
    ]
    graphJSON = json.dumps(data,cls=plotly.utils.PlotlyJSONEncoder)
    return graphJSON

中间步骤是在 flask route 函数中将 JSON 对象传递给 html，使用:

return flask.render_template("map.html",plot=graphJSON)

GPT 新协议评论生成

Max Woolf 为 OpenAI 的 GPT-2 的 Tensorflow 实现制作了一个真正好的包装器。他的软件包可以选择使用 117、124 或 355 百万参数模型。这个包使得定制您自己的语言生成模型变得非常容易:

import gpt_2_simple as gpt2model_name = "124M"
gpt2.download_gpt2(model_name=model_name)   # model is saved into current directory under /models/124M/

sess = gpt2.start_tf_sess()
gpt2.finetune(sess,
              'corpus.file_extension',
              model_name=model_name,
              steps=1000)   # steps is max number of training steps

gpt2.generate(sess)

注意，无论您在哪个目录中训练模块，都会有一个新形成的 checkpoint/run1 文件夹。如果您想要定制训练几个模型，将参数 run_name 插入到 gpt2.finetune()步骤中，如下所示:

import gpt_2_simple as gpt2sess = gpt2.start_tf_sess()
gpt2.finetune(sess,
             'corpus.file_extension',
              model_name=model_name, 
              run_name="your_custom_run_name",             
              steps=1000)

在设置此变量并训练模型之后，权重/模型将位于名为/check point/your _ custom _ run _ name/的新目录中。如果您将一列 CSV 传递给 finetune 函数，gpt-2-simple 模块将为您预处理文本，并且几乎不需要准备数据。如果您对自定义语料库的文本生成感兴趣，我强烈推荐使用 gpt-2-simple。

网站

我用 flask 做了这个网站，并在 AWS 的 uWSGI 服务器上托管。当有人选择“看看机器人在说什么”时，他们可以输入他们团队的名字，并设置温度，以控制波尔兹曼分布的随机性。温度越高，产生的输出越随机，越不像训练集。

用户可以输入一个团队的名称(只能做 22 个左右，因为我没有在每个团队中训练)，并设置所需的温度(建议 0.7-0.9)。请随意点击查看网站。

结论

这个项目在主题建模、情感分析和文本生成方面是一个很好的练习。我训练的一些初始团队进行了更多的迭代，所以如果你最终尝试了这个应用程序，这个模型对 Jets，Bears 和其他几个团队有更好的输出。gpt-2-simple 包是一个非常有趣的包，希望上面的启动代码对你有所帮助。感谢阅读！

来源

Trailhead，https://trail head . sales force . com/en/content/learn/modules/glove-and-word-vectors-for-sensation-analysis/use-glove-for-natural-language-processing

GPT-2-简单，【https://github.com/minimaxir/gpt-2-simple

项目源代码，【https://github.com/sethweiland/reddit_nfl_comments

BiDAF 中的建模和输出层——带 Minions 的图解指南

原文：https://towardsdatascience.com/modeling-and-output-layers-in-bidaf-an-illustrated-guide-with-minions-f2e101a10d83?source=collection_archive---------14-----------------------

BiDAF 的权威指南——第 4 部分，共 4 部分

BiDAF 是一个流行的机器学习模型，用于问答任务。本文借助一些可爱的小喽啰来解释 BiDAF 的建模层。

Source: Unsplash

他的文章是一系列四篇文章中的最后一篇，旨在说明 【双向注意力流】一种流行的问答机器学习模型(Q & A)的工作原理。

概括地说，BiDAF 是一个封闭领域，抽取的 Q & A 模型。这意味着为了能够回答一个查询**，BiDAF 需要查阅一个附带的文本，该文本包含回答查询所需的信息。这个附带的文本被称为上下文。 BiDAF 的工作方式是提取上下文中最能回答查询的子串——这就是我们所说的对查询的回答 。我有意将“查询”、“上下文”和“回答”这几个词大写，以表明我是在用它们的专业技术能力来使用它们。**

An example of Context, Query and Answer

对前面三篇文章的简要总结如下:

• 本系列的第 1 部分提供了 BiDAF 的高级概述。
• 第 2 部分解释了 BiDAF 如何使用 3 种嵌入算法来获得上下文和查询中单词的向量表示。
• 第 3 部分探讨了 BiDAF 的注意力机制，它结合了来自上下文和查询的信息。

前述注意步骤的输出是一个名为 G 的巨型矩阵。 G 是一个 8d -by- T 矩阵，它对上下文单词的查询感知表示进行编码。G 是建模层的输入，这将是本文的重点。

我们一直在做什么？G 实际上代表什么？喽啰来救援了！

好的，我知道在过去的三篇文章中我们已经经历了很多步骤。人们很容易迷失在无数的符号和方程式中，尤其是考虑到 BiDAF 论文中的符号选择并不是那么“用户友好”我的意思是，你还记得 H，u，ĥ和ũ分别代表什么吗？

因此，现在让我们后退一步，试着得到我们到目前为止所做的所有这些矩阵运算背后的直觉。

实际上，前面所有的步骤都可以分解为两个步骤集合: 嵌入步骤 和 注意步骤 。如上所述，所有这些步骤的结果是一个名为 G 的 8d -by- T 矩阵。

下面可以看到一个 G 的例子。 G 的每一列都是一个单词在上下文中的8dby-1向量表示。

An example of G. The length of the matrix, T, equals the number of words in the Context (9 in this example). Its height is 8d; d is a number that we preset in the word embedding and character embedding steps.

现在让我们玩一个小游戏，(希望)能帮助你理解之前文章中所有的数学术语。具体来说，让我们把上下文中的单词想象成一群有序的奴才。

Think of our Context as a bunch of Minions, with each Context word corresponding to one Minion

我们的每个小喽啰都有一个大脑，可以储存一些信息。现在，我们的奴才们的大脑已经很拥挤了。每个 Minion 的当前脑内容相当于 Minion 所代表的上下文单词的 8d 乘 1 列向量。在这里，我展示了“新加坡”奴才的大脑扫描图:

奴才们的脑子可不是一直这么满的！事实上，当它们出现时，它们的大脑几乎是空的。现在让我们回到过去，想想那些奴才们为了获得他们当前的知识状态而经历了什么样的“课程”。

小黄人上的前两节课是“单词嵌入”和“字符嵌入”在这些课程中，奴才们了解了自己的身份。“单词嵌入”课的老师 GloVe 教授教小喽啰们关于他们身份的基本信息。另一方面，“角色嵌入”课程是一门解剖学课程，在这门课程中，小喽啰们通过反复扫描来了解自己的身体结构。

Minions in the “Character Embedding” class

这是这两节课后“新加坡”奴才的脑部扫描。

The “Singapore” Minion understands his identity after attending the “Word Embedding” and “Character Embedding” lessons

紧接着，小喽啰们继续前进，参加了“情境嵌入”课程。本课是一堂会话课，在课上，小喽啰们必须通过一款名为 bi-LS 的 messenger 应用程序相互交流。bi-LS 协助的 convo 允许奴才们了解彼此的身份——这是他们在前两节课中学到的。很整洁，是吧？****

Two Minions having a fun conversation through bi-LS™, sharing information about themselves. Source: Giphy

我在“情境嵌入”课结束后，对“新加坡”的小跟班进行了另一次核磁共振扫描。如你所见，现在我们的小家伙知道更多的东西了！

Now “Singapore” knows both his and his neighbors’ identities!

我们的奴才们正在愉快地学习，突然一个人闯进了他们的学校😱原来他的名字是提问先生，他是一名记者。他在这里:

The inquisitive Mr. Query. He has an urgent question —”Where is Singapore situated” — and he knows some of our Minions hold relevant information for this question.**

Query 先生急需为他正在写的一篇文章收集一些信息。具体来说，他想知道"新加坡位于哪里【Query 先生知道我们的一些奴才在大脑中保存着这些信息。**

我们的奴才们虽然乐于助人，但还是想帮 Query 先生一把。为此，他们需要选择几名团队成员与 Query 先生会面，并提供他所寻求的信息。这一群拥有查询先生相关信息并将被派遣给他的爪牙被称为回答帮。

The Answer Gang, which collectively holds the answer to Mr. Query’s question. Only relevant Minions can join the Answer Gang!

现在，我们的跟班们有一项任务要做——他们需要集体决定谁应该、谁不应该加入回答帮。**他们这样做时需要小心！如果他们从答案帮中漏掉了太多本应包括在内的爪牙，查询先生就不会得到他需要的所有信息。这种情况被称为低回忆，Query 先生对此深恶痛绝。

另一方面，如果太多不必要的爪牙加入回答帮，查询先生就会被多余的信息淹没。他称这种情况为低精度，他也不喜欢这样！众所周知，Query 先生有一些愤怒管理问题👺所以给查询先生提供适量的信息对我们的下属最有利。

那么喽啰们怎么知道他们中的哪一个应该加入答题帮呢？

解决这个问题的方法是组织几次会面，统称为“关注”在这些会议中，每个仆人都可以单独与查询先生交谈，了解他的需求。换句话说，注意力会议允许下属衡量他们对提问者的问题的重要性。

这是“新加坡”奴才在注意力会议上闲逛时大脑的核磁共振扫描。这相当于我展示的第一张大脑扫描图像。

Singapore’s current brain content. He knows quite a bit— but he is still missing one thing!

正如我们所看到的，我们的奴才们的大脑现在已经很满了。以奴才们目前的知识水平，他们现在可以开始选择答案帮的成员了吗？不，不完全是！他们仍然缺少一条关键信息。我们的每个爪牙都知道他对查询先生的重要性。然而，在做出这个重要的决定之前，他们还需要意识到对方对 Query 先生的重要性。

正如你可能已经猜到的，这意味着奴才们必须第二次互相交谈！现在你知道这个对话是通过 bi-LS 应用完成的。

The Minions during the modeling step meeting. Here, they talk to each other through bi-LS™ and share their relative importance to Mr. Query. Source: Free PNG Logo

****这种由 bi-LS 促成的对话也被称为“建模步骤”,是我们当前文章的重点。现在让我们详细了解这一步骤！

第十步。建模层

好吧，让我们暂时离开我们的奴才，回到符号和方程上来，好吗？没有那么复杂，我保证！**

**建模层相对简单。它由两层双 LSTM 组成。如上所述，建模层的输入是**G。第一个双 LSTM 层将 G 转换为 2d -by- T 矩阵，称为 M1 。

M1 然后作为第二个双 LSTM 层的输入，第二个双层将其转换为另一个 2d -by- T 矩阵，称为 M2 。

从 G 到 M1 和 M2 的形成如下图所示。

Step 10. In the modeling layer, G is passed through two bi-LSTM layers to form M1 and M2

M1 和 M2 是上下文单词的另一种矩阵表示。 M1 和 M2 与之前的语境词表示的区别在于 M1 和 M2 在其中嵌入了关于整个语境段落以及查询的信息。

用奴才的话来说，这意味着我们的奴才现在拥有了所有他们需要的信息来决定谁应该加入回答帮。

The “Singapore” guy now has all he needs to decide if he should join the Answer Gang.

第十一步。输出层

好了，现在我们到了终曲！只要再走一步，我们就完成了！

对于上下文中的每个单词，我们都有两个数值向量来编码单词与查询的相关性。太棒了。我们需要做的最后一件事是将这些数字向量转换成两个概率值，这样我们就可以比较所有上下文单词的查询相关性。这正是输出层所做的。

在输出层中， M1 和 M2 先与 G 纵向串接形成G； M1 和克； M2 。两者都[G； M1 和克； M2 的尺寸为 10d -by- T 。

然后，我们通过以下步骤获得 p1 ，即整个上下文中开始索引的概率分布:

类似地，我们通过以下步骤获得最终指数的概率分布 p2 :

获取 p1 和 p2 的步骤如下图所示:

Step 11. The output layer, which converts M1 and M2 to two vector of probabilities, p1 and p2.

p1 和 p2 然后被用来寻找最佳答案跨度。最佳答案跨度只是具有最高跨度分数的上下文的子串。依次，跨度分数是该跨度中第一个单词的 p1 分数和该跨度中最后一个单词的 p2 分数的简单乘积。然后，我们返回具有最高 span 分数的 span 作为我们的答案。**

举个例子就能说明这一点。如您所知，我们目前正在处理以下查询/上下文对:

• 语境 : “新加坡是位于东南亚的一个小国。” ( T = 9)
• 提问 : “新加坡位于哪里？” ( J = 4)

在通过 BiDAF 运行这个查询/上下文对之后，我们获得了两个概率向量— p1 和 p2。

上下文中的每个单词都与一个 p1 值和一个 p2 值相关联。 p1 值表示单词成为答案范围开始单词的概率。以下是我们示例中的 p1 值:

我们看到模型认为我们的答案跨度最有可能的起始词是“东南”

p2 值表示单词成为答案范围的最后一个单词的概率。以下是我们示例中的 p2 值:

我们看到我们的模型非常确定，几乎 100%确定，我们的答案跨度最有可能的结尾词是“亚洲”

如果在原始上下文中，具有最高 p1 的单词在具有最高 p2 的单词之前出现，那么我们已经有了我们的最佳答案范围——它将简单地以前者开始，以后者结束。这就是我们例子中的情况。因此，模型返回的答案将简单地是“东南亚”**

就这样，女士们先生们——经过 11 个漫长的步骤，我们终于得到了我们问题的答案！

Here is Mr. Query with “Southeast” and “Asia”, both of whom have been selected to join the Answer Gang. It turns out that the information provided by “Southeast” and “Asia” is just what Mr. Query needs! Mr. Query is happy🎊

好吧，在我们结束这个系列之前有一个警告。假设具有最高 p1 的上下文单词在具有最高 p2 的上下文单词之后出现，我们还有一点工作要做。在这种情况下，我们需要生成所有可能的答案范围，并计算每个答案的范围分数。以下是我们的查询/上下文对的可能答案范围的一些示例:****

• ****可能答案跨度:“新加坡”；跨度分数 : 0.0000031
• ****可能答案跨度:“新加波是”；跨度分数 : 0.00000006
• ****可能答案跨度:“新加坡是 a”；跨度分数 : 0.0000000026
• ****可能答案跨度:“新加坡是个小”；跨度分数 : 0.0000316

然后，我们将具有最高跨度分数的跨度作为我们的答案。

这就是它 BiDAF 中每一步的详细说明，从开始到结束(撒上健康剂量的 Minion-joy)。我希望这个系列能帮助你理解这个迷人的 NLP 模型！

Source: Giphy

如果你对这篇文章有任何问题/评论，或者想联系我，请随时通过 LinkedIn 或 gmail DOT com 发送电子邮件至 meraldo.antonio。

词汇表

• ****上下文:包含查询答案的查询附带文本。
• ****查询:模型应该给出答案的问题。
• ****答案:包含可以回答查询的信息的上下文的子字符串。这个子串将由模型提取出来。
• ****跨度得分:一个答案跨度中第一个单词的 p1 值与该答案跨度中最后一个单词的 p2 值的乘积。
• T :上下文中的字数。
• J :查询的字数。
• G :一个大的、 8d -by- T 矩阵，包含查询感知的上下文表示。G 是建模层的输入。
• M1: 将 G 通过双 LSTM 得到的一个 2d -by- T 矩阵。 M1 包含上下文单词的向量表示，其具有关于整个上下文段落以及查询的信息。
• ****M2: 一个 2d -by- T 矩阵，通过一个 bi-LSTM 传递 M1 得到。 M2 ，就像 M1 一样，包含上下文单词的向量表示，其具有关于整个上下文段落以及查询的信息。
• p1: 一个长度为 T 的概率向量。每个上下文单词都有自己的 p1 值。该 p1 值表示该单词成为答案区间中第一个单词的概率。
• p2: 一个长度为 T 的概率向量。每个上下文单词都有自己的 p2 值。该 p2 值表示该单词成为答案区间中最后一个单词单词的概率。

参考

[1] 机器理解的双向注意力流(闵俊 Seo et。艾尔，2017)

如果你对这篇文章有任何意见或者想联系我，请随时通过 LinkedIn 给我发一个联系方式。此外，如果你能通过我的推荐链接支持我成为一名中级会员，我将非常感激。作为一名会员，你可以阅读我所有关于数据科学和个人发展的文章，并可以完全访问所有媒体上的故事。

累积影响建模—第一部分

原文：https://towardsdatascience.com/modeling-cumulative-impact-part-i-f7ef490ed5e3?source=collection_archive---------8-----------------------

创建累积影响的简单特征，使用健身疲劳模型预测运动表现

“一点一点，一点变成很多。”——坦桑尼亚谚语

欢迎来到建模累积影响，这是一个通过各种建模镜头查看运动训练对运动表现的累积影响的系列。旅程从第一部分开始，从生理学的一个模型“体能疲劳模型”开始，它使用指数衰减的训练强度的卷积。从第二部分开始，该系列将向更通用的技术过渡，包括样条、卡尔曼滤波器和 LSTM 驱动的神经网络。

纵向数据集，尤其是那些具有丰富事件历史的数据集，在公共领域很难找到。本文介绍了一个模拟链球运动员，该运动员模仿了 Busso、Candau 和 Lacour [1]在 1994 年发表的一篇论文中描述的一个真实的对手，在此称为“BCL94”我们称之为“H.T .”的掷链球运动员在研究时 20 岁，有 7 年的掷链球经验。这个模拟的所有代码都是 Github 上的。

The modern sport of hammer throwing. By US Army. Photo credit: Phil Sussman. — Fort Lewis Soldier…, Public Domain, Wikimedia

对于真实的和模拟的 H.T .，训练强度和表现数据是在比赛前 37 周(259 天)内收集的，包括两个一般阶段:

1. “积累”阶段持续到第 147 天，
2. 平均强度较低的“竞争”期持续到第 259 天，即“主要竞争”日

Simulated daily training intensities for our hammer thrower under a consistent weekly pattern of training intensity and different averages for the two periods.

作者的一个目标是使用训练强度的每日系列来模拟 H.T .的表现， w_n ， n = 1，…，259。(参见 BCL94 对强度的定义，强度是一个基于判断的公式，它将托举次数、练习投掷次数和跳跃次数合并为一个值。)这些强度在准备期平均为 34 个单位，在比赛期平均为 24 个单位，在最紧张的几周接近 50 个单位，至少低于 20 个单位。模拟大致模拟了这一点。

和训练强度一样， w_n ，H.T .每天测得的成绩值， p_n ，是艺术和科学的混合体。它们基于人类跑步表现的理论模型，该模型强化了一种直觉，即接近世界纪录的进步大于接近更普通水平的收获。就规模而言， p_n = 1000 相当于世界纪录，而 p_n = 0 相当于“任何身体健全的健康人都应该能够[做]”的基线表现

用来解释 H.T .在 BCL94 中表现的模型是“健身疲劳模型”，也以生理学家 Eric Banister 命名为“ Banister 模型”。在体能-疲劳模型中，一次运动训练引起两种对立效应，体能和疲劳、、T19、同时以相反的方向影响表现。乍一看，该模型表现为对两个对立协变量的线性回归:

e(p)*=496+0 . 069**体能 - 0.27 *疲劳。

注意疲劳对表现的负面影响几乎是健身正面影响的四倍。

到底什么是健身和*疲劳？*体能疲劳模型将两者表征为以下形式的指数衰减量:

g(t) = k * exp(-t / τ)，

其中 k 是初始量，τ是“衰减时间常数”(虽然健身疲劳是一个特定领域的模型，但读者可能会意识到由指数衰减近似的其他现象，如放射性衰减、记忆丧失，甚至啤酒泡沫衰退。)两者的关键区别在于，对于适应度，衰减时间常数要大得多；在 H.T .的情况下，适应性的τ_ 1 = 60相对于疲劳的τ_ 2 = 13。这导致脉冲响应的显著差异:

Under typical fitness-fatigue configurations, fatigue decays much more quickly than fitness. The impulse from fatigue at any time, however, is typically much greater than the impulse from fitness.

这些指数模型描述了初始水平的体能和疲劳的衰减，但 H.T .的体能和疲劳水平是基于过去许多训练课程的效果。应用叠加原理(适用于所有线性系统)，对一组刺激的总反应是对单个刺激反应的总和。因此，H.T .在任何时间点的体能和疲劳由所有过去训练期的时间递减效应的总和来定义。

该操作通过训练历史与指数衰减函数的离散卷积进行数学处理:

有无限数量的这样的特征可以用单个训练历史来构建，首先通过在上面的等式中调制τ，然后从除了纯指数衰减之外的模型中构建。

基于 BCL94 中给出的估计，H.T .的链球投掷成绩的个人体能-疲劳模型为:

Though it’s possible to code in a compact fashion, computing the fitness and fatigue features involves inner and outer loops through the data.

这导致了 h . t .疲劳、健康和相对(基线移除)性能的以下轨迹:

The plot above shows H.T.’s modeled antagonistic fitness (green) and fatigue (red) levels over time, as well as H.T.’s modeled performance sans the offset (i.e., 0.07 * fitness — 0.27 * fatigue).

最终的模拟性能还包括少量伪随机高斯噪声，以创建更真实的估计场景。

估计

为了估计体能疲劳模型中的参数，BCL94 首先创建时间常数τ_1 和τ_2 的值的网格。将网格中每个点的时间常数视为已知，通过时间计算健康度和疲劳度。用普通多元回归拟合 E(p _ n)=μ+k _ 1** fitness _ n**-k _ 2* fatigue _ n*，其中 n 随时间递增，记录每一步的残差平方和(RSS)。具有最小 RSS 的时间常数对与它们相关的回归估计相结合，完成了估计。

或者，也可以通过残差平方和的非线性最小化来直接恢复参数:

The fitness-fatigue parameters can be recovered using non-linear minimization of mean squared error, but it required some tweaking in R.

R 的通用优化器 optim 提供的“BFGS”方法在这个例子中运行良好，但是选择合理的初始值仍然是必要的。由于只有两个时变参数和每层的封闭解，BCL94 使用的“网格搜索”是一种有吸引力的估计方法。下面是模拟的实际性能(包括高斯噪声)和拟合模型的预测的全貌:

The points are the actual (simulated) performance measurements from H.T.; the blue line is the fitness-fatigue model’s predictions through time.

讨论和后续步骤

本文研究了一个运动表现模型，该模型利用两个对立的、基于卷积的特征来编码训练的效果。这两个特征都是指数衰减的训练强度的离散卷积，只是时间常数不同。该模型在整套参数中是非线性的，但当时间常数是固定的并被认为是已知的时是线性的，这使得网格搜索成为一种可行的估计方法。

虽然体能疲劳模型是一个来自生理学的领域特定模型，但它提出了一个更通用的特征工程策略:当存在一系列具有持续影响但也会随时间衰减的事件时，从具有指数衰减的卷积中构建特征。尝试不同的时间常数，或使用网格搜索交叉验证来探索不同的可能性。

本系列的下一部分将不再假设理论模型中的指数衰减，而是研究如何使用样条来模拟一般的衰减曲线。当手头的理论有限，或者你只是不信任存在的模型时，它会更有用。由于这位作者确实相信指数记忆衰退模型，请继续关注第二部分。

参考

[1] T. Busso，R. Candau 和 J. Lacour，疲劳和健康模拟训练对表现的影响 (1994)，欧洲应用生理学杂志

[2] R. Morton，R. Fitz-Clarke，E. Banister，模拟人类在跑步中的表现 (1985)，应用生理学杂志

累积影响建模—第二部分

原文：https://towardsdatascience.com/modeling-cumulative-impact-part-ii-2bf65db3bb98?source=collection_archive---------28-----------------------

超越指数衰减-用样条函数模拟累积影响

Welcome back, “H.T.” The sport of hammer throwing. By US Army. Photo credit: Phil Sussman. — Fort Lewis Soldier…, Public Domain, Wikimedia

在累积影响建模—第一部分中，体能疲劳模型通过期望值方程解释了我们模拟的链球运动员“H.T .”的平均表现:

在上面的预期中，“performance”p _ n经过校准，以使 0 到 1000 包括健全人到世界纪录[1]，而 w_i 是 H.T .的 performance、Busso、Candau 和 Lacour 的研究人员定义的强度加权训练量[2]。这个等式在概念上可以写成:

*预期成绩= 496 + .07 *体能-.27 疲劳，

其中体能和疲劳是训练强度 w_i 的卷积，呈指数衰减。第一部分对此进行了详细解释，但代码块 1 中的 R 代码完成了性能模拟，足以赶上:

Code Block 1. Simulating H.T.’s training and performance data using the fitness-fatigue model from Physiology.

另外，一个 Jupyter 笔记本有第一部分和第二部分的代码。为简明起见，用于绘图的 ggplot2 代码从代码块中省略，但可在笔记本或本要点中找到。

虽然体能疲劳模型在解释现实世界中许多运动员的表现方面做得很好，但我们怎么知道数据总是支持双重对抗指数呢？为了回答这个问题,《累积影响建模》的第二部分开发了一个灵活的基于样条的培训和绩效模型。为了进行测试，我们将使用它来恢复用于模拟 H.T .数据的理论模型，在下一部分中，我们将在一些真实的运动员身上进行尝试。

对于基函数 g_j(t)，t = 2，…， p ，参数 B 样条表示为η(t)=θ_ 1+σθ_ j g _ j(t)。使用 R 的样条包中的ns()函数，基函数完全由结位置的选择决定。给定一组合理的结，如果我们可以从数据中估计θ_1，…θ_j，我们将有一个具有灵活脉冲响应的新模型。

由于卷积在加法上是分布式的，所以用η(t)替换适应度和疲劳构造就可以了，最后我们得到了 p 基于卷积的变量，而不是性能回归中的两个变量。请参见下面的代数嵌入:

上面数学计算的结果是，我们不是用指数衰减来卷积训练强度，而是用样条基函数和常数偏移来卷积它们。这些新的卷积变量是:

而这些就像第一部分中的体能和疲劳一样进入了性能回归公式

为了获得我们的基函数 g_j，我们将不得不面对与参数样条相关的节点选择问题。纽结是两个近似三次多项式相遇的“断点”，应该选择它们来捕捉重要的非线性行为。但是，在一开始不知道要估计的函数的情况下，如何知道重要的非线性行为呢？

为了用鸡开始这个先有鸡还是先有蛋的问题，我们将为我们的模拟数据使用一个近乎理想的结放置。在健康-疲劳模型中，组合脉冲响应是健康和疲劳指数的线性组合:

用训练强度卷积这个单个函数与卷积两个独立的指数是一样的，我们的最终目标是用我们的样条方法从数据中估计ϕ(t。该函数反映了同时考虑健康和疲劳的训练事件的真实脉冲响应。

Figure 1. Theoretical system response (solid line) to one of H.T.’s training events. With three carefully chosen knots, our B-spline (dashed) is capable of approximating the function with minimal bias.

考察ϕ(t 的图表(图 1)与 H.T .的表现之间的关系，我们看到在培训活动后的最初几天，对 H.T .的表现有净负面影响。在训练活动后的第 25 天左右，净负面影响已经被中和。在比赛后的第 40 天左右，随着疲劳在体能训练中消散，H.T .仍能从最初的训练中获得最大收益。在第 40 天之后，体能指数缓慢下降，在大约 200 天后，一次锻炼的益处逐渐为零。

在这篇文章中，我们没有直接看到ϕ(t)，但将通过我们的 b 样条η(t)以灵活的方式估计它，它本身取决于估计的参数。一个只有三个内部节点的 B 样条曲线可以确定近似值，如果你把节点放在正确的位置，比如 14 天，40 天和 100 天，但是即使这样也会有偏差。图 1 还显示了η*(t)，或者我们的样条曲线通过线性回归直接拟合真实函数ϕ(t，这非常接近。代码块 2 中的回归代码如下所示:

Code Block 2. Comparing the true theoretical impulse response with B-spline approximation. Full example with plots shown in the Jupyter notebook.

除了直观地向我们展示η(t)与ϕ(t 的接近程度，η*(t)回归为我们提供了逼近ϕ(t).的最佳参数值θ_1，…，θ_p 如果我们可以使用训练和性能数据接近这些系数值，我们就在正确的轨道上。

print(eta_star_lm)Call:
lm(formula = plot_df$level ~ my_spline)Coefficients:
(Intercept)   my_spline1   my_spline2   my_spline3   my_spline4  
    -0.1646       0.1909       0.1014       0.3221       0.0756

是时候从实际训练和表现数据中用我们的 b 样条η(t)来估计ϕ(t 了。在这篇文章中，我们将带来另一件东西:我们完美选择的结。这是为了说明的简单性，考虑到我们对模拟世界的了解，这是“欺骗”。在本系列的下一部分中，我们将使用真正的游泳运动员的表现数据，而我们将没有这种奢侈。

回想一下第一部分中，通过对抗性适应度和疲劳结构来估算ϕ(t 是一项数值工作，需要直接非线性估算或 OLS 回归后的时间常数网格搜索。样条卷积方法提供的灵活性的另一个好处是它完全是线性的。

以下 20 行代码仅依赖于train_df中的模拟数据，完成基于样条的卷积估计:

Code Block 3: Introducing a spline-based method for modeling cumulative impact.

回归系数估计与上述η*(t)回归的估计相匹配，至少达到估计误差。请注意，变量z_1实际上对应于η*(t)回归的截距，因为该变量是通过一列 1 与训练强度的卷积创建的。新的截距只是性能基线，在概念上与体能疲劳模型相同。

> summary(spline_reg)Call:
lm(formula = perf ~ z_1 + z_2 + z_3 + z_4 + z_5, data = train_aug_df)Residuals:
     Min       1Q   Median       3Q      Max
-16.6576  -4.3120  -0.0274   4.1162  18.6338Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) 491.668936   2.947533  166.81   <2e-16 ***
z_1          -0.153090   0.012525  -12.22   <2e-16 ***
z_2           0.177312   0.011952   14.84   <2e-16 ***
z_3           0.096656   0.007135   13.55   <2e-16 ***
z_4           0.301346   0.026191   11.51   <2e-16 ***
z_5           0.066427   0.005505   12.07   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 6.965 on 253 degrees of freedom
Multiple R-squared:  0.9408,    Adjusted R-squared:  0.9396
F-statistic: 803.5 on 5 and 253 DF,  p-value: < 2.2e-16 

出于预测的目的，不需要重建ϕ(t 的估计)并重新与数据进行交互。基于卷积的特征已经在回归中，预测是存储在train_aug_df$perf_hat中的拟合值。从建模和预测的角度来看，我们用经验驱动的替代方法取代理论驱动的体能疲劳模型的工作已经完成。就像在第一部分中一样，我们根据过去的训练课程和重建总脉冲响应函数的方法进行性能预测。

为了完全关闭循环，代码块 4 执行以下操作:

• 使用上面获得的回归系数从数据中估计真实脉冲响应ϕ(t，即“η(t) hat”
• 用拟合的脉冲响应η(t)进行卷积训练，以证明除了端部效应，回归的拟合值是等价的。
• 卷积训练与真实脉冲响应ϕ(t)，即本文中的体能疲劳模型，以作比较。

Code Block 4. Closing the loop.

请注意，get_eta_hat函数丢弃了第一个系数，即性能基线，并将其添加回第 23 行。第二个系数也是截距，η(t)中的常数项，解释了附加到样条变量的一列。

上述计算导致遵循基于样条的方法的对比图和用于模拟数据的健身疲劳模型。

Figure 2. As in Figure 1, the solid line is still the theoretical system response to one of H.T.’s training events. The dashed line also looks similar, but here it is eta_hat, the spline fit from the training and performance data.

Figure 3. Performance predictions based on past training sessions for classical fitness-fatigue and presented spline-based method.

虽然指数衰减广泛适用，但还有许多其他有趣的衰减曲线。例如，时间序列中的移动平均模型认为，随机冲击在完全消失之前会影响固定数量的时间段。或者，也许一个事件的影响会衰减到非零水平，就像我们愿意相信的那样。本文展示的基于样条的方法适应了指数衰减的这些偏差。

这篇文章是累积影响建模的第二部分，仍然依赖于模拟数据，并能够根据对世界不切实际的了解来选择结。下一部分将在真实的运动员身上研究这种方法，其中存在不同时间尺度上发生的缺失数据和动态。

虽然该方法失去了作为模型构造的健康和疲劳的直观吸引力，但封闭形式的解决方案和对初始值的缺乏使其成为一种有吸引力的替代方案。该方法做出的性能预测非常类似于其理论驱动的、完全参数化的对应方法，并且提供了脉冲响应的估计，该估计继续具有直观的主题解释。

在下一篇文章中，我们“放下锤子”和五个真正的游泳者一起下水，在实践中测试这个方法。它将出现在累积影响建模的第三部分。

参考

[1] R. Morton，R. Fitz-Clarke，E. Banister，模拟人类在跑步中的表现 (1985 年)，《应用生理学杂志》

[2] T. Busso、R. Candau 和 J. Lacour，疲劳和健康，模拟训练对表现的影响 (1994 年)，《欧洲应用生理学杂志》

累积影响建模第三部分

原文：https://towardsdatascience.com/modeling-cumulative-impact-part-iii-1b216273b499?source=collection_archive---------12-----------------------

根据过去的训练强度，模拟五名优秀游泳运动员的表现

A swimmer performing freestyle, from Wikipedia

到目前为止，本系列已经研究了来自生理模型的模拟数据，即运动表现的“健康-疲劳”模型。虽然模型是特定的，但获得的关于累积影响的见解是普遍适用的。回想一下第一部分展示了指数衰减卷积的使用，而第二部分使用样条函数概括了卷积函数。

虽然模拟数据是理解方法的优秀工具，但它为我们提供了一些不自然的优势(例如，在选择参数初始值和样条节点放置方面)。本文将第一部分和第二部分中的方法应用于一个真实的数据集，涉及转换、试错和一点直觉。我们还将计算标准误差，因为不再可能“看到值已经恢复”

这些数据由克里斯蒂安·拉舍(约翰内斯·古腾堡大学美因茨**)体育科学研究所 Mark Pfeiffer 教授的工作组)提供，并由 Kolossa 等人[1]研究，这些数据来自一项持续 160 天的实验，其中五名游泳运动员在水中和陆地上进行训练，并进行每周一次的游泳速度测试。在水中，训练被量化为游泳强度(1、2、3、5 或 8)乘以每个强度下游泳的距离。在陆地上，1 小时的训练被转换为 2 公里的游泳距离，并根据训练类型进行加权(耐力 2，体能 5，力量 8)。水和陆地训练都累积到每日水平，并且该剂量在数据集中的范围从恰好零(即休息日)到 48.6 单位。**

原始成绩是游泳测试中三次重复的平均速度(m/s ),每次重复时增加阻力，游泳者必须在不爬墙的情况下游完 20 米。与[2]中的链球投掷器(在本系列的前两部分中研究过)不同，这里没有日常测试，因此在响应列中有“缺失值”。下面是游泳者 1 号的第一行数据，对它们的原始形式稍加处理。

day training_dose mean_velocity
1   16.02         1.241971153
2   19.31268      NA
3   9.03          1.250052057
4   12.91618      NA
5   0             NA
6   9.6657        NA
7   5.87          1.329180349

什么是绩效指标？投掷距离和游泳速度都是“越大越好”类型的变量，并且以它们的原始形式工作。然而，第一部分和第二部分使用了一个被称为“标准点”的转换后的性能指标对于标准点指标，0 代表身体健康的人的表现，1000 代表世界纪录的表现，接近人类极限的进步比更普通水平的进步更有价值。

克里斯蒂安·拉舍(提供了数据)提到，在当前背景下，使用标准点转换是“有争议的”。注意到他的谨慎，我们将在这里使用它。第一个原因是为了实现第一部分和第二部分的连续性，第一部分和第二部分使用了链球运动的标准点。第二个原因是，转换本身就是一个主题，是任何数据科学工具箱的一个有价值的补充。

在 1990 年的文章模拟人类在跑步中的表现 [3]中，一种量化表现测量的方法基于这样的观察，即随着时间的推移，世界记录近似遵循形式 y = L + a exp *(-x / b)，*以由时间常数 b 确定的速率接近渐近线 L。参数 a 可设为正值或负值，以在越大越好(如投掷)或越小越好(如计时跑步项目)类型评分之间切换。

Figure 1. World record progression for the men’s 200 m. By Ar558 — Own work, CC BY-SA 3.0

为了得到原始性能 y 的变换，作者建议求解上式中的时间，即x=bln(a/(y—L)。从字面上理解很奇怪:一名运动员的表现是人类历史上本应创造世界纪录的时间。作者称这个转换后的量为“临界点”，或“Cp”。

有一个程序先得到 a 、 b 和 L. ，设定 L 为人类表现的假设极限，留下两个未知数 a 和 b 。然后用x=bln(a/(y—L))创建两个方程。第一个将左侧设置为 1000，并用世界纪录性能代替 y 。第二个将左侧设置为 0，并用健全人的表演替换 y 。有了两个方程和两个未知数，这个系统可以解出 a 和*b。*它是一个非线性系统，然而，下面的代码块 1 包含一个解算器 R 函数，它将世界纪录、健全人和人类极限表现作为输入，并返回一个包括 a 和 b 的列表。

Code Block 1. An R function to solve for a and b given assumptions about the world record, able-bodied and human-limit performances. Supports types “running” (smaller is better) and “jumping” (larger is better).

没有办法知道甚至这个非常特殊的游泳测试的世界纪录，但是假设每秒 2.0 米的人类极限(记住阻力是加进去的)、每秒 1.8 米的世界纪录和每秒 0.5 米的健全人表现导致这样的规模，即我们的五名游泳运动员集体表现接近 0 到 1000 范围的中间(类似于之前研究的链球运动员)。使用这些输入和代码块 1 中提供的 R 函数，我们得到 a = -1.5 和 b = 496.3。

get_perf_params(world_record_perf = 1.8,
                able_bodied_perf = 0.5, limit = 2.0,
                type = "jumping")

在从五个游泳者的原始数据文件中读入数据后，下面的代码块 2 也使用了get_perf_params。它还将测量值转换为标准点(Cp)标度，并绘制数据。

Code Block 2. Reading in data, followed by transformation and plotting.

Figure 2. Transformed performance measurements and daily relative training intensities for the five swimmers.

当将上面的图表与我们在第一部分和第二部分中的链球投掷者进行比较时，有两点非常突出:

• 表现趋势比我们的链球运动员更微妙，他经历了相当剧烈的表现波动。游泳者 1 号在整个研究中的进步是最容易辨别的趋势。
• 与链球运动员相比，收集成绩数据的频率较低，一些游泳运动员比其他人更频繁。

T 体能-疲劳模型提供了一个有用的起点，因为它所强加的结构，我们将在本文后面考虑样条节点放置时从中受益。代码块 3 分别拟合每个游泳者的模型，产生参数估计值(表 1)、标准误差(表 2)和覆盖数据的拟合值图(图 3)。

Code Block 3. Fitting the fitness-fatigue model for each swimmer, and plotting the fitted values.

因为性能变量中有缺失值，所以在计算平方和时，rss函数(第 30 行)现在必须使用选项na.rm = TRUE。所有日常训练强度数据都是可用的，卷积变量的计算与之前完全一样。

第 36 到 60 行是计算非线性最小二乘标准误差的特殊函数。它们通过近似 nxp矩阵来工作

就像普通线性回归中的 X 矩阵一样进入方差公式。

和以前一样，我们需要初始值。在代码块 3 的第 72 行，前面文章中的起始值被注释掉，以支持第 73 行的值。当初始化为之前的值时，游泳运动员 1、2 和 4 收敛到一个比链球运动员具有更小的体能和疲劳时间常数的解，这表明游泳运动员总体上具有更快的动力学。另一方面，游泳者 3 和 5 收敛到奇怪的地方(例如，时间常数接近 2000)。更改为较小的开始时间常数不会改变游泳者 1、2 和 4 的解决方案，但会为游泳者 3 和 5 带来更直观的解决方案，因此将用于本文的剩余部分。这些解决方案以参数估计值(表 1)和近似标准误差(表 2)的形式显示如下。

Table 1. Fitness-fatigue model solutions for Swimmers 1–5.

Table 2. Fitness-fatigue model approximate standard errors for Swimmers 1–5.

在一封讨论他的数据集的电子邮件中，克里斯蒂安·拉舍提到这项研究中的游泳运动员被认为是“精英”他解释说，由于优秀运动员已经达到了高水平的表现，并且习惯于训练，他们可能不会表现出典型的表现差异。当有明显的性能改善时，他的经验法则是“疲劳延迟(时间常数)需要大于适应延迟(大约是适应延迟的两到三倍)，疲劳与适应的比率需要为> 1。”

对于游泳运动员 1，在整个研究中至少取得了一些明显的进步，估计的健康时间常数是游泳运动员中第二大的，为 34.5 天(仍然比链球运动员的 60 天短得多)。估计的疲劳时间常数非常小，只有 1.35 天(只有链球运动员的十分之一)，表明从高强度训练中恢复非常快。尽管衰减非常快，但疲劳重量是健康重量的 4.6 倍(相比之下链球运动员是 2.6 倍)。表 2 中的标准误差计算表明了这些估计值的一些质量，尽管疲劳参数的标准误差相对于其他参数相对较高。

游泳运动员 3 号的疲劳时间常数是所有运动员中最大的，仅略超过链球运动员的 13.4 天。疲劳重量仅比健康重量大 1.5 倍，减少了较长疲劳时间常数的后果。估计的健康时间常数是 23.5 天。除截距之外，游泳者 3 的参数的标准误差估计对于所有参数都相当大。

根据模型拟合，游泳运动员 4 的健身效果微乎其微，稍微积极的健身体重由疲劳体重主导(几乎是 30 倍！).疲劳再次消失得非常快，估计的健康时间常数是五名游泳运动员中最长的，为 46 天。除了健康时间常数之外，Swimmer 4 的标准误差相对较小，这表明各种参数估计值具有一定程度的精确度。

2 号和 5 号游泳运动员的体能和疲劳都非常短暂。对于游泳运动员 5，估计的体能和疲劳影响几乎相等，他也具有最少的测试数据点。这两位游泳运动员的标准误差估计值大得离谱。请特别注意游泳运动员 5 的健身和疲劳重量的标准误差，这说明了关于这两种影响的不同信息是多么少。

承认这里推论的局限性，这些优秀游泳运动员的体能和疲劳动态似乎比我们的大学生链球运动员发生得更快。在五分之四的情况下，疲劳效应消失得非常快。疲劳效应对健康效应的相对强度在模型拟合之间变化很大，但在所有情况下(甚至是游泳运动员 5 的情况)，疲劳重量大于健康重量。

下图 3 显示了上述模型拟合的预测图。该模型能够捕捉到游泳者 1 的积极趋势，但是没有捕捉到许多其他有趣的变化。游泳运动员 2 的模型预测，在第 60 天左右的大强度训练后，成绩会略有提高，这看起来可能已经发生了(也可能是巧合)。对于 3-5 岁的游泳者来说，没什么可看的。

Figure 3. Transformed performance measurements and daily relative training intensities for the five swimmers overlaid with fitness-fatigue model fitted values.

通过直接估计与训练数据卷积的滞后分布，plines 可以提供额外的灵活性。我们在第二部分中看到了这一点，但现在是时候将它应用到真实数据中了；代码块 4 为每个游泳者分别实现该方法。

Code Block 4. Implementing the spline-based convolution approach for the five swimmers and plotting both the estimated impulse responses and the fitted values.

由于回归过程默认情况下会弹出缺少值的行，因此该过程无需修改即可运行。注意第 6 行中结位置的选择。第 2 天和第 6 天的内部结捕获了由体能-疲劳拟合建议的短期动态，而第 25 天的结则适应了对游泳者 1 观察到的长期动态。拟合的脉冲响应如图 4 所示。

这些是由于第 0 天的一次训练而随着时间影响表现的冲动。在一个经典的设置中，我们期望它开始时为负，因为疲劳效应最初超过了健身效应，随着疲劳的消散达到正的最大值，随着健身效应的衰减慢慢衰减到零。详见第二部分。

Figure 4. The fitted impulse responses for the five swimmers. Knots are represented by the crosshairs on the green horizontal line.

在图 4 中，立即明显的是相对较弱的健身效果，考虑到游泳运动员的精英地位，这并不奇怪。游泳运动员 3 和 4 根本没有明显的健身效果，游泳运动员 3 甚至没有表现出太多的疲劳效果。游泳运动员 1 具有最类似于具有相对长的体能衰减曲线的经典体能-疲劳曲线的脉冲响应曲线(在 25 天时放置的结是为游泳运动员 1 准备的)。

游泳者 2 的估计脉冲响应具有一些伪像，这些伪像可以通过定制结放置(即，在 25 天时去除结)来减少。值得注意的是，游泳运动员 2 的健身效果是五名游泳运动员中最明显的，尽管它是短暂的。作为一种推测，也许这是一个不同于游泳者 1 的“健康”现象。

给定游泳者 5 的体能-疲劳拟合的尴尬参数估计，令人惊讶的是脉冲响应看起来与其他人(尤其是游泳者 2)相似。

下面，我们可以看到基于样条的卷积方法所做的性能预测(橙色)以及参数化健身疲劳模型的预测(蓝色)。两个模型的预测非常相似。

Figure 5. Predictions made by the spline-based convolution approach (in orange) together with the parametric fitness-fatigue model’s predictions (in blue).

在上面的绳结放置过程中，我们进行了反复试验。另一种方法是在最初几周和大约 160 天的最后边界之间打一个结。图 6 显示这导致了糟糕的结果。

Figure 6. Effects of a bad knot placement for the fitted impulse responses. Knots are represented by the crosshairs on the green horizontal line.

一种解释是，80 天的时滞在研究中几乎没有观察到两次，而从训练到表现的较小时滞观察得更频繁。在如此低的信息点放置一个结对我们的分析没有任何好处，并对这种基于样条的卷积方法提出了警告:在真正探索的情况下(没有参数模型指导我们)，一个不幸的结放置可能会导致我们产生不切实际的脉冲响应。明智的做法是关注时滞轴上可用的信息量。

真实数据总是会带来与模拟数据不同的体验，我们很幸运，这种体验和模拟数据一样相似。我们确实发现了游泳者的健康和疲劳效应，尽管这些效应的强度和时间尺度有很大的不同。标准误差表明在许多参数估计中存在高度的不确定性，我们对起始值和结选择的试错法增加了额外的不容易量化的手术变异。

我们将原始性能变量(米/秒的速度)转换为标准点(Cp)尺度。这引入了一个有趣的非线性转换，并创造了一种连续性的感觉，即使运动发生了变化。我们对值 a 、 b 和 L 做了强有力的假设，但是这些似乎不太可能从质量上改变分析。只需做一些小的调整就可以处理日常级别的“缺失”性能数据。

在第四部分中，累积影响的建模将最终超越基于卷积的特征，转向使用动态形式。这将使我们更接近标准的时间序列方法。在这个旅程的第一阶段，将通过使用潜在状态向量来适应累积影响，使得卡尔曼滤波器成为处理状态空间形式的健康-疲劳模型的优秀工具。当我们使用 statsmodels 工具包中的 MLEModel 类时，这也将激励我们将 R 留给 Python。

[1] D .科洛萨、M .宾.阿兹哈尔、c .拉什、s .恩德勒、f .哈纳卡姆、a .费劳提和 M .普发弗，使用带有卡尔曼滤波器反馈的健身疲劳模型进行性能评估 (2017)，《国际体育计算机科学杂志》

[2] T. Busso、R. Candau 和 J. Lacour，疲劳和健康，模拟训练对表现的影响 (1994 年)，《欧洲应用生理学杂志》

[3] R. Morton，R. Fitz-Clarke，E. Banister，模拟人类在跑步中的表现 (1985 年)，《应用生理学杂志》

用 Python 为一家电子商务公司建模客户流失

原文：https://towardsdatascience.com/modeling-customer-churn-for-an-e-commerce-business-with-python-874315e688bf?source=collection_archive---------5-----------------------

Photo by mohamad_hassan on Pixabay

如何使用寿命包来预测非合同流失风险

留住现有客户比获得新客户更具成本效益，这就是为什么跟踪高流失风险(流失)的客户并针对他们制定保留策略非常重要。

在这个项目中，我将基于巴西电子商务网站 Olist 的数据建立一个客户流失模型。我会用它来识别高风险客户，并告知保留策略和营销实验。

电子商务有一个复杂之处。虽然衡量合同业务(基于订阅的业务)的客户流失很简单，但在非合同业务(电子商务)中却看不到客户流失。在这些情况下，概率模型在估计客户死亡时间时会派上用场。我将使用的概率模型是寿命包中的 BG/NBD 模型。

加载数据

import pandas as pd
import numpy as np
import datetime as dt
import seaborn as sns
import matplotlib.pyplot as pltfrom lifetimes.utils import *
from lifetimes import BetaGeoFitter,GammaGammaFitter
from lifetimes.plotting import plot_probability_alive_matrix, plot_frequency_recency_matrix, plot_period_transactions, plot_cumulative_transactions,plot_incremental_transactions
from lifetimes.generate_data import beta_geometric_nbd_model
from lifetimes.plotting import plot_calibration_purchases_vs_holdout_purchases, plot_period_transactions,plot_history_aliveorders = pd.read_csv(‘brazilian-ecommerce/olist_orders_dataset.csv’)
items = pd.read_csv(‘brazilian-ecommerce/olist_order_items_dataset.csv’)
cust = pd.read_csv(‘brazilian-ecommerce/olist_customers_dataset.csv’)

生命周期包依赖于最近频率货币(RFM)分析来模拟客户流失和客户生命周期价值(CLV)。为了制作我们的模型，我们需要一个由最近、频率和货币列组成的数据框架。每个的定义如下。

• 最近:初次购买和最近(最后一次)购买之间的时间
• 频率:客户重复购买的次数(总购买次数— 1)
• 货币:购买的总支出

客户 ID 信息将来自cust。订单日期将从orders开始。价格将来自items。

构建数据集

print(cust.columns)

有两列用于标识客户。customer_id是为每个订单生成的客户 ID 令牌。如果同一个客户下了多个订单，他就有多个 customer_id 标识符。我们希望在这个分析中使用的是customer_unique_id，它对每个购买者都是唯一的，可以用来跟踪他们一段时间内的购买情况。

这是顾客购物的分布情况。

cust.groupby('customer_unique_id').size().value_counts()

大多数顾客只进行了一次购买。

orders = pd.merge(orders,cust[['customer_id','customer_unique_id']],on='customer_id')
orders.columns

在items数据集中，订单中的每一项都有单独的一行。price栏指的是累计订单购买量，而不是单项价格。因为我只需要订单价格，所以我会保留每个订单的第一个项目。

print(items.columns)

items.drop_duplicates('order_id',keep='first',inplace=True)

接下来，我将把订单与商品连接起来，以附加价格信息。

transaction_data = pd.merge(orders,items,'inner','order_id')
transaction_data = transaction_data[['customer_unique_id','order_purchase_timestamp','price']]## convert timestamp to date; only need the day
transaction_data['date'] = pd.to_datetime(transaction_data['order_purchase_timestamp']).dt.date
transaction_data = transaction_data.drop('order_purchase_timestamp',axis=1)
transaction_data.head()

现在我有了交易数据，我想用我在引言中提到的 RFM 变量将它转换成数据帧。生存期包具有将事务数据转换为 RFM 数据帧的功能。

summary = summary_data_from_transaction_data(transaction_data,'customer_unique_id','date',monetary_value_col='price',)
summary.describe()

了解 RFM 数据框架

summary.head()

summary 函数将客户交易转换为汇总表。许多客户的频率、新近度和货币= 0，如客户0000366f3b9a7992bf8c76cfdf3221e2。这是因为 Lifetimes 只考虑重复购买的顾客。

使用天数作为时间段(也可以定义为周或月)，生命周期模型的变量定义如下:

• frequency:客户重复购买的天数
• T:客户年龄天数
• recency:客户最近一次购买时的年龄(天数)
• monetary_value:客户购买的平均值，不包括第一次购买

frequency不包括客户的第一次购买，因为那是客户出生的那一天。之后，你可以开始质疑那个客户是否还活着。

summary[summary['frequency']>0].head()

transaction_data[transaction_data['customer_unique_id']=='004288347e5e88a27ded2bb23747066c']

请注意客户004288347e5e88a27ded2bb23747066c如何使用 Olist 进行了两次购买，但是根据频率和货币值的定义，他的frequency是 1 而monetary_value是 87.90 美元。

可视化 RFM 数据框架

我们将使用贝塔几何/NBD (BG/NBD)模型来分析客户流失。BG/NBD 模型是对帕累托/NBD 模型的改编。两个模型都描述了顾客流动不明显的企业中的重复购买模式；然而，BG/NBD 在计算上更加可行。

BG/NBD 模型的假设:

• 客户关系有两个阶段:“存活”一段时间，然后“死亡”
• 在活动期间，客户进行的交易数量遵循交易速率为λ的泊松分布
• λ的不均匀性遵循伽马分布
• 在任何交易之后，客户以概率 p 死亡；客户在多次交易后死亡的概率服从几何分布
• p 服从贝塔分布
• λ和 p 因客户而异

有关 BG/NBD 模型的更多信息，请查看彼得·法德尔的这篇论文和卡姆·戴维森-皮隆的这篇帖子。

bgf = BetaGeoFitter(penalizer_coef=0.0)
bgf.fit(summary['frequency'], summary['recency'], summary['T']);

plot_frequency_recency_matrix(bgf);

plot_probability_alive_matrix(bgf);

训练模型并评估模型性能

接下来，我们要评估该模型，看看它在未来的表现如何。我会将数据分为训练(校准)期和维持(观察)期，训练 BG/NBD 模型，并使用 Peter Fader 在本次演讲(@ 26:10)中概述的四幅图来评估性能。这些地块是:

1)校准周期直方图:模型是否符合训练数据？

2)累计成交图:模型对累计销售额的预测效果好吗？

3)增量交易图:模型是否捕捉到了交易的整体趋势？

4)条件期望图:模型能否根据训练数据预测客户的购买数量？

1)校准周期直方图

plot_period_transactions(bgf).set_yscale('log');

该模型相当好地代表了直到四次重复交易的真实数据。很少有顾客购买更多的东西。

2)累积交易图

summary_cal_holdout = calibration_and_holdout_data(transaction_data, 'customer_unique_id', 'date',calibration_period_end='2017-09-03', observation_period_end='2018-09-03' )

我们可以通过绘制两者来评估数据集的工作方式。

bgf.fit(summary_cal_holdout['frequency_cal'], summary_cal_holdout['recency_cal'], summary_cal_holdout['T_cal'])plot_cumulative_transactions(bgf, transaction_data, 'date', 'customer_unique_id', 730, 365);

红线代表左侧校准周期和右侧维持周期之间的边界。如你所见，BG/NBD 模型在预测累积交易方面做得非常出色。

3)增量交易图

plot_incremental_transactions(bgf, transaction_data, 'date', 'customer_unique_id', 730, 365);

该图表明，该模型在捕捉数据的总体趋势方面做得不错。

4)条件期望图

plot_calibration_purchases_vs_holdout_purchases(bgf, summary_cal_holdout);

该模型在三个校准周期内性能良好，但由于数据的分布，与维持数据存在偏差。

cust.groupby('customer_unique_id').size().value_counts()

只有不到 1%的顾客进行了四次或更多的购买，所以 BG/NBD 模型没有太多的数据来了解那些进行了多次重复交易的顾客。

实际上，如果我要继续对客户流失建模，我会考虑收集更多的数据。但是出于学习的目的，预测客户流失仍然是一个很好的练习。

客户流失风险建模

BG/NBD 模型假设只有在重复购买后才会发生死亡，因为顾客在购买过程中离开，而第一次购买被保留下来以表示顾客的出生。

正因为如此，只有一笔交易的客户，会有 100%的概率是活着的，这是值得怀疑的。考虑到这一限制，我们将只预测至少进行过一次重复交易的客户的流失风险。

df = summary[summary['frequency']>0]
df['prob_alive'] = bgf.conditional_probability_alive(df['frequency'],df['recency'],df['T'])
sns.distplot(df['prob_alive']);

从这里，我们可以根据客户“活着”的概率来可视化他们。利用领域知识，我们也许能够为可能已经流失的客户设定一个阈值，并且还能识别出处于流失风险中但尚未消失的客户。

接下来，我想为客户流失设定一个决策阈值。这是一个注入个人专业知识或与领域专家交流的机会。假设我与销售和营销经理交谈，我们同意考虑一位客户，感谢他阅读！如果你想跟着代码走，这里有GitHub repo。随时欢迎反馈和问题。下面评论！

df['churn'] = ['churned' if p < .1 else 'not churned' for p in df['prob_alive']]
sns.countplot(df['churn']);

A little over 92% of customers have churned, meaning that there’s a lot of opportunity for improvement regarding retention.

We can assume that the customers who have churned are already lost. But what is interesting in a business setting is the customers who are at high risk for churn, but haven’t churned yet. Later on, it might still be a good idea to apply different treatments to the churned group.

If I can identify them, maybe I can encourage the marketing team to target them with promotions.

sns.distplot(df[df['churn']=='not churned']['prob_alive']).set_title('Probability alive, not churned');

It seems reasonable to bucket customers with 80% or more churn risk to be considered high risk for churn.

df['churn'][(df['prob_alive']>=.1) & (df['prob_alive']<.2)] = "high risk"
df['churn'].value_counts()

Now that I have these churn groupings, I can move forward and apply special treatments to these groups. Ideally there would be more data and a bigger population of high-risk customers.

Conclusion

We’ve modeled churn risk in a non-contractual setting, and now have three customer segments — not churned, high risk, and churned. This could feed into a dashboard to give stakeholders a glimpse of “at-risk” customers. It also provides three different groups that we can run specific actions. Some ideas:

1. Reach out to churned customers to figure out why they left.

2. Send different types of targeted emails and special offers to the high risk group. If the sample size of high risk customers is large enough, you could split off a few small treatment groups and compare how their retention and CLV change with different promotional or customer relationship strategies.

3. Determine the the highest value customers in the non-churn group, and serve them additional benefits to ensure that they remain loyal customers.

Thanks for reading! If you want to follow along with the code, here’s the GitHub repo. Always open to feedback and questions. Comment below!

使用脸书预言家为推特和金融市场中的月亮周期建模

原文：https://towardsdatascience.com/modeling-lunar-cycles-in-tweets-and-financial-markets-using-facebook-prophet-d6ec0e9e20f?source=collection_archive---------26-----------------------

我最近读了一篇文章，认为满月时犯罪率会上升。出于对这种效应的好奇，我做了一点研究，似乎月亮周期是否影响人类行为的问题已经在学术界激烈争论了几十年(见塔库尔和夏尔马，1984 和罗顿和凯利，1985 反对观点的例子；这场争论的双方都有很多学者)。不满意，我决定通过检查一年的推文和财务数据来进行调查。

概观

这个项目的完整 Python 代码可以在这个 GitHub 库中找到。如果你想玩玩我的发现，我还做了一个交互式仪表盘。

为了这个项目，我使用各种关键词搜索短语在 2018 年的每一天收集了大约 1000 条推文。这一过程产生了每个搜索短语超过 365，000 条推文的数据集。然后，对产生的数据集进行情绪处理，并按日期汇总，以检查情绪随时间的趋势。此外，2018 年各种股票、货币和市场指数的每日金融数据都是从互联网上收集的。

然后，Twitter 和金融数据都使用脸书预言家进行建模，以确定与月亮可见阶段相关的任何季节性的程度(或缺失)。事实证明，月亮周期和人类行为之间确实存在关联，尽管影响的大小并不明显。

获取数据并探索它

当使用 twitterscraper 库时，抓取 Tweets 是一个相当简单的过程(警告，这可能需要几个小时才能完全执行):

**import** pandas **as** pd **import** json **from** twitterscraper **import** query_tweets
**import** twitterscraper def scrape_tweets(query, year=2018, num_tweets=1000):
    *"""scrapes X Tweets per day for a year.*
 *Scraping works backwards from midnight.*
 *'query' must be a string.*
 *------------------------------------*
 *Generates a JSON file for each day scraped.*
 *"""*
    dates_year = [str(date)[:10] for date in pd.date_range(start=f'1/1/{year}', end=f'12/31/{year}')]
    for i in range(len(dates_year)):
        begin_date = dates_year[i]
        if i == len(dates_year)-1:
            end_date = f'{year+1}-01-01'
        else:
            end_date = dates_year[i+1]
        day = dates_year[i]
        cmd = 'twitterscraper "{}" -l {} -o t{}.json -bd {} -ed {} --lang en'.format(query, num_tweets, day, begin_date, end_date)
        subprocess.run(cmd)
        if (i+1)%5 == 0:
            print(f"finished scraping {i+1} days of {year}")
    print("SCRAPING PROCESS COMPLETE!")
    pass

使用 alphavantage.co 的 API 搜集财务数据也相当容易，尽管这些数据需要一些额外的条件使其更有用。

**import** requests
**import** matplotlib.pyplot **as** plt
plt.style.use('ggplot')
%matplotlib inlinedef get_financial_data(symbol, stocks_apikey, column="close_24", year=2018, verbose=True):
    *"""Inputs:*
 *symbol         (string) Stock or Currency symbol*
 *stocks_apikey  (string) your API key for* 
 *https://www.alphavantage.co
    column         (string) column of output to use for plotting
    year           (int) the year you wish to examine*
 *=========================================*
 *Returns a DataFrame of daily financial information containing*
 *at least opening, closing, high, and low values.*
 *"""*
    valid_types = ['stock','index','currency','cryptocurrency']
    credentials = {'function':'TIME_SERIES_DAILY',
                       'symbol':symbol,
                       'outputsize':'full',
                       'apikey':stocks_apikey}
    r = requests.get('https://www.alphavantage.co/query', params=credentials)
    df = pd.DataFrame(r.json()["Time Series (Daily)"])
    df = df.T.reset_index()
    df.columns = ['date','open','high','low','close','volume']
    df.date = pd.to_datetime(df.date)
    df[['open','high','low','close','volume']] = df[['open',
            'high','low','close','volume']].astype(float) *# create a new column to account for after-hours trading*
    *# uses the next day's open value as the prior day's close value*
    cl24 = [df.loc[0].close]
    for val in df.open.values:
        cl24.append(val)
    cl24 = pd.DataFrame(cl24[:-1], columns=['close_24'])
    df = df.join(cl24)
   *# now account for after-hours trading exceeding high/low values*
    df['high_24'] = df[['high', 'close_24']].values.max(1)
    df['low_24'] = df[['low', 'close_24']].values.min(1)
    df['range'] = df['high'] - df['low']
    df['range_24'] = df['high_24'] - df['low_24']
    df['change_24'] = df['close_24'] - df['open']
    df.set_index('date', inplace=True)
    year_df = df[f'{year}':f'{year}']    *# getting 1 year's data* *# plotting the results* plt.figure(figsize=(15,5))
    plt.plot(df[column], label=column)
    plt.title(f"{symbol} Daily Performance for {year}", fontsize=16)
    plt.ylabel("Price (in USD)", fontsize=14)
    plt.legend(loc='best', borderpad=1, fontsize=14);
    return df

有了数据集，下一步是确定我们收集的所有推文的总体情绪。为此， TextBlob 库非常有用。TextBlob 可以快速解析文本以确定文本本质上通常是正面的还是负面的。我选择扩展“中性”极性值，以包含从-0.1 到 0.1 的所有极性分数，从而消除+/-类别中的一些噪声。

**from** **textblob** **import** **TextBlob**def get_tweet_sentiment(tweet): 
    analysis = TextBlob(tweet)
    polarity = analysis.sentiment.polarity
    subjectivity = analysis.sentiment.subjectivity
    if analysis.sentiment.polarity > 0.1: 
        sentiment = 'positive'
    elif analysis.sentiment.polarity < -0.1: 
        sentiment = 'negative'
    else: 
        sentiment = 'neutral'
    return sentiment, polarity, subjectivity

在所有这些处理之后，我最终可以使用 statmodels 的[seasonal_decompose](https://www.statsmodels.org/stable/generated/statsmodels.tsa.seasonal.seasonal_decompose.html)方法来确定收集到的 Twitter 数据中是否存在任何季节性。你瞧，事实证明 是7.4 和 29.5 天的季节性，这与月亮周期完全吻合。

Seasonal Decomposition of scraped tweets

和脸书先知一起做模特

脸书预言家，虽然很新，却是一个非常强大的适应性建模工具。幸运的是，Prophet()包含一个名为holidays=的参数，它允许用户输入一个格式正确的数据帧，FBProphet 会特别注意这个数据帧。其中，该功能模拟了holidays对所研究的时间序列的季节性影响。在下图中，满月发生在 1 月 2-3 日，2 月 1 日，3 月 1-2 日，等等。

**# `phases` was a pre-generated DataFrame of lunar phase dates****from** fbprophet **import** Prophet
**from** fbprophet.plot **import** plot_plotly, plot_cross_validation_metric
**from** fbprophet.diagnostics **import** cross_validation, performance_metrics*#to maintain pd.plotting functionality* pd.plotting.register_matplotlib_converters()def tweet_fbprophet(filename=''):
    tweets_df = pd.read_csv(filename)
    tweets_df.timestamp = pd.to_datetime(tw_df.timestamp, format='%Y%m%d')
    grouped = pd.DataFrame(tweets_df.groupby(['timestamp', 'sentiment'])['tally'].sum()).reset_index()

    *#prepare grouped sentiment data for FBProphet processing*
    *#here, positive sentiment only*
    grp_pos = grouped[grouped.sentiment == 'positive'].drop('sentiment', axis=1).reset_index(drop=True)
    grp_pos.columns = ['ds','y']

    m = Prophet(holidays=phases)
    m.fit(grp_pos)
    future = m.make_future_dataframe(periods=60, freq='D')
    forecast = m.predict(future)
    forecast[(forecast['fullmoon'] + forecast['lastquarter'] + 
             forecast['newmoon'] + forecast['firstquarter']).abs() > 
             0][['ds', 'fullmoon', 'lastquarter', 'newmoon', 'firstquarter']][:10]
    fig1 = m.plot(forecast);
    fig1.set_size_inches(15, 5)
    plt.show();
    fig2 = m.plot_components(forecast)
    fig2.set_size_inches(15, 10)
    plt.show();
    return m

Example of seasonality effect of lunar cycles found in tweets using FBProphet

FBProphet 的另一个好处是种类繁多的performance_metrics，很容易访问。我选择研究 RMSE(均方根误差),因为它很容易用原始数据集的单位来解释。实际上，在所有情况下，(推特和金融数据)RMSE 平均比月球周期的季节性影响高出 50%。虽然从纯预测的角度来看这并不太好，但需要注意的是，该项目的目标是确定是否存在相关性，而不是创建一个可以准确预测 twitter 情绪或股票价格的模型…这需要太多的额外变量。下面是性能指标的示例表和不同时间跨度的 RMSE 图。

RMSE plot from FBProphet cross-validation method

随着 Twitter 阴历季节性的确立，我对财务数据运行了同样的 FBProphet 程序。该数据还显示了与月亮周期相关的明显的季节性。然而，对于每个金融对象，相关性并不相同。一些股票/货币在满月时上涨，而另一些在满月时下跌。此外，似乎许多月相都有一个初始峰值，然后在第二天出现相应的修正(甚至过度修正)。例如，比较下面为比特币和微软生成的农历季节性图表。

Lunar seasonality comparison of BTC and MSFT (Microsoft does not trade on weekends, hence the gaps)

寻找信息的视觉化

最初，我认为生成单词云可能是从我的数据中提取更多信息的一种好方法。毕竟，单词云可以在视觉上提供大量信息，而其他图形则不能(因为它们缺乏单词云的莫邪)。不幸的是，由于我正在处理的语料库的规模(每个关键词搜索大约 5-6，000，000 个单词)，使用极其简单的 wordcloud 库被证明是不切实际的，因为执行时间延长到 10 多个小时。然而，这是我一次运行这个过程的代码和图表。

**from** wordcloud **import** WordCloud
**from** sklearn.feature_extraction.stop_words **import** ENGLISH_STOP_WORDSdf = pd.read_csv('tweets.csv')
df.timestamp = pd.to_datetime(df.timestamp, format='%Y%m%d')  
comment_words = ' '
stopwords = set(ENGLISH_STOP_WORDS)
stopwords.update(['com','ve','ll','just','don','really','00'])*#* ***this for loop took hours to run***
it = 0
for val in df.text: 
    val = str(val)          *# typecaste each val to string* 
    tokens = val.split()    *# split the value* 
    for i in range(len(tokens)): 
        tokens[i] = tokens[i].lower() *# Convert tokens to lowercase*
    for words in tokens: 
        comment_words = comment_words + words + ' '
    it += 1
    if (it)%10 == 0:
        print(f"Completed {i} of {len(df.text)} loops")

*# generate the word cloud*
wordcloud = WordCloud(width = 1200, height = 1200, 
                background_color ='white', 
                stopwords = stopwords, 
                min_font_size = 10).generate(comment_words)*# plot & save the WordCloud image  *                      
plt.figure(figsize = (20, 20), facecolor = None) 
plt.imshow(wordcloud) 
plt.axis("off") 
plt.tight_layout(pad = 0) 
plt.savefig('images/tweets_wordcloud.png')

Word cloud for the search “love OR hate OR peace OR war”

通过在一个查询语料库中比较正面/负面/中性，或者将一个语料库与另一个语料库进行比较，可以进行快速和简单的可视化来绘制每个搜索短语的情绪的每日变化。

Unsurprisingly, tweeters had the least positive things to say regarding politics.

Daily sentiment of Tweets with no search phrase given

一个更简单的方法是创建一个单词包，只过滤出现次数最多的第 n 个 T2，然后绘制一个直方图。这种简单的方法，当与日期过滤器结合使用来区分月相日期和非月相日期时，会为每个搜索短语产生一些有趣的结果(见下面的例子)。虽然没有出现什么过分令人担忧的情况，但有趣的是看到词频的顺序是如何根据月相变化的。

Comparison of most common words when no search phrase was entered

最后，也可能是我最喜欢的，是利用散点图库来制作一个漂亮的、交互式的、可搜索的 Twitter 数据图表。这个库实际上生成一个 HTML 文件来显示图表，对于像我这样大的语料库，每个 HTML 文件平均 30MB，这导致 web 浏览器的加载时间大约为 5-10 分钟。该库是为中小型数据集设计的，不是我的怪物，但结果仍然值得等待。可视化是交互式的(鼠标悬停在每个数据点上，会出现一个信息提示),并且是可关键字搜索的。最重要的是，选择一个关键字还会在图表下方显示该词的大量分类(肯定/否定)出现，以了解它在每个场景中的常用情况。我可以立即看到，这个图书馆有很大的潜力，有利于学术研究人员，特别是在文学，历史，政治学等领域。这里有一个精简后的例子，应该可以快速加载，还有一个下面功能的. gif 文件。

以下是我用于散文本的代码。

*!pip install scattertext
!pip install spacy
!python -m spacy download en_core_web_sm* **import** scattertext **as** st
**import** pandas **as** pd
**import** numpy **as** np **import** spacy
**import** en_core_web_smdf = pd.read_csv('tweets.csv')
df.timestamp = pd.to_datetime(df.timestamp, format='%Y%m%d')  
df = df[df['sentiment'].isin(['positive','negative'])]*# Turn DataFrame into a Scattertext Corpus*
nlp = en_core_web_sm.load()
corpus = st.CorpusFromPandas(data_frame=df,
                             category_col='sentiment',
                             text_col='text',
                             nlp=nlp).build()*# Create an HTML page for the interactive visualization*
html = st.produce_scattertext_explorer(corpus,
    category='positive',
    category_name='Positive',
    not_category_name='Negative',
    minimum_term_frequency=25,   *### high values reduce load time*
    minimum_not_category_term_frequency=25,  *### high values reduce load time*
    max_terms=5000,   *### aim for at least 1000 for a pretty graph*
    max_snippets=50,
    show_characteristic=True,
    width_in_pixels=1000)open("tweet_scattertext".html", 'wb').write(html.encode('utf-8'))

最终想法和建议

基于我的分析，似乎相当清楚的是，事实上人类行为和月相之间存在关联。然而，这必须留给未来的研究人员来确定这种相关性的本质和原因。这种影响相对较小，但仍然存在。就推文中的整体情绪而言，满月和上弦月似乎与积极情绪的增加相关，而新月和上弦月与消极情绪的增加相关。

金融市场也表现出月球季节性，包括股票、市场指数、货币和加密货币。然而，这种相关性并不是全面一致的；一些股票与满月正相关，其他的负相关。此外，大多数股票/货币在月相数据中表现出在一个方向上飙升的强烈趋势，然后在第二天立即在相反方向上修正(或过度修正)。任何试图使用这一信息的人都必须在个案的基础上检查季节性。尽管相关性很小(季节性占价格变化的大约 1%)，我怀疑从事算法交易的新兴实践的投资者将能够从将月球季节性纳入他们的投资算法中获利。

用 Python 建模新闻报道。第 2 部分:用有限的数据开始新闻预测

原文：https://towardsdatascience.com/modeling-news-coverage-with-python-part-2-starting-news-forecasting-with-limited-data-5c1092de3ea9?source=collection_archive---------36-----------------------

早些时候我写了一点关于获取新闻故事开始建模，现在我将玩一玩如何预测覆盖范围和衡量一个模型是否好。

我将使用 fast.ai 的现成神经网络来开始这篇文章，使用一组出版商的一些文章来生成一些出版商日常生产的模型，评估模型的拟合度，然后展示为什么这些模型对于我们正在处理的数据集的限制来说不是那么好(完全任意的限制以使代码易于运行，打乱粒度和添加更长时间的文章将会获得更好的结果)。然后，这将使用自动回归和移动平均，看看对于我们利用的有限数据，它如何优于神经网络，从而生成每个出版物的覆盖趋势，如下所示:

SARIMA models for select publications for select coverage with error on test predictions

我们开始吧！

检查我们是否还有数据

上次我们从 GDELT(【https://www.gdeltproject.org】T2)获取了 60 天的新闻报道。我们将开始确保我们仍然有它，如果我们没有，我们将下载它。

现在我们已经确认了我们需要的所有文件，我们可以像以前一样浏览它们，并为《纽约时报》、《华盛顿邮报》、《福克斯新闻频道》和《CNN》提取文章。

运行代码后，我们的 df 中就有了新闻文章。接下来，我们可以细分每个出版商每天在每个目标国家的覆盖范围。

首先，我们将确保该日期是一个好的日期时间，并将我们的索引设置为该日期。

df.DATE = df.DATE.apply(lambda x: str(x))
df.DATE = pd.to_datetime(df.DATE)
df.fillna("", inplace=True)
df.set_index("DATE", drop=True, inplace=True)

接下来，我们将对 df 进行一次性编码，以确定该文章是否被出版商和日期标记为与我们的某个国家或集团相关。

df["dprk"] = df["LOCATIONS"].apply(lambda x: x.find("North Korea") > -1)
df["ukraine"] = df["LOCATIONS"].apply(lambda x: x.find("Ukraine") > -1)
df["russia"] = df["LOCATIONS"].apply(lambda x: x.find("Russia") > -1)
df["iran"] = df["LOCATIONS"].apply(lambda x: x.find("Iran") > -1)
df["china"] = df["LOCATIONS"].apply(lambda x: x.find("China") > -1)loc_df = df.groupby(["SOURCES", "DATE"])[["dprk", "ukraine", "russia", "iran", "china"]].sum()

所以现在我们可以创建一个大的 df，我们称之为时间序列来保存每个出版商和国家的日期值，其中的列都有这样的标签。

time_series = pd.DataFrame()
for publisher in mySources:
  time_series = pd.concat([time_series, loc_df.ix[publisher].add_prefix(“{}_”.format(publisher))], axis=1)time_series.plot(figsize=(20,10))

很好，所以你现在应该有类似的东西，我们可以看到数据有每周的季节性。此外，乌克兰的东西大约在 9 月 23 日起飞。

因此，让我们通过神经网络运行这些时间序列，看看它会得出什么结果！

使用 Fast.ai 中的表格模型

所以现在我们将使用 fastai 的表格模型(https://docs.fast.ai/tabular.html)。我使用表格模型是因为它快速而简单，即使时间序列太短，不能期望一个很好的结果，我也不打算摆弄得到好的结果，但它能表达重点。

注意到季节性，我们希望模型查看最近几天的数据。由于是每周一次的季节性数据，我们可以只给出最近 7 天的数据，希望这能帮助它发现数据中有季节性成分。此外，我们可以给出文章是在一周的哪一天发表的，这可能会有所帮助。例如，对于第 X 天朝鲜的福克斯新闻频道，我们需要第 X-1 天、第 X-2 天的所有发布数据，依此类推。所以我们可以创建一个函数来完成这个任务。

然后，我们可以创建一种方法，让我们从我们的时间序列数据帧中选取一个序列，并训练一个神经网络来尝试对其建模。

为了演示的目的，我们可以选择一个时间序列，看看它生成的模型是什么样的。我们将以相当大的学习速度快速浏览一遍，看看它是什么样子。如果你不熟悉 fastai，去参加这个课程吧。否则，你基本上是想看看你用 fit_one_cycle 数打高尔夫能打多好，虽然这只是一个任意的例子。

learn = genNN(time_series, “foxnews.com_russia”)
learn.fit_one_cycle(1, 3e-02)

因此，为了可视化结果，我们需要从我们的学习者那里提取预测，并根据实际结果绘制图表。我们可以对训练集、验证集和测试集进行不同的着色，以了解模型不同部分的表现。

那么我们的模式是好是坏呢？10 月份的准确性可能是侥幸。如果你注意到了，我们还提取了训练集的平均值，这样我们就可以比较到目前为止我们是否比机会做得更好。我们将使用均方根误差来评估模型的表现。

>>> RMSE training  using mean: 7.255011530500299
>>> RMSE training using model: 7.2815201128101155
>>> RMSE testing using mean: 8.930425521776664
>>> RMSE testing using model: 6.544865479426577

所以这个模型实际上并不太糟糕。它在预测训练集将会发生什么方面做得非常糟糕，但它以某种方式计算出了 10 月初的报告高峰。(不知道为什么)。但是让我们重新运行模型生成，看看它是否会再次发生。我们应该得到不同的结果，因为模型开始时一切都是随机的，而我们只是遵循梯度。

learn = genNN(time_series, “foxnews.com_russia”)
learn.fit_one_cycle(1, 3e-02)

并重新运行代码以显示图形结果:

重新运行评估，我得到的结果是，通过取训练集的平均值，它的表现比之前的稍差。

>>> RMSE training  using mean: 7.255011530500299
>>> RMSE training using model: 7.314209362904273
>>> RMSE testing using mean: 8.930425521776664
>>> RMSE testing using model: 7.068090866517802

我们可以改变学习率，让它变得更好，

learn = genNN(time_series, "foxnews.com_russia")
learn.fit_one_cycle(1, 2e-02)
learn.fit_one_cycle(5, max_lr =3e-04)
...
>>> RMSE training  using mean: 7.255011530500299
>>> RMSE training using model: 9.782066217485792
>>> RMSE testing using mean: 8.930425521776664
>>> RMSE testing using model: 5.876779847258329

但总的来说，这并不是一个好的模型。但话说回来，这不应该是因为它只使用了几个出版商和几个普通国家的 60 天数据。

那很有趣。我们训练了一个神经网络，但它不是很聪明；尽管我们每天都给它喂食，但它似乎并不知道每周的季节性。让我们继续，如果我们使用更传统的模型来研究移动平均线的自回归，会发生什么。(我想我们应该从这里开始)。

季节性、自回归和移动平均线

用 python 分解时间序列非常容易。假设您仍然加载了时间序列数据框架，我们可以使用 stats models(【https://www.statsmodels.org/stable/index.html】)将每一列分解为其观察成分、趋势成分、季节成分和剩余成分，并绘制它们。

import statsmodels.api as sm
res = sm.tsa.seasonal_decompose(time_series["foxnews.com_russia"])
res.plot()

Foxnews Coverage including Russia Decomposed

它腐烂了！

我们还可以查看自相关，或者值与先前值的相关程度。

from statsmodels.graphics.tsaplots import plot_acf
plot_acf(time_series["foxnews.com_russia"])

我们可以看到季节性因素和趋势。

我们可以看看移动平均线。我们将取前两天的两天简单平均值。

f_ru = time_series[["foxnews.com_russia"]].copy()
f_ru.columns = ["actual"]
f_ru["sma"] = 0for i in range(0,len(f_ru)-2):
  f_ru.iloc[i+2,1] = np.round(((f_ru.iloc[i,0]+
  f_ru.iloc[i+1,0])/2),1)f_ru.plot(title="Fox Russia prior 2 days SMA")

我们可以检查一下 RSME

f_ru = f_ru[2:]
print("RMSE testing using SMA2: {}".format(sqrt(mean_squared_error(f_ru.actual, f_ru.sma))))>>> RMSE testing using SMA2: 7.666822836990112

所以我们有一个比上面的稍微差一点的 RSME，但是要少得多(除了我忘记了我必须使用它。复制()让熊猫不抱怨)。不算太寒酸，但其实比只取中庸还不如。大概吧。与这次相反，我实际上没有检查整个持续时间。因此，让我们使用 Python Statstools SARIMA 包来集成季节性、自回归和移动平均。

所以我们要关注一下这个的参数。基本上，我们将需要自回归移动平均组件和时间序列组件作为模型的订单和季节性订单。我们将假设它也是每周一次的季节性，这通常是正确的。

所以现在开始。首先，我们声明并拟合一个模型。

from statsmodels.tsa.statespace.sarimax import SARIMAX
issue = “foxnews.com_russia”
train_l = 55s_model = SARIMAX(endog = time_series[[issue]][:train_l][1:], 
 exog = time_series[[x for x in time_series.columns if x != issue]][:train_l].shift().add_suffix(“_l1”)[1:], order=(2,0,1), seasonal_order=(1,0,1,7)).fit()

现在我们想去得到我们的预测。我们可以使用以前计数的真实值，因为我们在现实世界中会有这些值，除非我们预测得非常远(但我们之前没有这样做)。

f_ru = time_series[[issue]].copy()[1:]
f_ru.columns = ["actual"]f_ru["predicted"] = s_model.predict(end=datetime.datetime(2019, 10, 6), endog = time_series[[issue]][-5:],exog = time_series[[x for x in time_series.columns if x != issue]].shift()[-5:], dynamic= False)#plotting
ax = f_ru["actual"].plot(title="Russia Mentions in Fox", figsize=(12,8))
f_ru["predicted"][:-5].plot(color="orangered")
f_ru["predicted"][-6:].plot(color="red")

所以我们有了一个看起来相当不错的 SARIMA 模型！RMSEs 怎么样了？

testing = f_ru[:-5].copy()
print("RMSE training using model: {}".format(sqrt(mean_squared_error(testing.actual, testing.predicted))))testing = f_ru[-5:].copy()
print("RMSE testing using model: {}".format(sqrt(mean_squared_error(testing.actual, testing.predicted))))>>> RMSE training using model: 3.098788909893913
>>> RMSE testing using model: 6.281417067978587

所以这里的 SARIMAX 比上面的神经网络做得更好。如果我们愿意的话，我们可以在一个循环中为我们正在寻找的所有组合重新运行代码(并且您可以添加您喜欢的发布者和主题的任何组合)。只是提醒一下，代码会发出很多警告。

It doesn’t do too well on the Ukraine forecasting but that should be expected because it didn’t fit the previous trends

尝试调整模型的参数，并围绕出版物和报道主题进行更改。或者更好的是，回去下载更长时间的文章，看看更多的数据如何改善训练和测试。

继续前进

在第一部分:介绍和这之间，你应该对我们能理解多少报纸之间的相互作用有一个很好的想法。通过对一段时间内的新闻报道进行建模，我们可以评估一个模型是否比其他模型更好地包含了影响新闻报道的因素。

这有什么关系？嗯，当人们跑来跑去说新闻在做这个或那个的时候，如果我们没有模型来评估他们的说法，真的很难分辨谁说的是真话还是假话。但是在仅仅几分钟内(希望如此),我们就能够获取数据，运行计算模型，并根据媒体报道预测未来的报道。所有这些都与关于议程设置的政治学/传播学/新闻学研究有关(例如:https://www . tandfonline . com/doi/full/10.1080/15205436.2014 . 964871)，但实际上预测新闻报道并不像我想象的那样广泛，我不确定这是为什么。

现在，这里无疑有很多数学问题，我把这些留给读者自己去解决(或者写下他们的结果:)数据也有问题；我们使用了 100，000 个可用新闻来源中的 4 个，从主题、地点和人物的可能组合的非常大的域中选择了 5 个地点，以及一年中只有 60 天的报道。我们也观察了这里的许多参数；最终，你不希望只关注一个主题，而是不断地监控所有的主题。然而，希望这足够简单，任何人都可以拿起它，搞乱他们感兴趣的任何东西。

我不完全确定我将在第 3 部分中涵盖什么。在谷歌搜索或各种 Twitter 数据的操作化问题上，可能包括公众利益的措施可能会很有趣，尽管一篇博客帖子的困难部分是确保数据易于访问。也可能做更长期的趋势。我刚刚开始写这些来拖延期中考试的学习，所以还没有想好。

用 Python 建模新闻报道。第 3 部分:新闻网站与 Google 搜索趋势的交互

原文：https://towardsdatascience.com/modeling-news-coverage-with-python-part-3-newspaper-coverage-and-google-search-trends-dca50f60a9c8?source=collection_archive---------22-----------------------

Fitting models of Google Search to Search Trends and News Articles

这篇文章整合了来自有限的报纸报道样本和谷歌搜索趋势的数据，以模拟两者之间的互动。在这些例子中，初步分析发现，新闻报道对于预测搜索趋势是有用的，但在另一个方向上，结果却是小而杂的。这些例子是在博客文章中用有限的资源和时间选择完成的，但是一切都足够通用，你可以交换你自己的主题或数据，这应该足以让你开始。

这是用 Python 演示新闻建模的系列博文中的第三篇。本系列的第一篇文章着眼于什么是时间序列，以及媒体中的交互可以被建模。的第二篇帖子查看了一些出版物和主题的新闻报道，用不同的方法预测报道。

装载物品

和以前一样，模型的所有必要代码都在本文中提供。它被格式化后可以复制粘贴到 Jupyter 笔记本文件中。该代码将获取文章数据，但是您可能希望按照说明和链接下载 Google 搜索趋势数据，或者如果您熟悉的话，自己下载。

我们将从使用 Python 包(【https://www.gdeltproject.org/data.html#rawdatafiles】T4)从 GDELT 获取 60 天的文章开始。这些文章都标有主题、人物和地点。这将检查文件是否已经下载，如果没有，则下载它们，并将其格式化为 pandas 数据帧，以便我们可以操作该系列。

现在我们将通过一个 SARIMA 模型来运行它，这个模型带有我们上次观察过的参数。在上一篇文章中，我们使用均方差来评估准确性，但这次我们将使用平均绝对百分比误差(MAPE)。这让我们可以对不同值的时间序列的误差进行标准化。我们将跟踪 mape_df 数据帧中的 MAPE 值，以便能够比较结果。

我们测试的第一个模型将同时使用出版物中其他国家的覆盖范围和其他出版物中该国家的覆盖范围，以符合之前的做法，确保一切都可以再次工作。

然后，我们可以更改外生系列，只从同一出版物中选择适合模型的文章，将结果保存到 mape_df dataframe 列“WithinNewspapersOnly_predicted”。

我们可以比较两者之间的 mape_df 差异:

mape_df["news_dif"] = mape_df["WithinNewspapersOnly_predicted"] - mape_df["NewspapersOnly_predicted"]for country in myCountries:
  print("{} mean hat change of {}".format(country, mape_df[mape_df.index.map(lambda x: x.find(country) > -1)]["news_dif"].mean()))>>> dprk mean hat change of -10.75
>>> ukraine mean hat change of -1.5
>>> russia mean hat change of 6.5
>>> iran mean hat change of -3.0
>>> china mean hat change of -13.75

看起来似乎大多数时候只使用同一份出版物比使用其他出版物的报道要好。也就是说，那是针对这个数据集和一小段时间的覆盖范围。好吧，那么谷歌搜索趋势会发生什么？下一步，我们将模型与搜索趋势相匹配，然后将新闻模型与新闻数据和搜索趋势相匹配。

谷歌搜索趋势

我们首先需要获得谷歌搜索数据。谷歌在 https://trends.google.com 发布他们的搜索趋势。如果你不熟悉的话，那里有很多有趣的功能。为了加快速度，我去此链接搜索了 2019 年 8 月 8 日至 2019 年 10 月 6 日相关国家(朝鲜、中国、俄国、伊朗和叙利亚)的趋势。

该页面应如下所示:

Google Trend Search Results

谷歌给出的搜索趋势是相对于搜索查询中最高搜索项的百分比，而不是原始搜索数字。在这里，乌克兰获得了最高数量的搜索结果，因此所有其他数据点都以乌克兰最大搜索量的百分比表示。我们可以点击“利息随时间变化”标签右侧的向下箭头下载 CSV 格式的文件。

然后，我们可以将 CSV 放入我们当前的工作目录中(或者其他位置，如果您想更改代码的话)，并使用 pandas 加载它。注意:Google 的 CSV 的第一行是某种元数据标签，所以我们跳过它，以便在加载时正确格式化列。

gdata = pd.read_csv("multiTimeline.csv", skiprows=2)
gdata.plot()

目测显示它看起来与网站上的原始数据相同，因此它似乎已经正确加载。现在，我们可以将其格式化为与我们之前的新闻媒体数据框相同的格式，并将它们放在一起，创建原始时间序列 df 的副本作为原始数据，以防我们搞砸了什么，需要返回。

gdata.columns = [x.replace(": (United States)", "").replace("North Korea", "dprk").lower() for x in gdata.columns]gdata.day = pd.to_datetime(gdata.day)gdata.set_index("day", inplace=True, drop=True)original = time_series.copy()time_series = time_series.join(gdata)

因此，现在我们可以继续前进，绘制报纸报道和谷歌搜索趋势，以直观地检查他们有多少排队。

我们可以先做一些观察。一些搜索趋势有明显的每周季节性(如中国)，但对于其他搜索趋势(如乌克兰、伊朗)，季节性不如那些随着新闻报道增加而出现的大峰值重要。这种差异会有点问题，因为搜索似乎忽略了朝鲜的新闻报道趋势，对乌克兰的第二次飙升做出了反应，对伊朗的第三次飙升做出了反应。

我们可以快速检查新闻对搜索是否有影响。为了建立基准错误率，我们可以继续在搜索结果上运行与新闻网站相同的 SARIMA 模型，只使用以前的 Google 搜索趋势和新闻媒体对模型相同主题的报道。

朝鲜和乌克兰在这一点上的误差最大。正如我们之前看到的，谷歌搜索偏离了新闻结果。果不其然，该模型从乌克兰的第一次高峰中吸取了经验，并预测到搜索趋势将再次出现高峰，而实际情况并非如此。此外，搜索趋势也忽略了朝鲜的新闻内容。我们真的没有任何东西可以预测这种数据中的行为，所以不去费心做分解，但更长时间的数据可能会让我们知道这是多么正常。

由于每个出版物每天的文章总数通常是相同的，覆盖一个国家可能会导致另一个国家的覆盖面减少(除非在同一篇文章中同时覆盖两个国家)，所以我们可以看看是否只包括所有内容会改进模型。

看起来，加入所有其他时间序列提高了对朝鲜和乌克兰的预测，而对俄罗斯没有影响，对伊朗和中国则有负面影响。更长的时间肯定有助于梳理出我们真正想要的数据类型，但我不想在这些博客文章中建立数据库，这已经足够说明问题了。

看起来新闻趋势帮助我们预测搜索趋势。反过来成立吗？

将谷歌搜索添加回新闻模式

到目前为止，我们已经了解了新闻如何预测谷歌搜索趋势。但是反过来呢？

因为我们已经将 Google 搜索趋势添加到 time_series 数据帧中，所以我们可以重新运行原始代码，稍微修改一下，只包含报纸，然后在 dataframe 列中保存 MAPE 值。

我们可以检查新模型在精确度上与原始模型有何不同，

mape_df["overall_change"] =mape_df["WithSearch_predicted"] - mape_df["NewspapersOnly_predicted"]for country in myCountries:
  print("{} mean hat change of {}".format(country, mape_df[mape_df.index.map(lambda x: x.find(country) > -1)]["overall_change"].mean()))>>> dprk mean hat change of 0.25
>>> ukraine mean hat change of -2.5
>>> russia mean hat change of 7.0
>>> iran mean hat change of 5.0
>>> china mean hat change of -2.25

有了这些参数和数据，加入谷歌搜索趋势对模型有混合但小的影响。不过这也没什么，对于 60 天的数据，我们也不指望会有很好的结果。

最后

我们看到，有时新闻可以帮助预测谷歌搜索趋势，但并不总是如此。SARIMA 模型在搜索趋势偏离新闻报道趋势的主题上表现更差，这是意料之中的。如果没有之前乌克兰峰值的例子，我们不会指望这个模型真的能猜到覆盖范围会有这么大的差异。

我们没有从将搜索结果添加回新闻模型中得到很好的结果。那很好。这些模型的参数并没有针对任何搜索词的细节进行调整，而且涵盖的时间也很短。看看新闻->搜索和搜索->新闻之间的关系如何维持更长的时间会很有趣。此外，这些都是相当广泛的主题。在下一篇文章中，我们将看看每一期子主题出版物的趋势。

用正则化线性模型& Xgboost 对价格建模

原文：https://towardsdatascience.com/modeling-price-with-regularized-linear-model-xgboost-55e59eae4482?source=collection_archive---------5-----------------------

Photo Credit: Pixabay

开发用于预测个人房价的统计模型

我们想模拟房子的价格，我们知道价格取决于房子的位置、房子的面积、建造年份、翻新年份、卧室数量、车库数量等。因此，这些因素促成了这种模式——优质的位置通常会导致更高的价格。然而，在相同的面积和相同的平方英尺内，所有的房子并没有完全相同的价格。价格的变化就是噪音。我们在价格建模中的目标是对模式建模，忽略噪声。同样的概念也适用于建模酒店房间价格。

因此，首先，我们将为房价数据的线性回归实现正则化技术。

数据

有一个优秀的房价数据集可以在这里找到。

import warnings
def ignore_warn(*args, **kwargs):
    pass
warnings.warn = ignore_warn
import numpy as np 
import pandas as pd 
%matplotlib inline
import matplotlib.pyplot as plt 
import seaborn as sns
from scipy import stats
from scipy.stats import norm, skew
from sklearn import preprocessing
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.linear_model import ElasticNetCV, ElasticNet
from xgboost import XGBRegressor, plot_importance 
from sklearn.model_selection import RandomizedSearchCV
from sklearn.model_selection import StratifiedKFold
pd.set_option('display.float_format', lambda x: '{:.3f}'.format(x))
df = pd.read_csv('house_train.csv')
df.shape

Figure 1

(df.isnull().sum() / len(df)).sort_values(ascending=False)[:20]

Figure 2

好消息是我们有很多特性可以玩(81)，坏消息是 19 个特性有缺失值，其中 4 个缺失值超过 80%。对于任何特性，如果它丢失了 80%的值，那么它就不会那么重要，因此，我决定删除这 4 个特性。

df.drop(['PoolQC', 'MiscFeature', 'Alley', 'Fence', 'Id'], axis=1, inplace=True)

探索功能

目标特征分布

sns.distplot(df['SalePrice'] , fit=norm);# Get the fitted parameters used by the function
(mu, sigma) = norm.fit(df['SalePrice'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))#Now plot the distribution
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
            loc='best')
plt.ylabel('Frequency')
plt.title('Sale Price distribution')#Get also the QQ-plot
fig = plt.figure()
res = stats.probplot(df['SalePrice'], plot=plt)
plt.show();

Figure 3

目标特征 SalePrice 是右偏的。由于线性模型喜欢正态分布的数据，我们将对 SalePrice 进行转换，使其更符合正态分布。

sns.distplot(np.log1p(df['SalePrice']) , fit=norm);# Get the fitted parameters used by the function
(mu, sigma) = norm.fit(np.log1p(df['SalePrice']))
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))#Now plot the distribution
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
            loc='best')
plt.ylabel('Frequency')
plt.title('log(Sale Price+1) distribution')#Get also the QQ-plot
fig = plt.figure()
res = stats.probplot(np.log1p(df['SalePrice']), plot=plt)
plt.show();

Figure 4

数字特征之间的相关性

pd.set_option('precision',2)
plt.figure(figsize=(10, 8))
sns.heatmap(df.drop(['SalePrice'],axis=1).corr(), square=True)
plt.suptitle("Pearson Correlation Heatmap")
plt.show();

Figure 5

一些特征之间存在很强的相关性。比如 GarageYrBlt 和 YearBuilt，TotRmsAbvGrd 和 GrLivArea，GarageArea 和 GarageCars 都是强相关的。它们实际上或多或少表达了相同的东西。我以后会让elastic netcv帮助减少冗余。

销售价格和其他数字特征之间的相关性

corr_with_sale_price = df.corr()["SalePrice"].sort_values(ascending=False)
plt.figure(figsize=(14,6))
corr_with_sale_price.drop("SalePrice").plot.bar()
plt.show();

Figure 6

销售价格与总体质量的相关性最大(约 0.8)。此外，GrLivArea 的相关性超过 0.7，GarageCars 的相关性超过 0.6。让我们更详细地看一下这 4 个特性。

sns.pairplot(df[['SalePrice', 'OverallQual', 'GrLivArea', 'GarageCars']])
plt.show();

Figure 7

特征工程

• 具有高度偏斜分布(偏斜> 0.75)的对数变换要素
• 虚拟编码分类特征
• 用该列的平均值填写 NaN。
• 训练集和测试集分离。

feature_engineering_price.py

弹力网

• 岭和套索回归是正则化的线性回归模型。
• ElasticNe t 本质上是一个套索/脊混合体，它需要一个目标函数的最小化，这个目标函数包括 L1 (套索)和 L2 (脊)规范。
• 当存在多个相互关联的特征时，ElasticNet 非常有用。
• 类 ElasticNetCV 可用于通过交叉验证设置参数alpha (α)和l1_ratio (ρ)。
• elastic net cv:elastic net模型，通过交叉验证选择最佳模型。

让我们看看elastic tcv将为我们选择什么。

ElasticNetCV.py

Figure 8

0 < The optimal l1_ratio <1 , indicating the penalty is a combination of L1 and L2, that is, the combination of 拉索和脊。

模型评估

ElasticNetCV_evaluation.py

Figure 9

这里的 RMSE 其实就是 RMSLE (均方根对数误差)。因为我们取了实际值的对数。这是一篇很好的文章，解释了 RMSE 和 RMSLE 的不同之处。

特征重要性

feature_importance = pd.Series(index = X_train.columns, data = np.abs(cv_model.coef_))n_selected_features = (feature_importance>0).sum()
print('{0:d} features, reduction of {1:2.2f}%'.format(
    n_selected_features,(1-n_selected_features/len(feature_importance))*100))feature_importance.sort_values().tail(30).plot(kind = 'bar', figsize = (12,5));

Figure 10

减少 58.91%的功能看起来很有效。ElasticNetCV 选出的前 4 个最重要的特性分别是 Condition2_PosN 、 MSZoning_C(all) 、exterior 1st _ brk comm&GrLivArea。我们将看看这些特性与 Xgboost 选择的特性相比如何。

Xgboost

第一个 Xgboost 模型，我们从默认参数开始。

xgb_model1.py

Figure 11

它已经比 ElasticNetCV 选择的型号好得多了！

第二个 Xgboost 模型，我们逐渐添加一些参数，以增加模型的准确性。

xgb_model2.py

Figure 12

又有进步了！

第三个 Xgboost 模型，我们添加了一个学习率，希望它会产生一个更准确的模型。

xgb_model3.py

Figure 13

可惜没有什么起色。我的结论是 xgb_model2 是最好的型号。

特征重要性

from collections import OrderedDict
OrderedDict(sorted(xgb_model2.get_booster().get_fscore().items(), key=lambda t: t[1], reverse=True))

Figure 14

Xgboost 选择的前 4 个最重要的特性是 LotArea ， GrLivArea ，overall qual&TotalBsmtSF。

只有一个特性 GrLivArea 被 ElasticNetCV 和 Xgboost 同时选中。

因此，现在我们将选择一些相关的功能，并再次适应 Xgboost。

xgb_model5.py

Figure 15

又一个小改进！

Jupyter 笔记本可以在 Github 上找到。享受这周剩下的时光吧！

基于变分自动编码器的电信客户流失建模

原文：https://towardsdatascience.com/modeling-telecom-customer-churn-with-variational-autoencoder-4e5cf6194871?source=collection_archive---------10-----------------------

Photo credit: Pixabay

如何应用深度卷积神经网络和自动编码器建立流失预测模型。

An autoencoder 是深度学习对降维的回答。这个想法非常简单:通过一系列隐藏层来转换输入，但要确保最终的输出层与输入层具有相同的维度。然而，插入的隐藏层具有逐渐减少的节点数量(因此，减少了输入矩阵的维数)。如果输出与输入紧密匹配或编码，则最小隐藏层的节点可以作为有效的降维数据集。

变分自动编码器 (VAE)类似于经典自动编码器，是由编码器、解码器和损失函数组成的神经网络。它们让我们设计复杂的数据生成模型，并使它们适合大型数据集。

在阅读了一篇关于使用卷积网络和自动编码器洞察用户流失的文章后。我决定对电信客户流失数据集实现 VAE，该数据集可以从 IBM 样本数据集下载。将 VAE 应用于这样一个相对较小的数据集有点矫枉过正，但为了学习 VAE，我还是要这么做。

数据

每行代表一个客户，每列包含列元数据上描述的客户属性。

数据集包括以下信息:

• 上个月内离开的客户—这一列称为流失
• 每位客户已注册的服务—电话、多条线路、互联网、在线安全、在线备份、设备保护、技术支持以及流媒体电视和电影
• 客户账户信息——他们成为客户的时间、合同、支付方式、无纸化账单、每月费用和总费用
• 客户的人口统计信息—性别、年龄范围，以及他们是否有伴侣和家属。

import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler,MinMaxScaler
import collections
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import preprocessing
from sklearn.metrics import (confusion_matrix, precision_recall_curve, auc,
                             roc_curve, recall_score, classification_report, f1_score,
                             precision_recall_fscore_support)
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler,MinMaxScaler
from keras.layers import Input, Dense, Lambda
from keras.models import Model
from keras.objectives import binary_crossentropy
from keras.callbacks import LearningRateScheduler
from keras.utils.vis_utils import model_to_dot
from keras.callbacks import EarlyStopping, ModelCheckpoint
import keras.backend as K
from keras.callbacks import Callback
import matplotlib
matplotlib.rcParams['figure.figsize'] = (10.0, 6.0)df = pd.read_csv('WA_Fn-UseC_-Telco-Customer-Churn.csv')
df.info()

Figure 1

总费用应转换成数字。

df['TotalCharges'] = pd.to_numeric(df['TotalCharges'], errors='coerce')

数据集中的大多数特征都是分类的。我们将首先可视化它们，然后创建虚拟变量。

可视化并分析分类特征

性别

gender_plot = df.groupby(['gender', 'Churn']).size().reset_index().pivot(columns='Churn', index='gender', values=0)
gender_plot.plot(x=gender_plot.index, kind='bar', stacked=True);
print('Gender', collections.Counter(df['gender']))

Figure 2

性别似乎对客户流失没有影响。

合作伙伴

partner_plot = df.groupby(['Partner', 'Churn']).size().reset_index().pivot(columns='Churn', index='Partner', values=0)
partner_plot.plot(x=partner_plot.index, kind='bar', stacked=True);
print('Partner', collections.Counter(df['Partner']))

Figure 3

客户是否有合作伙伴似乎对客户流失有一些影响。

受赡养者

dependents_plot = df.groupby(['Dependents', 'Churn']).size().reset_index().pivot(columns='Churn', index='Dependents', values=0)
dependents_plot.plot(x=dependents_plot.index, kind='bar', stacked=True);
print('Dependents', collections.Counter(df['Dependents']))

Figure 4

没有家属的客户比有家属的客户更容易流失。

电话服务

phoneservice_plot = df.groupby(['PhoneService', 'Churn']).size().reset_index().pivot(columns='Churn', index='PhoneService', values=0)
phoneservice_plot.plot(x=phoneservice_plot.index, kind='bar', stacked=True);
print('PhoneService', collections.Counter(df['PhoneService']))

Figure 5

没有多少客户没有注册电话服务，客户是否有电话服务似乎对客户流失没有影响。

多线

multiplelines_plot = df.groupby(['MultipleLines', 'Churn']).size().reset_index().pivot(columns='Churn', index='MultipleLines', values=0)
multiplelines_plot.plot(x=multiplelines_plot.index, kind='bar', stacked=True);
print('MultipleLines', collections.Counter(df['MultipleLines']))

Figure 6

客户是否签署了多线似乎对客户流失没有影响。

互联网服务

internetservice_plot = df.groupby(['InternetService', 'Churn']).size().reset_index().pivot(columns='Churn', index='InternetService', values=0)
internetservice_plot.plot(x=internetservice_plot.index, kind='bar', stacked=True);
print('InternetService', collections.Counter(df['InternetService']))

Figure 7

似乎注册光纤的客户最有可能流失，几乎 50%的客户流失了。

在线安全

onlinesecurity_plot = df.groupby(['OnlineSecurity', 'Churn']).size().reset_index().pivot(columns='Churn', index='OnlineSecurity', values=0)
onlinesecurity_plot.plot(x=onlinesecurity_plot.index, kind='bar', stacked=True);
print('OnlineSecurity', collections.Counter(df['OnlineSecurity']))

Figure 8

没有注册 OnlineSecurity 的客户最有可能流失。

在线备份

onlinebackup_plot = df.groupby(['OnlineBackup', 'Churn']).size().reset_index().pivot(columns='Churn', index='OnlineBackup', values=0)
onlinebackup_plot.plot(x=onlinebackup_plot.index, kind='bar', stacked=True);
print('OnlineBackup', collections.Counter(df['OnlineBackup']))

Figure 9

没有注册 OnlineBackUp 的客户最有可能流失。

设备保护

deviceprotection_plot = df.groupby(['DeviceProtection', 'Churn']).size().reset_index().pivot(columns='Churn', index='DeviceProtection', values=0)
deviceprotection_plot.plot(x=deviceprotection_plot.index, kind='bar', stacked=True);
print('DeviceProtection', collections.Counter(df['DeviceProtection']))

Figure 10

没有注册 DeviceProtection 的客户最有可能流失。

技术支持

techsupport_plot = df.groupby(['TechSupport', 'Churn']).size().reset_index().pivot(columns='Churn', index='TechSupport', values=0)
techsupport_plot.plot(x=techsupport_plot.index, kind='bar', stacked=True);
print('TechSupport', collections.Counter(df['TechSupport']))

Figure 11

没有注册技术支持的客户最有可能流失。

流媒体电视

streamingtv_plot = df.groupby(['StreamingTV', 'Churn']).size().reset_index().pivot(columns='Churn', index='StreamingTV', values=0)
streamingtv_plot.plot(x=streamingtv_plot.index, kind='bar', stacked=True);
print('StreamingTV', collections.Counter(df['StreamingTV']))

Figure 12

流媒体电影

streamingmovies_plot = df.groupby(['StreamingMovies', 'Churn']).size().reset_index().pivot(columns='Churn', index='StreamingMovies', values=0)
streamingmovies_plot.plot(x=streamingmovies_plot.index, kind='bar', stacked=True);
print('StreamingMovies', collections.Counter(df['StreamingMovies']))

Figure 13

从以上七个图中，我们可以看到，没有互联网服务的客户流失率非常低。

合同

contract_plot = df.groupby(['Contract', 'Churn']).size().reset_index().pivot(columns='Churn', index='Contract', values=0)
contract_plot.plot(x=contract_plot.index, kind='bar', stacked=True);
print('Contract', collections.Counter(df['Contract']))

Figure 14

很明显，合同条款确实对客户流失有影响。当客户有一份两年的合同时，很少有不愉快的事情发生。大多数客户都有逐月合约。

无纸账单

paperlessbilling_plot = df.groupby(['PaperlessBilling', 'Churn']).size().reset_index().pivot(columns='Churn', index='PaperlessBilling', values=0)
paperlessbilling_plot.plot(x=paperlessbilling_plot.index, kind='bar', stacked=True);
print('PaperlessBilling', collections.Counter(df['PaperlessBilling']))

Figure 15

支付方式

paymentmethod_plot = df.groupby(['PaymentMethod', 'Churn']).size().reset_index().pivot(columns='Churn', index='PaymentMethod', values=0)
paymentmethod_plot.plot(x=paymentmethod_plot.index, kind='bar', stacked=True);
print('PaymentMethod', collections.Counter(df['PaymentMethod']))

Figure 16

支付方式似乎对流失率有影响，尤其是电子支票支付的流失率最高。

老年人

seniorcitizen_plot = df.groupby(['SeniorCitizen', 'Churn']).size().reset_index().pivot(columns='Churn', index='SeniorCitizen', values=0)
seniorcitizen_plot.plot(x=seniorcitizen_plot.index, kind='bar', stacked=True);
print('SeniorCitizen', collections.Counter(df['SeniorCitizen']))

Figure 17

我们的数据中没有很多老年人。看起来顾客是否是老年人对客户流失率没有影响。

探索数字特征

任期

sns.kdeplot(df['tenure'].loc[df['Churn'] == 'No'], label='not churn', shade=True);
sns.kdeplot(df['tenure'].loc[df['Churn'] == 'Yes'], label='churn', shade=True);

Figure 18

df['tenure'].loc[df['Churn'] == 'No'].describe()

Figure 19

df['tenure'].loc[df['Churn'] == 'Yes'].describe()

Figure 20

没有流失的客户比流失的客户拥有更长的平均任期(20 个月)。

月费

sns.kdeplot(df['MonthlyCharges'].loc[df['Churn'] == 'No'], label='not churn', shade=True);
sns.kdeplot(df['MonthlyCharges'].loc[df['Churn'] == 'Yes'], label='churn', shade=True);

Figure 21

df['MonthlyCharges'].loc[df['Churn'] == 'No'].describe()

Figure 22

df['MonthlyCharges'].loc[df['Churn'] == 'Yes'].describe()

Figure 23

流失客户支付的月平均费用比非流失客户高 20%以上。

总费用

sns.kdeplot(df['TotalCharges'].loc[df['Churn'] == 'No'], label='not churn', shade=True);
sns.kdeplot(df['TotalCharges'].loc[df['Churn'] == 'Yes'], label='churn', shade=True);

Figure 24

数据预处理

用 0 到 1 之间的值对标签进行编码。

le = preprocessing.LabelEncoder()
df['Churn'] = le.fit_transform(df.Churn.values)

用该列的平均值填写 nan。

df = df.fillna(df.mean())

编码分类特征。

categorical =  ['gender', 'Partner', 'Dependents', 'PhoneService', 'MultipleLines', 'InternetService', 'OnlineSecurity', 'OnlineBackup', 'DeviceProtection', 'TechSupport', 'StreamingTV', 'StreamingMovies', 'Contract', 'PaperlessBilling', 'PaymentMethod']
for f in categorical:
    dummies = pd.get_dummies(df[f], prefix = f, prefix_sep = '_')
    df = pd.concat([df, dummies], axis = 1)
# drop original categorical features
df.drop(categorical, axis = 1, inplace = True)

将数据分为训练集、验证集和测试集，并创建批处理通过我们的网络发送。

autoencoder_preprocessing.py

VAE 在喀拉斯的实施

以下代码脚本主要来自 Agustinus Kristiadi 的博客文章:变化的自动编码器:直觉和实现。

• 定义输入层。
• 定义编码器层。
• 编码器模型，将输入编码为潜在变量。
• 我们使用平均值作为输出，因为它是中心点，代表高斯分布。
• 我们从 2 个密集层的输出中取样。
• 定义 VAE 模型中的解码层。
• 定义总体 VAE 模型，用于重建和培训。
• 定义生成器模型，生成给定潜在变量 z 的新数据。
• 将我们的损失转化为 Keras 代码。
• 开始训练。

VAE.py

该模型在批量为 100 个样本的 55 个时期后停止训练。

估价

plt.plot(vae_history.history['loss'])
plt.plot(vae_history.history['val_loss'])
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='upper left')
plt.show();

Figure 25

从上面的损失图中，我们可以看到，该模型在训练和验证数据集上都具有可比的性能，并且在最后似乎收敛得很好。

我们使用重构误差来衡量解码器的性能。训练自动编码器以减少重建误差，如下所示:

x_train_encoded = encoder.predict(X_train)pred_train = decoder.predict(x_train_encoded)
mse = np.mean(np.power(X_train - pred_train, 2), axis=1)
error_df = pd.DataFrame({'recon_error': mse,
                        'churn': y_train})plt.figure(figsize=(10,6))
sns.kdeplot(error_df.recon_error[error_df.churn==0], label='not churn', shade=True, clip=(0,10))
sns.kdeplot(error_df.recon_error[error_df.churn==1], label='churn', shade=True, clip=(0,10))
plt.xlabel('reconstruction error');
plt.title('Reconstruction error - Train set');

Figure 26

x_val_encoded = encoder.predict(X_val)pred = decoder.predict(x_val_encoded)
mseV = np.mean(np.power(X_val - pred, 2), axis=1)
error_df = pd.DataFrame({'recon_error': mseV,
                        'churn': y_val})plt.figure(figsize=(10,6))
sns.kdeplot(error_df.recon_error[error_df.churn==0], label='not churn', shade=True, clip=(0,10))
sns.kdeplot(error_df.recon_error[error_df.churn==1], label='churn', shade=True, clip=(0,10))
plt.xlabel('reconstruction error');
plt.title('Reconstruction error - Validation set');

Figure 27

潜在空间可视化

我们可以在 2D 潜在空间中聚类客户，并可视化搅动和未搅动的客户，它们可以在潜在空间中分离，并揭示不同聚类的形成。

x_train_encoded = encoder.predict(X_train)plt.scatter(x_train_encoded[:, 0], x_train_encoded[:, 1], 
            c=y_train, alpha=0.6)
plt.title('Train set in latent space')
plt.show();

Figure 28

x_val_encoded = encoder.predict(X_val)plt.scatter(x_val_encoded[:, 0], x_val_encoded[:, 1], 
            c=y_val, alpha=0.6)plt.title('Validation set in latent space')
plt.show();

Figure 29

对验证集的预测

x_val_encoded = encoder.predict(X_val)
fpr, tpr, thresholds = roc_curve(y_val, clf.predict(x_val_encoded))
roc_auc = auc(fpr, tpr)plt.title('Receiver Operating Characteristic')
plt.plot(fpr, tpr, label='AUC = %0.4f'% roc_auc)
plt.legend(loc='lower right')
plt.plot([0,1],[0,1],'r--')
plt.xlim([-0.001, 1])
plt.ylim([0, 1.001])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.show();

Figure 29

print('Accuracy:')
print(accuracy_score(y_val, clf.predict(x_val_encoded)))
print("Confusion Matrix:")
print(confusion_matrix(y_val,clf.predict(x_val_encoded)))
print("Classification Report:")
print(classification_report(y_val,clf.predict(x_val_encoded)))

Figure 30

对测试集的预测

x_test_encoded = encoder.predict(X_test)
fpr, tpr, thresholds = roc_curve(y_test, clf.predict(x_test_encoded))
roc_auc = auc(fpr, tpr)
plt.title('Receiver Operating Characteristic')
plt.plot(fpr, tpr, label='AUC = %0.4f'% roc_auc)
plt.legend(loc='lower right')
plt.plot([0,1],[0,1],'r--')
plt.xlim([-0.001, 1])
plt.ylim([0, 1.001])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.show();

Figure 31

print('Accuracy:')
print(accuracy_score(y_test, clf.predict(x_test_encoded)))
print("Confusion Matrix:")
print(confusion_matrix(y_test,clf.predict(x_test_encoded)))
print("Classification Report:")
print(classification_report(y_test,clf.predict(x_test_encoded)))

Figure 32

就是这样！Jupyter 笔记本可以在 Github 上找到。周一快乐！

参考资料:

[## 可变自动编码器。

在我关于自动编码器的介绍性文章中，我讨论了各种模型(欠完整、稀疏、去噪、收缩)…

www.jeremyjordan.me](https://www.jeremyjordan.me/variational-autoencoders/) [## 教程-什么是变分自动编码器？-贾恩·阿尔托萨尔

从两个角度理解变分自动编码器:深度学习和图形模型。

jaan.io](https://jaan.io/what-is-variational-autoencoder-vae-tutorial/) [## 解释了各种自动编码器

在我上一篇关于生成对立网络的文章中，我介绍了一个简单的方法来训练一个可以…

kvfrans.com](http://kvfrans.com/variational-autoencoders-explained/) [## 使用 Keras 中的自动编码器检测信用卡欺诈——面向黑客的 TensorFlow(第七部分)

信用卡交易中的异常检测是如何工作的？

medium.com](https://medium.com/@curiousily/credit-card-fraud-detection-using-autoencoders-in-keras-tensorflow-for-hackers-part-vii-20e0c85301bd) [## 基本自动编码器-使用重构误差进行异常检测| Deeplearning4j

下载此笔记本请查看自述文件，了解安装、设置依赖关系和导入…

deeplearning4j.org](https://deeplearning4j.org/tutorials/05-basic-autoencoder-anomaly-detection-using-reconstruction-error) [## Naomi fridman/深度 VAE 客户流失预测

预测客户流失的可变深度自动编码器-Naomi fridman/深度 VAE-客户流失预测

github.co](https://github.com/naomifridman/Deep-VAE-prediction-of-churn-customer)