使用神经网络(AI)预测西雅图警察拦截期间逮捕和搜身的可能性
面对 2008-2009 年大衰退期间全球经济的不确定性，技术世界继续以惊人的速度前进。苹果推出了 iPhone 3G，这是他们的第二款 iPhone 产品，也是第一款包含现在广受欢迎的 App Store 的产品，它将在 2020 年超过每周$1B 的收入，并独自超越世界上大多数国家的国内生产总值。几年前收购安卓之后，谷歌推出了安卓操作系统，后来占据了手机市场近 75%的市场份额。伴随这些重大技术的是“物联网”的形成，这是全球数十亿设备的集合，这些设备现在连接到互联网，并持续收集和共享数据。物联网诞生于 2008 年或 2009 年，当时联网设备的数量超过了联网人数。
自 2008 年以来，物联网(IoT)的存在在当今世界的各个方面都变得更加明显。从这些物联网设备中创建的大量数据是巨大的 IBM 在 2017 年估计，曾经存在的 90%的数据是在 2017 年之前的两年中创建的！我们每天创建 2.5 万亿字节的数据，随着技术的进步和数字鸿沟的缩小，这个数字还在不断增长。许多组织已经学会了利用他们的数据来快速发展他们的策略，以满足他们的客户的经常不明确的需求。这样做为这些公司维持了巨大的增长和盈利水平——例如，脸书、Alphabet、亚马逊、网飞和沃尔玛，为他们自己的市场价值增加了数十亿美元，并使自己成为世界上最有价值的公司之一。
物联网凸显了数据在当今数字化社会中的重要性。没有采用物联网传感器的组织至少正在改变他们的标准操作程序，以从常规业务流程中捕获数据。由于大规模数据管理和仓储的日益复杂性，以及缺乏这样做的短期激励，许多政府组织在采用和促进这一新的数字时代方面行动迟缓。西雅图的亚马逊和雷蒙德的微软总部所在地金县是一个反驳这种观点的地方政府。金县有一个开放的数据政策，定期将最新数据发布到他们位于 https://data.kingcounty.gov/的公民门户网站。数据主题包括执法、环境质量、交通等。
我使用人工智能(AI)预测西雅图警察局在“特里”拦截(临时拘留或涉嫌犯罪活动的主体)期间进行逮捕或搜身的可能性。西雅图市提供了一个非常棒的数据集，我们将用来训练一个机器学习模型:https://data . Seattle . gov/Public-Safety/Terry-Stops/28ny-9ts 8/data。
数据转换
根据对数据集的检查，并不是该属性中提供的每个属性都有助于进行这些预测。例如，我们可以自信地假设主体 ID 和 Terry Stop ID 号是任意值，与逮捕或搜查主体的决定无关。也就是说，我们将用来完成这项任务的数据集的属性包括在下面。
我们将用来完成这项任务的人工智能概念是一个神经网络。神经网络是受构成动物大脑的生物神经元启发并模仿的计算系统。神经网络接受数字阵列作为输入，然后对这些输入中的每一个应用加权来计算输出。通过使用“反向传播”,神经网络将在训练过程中“学习”,并调整其内部权重，使其越来越接近所需的输出。我将特别使用一个神经网络系统，这个系统是我用微软的编程语言 C#从头开始创建的。NET 框架。我的神经网络使用双曲正切(TANH)激活函数，我有一个具有 50 个隐藏神经元的隐藏层，每层拥有一个偏置神经元，我使用 0.018 的学习率进行反向传播。
神经网络本质上是复杂的矩阵计算器。计算器只接受数值，不接受数据集中存在的字符串(文本)值。这意味着我们需要对数据进行一些翻译，以便计算机能够理解输入。除了军官的出生年份之外，该数据集中提供的所有属性都需要从字符串值转换为数字值。
那么我们如何将这些字符串值转换成数字呢？对于每个属性中的每个唯一值，我们将创建另一个属性，如果原始属性与该唯一值匹配，则包含值 1，如果不匹配，则包含值 0。例如，如果我们要对主题感知种族属性进行转换，并且该属性的唯一值只有“白人”、“黑人”和“西班牙裔”，我们将创建新属性“主题 _ 是 _ 白人”、“主题 _ 是 _ 黑人”和“主题 _ 是 _ 西班牙裔”。对于每条记录，如果受试者感知的种族与这些新属性相匹配，我们将为每条新属性填入 1，如果不匹配，则填入 0。
我将这种转换应用于数据集中的许多属性，以创建新的输入列表，我们将输入到我们的神经网络中:
对于只包含两个唯一值(即性别)的属性，我们不必应用这种转换。对于这些属性，我只是将其中一个值替换为 0，将另一个值替换为 1。对于性别，我将所有“男性”值替换为 0，将所有“女性”值替换为 1。
在我们准备好训练我们的神经网络之前，还有一个转换必须应用。军官年龄字段包含范围从 22 到 69 的值。神经网络将一个激活函数应用到它们的每个内部计算中，以模拟生物神经元的阈值“触发”。TANH(以及许多常用的激活函数)的美妙之处在于，输出将始终在–1 和 1 的范围内，防止内部权重在训练时膨胀到无穷大。
双曲正切函数是非线性的，这意味着输出的梯度(导数)不一致。这种非线性赋予每个神经元**“激发或不激发”的特性，这在动物的生物神经元中是可以找到的。此外，激活函数的非线性使得我们的神经网络能够发现数据中的非线性趋势。每个神经元由动作电位(输入值)刺激，然后通过双曲正切函数传递。如果动作电位足够高，神经元就会“点火”，将正数传递给下一层的每一个神经元。这个输出然后被连接神经元和下一层神经元的轴突放大或沉默**。
由于双曲正切函数的梯度在–1 和 1 之间最为显著(如上图中的紫色框所示)，我们将对数据集的连续属性应用最小-最大归一化，以将其向下“归一化”(拉伸)至–1 和 1 之间。我们需要应用的唯一属性是“军官年龄”属性。
在对上面详细描述的数据进行转换之后，我们现在准备好训练我们的神经网络。我们的数据集包含 13，573 条完整且可用的记录。我决定用 75%的数据来训练我们的模型，剩下的 25%用来测试我们神经网络预测的准确性。在初始化了具有 38 个输入节点、50 个隐藏节点和 2 个输出节点的新神经网络之后，我们准备好训练我们的模型。像这样大规模的神经网络训练可能是计算密集型的，需要一些时间，所以我将工作负载外包给了微软 Azure 的计算优化 FSV2 虚拟机的一个实例，该虚拟机由英特尔的至强白金处理器提供支持。
下面是我们正在使用的神经网络图。原谅我的小文本，但它必须小，以适应。
10，180 次训练迭代和 38 小时后，我们有了一个完全训练好的神经网络，可以做预测了。大约有 3393 个预测，所以我不能在这里全部显示出来。但是，以下是对结果的高级分析。
预测精度
神经网络正确预测逮捕了吗？
神经网络计算了每个事件记录中发生逮捕的概率。对于网络预测最有可能逮捕的 100 条记录，在这 100 起事件中有 47 起真的发生了逮捕。相反，网络预测的 100 条记录中最不可能被逮捕，事实上，其中 73 条没有导致逮捕。结合这两个指标，这意味着在前 200 个最极端的案例中，预测正确 120 次，或成功率为 60%。比扔硬币要好，但还不足以非常自信！这是否意味着我们的神经网络失败了？不一定。如果我们后退一步，我们可以认识到输入的数据可能与逮捕的机会没有很大的相关性。我敢肯定，有更多的逮捕发生的决定因素没有在这个数据集中进行测量——对象的举止，药物和酒精的影响，交谈中的尊重，以及以前的犯罪历史都是警察在拦截中可能的决定因素，没有作为输入纳入我们的模型。因此，由于数据的性质，我们可能要求我们的神经网络执行一项不可能的任务。
神经网络是否正确预测了搜身？
与逮捕类似，我们的神经网络也计算了每起事件发生搜身的概率。对于网络预测最有可能被搜身的 100 条记录，实际上有 82 条被搜身了。对于网络预测最不可能被逮捕的 100 条记录来说，事实上有 87 条记录没有被搜身，而不是被 T2 搜身。总的来说，这意味着我们的模型在 200 个案例中有 169 个是正确的，或者说成功率大约为 85%。这相对于停搏预测的准确性是一个很大的改进，并且足够高，让我相信我们的模型在数据中发现了很强的相关性。为什么不是 100%呢？请记住，这是不完善和不完整的信息。这并没有(也可能永远不会)完全捕捉到在每个警察停车时发生的每一件事。因此，期望从我们提供的有限数据中获得近乎完美的预测精度是不合理的。理论上，随着我们提供更多相关数据，这种准确性应该会提高。
现在将我们的范围改为整个测试数据集，神经网络在 66%的情况下正确预测了逮捕，在 75%的情况下正确预测了搜查。虽然这些结果并不完美，但它们确实表明我们的神经网络能够从数据的相关性中提取和学习。那么神经网络是怎么做到的呢？让我们进一步分析结果，尝试阐明相关性。
属性关联
武器出现
武器的存在如何影响逮捕或搜身的预测？在我们的网络最有信心逮捕的前 100 个案例中，有 70 个案例确实有武器。在我们的网络最有信心进行搜身的前 100 个场景中，武器出现了 94 次。另一方面，在最不可能被逮捕或搜身的 100 个案例中，没有一个案例带有武器。这非常表明，这种输入是神经网络预测结果的强有力指标。
最有可能搜身时携带武器的案件有 94 起，最有可能逮捕时携带武器的案件有 70 起，这一事实表明，武器的存在与搜身而不是逮捕更有关联。这可能与常识相反，但实际上可能是有道理的——也许如果有武器，搜身几乎总是会发生，但考虑到武器可能是合法拥有的，并不总是导致逮捕。
官龄
军官的年龄如何影响我们的预测？以下是结果:
最有可能被捕的平均年龄是 29.3 岁，而最不可能被捕的平均年龄是 38.8 岁。这是一个巨大的差异！这意味着我们的神经网络已经发现了警官年龄和被捕几率之间的相关性，并表明年轻警官比年长警官更有可能实施逮捕。也许年纪较大的警官成熟到不需要进行逮捕，或者积极地试图完全避免这种情况。也许年纪较大的官员通常更多地担任管理角色，因此不太可能实施逮捕。
最有可能搜身和最不可能搜身的平均年龄非常相似——29.9 岁和 28.2 岁，所以不足以进行深入研究和推断。
军官性别
军官的性别如何影响我们的预测？
这里没有大的区别，但有些值得注意。如果警察的性别是男性而不是女性，被逮捕和搜身的可能性更大。尤其是考虑到被捕的可能性。
主题性别
受试者的性别如何影响我们的预测？
我们的神经网络已经注意到了这种明显的相关性。无论是最有可能逮捕的案件还是最有可能搜身的案件，研究对象都是 100%的男性。当观察最不可能逮捕和最不可能搜身的案件时，这一比率变化很大。最不可能被逮捕的案件中只有 77%是男性(23%是女性)，而最不可能被搜身的案件中，男性比例骤降至 28%(72%是女性)。
最有可能搜身和最不可能搜身的男女比例发生了巨大变化。最有可能的情况下是 100%男性，最不可能的情况下是 72%女性！为什么会这样？我们可以做出许多假设，并用数据来检验每一个假设，但我的最佳猜测是，这是因为由于潜在的性骚扰反弹或诉讼，警察不太可能对女性进行搜身。
这是否意味着官员在工作中有意识或无意识的偏见？也许吧！但是匆忙得出这样的结论是危险的…这可能是许多事情的结果。从理论上讲，这也可能是男性从事的活动的违法程度比女性平均水平更严重的结果。
学科年龄
受试者的年龄如何影响我们的预测？让我们来看看结果:
最有可能被捕:
最有可能搜身的:
最不可能被逮捕:
最不可能搜身的:
请注意，在被标记为最有可能被捕的前 100 个案例中，有 66 个案例的年龄在 26 至 35 岁之间。这意味着我们的神经网络表明这个年龄范围特别容易被逮捕。再问一次，这是否意味着西雅图的警察对这个年龄范围有偏见？也许，但不一定。这可能是该年龄段的人从事的非法活动比其他年龄段的人多得多的结果。
另一方面，最不可能被逮捕的人有更多的年轻和年长的人，这表明这些年龄的人确实不太可能被逮捕。类似的趋势也适用于弗里斯克的预测。
军官竞赛
军官的种族作为一个属性如何影响我们模型的预测？
最有可能被捕:
最有可能被搜身的:
最不可能被逮捕:
最不可能搜身:
这里的比率变化不大，不能提供有用的信息来推断模式或得出结论。
学科竞赛
受试者的种族如何影响我们的预测？
最有可能被捕:
最有可能搜身的:
最不可能被逮捕:
最不可能搜身的:
与警官的种族相似，受试者的种族似乎不是逮捕或搜身结果的有影响的决定因素，但也许在某些地区有一些小的关联。然而，请注意最有可能搜身的案例和最不可能搜身的案例之间的组合变化。最不可能搜身的案件非常严重地偏向黑人。为什么警察比其他种族的警察更少对黑人搜身？没有办法用我们现有的数据来判断，但神经网络现在已经警告我们这一趋势——如果我们有这些资源可供支配，我们可以选择用其他方法和数据来进一步挖掘。
最有可能和最不可能的逮捕组合非常相似。这是不是意味着我们的神经网络坏了？不，这表明这一指标与逮捕或搜身的可能性之间没有很大的关联。
星期几
星期几如何影响预测？
最有可能被逮捕:
最有可能被搜身的是:
最不可能被逮捕:
最不可能搜身的是
除了周二被捕的可能性最小之外，这里似乎没有任何值得注意的趋势。
一天中的时间
我们把一天的时间分成四个不同的组:午夜到早上 6 点，早上 6 点到中午，中午到下午 6 点，下午 6 点到午夜。我们完全可以把它分成 24 小时，甚至更长时间，但是这样效率很低，而且可能不会有什么深刻的发现。训练我们的神经网络需要更多的时间，所以我妥协了，把一天中的时间分成四组。以下是结果:
最有可能被捕:
最有可能搜身的:
最不可能被捕:
最不可能搜身的:
结果表明了你的预期:逮捕和搜身更有可能发生在深夜或凌晨。
辖区
最后，辖区如何影响我们神经网络的预测？数据集的 terry 停靠点出现在五个不同的区域:北部、东部、南部、西部和西南部。
最有可能被捕:
最有可能搜身的:
最不可能被逮捕:
最不可能搜身的:
这里最清楚的指标在于逮捕数据:请注意，我们的网络预测的最有可能逮捕的前 100 起案件中有 64%发生在西区。进一步重申这一点，在我们的神经网络最有信心逮捕不会不会发生的前 100 条记录中，北部似乎是一个相对安全的区域，84%的记录在北部，而西部没有记录发生！
这让我们清楚地看到西雅图的哪些选区可能是安全的，哪些选区可能会变得更加危险。从神经网络做出的预测中，我们可以推断，与西部或者南部相比，北部和西南部是城镇中相对安全的部分。
但是这个假设有保证吗？绝对不行！也许数据是扭曲的——可能西雅图警察局有更多的警察在某些地区巡逻，比如西部，这就是那里有更多逮捕的原因。然而，这一反驳点可能被证明是错误的，因为北部是可能发生搜身的地区之一，这表明那里确实有警察存在。
总结
如前所述，我们的神经网络正确预测逮捕的准确率为 66%，正确预测搜身的准确率为 75%。虽然这些数字并不完美，但它们比扔硬币要好，表明我们的模型能够发现并学习输入数据与逮捕或搜身可能性之间的相关性。
这对我们有什么用？回想起来，这种学习模型可以帮助警察或政府官员发现数据中的异常和趋势。这些见解可用于帮助规划资源分配。例如，如果神经网络模型发现一种趋势，即逮捕很可能发生在城镇的特定地区，在一周的特定日期，在一天的特定时间，官员可以将更多的警力资源分配到该地区，以在犯罪活动发生之前主动预防犯罪活动。在人工分析的监督下，像这样的多层趋势通常很难发现；多亏了神经网络，我们现在可以揭示这些趋势。
我为处理这个数据集而创建的神经网络花了大约 38 个小时来训练 10，180 个例子，或者每次迭代大约 13.5 秒。训练时间长的原因是计算密集型顺序反向传播算法，该算法调整模型的内部权重以更好地逼近每个输出。幸运的是，正向传播(要求模型做出预测)要快得多。即使比我们在这里使用的模型更深的模型也可以在不到十分之一秒的时间内做出预测，即使是在像 5 美元的 Raspberry Pi Zero 这样的低功耗廉价硬件上。
在一个神经网络被创建并在至少性能良好的硬件上训练之后，该模型可以被部署到低性能硬件上，这样它的预测能力就可以在“边缘”(前线)使用——通过物联网设备在警察手中使用。也许在未来，警察将佩戴一系列小型物联网传感器，这些传感器可以实时捕获比该数据集中列出的更多的数据:对话对话、面部表情图像、与对象的接近程度等等。这些传感器可以通过应用情绪分析(可能也由另一个神经网络驱动)来判断对话的情绪，并将其对对话中紧张程度的估计作为输入传递给我们的模型。计算机视觉技术，如微软的 Azure Face 认知服务，可以分析实时图像，并对军官和受试者的压力水平进行预测，并将这些作为输入传递给我们的模型。你甚至可以给警察配备心率监测器或皮质醇监测设备来检测压力或不适，并以此作为输入。
所有这些输入都可以输入到每个官员佩戴的小型离散物联网设备中。这种设备将使用最新版本的经过训练的神经网络，随着情况的发展，不断预测逮捕或搜身的可能性。如果逮捕的可能性超过了预定的阈值，连接的物联网设备可以自动向该地区的另一名官员发出支持信号。军官们将不再需要自己通过无线电寻求帮助，这可能会通过诱导受试者的“战斗或逃跑”反应来使局势进一步升级，并且还会使他们自己处于这样做的潜在反弹的危险中。
随着物联网传感器在各个经济部门的各个业务领域变得越来越普遍，高质量数据的可用性将会增加。幸运的是，希望参与这场数据驱动的人工智能革命的组织可以通过价格合理的现成技术来实现。这样做将极大地改进学习模型，因为引入了高度相关的数据，增强了模型的预测准确性。我鼓励任何对企业数字化转型感兴趣的组织领导者进行投资——与未来的收益或成本节约相比，今天的成本可能微不足道。
预测网址的恶意性
创建用于分离恶意和良性 URL 的特征向量
在本文中，我们将逐步开发一个用于识别恶意 URL 的简单特性集表示。我们将为 URL 创建特征向量，并使用这些向量来开发用于识别恶意 URL 的分类模型。为了评估这些特征在区分恶意网址和良性网址方面有多好，我们建立了一个基于决策树的机器学习模型来预测给定网址的恶意性。
恶意网站是网络安全中众所周知的威胁。它们是在线传播病毒、蠕虫和其他类型恶意代码的有效工具，并且是大多数网络攻击的 60%以上的罪魁祸首。恶意 URL 可以通过电子邮件链接、文本消息、浏览器弹出窗口、页面广告等方式发送。这些 URL 可能是指向可疑网站的链接，或者很可能嵌入了“可下载内容”。这些嵌入的下载可能是间谍软件、键盘记录器、病毒、蠕虫等。因此，及时检测和缓解恶意代码在其网络中的传播已成为网络防御人员的当务之急。恶意 URL 检测器的各种技术以前主要依赖于 URL 黑名单或签名黑名单。这些技术大多提供“事后”解决方案。为了提高恶意 URL 检测方法的及时性和抽象性，机器学习技术正日益被接受。
为了开发机器学习模型，我们需要一个特征提取框架来表征 URL 或将 URL 转换为特征向量。在本文中，我们将收集已知恶意网址和已知良性网址的样本。然后，我们开发一个指纹识别框架，并为样本中的所有 URL 提取一组给定的***【M】***特征。我们通过开发一个具有这些特征的简单预测模型来测试这些特征在区分恶意 URL 和良性 URL 方面的有效性。最后，我们测量模型预测恶意网址的能力，作为区分恶意网址和良性网址的特征的有效性。
下图是本文中方法流程的概述。
分析处理
数据
我们从两个来源收集数据:Alexa 1000 强网站和 phishtank.com。从 Alexa 的前 1000 个网站中抓取了 1000 个假定的良性 URL，从phishtank.com中抓取了 1000 个可疑的恶意 URL。由于 virustotal API 的限制率，我们随机抽取了 500 个假定的良性 URL 和 500 个假定的恶意 URL。然后通过病毒扫描这些网址。检测到 0 个恶意内容的 URL 被标记为良性(b_urlX)，检测到至少 8 个恶意内容的 URL 被标记为恶意(m_urlX)。我们将每次扫描的 JSON 结果转储到相应的文件‘b _ urlx . JSON’、‘m _ urlx . JSON’中。您可以在这里找到这些文件。
from requests import get
from os import listdir
import pandas **as** pd
import numpy **as** np
from pandas.io.json import json_normalize
import seaborn **as** sns
import matplotlib.pyplot **as** plt
import math
from datetime import datetime
plt**.**rcParams["figure.figsize"] **=** (20,20)
处理应用编程接口速率限制和 IP 阻塞
为了确认示例中的恶意 URL 是恶意的，我们需要向 VirusTotal 发送多个请求。病毒扫描提供多个病毒扫描引擎的汇总结果。此外，我们通过(Shodan)[shodan.io]传递 URL。Shodan 是所有连接到互联网的设备的搜索引擎，提供基于服务的 URL 服务器功能。VirusTotal 和 Shodan 目前的 API 速率限制分别是每分钟 4 个请求和每个 API 密钥每月至少 10，000 个请求。虽然对数据的 URL 请求数量在 Shodan API 请求限制之内，但事实证明，病毒总数有点困难。这是通过创建几个虚拟磁带库应用编程接口密钥(最多 4 个)并在每个请求中随机采样来解决的。除了限制 API 请求的数量之外，在短时间内发送多个请求还会导致来自 VT 和 Shodan 服务器的 IP 阻塞。我们编写了一个小爬虫，从https://free-proxy-list.net/获取最新的精英 IP 地址集，并为每个请求创建一个新的代理列表，因为自由代理的寿命非常短。除了 IP 池之外，我们还使用 Python 的 FakeUserAgent 库在每个请求上切换用户代理。
最后，对于每个请求，使用新的代理和用户代理，我们可以每分钟发送 16 个请求，而不是之前的 4 个。每个请求都有以下请求参数:
- 1 病毒总密钥:来自 VT API 密钥池的样本。
- 1 Shodan 密钥:来自 Shodan API 密钥池的示例。
- 1 IP:向https://free-proxy-list.net/发送请求，获取最新的免费精英代理。
- 1 用户代理:Python 的可用用户代理示例(假冒用户代理)[https://pypi . org/project/假冒用户代理/]
Shodan 和 VT 的扫描产生了以下数据集。从 shodan 中，我们提取了以下特征:
- numServices:主机上运行的服务(开放端口)总数
- robotstxt:网站是否启用了 robots txt
扫描后的最终数据集在这里可用。您可以下载这些数据并运行您的分析。
指纹 URL(用于恶意软件 URL 检测的特征 URL)
目标是提取 URL 特征，这些特征对于区分恶意 URL 和良好 URL 非常重要。首先，让我们看看一个 URL 的结构中的相关部分。
典型的 URL
URL(统一资源定位器的缩写)是一种引用，它指定了 web 资源在计算机网络上的位置以及检索它的机制。如下图所示，URL 由不同的组件组成。该协议或方案规定了如何(或需要什么)传输信息。主机名是计算机网络上计算机 IP 地址的人类可读的唯一引用。域名服务(DNS)命名层次结构将 IP 地址映射到主机名。被攻破的网址被用来实施网络攻击。这些攻击可能是任何或多种形式的网络钓鱼电子邮件、垃圾邮件和驾车下载。
关于域名，所有者购买人们觉得更容易记住的域名。所有者通常希望他们提供的品牌、产品或服务的名称是特定的。URL 的这一部分(域)一旦设置就不能更改。恶意域名所有者可能会选择多个廉价域名，例如“xsertyh.com”。
免费 URL 参数是 URL 的一部分，可以通过更改来创建新的 URL。这些包括目录名、文件路径和 URL 参数。这些免费的 URL 参数通常被攻击者操纵来创建新的 URL、嵌入恶意代码并传播它们。
有许多用于恶意 URL 检测的技术，两种主要技术是 a)黑名单技术，和 b)机器学习技术。黑名单包括维护已知恶意域的数据库，并将新 URL 的主机名与数据库中的主机名进行比较。这是一个“事后”问题。它将无法检测新的和看不见的恶意网址，只有当它被认为是来自受害者的恶意网址时，它才会被加入黑名单。另一方面，机器学习方法提供了一种可跨平台推广且独立于已知签名的先验知识的预测方法。给定恶意和良性恶意软件样本的样本，ML 技术将提取已知好的和坏的 URL 的特征，并概括这些特征以识别新的和看不见的好的或坏的 URL。
URL 指纹识别过程针对 3 种类型的 URL 功能:
- URL 字符串特征:源自 URL 字符串本身的特征。
- URL 域特征:URL 域的域特征。其中包括 whois 信息和 shodan 信息。
- 页面内容特征:从 URL 的页面提取的特征(如果有的话)
下表显示了提取的所有特征的摘要:
URL 的特征向量摘要
运行上面的脚本会产生以下带有 23 个特征的数据。我们将把整数、布尔值和对象列名放在不同的列表中，以便于数据访问。
objects **=** [i **for** i **in** data**.**columns **if** 'object' **in** str(data**.**dtypes[i])]
booleans **=** [i **for** i **in** data**.**columns **if** 'bool' **in** str(data**.**dtypes[i])]
ints **=** [i **for** i **in** data**.**columns **if** 'int' **in** str(data**.**dtypes[i]) **or** 'float' **in** str(data**.**dtypes[i])]
移除高度相关的要素
最线性分析假设预测变量之间不存在多重共线性，即预测要素对不得相关。这种假设背后的直觉是，没有额外的信息添加到具有多个相关特征的模型中，因为相同的信息被其中一个特征捕获。
多重相关特征也表示数据中的冗余特征，删除它们是数据降维的良好开端。通过移除相关特征(仅保留一组观察到的相关特征)，我们可以解决预测值之间的特征冗余和共线性问题。
让我们创建一个简单的相关性网格来观察恶意和良性 URL 的派生特征之间的相关性，并删除一个或多个高度相关的特征。
corr **=** data[ints**+**booleans]**.**corr()
*# Generate a mask for the upper triangle* mask **=** np**.**triu(np**.**ones_like(corr, dtype**=**np**.**bool))*# Set up the matplotlib figure* f, ax **=** plt**.**subplots(figsize**=**(20, 15))*# Generate a custom diverging colormap* cmap **=** sns**.**diverging_palette(220, 10, as_cmap**=**True)*# Draw the heatmap with the mask and correct aspect ratio* sns**.**heatmap(corr, mask**=**mask, cmap**=**cmap, vmax**=**.3, center**=**0,
square**=**True, linewidths**=**.5, cbar_kws**=**{"shrink": .5}, annot**=**True)
互相关系数
然而，我们并不希望删除所有相关的变量——只删除那些相关性非常强、不会给模型增加额外信息的变量。为此，我们为观察到的正相关和负相关定义了某个“阈值”(0.7)。
我们看到大部分高度相关的特征都是负相关的。例如，在 URL 中的字符数和 URL 的熵之间存在 0.56 的负相关系数，这表明较短的 URL 具有
这里我们将创建一个函数来识别和删除多个相关特征中的一个。
预测 URL 的恶意性(决策树)
建模根据以前观察到的数据模式构建了解释数据的蓝图。建模通常是预测性的，因为它试图使用这种开发的“蓝图”来预测基于过去观察的未来或新观察的值。
基于提取的特征，我们想要最好的预测模型来告诉我们一个看不见的 URL 是恶意的还是良性的。因此，我们寻求有用功能的独特组合，以准确区分恶意和良性 URL。我们将经历两个阶段，特征选择，其中我们只选择对预测目标变量有用的特征，并用决策树建模以开发恶意和良性 URL 的预测模型。
特征选择
什么样的变量对识别一个网址是“恶意的”还是“良性的”最有用？在计算上，我们可以通过测试哪些变量“改善”或“未能改善”预测模型的整体性能，来自动选择哪些变量最有用。这个过程被称为“特征选择”。特征选择还有助于降低数据的维数，解决计算复杂性和模型性能的问题。特征选择的目标是获得原始数据的有用子集，该子集以不丢失有用信息的方式预测目标特征(一起考虑所有预测)。尽管特征选择超出了简单的相关性消除，但是对于本文，我们将我们的特征选择方法限制为仅仅保留这些特征。让我们创建一个只包含不相关要素的原始数据子集。
predictor_columns **=** data2**.**columns
d **=** data2[predictor_columns]
x, y **=** d[predictor_columns], data['vt_class']
我们只保留对模型有独特贡献的特征。我们现在可以开始用原始样本的 70%和这 14 个特征来开发模型。我们将保留 30%的样本来评估模型在新数据上的性能。
- 数字服务
- 熵
- 数字
- 数字图表
- 体长
- 数字标题
- 数字图像
- 数字链接
- 每日业务报告(Daily Service Report)
- data-storageequipment 数据储藏设备
- sscr
- 小触须（同 small bristles）
- 机器人
- 哈希 Http
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test **=** train_test_split(x, y, test_size **=** 0.3, random_state **=** 100)
决策树
from sklearn import tree
from sklearn.metrics import accuracy_scorefrom sklearn.externals.six import StringIO
from IPython.display import Image
from sklearn.tree import export_graphviz
import pydotplus
有多种机器学习算法(分类)算法可用于识别恶意 URL。在将 URL 转换成代表性的特征向量之后，我们将“恶意 URL 识别问题”建模为二元分类问题。二元分类模型为只有两种结果“恶意”和“良性”的类训练预测模型。批量学习算法是在以下假设下工作的机器学习算法:
- the entire training data is available before model development and
- the target variable is known before the model training task.
批处理算法是理想和有效的，因为它们是可解释的区分学习模型，使用训练数据点之间的简单损失最小化。决策树是机器学习中的一种批量学习算法。
在决策分析中，决策树是模型得出某些结论的决策过程的可视化表示。决策树背后的基本思想是试图理解什么因素影响类成员，或者为什么一个数据点属于一个类标签。决策树清楚地显示了构成类成员的条件。因此，它们是决策过程的可视化表示。
决策树通过将数据集分解成越来越小的部分来构建预测模型。分割子集的决定基于最大化分割的信息增益或最小化分割的信息损失。从根节点(没有不确定性的最纯粹的特征)开始，通过基于子集的纯粹性创建各种叶节点来形成树。
在这种情况下，决策树将解释每个特征的类别边界，以将 URL 分类为恶意或良性。构建决策树时，有两个主要因素需要考虑:
`- a) What criteria to use in splitting or creating leaf nodes and
- b) tree pruning to control how long a tree is allowed to grow to control the risk of over-fitting.`
决策树算法的 criteria 参数指定控制什么标准(Gini 或熵),而 max_depth 参数控制允许树增长多远。基尼系数是指如果我们根据分支中的分布随机选择一个标签，随机样本被错误分类的概率。熵是对信息(或者说缺乏信息)的一种度量。
不幸的是，由于事先不知道标准和树深度的正确组合，我们将不得不反复测试这两个参数的最佳值。我们针对两个标准测试了 50 次迭代的 max_depth，并可视化了模型准确度分数。
参数调谐
似乎最好的模型是最简单的模型，基尼指数和最大深度为 4，样本外准确率为 84%。此外，最大化熵似乎不会产生好的结果，这表明添加到模型中的新参数不一定给出新的信息，但是可以产生改进的节点概率纯度。因此，我们可以用 max_depth = 4 和 Gini 标准来拟合和可视化该树，以识别哪些特征在分离恶意和良性 URL 中最重要。
建立模型…
*###create decision tree classifier object* DT **=** tree**.**DecisionTreeClassifier(criterion**=**"gini", max_depth**=**4)*##fit decision tree model with training data* DT**.**fit(X_train, y_train)*##test data prediction* DT_expost_preds **=** DT**.**predict(X_test)
想象这棵树…
dot_data **=** StringIO()
export_graphviz(DT, out_file**=**dot_data,
filled**=**True, rounded**=**True,
special_characters**=**True,feature_names**=**X_train**.**columns, class_names**=**DT**.**classes_)
graph **=** pydotplus**.**graph_from_dot_data(dot_data**.**getvalue())
Image(graph**.**create_png())
树
预测模型的准确性对 max_depth(树修剪)和分裂质量标准(节点分裂)的参数调整非常敏感。这也有助于实现最简单的简约模型，防止过度拟合，并在看不见的数据上表现得一样好。这些参数特定于不同的数据问题，最好测试不同参数值的组合。
该模型显示，恶意 URL 的脚本与特殊字符比率(sscr)较低，URL 字符相对更“有序”或更单调。此外，恶意网址的域名可能在 5-9 个月前已经过期。我们还知道“恶意广告”的问题，即骗子获得过期合法域名的所有权，以分发可下载的恶意代码。最后，大概良性 URL 最显著的特点就是长寿。他们似乎缓和了 HTML 正文内容中脚本与特殊字符的比例，使域名寿命延长到 4-8 年。
用深度学习预测一首歌的音乐情绪。
使用 Python 上的 Keras 和 Tensorflow 库，通过神经网络模型预测音乐曲目情绪的酷方法。
乌列尔·索伯兰斯在 Unsplash 上拍摄的照片
音乐是一种表达我们情感的强大语言，在很多情况下，它被用来治疗我们生活中的困难时刻。情绪和心情可以很容易地反映在音乐中，当我们做运动时，我们倾向于听充满活力的音乐，同样，当我们焦虑或疲劳时，一首优美放松的歌曲可以帮助我们平静下来。这就是为什么我试图弄清楚分类模型如何帮助确定特定曲目的情绪。
在这篇文章中，我将展示并解释我如何使用一个用于分类的多类神经网络和 Spotify 提供的一个很酷的数据集来实现这个想法。所以不要再说了，让我们开始工作吧！！。
*我的 Github 存储库中的全部代码、脚本、笔记本和数据(点击此处)
所需工具:
- 熊猫和 Numpy 进行数据分析。
- Keras 和 Tensorflow 构建深度学习模型。
- Sklearn 验证模型。
- Seaborn 和 Matplotlib 绘制一个漂亮的图形。
- Spotipy Python 库(点击这里了解更多信息)。
- Spotify 凭证访问其 API 和数据采集(点击此处了解更多信息)。
Spotify 音频功能:
Spotify 使用一系列不同的功能对歌曲进行分类。我从 Spotify 网页复制/粘贴信息。
- 声音:一种置信度，从 0.0 到 1.0，表示音轨是否是声音的。1.0 表示音轨是声学的高置信度。
- 可跳舞性:可跳舞性描述了一个曲目在音乐元素组合的基础上适合跳舞的程度，包括速度、节奏稳定性、节拍强度和整体规律性。值 0.0 最不适合跳舞，1.0 最适合跳舞。
- 能量:能量是一个从 0.0 到 1.0 的度量，代表强度和活动的感知度量。通常，高能轨道感觉起来很快，很响，很嘈杂。例如，死亡金属具有高能量，而巴赫前奏曲在音阶上得分较低。对该属性有贡献的感知特征包括动态范围、感知响度、音色、开始速率和一般熵。
- 乐器性:预测音轨是否不包含人声。“Ooh”和“aah”在这种情况下被视为乐器。Rap 或口语词轨道明显是“有声的”。乐器度值越接近 1.0，轨道不包含人声内容的可能性就越大。高于 0.5 的值旨在表示乐器轨道，但随着该值接近 1.0，置信度会更高。
- 活跃度:检测录音中是否有观众。较高的活跃度值表示音轨被现场执行的概率增加。高于 0.8 的值很有可能表示该音轨是实时的。
- 响度:轨道的整体响度，以分贝(dB)为单位。响度值是整个轨道的平均值，可用于比较轨道的相对响度。响度是声音的质量，是与体力(振幅)相关的主要心理因素。值通常在-60 和 0 db 之间。
- 语速:语速检测音轨中是否存在口语单词。越是类似语音的录音(例如脱口秀、有声读物、诗歌)，属性值就越接近 1.0。高于 0.66 的值描述可能完全由口语单词组成的轨道。介于 0.33 和 0.66 之间的值描述可能包含音乐和语音的轨道，可以是分段的，也可以是分层的，包括说唱音乐。低于 0.33 的值很可能代表音乐和其他非语音类轨道。
- 配价:从 0.0 到 1.0 的一种量度，描述音轨所传达的音乐积极性。高价曲目听起来更积极(例如，快乐、愉快、欣快)，而低价曲目听起来更消极(例如，悲伤、沮丧、愤怒)。
- 速度:轨道的整体估计速度，单位为每分钟节拍数(BPM)。在音乐术语中，速度是给定作品的速度或步调，直接来源于平均节拍持续时间。
1.解释和分析数据:
为了获得数据，我必须使用 Spotipy 库创建一系列函数。这个库有助于 Spotify 服务自动下载更多关于播放列表、歌曲、艺术家音乐等的技术信息(如上所述)。出于本文的主要目的，我不会提及我是如何获得数据的，但我会解释数据由什么组成。
您可能知道，分类问题使用带标签的数据，所以我必须创建这些标签。我决定创建 4 个类别来标记曲目，这些类别是“充满活力的”、“冷静的”、“快乐的”和“悲伤的”。我根据下面这篇文章选择这些类别，这篇文章解释了按情绪对音乐进行分类的最佳方式。
本文将涵盖使用各种 DSP 和音乐理论技术分析音乐，包括节奏，和声…
sites.tufts.edu](https://sites.tufts.edu/eeseniordesignhandbook/2015/music-mood-classification/)
然后，我在 Spotify 上搜索了一些基于这 4 个标签的不同音乐曲目的播放列表(每个标签 200 首曲目)，最后，我将所有这些曲目连接到由每个情绪标记的主数据框中。你可以在我的 GitHub 知识库上查看这个数据集(点击这里)
主要数据有 800 行和 18 列，但是为了信息还原的目的，我决定使用长度、舞蹈性、声音性、能量、I 乐器性、活跃度、化合价、响度
我通过标签对数据框进行分组，计算轨迹特征的平均值。我得到了以下结果:
使用平均统计分组的数据帧。(图片由作者提供)
做这个简单的分析，我很快注意到最流行的歌曲是快乐的，悲伤的歌曲往往有很长的长度，充满活力的歌曲节奏最快，平静的歌曲往往是有声的。
2.构建模型:
2.1-数据预处理:
为了标准化特征，我使用了最小最大缩放器来缩放范围为[0，1]的值，并保持原始分布的形状。我还对 4 个标签进行了编码，因为神经网络使用数值来训练和测试。最后，我将数据分成 80%用于训练和 20%用于测试。
*#Libraries to pre-process the variables*
from sklearn.preprocessing import LabelEncoder,MinMaxScaler
from sklearn.model_selection import train_test_split
我用来处理数据的代码:
*#Define the features and the target*
col_features = df.columns[6:-3]
X = df[col_features]
Y = df['mood']*#Normalize the features*
X= MinMaxScaler().fit_transform(X)*#Encode the labels (targets)*
encoder = LabelEncoder()
encoder.fit(Y)
encoded_y = encoder.transform(Y)*#Split train and test data with a test size of 20%*
X_train,X_test,Y_train,Y_test = train_test_split(X,encoded_y,test_size=0.2,random_state=15)
标签编码如下:
标签及其编码。(图片由作者提供)
2.2 创建模型:
为了建立这个模型，我使用了 Keras 库，这个库被设计成能够快速实验深度神经网络**，专注于用户友好。我的主要目标是将轨迹分为 4 类情绪(平静、精力充沛、快乐和悲伤)，因此我的模型由一个多类神经网络组成，输入为 10 个特征， 1 层有 8 个节点，4 个输出有输出层。我还需要使用一个分类器作为估计器，在这种情况下，分类器是 KerasClassifier ，它将我之前使用定义的神经网络模型创建的函数作为参数。激活函数对应于一个整流线性单元 ( Relu) ，损失函数是一个 L 逻辑函数和亚当梯度下降算法是优化器。**
*#Libraries to create the Multi-class Neural Network*
from keras.models import Sequential
from keras.layers import Dense
from keras.wrappers.scikit_learn import KerasClassifier
from keras.utils import np_utils*#Import tensorflow and disable the v2 behavior and eager mode*
import tensorflow as tf
tf.compat.v1.disable_eager_execution()
tf.compat.v1.disable_v2_behavior()
- 重要提示:我禁用了 TensorFlow 的急切执行和 v2 行为，因为我一直试图理解和学习库在那些模式下是如何工作的(抱歉我是 Tensorflow 的新手呵呵)。
我用来建立神经网络的代码:
*#Function that creates the structure of the Neural Network*
def base_model():
*#Create the model*
model = Sequential()*#Add 1 layer with 8 nodes,input of 4 dim with relu function*
model.add(Dense(8,input_dim=10,activation='relu'))*#Add 1 layer with output 3 and softmax function*
model.add(Dense(4,activation='softmax'))*#Compile the model using logistic loss function and adam * optimizer, accuracy correspond to the metric displayed
model.compile(loss='categorical_crossentropy',optimizer='adam',
metrics=['accuracy'])
return model*#Configure the estimator with 300 epochs and 200 batchs. the build_fn takes the function defined above*.
estimator = KerasClassifier(build_fn=base_model,epochs=300,
batch_size=200)
3.3 评估模型:
使用 K 倍交叉验证我评估了使用训练数据的估计器。拆分的次数是 K= 10 洗牌所有的值。
*#Library to evaluate the model*
from sklearn.model_selection import cross_val_score, KFold
我用来评估模型的代码:
*#Evaluate the model using KFold cross validation*kfold = KFold(n_splits=10,shuffle=True)results = cross_val_score(estimator,X,encoded_y,cv=kfold)print("%.2f%% (%.2f%%)" % (results.mean()*100,results.std()*100))
模型的精确度是每个折叠精确度的平均值，在本例中，精确度为 72.75 %。
3.4 训练模型:
是时候训练模型了！所以让我们开始编码:
*#Train the model with the train data*
estimator.fit(X_train,Y_train)*#Predict the model with the test data*
y_preds = estimator.predict(X_test)
需要特别提到的是，模型是用 640 个样本(主数据的 80%)训练的。
模型训练时最后一个时期的一些输出。
历元、时间、损失和训练过程中模型的准确性。(图片由作者提供)
3.多类神经网络的精度；
最后，为了评估模型的准确性，我使用 Seaborn 库和 Matplotlib 绘制了一个混淆矩阵。我还计算了 Sklearn 库提供的准确率分数。
*#Create the confusion matrix using test data and predictions*
cm = confusion_matrix(Y_test,y_preds)*#plot the confusion matrix*
ax = plt.subplot()
sns.heatmap(cm,annot=True,ax=ax)
labels = target['mood'].tolist()
ax.set_xlabel('Predicted labels')
ax.set_ylabel('True labels')
ax.set_title('Confusion Matrix')
ax.xaxis.set_ticklabels(labels)
ax.yaxis.set_ticklabels(labels)
plt.show()*#Show the accuracy score*
print("Accuracy Score",accuracy_score(Y_test,y_preds))
音乐情绪的混乱矩阵
最终的准确度分数为 76% ，看一下混淆矩阵，我注意到我的模型可以很好地分类平静和悲伤的歌曲，但在处理充满活力和快乐的歌曲时有一些问题。我可以修改一些参数，如批次大小、时期，或者聚合或删除一些追踪要素来训练我的模型，从而帮助提高模型的准确性。
4.享受按情绪分类音乐的乐趣:
我想展示如何预测一首歌的情绪，这首歌可能你懒得完全听，但你想知道这首歌会让你跳舞还是哭。
我的 Github 存储库(点击此处)上有一个名为 helpers.py、的脚本，使用它你只需要在 Spotify 上为开发者创建一个应用程序(点击此处了解更多信息)并获得一个 Client_id 、 Client_secret 和重定向 URL 。使用这个脚本，你可以通过 这个多类分类模型的一点帮助，下载预测任何歌曲 情绪所需的特性。(没错，就像披头士的歌)。
首先，我们需要获得 Spotify 应用程序提供的任何歌曲的 Spotify URI。比如，我会预测周末的炫目灯光的心情。
眩目灯光下的 Spotify URI:Spotify:track:0 vjijw 4g luzamy D2 vxmi 3 b
Spotify 应用面板(作者截图)。
然后，我会将歌曲的 Spotify URI 传递给一个我定义的名为 predict_mood 的函数。该函数将歌曲的 Id 作为参数，并包含在创建的神经网络模型中
预测情绪函数的代码是:
*#Import the Script helpers.py*
from helpers import *def predict_mood(id_song):
*#Join the model and the MinMaxScaler in a Pipeline*
pip = Pipeline([('minmaxscaler',MinMaxScaler()),('keras',
KerasClassifier(build_fn=base_model,epochs=300,
batch_size=200,verbose=0))])*#Fit the Pipeline*
pip.fit(X2,encoded_y)
*#Obtain the features of the song (Function created on helpers.py)*
preds = get_songs_features(id_song)*#Pre-processing the input features for the Model*
preds_features = np.array(preds[0][6:-2]).reshape(-1,1).T
*#Predict the features of the song*
results = pip.predict(preds_features)
mood = np.array(target['mood'][target['encode']==int(results)])
#Obtain the name of the song and the artist
name_song = preds[0][0]
artist = preds[0][2]*#Store the name,artist and mood of the song to print.* result_pred=print("{0} by {1} is a {2} song".format(name_song,
artist,mood))return result_pred
《致盲之光》的 Spotify URI 是“spotify:track:0 vjijw 4g luzamy D2 vxmi 3 b”但是predict_mood 函数取 Id，所以我只需要取“Spotify:track”之后的代码:在这种情况下，Id 是 0VjIjW4GlUZAMYd2vXMi3b。
结果显示，周末的强光使人精神振奋。很准确，你不这样认为吗？。
5.结论
深度学习算法在实现与你喜欢的事情相关的想法或项目方面有很多乐趣。就我而言，我非常喜欢音乐，所以我可以用这些知识创造一些很酷的方法来帮助我自动完成一项可能需要很长时间才能完成的任务。我还可以更多地了解这个神奇的数据科学世界，以及我对音乐品味的倾向。
我的文章:
- 一个友好的公共汽车和分析人员到达网页上的部门
- 使用 Python 和 R 对音乐进行聚类，在 Spotify 上创建您的播放列表
- 一种在网络上进行公共汽车和分析工作的交互形式。
- 从你电脑的音乐文件夹中自动播放 Spotify 上的播放列表
参考资料:
- https://towards data science . com/cross-validation-explained-evaluating-estimator-performance-e51e 5430 ff 85
- https://towards data science . com/clustering-music-to-create-your-personal-playlists-on-Spotify-using-python-and-k-means-a39c 4158589 a
- https://machine learning mastery . com/multi-class-classification-tutorial-keras-deep-learning-library/
使用人工神经网络预测下一个音高
太空人真的需要使用照相机吗…
来源:来自 Pixabay 的 Vic Cherubuni
欺骗和猜测:在棒球中获得优势
休斯顿太空人队在 2017 年赢得了世界大赛。他们有一个令人难以置信的运行，击败了波士顿红袜队，纽约洋基队，和(谢天谢地)洛杉矶道奇队。然而，美国职业棒球大联盟最近的调查发现，太空人作弊。尽管休斯顿(和旧金山)的球迷不关心太空人队如何击败道奇队赢得世界大赛，但太空人队被发现利用技术窃取标志。他们从一个非法的中外场摄像机转播了对方接球手的手势。在俱乐部会所隧道里的教练会收到摄像头的反馈，并通过敲打垃圾桶来发出超速信号，而不是敲打快速球。这里有一个视频分解，原谅NSFW语言。
在棒球比赛中，知道什么球会到来是一个不可思议的优势。捕鼠人试图不被发现是有原因的。作为一名前球员，我可以证明知道接下来会发生什么是非常有用的。此外，有一个复杂的方法和预测投球通常是区分伟大击球手的特征。这里有一个技安进攻方式的例子:
“我通常打区域联防，大部分时间我会坐快速球。但是如果有一个人会给我一个 2-1 的变速球，我会坐以待毙。或者(与)一个投滑球的左撇子，我会坐滑球。我会坐在这些球场上，然后让一切都过去。”
此外，名人堂外野手和前全垒打王汉克·亚伦曾经说过:
“猜测投手会投出什么是成功击球手的 80%。另外 20%只是执行力。”
在一个复杂的击球方法中考虑了多个变量后，我意识到猜测或坐下来，下一次投球是机器学习的完美应用，特别是香草前馈神经网络。
先决条件:
我是一名数据科学家新手，但这篇博文假设我熟悉基本的深度学习和 python 概念。
数据收集和清理:
Statcast 是棒球领域最近的一项创新，它使整个项目成为可能。Statcast 给出了每一个音高的 CSV 文件，每个音高都有大量的数据，如速度、音高位置、音高结果、位置等。我想把这种算法用在拥有高级混合音高的投手身上，也就是说，以相对高的频率投出至少三种主要音高。米科拉斯在 2019 年投掷了大约 2800 个跟踪球，斯特罗曼投掷了大约 3000 个球。
在推特上关注达伦·威尔曼，替我谢谢他。他在 Statcast 建立了基础设施来实现这一点。推特:@darenw
斯特罗曼的球场阿森纳/来源:Statcast
米科拉斯的球场/来源:Statcast
我在 excel 中做了最初的数据清理，因为无可否认，作为一个以金融为导向的人，我对 excel 比对熊猫更有信心。清理数据后，我根据自己的经验和职业球员的陈述选择了以下相关变量来预测投球:以前的投球类型(string)、以前的投球结果(string)、击球手的惯用手(string)、投球数(ints)、跑垒员(boolean)、出局数(int)、局数(int)和当前跑垒差(int)。
初始模型:
以下是我用过的库。如果有疑问，就扔进去:
输入库后，我读入 CSV 文件，然后去熊猫工作。我们的第一步是一次性编码所有的字符串值。一种热门的编码方式是将分类数据转换为数字数据，以便机器学习模型可以用它来执行逻辑回归。
所有数据都整齐地组织在 NumPy 数组中，我使用 Sklearn 库创建了一个随机训练/测试分割。在试验了不同的优化器之后，我使用 Keras API 构建了一个带有一个隐藏层的前馈模型。
对于斯特罗曼来说，这在预测超过 600 个音高时产生了 46%的准确率。这个结果比随机预测的概率要高得多。
然而，46%让我深感不满。我很反感。我想做得比抛硬币更好，即使硬币真的是六面的，形状也很奇怪。我想到了一个更准确的方法来建立这个模型。
解决问题:
当检查 Astros 中继系统时，我意识到他们没有中继准确的音高。他们转播快速球或慢速球。猜测准确的下一个音高不像猜测音高的特征那样实用，比如移动和速度。我回到 CSV，将相似的球种归类为硬球(两缝线快速球(伸卡球)和四缝线快速球)、破球(切球和滑球)和软球(曲球和变速球)。这更类似于宇航员使用的撞击系统。
通过稍微改变模型和输入，我能够获得接近 60%的准确率。
应用和结论:
在过去的 20 年里，棒球一直是体育界应用定量方法的领导者。从应用上来说，我觉得这个模型可以作为打者或者任何东西的学习工具来真正打败道奇队。击球手可以使用这个模型来策划球探报告，并在第二天面对先发投手的击球上进行模拟。此外，我知道这个模型还可以改进。我使用的变量并不是猜测音高艺术的全部。我很想从一些玩家和捕手(他们称之为游戏)的头脑中了解如何改进这个模型。请在下面找到我的 git 库，里面有模型和数据。
[## danielajk99/pitch_predictor
用于基音预测的人工神经网络。为 danielajk99/pitch_predictor 开发做出贡献，创建一个…
github.com](https://github.com/danielajk99/pitch_predictor)
特别感谢布兰登·威廉姆斯和诺亚·鲁本帮助我完成了这个项目。还要特别感谢 Ulk Aslak 教授和 Lucian Leahu 教授，他们向我介绍了人工神经网络。
数据和球场画面来自棒球专家的 Statcast 搜索工具。按照这个 链接 在 Statcast 中创建自己的数据集。
丹尼尔·丹·金是乔治敦大学的本科生。这个博客是他在丹麦留学项目的最终项目的一部分(由于新冠肺炎的原因，现在实际上是一个在家学习的项目)。在推特上关注我:@danielajk99。
利用经济指标预测名义 GDP:一种数据科学方法
利用流行的分类模型预测国内生产总值
凯文·Ku 在 Unsplash 上的照片
这篇文章讨论了使用网络上的数据预测名义 GDP(国内生产总值)。问题陈述是从真实来源收集数据，执行探索性数据分析(EDA)，训练模型，以及预测名义 GDP(加拿大)。作为一名 IT 人员，我对这些经济学相关术语知之甚少。这就是为什么首先要做的是熟悉这些概念。接下来是用 R 编程语言学习和执行这个项目，并熟悉流行的 R 包，如 tidyverse、ggplot、caret 等。在实践中，很好地利用在线资源可以很容易地习惯 R 及其语法。
步入任务，遇到各种经济统计又名指标，分类为滞后和领先指标。经济指标可能对一个国家的名义 GDP 有直接影响。我们的方法是用这些指标来预测加拿大的名义 GDP。在研究了各种在线博客、文章和论文中的这些指标后，最终确定了以下指标:
人口
难民人口
实际利率
国内公司数量
旅游服务(占进口服务-国际收支的百分比)
税收
房地产市场
劳动生产率
政府债券收益率 10 年期利率
个人汇款(收到&以美元支付)
铁路运送旅客
航空运送旅客
通货膨胀
收入增长
失业率
政府赤字
消费价格指数
商品价格指数
欧元对加元的兑换率
美元对加元的兑换率
多伦多证券交易所交易价值
多伦多证券交易所交易量
入境国际游客
数据来自加拿大统计局、T2、世界银行、T4 统计局和其他来源。这些指标被假定为特征，这些特征的数据是在 2009 年至 2018 年间收集的。数据集分为两部分，训练集和测试集。
T
探索性数据分析
在数据集组装、清理和整齐排列后，开始对数据集进行探索性数据分析。使用标准正态分布的概率密度 dnorm 函数进行数据标准化。为了检查数据是否正常，对每个特征应用夏皮罗-维尔克正态性检验**。除 3 个变量外，所有变量的 p 值都大于 0.05，这表明它们是正常的。为了分析数据集特征之间的相关性，使用了 caret 包的find correlation函数，去除了高度相关的属性。由于数据点较少，即 8 个(2009 年至 2016 年)，必须降低数据的维度。这是通过使用 prcomp 函数应用主成分分析(PCA)** 完成的。至此，EDA 结束了。现在让我们开始建模吧！参考 GitHub 链接… 点击这里！
线性回归建模
在这种类型的建模中，PCA 值用于使用 lm() 函数训练简单的线性回归模型。使用汇总功能对生成的模型进行汇总。关键观察是 p 值小于显著性水平(< 0.05)。这意味着我可以安全地拒绝零假设，因为预测系数的系数β为零。此外，模型的 F 统计值为 428.6(越多越好)。因此，可以得出结论，该模型具有统计学意义。为了验证该模型，我对 2017 年和 2018 年的值进行了测试。模型的总体准确率为 98.31%。文章末尾的表格提供了关于实验的更多细节。
随机森林造型
对于这个建模，R 语言的 randomForest 库使得选择特征和创建随机森林模型变得更加容易。 importance() 函数根据数据叶节点杂质找到并显示每个特征的重要性。这些数据存储在一个新的数据集中，该数据集根据重要性值进行排序。这里的技巧不是选择重要性值仅较高或较低的特征，而是选择中和频谱的值，即较高、较低和中等重要性值的组合，以保持平衡。使用试错法挑选了六个特征。
个人汇款(已收)
实际利率
旅游服务进口
国债收益率 10 年期利率
政府赤字
国内公司数量
使用随机森林函数对模型进行训练。 mtry 参数设置为 6(特征数量)，树的数量( ntree )设置为 1000。该模型对 2017 年和 2019 年 NGDP 的预测平均准确率为 95.68%。
支持向量机建模
这是我们实验清单上的最后一个。为了在我们的使用中利用 SVM 建模，我们使用了提供 svn 函数的 R 的 e1071 包来建模。因为我们的用例不是一个分类用例，所以必须应用回归选项。为此， svm() 函数提供了一个类型参数，使我们能够进行选择。相应地，选择EPS-回归类型，并将内核参数设置为径向。关于内核的更多细节，点击这里。这里预测的平均准确率为 98.18%。下面的图表提供了实验的更多细节。
显示实际值与预测值的图表
总之，尽管数据集很小，但所有模型的预测准确率都在 95%以上。当使用 PCA 值训练线性模型时，得到的输出比其他的更好。然而，随机森林可以被认为是最不可靠的，因为它使用数据集的最少特征。如果不是准确性或可靠性，该结果至少倡导了经济指标在预测名义 GDP 中的意义。
PS:这个实验是提交给加拿大统计局商业数据科学家挑战 2019/2020 的一部分。
使用优先机器学习预测奥斯卡
奥斯卡和他们的优先投票让我创造了一种新颖的机器学习方法来模仿这种投票系统
去年对电影来说是很棒的一年，如果你像我一样，沐浴在电影票热潮的余晖中，仍然在电影院看很多电影，你知道好莱坞往事、寄生虫、 1917 和更多的电影提供了独特的电影体验。每年的奥斯卡周日，好莱坞都会聚在一起，给自己一个大大的鼓励。当晚最大的奖项是最佳影片奖，它可以将一部电影载入电影史。与奥斯卡之夜颁发的其他 23 个奖项不同，梦寐以求的最佳影片奖是通过一种称为优先投票的方法选出的，这种方法比传统投票更复杂。为今年的奥斯卡做准备，并了解优先投票让我写了一些程序，用机器学习来模拟这个投票系统。
Reddit 用户 u/Tillmann_S 的 2020 奥斯卡艺术
在这篇文章中，我:
- 选择用于预测奥斯卡的数据
- 从数据科学的角度探索优先投票的工作原理
- 演示一个我自己设计的方法，我称之为优先投票随机森林
- 模拟最佳影片投票的幕后发生的事情
- 预测今年的最佳影片得主
我没有在本文中包含我的任何代码，但是这里是存储库，包含我在这个分析中使用的笔记本
如何预测奥斯卡:数据集
为了使用机器学习来预测任何事情，我们需要一个有意义的数据集来训练我们的模型。在最佳影片竞赛的情况下，我们有 9 部 2019 年提名的电影。尽管我对奥斯卡充满敬意(毕竟我有足够的兴趣写这篇文章)，但我毫无保留地认为，今年的最佳电影将赢得奥斯卡最佳影片奖。该学院由电影行业各个领域的数千名成员组成，他们每个人都有自己的偏见，这导致了他们的投票。因为投票背后有真实的人，我们不能依赖票房利润或总评论家分数等电影质量的数字指标。但是你知道什么和电影制作人的投票密切相关吗？其他电影人投票。
还有许多其他的颁奖典礼组成了“颁奖季”，像美国演员工会奖和导演工会奖这样的活动的投票者通常也是奥斯卡奖投票者。利用早期奖项如 SAGs、DGAs、PGAs、金球奖和 BAFTAs 的结果，并结合奥斯卡信息如提名数，我可以训练一个模型来预测今年的最佳影片获奖者。为了获得一致的电影数据和命名惯例，我从维基百科上收集了每个颁奖典礼的提名者和获奖者的数据，并使用熊猫和美丽的汤包将它们合并到 Python 中的一个数据集。
优先投票的工作原理
优先投票，也称为即时决胜投票，通常用于有许多候选人争夺一个席位的情况。自 2009 年以来，奥斯卡一直使用这种计票系统来决定最佳影片竞赛，当时该领域从五个提名人扩大到十个。在优先投票中，选民不是为一部电影投票，而是提交一份对所有选项进行排名的投票，排名第一的选项将作为该电影的投票。然后开始一个迭代过程，其中最不受欢迎的电影被淘汰，所有的投票被重新排序，直到一部电影获得超过 50%的#1 选票。在一部电影从所有选票中被淘汰后，之前将被淘汰的电影放在第一位的选票现在将第二位移到第一位，这增加了剩余电影的票数。这一过程一直持续到一部电影获得超过 50%的#1 选票，然后它被宣布为获胜者。这种消除过程的模拟如下所示。
图 1:优先投票淘汰的模拟。由我的偏好随机森林生成
优先投票法的批评者声称，这种方法奖励容易喜欢或无争议的电影，因为无争议的电影将在人们的排名中处于中间位置，而有争议的电影可能在一些人的投票中名列前茅，但在其他人的投票中垫底，因此它们很容易被淘汰。去年，当更具艺术性的电影《T2》《罗马》《T3》输给更具大众吸引力的电影《T4》《绿皮书》时，这种效应就显现出来了。
图片来源:左——阿方索·卡隆/网飞，右——环球影业/参与者/梦工厂
优先抽签随机森林
我们过去已经看到，优先投票可以改变最佳影片竞赛的结果，因此我创建了一个模型来反映这种独特的计票方法。随机森林分类器模型通过使用多个去相关的决策树分类器进行预测。这里有一篇文章更侧重于传统随机森林如何工作的细节。一般来说，一个随机森林会根据树叶的大小将每棵树的“投票”作为一个分数，并根据哪个类别在所有树中拥有最多“投票”来选择最终标签。对于这个优先投票随机森林，我们改为使用测试集中每部电影的 ProbA 值，并使用它们来创建电影的第 1-9 名排名。ProbA 值是该项在“赢家”类别中的可能性，代表比二元“赢家”或“输家”分类更软的预测。这种更软的预测允许我们将预测从布尔分类改为范围。每个决策树产生一张选票，一旦整个森林产生了他们的选票，优先选票淘汰的迭代过程开始确定森林对获胜者的选择。通过使用排名而不是挑选一个类别，我的优先投票随机森林保存了传统随机森林会丢弃的信息，并在优先投票的消除和重新排名阶段再次使用。
图 2:单个决策树对测试集的投票
模拟奥斯卡颁奖典礼
利用我的优先投票随机森林，我模拟了今年的最佳图片竞赛。为了消除每棵决策树的相关性，我改变了每棵树看到的奖励，类似于随机森林的 max_features 超参数。在这个模拟中， max_features 代表投票的学院成员可能属于哪个公会，或者他们对该季其他奖项的关注程度。我还为每个决策树添加了一个随机噪声特征来进行训练，代表每个投票者对某些电影的天生偏见。该学院由大约 7000 名不同的投票者组成，所以我启动了我的森林，很快就产生了 7000 张选票。经过 6 轮淘汰最后一部电影后，第一部电影获得了超过 50%的第一票，我的模型选出了最佳影片奖…
图 3:6 轮优先投票淘汰后的最终排名。当电影《1917》获得超过 50%的票数时，这一过程就停止了。
最终预测
我的优先投票随机森林是一种模拟奥斯卡的新颖方法，我希望它能帮助你了解最佳影片投票和随机森林分类器的一些内容，但除了优先投票，让我们进入正题，真正预测这些坏男孩。利用我搜集的获奖电影数据集，我实现了 H2O 强大的 AutoML 工具来训练 100 个不同的随机森林、XBGT 和深度学习模型，用各种参数来预测今年的奥斯卡奖。AutoML 选择了 XGBoost 模型，该模型在交叉验证中正确预测了 159 部电影中的 147 部电影的奥斯卡奖结果。这个模型的漩涡预测了哪部电影？还有 1917 ！看起来这部电影的前景不错，因为优先投票随机森林和我的 AutoML 模型都预测到了。
图片来源:环球影业，弗朗索瓦·杜哈梅尔
链接和推销:
抓取代码灵感来自 Github 用户 Buzdygan
旧金山大学 MSDS 分校
预测过去:验证的概念
如何确定他/她的模型是否准备好分析看不见的数据？
首先，我们必须明白，虽然传统统计学更侧重于理解给定数据集的复杂性，但数据科学旨在理解所述复杂性，对未来情景做出预测或处方。
因此，评估模型可靠性的主要部分是通过验证。验证的核心思想是:既然我们无法拥有未来事件的数据，我们就必须“预测过去”。
那是什么意思？
因为我们已经知道过去事件(X)的结果(Y)，在验证中，我们评估该模型预测事件的结果有多准确。
因此，我们必须从数据集中删除一条信息，我们将把它伪装成未来的数据，模拟生产环境。这个拆下来的零件不是用于模型训练，而是作为一个测试。
这可以通过几种方式来实现。最简单的是将数据分成两组，一组用于建模，另一组用于测试。另一方面，这种方法的演变是交叉验证。
与其把数据只分成两组，为什么不分成五组呢？四个用于模型训练，一个用于精度测试？如果您错过了窍门，这可能会适得其反:如果您将数据集分成五部分，您将进行五种不同的评估，并更改哪一部分是验证集。
通过这种方式，连续的测试可以让你对你的模型的可靠性有一个更广阔的视野。
不过，要记住的重要一点是，如果你没有使用时间序列，那么一定要随机划分样本。
如果你有一个“时间序列”数据集呢？
我们称时间序列为按时间顺序排列的数据集。这个小定义已经表明，传统的随机分裂方法没有用，因为我们可能会制造混乱，甚至使模型不适合。
因此，有两种主要方法和一种非常模糊的方法来验证时间序列数据集:扩展和滑动是主流，交替阻塞是一种非常非传统但非常强大的方法。
这些方法依赖于根据时间非随机地划分数据集。他们每个人都使用不同的分裂组合来验证。例如，2018 年的数据可用于训练模型，以预测 2019 年第三季度的销售额，等等。
我打算稍后再写关于时间序列验证的文章。这篇短文的主要目的是说:验证旨在通过向模型隐藏数据集的一部分，并使其预测隐藏的部分，来模拟生产和新数据的条件。
然后，通过计算误差量来评估可靠性。通过这种方式，可以检查对模型所做的更改在统计上是否显著，以及预测未来数据是否可靠。
预测二手车的价格
帕克·吉布斯在 Unsplash 上的照片
具有特征选择和特征工程的 XGBoost 实现
本文旨在分析二手车数据集的属性，以便在使用 XGBoost 模型预测这些汽车的价格之前执行特征选择和特征工程。数据处理、清理、特征选择和建模的整个过程在下面的部分中示出。术语“属性”、“变量”和“特征”在本文中可互换使用，具有相同的含义，即它们指的是数据集的列。
我们先来了解一下特征选择和特征工程是什么意思。
特征选择
选择那些与输出变量预测最相关的输入特征的过程称为特征选择。这有助于减少数据的维度，并确保模型不会学习那些可能降低其准确性的不相关特征。
有关功能选择的更多信息，请点击此处。
特征工程
在使用或不使用现有要素的情况下，为给定数据集获取新要素以提高该数据集模型性能的过程称为要素工程。
有关特征工程的更多信息点击此处。
现在让我们深入到实现部分。
履行
导入必要的库
我们导入了预测过程所需的以下库。
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from sklearn.preprocessing import StandardScaler, LabelEncoder
from sklearn.model_selection import train_test_split
import math
import matplotlib
import seaborn as sns%matplotlib inline
加载数据集
这里使用的数据集是加载的二手车数据集。它拥有从易贝收集的 37 万行数据和描述每辆二手车细节的 28 个属性。但是使用了具有 50000 行的总数据集的子集。这些数据的内容是德语，然后翻译成英语。这个原始数据集的翻译子集可以在这里找到。数据集属性的描述如下所示:
df = pd.read_csv(‘../input/data.csv’, sep=’,’, header=0, encoding=’cp1252')
df.describe()
从数据描述中我们可以看到，有些数值变量有不能存在的值。例如， yearOfRegistration 变量的最小值为 1000，最大值为 9999，而这两个值是不可能的。这些值将在接下来的代码块中处理。
删除不必要的列和不允许的值
然后，我们删除数据集中对预测汽车成本没有意义的属性。
df.drop([‘name’,’seller’, ‘offerType’, ‘dateCrawled’, ‘lastSeen’], axis=’columns’, inplace=True)
然后删除属性值不可能(不允许)的行。
df=df[(df.price>=100) & (df.price<=200000) & (df.yearOfRegistration>=1950) & (df.yearOfRegistration<=2019) & (df.monthOfRegistration>0) & (df.monthOfRegistration<=12) & (df.powerPS>10) & (df.powerPS<1000)]df.describe()
从上面的数据描述中我们可以看到，属性现在只有可能的值。
NA 值的插补
我们首先需要找出哪些属性具有 NA 值，即空字段。
df.isna().sum()
所以 NA 值只存在于分类变量中。为了填充每个空字段(插补，使用该属性的模式值。
var=[‘vehicleType’,’gearbox’,’model’,’fuelType’,’notRepairedDamage’]
for i in var:
df[i].fillna(df[i].mode()[0],inplace=True)
df.isna().sum()
使用目标变量的特征可视化
根据价格变量绘制特征注册年份、注册月份、和邮政编码。
for i in [‘yearOfRegistration’, ‘monthOfRegistration’,’postalCode’]:
sns.jointplot(x=i, y=’price’,data=df[[i,’price’]],size=7)
注册年份与价格
注册与价格
邮政编码与价格
为每对变量显示的图以及皮尔逊相关性和 p 值可用作特征选择度量。第二个图显示月注册与*价格的相关值非常低，而 p 值很高。*所以这个变量被删除了。
df.drop([‘monthOfRegistration’],axis=’columns’,inplace=True)
特征工程——创造新的特征“ageOfCar”
这个数据集有一个属性 dateCreated ，这是在易贝上创建汽车广告的日期。通过使用这个属性和 yearOfRegistration 属性，我们可以创建一个告诉我们汽车年龄的新特性( ageOfCar )。
df[‘ageOfCar’]=pd.DatetimeIndex(df[‘dateCreated’]).year-df[‘yearOfRegistration’]
sns.jointplot(x=’ageOfCar’,y=’price’,data=df[[‘ageOfCar’,’price’]],size=7)
车龄与价格
从上面的图中我们可以看到，皮尔逊相关性的大小与注册年份和价格的皮尔逊相关性的大小相同，但是为负。这意味着随着车龄的增加，价格会下降，这在现实生活中也是如此。
现在我们删除变量 yearOfRegistration 和 dateCreated。
df.drop([‘yearOfRegistration’,’dateCreated’],axis=’columns’,inplace=True)
用箱线图识别和去除异常值
在统计学中，异常值是与其他观察值显著不同的数据点。异常值可能是由于测量中的可变性造成的，或者它可能表示实验误差；后者有时被排除在数据集之外。异常值会在统计分析中引起严重的问题。(来源:维基百科)
所以基本上离群点就是不符合一般数据分布的数据点。数据中存在的异常值会扭曲模型并扰乱其学习过程。因此，现在我们将检测数据集中的这些异常值，并使用箱线图(四分位数间距)方法移除它们。要了解不同的异常值检测方法请点击这里。
sns.boxplot(x=df[‘price’])
价格箱线图
因此价格变量的异常值位于值 25000 之后。
sns.boxplot(x=df[‘ageOfCar’])
ageOfCar 箱线图
sns.boxplot(x=df['powerPS'])
powerPS 的箱线图
sns.boxplot(x=df['kilometer'])
公里箱线图
同样，对于特征而言，ageOfCar、powerPS、千米值分别高于 30、高于 280 和低于 25000 都是异常值。这些异常值现在已被移除。
df=df[(df.price<=25000) & (df.ageOfCar>=0) & (df.ageOfCar<=30) & (df.powerPS<=280) & (df.kilometer>=25000)]
编码分类变量
现在，我们需要对数据集的分类变量进行编码，以便模型可以使用它进行训练。分类变量有几种编码技术。在这种情况下，我们将对具有两个类的变量使用标签编码，对具有两个以上类的变量使用目标编码。目标编码优于一键编码，因为后者会导致数据维数非常高。
x=df.drop([‘price’],axis=’columns’,inplace=False)
y=df[‘price’]le1=LabelEncoder().fit(x['gearbox'])
x['gearbox'] =le1.transform(x['gearbox'])
le2=LabelEncoder().fit(x['notRepairedDamage'])
x['notRepairedDamage'] =le2.transform(x['notRepairedDamage'])
le3=LabelEncoder().fit(x['abtest'])
x['abtest']=le3.transform(x['abtest'])from category_encoders import TargetEncoderte=TargetEncoder(cols=['brand','model','vehicleType','fuelType','postalCode']).fit(x,y)
x=te.transform(x)
使用相关热图
通过使用关联热图，我们可以了解变量之间的关系。从热图中可以得出两个主要观察结果(不包括对角线元素):
- 与其他输入变量高度相关的输入变量。
- 与输出变量相关性低的输入变量。
现在让我们生成热图。
f, ax = plt.subplots(figsize=(11, 9))
cmap = sns.diverging_palette(220, 10, as_cmap=True)
x[‘price’]=y
sns.heatmap(x.corr(),cmap=cmap,square=True,annot=True)
从热图中，我们可以看到输入变量之间没有太大的相关性，但是， abtest 与 *price 的相关性很差。*因此 abtest 被移除。
x.drop([‘abtest’],axis=’columns’,inplace=True)y=x['price']
x.drop(['price'],axis='columns',inplace=True)
标准化，分为训练集和验证集
我们现在将变量标准化，并将其分为训练集(80%)和验证集(20%)。
sc=StandardScaler()
x=sc.fit_transform(x)X_train, X_val, y_train, y_val = train_test_split(x, y, test_size=0.2, random_state = 42)
使用 XGBoost 模型进行训练和测试
我们使用一个 XGBoost 模型来训练数据并预测输出。我们还使用 sklearn 模块的 GridSearchCV 来确定模型的最佳参数。使用 GridSearchCV 调优的参数有:
- n_estimators: 定义 XGBoost 模型使用的估计器(树)的数量。
- max_depth: 定义每棵树允许的最大深度，用于控制过拟合。
- eta: 这是模型的学习率，它缩小了节点的权重，使其更加稳健。
为了获得模型性能的无偏估计，我们使用了带有cv=5的 K 倍交叉验证方法。通过使用 RMSE 作为误差度量来识别最佳参数。因为我们需要找到产生最小 RMSE 的模型，我们将 RMSE 的负数传递给scoring参数。
import xgboost as xgb
from sklearn.model_selection import GridSearchCVparam={‘n_estimators’:[10,20,50,100,150,200],
‘max_depth’:range(3,11),
‘eta’:[0.05,0.1,0.15,0.2,0.25,0.3]}xgr=xgb.XGBRegressor()
gs=GridSearchCV(estimator=xgr,param_grid=param,scoring=’neg_root_mean_squared_error’,cv=5,verbose=3)
gs.fit(X_train,y_train)
参数搜索完成后，我们现在可以查看确定的最佳参数及其相应的指标得分。
print(gs.best_params_)
print(gs.best_score_)
eta=0.05、n_estimators=200、max_depth=7被确定为产生最低 RMSE 分数的参数。使用 GridSearchCV 确定的这些参数，我们构建最终模型并将其用于预测。对于回归问题，比如这个问题，误差度量比如 RMSE、MSE、MAE、R 等。被使用。这里，我们使用 RMSE 和 T21(决定系数)作为误差度量。XGBoost 的score()函数返回回归问题的 R 值。
from sklearn.metrics import mean_squared_errorxgr=xgb.XGBRegressor(eta=0.05,max_depth=7,n_estimators=200)
xgr.fit(X_train,y_train)
pred=xgr.predict(X_val)
print('RMSE: ',mean_squared_error(y_val,pred,squared=False))
print('R2Score: ',xgr.score(X_val,y_val))
RMSE 和 R 值表明模型的结果相当准确。现在，为了评估模型的性能，并了解数据分布的捕捉情况，我们从残差中生成图。
在回归分析中，因变量的观测值(y)与预测值(ŷ)之差称为残差 (e)。每个数据点有一个残差。(来源: Stattrek )
我们使用标准化残差(残差除以其标准差)来创建散点图和直方图。
def residual_plot(y_test, y_pred):
res = y_test — y_pred
fig, axs = plt.subplots(1,2,figsize=(30,10))
std_res = res/np.std(res)
axs[0].title.set_text(‘Scatter Plot of residuals’)
axs[0].set_xlabel(‘True Output’)
axs[0].set_ylabel(‘Residuals’)
axs[0].scatter(y_test,std_res)
axs[1].title.set_text(‘Histogram of residuals’)
sns.distplot(std_res, ax = axs[1]);
plt.show()residual_plot(y_val,pred)
残差的散点图和直方图
散点图显示存在一些异常值，并且残差不具有恒定的方差。因此，它们不遵循正态分布。直方图显示大部分残差集中在零附近。
特征重要性
为了理解输出预测过程中每个特征的重要性并得出未来的推论，我们绘制了直接从 XGBoost 模型中获得的特征重要性值。
feature_important = xgr.feature_importances_
values = list(feature_important)data = pd.DataFrame(data=values, index=feat, columns=[“score”]).sort_values(by = “score”)
ax=data.plot(kind=’barh’,title=’Feature Importance’,legend=None)
ax.set_xlabel(‘Importance value’)
ax.set_ylabel(‘Features’)
plt.show()
特征重要性
未修复损坏、燃油类型、和变速箱对价格的预测作用不大。
结论
总之，我们利用二手车数据集来理解数据清理、插补、特征选择和特征工程的概念。我们还学习了如何使用 XGBoost 模型，并通过绘制残差和特征重要性值对其进行回归分析。
要获得 IPython 笔记本的完整代码请点击这里。
用 R 和 Python 中的多元逻辑回归预测 SARS CoV-2 结果的概率
对 SARS 冠状病毒 2 型患者进行分类以及哪些变量会影响结果。
来源:维克斯
编者注: 走向数据科学 是一份以研究数据科学和机器学习为主的中型刊物。我们不是健康专家或流行病学家，本文的观点不应被解释为专业建议。想了解更多关于疫情冠状病毒的信息，可以点击 这里 。
简介
大约两个月前，巴西的以色列阿尔伯特·爱因斯坦医院发布了关于 Kaggle 的公开数据，以便更好地了解这种新型冠状病毒。因此，我试图使用这个数据集来预测患 SARS CoV-2 的概率，并检查哪些变量影响最具预测性的 SARS CoV-2 结果。
背景
截至 3 月 27 日，圣保罗州记录了 1223 例新冠肺炎确诊病例，68 例相关死亡，而拥有约 1200 万人口的圣保罗州，以及以色列阿尔伯特·爱因斯坦医院所在地，截至 3 月 23 日，有 477 例确诊病例和 30 例相关死亡。该州和圣保罗县都决定建立隔离和社会距离措施，这将至少持续到 4 月初，以减缓病毒的传播。
数据集
该数据集包含在巴西圣保罗的以色列阿尔伯特爱因斯坦医院就诊的患者的匿名数据，这些患者在就诊期间收集了样本以进行新型冠状病毒 RT-PCR 和其他实验室测试。
数据集总共有 5644 行和 111 列。
分析中使用的变量(预测值)
这篇文章不足以解释所有使用的变量。因此，我将即时讨论一些重要的预测因素。
为了准备数据，我使用了 Python 中的 Pandas。对于逻辑回归分析，我使用了 r。现在让我们跳到数据准备。
数据准备
第一眼看数据集，我发现了许多 NaN 值。在前进之前，必须纠正缺少值的问题。
按降序排列的每一列的值计数
从这里开始，我将数据集分成两个子集，命名为条件和疾病。(它们只是名字，不代表任何东西)橙色标记的是疾病子集，绿色标记的是病情子集。
变量中性粒细胞和蛋白质 C reativa 中的缺失值必须进行估算。我已经使用 KNNImputer 来估算丢失的值。
对于这两个数据集，我们的目标变量是 *SARS_CoV2_exam_result。*数据集疾病具有所有分类预测值，数据集条件具有所有连续数值预测值。
对于分类变量，0 表示未检测到疾病，1 表示检测到疾病。
分析
从这里开始，我将使用 R 进行分析。让我们从疾病数据集开始。
疾病数据集分析
疾病数据集摘要
我不得不跳过副流感 2 号，因为没有病例。现在，让我们绘制克拉姆夫相关图，以检查分类变量之间的相关性。
# CramerV Correlation to check for any correlation between catagorical variables
# Except Patient's age, all other variables are catagorical(binary)
>disease.corr = PairApply(disease[,names(disease) != "Patient_Age_Quantile"],cramerV, symmetric = TRUE)
# Displaying correlation with variable SARS_COV2_Result
>disease.corr[,2]
# Correlation plot
>corrplot(disease.corr, method = "square", type = "lower")
疾病数据集相关图
这里最重要的预测因子似乎是我们目标上的鼻病毒或肠道病毒，相关值为 0.1517。(不太重要)现在让我们将数据集分成 80%的训练和 20%的测试。
# Dividing Train/Test data with 80% training and 20% test
>sample_size <- floor(0.8 * nrow(disease))
>train <- sample(nrow(disease), size = sample_size)
>disease.train <- as.data.frame(disease[train,])
>disease.test <- as.data.frame(disease[-train,])
我们的数据集还没有准备好进行逻辑回归分析。
# Logistic regression considering top predictors
>disease.function = paste("SARS_COV2_Result", "~", "Patient_Age_Quantile + Rhinovirus_OR_Enterovirus")
disease.glm = glm(as.formula(disease.function), data = disease.train , family = binomial)
>summary(disease.glm)
疾病概述
在这里，预测因子患者年龄分位数和鼻病毒或肠道病毒的 P 值显示了巨大的希望。
根据上面的等式，我们可以说，如果一个人患有鼻病毒，他/她患 SARS CoV2 的可能性较小。
我使用了 10 重交叉验证来验证我们的模型。下一步是预测测试数据集上的值并准备混淆矩阵。
# 10 fold cross-validation to verify the model
>cv.glm(disease.train,disease.glm,K=10)$delta[1]# Predicting on test data based on training set
>disease.glm.predict <- predict(disease.glm,disease.test,type = "response")
>summary(disease.glm.predict)
# Mean predict results for SARS Cov2 Results diagnosis
>tapply(disease.glm.predict, disease.test$SARS_COV2_Result, mean)# Confusion matrix for threshold of 1%
>disease.confusion = table(disease.test$SARS_COV2_Result, disease.glm.predict > 0.01)
>disease.confusion# False negative error(Type II error)
>disease.type2error = disease.confusion[1,1]/
(disease.confusion[1,1]+disease.confusion[2,2])
>disease.type2error
我设置了一个严格的阈值(1%)来最小化第二类错误。该模型必须最大限度地减少误检患者的数量。60%的患有新冠肺炎的患者将从模型中被错误地检测为阴性，这是相当高的。
正如预测的那样，ROC 曲线是 okayish 的。
# Plotting ROCR curve
>disease.ROCRpred = prediction(disease.glm.predict, >disease.test$SARS_COV2_Result)
>disease.ROCRperf = performance(disease.ROCRpred, "tpr", "fpr")
>par(mfrow=c(1,2))
>plot(disease.ROCRperf, colorize=TRUE, >print.cutoffs.at=seq(0,1,by=0.1), text.adj=c(-0.2,1.7))# Probability of SARS vs Rhinovirus plot
>disease.predicted.data <- data.frame(
probability.of.having.SARS=disease.glm.predict,
Rhinovirus=disease.test$Rhinovirus_OR_Enterovirus)
>disease.predicted.data <- disease.predicted.data[
order(disease.predicted.data$probability.of.having.SARS, decreasing=FALSE),]
>disease.predicted.data$rank <- 1:nrow(disease.predicted.data)>plot(probability.of.having.SARS ~ Rhinovirus,disease.predicted.data, type = "l", lwd = 2)
ROC 图和 SARS 与鼻病毒的概率图
鼻病毒的症状接近 SARS CoV-2 病毒，多见于鼻部(鼻病毒)。因此，随着鼻病毒的检测，患 SARS 的概率降低。
条件数据集分析
条件数据集包含所有数字预测值，即从血液样本中获取的值。从摘要来看，数据看起来正常。
条件数据集摘要
重复相同的步骤，对目标变量的较高相关预测值应用逻辑回归。
# attaching the final dataset
>attach(condition)# Correlation between varables
>condition.corr = cor(condition)
>condition.corr[,2]
# Correlation plot
>corrplot(condition.corr, method = "square", type = "lower")# Dividing Train/Test data with 80% training dataset
>sample_size <- floor(0.8 * nrow(condition))
>train_ind <- sample(nrow(condition), size = sample_size)
>condition.train <- as.data.frame(condition[train_ind,])
>condition.test <- as.data.frame(condition[-train_ind,])# Logistic regression considering all the variables on the targer variable SARS_COV2_Result
>condition.function = paste("SARS_COV2_Result", "~", "Patient_Age_Quantile + Leukocytes + Eosinophils + Red_blood_cell_distribution_width_RDW + Platelets + Proteina_C_reativa_mg_dL")
>condition.glm = glm(as.formula(condition.function), data = >condition.train , family = binomial)
>summary(condition.glm)# 10 fold cross-validation to verify the model
>cv.glm(condition.train,condition.glm,K=10)$delta[1]# Predicting on test data based on training set
>condition.glm.predict <- predict(condition.glm,condition.test,type = "response")
>summary(condition.glm.predict)
>tapply(condition.glm.predict, condition.test$SARS_COV2_Result, mean)# Confusion matrix for threshold of 1%
>condition.confusion = table(condition.test$SARS_COV2_Result, condition.glm.predict > 0.01)
>condition.confusion# False negative error rate (Type II error)
>condition.type2error = condition.confusion[1,1]/ (condition.confusion[1,1]+condition.confusion[2,2])
>condition.type2error
条件数据集相关图
这里我们看到一些高度相关的变量。然而，我们只考虑与目标变量高度相关的预测值。因此，为了进行分析，我将使用患者年龄分位数、白细胞、嗜酸性粒细胞、红细胞分布宽度 RDW、血小板和蛋白质 a_C_reativa_mg_dL 预测值。我们来看看 glm 总结。
条件数据集 glm 摘要
让我们预测一下测试数据集的误差。使用该模型，121 例中的 30 例给出错误的结果，并且 ROC 曲线优于先前的模型。
# Plotting ROCR curve
>condition.ROCRpred = prediction(condition.glm.predict, condition.test$SARS_COV2_Result)
>condition.ROCRperf = performance(condition.ROCRpred, "tpr", "fpr")
>plot(condition.ROCRperf, colorize=TRUE, print.cutoffs.at=seq(0,1,by=0.1), text.adj=c(-0.2,1.7))
条件数据集 ROC 曲线
绘制结果，
# Creating a dataframe with variables and predicted values of SARS results
>condition.predict.dataframe <- data.frame(
probability.of.having.SARS=condition.glm$fitted.values,
Leukocytes=condition.train$Leukocytes,
Patient_Age_Quantile = condition.train$Patient_Age_Quantile,
Eosinophils = condition.train$Eosinophils,
Red_blood_cell_distribution_width_RDW = >condition.train$Red_blood_cell_distribution_width_RDW,
Platelets = condition.train$Platelets,
Proteina_C_reativa_mg_dL = >condition.train$Proteina_C_reativa_mg_dL)>plot1 = ggplot(data=condition.predict.dataframe, aes(x=Patient_Age_Quantile, y=probability.of.having.SARS)) +
geom_point(aes(color=Patient_Age_Quantile), size=4)
>plot2 = ggplot(data=condition.predict.dataframe, aes(x=Leukocytes, y=probability.of.having.SARS)) +
geom_point(aes(color=Leukocytes), size=4)
>plot3 = ggplot(data=condition.predict.dataframe, aes(x=Red_blood_cell_distribution_width_RDW, y=probability.of.having.SARS)) +
geom_point(aes(color=Red_blood_cell_distribution_width_RDW), size=4)
>plot4 = ggplot(data=condition.predict.dataframe, aes(x=Eosinophils, y=probability.of.having.SARS)) +
geom_point(aes(color=Eosinophils), size=4)
>plot5 = ggplot(data=condition.predict.dataframe, aes(x=Platelets, y=probability.of.having.SARS)) +
geom_point(aes(color=Platelets), size=4)
>plot6 = ggplot(data=condition.predict.dataframe, aes(x=Proteina_C_reativa_mg_dL, y=probability.of.having.SARS)) +
geom_point(aes(color=Proteina_C_reativa_mg_dL), size=4)
# Plotting the values
>grid.arrange(plot1, plot2, plot3, plot4, plot5, plot6, ncol=3 , nrow = 2)
患 SARS 的概率与预测值
我们可以看到，患者的年龄分位数和蛋白 a C reativa 呈正相关，其余与目标变量的概率呈负相关。
预测器的单独调查
一些预测因素在单独分析时显示出巨大的潜力。这些预测是，
我现在将创建一个函数来分析这些变量，以查看它们的效果。
# Plotting Seperate probability graphs
>plotting.function <- function(var,variableORrank){
condition.sep.function = paste("SARS_COV2_Result", "~", as.character(var))
condition.sep.glm = glm(as.formula(condition.sep.function), data = condition.train , family = binomial)
print(summary(condition.sep.glm))
cv.glm(condition.train,condition.sep.glm,K=10)$delta[1]
condition.predicted.data <- data.frame(
probability.of.SARS=condition.sep.glm$fitted.values,
variable=condition.train[,as.character(var)])
condition.predicted.data <- condition.predicted.data[
order(condition.predicted.data$variable, decreasing=FALSE),]
condition.predicted.data$rank <- 1:nrow(condition.predicted.data)
ggplot(data=condition.predicted.data, aes(x= variable, y=probability.of.SARS)) +
geom_point(aes(color=variable), size=3) +
xlab(as.character(var)) +
ylab("Probability of having SARS CoV-2") +
scale_colour_gradient(low = "darkgreen", high = "darkred", na.value = NA) +
ggtitle(coef(summary(condition.sep.glm))[,'Pr(>|z|)'])
}>plotfun1 = plotting.function("Platelets")
>plotfun2 = plotting.function("Monocytes")
>plotfun3 = plotting.function("Hemoglobin")
>plotfun4 = plotting.function("Red_blood_cells")
>plotfun5 = plotting.function("Mean_platelet_volume")
>grid.arrange(plotfun1, plotfun2, plotfun3, plotfun4, plotfun5, ncol=3 , nrow = 2)
结果是，
患 SARS 的概率与单一预测因子分析(图表标题为 Pr(>|z|)值，越低越好)(单核细胞是一种白细胞，血红蛋白是一种含铁的氧运输金属蛋白)
除了血小板之外，所有预测因子都与患 SARS CoV-2 的概率呈正相关。
结论
对于医院来说，区分鼻病毒携带者和新冠肺炎病毒携带者是非常重要的。将鼻病毒和新冠肺炎病毒患者分开会让这个过程快很多。
高%vol 的白细胞表明身体正在对抗病毒。
当携带病毒时，身体不能产生足够的血小板来帮助止血。
利用给定的患者血液数据，该模型可以预测患 SARS CoV-2 的概率。
分析中使用的代码
介质物品分析代码。通过创建一个……为零熊猫/新冠肺炎回归分析的发展做出贡献
github.com](https://github.com/ZeroPanda/COVID-19-Regression-Analysis)
参考文献
统计学习介绍及在 R 中的应用；加雷斯·詹姆斯，丹妮拉·威滕，特雷弗·哈斯蒂，罗伯特·蒂布拉尼
支持临床决策的人工智能和数据科学(3 月 28 日至 4 月 3 日)
www.kaggle.com](https://www.kaggle.com/einsteindata4u/covid19) [## 维基百科(一个基于 wiki 技术的多语言的百科全书协作计划ˌ也是一部用不同语言写成的网络百科全书ˌ 其目标及宗旨是为全人类提供自由的百科全书)ˌ开放性的百科全书
维基百科是一个免费的在线百科全书，由世界各地的志愿者创建和编辑，由维基媒体托管…
www.wikipedia.org](https://www.wikipedia.org/)
预测特许学校的弹性
我开发工具帮助家庭在教育选择中感到安全的经历
图片来源:Neonbrand，via Unsplash ( Unsplash 许可)
特许学校对学生来说是一个很好的选择，但是关于特许学校稳定性的问题可能会阻止一些家庭给他们机会。在 Insight Data Science，我开发了一个工具来帮助加州的家庭了解上一所新开的特许学校的风险。
特许学校是公共资助但私人经营的，所以它们不是由传统的公立学区经营的，但它们不像私立学校那样收取学费。特许学校在过去 25 年中经历了大幅增长，2016 年全国约有 300 万学生入学(占公立学校总入学人数的 6%)。特许学校通常专注于特定的课程，如 STEM(科学、技术、工程和数学)或艺术，并且通常比传统的公立学校规模小。这些特点可能会吸引那些寻找更传统教育环境的替代方案的学生，并在一定程度上解释了特许经营行业的增长。
来源:国家教育统计中心
上特许学校的好处和风险
当我在亚利桑那大学攻读经济学博士学位时，我有机会为 BASIS Schools(一个特许学校网络)进行定量研究。在 BASIS 学校，我观察了课堂教学和其他各方面的管理，对网络提供的教育体验印象深刻。
BASIS Schools 从不起眼的地方起步:它的第一所学校于 1998 年在一个租赁的地方开办，那里以前是一个杂货店。BASIS Schools 在《美国新闻与世界报道》2019 年全国最佳公立高中排名中增加了两所高中，包括排名第三的学校 BASIS Scottsdale。
对于一个从以前的杂货店起步的学校系统来说，这已经不错了！来源:美国新闻与世界报道
不幸的是，并不是所有特许学校的故事都以同样的成功告终。由于特许学校通常缺乏传统公立学校管理的学校的运营和财政支持，它们可能面临更大程度的不稳定。
自 2010 年以来，加州开设了 831 所特许学校，其中 73.5%仍在运营(相比之下，同期开设的加州传统公立学校的比例为 90.4%)。意外的学校关闭可能会给家庭带来严重的痛苦，让家长和学生争相寻找合适的教育替代方案。
一位母亲的女儿就读于一所意外关闭的新特许学校，她说她会考虑再次就读一所特许学校，但想知道:
“…我怎样才能更好地评估每个人的数字？…我们需要问的难题是什么？因此，如果我们能够在未来防止这种情况发生，我们绝对会这样做。”
八岁的麦肯齐·惠斯勒的圣诞假期开局不利。她刚刚发现她不会…
www.wunc.org](https://www.wunc.org/post/what-happens-when-charter-school-closes)
介绍… CharterPredictor
了解到进入一所新的特许学校的潜在积极和消极方面，我创建了 CharterPredictor。在 charterpredictor.com 的，加州的家庭可以找到一所 2017 或 2018 学年新开的特许学校，并在学校上运行 CharterPredictor。CharterPredictor 随后为家长提供了一个衡量学校是“风险”还是“安全”选择的标准，以及学校在最初开放后继续运营三年的可能性。
CharterPredictor 还显示学校的初始注册人数，这是一个预测学校是否会继续运营的重要指标(学校资金与注册学生人数挂钩，如果注册学生人数不足，就很难覆盖固定成本)。CharterPredictor 还提供与感兴趣的学校相邻的新特许学校的附加信息，以便有兴趣支持新特许学校的家庭可以快速了解其他相关选项。
背景图片来源:Moren Hsu，via Unsplash ( Unsplash 许可)
虽然最终我们对所有特许学校的稳定性感兴趣，但对于这个项目的迭代，我想重点关注新的特许学校，因为这些学校特别容易关闭，因为它们往往缺乏更成熟的学校的资源。
CharterPredictor 背后的数据
CharterPredictor 背后的数据来自加州教育部，其中包括学校开学和关闭的信息，入学和人口统计信息。我收集了 2004 至 2018 学年的记录，然后将它们与美国社区调查中关于失业和教育水平的邮政编码水平数据结合起来。最后，我获得了与每所特许学校相邻的传统公立学校的家长评级(来自 great schools——great schools 就像 Yelp for schools，家长可以在这里留下对学校的评论和评级),以便衡量传统公立学校的竞争。
为了预测哪些新的特许学校有关闭的风险，我收集了 2004 年至 2016 年期间加州所有新的特许学校(包括 1017 所新学校)，并使用 60/40 的拆分将数据分为训练集和测试集。CharterPredictor 然后使用 2017 年和 2018 年开设的学校的信息进行预测。
模拟学校关闭
我使用 XGBoost 来模拟特许学校的关闭，因为 XGBoost 给我的预测比逻辑回归模型更准确。XGBoost 模型导致曲线下面积(“AUC”)测量值为 0.95，而逻辑回归导致 AUC 为 0.94。虽然逻辑回归更容易理解，但我继续使用 XGBoost，因为我希望它能更好地推广到加利福尼亚州以外的新的州和新的功能。
我在几个超参数上使用交叉验证和网格搜索，包括最大树深度、收缩率和提升轮数来调整模型(您可以在 CharterPredictor GitHub 资源库查看构建分析数据集的代码)。预计关闭概率为 20%或更高的特许学校被归类为“有风险的”，因为如果我是一名家长，考虑到学校的稳定性，这似乎是我希望得到通知的水平，尽管这个截止值无疑有些武断。
在评估测试集时，经过训练的模型具有 94.1%的准确性，70.3%的召回率和 70.3%的精确度。在我们测试的 38 所实际关闭的学校中，CharterPredictor 可以正确地将其中的 26 所归类为风险学校。如果 CharterPredictor 的主要目标是保护学生免受不稳定因素的影响，我们可能会选择更低的关闭概率作为阈值来定义一所学校有风险。然而，重要的是避免过度热衷于将学校分类为风险学校，因为其中一些学校可能是我们希望家庭考虑的好机会。
结论和下一步措施
由于教育数据和特许学校政策因州而异，我想先从加州开始。下一步要使它成为更有用的工具，就要把它扩展到加州以外的地方。到目前为止，我相信 CharterPredictor 可以成为一个有用的工具，帮助加州家庭对他们的教育选择更有信心。
总的来说，我喜欢与跨越私人和公共领域的行业合作(如教育、卫生、能源)，因为它们往往有大量可用的数据(我可能很难获得完全私人市场的企业细节和开业/倒闭)。此外，这些公共/私人市场往往与我们从社会角度非常关心的结果联系在一起。
关于作者
帕特里克·西泽克在亚利桑那大学获得经济学博士学位，专攻经验微观经济学。Patrick 目前是 Esurance 的数据科学家。
在 Esurance 工作之前，Patrick 是 Insight Data Science 的数据科学研究员。访问 Insight Data Science 的网站。
预测推迟的 2020 年橄榄球六国赛的比分
访问和分析六国橄榄球数据，以及预测剩余延期比赛的分数
介绍
六国橄榄球锦标赛是一年一度的橄榄球联盟比赛，由英格兰、法国、爱尔兰、意大利、苏格兰和威尔士的国家队参加。锦标赛定于 2020 年 3 月 14 日结束；然而，由于新冠肺炎疫情，意大利对爱尔兰的倒数第二场比赛和所有的决赛都被推迟，意在改期。这是自 2001 年爆发口蹄疫以来，首次有超过一场比赛被推迟。鉴于此，威尔士橄榄球球迷and 的数据科学家对冠军赛的相关数据以及利用这些数据可能实现的目标感到好奇。这是我写的一篇关于我对六国橄榄球数据的分析，以及预测 2020 年推迟比赛的分数的文章。
六国橄榄球数据
首先要克服的问题是获得一些我可以用来模拟橄榄球比赛的数据。鉴于锦标赛始于 1882 年，我很早就决定要调查每一场可能的比赛，这需要花费相当多的时间和精力。随着锦标赛在 2000 年正式成为六国赛，我转向谷歌，开始寻找 10 年的橄榄球数据。
在网上查看后，我偶然发现了 ESPN Scrum 网站，它包含了从 1882 年以来所有六个国家比赛的基本数据。
ESPN Scrum 网站
最突出的一点是数据是表格形式的。当我看到这张写着“的桌子”时，我不禁暗自高兴。表格是我们的朋友，它让处理数据变得更干净，更轻松。也就是说，只有当数据源结构良好时，这才是正确的。
让我们看看这个网站还能提供什么。当你点击冠军年，在这种情况下，2020 年，你会被重定向到另一个页面，其中包含了迄今为止的结果表。我们可以看到，爱尔兰和意大利、威尔士和苏格兰、意大利和英格兰、法国和爱尔兰之间的最后四场比赛由于延期，因此没有与之相关联的比分。这些是我们想要预测的。
2020 年橄榄球六国记分牌至今
这个阶段需要考虑的重要的是可用的数据；也就是说，这些变量是否对匹配建模有用。我们目前知道的是，对于每支球队来说，比赛是在主场还是客场进行的，比分，以及球队是赢了还是输了比赛。在这个阶段，这些属性可能足以做出合理的预测。然而，关注球队层面的变量，比如每支球队的出场总数，以及球员层面的变量，比如他们的身高、体重、最近的表现等等，会很有趣。如果我们可以找到一个可靠的数据资源，并且数据容易访问，我们可能会在以后添加这些。但是首先，让我们从 ESPN 网站上摘录这些信息。
网页抓取
ESPN 网站不提供任何类型的可下载数据，所以是时候深入网页的来源，提高我的网络抓取技能了！
这个网站的好处在于，2000 年至 2020 年之间的每个记分牌的 HTML 表格都使用相同的 ID。这意味着，给定所有 20 个页面，我可以遍历 URL 并应用相同的预处理步骤来访问每个表。
首先，我复制了 2000 年到 2020 年之间的所有记分牌 URL，并将它们存储在一个数组中。然后，我遍历数组，发出 GET 请求从每个页面获取原始 HTML 内容。使用 BeautifulSoup，HTML 内容被解析为 Python，这样我就可以使用 ID scrumArticleBoxContent找到特定的记分板表。
现在我有了所有字符串形式的 HTML 表，我可以使用read_html函数将其解析为 Pandas。因为这个函数输出一系列数据帧，所以我可以使用concat函数将它们连接成一个主数据帧。然后，我删除了通过 HTML 潜入的不相关的列。
一些行的日期为空，这也反映在 ESPN 网站上显示的表格中。通过使用ffill函数，我可以用该列中最后一个已知的日期来填充那些缺失的日期。
为了预测每个国家在推迟的比赛中取得的分数，我们将不得不分开每个比赛记录，以便每个队在单独的一行中得到代表。在此之前，我们还必须考虑到国名的顺序代表了哪个国家是主场还是客场比赛。由于球队的名字有一个破折号作为固定的分隔符，我们可以根据破折号分割字符串，将左边的名字解析到home_team列，将右边的名字解析到away_team列。
现在队伍被分开了，我们也必须从队伍名称中删除分数。我们可以使用 Regex 从团队名称中提取数字，并将它们解析到home_team_score和away_team_score列中。
我们还想知道比赛场地的名称。鉴于这些数据在另一页的另一个表格中，我将打一张懒牌，假设在过去 10 年里，这些比赛在每个国家的同一个体育场进行。我可以通过迭代球场名称并从我的球场字典中分配相应的值，将球场名称与主队联系起来。
基于这最后一个数据帧，我们知道如果home_team_score不为空，游戏就开始了，如果为空，那么这些就是被推迟的。我们可以根据这个规则将数据分为训练和测试。
让我们先关注一下训练数据。在这个博客的下一部分，我需要知道哪个国家赢了这场比赛。这不会作为一个变量包含在游戏比分预测中。为此，我遍历训练数据帧，并将获胜的团队关联到一个winner列。
最后，为了分割数据，使每个国家都是数据框中的一行，我可以根据国家是主场还是客场来选择这些国家，然后将它们附加到另一个国家的顶部。
快速分析
这篇博客的主要焦点是看我们是否能做出合理的比赛分数预测。然而，我认为快速查看数据中是否存在任何潜在模式也很有趣。我感兴趣的一件事是，与客场比赛相比，一个球队在主场比赛时是否会赢得更多比赛，这两者之间是否存在某种关联。
围绕主场优势有很多卓有成效的研究。也就是说，有没有据说主队比客场队获得的好处？这些好处可以归因于支持球迷对竞争对手或裁判的心理影响、时区或气候的差异、旅行后的疲劳以及许多其他因素。所以，我想知道 10 年来 6 个国家的数据是否反映了这一点。
熊猫的correlation函数用于查找数据帧中各列之间的成对相关性。因为这个函数忽略任何非数字数据，所以必须首先将它们映射到某种数字表示。熊猫有另一个内置功能可以为你处理这些。当您将列的数据类型转换为category时，该列会将非数字数据转换为分类表示。因此，在这种情况下，如果我们要查看一支球队在主场或客场比赛时与他们输赢之间是否存在相关性，例如，这些列将把win转换为 1，把lose转换为 0。
如果我们对整个数据集应用correlation函数，我们会得到 12%的相关分数。相关性旨在最大化其价值；也就是说，较高的分数意味着属性之间有很强的关系。所以在这种情况下，12%的低得分表明没有太大的主场优势。
现在，我不是 100%确信这是真的。当我的祖国威尔士在卡迪夫千年体育场主场比赛时，我确信当威尔士观众放声高唱卡隆·兰时，这点燃了威尔士橄榄球队的某种激情。我自己想想都起鸡皮疙瘩！所以，我们来看看这个。
如果我们关注威尔士在 10 年六国赛中的表现，我们可以绘制出一些结果。首先，我们来看看他们的总体表现。10 年来，威尔士队参加了 100 场锦标赛。在这 100 场比赛中，威尔士赢了 52 场，输了 45 场，平了 3 场。
如果我们看看威尔士主场 49 场比赛的结果，他们赢了 28 场，输了 20 场，平了 1 场。
最后，如果我们看看威尔士 51 场客场比赛的结果，他们赢了 24 场，输了 25 场，平了 2 场。
鉴于这些结果，可以有把握地说correlation指标反映了相同的输出。与客场作战相比，威尔士主场作战时只多赢了 4 场比赛。但这并没有抵消我们的主场优势。这只告诉我们这个数据集内的相关性。我们可以将其他几场橄榄球比赛纳入分析，这些比赛可能会改变这一结果。
进行分数预测
将主场优势放在一边，我想看的下一件事是预测剩下的延期比赛的比分。有两种监督学习方法。回归用于预测连续值，而分类预测离散输出。例如，预测房子的英镑价格是一个回归问题，而预测肿瘤是恶性还是良性是一个分类问题。由于橄榄球比赛的最终得分在技术上可以是任何正数(甚至是零)，我们将研究回归方法。
首先，让我们为训练和测试添加一个标签，这样我们就可以知道哪个数据集是哪个数据集。然后我们将它们组合起来，并对team, stadium和home_or_away列进行编码。这是为了确保在我们再次分割它们之后，我们在每个集合中有相同的列和编码。这里需要注意的是，我没有将date作为建模中的一个变量。现阶段我只是希望能够根据前面几场的比分来预测剩下几场的比分。此外，由于锦标赛通常在每年的同一时间举行，我不认为这对造型有多大影响。如果是季节性的，这可能是一件有趣的事情。
一旦我们对数据进行了编码并再次将其拆分，这样我们就有了训练集和测试集，我们将需要再次拆分我们的训练数据，这样我们就有了验证集。也就是说，我们将从训练集中抽取一个数据样本，这样我们就可以评估我们的模型在预测这些样本的分数方面的性能，因为我们已经知道答案。为此，我们将抽取 30%的训练集样本。所以在 610 个记录中，有 183 个样本。我们还将确定哪些列是属性，哪些是标签。属性是自变量，而标签是我们想要预测其值的因变量。现在我们已经对数据进行了编码，我们有 14 个属性，我们希望根据这些属性来预测分数。因此，我们的属性被设置为 X 变量，而score列被设置为 y 变量。
现在来训练一个模型。这个问题虽然是回归问题，但不是线性问题。并不指望每支球队在每场比赛后得分都会增加。这就抵消了那些方法。在这种情况下，我决定使用一个 RandomForestRegressor ，因为该算法易于使用，相对准确，而且与标准决策树相比，它可以减少过拟合。随机森林算法创建了几个决策树，在特征权重中注入了一些随机性。这些决策树然后被组合以创建一个用于最终分析的森林(因此是决策树 ) 的随机森林)。该算法既支持分类，也支持回归，使其非常灵活地用于各种应用。
在构建和训练我们的模型之前，我们首先需要设置一些超参数。这些参数通常是最难设置的，因为这些设置通常没有完美的值。一般的经验法则是最初坚持使用默认值，然后一旦模型被训练和测试，就开始使用试错法来调整这些值，直到获得最佳结果。对于这个模型，我发现 50 个 n 估计值和 4 个 max_depth 是提供最佳结果的一组很好的值。有关这些特定设置的更多详细信息，可在官方 scikit-learn 文档中找到。
我将在X_train(属性)和y_train(标签)上拟合模型，然后预测我们从训练中移除的验证样本X_val 的得分y_pred。然后，我会通过测量均方根误差(RMSE)来比较预测值和实际得分(T4)。
RMSE 是残差方差的平方根。它表示模型与数据的绝对拟合。也就是说，观察到的数据点与模型的预测值有多接近。它还有一个有用的特性，就是和你要预测的变量在同一个单位里；所以在这种情况下。
模型的 RSME 结果是 10.1 (10)分…好吗？不好？绝对不完美。由于大多数球队每场比赛得分不到 40 分，10 分的误差并不是特别准确。正如这里所展示的，该模型在预测所有结果的分数方面做了相当多的工作。在预测得分较低的游戏时，该模型似乎表现得更好，但不是明显更好。也许如果我们把那些高分匹配从图片中去掉，我们可能会看到模型中的错误减少。但我们希望捕捉这些可能性，不会把它们视为离群值。这是需要进一步调查的事情。
让我们看看这个模型对延期比赛的预测得分。我们可以在X_test上应用模型，并将预测作为一列附加到测试数据帧上，这样更容易阅读。
好的。模型还没有给出不切实际的预测！预计英格兰会以 36 比 13 战胜意大利。虽然这么说让我很痛苦，但考虑到英格兰目前领先赢得整场比赛，这似乎是一个令人尊敬的预测。爱尔兰看起来像预测他们会输给法国 19 比 25。同样，这似乎不是一个不切实际的预测。在每场比赛以大约 10 分的优势击败英格兰和意大利之后，法国和英格兰并列榜首。预测威尔士 25 比 16 赢苏格兰(耶！🏴).在只赢了意大利之后，这意味着我们将会避免那个木勺！最后，预计爱尔兰将以 25 比 11 战胜意大利。
结论
那么，我从这个分析中学到了什么？
非线性回归问题并不总是容易的。虽然模型的 RSME 没有我希望的那么低，但我怀疑这是因为数据的可变性。有时威尔士在球场上全力以赴，以超过 40 分的比分获胜，但其他时候，他们真的没有！这是我想进一步调查的事情。我还想在预测中加入其他变量，比如天气状况，裁判是谁，因为他们总是为威尔士的失利负责…还有更多！
当前橄榄球六国 2020 领先委员会
然而，该模型预测的分数并没有那么不切实际。根据 2020 年锦标赛的分数，预测的产出遵循类似的趋势。总体而言，如果没有奖励积分，法国和英格兰可能会成为 2020 年橄榄球六国赛的共同赢家。
完整的笔记本，请查看下面我的 GitHub 回购:https://github.com/LowriWilliams/Rugby_Six_Nations_2020
如果你喜欢这篇文章，别忘了点赞和分享。
利用网络预测 R-1 中的新冠肺炎病例
激励和创建网络
第二部分可以在这里找到。
像所有传染病一样，新冠肺炎病毒的传播发生在非常局部的水平上。它不只是在世界各地随机跳跃，它必须在人与人之间传播，这些人必须旅行到一个新的地方传播病毒。正是因为这个原因，网络经常被用来帮助模拟疾病的传播，无论是通过使用社交网络来模拟疾病在人与人之间的传播，还是通过地理网络来帮助我们根据邻近地方发生的事情预测疾病的爆发。
在这个博客系列中，我们将关注地理网络，其中节点是地点，如果两个地点相邻，则边连接这两个地点。最终，我们将尝试使用这样一个网络，根据邻近地区的 Covid 病例来预测一个地区的 Covid 病例。
这种方法的准确性将受到限制，因为它没有考虑流行病学因素(如群体免疫和社会距离)或其他地区联系方式(如火车)。然而，这个博客系列仍然为可视化网络和使用它们创建预测分析的特征提供了有益的指导。它还表明网络分析可能是更大的 Covid 模型的有用部分。
(图片由作者创作)
在英国，发布 Covid 案例的最小地理区域被称为中间层超级输出区域(MSOAs ),因此我们将重点关注这些区域。你可以在这里的一个互动地图里看到这个数据，在这里下载。
在这第一部分中，我们将创造一些形象化的东西来激发为什么这种方法可能有效，然后我们将创造出上图中看到的网络。我们将展示代码，同时，所有代码都可以在 github 上找到。
这些都是用来创建本帖中显示的可视化的包。
# Set up ----
# Load packages
library(dplyr)
library(readr) # Loading the data
library(tidyr)
library(sf) # For the maps
library(sp) # Transform coordinates
library(ggplot2)
library(viridis)
library(gganimate) # For the animated map# Network packages
library(igraph) # build network
library(spdep) # builds network
library(tidygraph)
library(ggraph) # for plotting networks
动机
为了了解这种方法可能奏效的原因，我们将看看莱斯特，该市 6 月份的 Covid 病例激增，导致英国首次局部封锁。我们将通过使用{ggplot2}和{sf}包创建地图来可视化 Covid 在该区域的分布，并使用{gganimate}包制作动画。
我们需要下载案例数据(可以直接从 URL 下载)，并将其转换成正确的格式。初始表很宽(每个 MSOA 一行，每周一列)，有些周没有数据，所以我们利用相对较新的{tidyr}和{dplyr}函数where、pivot_longer和across来清理这些数据。
注意:自从这篇文章发表后，数据和访问数据的位置都发生了变化。这里的新网址是。我们在第 2 部分中使用这些新数据，你可以在 github 上访问新代码。
url <- '[https://c19downloads.azureedge.net/downloads/msoa_data/MSOAs_latest.csv'](https://c19downloads.azureedge.net/downloads/msoa_data/MSOAs_latest.csv')
msoa <- readr::read_csv(url, col_types = cols()) %>%
# Drop columns that are all NA
select(where(~!all(is.na(.)))) %>%
# Pivot longer
pivot_longer(dplyr::starts_with("wk_"),
names_to = "week",
values_to = "cases",
names_prefix = "wk_") %>%
# Turn week to numeric
mutate(week = as.integer(week)) %>%
# Turn NAs into 0s
mutate(across(c(last_7_days, cases),
.fns = ~ifelse(is.na(.), 0, .)))
生成的表如下所示(为了简化，我去掉了前几列)。前两列告诉我们地方当局辖区(LAD ),其面积比 MSOA 还大。接下来的两列告诉我们 MSOA。最后两列告诉我们那一周，以及那一周有多少阳性病例。“last_7_days”列告诉我们从上一整周到今天的阳性病例数。
我们应该做的第一件事是检查莱斯特峰值是否真的出现在数据中，这可以通过绘制简单的线图来完成。
msoa %>%
group_by(lad19_nm, week) %>%
summarise(cases = sum(cases))%>%
filter(lad19_nm == "Leicester") %>%
ggplot(aes(x = week, y = cases)) +
geom_line(size = 1) +
theme_minimal() +
ggtitle("Time series gaph of cases for Leicester")
(图片由作者创作)
我们看到数据中有一个很大的峰值，在第 25 周达到 500 例。
下一步是看看这是如何在 MSOAs 之间分配的。我们可以将它作为另一个线图来做(事实上，我们在 github 代码中也是这样做的)，但这不会给我们与在地图上绘制相同的洞察力。
要绘制地图，我们必须下载英国 MSOAs 的形状文件，你可以在这里做。保存并解压缩整个文件夹后，我们可以通过以下方式加载它:
msoa_sf <- sf::st_read("data/maps/Middle_Layer_Super_Output_Areas__December_2011__Boundaries_EW_BFE.shp")
绘制莱斯特的 MSOAs 非常简单，只需将形状文件限制在莱斯特的文件中(我们通过内部连接来实现)，然后向 ggplot 添加一个geom_sf()。为了绘制病例数据，我们限制在特定的一周(让我们选择第 25 周的峰值)并加入病例数据。现在，我们可以根据第 25 周的阳性病例数来填充 MSOA 地区。
msoa_data_wk25 <- msoa %>%
filter(lad19_nm == "Leicester",
week == 25)# We see that covid cases are clustered in several MSOAs
msoa_sf %>%
inner_join(msoa_data_wk25 %>% distinct(msoa11_cd, cases),
by = c("MSOA11CD" = "msoa11_cd")) %>%
ggplot()+
geom_sf(aes(geometry = geometry, fill = cases)) +
# Everything past this point is just formatting:
scale_fill_viridis_c() +
theme_void() +
ggtitle(label = "MSOA map for Leicester",
subtitle = "The colour shows Covid cases during week 25")
(图片由作者创作)
这是我们第一次看到网络有用的迹象。我们看到邻近的澳门特别行政区也有类似的情况。这是有道理的。如果你在一个地方爆发疫情，你预计它会蔓延到邻近地区。
像这样可视化数据的一个缺点是，很难看到案例如何随时间变化。这就是{gganimate}的用武之地。{gganimate}是一个很棒的软件包，可以让您根据特定的变量来制作图形动画。它允许我们为每个星期创建上面的地图，然后为我们把它组合成一个 gif。
# gganimate
# Need to install transformr for this to work
p <- msoa_sf %>%
inner_join(msoa %>% filter(lad19_nm == "Leicester"),
by = c("MSOA11CD" = "msoa11_cd")) %>%
ggplot(aes(group = week)) +
geom_sf(aes(fill = cases)) +
scale_fill_viridis_c() +
transition_time(week) +
labs(title = paste0("New covid cases for MSOAs within Leicester"),
subtitle = "Week = {frame_time}") +
theme_void()# We define this to help us pick the right number of frames
num_weeks <- n_distinct(msoa$week)animate(p, nframes = num_weeks, fps = 1, end_pause = 4)
我们看到，到transition_time()为止，这段代码看起来与我们上一张图表非常相似。唯一的区别是我们没有将数据过滤到特定的一周，我们在图的aes()中包含了一个group = week。transition_time()是我们如何指示它使用week作为变量来制作动画。
最后一点代码只是帮助我们指定一个好的帧数，这样每周都会得到相同的帧数。结果如下。
邻近地区的重要性再次凸显出来。高病例数从一个地区蔓延到邻近地区。此外，就在单个 MSOA 病例数达到峰值之前(发生在第 27 周)，我们看到邻近地区的病例数也很高。网络分析能让我们预测到这个特殊的区域在第 27 周会有这样的病例高峰吗？只有一个办法可以知道。
网络
现在我们知道为什么网络可以帮助我们，让我们实际上建立一个。要定义一个网络，知道所有的节点是什么，以及它们如何连接在一起(即边)就足够了。
在我们的例子中，节点很简单，这些只是 MSOAs。我们知道，如果两个节点所代表的两个 MSOAs 在我们的地图上是邻居，那么我们希望这两个边是连接在一起的。幸运的是，来自{spdep}包的poly2nb()函数(读“多边形到邻居”)帮助我们做到了这一点。它可以生成一个形状文件，并告诉我们哪些区域是邻居。
# First we just need to restrict to Leicester MSOAs for one week.
leicester_sf <- msoa_sf %>%
inner_join(msoa %>% filter(lad19_nm == "Leicester", week == 27) ,
by = c("MSOA11CD" = "msoa11_cd"))# Use the poly2nb to get the neighbourhood of each area
leicester_msoa_neighbours <- spdep::poly2nb(leicester_sf)
然后我们可以使用nb2mat()函数(读作“矩阵的邻居”)把它变成一个邻接矩阵。对于网络，邻接矩阵是这样的:如果存在连接节点 i 到节点 *j，*的边，则单元 i，j 等于 1，否则等于 0。因此邻接矩阵定义了一个网络，事实上，我们可以使用{igraph}包中的graph_from_adjacency_matrix函数来创建一个网络对象。
# Use nb2mat to turn this into an adjacency matrix
adj_mat <- spdep::nb2mat(leicester_msoa_neighbours, style = "B")
rownames(adj_mat) <- leicester_sf$msoa11_hclnm
colnames(adj_mat) <- leicester_sf$msoa11_hclnm# Use graph_from_adjacency_matrix to create a graph object from the adjacency matrix
leicester_network <- igraph::graph_from_adjacency_matrix(adj_mat, mode = "undirected")
您可以就此打住，但我发现{tidygraph}包提供了一种使用标准{dplyr}动词(我们将在后面看到)操作网络对象的惊人方法。所以我们用as_tbl_graph()把这个网络对象变成一个整齐的图形对象。
# Use as_tbl_graph to turn this into a tidygraph
leicester_network <- igraph::graph_from_adjacency_matrix(adj_mat, mode = "undirected") %>%
tidygraph::as_tbl_graph()
现在我们已经创建了这个对象，我们如何与它交互？
打印对象显示我们的网络有 37 个节点和 91 条边。我们还看到它是由两个 tibbles 定义的，一个定义节点数据，另一个定义边缘数据。
> leicester_network
# A tbl_graph: 37 nodes and 91 edges
#
# An undirected simple graph with 1 component
#
# Node Data: 37 x 1 (active)
name
<chr>
1 Beamont Park
2 Rushey Mead North
3 Stocking Farm & Mowmacre
4 Bradgate Heights & Beaumont Leys
5 Rushey Mead South
6 Belgrave North West
# ... with 31 more rows
#
# Edge Data: 91 x 2
from to
<int> <int>
1 1 3
2 1 4
3 2 3
# ... with 88 more rows
注意，节点数据只包含 MSOA 的名称，边数据根据节点在节点数据中的位置告诉我们哪些节点是连接的。例如，边缘数据的第一行告诉我们“Beamont Park”和“Stocking Farm & Mowmacre”是邻居。
如果我们想可视化网络，我们可以使用{ggraph}，它是{ggplot2}的网络特定扩展，工作方式大致相同。下面是我们如何创建一个简单的情节。节点位置由算法确定。
ggraph(leicester_network) +
geom_edge_link(colour = 'black', width = 2) +
geom_node_point(size = 5, colour = 'steelblue') +
theme_void() +
ggtitle("Network plot of MSOAs in Leicester",
subtitle = "Each node is an MSOA, and an edge joins two
nodes if they are neighbours")
(图片由作者创作)
如果我们能在之前的地图上覆盖这个网络，那就太好了。这将允许我们感觉检查网络，并使它更容易在网络和现实之间转换。
为了实现这一点，我们首先需要获得每个 MSOA 的中心坐标。幸运的是，我们下载的形状文件为每个 MSOA 提供了一对经纬度坐标。不幸的是，多边形本身处于不同的空间投影下。我们使用以下代码提取坐标并转换投影。
coords <- leicester_sf %>%
as.data.frame() %>%
select(LONG, LAT) %>%
sp::SpatialPoints(proj4string=CRS("+init=epsg:4326")) %>% # LAT LONG code
sp::spTransform(CRS("+init=epsg:27700")) %>% # UK grid code
as.data.frame() %>%
bind_cols(msoa11_hclnm = leicester_sf$msoa11_hclnm) # to allow us to bind on
coords看起来是这样的:
LONG LAT msoa11_hclnm
<dbl> <dbl> <chr>
1 457176\. 309094\. Beamont Park
2 461975\. 307767\. Rushey Mead North
3 458317\. 307856\. Stocking Farm & Mowmacre
4 456550\. 306747\. Bradgate Heights & Beaumont Leys
5 461046\. 307294\. Rushey Mead South
6 459732\. 307380\. Belgrave North West
7 460299\. 306733\. Belgrave North East
8 458059\. 305765\. Abbey Park
9 462606\. 306391\. Hamilton & Humberstone
10 459885\. 305741\. Belgrave South
# ... with 27 more rows
通过绑定 MSOA 名称，很容易将它加入到我们的网络对象中。我们本质上想要将left_join放到我们的节点表中。{tidygraph}让这变得非常简单。我们所要做的就是通过使用activate()告诉它我们想对节点表做些什么，然后我们就可以使用我们都知道并且喜欢的经典动词，比如left_join()。
# Join on cordinates and sf geometry
leicester_network <- leicester_network %>%
activate("nodes") %>%
left_join(coords, by = c("name" = "msoa11_hclnm"))
任何时候，您都可以使用as_tibble()提取激活的 tibble(即节点或边)。如果我们现在这样做，我们会看到坐标已经成功连接。
leicester_network %>%
as_tibble() %>%
head name LONG LAT
<chr> <dbl> <dbl>
1 Beamont Park 457176\. 309094.
2 Rushey Mead North 461975\. 307767.
3 Stocking Farm & Mowmacre 458317\. 307856.
4 Bradgate Heights & Beaumont Leys 456550\. 306747.
5 Rushey Mead South 461046\. 307294.
6 Belgrave North West 459732\. 307380.
现在我们创建一个图形，在这里我们使用 MSOA 形状文件作为一个geom_sf层，然后覆盖我们的网络。
ggraph(leicester_network, layout = "manual", x = LONG, y = LAT) +
geom_sf(data = leicester_sf, fill = "white") +
theme_void() +
geom_edge_link(colour = 'black', width = 2) +
geom_node_point(size = 5, colour = 'steelblue') +
ggtitle("MSOA network overlaying shapefile")
(图片由作者创作)
就这样，我们展示了网络正是我们想要的样子；每个节点都是一个 MSOA，如果两个节点是邻居，那么一条边连接两个节点。
下次
我们已经成功展示了为什么我们可以依靠网络来帮助我们理解 COVID 是如何传播的，以及如何构建和可视化这样一个网络。
在第 2 部分中，我们将看到网络结构如何让我们直观快速地创建强大的功能，然后使用这些新功能创建预测模型。
这篇文章中使用的所有代码都可以在 GitHub 上找到，这里的renv被用来处理包版本和依赖关系。
使用 R-Part 2 中的网络预测新冠肺炎的传播
生成网络要素并进行预测
在第一部分中，我们展示了为什么网络可以用于预测新冠肺炎的传播，我们创建并可视化了一个名为 MSOAs 的英国小区域地理网络。我们最后展示了莱斯特这些地区的网络，覆盖在莱斯特的地图上(如下所示)。我们关注莱斯特，因为这是英格兰第一个实行局部封锁的地区。
(图片由作者创作)
在第二部分也是最后一部分，我们将使用该网络来生成每个区域的要素。最后，我们将使用这些功能创建一个模型来预测一个地区在未来一周将获得的 Covid 病例数。
除了我们在第一部分中涉及的工作，我们将展示代码，你可以在 github 上找到所有的代码。
我们从两个 R 对象开始。第一个是msoa，它是一个包含 MSOA 案例数据的 tibble，第二个是msoa_network，它是英国每个 MSOA 的{tidygraph}网络对象，其中的节点是 msoa，如果它们是邻居，则一条边连接两个 msoa。上面的图像向我们展示了这个网络，但只是针对莱斯特。
这些都是用来创建本帖中显示的数据和可视化的包。
library(dplyr)
library(readr) # Loading the data
library(tidyr)
library(sf) # For the maps
library(ggplot2)
library(viridis)
library(cowplot)
library(glue)# Network packages
library(igraph) # build network
library(spdep) # builds network
library(tidygraph)
library(ggraph) # for plotting networks
使用我们的网络创建要素
有两种方法可以使用我们的网络来创建功能，这有助于将它们分开考虑。
第一类是基于网络结构的特征，例如，一个区域的邻居数量(也称为其度数)。这里的想法是，LAs 的分组和连接方式可能是 Covid 如何在网络中传播的良好指示。
第二个是基于邻近地区正在发生的情况的特征，例如，在邻近的 MSOAs 中有多少病例。这里的想法是，一个被高病例数包围的地区可以预期在接下来的一周内看到自己的上升。
基于网络结构的特征
我们将研究其中一些特性，展示我们是如何为我们的网络推导出这些特性的，然后展示这些特性在莱斯特每个 MSOAs 中的价值。
度数=邻居数量。你可以从下面的图片中看到，我们有两种不同的方式来观想这一点。
首先是使用网络，这样可以更容易地看到连接每个 MSOA 的边。使用这个网络，您可以很容易地计算任何中央 MSOAs 的所有邻居，以查看它是否与它们的度数相匹配。我们不能对接壤的 MSOAs 这样做，因为我们没有标绘他们所有的邻居，因为他们不都在莱斯特。
第二种方法是使用地图来可视化。这可以更容易地看到整个洛杉矶的价值分布。
(图片由作者创作)
**介数=一个节点有多少条最短路径。**这是一种节点重要性的度量。如果你想从一个 MSOA 旅行到另一个，总有至少一条最短的路线可以减少你必须经过的机场数量。因此，您可能会认为介数高的节点会有更多的人通过它。在现实中，这种方法不太可能代表这一点，而且现实生活中的最短旅行距离与您旅行经过的 MSOAs 最小化没有任何关系。然而，它确实提供了一些额外的背景，模型可能会发现有帮助。
同样，由于我们只是绘制了整个网络的一个子集，我们不能仅仅从莱斯特网络来计算介数。
(图片由作者创作)
三角形=该节点所属的三角形数量。三角形就是有三个节点相互连接的情况。在我们的例子中，当一个 MSOA 与两个也彼此相邻的 msoa 相邻时，就会发生这种情况。我们可能期望 Covid 在有许多三角形的区域中传播得更多，因为 MSOAs 连接得更多。
(图片由作者创作)
**传递性=有多少局部聚类。**这是对节点的邻居中有多少也是邻居的度量，通过将邻居对的数量除以该节点所在的三角形的数量来计算。社交网络具有很高的传递性，因为朋友的朋友很可能自己也是朋友。我们可能会认为更高的值会影响 Covid 在网络中的传播方式，因为这意味着区域之间的互联程度更高。
(图片由作者创作)
通过将方便的{igraph}函数与{tidygraph} mutate方法相结合，创建这些功能非常容易。
我们从网络对象开始。我们activate了节点，所以它知道下面的 mutate 将应用于底层节点数据。然后我们应用 mutate，将我们的新特性定义为应用于我们的网络对象的相应的{igraph}函数，它是使用.G()获得的。
msoa_features <- msoa_network %>%
activate("nodes") %>%
mutate(degree = igraph::degree(.G()),
betweeness = igraph::betweenness(.G()),
triangles = igraph::count_triangles(.G()),
transitivity = igraph::transitivity(.G(), type = "local")))
唯一的区别是对于传递性你必须指定你想要使用局部传递性(即每个节点的传递性)而不是网络传递性(这是整个网络的单一度量)。
基于邻域的特征
现在，我们希望根据相邻 MSOAs 中正在发生的情况(例如，病例数)来生成要素。与上面只需要生成一次的特性不同，我们需要为每周的案例数据生成这个特性。首先，我们将只看一个特定的星期，即第 25 周，这是莱斯特的第一次高峰。
首先，我们将案例数据连接到我们的网络对象。同样，我们利用了在{tidygraph}对象上使用像left_join()这样的{dplyr}动词的方式。
msoa_cases_network <- msoa_network %>%
activate("nodes") %>%
left_join(msoa_features %>% filter(week == 25),
by = c("name" = "msoa11_cd"))
现在有了一整套的{tidygraph}函数，可以让你在邻居级别应用自定义函数(例如，参见 map_local )然而，我们在这里想要的非常简单，最有效的计算方法(很远)是使用邻接矩阵。
你可能还记得上次提到的邻接矩阵定义了一个网络。对于网络中的每个节点，它都有一行和一列，如果一条边将节点 i 连接到节点 j ，则单元 (i，j) 的值为 1，否则为 0。
具有相应邻接矩阵的网络示例。Jesus Najera 创作的图片，作为图论博客系列的一部分。
我们首先使用邻接矩阵来创建我们的网络，然而，通过以下方式从现有网络中获得邻接矩阵非常容易:
adj_matrix <- igraph::get.adjacency(msoa_cases_network)
现在，如果我们有网络中每个节点的值的向量，例如 covid 情况，让我们想想如果我们将邻接矩阵乘以这个向量会发生什么。结果将是一个相同长度的向量。位置 i 将只是原始向量的所有值的总和乘以它们在邻接矩阵的行 i 中的位置。但是我们知道第 i 排除了第 i 排的邻居都是 0，而在那些地方，正好是 1 *。*因此，我们看到结果只是所有节点中所有 Covid 案例的总和，而这正是我们想要的！因为这只是矩阵乘法，所以它发生得非常快。
因此，我们现在知道，为了生成这种“相邻 MSOA 中的 covid 案例”特征，我们只需将邻接矩阵乘以每个 MSOA 中的 Covid 案例的向量。这可以通过以下方式轻松实现:
msoa_cases_network <- msoa_cases_network %>%
mutate(cases_in_neighbourhood =
as.vector(adj_matrix %*% .N()$cases),
cases_in_neighbourhood_average =
cases_in_neighbourhood/degree)
同样，我们可以将mutate应用于 tidygraph 对象。我们通过使用访问网络底层节点数据的.N()得到案例向量，然后使用$cases得到案例向量。我们使用矩阵乘法将我们的邻接矩阵与此相乘，在 R 中用%*%表示。最后，我们将结果包装在as.vector()中，这样 R，更确切地说是 mutate，就知道把它当作一个向量。我们通过除以节点的度来计算平均案例数。
为每周创建此内容
现在我们已经这样做了一周，我们只需要每周都这样做。我以下面的方式使用lapply来做这件事。如果你不熟悉应用系列，它类似于for循环，但更有效。我们每周应用的代码应该看起来非常类似于我们上面所做的。这里的关键区别是，我们在数据中每周都这样做，然后在最后将所有结果绑定在一起。
data_for_predictions <- lapply(msoa_features %>%
distinct(week) %>%
pull(week),
function(wk){
message(paste0("Creating features for week ", wk))
# Filter to the correct week
weekly_data <- msoa_features %>%
filter(week == wk)
# join onto the network and create network features
weekly_network <- msoa_network %>%
activate("nodes") %>%
left_join(weekly_data, by = c("name" = "msoa11_cd")) %>%
mutate(cases_in_neighbourhood =
as.vector(adj_matrix %*% .N()$cases),
cases_in_neighbourhood_average =
cases_in_neighbourhood/degree)
# Extract as tibble
weekly_network %>%
as_tibble()
}) %>% # Bind all these tibbles together
bind_rows()
这样，我们就创建了一个数据集，对于进行预测建模来说，它应该看起来更熟悉。我们每个 MSOA 每周有一行，有很多功能可以帮助我们预测下周。我们在下面提供了一个摘录，这样你就可以看到我们为某个 MSOA 的几个星期的子集创建的这些新功能。注意，我们已经创建了一个“cases_next_week”变量，因为这将是我们的结果变量。换句话说，我们将使用所有其他变量来尝试和预测下周会有多少病例。
(图片由作者创作)
像其他功能一样，我们可以在地图上看到任何一周的这些功能。下面的可视化显示了第 25 周的情况。换句话说，一个区域越亮，其相邻区域的病例越多。我们称这个变量为“案例 nb”。
(图片由作者创作)
创建预测病例的模型
我们将使用上述数据集，尝试使用其他变量来预测“cases_next_week”。这篇博文不打算深究预测建模或机器学习方法或最佳实践。相反，我们将向您展示如何使用上面创建的数据集创建一个简单的线性模型。如果你正在使用这个博客系列从网络中创建你自己的预测模型，你可以用你最喜欢的模型来代替。
训练和测试数据集
为了了解模型结果的概括程度，我们将数据分为训练数据集和测试数据集。我们这样做是为了让大约 70%的数据用于训练，30%用于测试。
我们使用 week 变量进行划分，因此我们对某个时间点之前的所有数据进行训练，并对其余数据进行测试。这是为了避免我们对过去的数据进行测试，这可能会导致数据泄漏。如果我们试图将这个模型付诸实践并预测下周发生的事情，我们将永远不会有来自未来的训练数据。
processed_df <- data_for_predictions %>%
select(name, week, degree:cases_in_neighbourhood_average) %>%
drop_na() # This is a bit of a cheat to get results quickly
training_threshold <- 0.7
training_weeks <- processed_df %>%
distinct(week) %>%
slice(1:round(nrow(.) * training_threshold)) %>%
pull()train_df <- processed_df %>%
filter(week %in% training_weeks)test_df <- processed_df %>%
filter(!week %in% training_weeks)
模型创建
现在我们使用lm函数创建一个模型。
linear_model <- lm(cases_next_week ~ .,
data = train_df %>% select(-name, -week))
在这里，我们可以检查模型的摘要，并放心地看到“cases _ in _ neighbourhood _ average”与您认为极其重要的“cases”具有相似的大小系数。这意味着 MSOA 多了一个病例，对下周的病例会产生类似的联想效应，而不是我们在每个邻近地区都多了一个病例。换句话说，模型依赖于我们的邻近特征。
(图片由作者创作)
有趣的是，看起来基于网络结构的特征与下周的案例有更不确定的关系，并且没有我们预期的效果。比如他们的系数大部分都是负的。这可能与以下事实有关:如果区域中没有太多 Covid，那么网络结构容易受到 Covid 传播的影响就没有多大意义。
我们可以用下面的方法为我们的测试数据创建预测。
test_pred <- stats::predict(linear_model,
test_df %>%
select(-name, -week, -cases_next_week))
这让我们能够比较对未知数据的预测和观察。例如，在这里，我们计算测试数据的 RMSE，看到它是 8.36，比训练数据达到的 0.35 高得多，显示模型遭受过拟合。
pred.obs <- tibble(pred = test_pred,
obs = test_df$cases_next_week)rmse <- pred.obs %>%
summarise(sqrt(mean((obs-pred)^2)))
使用模型
显然，像这样的模型有数百种不同的使用方式。一种方法是将输出可视化在地图上。下面显示了我们的模型根据第 40 周(在我们的测试数据中)的数据对第 41 周每个 MSOA 中的病例进行的预测，并将其与第 41 周实际发生的情况进行了比较。
左边的地图显示实际案例，右边的地图显示根据我们的模型预测的案例。(图片由作者创作)
模型性能和改进
现在你可能可以从上面的视觉化图像中看出，这个模型并不是特别好。它在很多领域都是错误的，即使它是正确的，它也严重依赖于该领域已经有多少案例，这不是很有帮助。为了让你了解这一点，如果我们用一个基线模型来代替这个模型，这个基线模型只是假设病例数不变，我们将得到 8.18 的 RMSE，比我们的线性模型好**。**
如果我们希望在现实世界中使用这个模型，那么需要做大量的工作来理解如何使用它，并围绕它设计一个模型管道。**这只是预测建模的本质。**如果你自己要做这件事，这里有一些你可能想探索的事情的想法:
- 将结果变量改为“病例数的变化”,以尝试消除对使用当前病例数的过度依赖。
- 使用分类结果变量，如“下周病例会增加 x%吗”或“病例会超过此固定阈值吗”。
- 考虑到每个 MSOA 的人口，所以较小的地区不会有偏见。
- 使用可以减少过度拟合的更复杂的模型(例如，random forest 或 xgboost)。
- 添加更多数据或特征，例如，考虑交通连接区域的方式或额外背景，如该区域的城市化程度。
- 更多的特征工程，例如标准化特征，或者向特征添加噪声(这可以通过减少过拟合来提高模型性能)。
结论
这篇文章有效地展示了你如何从关于事物如何联系的信息(如地理上的邻居)开始，并创建许多可用于预测建模的有用特征。
新冠肺炎继续对我们的生活和社会产生巨大影响。由于疾病传播的指数性质，得到任何关于哪里可能爆发疾病的警告会在应对的选择和影响方面产生巨大的差异。已经有很多令人印象深刻的建模工作在进行，试图给出这个额外的警告(例如帝国理工学院的这个工作)。这篇博客文章表明，网络可以用来提供额外的功能，使模型更加准确，从而帮助我们对抗这种疾病。
这篇文章中使用的所有代码都可以在 GitHub 上找到，这里的renv被用来处理包版本和依赖关系。
预测股市很难:创建一个机器学习模型(可能)不会有帮助
马库斯·温克勒在 Unsplash 上的照片
来自《走向数据科学》编辑的提示: 虽然我们允许独立作者根据我们的 规则和指导方针 发表文章，但我们并不认可每个作者的贡献。你不应该在没有寻求专业建议的情况下依赖一个作者的作品。详见我们的 读者术语 。
有许多关于制作机器学习模型来预测股票价格的文章和方法。我不是在这里宣布他们是错误的或被误导的。事实上，我认为这么多人对试图解决一个古老的问题感兴趣是件好事，对 Python 等编程语言的贡献在一定程度上帮助了投资的“民主化”，人们现在可以在几分钟内建立一个相对复杂的机器学习模型。相反，本文想要展示的是支持股票价格可能无法预测这一观点的证据。具体来说，本文将使用 Python 来强调股票价格或多或少遵循随机游走的思想。文章的主旨是，随机漫步并不排除击败市场的可能性，也不主张将投资模型扔到一边，但它应该发出信号，这些模型应该以极度的谨慎和勤奋来构建，并不断进行重新评估。建立交易算法的过程在很多情况下是一个终生的过程。
随机游走理论综述
股票价格的随机游走理论有很长的历史，最早是在 19 世纪中期提出的，然后在 20 世纪中期由尤金·法马和伯顿·马尔基尔等人推广，后者以他的经典之作《华尔街的随机漫步》使随机游走理论在投资界变得常见。该理论最简单的形式是，股票价格无法预测，因为股票价格的变化是随机的。该理论假定金融市场是有效的。也就是说，由于市场参与者是理性的利润最大化者，所有公开信息都包含在股票的当前价格中。如果存在任何异常，它将被迅速利用和消除，从而导致更高效的状态。例如，假设市场认为苹果的股票被低估了。随机游走理论假设市场参与者会立即购买股票，这反过来会导致股票不再被低估。这种“完美”的效率使得股票价格被相应地定价，并反映所有可用的信息。新消息是唯一能改变股票价格的东西，因为消息周期是不可预测的，所以股票价格会随机波动。
你可能已经注意到，随机漫步假说(RWH)做了几个假设，其中最主要的是金融市场是有效的。然而，泡沫存在的事实似乎推翻了这一假设，从而对随机漫步理论提出了质疑。此外，该理论催生了行为金融学领域，理查德·塞勒是其主要支持者之一，该理论旨在表明投资者在许多情况下远非理性行为者，这反过来意味着当行为者经常非理性行为时，市场很难有效。
尽管有这些(合理的)反对意见，随机漫步假说得到了大量经验证据的支持，本文的剩余部分将借助 Python 强调其中的一些。代码可以在这里找到
用 Python 演示随机漫步
根据 RWH 的说法，因为股票价格遵循随机游走，所以股票今天的价格是对其明天价格的最佳预测。为了检验这种说法，我们可以将某一特定股票在不同时间间隔内的滞后价格与该股票的最近价格进行比较，以确定它们是否代表今天的价格。
滞后 1、2、3、4 和 5 天的 DOCU
同样，这里的想法是评估今天的价格是否是明天价格的最佳指标。如果是这样的话，那么股票会相应地定价，因此会随机波动。另一方面，如果滞后价格与今天的价格几乎没有关系，那么今天的价格就不是股票价格的最佳预测者——也许过去的价格是——因此市场是无效的。使用线性回归来评估这一点，我们可以看到滞后价格与今天的价格有着极其密切的关系。在这个特殊的例子中，DocuSign 的滞后价格代表了今天的价格。
将这一分析扩展到纳斯达克 100 指数中的每一只股票，这一现象被证明是普遍的。纳斯达克 100 指数中每只股票的每个线性回归模型的 r 平方得分为 0.99，表明该模型几乎解释了数据中的所有变化。也就是说，滞后价格系列，尤其是 lagged_1，完全解释了最近的价格:
这是否意味着投资者不能选股
RWH 的一个含义是，传统的选股方法，如技术和基本面分析，用处不大。这两种技术都意味着投资者可以利用这些方法来制定有利可图的交易策略，但 RWH 认为这通常会弄巧成拙，因为交易者会利用并因此抵消这些异常，使市场变得有效。然而，该假说允许在某些情况下选股，分析师会比遵循简单的买入并持有政策的投资者做得更好，只要他能比其他分析师和投资者更快地识别实际价格和内在价值之间存在不可忽略的差异的情况，并且如果他能更好地预测重要事件的发生并评估它们对内在价值的影响。
结论
是的，随机漫步理论享有实证支持，并为算法设定了很高的标准，以持续不断地预测股价走势。但这并不自动意味着建立投资模型的雄心被抛弃。甚至 RWH 也允许 T2 的一些交易者表现优于其他交易者。吉姆·西蒙斯是一个克服困难的完美例子。他花了数年时间试图“解决市场问题”，一个方法接一个方法都失败了，最终取得了突破，经过了数年和无数个小时来完善他的方法。即使是现在，他的公司也只有一半的交易是正确的。重点不是阻止你建立下一个机器学习模型，而是现实地设定期望，并鼓励你长期坚持下去。
进一步阅读
股票市场价格的随机波动
用机器学习预测股市。标杆管理。
亨尼投资
用随机森林分类器绘制基线。
本文描述了一个机器学习模型的实现细节和基线性能，该模型用于预测过去 **8 年中 S & P500 公司下一年的收入。**如果你还没有看到介绍，请阅读用机器学习预测股市。引言。熟悉背景、问题陈述和我选择的解决方法。我们开始吧！
什么是基线？
基线模型是一个简单的模型，用于建立一个参考点，并有效地跟踪一个人的进展。我如何知道我在开发机器学习管道的各个方面时取得了进展？当新型号比以前型号表现得更好时，我知道我取得了进步。在本文中，我将构建这个先前的模型，它将成为下一个更复杂模型的基准。此外，它将帮助我对重要的特性建立一些直觉，并帮助确定进一步开发的优先级。
偏见意识
在机器学习中最危险的两个问题，在我看来就是未来泄漏和过拟合。
未来泄漏是指将表示在预测时不可用的信息的特征包括到训练数据集中。这是一个时间序列问题。因此，防止未来泄漏的第一步是确保在每个时间戳使用的功能不来自该时间戳之后发布的报告——这是一个快速验证，但我们最好确定。
过度拟合指的是以这样一种方式训练模型，它可以很好地预测过去(在训练集上做得很好)，但不能预测未来。在这项研究中，这个问题可能来自于财务特征的概念不会随着时间的推移而发生巨大的变化。如果我们在训练和验证集合中的不同时间戳都包括同一家公司，并且这家公司碰巧持续增加其收入，那么模型将最有可能了解什么特征组合代表这家公司(这将产生正确的预测)，而不是什么特征组合驱动实际增长。我希望我的模型能够对以前从未见过的公司做出正确的预测。因此，我需要确保它学习功能的趋势，而不是公司的趋势。
培训—验证—保持
有了 8 年的历史和 500 家公司，我将使用每个时间戳每个公司 (2)年的历史数据。由于每个时间戳我需要 2 年的历史数据，所以可预测的最早年份是 2014 (其中 2012 和 2013 构成特征)，最晚是 2019 (其中 2017 和 2018 构成特征)。我使用 2018 年和 2019 年作为样本外的维持，它不会用于训练模型。相反，我将使用它来评估和报告模型性能。
该模型在 2014 年至 2017 年使用群组 K 倍验证进行训练。使用这种方法，训练包括使用较晚的年份来预测较早的年份，这在偏差意识一节中被警告过。但是，在这种情况下，这不是一个问题，因为我使用 Group K-Fold 抽样方法，通过公司 ID 而不是随机地在培训和验证之间区分公司。如果某个公司出现在定型集中，它将根本不会出现在验证集中，无论是哪一年。同时，正如偏见意识一节所警告的，培训(A 公司，2016)和验证(B 公司，2014)不会引入过度拟合，因为 B 公司对 A 公司没有任何洞察力，除非这两家公司的目标高度相关。
选择模型
在这个阶段，我想要一个实现简单并且不需要太多调整的模型。逻辑回归可能是一个不错的选择，但它需要仔细考虑处理空值以及相关特征，以获得有意义的系数。 XGBoost 可能是解决这个问题的最佳模型，但是它需要相当大的努力来调整参数。我将改为实现一个随机森林。与 XGBoost 不同，它本身不能处理空值，但是它具有良好的开箱即用性能，并且它的特性重要性在某种程度上代表了它预测高概率的原因。
基准性能
评估分类模型的一个好方法是查看其提升曲线。提升曲线允许我回答以下问题:“如果我从模型中选择了前 N 个预测，我将在 P 个案例中的 K 个案例中是正确的”。
下面是第 5 级的提升曲线，预测公司在下一年的收入是否会翻倍。
该图描绘了三(3)个场景，向我们展示了我们的模型的性能，随机猜测将如何执行，以及完美模型将如何在维持集上执行。在这个类中，只有 ~12%的阳性，所以我们的模型必须相当有选择性。假设我们画出了概率的前 8%(垂直虚线)。最佳模型会将前 8%的概率分配给所有正样本，这将在我们的维持集中提取所有正样本的~ 75%(顶部虚线)。我们的模型提取了所有阳性(中间虚线)的 40% ，大约是完美模型的一半。随机猜测只有 8%的正面结果(底部虚线)，比我们的模型差五(5)倍。
通过提高概率等级，我的模型的性能更接近完美模型的性能。如果我画出预测的前 4%,我的模型的性能将几乎是完美的，换来的是给我更少的选择。虽然这听起来很诱人，但应该谨慎，因为数据量不够大，而且前 4%只产生了几十个样本。
也就是说，要成为一个成功的投资者，你可能只需要十几家公司。如果做得正确，这些人会排在前百分之几。
结论
我现在已经用一个开箱即用的随机森林模型对我在第 5 类上的表现进行了基准测试。虽然这个性能看起来不错，但我仍然不知道这个模型依靠什么来做出决策。为了找到答案，我将打开盒子，通过查看特征重要性，看看是什么驱动了这些预测。这将有助于我对推动下一年收入增长的因素建立一些直觉。下一篇见。
财务数据
在这项研究中，我使用了 www.vhinny.com 的历史金融数据。Vhinny 的 Alpha 数据集为 S & P500 公司提供从 2011 年开始的 8 年多的基础数据，如资产负债表、损益表和现金流量表。
我们连线吧！
我很高兴和与我志同道合的人联系，这是对财务独立的追求。如果你也在寻求经济独立，或者你想合作，交流想法，请随时联系我们！以下是我管理的一些资源:
- www.vhinny.com—投资研究平台，为您自己的分析提供金融数据
- https://www.linkedin.com/company/vhinny——加入我们 LinkedIn 上的社区，在那里我和其他贡献者分享与投资相关的内容
干杯！
来自《走向数据科学》编辑的注释: 虽然我们允许独立作者根据我们的 规则和指导方针 发表文章，但我们不认可每个作者的贡献。你不应该在没有寻求专业建议的情况下依赖一个作者的作品。详见我们的 读者术语 。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
