利用卫星图像预测赞比亚的粮食不安全状况
通过可视化我们的模型预测，我们能够识别作物和水可及性的高风险区域。
气候变化预计将通过更高的温度、更大的降雨量波动以及可能对作物和淡水供应产生毁灭性影响的极端天气，在 21 世纪增加营养不良的发生率。在义神星，我们专注于将机器学习应用于卫星图像，并开发了一个名为 EarthAI 的地理空间机器学习平台，以处理遥感图像的数量和复杂性。在我们的作物分类项目中，我们训练了一个模型，使用遥感图像对作物和水进行分类，并使用这些预测来研究赞比亚食物和淡水可获得性的变化，赞比亚在 2018 年全球饥饿指数中排名第五，属于“令人担忧”的类别。⁵ 60%的赞比亚人口生活在贫困中，40%的儿童因营养不良而发育不良。⁴通过这项研究，我们能够发现自然资源可获得性的变化，并跟踪作物可用性和轮作情况，从而了解该地区的粮食不安全状况。
“作为一家福利公司，将个人和组织与他们改善周围世界所需的信息联系起来对我们来说非常重要。几乎有无限多的方法可以做到这一点，但在无法获得信息或技术的项目中工作真的很有意义。让地球关注赞比亚的食品短缺对我们公司来说是一个鼓舞人心的挑战。”
~首席执行官布兰登·理查森
数据
由于它的全球覆盖范围和长期的免费和公开可用数据的存档，我们决定在这个项目中使用 Landsat 8 图像。这颗卫星大约每 16 天重访地球上的同一个地点，并在可见光、近红外和短波红外光谱中产生空间分辨率为 30 米的图像。
对于世界上的大部分地区，很难找到地面实况地理参考作物数据集，但幸运的是，自 2008 年以来，美国农业部(USDA)每年都会为美国大陆制作一个地面实况地理参考作物数据集。这个名为 CropScape 的数据集以 30 米的分辨率标注了 92 种不同的作物和 25 种土地覆盖类型。我们使用 CropScape 数据集在美国训练我们的模型，然后将这个模型推广到赞比亚。
图一。来源: Landsat 8 影像和作物景观农田数据层。
感兴趣区域
赞比亚享有长期的和平与稳定，但受到气候变化的严重影响。⁴我们特意选择了一个政治稳定的国家来控制混淆变量，如战争或经济制裁，以便粮食不安全问题是由气候变化引起的极端天气造成的，随着时间的推移，气候变化会慢慢影响作物，而不是由动荡的政治环境或其他超出遥感数据分析能力的无形变量造成的。
在赞比亚，农业支撑着 population,⁶85%的人口的生计，但赞比亚是一个内陆国家，大多数农民依靠降雨来灌溉作物，因此特别容易受到气候变化引起的天气波动的影响。⁴由于这些因素以及赞比亚种植的作物与美国的作物非常匹配，赞比亚是模型外推的一个简单选择。
利用赞比亚作物生产的一些基本统计数据(非洲数据、粮农组织国家报告和赞比亚数据门户)，我们根据 2015 年和 2016 年种植的公顷数和每种作物的产量确定了赞比亚最常见的作物。总之，我们在模型中包括了九种作物类型:玉米、小麦、小米、马铃薯、杂豆、花生、甘蔗、烟草和棉花。
光谱特征工程
用于训练我们的分类模型的特征是来自 Landsat 8 的原始波段(红色、绿色、蓝色、近红外、两个短波红外和一个沿海波段)和一组派生的水、植被和土壤指数。由于不同作物的植被指数在不同时间达到峰值，因此包含一个时间序列的特征值以区分不同作物非常重要。
归一化差异植被指数(NDVI)是用于测量植物健康的近红外和红色波段的比率。⁷较高的 NDVI 值表明植物更健康。⁷图 2 向您展示了 NDVI 时间序列如何区分作物。玉米的 NDVI 在 6 月和 7 月初达到高峰；棉花、花生、小米和烟草在 8 月达到高峰；干豆在九月达到高峰。
图 2:按作物和月份划分的平均 NDVI 时间序列图。
建模
义神星 EarthAI 平台建立在 Apache Spark 和亚马逊网络服务(AWS)基础设施之上。因为我们需要为这个建模工作处理的图像量几乎是 1tb，所以我们使用 EarthAI 以及时和分布式的方式来执行建模。我们的流程概述如下。
- 读入 Landsat 8 影像和 ground truth CropScape 数据，并通过重新投影和重新采样 ground truth 数据以匹配遥感数据，在空间上连接这两个栅格数据集。
- 从原始的 Landsat 8 波段导出水、植被和土壤指数。
- 将栅格数据分解为像素，因此数据集中的每一行都是一个像素，每一列都是不同的波段。
- 使用 Landsat 8 质量评估(QA)波段过滤掉不良数据，如被云覆盖的像素。
- 通过按月平均每个光谱特征来创建时间序列特征。
- 将数据分成训练集和测试集。
- 训练和评估分类模型。
- 在 2015 年和 2016 年对赞比亚的模型进行评分。
- 对模型的预测执行分类后更改检测。
图 3:使用 EarthAI 平台建模工作流。
因为我们选择的九种作物的作物日历并不同步，所以我们决定为夏季和冬季作物建立两个独立的模型。在美国，甘蔗、小麦和马铃薯是冬季作物，在夏末秋初种植，第二年夏天收获；而玉米、棉花、花生、小米、干豆和烟草是夏季作物，在早春种植，秋季收获。在赞比亚的⁸，由于半球不同，作物日历是相反的。⁹
对于夏季和冬季生长季节，我们表现最好的模型是梯度促进树木(GBT)模型。GBT 夏季模型的总体准确度为 0.876，F1 值为 0.792。类别度量如图 4 所示。与其他类相比，小米的召回率很低，因为我们的训练集不平衡，而且我们的此类示例太少。
GBT 冬季模型的总体准确度为 0.956，F1 值为 0.895。较少的班级和贫瘠的冬季景观使得这些班级比夏季班级更容易建模。类别度量如图 5 所示。
图 4:测试集的夏季模型度量。
图 5:测试集的冬季模型度量。
调查结果
我们对赞比亚两年的夏季和冬季模型进行了评分，并对赞比亚境内的区域进行了分类后变化检测。我们能够检测自然资源可用性的变化，如图 6 所示，2015 年至 2016 年间，塔米山淡水水库水位显著下降。
图 6:2015 年和 2016 年叠加在赞比亚姆库西区米塔丘陵大坝上的模型预测。
我们还能够跟踪作物的可用性和轮作时间，如图 7 所示的例子，显示了粮食作物向经济作物的转移。在需要的时候，经济作物会产生更多的收入，因为它们是为了销售而不是消费而种植的。
图 7:2015 年和 2016 年赞比亚农业区的模型预测。
通过可视化我们的模型预测，我们能够识别粮食不安全和水可及性的高风险区域。及时分享这些信息可以加快预防饥饿和营养不良的干预努力。
挑战
当试图将一组特定数据上训练的模型外推至一组新数据时，概化将是一个挑战。我们的模型是根据美国的影像进行训练的，因为我们在赞比亚找不到地面实况地理参考作物数据。美国幅员辽阔，并非所有州都表现出与赞比亚相似的特征，因此我们分析并选择了在气候、天气、土壤、生物群落、生态区域和作物可用性方面与赞比亚最相似的美国各州，以创建最具代表性的训练集。
这个项目的另一个挑战是需要处理的大量遥感数据以及我们完成特征工程的速度。幸运的是，义神星 EarthAI 的使用缓解了这些问题以及处理地理空间数据源的所有隐藏的复杂性。
“将地球观测(EO)数据纳入数据科学的主要挑战之一是处理不同的地图投影和坐标系统。另一个挑战是数据的巨大规模。我们从头开始设计义神星 EarthAI 来处理这两个挑战，对数据科学家几乎没有影响。EarthAI 的“栅格连接”功能就是一个例子，通过该功能，来自全球各地的大型数据集可以在空间和时间上与平台管理的复杂性相连接。”
~研发副总裁 Simeon Fitch
我们在这个项目中使用的遥感数据的大小约为 1tb，在 EarthAI 中聚合和过滤数据后，数据减少到 11GB。尺寸显著减小的部分原因是我们使用光学数据，并且必须应对云。如果一朵云在整个生长季节的任何时候覆盖了一个特定的像素，那么这个像素就会从我们的数据中删除，这对于我们的数据中已经很少见的作物(如豆类和小米)来说是不可取的。因此，我们开始研究减少缺失数据的不同方法，最有前途的方法之一是将合成孔径雷达(SAR)数据源整合到我们的模型中。
合成孔径雷达是一种雷达，它通过向地球发射微波脉冲来工作。当脉冲与物体接触时会被反射或散射，然后返回到卫星，在卫星上，传感器根据返回的回波生成图像。与需要阳光照射物体的被动光学传感器不同，SAR 可以在白天和晚上进行感知，从而产生更多的图像。SAR 的另一个优势是微波脉冲能够穿透云层，因此 SAR 即使在阴天也能产生图像。我们已经有了一些有希望的初步结果，将 SAR 和光学图像结合起来对作物进行分类，如图 8 所示。
图 8:有和没有 SAR 数据的模型预测的比较。
通过利用员工的数据科学专业知识和义神星 EarthAI 平台的可扩展处理能力，我们开发了一个模型，能够跟踪一个面临粮食不安全风险的国家的水和作物可用性。随着气候变化继续改变天气模式，受影响的国家将不得不适应不同的生长季节和水模式，因此，一个允许人们可视化作物和水分布随时间变化的模型将有助于各国适应并确定粮食不安全的高风险地区。
如有疑问或意见，请访问我们的网站:www . astrea . earth。
由 撰写考特尼 ，资深数据科学家
非常感谢义神星大学参与这个项目的其他数据科学家:
金伯利·斯科特博士，联合创始人&数据科学副总裁
杰森·布朗，高级数据科学家
埃里克·库伯特森，数据科学家
参考文献
- 【http://agris.fao.org/agris-search/search.do? recordID=GB2013200329
- https://en.wikipedia.org/wiki/Landsat_8
- https://www . NASS . USDA . gov/Research _ and _ Science/wh 田/sarsfaqs 2 . PHP # section 1 _ 1.0
- https://www . concern USA . org/story/worlds-ten-hungriest-countries/
- https://www . concern USA . org/WP-content/uploads/2018/10/GHI-2018 . pdf
- 【http://www.yieldgap.org/zambia
- https://medium . com/regen-network/remote-sensing-indexes-389153 e3d 947
- https://downloads . USDA . library . Cornell . edu/USDA-esmis/files/VM 40 xr 56k/dv 13 zw 65 p/w 9505297d/planning-10-29-2010 . pdf
- http://www.fao.org/giews/countrybrief/country.jsp?code=ZMB
用张量流预测 GCP 森林覆盖类型及模型部署
通过 Kaggle 竞赛开始使用 Tensorflow 并了解如何在 GCP 部署模型
迪皮卡·巴德利用张量流进行森林覆盖类型分类
在这篇文章中，我将分享:
- 我如何开始使用 Tensorflow
- 用深度学习解决一个卡格竞赛
- 在 GCP 部署模型
- 在 GCP 为曼梯·里建一条管道
本项目使用的数据来自森林覆盖类型的 Kaggle 竞赛。虽然这不是 Kaggle 上的一个活跃的竞争，但这正好符合我的数字/分类数据标准，可以很容易地进行预测，因此我们可以专注于在 Tensorflow 中建立模型，并为 GCP 建立一个小管道。数据的详细信息在其数据描述页面上提供。数据由美国地质调查局和 USFS(林务局)提供。本问题将预测七种森林覆盖类型:
1 —云杉/冷杉
2 —黑松
3 —黄松
4 —杨木/柳树
5 —白杨
6 —花旗松
7 —克鲁姆霍尔茨
我将通过以下步骤来解决这个问题:
- 正在加载数据集
- 预处理数据集
- Tensorflow 入门
- 创建张量流数据集
- 用 Keras 建模
- 训练模型
- 测试模型
- 向 Kaggle 提交结果
- 在 GCP 部署模型
1.正在加载数据集
从上述 Kaggle 竞赛中下载数据，并将其存储在您的 google drive 或本地的适当文件夹中。将文件夹路径设置为FOLDER_NAME变量。在我的例子中，我将数据存储在 google drive 中，并使用 Google Colab 读取输入文件。如果你想在 Kaggle 笔记本上写代码，你可以按照我在 Kaggle和本文一起发表的代码来写。这些变化仅发生在 Kaggle 笔记本电脑加载和存储数据的情况下。
让我们从安装驱动器开始使用 Google Colab:
这将给出一个链接，以获取代码，你需要输入到输入框出现。一旦完成，我们就可以在 dataframe 中加载数据了。
输出:
部分输出
训练数据集有 15120 行。从描述表中可以看出，土壤类型 7 和 15 具有恒定的 0 值，因此应将其移除。除了分类列(Soil_Types和Wilderness_Areas)之外，一些具有数值的列应该被规范化以获得更好的结果。在下一步中，我们将执行所有的预处理步骤，以便为预测准备好数据。
2.预处理数据集
所有土壤类型的类别数据列将从一个热编码表合并为一列，类似的荒野区域列也将被转换。
输出:
土壤类型 8 和 25 每种类型只有一行，因此它们被转换成另一列，标有 NA 下标。这是可选的，您也可以删除这些列。
对于数字列，MinMaxScaler是将被应用来获得标准化列的转换器。在此之前，我们需要将数据分为训练、测试和验证，以便对这些数据进行标准化。
将数据分为训练、val 和测试
数据将分为训练、验证和测试。培训占 60%，测试占 20%，验证占 20%。
输出:
一旦数据被分割，归一化可以应用如下:
3.Tensorflow 入门
Tensorflow 是谷歌基于 Theano (Python 库)的开源深度学习库，用于研究和生产。你需要知道的一些核心组件是张量和图形。张量是一个 n 维向量或矩阵，用于存储数据和执行功能。Graph 是描述节点之间所有操作和连接的地方，可以在多个 CPU 或 GPU 上运行。
在这篇文章中，我将解释我如何使用 Tensorflow 创建我的第一个模型，而不是 tensorflow 的基础知识。如果你有兴趣学习，可以从这里过基础教程。我必须承认浏览 Tensorflow 文档不像 PyTorch 那么容易。使用 PyTorch，我能够在一周内构建我的第一个模型，但使用 Tensorflow 文档，很难理解加载数据本身的正确方式，此外，函数格式之间存在冲突，并且它与其他版本的向后兼容性也存在问题。因此，接下来的教程将会导致其他不容易解决的错误。
在你的机器或谷歌 Colab 中安装 tensorflow。我用过 Google Colab 搭配 GPU 运行时。这个问题用的是 Tensorflow 2。
!pip install tensorflow
4。创建张量流数据集
在所有的处理之后，我们将把数据放入 tensorflow 数据集。这将有助于构建模型的输入管道。对于数据集，必须定义 batch_size，这是可以成批而不是一次访问数据的大小。你也可以随意调整行。
输出:
输出的一部分(水平方向的列名无法放入一个屏幕截图中)
5.用 Keras 建模
将使用具有前馈神经架构的 tensorflow Keras 构建模型。使用 tensorflow Keras layers.DenseFeatures构建要素图层。分类列和数字列被分别处理，以创建代码中所示的输入层。模型在build_model 函数中定义，具有 100 和 50 两个隐藏层，后面是包含输出神经元数量为 8 的输出层。它是 8，因为类型是从 1 到 8 的整数值，而不是 0 到 7。在混淆矩阵中，我们忽略了类 0。使用优化器adam和激活功能relu。建筑模型和管线的代码如下:
6.训练模型
对于训练，我们为模型的拟合函数提供训练和验证数据集。验证损失使得在训练期间更容易发现模型的过度拟合。模型被训练 100 个纪元。
输出:
部分输出
在输出中，很明显，网络是基于随着验证精度增加的训练精度来学习的。如果在某个点上，验证准确性在增加，而训练准确性在增加，那么这就是模型过度拟合的点。这样你就能知道在那个时候停止纪元。如果它是随机的，那么这个模型没有学到任何东西。有了设计好的架构，我在训练集和验证集上的准确率分别达到了 85%和 84%。
在模型摘要中，您可以看到每一层的参数数量。在第一个隐藏层中，51 个输入特征连接到 100 个隐藏节点，对于完全连接的网络具有 5100 个权重，对于每个节点具有 100 个偏差参数，总计 5200 个参数。在下一层中，100 个节点连接到 50 个节点，具有 5000 个连接，第二隐藏层中的每个节点具有 50 个偏置参数，构成 5050 个参数。类似地，计算下一批参数。这就是你阅读模型摘要的方式。它显示最后学习了多少个参数。最后模型被保存到一个目录中。
7.测试模型
为了测试模型如何处理测试数据，我们将对 sklearn 的报告进行分类。分类报告将显示每种森林覆盖类型的精确度、召回率和 f1 分数以及平均精确度等。
输出:
该模型的平均准确率达到了 80%左右。对于前馈神经网络的简单结构，这是足够好的结果。现在，我们可以使用模型获得 kaggle 测试数据集的结果，并将其提交给 kaggle。
8.向 Kaggle 提交结果
我们需要预测 kaggle 上给出的测试数据的封面类型。该数据与第一步中下载的数据位于同一文件夹中。预期的输出是包含Id和Cover_Type列的 csv 文件。对训练数据进行的类似转换必须对测试数据进行。以下代码显示了如何获取测试数据的结果:
输出:
一旦你准备好了文件，你可以把它上传到 Kaggle 竞赛页面的我的提交部分。一旦提交，你可以在一段时间后看到分数。提交后，我在测试数据上得到了~62%的准确率。你可以用不同的实验做尽可能多的提交，努力提高这个分数。就这些，你可以开始参加不同的比赛，用不同种类的数据集做实验。我从有数字/分类数据的简单预测问题开始。我的目标是通过一个真实世界的例子来学习张量流，所以我从 Kaggle 竞争数据集开始，该数据集并不活跃，但有简单的问题需要处理。
9.在 GCP 部署模型
谷歌人工智能平台
前一节中保存的模型可用于部署在 google cloud 中，以便您拥有的任何应用程序都可以访问它。我认为你对谷歌云有基本的了解，并且在这方面有所建树。正如我不会解释的那样，如何开始学习 GCP(你可以在 coursera / GCP Qwiklabs 上学习很多入门课程)。
先决条件:
你可以创建一个免费的谷歌云账户，如果你还没有的话，它有 300 美元的免费信用。在接下来的步骤中，我们需要这一点。同时在你的电脑上安装 GCP 软件开发工具包。如果没有项目，在 Google Cloud 中创建一个项目。确保在 IAM 中下载一个服务帐户密钥，并存储在环境变量中。(关于通过命令控制台与 GCP 交互的基本设置，请参考谷歌云文档)
export GOOGLE_APPLICATION_CREDENTIALS=<PATH TO KEY JSON FILE>
运行以下命令，使用 Google Cloud 帐户进行身份验证，并根据输出中的给定说明设置项目。
gcloud auth login
第一步，将前面章节中的文件夹forest_model_layer_100_50_epoch_100上传到谷歌云存储中。我创建了一个 bucket forest-cover-model并在那个位置上传了文件夹。
一旦完成，您就可以部署模型了。
gcloud components install beta
这在部署模型之前是必需的。
gcloud beta ai-platform versions create v7 \--model ForestPredictor \--origin gs://forest-cover-model/forest_model_layer_100_50_epoch_100/ \--runtime-version=2.1 \--python-version=3.7
v7 版本是在我做了一些实验后有效的版本。可以以 v1 作为版本名。这样，您可以保留模型的不同版本。合适的运行时版本可以在这里找到。Origin 参数是存储模型的 google 存储路径的路径。
最初我计划用预处理类等定制预测例程。但不幸的是，在创建了所有这些之后，在部署的时候，我才知道它只适用于
tensorflow>=1.13,<2。它正在发展，所以它可能会在未来支持，检查这里的更新。(自定义管道的代码存在于我的库中，我已经分享给感兴趣的人)
要测试部署的模型，可以浏览到 AI Platform > Models，在模型名称下点击想要测试的版本号。有测试和使用选项，你可以给自定义输入。输入格式如下，使用下面的例子进行测试。
{
"instances":[
{
"Elevation":0.4107053527,
"Aspect":0.9833333333,
"Slope":0.2121212121,
"Horizontal_Distance_To_Hydrology":0.0,
"Vertical_Distance_To_Hydrology":0.2235142119,
"Horizontal_Distance_To_Roadways":0.3771251932,
"Hillshade_9am":0.7716535433,
"Hillshade_Noon":0.842519685,
"Hillshade_3pm":0.6141732283,
"Horizontal_Distance_To_Fire_Points":0.9263906315,
"Wilderness_Areas":"Wilderness_Area1",
"Soil_Types":"Soil_Type33"
}
]
}
您可以看到如下输出:
输出给出了输入的每种覆盖类型的概率。
一旦这些都工作正常，您就可以使用该模型对任何输入数据进行预测。我已经在 Github 资源库中分享了代码。文件forest_classification.py包含调用模型和输入管道的代码。
我将输入数据分成小块，因为它不能一次返回所有行的结果。这些是你需要在应用程序中管理的东西。如果您有兴趣探索如何为部署在 GCP 的模型创建管道，您可以参考该文档。
Tadaa！您已经准备好在实际解决方案中使用该模型，并为各种实验维护版本控制，并在此基础上跟踪性能。在现实世界中，除了只是建立一个模型，这些技能是非常重要的。有趣的是，如果您愿意的话，您可以尝试看看这些管道在其他云环境中是如何构建的。如果您希望尝试如何使用各种参数或其他优化函数来优化模型，那么您可以参考我以前的文章，在那里我提出了在深度神经网络中构建模型的改进/其他机会。尽管这是针对 PyTorch 的，但是无论您使用哪种框架，更改架构和参数的基础都是相同的。你必须找到在其他库中定义和做同样事情的方法。希望这能帮助你开始学习 GCP 和 Kaggle 以及 Tensorflow 框架！
一如既往—愉快的实验和学习:)
用机器学习工作流预测森林覆盖类型
机器学习
又名神奇的树:在哪里找到，如何检测它们🌲
在本文中，我将解释如何使用端到端的工作流来处理多类分类监督机器学习(ML)问题(或项目)。
- **监督:**特征(生成预测的变量)和目标(待确定的变量)在数据集中可用。
- **多类分类:**有七个离散的类别来区分目标。
该项目基于机器学习社区中一个著名的数据集，称为森林覆盖类型，可在 UCI 机器学习知识库中下载。
来自原始数据集的分层样本用于应用工作流和独立测试集以生成最终预测，被用作 Kaggle 中初学者友好型竞赛的一部分。
**项目(和竞赛)的目标:**以最佳精度预测北科罗拉多州罗斯福国家森林四个不同荒野地区的七种不同覆盖类型。
四个荒野地区是:
- 1:拉瓦
- 2:纽塔
- 3:科曼奇峰
- 4: Cache la Poudre
Cache la Poudre 荒野地区(来源)
在Cover_Type栏中从 1 到 7 的七个类别，待分类:
- 1:云杉/冷杉
- 2:黑松
- 3:美国黄松
- 4:杨木/柳木
- 5:阿斯彭
- 6:道格拉斯冷杉
- 7:克鲁姆霍尔茨
目标(封面类型)的名字让我想起了神奇的野兽，所以我称它们为神奇的树，以增加项目的想象力🙃。
**端到端机器学习工作流程步骤:**为了对封面类型进行分类并回答发起的问题，在哪里找到奇异的树以及如何检测它们，将遵循以下步骤:
- 理解、清理和格式化数据
- 探索性数据分析
- 特征工程和选择
- 比较几种机器学习模型
- 对最佳模型执行超参数调整
- 解释模型结果
- 用测试数据评估最佳模型(回答初始问题)
- 总结和结论
我将在本文中提供该项目的亮点，项目背后的全面分析和完整代码可在 Kaggle 笔记本和 GitHub 中获得。
对于数据角力和可视化 numpy ， pandas ， matplotlib ，seaborn；为了建立机器学习模型 xgboost ，lright GBM和scikit-learn；为了执行预处理步骤，将使用 scikit-learn 。
1.理解、清理和格式化数据
让我们加载训练数据并创建trees数据框:
trees = pd.read_csv("/kaggle/input/learn-together/train.csv")
我总是发现查看第一行和最后一行很有用:
数据集的第一列和第一行
数据集的最后一列和最后一行
第一个观察结果是，Vertical_Distance_To_Hydrology列中有一些负值。我将在检查异常&异常值部分对此进行更详细的检查。
为了理解trees数据框架，让我们看看数据类型和描述性统计。使用 pandas info方法，我们可以列出非空值和数据类型:
**Data columns (total 56 columns):
Id 15120 non-null int64
Elevation 15120 non-null int64
Aspect 15120 non-null int64
Slope 15120 non-null int64
Horizontal_Distance_To_Hydrology 15120 non-null int64
Vertical_Distance_To_Hydrology 15120 non-null int64
Horizontal_Distance_To_Roadways 15120 non-null int64
Hillshade_9am 15120 non-null int64
Hillshade_Noon 15120 non-null int64
Hillshade_3pm 15120 non-null int64
Horizontal_Distance_To_Fire_Points 15120 non-null int64
Wilderness_Area1 15120 non-null int64
Wilderness_Area2 15120 non-null int64
Wilderness_Area3 15120 non-null int64
Wilderness_Area4 15120 non-null int64
Soil_Type1 15120 non-null int64
Soil_Type2 15120 non-null int64
Soil_Type3 15120 non-null int64
Soil_Type4 15120 non-null int64
Soil_Type5 15120 non-null int64
Soil_Type6 15120 non-null int64
Soil_Type7 15120 non-null int64
Soil_Type8 15120 non-null int64
Soil_Type9 15120 non-null int64
Soil_Type10 15120 non-null int64
Soil_Type11 15120 non-null int64
Soil_Type12 15120 non-null int64
Soil_Type13 15120 non-null int64
Soil_Type14 15120 non-null int64
Soil_Type15 15120 non-null int64
Soil_Type16 15120 non-null int64
Soil_Type17 15120 non-null int64
Soil_Type18 15120 non-null int64
Soil_Type19 15120 non-null int64
Soil_Type20 15120 non-null int64
Soil_Type21 15120 non-null int64
Soil_Type22 15120 non-null int64
Soil_Type23 15120 non-null int64
Soil_Type24 15120 non-null int64
Soil_Type25 15120 non-null int64
Soil_Type26 15120 non-null int64
Soil_Type27 15120 non-null int64
Soil_Type28 15120 non-null int64
Soil_Type29 15120 non-null int64
Soil_Type30 15120 non-null int64
Soil_Type31 15120 non-null int64
Soil_Type32 15120 non-null int64
Soil_Type33 15120 non-null int64
Soil_Type34 15120 non-null int64
Soil_Type35 15120 non-null int64
Soil_Type36 15120 non-null int64
Soil_Type37 15120 non-null int64
Soil_Type38 15120 non-null int64
Soil_Type39 15120 non-null int64
Soil_Type40 15120 non-null int64
Cover_Type 15120 non-null int64
dtypes: int64(56)
memory usage: 6.5 MB**
用describe方法，我们可以观察到描述性统计:
第一列的描述性统计
info法提供了一些有价值的信息:
检查异常&异常值:
观察到的第一个异常是Vertical_Distance_To_Hydrology列中的负值。定义是:
到最近地表水特征的垂直距离
通过一些研究和使用逻辑，负值表明最近的地表水低于该数据点或低于海平面。两种情况都有道理，所以我要保留负值。
为了帮助未来的 ML 模型掌握数据中的模式，我将搜索异常值并使用极端异常值方法来确定它们。
如果数据点位于第一个四分位数以下或第三个四分位数以上的四分位数间距的 3 倍以上，则数据点将被删除。
*# loop through all columns to see if there are any outliers*
for column in trees.columns:
if outlier_function(trees, column)[2] > 0:
print("There are {} outliers in {}".format(outlier_function(trees, column)[2], column))
知道了荒野面积和土壤类型列是一次性编码的，我们就可以专注于其余部分:
**There are 53 outliers in Horizontal_Distance_To_Hydrology
There are 49 outliers in Vertical_Distance_To_Hydrology
There are 3 outliers in Horizontal_Distance_To_Roadways
There are 7 outliers in Hillshade_9am
There are 20 outliers in Hillshade_Noon
There are 132 outliers in Horizontal_Distance_To_Fire_Points**
山体阴影列是特定时间阴影的 RGB 颜色表示，因此范围已经固定在 0 到 255 之间。
考虑到Horizontal_Distance_To_Firepoints具有最大数量的异常值和最宽的数据范围【0，6993】，我将只删除该列中的异常值。
trees = trees[(trees['Horizontal_Distance_To_Fire_Points'] > outlier_function(trees, 'Horizontal_Distance_To_Fire_Points')[0]) &
(trees['Horizontal_Distance_To_Fire_Points'] < outlier_function(trees, 'Horizontal_Distance_To_Fire_Points')[1])]
2.探索性数据分析
EDA 是该工作流中的第一步，在该步骤中，针对特征选择启动决策过程。通过观察目标的分布、特征与目标的关系以及特征之间的联系，可以获得一些有价值的见解。
我的偏好是从查看目标开始，然后检查特征及其与目标的关系。
目标分布:
奇异的树/标签/覆盖类型的分布
数据集具有平衡的标签，导致每个类别的封面类型数量几乎相等。当我们应用 ML 算法时，这将是一个优势，因为模型将有很好的机会学习所有类的模式，而不需要进一步的平衡策略。
荒野地区——覆盖类型:
为了了解这种关系，荒野区域列将被反向热编码。
trees['Wilderness_Area_Type'] = (trees.iloc[:, 11:15] == 1).idxmax(1)
不同荒野地区覆盖类型的 KDE 样地
荒野区域是确定Cover_Type的显著特征。
荒野地区—土壤类型—覆盖类型:
反向一热编码将使用以下函数应用于土壤类型列。因此，将添加一个具有 1 到 40 之间的离散数字的列Soil_Type。
不同荒野地区土壤类型和覆盖类型散点图
不同的荒野地区由一些特定的树木组成。有趣的是，有一种奇妙的树，棉白杨/柳树，特别喜欢生长在 Cache la Poudre(荒野区域 4)。虽然云杉/冷杉、黑松、白杨和花旗松可以在任何土壤类型中生长，但其他覆盖类型可以在特定的土壤类型中生长。
连续特征之间的关系:
土壤类型和荒野区域列是离散的，并且都是一个分类要素的一个热编码版本。
其余要素被视为连续要素:高程、坡向、坡度、水平距离水文、垂直距离水文、水平距离道路、山体阴影上午 9 点、山体阴影中午 3 点、水平距离火灾点。
为了用一个函数将它们可视化，将绘制 Seaborn 的PairGrid,以提供用不同的图调整上下对角线的灵活性:
连续特征的成对网格
上半部分显示了具有皮尔逊系数的 KDE 图，下半部分显示了散点图。对角线是特定特征的直方图。正如所料，山体阴影特征是共线的:
- 山体阴影中午-山体阴影下午 3 点
- 山体阴影下午 3 点—山体阴影上午 9 点
这些对为模型提供相同的输入，为了更好的可解释性，其中一个将在特征工程&选择中被删除。
特性和目标的皮尔逊系数:
作为 EDA 的最后一步，当观察特征和目标的 Pearson 系数时，只有 1%的 one-hot-encoded 土壤类型列在确定Cover_Type时是有效的(本文未显示，但在笔记本这里)。因此，它们将被排除，皮尔逊系数将被重新使用连续特征、一个热编码的荒野区域、Soil_Type和Cover_Type。
continuous_variables = trees.columns[1:11].tolist()wilderness_areas = sorted(trees['Wilderness_Area_Type'].value_counts().index.tolist())all_features_w_label = continuous_variables + wilderness_areas + ["Soil_Type"] + ["Cover_Type"]trees_w_numeric_soil = trees[all_features_w_label]
cover_type 数据集的皮尔逊系数热图
深色代表强烈的相关性。不幸的是，最后一列由浅色组成，导致特征相对于目标的弱皮尔逊系数，在[-0.22，0.12]的范围内。
用标签编码的 Soil_Type列与Cover_Type有更强的相关性。
Hillshade_9am在确定Cover_Type时最不重要。因此，它将在下一节中被删除。
3.特征工程和选择
在已有特征的基础上提取新特征，用一些方法和算法消除特征，称为特征工程&选择。
水文要素有水平距离和垂直距离，两者的欧几里德距离相加闪烁。
trees_w_numeric_soil['Euclidian_Distance_To_Hydrology'] = (trees_w_numeric_soil['Horizontal_Distance_To_Hydrology']**2 + trees_w_numeric_soil['Vertical_Distance_To_Hydrology']**2)**0.5
此外，添加数字特征的线性组合是特征工程中的常见做法。对于一些数字特征，两个变量的平均值相加:
trees_w_numeric_soil['Mean_Elevation_Vertical_Distance_Hydrology'] = (trees_w_numeric_soil['Elevation'] + trees_w_numeric_soil['Vertical_Distance_To_Hydrology'])/2trees_w_numeric_soil['Mean_Distance_Hydrology_Firepoints'] = (trees_w_numeric_soil['Horizontal_Distance_To_Hydrology'] + trees_w_numeric_soil['Horizontal_Distance_To_Fire_Points'])/2trees_w_numeric_soil['Mean_Distance_Hydrology_Roadways'] = (trees_w_numeric_soil['Horizontal_Distance_To_Hydrology'] + trees_w_numeric_soil['Horizontal_Distance_To_Roadways'])/2trees_w_numeric_soil['Mean_Distance_Firepoints_Roadways'] = (trees_w_numeric_soil['Horizontal_Distance_To_Fire_Points'] + trees_w_numeric_soil['Horizontal_Distance_To_Roadways'])/2
另一种常见的做法是对数字特征执行对数和平方根变换。再添加 5 个特征后，对正特征应用平方根变换:
for col **in** trees_w_numeric_soil.columns:
if trees_w_numeric_soil[col].min() >= 0:
if col == 'Cover_Type':
next
else:
trees_w_numeric_soil['sqrt' + col] = np.sqrt(trees_w_numeric_soil[col])
在重新访问皮尔逊系数之后，如果新添加的特征显示出与目标更强的相关性，则保留该新添加的特征，而丢弃原始特征。此外，Hillshade_9am因与Hillshade_3pm密切相关而被剔除。
特征工程后的皮尔逊系数
最后的特点是:
*# final list of features*
transformed_features = ['sqrtHorizontal_Distance_To_Hydrology', 'sqrtMean_Distance_Hydrology_Roadways', 'sqrtEuclidian_Distance_To_Hydrology',
'Mean_Elevation_Vertical_Distance_Hydrology', 'Mean_Distance_Firepoints_Roadways', 'Mean_Distance_Hydrology_Firepoints']all_features = (['Elevation', 'Aspect', 'Slope', 'Vertical_Distance_To_Hydrology', 'Horizontal_Distance_To_Roadways',
'Hillshade_Noon', 'Hillshade_3pm', 'Horizontal_Distance_To_Fire_Points' ] +
wilderness_areas + ['Soil_Type'] + transformed_features)trees_training = trees_w_numeric_soil[all_features]
labels_training = trees_w_numeric_soil["Cover_Type"].as_matrix()
为了比较 ML 模型并建立基线，trees_training和labels_training数据帧被分成训练集和验证集。
X_train, X_valid, y_train, y_valid = train_test_split(trees_training, labels_training, test_size=0.2, random_state=1)
- 训练集用作输入，以便机器学习模型可以捕捉特征中的模式，并利用它们来区分目标。
- 验证集用于评估 ML 模型的性能，并量化其将模式归纳到新数据集的能力。
创建基线度量:
在深入研究 ML 分类算法之前，我将计算一个常识基线。常识基线可以定义为了解该领域的人如何在不使用任何 ML 技巧的情况下解决问题。它可以用人类的直觉以及一个虚拟或简单的算法来计算，只需要几行代码。
我将使用 scikit-learn 库中的一个虚拟算法。使用该算法，我将建立一个具有准确性的基线度量，即所有覆盖类型中正确预测的覆盖类型的百分比。准确性是本次竞赛的评估标准，将在整个项目中使用，记住它对于某些分类问题来说不是最有效的标准。
基线指标在某种程度上很重要，如果机器学习模型不能击败人的简单直观预测或算法的猜测，那么原始问题需要重新考虑或训练数据需要重新构建。
*# Create dummy classifer*
dummy = DummyClassifier(strategy='stratified', random_state=1)*# train the model*
dummy.fit(X_train, y_train)*# Get accuracy score*
baseline_accuracy = dummy.score(X_valid, y_valid)
print("Our dummy algorithm classified {:0.2f} of the of the trees correctly".format(baseline_accuracy))
**Our dummy algorithm classified 0.14 of the of the trees correctly**
现在，我期望以下 ML 模型击败 0.14 的准确度分数！
4.比较几种机器学习模型
有时很难知道哪种机器学习模型将有效地解决给定的问题。所以，我总是尝试几种机器学习模型。
一项研究显示，基于树和基于距离的算法在分析的 165 个数据集上优于其他最大似然算法。
机器学习算法性能比较(来源)
我将比较一个基于距离的算法和四个基于树的算法的准确性。
基于距离:
基于树:
- llight GBM库中的 light Gradient Boosting Machine(light GBM)分类器
- 来自 xgboost 库的额外梯度增强(XGBoost)分类器
- 来自 scikit-learn 库的随机森林分类器
- 来自 scikit-learn 库的额外树(随机森林)分类器
(型号后面的代码是这里是。)
与研究结果不同，额外树分类器优于其他分类器。所有模型都优于基线度量，表明机器学习适用于森林覆盖类型的分类。
我们最好的模型，额外树分类器，是一个基于集合树的模型。sci kit-learn 库中的定义如下:
该类实现了一个元估计器，它在数据集的各个子样本上拟合多个随机决策树(也称为额外树),并使用平均来提高预测精度和控制过度拟合。
与随机森林算法的主要区别在于:
- 不是为节点中的特征寻找最有区别的分裂值，而是完全随机地选择分裂作为阈值。
- 子样本取自整个训练集，而不是训练集的引导样本。
因此，该算法对于过拟合变得更加鲁棒。
5.对最佳模型执行超参数调整
寻找模型参数的最佳组合被称为超参数调整，它可以极大地提高模型的性能。我将使用交叉验证的随机搜索算法进行超参数调整:
- **随机搜索:**在一个范围内定义一组 ML 模型的参数，并输入到 sklearn 的T0。该算法随机选择一些参数组合，并将定义的
score(精确度，对于这个问题)与迭代进行比较。可以用参数n_iter控制随机搜索运行时间和迭代次数。 - **K-Fold 交叉验证:**一种用于评估超参数在整个数据集上的性能的方法。(目测此处为)。不是将数据集分成训练集和验证集的两个静态子集，而是对于给定的 K，数据集被均等地划分，并且随着迭代，不同的 K-1 个子集被训练，并且用不同的子集测试模型。
使用以下一组参数和 5 重交叉验证，RandomizedSearchCV将寻找最佳组合:
*# The number of trees in the forest algorithm, default value is 100.*
n_estimators = [50, 100, 300, 500, 1000]
*# The minimum number of samples required to split an internal node, default value is 2.*
min_samples_split = [2, 3, 5, 7, 9]
*# The minimum number of samples required to be at a leaf node, default value is 1.*
min_samples_leaf = [1, 2, 4, 6, 8]
*# The number of features to consider when looking for the best split, default value is auto.*
max_features = ['auto', 'sqrt', 'log2', None]
*# Define the grid of hyperparameters to search*
hyperparameter_grid = {'n_estimators': n_estimators,
'min_samples_leaf': min_samples_leaf,
'min_samples_split': min_samples_split,
'max_features': max_features}*# create model*
best_model = ExtraTreesClassifier(random_state=42)
*# create Randomized search object*
random_cv = RandomizedSearchCV(estimator=best_model, param_distributions=hyperparameter_grid,
cv=5, n_iter=20,
scoring = 'accuracy',
n_jobs = -1, verbose = 1,
return_train_score = True,
random_state=42)*# Fit on the all training data using random search object*
random_cv.fit(trees_training, labels_training)
random_cv.best_estimator_
以下是参数的最佳组合:
n_estimators= 300max_features=无min_samples_leaf= 1min_samples_split= 2
当我将它们输入额外的树分类器时:
**Accuracy score in the previous extra random forests model: 0.8659106070713809
Accuracy score after hyperparameter tuning: 0.885923949299533**
导致精确度增加 2 个点。
另一种搜索方法是GridSearchCV，与RandomizedSearchCV相反，搜索是在给定参数的每一个组合上进行的。(此处应用和讨论【原笔记本)
6.解释模型结果
混乱矩阵:
可视化分类模型结果的最常见方法之一是混淆矩阵。
奇异树混淆矩阵将是一个 7x7 矩阵。我将使用标准化的混淆矩阵，因此在该特定类别的所有猜测中，正确猜到的实际封面类型的百分比将出现在矩阵的对角线上。
非对角线元素将显示模型错误标记的元素。混淆矩阵对角线上的百分比越高、颜色越深越好，表明有许多正确的预测。
我将使用 scikit 中的函数——学习绘制它:
超参数调整的超随机森林分类器的混淆矩阵
该模型很好地检测了美国黄松、三叶杨/柳树、白杨、花旗松、克拉姆霍尔茨等奇异的树木，但似乎对检测云杉/冷杉和黑松(覆盖类型 1 和 2)有些混淆。
特征重要性:
另一种方法是查看具有feature_importances_的特征重要性:一个介于 0 和 1 之间的数字，显示每个特征对预测的贡献。数字越高，说明贡献越大。
利用当前选择的功能、额外的树分类器和参数，前 10 个功能是:
其中 5 个是在本项目范围内创建的。这个列表强调了特征工程和选择的重要性。
7.用测试数据评估最佳模型(回答初始问题)
如何检测奇异的树(生成最终预测):
Kaggle 为比赛提供单独的训练和测试设备。直到现在，我都在训练台上工作。尽管如此，测试集用于最终预测，所以我将测试集和训练集在这里对齐，并将其输入超参数调整的额外树分类器
成功检测奇异树，准确率 78%。
哪里可以找到神奇的树:
云杉/冷杉、Lodgepole Pine 和 Krummholz 喜欢在 Rawah、Neota 和 Comanche Peak 荒野地区出没。
Cache la Poudre 荒野区是种植黄松和棉白杨/柳树的最佳地点。
如果你看到一个白杨，怀疑你是在 Rawah 或 Comanche。
花旗松是一种容易相处的物种，适合任何荒野地区。
8.总结和结论
在本文中，我重点介绍了端到端的机器学习工作流在有监督的多类分类问题中的应用。我从使用 EDA 理解和可视化数据开始，并形成了对 cover 类型数据集的见解。利用 EDA 的输出，我进行了特征工程，在其中我转换、添加和删除了特征。
额外树分类器很好地匹配了这个分类问题的准确性度量。随着超参数的调整，通过调整n_estimators参数，模型的准确度分数增加了 2 个点。解释模型结果显示了最重要的特征以及如何进一步提高准确性(通过更好地区分封面类型 1 和 2)。
项目背后的综合代码:
【https://github.com/cereniyim/Tree-Classification-ML-Model】Github 回购:T4
Kaggle 笔记本:https://www . ka ggle . com/cereniyim/fantastic-trees-where-to-find-how-detect-them/Notebook
感谢阅读，这种端到端的工作流程可以应用于任何机器学习问题，我鼓励你去尝试一下！
对于任何问题、评论或建设性反馈，您可以通过回复联系我， Twitter 或 Linkedin ！
使用 Elo 评级预测一级方程式赛车结果
体育分析
通过模型、数据管道和蒙特卡洛模拟建立预测
一级方程式赛车比赛通常感觉是可预测的:梅赛德斯在 2020 赛季赢得了过去的 5 个冠军和前 9 个 gp 中的 7 个。网上当然不乏 f1 预测——从专业专家到自封的纸上谈兵的“专家”，任何人都经常热衷于分享他们的预测。我的目标是对这个问题进行更深入的分析:一级方程式比赛能被算法预测吗？
在过去的两年里，我一直致力于预测 f1 资格赛的结果。从一开始，我的目标也是预测比赛结果，但这比预期的要复杂。排位赛结果相对更容易建模:在一级方程式排位赛中，获胜者是创下最快时间的车手。因此，司机的表现可以被认为是独立于其他人的:他们与时间竞争，只是间接的相互竞争。这使得预测结果更容易:每个司机的预测结果可以直接与他们的成绩进行比较。一级方程式比赛的建模更加复杂，因为车手之间直接竞争。但是随着一个封闭的夏天的到来，我决定是时候尝试一下赛车模型了。两位伟大的开发人员 Raiyan 和 Philip 加入了我的行列，共同探索建立一个预测 f1 比赛结果的模型。
在 F1 比赛中，完成时间并不重要:无论你以 20 秒还是 0.1 秒的优势赢得比赛，你都可以获得相同的分数。我们决定将比赛建模为成对参与者之间的一系列独立的、面对面的比赛，最终一名参与者获胜或失败(或者未解决，如果一方或双方没有完成比赛)。该模型不考虑完成时间或获得的分数，而是检查参与者赢得了多少场比赛，以及它与他们的预期分数如何匹配。这种方法背后的一个主要好处是它在“反常”比赛中的相对弹性:例如，在一个只关注终点位置的模型中，大多数领先者没有完成的比赛将导致较弱的车手获得不合理的大收益。一个关于头对头比赛的推理模型不会将退役的车手考虑在内，从而导致对车手分数的现实调整。
这些车手之间的竞争可以用一个叫做 Elo ratings 的概念来模拟。Elo 评级最初是为国际象棋选手排名而开发的。在 Elo 系统中，战胜一个更强的对手比战胜一个统计上更弱的对手获得更多的分数。每位车手的 Elo 评分在比赛结束后会根据他们在每场比赛中的表现与预期分数的对比情况进行上调或下调。预期分数是通过比较两个参与者的 Elo 评级计算出来的，它是一个介于 1 和 0 之间的值。预期分数 1 意味着车手将永远赢得比赛，而预期分数 0.5 意味着他们将赢得一半的时间。一级方程式赛车比赛的“真实”分数总是 1 或 0，取决于车手是赢还是输(如果一名或两名车手没有完成比赛，比赛不计算分数)。实际分数和预期分数之间的差异将应用于每位驾驶员的分数，K 系数为 4。
f1 比赛因排位赛而变得更加复杂，排位赛规定了比赛的发车顺序。因为车手的发车位置对他们的预期完成位置有着重大影响，所以参与者的 Elo 分数会根据他们的发车位置上下调整。网格调整后的 Elo 分数(gaElo)是用于预测车手比赛成绩的分数。与排位赛模型一样，预测是通过蒙特卡洛模拟生成的，在模拟过程中，每个车手、建造商和引擎的 gaElo 分数都会随机上调或下调。这里还模拟了一名车手退出比赛的可能性。通过多次运行模拟，可以计算出每个车手完成位置的百分比概率。
如前所述，比赛模型的预测取决于给定的起始网格。除了在周六的排位赛中决定的发车，我们还可以使用排位赛模型的模拟排位赛结果来生成比赛预测。这被称为“资格前预测”，它往往比“资格后预测”具有更大的分数方差。这是可以理解的，因为定性前的预测包含两倍的不确定性。构建允许两个模型顺序运行的管道并不容易，但是结果是非常值得的！
比赛模型产生的预测与现有的排位赛模型一起出现在项目的网站上。对 2020 年德国 Eifel GP 的预测是由当前的模型生成的，而早期的预测是由实验模型生成的。
马克斯·维斯塔潘的比赛预测截图。图片作者:Ville Kuosmanen
该模型的性能是根据两个基准来衡量的:一个随机模型，随机预测结果；一个模型，总是预测结果与起始网格相同。我们的评估显示，该模型的性能优于基准测试，尽管与基于网格的模型相比，它的性能仅略胜一筹，这表明该模型在准确性方面还有很长的路要走。比赛模型尤其难以预测由于发车区处罚或排位赛中的问题而在场地后部起步的快速车手的成绩，这是我们希望在未来解决的问题。
比赛模式仍然有很多需要改进的地方。尽管如此，我不相信比赛模型能够匹配排位赛模型的表现:f1 比赛从本质上来说比排位赛更随机(这是否是真的还有待证明)。我们想做的下一个功能包括通过模拟未来的比赛来预测整个锦标赛。就像连接排位赛模型和正赛模型一样，这是一个工程问题，需要建设新的数据管道。我还想通过在网站上展示对模型内部工作的一些见解来提高模型的可解释性(这是我在推荐系统中写过的一个主题)。调查赛季前测试的数据也很有趣:比赛和排位赛车型往往在每个赛季开始时表现最差，因为制造商之间的力量平衡被打破了。
如果你有兴趣了解我们的一级方程式预测项目，你可以查看 GitHub 上的项目页面(完整的开源代码)。虽然作为一名工程师，我希望看到一个可以完美预测每场比赛结果的模型，但作为一名观众，每隔一段时间就有一场意想不到的比赛获胜者被证明是错误的，这令人耳目一新！
使用 NLP +深度学习预测欺诈性新闻文章
利用 NLP、ML 和深度学习分析和预测欺诈性新闻文章。
生活在一个网上信息泛滥的时代，让我们容易相信欺骗性的新闻报道。这个问题在互联网上大量存在，其需求要求我们用机器学习和深度学习模型进行实验，以对潜在虚假的新闻进行分类。
这篇文章希望阐明“如何”利用自然语言处理、机器学习和深度学习来预测和识别欺诈性新闻文章。它还分享了从分析中得出的方法、结果和结论。
获取数据:
包含真假新闻文章的数据集是从 Kaggle 获得的。数据集的形状是 20800 行和 5 列。数据集中有 10413 篇真实文章和 10387 篇虚假文章，表明数据集是平衡的。“**标签“**栏表示该商品是真的还是假的； 1 表示假， 0 表示真。
图 1:数据快照。
清理数据:
由于文本数据是非结构化和杂乱的，因此需要进行文本预处理。一个助手函数被应用于数据，以移除换行符、新行、超链接、与符号、大于/小于符号、不间断空格、电子邮件、新行字符和分散注意力的单引号。创建了一个包含文章文本长度的新列“length”。
图 2:创建长度列后的数据集。
数据集还在“author”和“article_title”列中包含由“Unknown”替换的空值。长度少于 50 个字符的文章也被删除，理由是少于 50 个字符的文章不是文章。
探索性数据分析
名为**‘text _ polarity’的列被添加到数据集中。使用 TextBlob 库的情感极性函数对‘文本’**列返回文本的极性，从 1 到 1 表示从负到正。
观察真实文章和虚假文章之间文本极性的分布如何不同将是有趣的。下图 3A 和 3B 分别展示了真文和假文的文本极性分布。
左图 3A，右图 3B
这两类文章中文本极性的分布似乎是均匀的，有些文章是完全负面的，有些是完全正面的。
图 4 —相关矩阵
为了检查数据集中不同特征之间的相关性，使用了热图。数据集中不存在强相关，数据集中长度和标签之间的最大相关值为-0.12，表明非常弱的负相关。
使用 CountVectorizer 方法，绘制了去除和不去除停用词的前 20 个单字、双字和三元字。停用词指的是语言中最常见的词。在下面的例子中，图 5A 和 5B 显示了去除停用词后真实文章和虚假文章中的前 20 个三元模型。
左侧—图 5A，右侧—图 5B
推断统计分析:
接下来，进行单样本 t 检验，以检验总体平均值是否与某些假设值有显著差异。进行统计 t-检验以确定假文章的平均文本极性是否不同于所有文章的平均文本极性。
**假设检验:**是假的文章和全文字文章在文字极性的均值上有显著差异吗？
**无效假设:**无效假设将是假文章和所有文本文章之间的文本极性没有差异。
**备选假设:**备选假设将是假文章和所有文章之间的文本极性存在差异。
图 6
基于图 6 中的单样本 t-检验结果，所有文章的平均文本极性与虚假文章的平均文本极性之间存在显著差异。在 5%置信区间 0.0025 的低 P 值是拒绝零假设的良好指标。
建模
下一步，也是我们分析的关键，需要训练模型来预测一篇新闻文章是真是假。在训练模型之前，我们需要将文本转换成适合这些算法的输入。为此，我们将使用两种方法，称为计数矢量器和 TF-IDF(术语频率-逆文档频率)矢量器。
CountVectorizer 用于通过将文本语料库转换为表示相应语料库中字数的向量来对其进行标记化。它还允许我们删除文本中的停用词，并检查最流行的 N ’ unigrams，bigrams 和 trigrams。
相反，TF-IDF 是通过检查语料库并计算“术语频率”和“逆文档频率”获得的词频分数。术语频率是术语在文档中出现的频率，逆文档频率为文档中常见的词分配低权重。因此，TfidfVecotrizer 获取一组文档，并将它们分配给一个 TF-IDF 特性矩阵。然后，这种矢量化形式可以用作训练模型的输入。
这两种方法都将作为多项式朴素贝叶斯分类器和逻辑回归分类器的输入。逻辑回归被用作基准模型。基准模型是一种用于参考的模型，用来比较其他模型相对于它的表现。
没有为逻辑回归执行超参数调整，并且使用 0.01 的 C 值 对其进行训练。或者，对于多项式 NB 模型，使用交叉验证执行超参数调整，并发现 alpha = 0.1 是计数矢量器和 TF-IDF 矢量器模型中的最佳值。
LSTM(长短期记忆)神经网络是一种特殊类型的递归神经网络。简单来说，LSTM 网络有一些内部情境状态细胞，充当长期或短期记忆细胞。LSTM 网络的输出由这些单元的状态调制。当我们需要神经网络的预测依赖于输入的历史背景，而不是仅仅依赖于最后一个输入时，这是一个非常重要的属性。”[1] 阿萨德·莫瓦德，LSTM 神经网络的魔力 (2018)
创建 LSTM 模型时，每个文本中使用的最大字数等于 50，000， 最大字数 为 250。LSTM 模型的spatial drop out 1d值为 0.2， 脱落值 为 0.1， 经常性脱落值 为 0.1。
结果
为了检验结果，计算了不同输入(计数矢量器和 TF-IDF)和不同模型的准确性和 AUC-ROC 评分。在评估模型时，准确性本身可能是一个误导性的指标，因此还查看了 AUC-ROC 评分和 ROC 图。AUC-ROC 是一个很好的评估指标，因为它考察了模型正确分类的能力，也为我们提供了假阳性和假阴性的概念。如果 AUC 分数高，则模型在预测类别标签方面做得很好，即 0 表示 0，1 表示 1。在这种特殊情况下，如果 AUC 很高，这意味着模型预测真实的文章是真实的，而假文章是假的。
图 7A 显示了使用计数矢量器和 TF-IDF 转换数据的不同分类器的结果。图 7B 显示了相同分类器的 ROC 曲线。
**
左—图 7A，右—图 7B
用计数矢量器来表征数据的逻辑回归给了我们最好的结果。它实现了 95%的准确性和 95%的 AUC-ROC 评分。尽管 LSTM 模型获得了更高的 AUC-ROC 评分，但使用它可能会矫枉过正，因为逻辑回归在准确性方面表现更好，并且具有更高的 AUC-ROC 评分。此外，逻辑回归计算成本较低。
最差的模型性能也是通过用 TF-IDF 矢量器的逻辑回归实现的，准确度和 AUC-ROC 得分为 0.88。
模型评估
图 8A
图 8B
图 8C
除了仅使用准确性作为评估指标之外，我们还使用了 AUC-ROC 评分，它绘制了真阳性率与假阳性率的关系。我们还应该考虑假阳性和假阴性误差来评估我们的模型的性能。
假阳性(I 型错误):你预测文章是假的但却是真的。
假阴性(ⅱ型错误):你预测文章是真的但却是假的。
图 8A、8B 和 8C 显示了每个分类器的混淆矩阵和分类报告。
解释结果
重要的是，我们不仅要预测新闻文章是真是假，还要估计新闻文章是真是假的概率。这将允许我们标记假新闻文章，并显示它是假的概率。
用这种方法我们会遇到两种问题。首先，考虑新闻文章被模型预测为假的，但实际上，新闻文章是真实的。这被称为误报，这种错误可能代价高昂，因为它可能导致公司尴尬，还可能导致用户流失。不过，这可以通过建立一个机制来报告被错误标记为假的文章来快速纠正。这可能是这个项目发展的许多方向之一。
相反，考虑新闻文章被预测为真实但实际上是假的。随着越来越多的人接触到这篇假新闻，同时相信它是真实的，这篇文章将继续不受限制地被分享。这被称为假阴性，它有能力在敏感问题上播种不信任和分化在线用户。这对我们模型的有效性以及提供这些服务的公司的声誉更有害，因为它会导致人们不信任这个平台及其辨别真相和谎言的能力。为了缓解这个问题，可以建立一种机制，允许人们报告文章是否是假的，然后一个负责检查文章事实的团队根据他们的分析做出最终决定。这些只是一些可以实现的解决方案，它们既不意味着完美，也不声称是完美的。
你可以在这里找到代码为的完整笔记本。如果您有问题、可以提出的改进建议或想要讨论想法，请随时联系我们。
用 Python 和机器学习预测未来股市趋势
利用集成建模和交叉验证准确预测股票市场趋势
注来自《走向数据科学》的编辑: 虽然我们允许独立作者根据我们的 规则和指导方针 发表文章，但我们不认可每个作者的贡献。你不应该在没有寻求专业建议的情况下依赖一个作者的作品。详见我们的 读者术语 。
随着最近由于新冠肺炎疫情引起的股票市场的波动，我认为尝试和利用机器学习来预测股票市场的近期趋势是一个好主意。我对机器学习相当陌生，这是我的第一篇媒体文章，所以我认为这将是一个很好的开始和展示项目。
这篇文章讨论了关于数据科学的不同主题，即:数据收集和清理、特征工程，以及创建机器学习模型来进行预测。
注意:我之前在这篇文章的代码中使用了前瞻偏差，这产生了一些非常好的结果(非常好)。然而，它们是误导性的，我现在的目标是尝试用交叉验证来解决这个问题。
作者声明:本项目非财务或投资建议。这并不能保证它在大多数情况下都能提供正确的结果。因此，你应该非常小心，不要把它作为交易洞察力的主要来源。
你可以在我的 github 上的 jupyter 笔记本上找到所有代码:
此时您不能执行该操作。您已使用另一个标签页或窗口登录。您已在另一个选项卡中注销，或者…
github.com](https://github.com/lucasrea/StockForecast)
1.导入和数据收集
首先，我们包括了这个项目使用的所有库。我使用 yfinance API 收集了所有的历史股票市场数据。它直接取自雅虎财经网站，所以是非常可靠的数据。
然后我们定义一些用于数据检索和数据处理的常数。带有指示符号的列表有助于为我们的模型生成更多特征。
这里有一个链接，你可以在这里找到这些功能的实际名称。
熊猫实施的通用财务技术指标。这是正在进行的工作，错误是意料之中的，结果…
github.com](https://github.com/peerchemist/finta)
现在我们从 yfinance 中提取历史数据。我们没有太多可以使用的特性——除非我们找到一种方法来至少规范化它们或者从它们派生出更多的特性，否则这些特性不是特别有用。
“关闭”列中的数据
2.数据处理和特征工程
我们看到上面的数据是粗略的，包含了大量的时间序列峰值。它不是很平滑，模型很难从中提取趋势。为了减少这种情况的出现，我们希望在计算技术指标之前指数平滑我们的数据。
来自“Close”的数据，但已被平滑
我们可以看到数据平滑了很多。当计算技术指标时，有许多波峰和波谷会使它难以近似，或者难以提取趋势。它会让模型失控。
现在是时候计算我们的技术指标了。如上所述，我使用 finta 库结合 python 内置的 eval 函数来快速计算指标列表中的所有指标。除了正常的成交量，我还计算了一些不同平均长度的均线。
我去掉了像“开盘”、“高”、“低”和“收盘调整”这样的栏，因为除了指标之外，我们还可以用均线得到足够好的近似值。成交量已经被证明与价格波动有相关性，这也是我将其归一化的原因。
Index(['close', 'RSI', 'MACD', 'SIGNAL', '14 period STOCH %K','MFV', '14 period ATR', 'MOM', '14 period MFI', 'ROC', 'OBV_x', 'OBV_y', '20 period CCI', '14 period EMV', 'VIm', 'VIp', 'ema50', 'ema21', 'ema14', 'ema5', 'normVol'], dtype='object')
就在我们收集预测之前，我决定保留一点数据来预测未来的值。这一行捕获了与 7 月 27 日这一周的 5 天相对应的 5 行。
live_pred_data = data.iloc[-16:-11]
现在是这个项目最重要的部分之一——计算真值。如果没有这些，我们甚至无法训练一个机器学习模型来进行预测。
我们如何获得真理价值？这很直观。如果我们想知道一只股票何时上涨或下跌(希望赚一百万美元！)我们只需要展望未来，观察价格，以决定我们现在是应该买入还是卖出。嗯，有了这些历史数据，这正是我们能做的。
回到我们最初提取数据的表格，如果我们想知道 1993 年 3 月 29 日当天(收盘价为 11.4375)的买入(1)或卖出(0)决定，我们只需要向前看 X 天，看看价格是高于还是低于 1993 年 3 月 29 日的价格。因此，如果我们向前看 1 天，我们会看到价格上升到 11.5。因此，1993 年 3 月 29 日的真实价值是买入(1)。
由于这也是数据处理的最后一步，我们删除了指标和预测生成的所有 NaN 值，并删除了“close”列。
因为我们使用了 Pandas 的 shift()函数，我们从数据集的末尾丢失了大约 15 行(这就是为什么我在这一步之前捕获了 7 月 27 日这一周)。
3.模型创建
在训练我们的模型之前，我们必须将数据分成训练集和测试集。然而，由于时间序列的性质，我们需要小心处理这部分。如果我们随机化我们的训练测试集，我们可能会遇到前瞻偏差，这对于预测股票市场是不利的。这是因为你用模型已经看过的数据训练它。
为了防止这种情况，我们将使用一种称为交叉验证的不同技术来训练该模型。下图说明了我们将如何划分数据并测试模型的准确性。
时间序列中交叉验证的滑动窗口方法
首先，我们将使用多个分类器来创建一个集成模型。这里的目标是结合几个模型的预测，尝试提高可预测性。对于每个子模型，我们还将使用 Sklearn 的一个特性 GridSearchCV 来优化每个模型，以获得最佳结果。
首先，我们创建随机森林模型。
然后是 KNN 模式。
最后，我们创建投票分类器
一旦我们建立了模型，我们就可以把它和交叉验证放在一起。我们从编写一些代码开始，这些代码将允许我们用许多大小均匀的块来迭代我们的数据。
打印出我们的数据帧的索引，我们已经成功地分割了我们的数据，就像在滑动窗口图像中一样。
10 50
20 60
30 70
40 80
50 90
60 100
.
.
.
6820 6860
6830 6870
现在我们添加代码，将分区分成训练集和测试集。在 train_test_split 函数中设置 shuffle=False 非常重要——这是避免前瞻偏差的方法。
最后，我们合并我们的模型。
这样的样本运行看起来像…结果被夸大了。
rf prediction is [1\. 1\. 1\. 0\. 0\. 1\. 1\. 1\. 1\. 1\. 1\. 1.]
knn prediction is [1\. 1\. 1\. 1\. 1\. 1\. 1\. 1\. 1\. 1\. 1\. 1.]
ensemble prediction is [1\. 1\. 1\. 1\. 1\. 1\. 1\. 1\. 1\. 1\. 1\. 1.]
truth values are [1\. 1\. 1\. 0\. 1\. 1\. 1\. 1\. 1\. 1\. 1\. 1.]0.9166666666666666 0.9166666666666666 0.9166666666666666
最后的结果给了我们
RF Accuracy = 0.6732649071358748
KNN Accuracy = 0.67460899315738
Ensemble Accuracy = 0.6925708699902248
我们可以看到，通过使用集合建模，我们获得了更高的精度。与不使用交叉验证相比，结果的准确性要低得多。然而，这是解决这个问题的更正确的方法……大约 70%的准确率也不错！
4.结果验证
下一步，我们将使用预测模型预测 S&P500 的行为。我在 8 月 17 日的周末写这篇文章。因此，为了看看这个模型是否能产生准确的结果，我将使用本周的收盘数据作为预测的“真实”值。由于这个模型被调整为 15 天的窗口，我们需要输入 7 月 27 日这一周的数据。
七月二十七日->八月十七日
七月二十八日->八月十八日
七月二十九日->八月十九日
七月三十日->八月二十日
7 月 31 日-> 8 月 21 日
我们保存了将要在 live_pred_data 中使用的周。
live_pred_data.head()
以下是我们将要进行预测的五个主要日期。看起来模型预测价格每天都会上涨。
让我们用实际结果来验证我们的预测。
del(live_pred_data['close'])
prediction = ensemble_model.predict(live_pred_data)
print(prediction)[1\. 1\. 1\. 1\. 1.]
结果
7 月 27 日:322.78 美元——8 月 17 日:337.91 美元
7 月 28 日:321.74 美元——8 月 18 日:338.64 美元
7 月 29 日:323.93 美元——8 月 19 日:337.23 美元
7 月 30 日:323.95 美元——8 月 20 日:338.28 美元
7 月 31 日:325.62 美元——8 月 21 日:339.48 美元
正如我们从实际结果中看到的，我们可以确认模型的所有预测都是正确的。然而，决定股票价格的因素有很多，所以说这个模型每次都会产生相似的结果是幼稚的。然而，在相对正常的时间段内(没有导致股票市场波动的重大恐慌)，该模型应该能够产生良好的结果。
5.摘要
总结一下我们在这个项目中所做的，
- 我们已经收集了用于分析和特征创建的数据。
- 我们已经使用熊猫来计算许多模型特征，并产生干净的数据来帮助我们进行机器学习。使用熊猫创建预测或真实值。
- 利用交叉验证来避免前瞻偏差。训练许多机器学习模型，然后使用集成学习将它们结合起来，以产生更高的预测精度。
- 确保我们的预测与真实世界的数据相符。
通过这个项目，我学到了很多关于数据科学和机器学习的知识，希望你也一样。作为我的第一篇文章，我喜欢任何形式的反馈来帮助提高我作为程序员和数据科学家的技能。
感谢阅读:)
预测未来战争
在 Unsplash 上由 Stijn Swinnen 拍摄的照片。
来自开放数据和机器学习的见解
我知道你在想什么:战争是罕见而复杂的事件，人们不能指望考虑到它们的全部复杂性。你是对的，它们源于错综复杂的政治、经济和历史原因，但不要忘记厚厚的随机性外衣，因此应该谨慎对待。甚至政治专家也一直在争论是什么导致了冲突，有时是在事件发生几十年后。冲突可能是黑天鹅，本质上不可预测的改变游戏规则的事件(即使专家经常声称事后发现了明显的警告信号)。但这并不意味着我们不能通过使用机器学习来评估导致战争的氛围。
用经济指标来预测冲突并不是什么新想法，关于这个问题已经做了很多研究 。结果是有希望的，但他们使用给定年份的数据来预测同一年一个国家的状态，这有点违背直觉。如果你已经知道了国家的经济状况，那么为什么还需要“预测”它是否处于战争状态呢？信息已经有了。这迫使他们使用联合国等组织提供的预测来应用他们的模型。这使得他们的预测依赖于他们。在这里，只有过去和现在的数据被用来预测一个国家在不久的将来的状态。
使用的功能
为了将这些指标与冲突联系起来，使用了 UCDP/PRIO 武装冲突数据集 ，其中列出了 1946 年至 2018 年的所有武装冲突。根据冲突的致命程度和所反对的群体类型，冲突被分为不同的类型。就特征而言，采用了主要来自世界银行的年度时间序列形式的经济和社会指标。目标特征是通过增加表明一个国家在未来五年内是否会陷入冲突的变量而获得的。在这里，我们跳过预处理阶段，更详细的描述见报告。最终数据集包含 12 900 行(60 年* 215 个国家)。
分析因果关系
对于多变量时间序列，皮尔逊相关是不合适的，当前冲突状态可以是其他变量的过去值的函数。这就是应用格兰杰因果关系检验的原因。它决定了一个时间序列如何影响另一个时间序列。每年的汇总版本用于计算测试。
经济指标和冲突的格兰杰因果关系测量(X 影响 Y)。鸣谢:作者。
小冲突受到债务、账户余额和自然资源租金的强烈影响。总的来说，冲突受到婴儿死亡率的强烈影响，婴儿死亡率是贫困引发冲突的一个指标。国家间的冲突是由国内生产总值和账户余额引发的。外国直接投资与新冲突的爆发或出现之间也有着密切的联系。这些观察表明，经济是有影响的，但是经济的不同方面与不同类型的冲突有关。
模型
最好的结果不是通过特定于时间序列的模型(VAR 模型)获得的，而是通过将每个国家的每一年作为一个数据点。一个有 100 棵树并且最大深度为 15 的随机森林获得了如下所示的结果。
随机森林模型的结果:准确率为 92%，F1 值为 74%，精确度为 71%，召回率为 76%。鸣谢:作者。
用于区分的特征如下所示。
根据辨别能力对每个特征进行排序。鸣谢:作者。
对该模型最有影响的特征是腐败、婴儿死亡率和自然资源租金。这意味着该模型的预测基于政权类型(腐败)、贫困(婴儿死亡率)和对资源的依赖(自然资源)。
预言
最后，有趣的部分是，最近几年(2018 年和 2019 年)的数据被输入到模型中，以查看世界上哪个国家将在未来五年(2019 年至 2024 年)陷入新的冲突。结果如下图所示。
我们的模型给出的每个国家在 2019-2024 年经历新冲突的概率。鸣谢:作者。
许多更有可能面临新冲突的国家已经经历了动荡。但像巴基斯坦、印尼这些相对稳定的国家，近期爆发新冲突的概率很大。伊朗最近因制裁和最近的政治气氛而承受巨大压力，与其他国家相比，它也很有可能陷入新的冲突。
如果只考虑主要的冲突，即战争，我们会得到这些预测。
我们的模型给出的每个国家在 2019-2024 年面临新的重大冲突的概率。鸣谢:作者。
在这种情况下，只有叙利亚和索马里这两个已经处于非常不稳定状态的国家很有可能在未来几年面临新的战争。
尽管战争远不是一个确定性的系统，但机器学习是评估一个国家历史上关键时期的有用工具。
完整的项目可以在这里访问。
从好书到名著
使用 python 来预测是什么让书籍变得伟大
什么
这是 2020 年 7 月 30 日从 Goodreads 的 API 中提取的前 50，000 个图书 id 的数据集。几千个 id 没有通过，因为图书 id 被更改，URL 或 API 被破坏，或者信息以非典型格式存储。
为什么
从读者的角度来看，书籍是一种多小时的学习和休闲的承诺(他们不会称之为好的免费阅读)。从作者和出版商的角度来看，书是一种生活方式(也有一些学习和休闲)。在这两种情况下，知道哪些因素解释和预测伟大的书籍将节省您的时间和金钱。因为虽然不同的人有不同的品味和价值观，但了解一本书的总体评价是一个明智的起点。你可以随时更新它。
环境
在虚拟环境中工作是很好的实践，一个有自己的库和版本的沙箱，所以我们将为这个项目制作一个。有几种方法可以做到这一点，但是我们将使用 Anaconda 。要使用 Python 3.7 创建并激活名为“gr”(Goodreads)的 Anaconda 虚拟环境，请在您的终端或命令行中运行以下命令:
装置
您应该在提示符的左侧看到“gr”或您命名的环境。如果是这样，请运行以下命令。Anaconda 将自动安装这些包的任何依赖项，包括 matplotlib、numpy、pandas 和 scikit-learn。
进口
数据收集
我们使用由 Michelle D. Zhang (代码和文档在这里)制作的 Goodreads API 的轻量级包装器提取前 50，000 个图书 id 及其相关信息，然后将每个 id 作为字典写入名为book_data的 JSON 文件。
数据清理
我们将在下面定义和描述一些关键函数，但是我们稍后将在一个大的 wrangle 函数中运行它们。
威尔逊下界
基于 20 条评论的 4 星评价和基于 20000 条评论的 4 星评价是不相等的。基于更多评论的评级具有较少的不确定性，因此是对“真实”评级的更可靠的估计。为了正确定义和预测伟大的著作，我们必须通过惩罚不确定性来改变average_rating。
我们将通过计算威尔逊下限来实现这一点，其中我们估计特定评级的置信区间，并将其下限作为新的评级。基于数万条评论的评级几乎不会受到影响，因为它们的置信区间很窄。然而，基于更少评论的评级具有更宽的置信区间，并且将被降低更多。
注:我们修改了公式，因为我们的数据是从 5 分制计算的，而不是 Wilson 描述的二进制。具体来说，我们将average_rating减 1，作为真实非膨胀评级的保守估计，然后将其归一化。如果这个惩罚太重或太轻，随着时间的推移，更多的评级将分别提高或降低该书的评级。换句话说，随着更多的信息，这种调整是自我纠正的。
流派
Goodreads 的 API 返回了“书架”,它包含了像“科幻小说”这样的真实类型和像“阅读”这样的用户创建的类别。在提取数据时，我们只提取了 5 个最受欢迎的货架，以限制这种清理；在这里，我们会完成任务。
经过一些检查，我们看到这些子串代表了非流派书架的主体。我们将使用正则表达式过滤掉它们。注意:我们在正则表达式中使用了两个字符串，所以这一行不会被切断。圆括号内的相邻字符串在编译时连接在一起。
多功能一体清洗
现在，我们将构建并运行一个函数来处理数据集。这样，清理更具可重复性和可调试性。
比较未调整和调整后的平均评分
从数字上看，平均值(蓝色)和中值(绿色)趋势的中心度量值略有下降，方差显著下降。
视觉上，我们可以看到评级调整在更平滑和更宽的分布中(尽管注意 x 轴被截断)。这是因为消除了没有评级或评级很少的异常书籍，并缩小了具有高度不确定性的评级。
作者图片
Unadjusted mean: 3.82
Unadjusted median: 3.93
Unadjusted variance: 0.48
作者图片
Adjusted mean: 3.71
Adjusted median: 3.77
Adjusted variance: 0.17
数据泄露
因为我们的目标来源于评级，使用评级来训练我们的模型实际上是用目标来训练。为了避免扭曲模型，我们必须删除这些列。
也有可能review_count有点漏，但它似乎更像是流行度的代理，而不是伟大，就像流行(ular)歌曲通常不被认为是经典一样。当然，如果它的排列重要性高得可疑，我们会重新考虑这个问题。
分割数据
我们将进行 85/15 的训练测试分割，然后重新分割我们的训练集，使验证集与测试集的大小大致相同。
(20281, 12) (20281,) (4348, 12) (4348,) (4347, 12) (4347,)
评估指标
对于这种不平衡的分类，准确性(正确预测/总预测)可能会产生误导。对于这个分数来说，没有足够的真阳性作为模型性能的最佳度量。因此，我们还将使用 ROC AUC，曲线下的受试者操作者特征面积。这是一张彩图，由马丁·托马斯提供。
Martin Thoma 绘制的 XKCD 风格的 ROC AUC 图，在知识共享 CC0 1.0 通用公共领域专用下提供
ROC 曲线是分类模型的真阳性率(TPR)与其假阳性率(FPR)的关系图。ROC AUC 是从[0，1]到该曲线右下方的面积。因为最佳模型性能使真阳性最大化，假阳性最小化，所以这个 1x1 图中的最佳点在左上角，曲线下面积(ROC AUC) = 1。
对于像great这样的不平衡类，ROC AUC 作为一个指标优于准确性，因为它更好地反映了真阳性和假阳性之间的关系。它还描述了分类器在其所有值上的性能，为我们提供了更多关于模型何时何地改善、停滞或受损的信息。
适合模型
预测名著是一个二元分类问题，所以我们需要一个分类器。下面，我们将对一个线性模型(逻辑回归)和两个基于树的模型(随机森林和 XGBoost)进行编码、估算和拟合，然后将它们相互比较并与多数基线进行比较。我们将计算它们的准确性和 ROC AUC，然后进行可视化。
多数类基线
首先，从结构上来看，great的书是 Wilson 调整后排名前 20%的书。这意味着我们的多数类基线(没有书是伟大的)有 80%的准确性。
第二，这个“模型”不会改进、停滞或受损，因为它从一开始就没有洞察力。随机选择的阳性和随机选择的阴性将被同等对待。在其他研究中，其 ROC AUC = 0.5。
Baseline Validation Accuracy: 0.8
Baseline Validation ROC AUC: 0.5
逻辑回归
现在，我们将使用交叉验证拟合线性模型，重新计算评估指标，并绘制混淆矩阵。
Baseline Validation Accuracy: 0.8
Logistic Regression Validation Accuracy: 0.8013
Baseline Validation ROC AUC: 0.5
Logistic Regression Validation ROC AUC: 0.6424
逻辑回归混淆矩阵
作者图片
随机森林分类器
现在我们将使用基于树的模型和 bagging 做同样的事情(BootstrapAGGregation)。
Baseline Validation Accuracy: 0.8
Logistic Regression Validation Accuracy: 0.8013
Random Forest Validation Accuracy: 0.8222
Majority Class Baseline ROC AUC: 0.5
Logistic Regression Validation ROC AUC: 0.6424
Random Forest Validation ROC AUC: 0.8015
随机森林混淆矩阵
作者图片
XGBoost 分类器
现在，我们将使用另一个基于树的模型做同样的事情，这次使用 boosting。
Baseline Validation Accuracy: 0.8
Logistic Regression Validation Accuracy: 0.8013
Random Forest Validation Accuracy: 0.8245
XGBoost Validation Accuracy: 0.8427
Majority Class Baseline ROC AUC: 0.5
Logistic Regression Validation ROC AUC: 0.6424
Random Forest Validation ROC AUC: 0.8011
XGBoost Validation ROC AUC 0.84
XGBClassifier 在准确性和 ROC AUC 方面表现最佳。
绘制并比较模型的 ROC AUC
下面，我们看到逻辑回归在获得高 ROC AUC 方面远远落后于 XGBoost 和 Random Forests。在前两者中，XGBoost 最初优于 RandomForest，然后两者大致收敛在 FPR=0.6 左右。然而，我们在右下方的图例中看到，XGBoost 的 AUC 最高，为 0.84，其次是随机森林，为 0.80，逻辑回归为 0.64。
用不太专业的语言来说，XGBoost 模型最擅长将好书归类为好书(真阳性)，而不将不太好的书归类为好书(假阳性)。
作者图片
排列重要性
一个直观的识别某样东西是否重要以及重要到什么程度的方法是看当你拿走它时会发生什么。在不受时间和金钱约束的情况下，这是最好的。
但是在有实际约束的现实世界中，我们可以用置换来代替。我们不是通过删除列值来消除它们，而是通过随机化来消除列的信号。如果该列真的是一个预测特性，那么它的值的顺序就很重要，如果不破坏关系的话，将它们移动会大大削弱关系。因此，如果特征的预测能力没有真正受到伤害，或者随机化甚至有所帮助，我们可以得出结论，它实际上并不重要。
让我们仔细看看 XGBoost 模型的排列重要性。我们将不得不改装它以与 eli5 兼容。
排列重要性分析
作者图片
正如我们在开始时假设的那样，review_count很重要，但并没有高到令人怀疑的程度。这似乎没有上升到数据泄露的程度。这意味着，如果你想知道下一本书该读什么，一个有用的指标是它有多少评论，代表有多少人读过它。
我们看到genres是 XGBoost 模型中 ROC AUC 的第二个最重要的特征。
author排在第三，这令人惊讶，或许也有点令人担忧。因为我们的测试集并不大，这个模型可能只是识别出按 Wilson-adjusted 术语评价最高的作者，比如 J.K .罗琳和苏珊·科林斯。更多的数据将有助于检验这一理论。
第四个是num_pages。我认为这个数字会更高，原因有二:
- 很长的书的评分似乎有一点向上的倾向，因为愿意开始和完成它们的人会给它们更高的评分。冗长的篇幅筛选掉了不太感兴趣的边缘读者，他们可能一开始就不会给这本书很高的评价。
- 阅读并炫耀你正在阅读或已经阅读了长篇大论是社会地位高的标志。典型的例子:无限的玩笑。
外卖食品
我们已经看到了如何收集、清理、分析、可视化和建模数据。一些可操作的要点是，什么时候和谁出版一本书并不重要，但它的评论数很重要——评论越多越好。
为了进一步分析，我们可以分解genres和authors来找出哪些被评为最高。现在，快乐阅读。
用随机森林预测幸福
使用机器学习探索影响人们幸福水平的因素
迈克·彼得鲁奇在 Unsplash 上的照片
为什么我们试图预测幸福？能够预测幸福意味着我们能够操纵或试图改善某些因素，以增加我们自己的幸福，也可能增加政府的国民幸福。我发现随机森林(RF)是最简单和最有效的软件包，所以让我们开始吧！
内容:
- 数据
- 随机森林模型
- 数据清理
- 培训和测试
- 特征重要性
- 修改变量的数量
- 评估模型
数据:
从 #WorldValuesSurvey 获得的数据包含了> 290 个问题&在去除了缺失的幸福水平数据后，由大约 69k 个回答组成。是跨年度的跨国调查，问卷可以在网站上找到。特别是，我们将关注 2017 年至 2020 年的数据集。数据集的大小使其成为机器学习的最佳选择。
随机森林模型:
首先，我们将使用 RF 分类器*，因为我们希望机器能够预测群体的幸福水平(非常幸福、相当幸福、不太幸福、一点也不幸福)。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn import metrics
数据清理:选择数据
让我们先从问题栏开始，去掉 Q46 关于幸福水平的回答中的负值。
var="Q46"
df=df[df.columns[32:349]]
df=df[df[var]>0]
*负值表示受访者表示不知道、没有答案、没有被问到或回答缺失。这些值将使机器更难对它们进行分类，因为它增加了类别的数量，而不是我们所寻找的。
剩余的数据集如下所示:
进一步的数据清理:
下一个问题是，我们必须处理其他列中缺失的值。有 3 个选项可供考虑:
- 用 0 替换丢失的值
- 用平均值替换缺失的值
- 删除缺少值的行(数据集变为空)。
由于第三个选项不可行，我们必须检查选项 1 或 2 中哪一个会给出最高的精度。在这种情况下，我发现用 0 替换会更准确。
df.fillna(0, inplace=True)
准备列车标签:
现在，我们为机器设置“标签”,以识别我希望它预测的特征，并将数据分为训练集和测试集。
train_labels = pd.DataFrame(df[var])
train_labels = np.array(df[var])
train_features= df.drop(var, axis = 1)
feature_list = list(train_features.columns)
train_features = np.array(train_features)
train_features, test_features, train_labels, test_labels = train_test_split(train_features, train_labels, test_size = 0.25, random_state = 42)
训练并测试模型:
培训和测试的过程很简单。为了提高预测能力和/或模型速度，我们可以简单地修改 RF 分类器中的参数。
增加精度:
n_estimators —多数表决前算法构建的树的数量
max _ features 随机森林考虑分割节点的最大要素数
min _ sample _ leaf 分割内部节点所需的最小叶子数。
增加速度:
n_jobs —允许使用的处理器数量。如果= 1，只使用一个处理器。If =-1，无限制
random _ state 使模型的输出可复制，即给定相同的超参数和训练数据，总是产生相同的结果
oob_score:随机森林交叉验证方法
rf=RandomForestClassifier(n_estimators = 1000, oob_score = True, n_jobs = -1,random_state =42,max_features = “auto”, min_samples_leaf = 12)
rf.fit(train_features, train_labels)
predictions = rf.predict(test_features)
print(metrics.accuracy_score(test_labels, predictions))
该模型需要 1.3 分钟来训练~52k 训练行和> 290 列，并需要 1 秒来测试。准确率为 63.70% 。如果我们选择用平均值来填充缺失值，精确度将是 63.55% 。但重要的是找出是什么影响了机器的预测，因为这些是我们想要研究的变量。我们当然不能指望每个人都回答 290 多个问题，或者试图在所有 290 个方面努力提高幸福感(这将花费很多)。因此，我们将着眼于特性的重要性。
特性重要性:
如果您还记得，feature_list 包含除 Q46 之外的所有其他变量的列。目标是了解影响预测的变量。
importances = list(rf.feature_importances_)
feature_importances = [(feature, round(importance, 2)) for feature, importance in zip(feature_list, importances)]
feature_importances = sorted(feature_importances, key = lambda x: x[1], reverse = True)
[print('Variable: {:20} Importance: {}'.format(*pair)) for pair in feature_importances]
x_values = list(range(len(importances)))# Make a bar chart
plt.bar(x_values, importances, orientation = 'vertical', color = 'r', edgecolor = 'k', linewidth = 1.2)
# Tick labels for x axis
plt.xticks(x_values, feature_list, rotation='vertical')
# Axis labels and title
plt.ylabel('Importance'); plt.xlabel('Variable'); plt.title('Variable Importances');
特征重要性总和为 1，我们注意到，与其他变量相比，某些变量对预测的影响更大，几乎每个变量都有某种形式的影响，尽管由于变量太多而非常小。接下来的事情是继续改进我们的模型，让我们更好地理解幸福。
修改变量数量:
让我们来看看前 20 个特征，并用这 20 个变量(+ var 本身)建立一个新模型。我们将重复数据清理和相同的 RF 模型。我得到了 **64.47%的准确率。**如果我们选择用平均值替换缺失值，则精确度将为 64.41% 。这里令人惊讶的是，随着变量数量的减少，模型变得更加精确(从 63.70% 到 64.47% )。这可能是因为其他变量在模型中产生噪声，导致模型不太准确。
让我们再来看看特性的重要性:
这一次，可以更清楚地看出哪些变量更重要。你可以参考 WVS 网站上的调查问卷了解更多详细信息。我将总结一下这些问题所涉及的主题。
评估模型:
让我们看看前 200 个测试值的实际值与预测值的图表。为了更好地了解整个测试集，让我们简单地计算一下预测值和实际值之间的差异(预测值减去实际值)。
该模型在预测幸福水平方面似乎是消极的多于积极的，但在其他方面仍然被认为是平衡的！
感悟:
我所做的是检查 WVS 中 290 多个与幸福水平更相关的关键问题。这将意味着我们在研究幸福时可以试着特别关注这些方面。
在调查问卷中，我们还会注意到 Q261 和 Q262 是同一个东西(年龄和出生年份)，因此我们可以删除其中一个以包含另一个特征。对于问题 266，267，268(回答者和父母的出生国)，它们似乎是重复的，但并不完全相同，因为移民/跨文化婚姻可能会发生。尽管如此，我们可以考虑删除其中两个，因为这种情况很少发生。
一般话题有:
个人层面:
生活满意度、健康、财务、自由、年龄、安全、宗教、婚姻、家庭。
国家层面:
国家、对腐败的看法、民主/政治影响力、民族自豪感
特别是，健康、财务和年龄是机器认为最重要的特征。从这个意义上说，个人层面的因素比国家层面的因素对一个人的幸福水平有更大的影响。
然而，我注意到 WVS 没有关于睡眠时间的数据，这是我在之前的文章中观察到的一个关键因素。尽管如此，它仍然非常有用，因为我们可以考虑这些方面进行进一步的分析！我会带着对这些方面和幸福之间的相关性的更多见解回来，以确定我们如何才能提高我们的幸福水平。在那之前，记得保持快乐！
预测危险的地震颠簸第二部分:训练和调整监督最大似然分类器和模型性能分析
本文演示了使用不同的监督分类器预测危险的地震颠簸、调整模型超参数、准确性悖论以及理解“业务问题”对性能评估的重要性
肖恩·麦克伦登在 Unsplash 上的照片
介绍
在我的上一篇文章中，关于来自 UCI 数据档案库的地震颠簸数据集，我应用了特征工程的基本数据分析技术和不平衡数据集的测试序列分割策略。在本文中，我展示了如何应用监督机器学习算法(KNN、随机森林和 SVM)进行预测，并调整超参数(使用 GridSearchCV)以获得最佳结果。性能评估的结果也清楚地展示了准确性悖论，我将在下面的章节中详细阐述。此外，它还展示了为什么理解“业务问题”对于选择最佳模式至关重要。
完整的笔记本可以在我的 GitHub 库中找到
我从…
创建一个文件夹来保存我的模型。将多种算法用于数据时，最好在训练和调整后保存模型。为此，我使用 datetime python 模块创建了一个文件夹。
**import** **datetime****def** model_store_location():
**return** "model-store-**{}**".format(datetime.date.today())model_store_location = model_store_location()
print(model_store_location)!mkdir {model_store_location}
每个性能最好的模型都将存储在该文件夹中。
接下来，我建立了一个基线
这是一个二元分类任务，我选择了最简单的分类模型，K-最近邻(KNN)来开始预测建模。KNN 是最简单的分类器之一，它基于前 K 个相似的可见实例中标签的最大数量来将标签分配给不可见实例。
为了改进模型，需要一个基线，随后的模型将与该基线进行比较。因此，它被称为“基线模型”。因此，对于基线模型，我使用默认参数初始化了 KNN 分类器:
model = KNeighborsClassifier()
接下来，我首先检查了性能，使用 10 次分割的 StratifiedKFold，保持评分标准“准确性”。
skf = StratifiedKFold(n_splits=10)
作者图片
对于所有的折叠，精确度在 0.92 到 0.93 之间，并且数值的接近性通过分数的标准偏差来显示。
这里的准确率真的很高。现在，让我们在训练数据中拟合模型，并检查模型的性能—
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
我添加了一个数据框来保存每个模型的超参数和性能，如下所示:
df_results = pd.DataFrame(index = ['scoring_technique',
'algorithm', 'n_neighbors',
'weights', 'leaf_size',
'accuracy', 'precision',
'recall', 'f1-score'])
并添加了我的基线模型的性能:
**from** **sklearn.metrics** **import** accuracy_score, precision_score, recall_score, f1_scoredf_results['Baseline'] = ['None', model.algorithm,
model.n_neighbors,
model.weights,
model.leaf_size,
accuracy_score(y_test, y_pred),
precision_score(y_test, y_pred),
recall_score(y_test, y_pred),
f1_score(y_test, y_pred)]
作者图片
这就是 KNN 的表演看起来超越了准确性。这是一个典型的准确性悖论的例子。由于大多数地震凸起是无危险的，因此真负值(TN)比 TP、FP 和 FN 的总和要大得多， 从而在数字上提高了准确性，并创造了预测正确的假象 ，但正如精度和召回率所示，危险地震凸起识别率低于 30%。
二元分类性能指标—
= TP+TN/(TP+FP+TN+FN)
精度 = TP / (TP + FP)
回忆 = TP / (TP + FN)
F1-Score= 2 * Precision x Recall/(Precision+Recall)= 2TP/(2TP+FP+FN)**
因此，基于基线性能度量，只有 25%的预测危险隆起是正确的，只有 5%的实际危险地震隆起被模型正确识别。93%的数据实际上包含无危险的地震隆起，因此cross_val_score 中的精确数字是有意义的，也在 93%左右。
然而，正如您所看到的，精确度和召回率都受到了影响，让我们看看如何使用GridSearchCV来调优超参数
GridSearchCV 用于调整 KNN 超参数
马库斯·温克勒在 Unsplash 上拍摄的照片
GridSearchCV 是一个函数，它接受超参数的可能值或范围，运行指定的超参数值的所有组合，并根据提到的评分标准计算性能。这个评分标准应该基本上与业务价值或您的目标一致。
要查看提供了哪些评分标准，此代码有助于:
***import** **sklearn**
sorted(sklearn.metrics.SCORERS.keys())*
该代码将生成可用评分方法的完整列表。
接下来，我为 KNN 定义了超参数，如下面的代码所示，并创建了一个名为param_grid 的字典，其中包含 KNN 超参数的参数，如 Scikit-Learns 关于 KNN 的文档中所定义的键，而值包含相应的值列表。
*n_neighbors = [1, 2, 3, 4, 5]
weights = ['uniform', 'distance']
algorithm = ['ball_tree', 'kd_tree', 'brute']
leaf_size = [10, 20, 30 , 40, 50] #param_grid = dict(n_neighbors=n_neighbors,
weights=weights,
algorithm=algorithm,
leaf_size=leaf_size)*
接下来，我使用模型、最佳超参数搜索的网格、交叉验证类型、作为精度的评分度量来初始化 GridSearchCV，并且我还使用了refit参数来确保最佳评估的估计值作为best_estimator_可用，这是当训练集可用于进行预测时的拟合模型。Verbose决定您希望在屏幕上显示多少日志信息。
*grid = GridSearchCV(estimator=model,
param_grid=param_grid,
cv=StratifiedKFold(shuffle=**True**),
scoring=['precision'],
refit='precision',
verbose=10)*
我调用 fit 方法来启动超参数调优过程。
*grid_result = grid.fit(X_train, y_train)*
基本上，grid_result包含了所有的输出，从拟合时间到个人得分，以及基于传递的参数评估的最佳模型。grid_result.best_estimator_包含最佳拟合模型的超参数的选定值，根据评分标准进行优化。
我试验了不同的scoring
型号 1 — scoring:精度
型号 2 — scoring:召回
型号 3 — scoring : f1
以下是模型性能的总结:
作者图片
如果您在此网格中注意到，每个模型的评分指标的最佳可能结果属于在整个网格中指定的相应scoring指标，即，当参数为“precision”时，任何其他模型的“precision”性能都不如模型 1。同样，对于“召回”，网格中的最佳由模型 2 和模型 3 共享。其实他们的表现是一样的。
在模型 1 中，其精度为 50%，即 50%的预测的危险地震颠簸是正确的，有可能评论说该模型只不过是掷硬币。但是抛硬币是公平的，而我们的数据集却不是，它是一个不平衡的数据集。在一次公平的掷硬币中，正面和反面同样可能发生，但是对于这个数据集，危险的和非危险的地震颠簸不会发生。
在模型 2 和模型 3 中，召回率最高，可以解释为正确识别的实际危险地震碰撞的 14%。精度也是 16%。f1 的成绩远远好于 Model 1，相当不错。模型 1 的低 f1 分数是因为 2%的不良回忆，这意味着该模型只能预测实际危险地震颠簸的 2%,这是 no bu eno。
随机森林分类器能赢得这场挑战吗？
我选择的下一个模型是随机森林分类器，竞赛获奖的 ML 模型。这是一个健壮的 ML 模型，其中构建了多个决策树来集成它们的输出。最终预测是所有单个决策树预测的函数。这就是随机森林模型表现更好的原因。
我从初始化随机森林分类器开始:
*from sklearn.ensemble import RandomForestClassifier
model = RandomForestClassifier()*
接下来，我构建了一个网格来寻找最合适的参数:
*param_grid = {
‘bootstrap’: [True],
‘max_depth’: [80, 90, 100, 110],
‘max_features’: [2, 3, 4],
‘min_samples_leaf’: [3, 4, 5, 6],
‘min_samples_split’: [8, 10, 12],
‘n_estimators’: [100, 200, 300, 500]
}*
我还再次为“得分”选择了“精确度”。以下是完整的代码:
*from sklearn.ensemble import RandomForestClassifier
model = RandomForestClassifier()param_grid = {
‘bootstrap’: [True],
‘max_depth’: [80, 90, 100, 110],
‘max_features’: [2, 3, 4],
‘min_samples_leaf’: [3, 4, 5, 6],
‘min_samples_split’: [8, 10, 12],
‘n_estimators’: [100, 200, 300, 500]
}grid = GridSearchCV(estimator=model,
param_grid=param_grid,
cv=StratifiedKFold(shuffle=True),
scoring=['precision'],
refit='precision',
verbose=10)
grid_result = grid.fit(X_train, y_train)file_name = 'seismic_bump_rf_model.sav'
joblib.dump(model, model_store_location + file_name)print(grid_result.best_estimator_)*
由此得出的模型是—
然而，在预测测试集时，精度为零。
作者图片
训练和调整支持向量机
由于随机森林模型真的让我失望，我想到了使用支持向量机(SVM)。支持向量机属于分类算法的核方法组，用于拟合决策边界。分类基于决策边界的哪一侧是数据点。
支持向量机的计算开销很大，可以通过以下网格搜索设置进行演示:
*from sklearn.svm import SVC
model = SVC()param_grid1 = {'C': [0.1, 1, 10, 100, 1000],
'gamma': [1, 0.1, 0.01, 0.001, 0.0001],
'kernel': [ 'linear', 'poly', 'rbf', 'sigmoid']}*
安装这个网格花了超过 10 个小时，实际上我不得不中断。我检查了多边形内核有更好的结果。所有线性核都很快，但精度为零。总的来说，这个网格相当于 500 次拟合，每 100 个候选对象 5 次折叠(5 * 5 * 4)。
因此，我构建了一个新网格，如下所示:
*param_grid2 = {'C': [1],
'gamma': [1],
'kernel': ['poly', 'rbf', 'sigmoid']}*
这个网格总共有 15 次拟合，花了大约 3.5 个小时来完成调整过程。
现在，按照同样的过程，我看了看最适合的模型，如下所示:
*grid_result.best_estimator_***Output:
SVC(C=1, gamma=1)****
默认的核是“rbf”或径向基函数(RBF)。要了解更多关于内核函数的信息，请参考 Scikit-learn 文档。
*grid_result.best_estimator_.kernel***Output:
rbf****
以下是 SVM 的最终性能结果:
作者图片
就精度而言，这看起来比迄今为止所有其他模型都要好。然而，召回只是稍微好一点。
现在，我们来对比总结一下…
每个分类器的最佳精度(图片由作者提供)
在上表中，SVM 的精确度略好，但精确度有了很大的提高，67%的危险地震隆起的实际数量被正确预测。然而，召回率是 6%，即实际危险的地震颠簸中只有 6%是正确的。KNN 车型的 f1 得分低于 SVM 车型。
总的来说，由于无法识别危险的地震碰撞带来的风险因素，该问题陈述需要更好的回忆。召回中唯一表现较好的模型是 KNN，其权重一致，n _ neighbours 为 1，leaf_size 为 10。
业务需求的重要性因情况而异。有时回忆是首选(像这里)，有时精确是首选。如果对预测的结果采取行动是有成本的，那么精确度就很重要，也就是说，你会想要更少的假阳性。当成本与每个假阴性相关联时，需要更高的召回率。在两者都重要的情况下，f1 分数，即精确度和召回率的调和平均值，变得重要。
参考:
[1]西科拉·m .、弗罗贝尔·l .:规则归纳算法在煤矿地震灾害监测系统收集的数据分析中的应用。采矿科学档案，55(1)，2010，91–114。
感谢您的来访。我希望你喜欢阅读这个博客！
有关使用 SMOTE 处理不平衡类的更多信息，请参考这篇由 CathL 及其团队撰写的博客:
* [## 预测危险的地震颠簸第三部分:改善不平衡数据集的模型性能
新手评估试用机器学习模型进行分类和数据扩充，以更好地支持…
medium.com](https://medium.com/womeninai/predicting-hazardous-seismic-bumps-part-iii-improving-model-performance-for-imbalanced-datasets-88fa64b4d622)*
本笔记本的 GitHub 链接:
* [## royn5618/Medium_Blog_Codes
permalink dissolve GitHub 是超过 5000 万开发人员的家园，他们一起工作来托管和审查代码，管理…
github.com](https://github.com/royn5618/Medium_Blog_Codes/blob/master/Predicting%20Hazardrous%20Seismic%20Bumps/Predicting_Seismic_Bumps.ipynb)*
预测危险的地震碰撞第一部分:EDA，不平衡数据集的特征工程和训练测试分裂
本文展示了探索性数据分析(EDA)、特征工程以及使用 UCI 数据档案中的地震颠簸数据集对不平衡数据进行分割的策略。
简介:
地震冲击数据集是一种鲜为人知的二进制分类数据集，它使用长壁煤矿中的地震和地震声学系统来捕捉地质条件，以评估它们是否容易发生导致地震灾害的岩爆。
【https://archive.ics.uci.edu/ml/datasets/seismic-bumps】链接到数据集:
这是一个很好的数据集，它实际展示了不平衡数据集，处理不同类型的数据分割，并评估分类器的性能指标，包括展示准确性悖论。
关于该数据集的另一件事是，它既有分类又有数字特征，这为争论和尝试不同的特征转换方法提供了一个平台。
这篇文章代码不多，但是在理解什么是对的和错的方面更直观一点！代码可以在我的 GitHub 库中找到。
注意——我没有详细说明 EDA 和功能工程步骤中的每个功能，因为它们是重复的。我只有这个博客和 中的样本，完整的代码可以在 GitHub 的这个链接 中找到。
探索性数据分析
数据事实:
该数据集有 19 列 2584 个实例，其中有 4 个分类特征、8 个离散特征和 6 个数值特征。最后一个是标签列，包含 0 表示无危险，0 表示无危险的地震碰撞。为了便于使用，我将特性名称分类并保存如下:
**col_list_categorical** = ['seismic', 'seismoacoustic', 'shift', 'ghazard']
**col_list_numerical** = ['genergy', 'gpuls', 'gdenergy', 'gdpuls', 'energy', 'maxenergy']
**col_list_discrete** = ['nbumps', 'nbumps2', 'nbumps3', 'nbumps4', 'nbumps5', 'nbumps6', 'nbumps7', 'nbumps89']
**label** = 'class'
属性信息[ 来源:
1 .地震:通过地震方法获得的矿山工作中的移位地震危险评估结果(a-无危险，b-低危险，c-高危险，d-危险状态)；
2。地震声学:通过地震声学方法获得的矿山开采中的移位地震危险性评估的结果；
3。班次:有关班次类型的信息(W —采煤，N—准备班次)；
4。地震能量:监测长壁的
地震检波器中最活跃的地震检波器(GMax)在前一个班次内记录的地震能量；
5。gpuls:GMax 在前一个班次内记录的脉冲数；
6。gdenergy:GMax 在前一个班次中记录的能量与前八个班次中记录的平均能量的偏差；
7。GD puls:GMax 在前一次移位中记录的脉冲数与前八次移位中记录的平均脉冲数的偏差；
8。ghazard:根据仅来自 GMax 的注册，通过地震声学方法获得的矿山作业中的移位地震危险评估结果；
9。nbumps:前一个班次内记录的地震颠簸次数；
10。nbumps2:前一个班次内记录的地震颠簸次数(能量范围[1⁰,1⁰])；
11。nbumps3:前一个班次内记录的地震颠簸次数(能量范围[1⁰,1⁰⁴])；12。nbumps4:前一次移动中记录的地震颠簸次数(能量范围[1⁰⁴,1⁰⁵])；
13。nbumps5:在最后一次移动中记录的地震颠簸次数(能量范围[1⁰⁵,1⁰⁶])；14。nbumps6:前一个班次内记录的地震颠簸次数(能量范围[1⁰⁶,1⁰⁷])；15。nbumps7:前一个班次内记录的地震颠簸次数(能量范围[1⁰⁷,1⁰⁸])；16。nbumps89:前一个班次内记录的地震颠簸次数(能量范围[1⁰⁸,1⁰⁰])；17。能量:前一次移动中记录的地震冲击的总能量；18。最大能量:前一次移动中记录的地震冲击的最大能量；19。类别:决策属性—“1”表示高能地震冲击发生在下一个班次(“危险状态”)，“0”表示在下一个班次(“非危险状态”)没有高能地震冲击发生。
目标类别分布
在 2584 份记录中，只有 6.5%的危险事件。换句话说，50 次地震中有 3 次是危险的。
sns.countplot(x=label, data=df, palette=colors)
plt.xlabel('CLASS')
plt.ylabel('COUNT')
分类特征
首先，我查看了分类特征，以了解它们与标签类别之间是否存在关系偏差。为此，我使用 pandas.crosstab 函数构建了列联表，并评估了分类特征中每个类别的内容相对于类别标签的比率和比例，以确定是否存在偏差。
data_crosstab = pd.crosstab(df['seismoacoustic'], df[label], colnames=['class'])
类别 1 到类别 0 的每个地震声学类别的比率约为 0.06，没有表现出任何强偏差。特征“地震”和“移动”包含一些分布偏差。“地震”特征中的“b”类包含更多的危险地震碰撞，而“shift”特征中的“W”类包含的地震碰撞比“N”类多。这些分类特征的列联表如下:
左侧——针对“班次”与“等级”的列联表；中心——地震与等级的应变表；右侧——针对“ghazard”与“Class”的列联表
数字特征
接下来，我查看数据并挑选出数字列，以理解它们的描述性统计、相关性和分布图。
- 描述性统计
在 Pandas 中， 描述 方法仅提供计数、平均值、最小值、最大值、标准值和百分位数，我们可以从这些值中评估数据在某种程度上的偏斜度。为了更好地理解它们的分布，视觉效果是必不可少的，如分布图和直方图。
df[col_list_numerical].describe()
- 相关性&热图
我用熊猫评估了这些特征之间的相关性。DataFrame.corr 并使用 seaborn.heatmap 可视化。
df_corr = df[col_list_numerical].corr()
plt.figure(figsize=[8, 8])
sns.heatmap(data=df_corr, vmin=-1, vmax=1, cmap=’gist_earth_r’, annot=True, square=True, linewidths=1)
plt.xticks(rotation=90)
plt.yticks(rotation=0)
- 来自描述性统计和相关性的观察:
能量和最大能量的分布看起来很相似。genergy 看起来右倾。为了更好地理解分布，我使用了 seaborn.distplot 来可视化每个特性，如下所示。
在热图中，很明显 genergy 和 gpuls、gdenergy 和 gdpuls 之间有很强的相关性。能量和最大能量几乎是完全相关的。在 EDA 的后面部分，我试图为这些高度相关的对构建散点图，并可视化它们在每个目标类中的分布。
- 更多 EDA 可视化:分布图&散点图
就像我说的，我想将分布可视化如下，很明显 genergy 是右偏的，就像所有其他数字特性一样(如 Github Notebook 所示)。
sns.distplot(df['genergy'], hist=True)
genergy 的更多实例在 0 到 200000 的范围内，并且分布随后向更高的能量值逐渐减少。所有其他数值特征的模式相同，但 gdenergy 和 gdpuls 的范围小于其他数值特征。
此外，我使用散点图来可视化相关系数大于 0.70 的相关特征。我构建了散点图来查看这些数值是如何分布在每个目标类中的。
plt.figure(figsize=[10, 8])
sns.scatterplot(x='genergy', y='gpuls', hue='class', data=df)
左侧 genergy 与 gpuls 的散点图；右侧—能量与最大能量的散点图
在左图中，地震凸起更集中在较高的 gpuls 和 genergy 值处。右侧散点图中观察到的能量与最大能量的线性关系也很强，如描述性统计数据、它们的散点图和此散点图所示。
离散特征
“nbumpsX”是包含从 0 到 9 的整数的离散特征。这里可以看到，‘nbumps 6’全是零，与‘nbumps 7’和‘nbumps 89’相同。此外，它们都是右偏或正偏分布。
我使用下面的代码删除了“nbumps6”、“nbumps7”和“nbumps89”。其他一切看起来都很好，但是任何 ML 模型都很有可能将‘nbumpsX’信息视为序数。
df.drop(columns=['nbumps6', 'nbumps7', 'nbumps89'], inplace=True)
我还使用以下代码再次构建了列联表，以查看 nbumpsX 如何与目标类联合分布:
for each_col in col_list_discrete:
data_crosstab = pd.crosstab(df[each_col], df[label], colnames=['class'])
print(data_crosstab)
print('-----'
在表中，危险和非危险情况的分布没有明显的差异。在所有这些分立特性中，它们在 nbumpsX 计数中按比例出现。
因为所有这些看起来都很好，所以我转到了下一节“特性工程”，在这里我处理了分类和数字特性。
特征工程
我一次性编码了分类变量，并转换了数值列以缩小范围。我保持了独立的特征不变。
- 分类特征的一键编码
我将所有两个类别的分类特征编码成一个包含 0 和 1 的特征，同时将两个以上类别的特征转换成每个类别的二进制列。这个 one-hot encoded feature 列表示该实例是否包含该特定类别。
例如，地震声学特征包含三个类别——a、b 和 c(如左图所示),它们被转换为扩展的三个二元特征(如右图所示)。我使用 drop='first '来避免伪变量陷阱。
label_encoder = LabelEncoder()
onehot_encoder = OneHotEncoder(drop='first', sparse=False)
encoded_array = label_encoder.fit_transform(df[col_name])
encoded_array_reshaped = encoded_array.reshape(len(encoded_array),1)
one_hot_encoded_array = onehot_encoder.fit_transform(encoded_array_reshaped)
num_features = one_hot_encoded_array.shape[1]
new_enc_col_names = [col + '_enc_' + str(num) for num in range(0, num_features)]
df_enc = pd.DataFrame(one_hot_encoded_array)
df_enc.columns = new_enc_col_names
df = pd.concat([df, df_enc], axis=1)
df.drop(columns=col, inplace=True)
- 数字列的转换
然后，我使用对数变换来转换所有的数字特征。z 分数变换并比较原始分布的变换。分布现在看起来更“正常”，范围也更好了。参考下面的代码和可视化进行比较。
sns.distplot(**np.log**(df['genergy']), hist=True)
# Results in the figure on the right
左 genergy 的实际分布；右——遗传的对数变换分布
我基本上创建了一个字典来定义我希望在这些列上应用什么转换，然后循环遍历字典条目，如下所示:
def shifted_log_func(df_col):
return np.log(1 + df_col)dict_num_cols_trnsfm = {'genergy': np.log,
'gpuls' : np.log,
'gdenergy': stats.zscore,
'gdenergy': stats.zscore,
'energy': shifted_log_func}for col_names, transfm_func in dict_num_cols_trnsfm.items():
df['scaled_' + col_names] = transfm_func(df[col_names])
df.drop(columns=col_list_numerical, inplace=True)
df[[col for col in df.columns if 'scaled_' in col]].describe()
由于两者的数据相似，删除 maxenergy 后，数值的最终描述性统计数据。
左侧-缩放数字要素的描述性统计数据；右-原始数值数据的描述性统计信息
左侧表格中的范围似乎不错。右边的表格包含一些能量和最大能量的零值。为了解决这个问题，我在“能量”栏中加了 1(因为我去掉了 maxenergy ),然后应用对数变换。
分割培训和测试数据
由于这是一个不平衡的数据集，我使用分层混洗分割来分割和按比例分配目标类。
stratified_split = StratifiedShuffleSplit(n_splits=1, test_size=0.20)for train_idx, test_idx in stratified_split.split(X, y):
y_train= y[train_idx]
X_train = X[train_idx]
X_test, y_test = X[test_idx], y[test_idx]print("Training Set Target Class Distribution:")
print(y_train.value_counts()/len(y_train))
print("Test Set Target Class Distribution:")
print(y_test.value_counts()/len(y_test))
如果我没有使用这种技术，我可能会以稍微不平衡的比例结束危险的地震颠簸，如下图所示。机器学习算法将能够从更公平比例的类别标签中学习，而不是从随机分布中学习，在随机分布中，目标类别的较小类别可能太少而无法学习，因此降低了模型的性能。
(左)-使用 train_test_split with shuffle 进行拆分时的分布= False(右)-使用 train_test_split with shuffle 进行拆分时的分布=True
这些微小差异的影响可能会影响机器学习算法的性能，因此对于二进制和多类分类问题都是如此。因此，建议对分类问题使用分层分割，尤其是当存在不平衡数据集时，以保留每个类别的样本百分比并避免采样偏差。
感谢您的来访。我希望你喜欢阅读这篇博客。 下一篇文章 是这篇探索模型构建和性能指标比较的博客的继续。下面链接-
[## 预测危险的地震颠簸第二部分:训练和调整监督 ML 分类器和模型…
本文展示了预测危险的地震颠簸使用二元分类器，调整模型超参数…
towardsdatascience.com](/predicting-hazardous-seismic-bumps-part-ii-training-supervised-classifier-models-and-8b9104b611b0) [## 预测危险的地震颠簸第三部分:改善不平衡数据集的模型性能
新手评估试用机器学习模型进行分类和数据扩充，以更好地支持…
medium.com](https://medium.com/womeninai/predicting-hazardous-seismic-bumps-part-iii-improving-model-performance-for-imbalanced-datasets-88fa64b4d622)
本笔记本的 GitHub 链接:
[## royn5618/Medium_Blog_Codes
permalink dissolve GitHub 是超过 5000 万开发人员的家园，他们一起工作来托管和审查代码，管理…
github.com](https://github.com/royn5618/Medium_Blog_Codes/blob/master/Predicting%20Hazardrous%20Seismic%20Bumps/Predicting_Seismic_Bumps.ipynb)
基于神经网络的 HDB 房价预测
实践教程
除了问你的房屋中介，你有没有想过你的 HDB 值多少钱？或者想知道对你的 HDB 有很大影响的主要特征是什么？看到超过 92.2%的新加坡人拥有 HDB，我相信这些是大多数新加坡人都有过的常见问题，如果我能使用数据分析来获得对新加坡政府住房市场的任何见解，那将是很好的。序列预测模型给出的平均绝对误差为 22，429 美元(+/- 571.83)
关于可视化的更多信息，请查看我的 tableau 公共简介。有关代码的更多信息，请查看我的 GitHub。所有数据均在data.gov.sg获得。
概观
data.gov 的所有数据都有。以下数据分析项目遵循典型的方法(CRISP-DM)。
- 数据清洗+探索性数据分析
- 顺序神经网络(3 层)
- 估价
- 结论
- 未来的工作
在接下来的 EDA 阶段，我将致力于数据可视化方面的深入研究，为我的读者提供对新加坡 HDB 房地产市场的更深入的见解。
数据清理
该数据集包括从 1990 年到 2020 年的数据。总共有 677，114 行和 10 列。虽然数据集中没有空值，但也有一些异常值(例如，封顶与未封顶的要素等。).首先，由于我们考虑的是在 MOP(最短入住期)后出售的公寓，任何超过 95 年租约(剩余租约)的公寓都将被删除。
然而，在数据集中有大量的重复，它们已被删除。
目标功能概述
我们整个数据集的中值价格约为 35 万美元。另外我们可以看到，80 万以上的房子只有少数。有趣的事实:根据我们的数据，这些年来卖出的最贵的房子是 123.2 万英镑。
图一。房价分布
探索性数据分析
相关矩阵
图二。所有定量特征的相关矩阵
公寓类型和建筑面积
正如我们所知，在新加坡有多种类型的公寓，一室，两室，行政等。直觉上，我们认为更大的房子(更大的建筑面积)是价格的有力指标。我把两张图拼在一起，这样你就可以看到公寓类型/建筑面积和转售价格之间的强烈相关性。正如我们所想象的，房子越大，价格就越高。
图 3。建筑面积(平方米)与转售价格
图 4。建筑面积与转售价格中位数
越大=越好？简短的回答；是的。
故事范围
在数据集中，地板水平被给定为一个范围并被注册为一个字符串。(例如“01 至 03”、“03 至 05”等。).此外，还有一些从 01 到 05 的楼层范围，跨越两个特征。我根据这些特征的平均楼层对它们进行了编码，并相应地将它们从 1(最低楼层范围)到 17(最高楼层范围)进行了排序。
图 5。转售价格与楼层范围
根据转售价格对它们进行排序，我们注意到，与其他功能相比，故事范围从 01 到 05 具有更高的中值价格。而总的趋势显示，地板越高，转售价格越高。
同样，对于故事范围 15、16 和 17，它们转换为 43 到 45、46 到 48 和 48 到 51，差异是< $200,000 (1 Sigma). There aren’t too many HDB in Singapore that has floor that are higher than 43 floors, deeper analysis shows that most apartments sold there are from here.
Flat Model/Flat Layout
Flat model essentially refers to the layout of the house. As HDB have a long history in Singapore housing market. some of these models are obsolete, as the land are getting scarcer, some flat models will never return (e.g. Masionette).
Figure 6. Resale Price VS Model Type, ranked according to mean
Date (Month & Year)
The date was split into month and years as I firmly believe that the months do not matter as much as compared to the years. Perhaps changing them into Q1Y19 would give a deeper resolution into how HDB resale price changes overtime.
月
显然，这几个月没有太大的区别。月份从最小的中间转售价格开始排序。
图 7。转售价格中位数与月数
图 8。几个月内售出的房产总数
我把几个月作为一个特征，因为它们不会对转售价格的变化产生很大影响。有趣的是，计数图显示，多年来，二月一直是不受欢迎的购房月份。我认为这是因为二月的天数最少。
年份(HDB 转售指数)
通过对历年(1990 年至 2020 年)转售价格的排名，我们可以看到 HDB 转售价格多年来有明显的趋势。
图 9。转售价格与年销售量
HDB 转售指数随时间波动，在 2013 年达到峰值，HDB 的平均售价约为 40 万美元。从 1997 年到 1998 年，一年中最大的年同比降幅约为 50，000 美元。
塑造转售指数的一些值得注意的事件是 1997 年和 2008 年房地产泡沫期间发生的亚洲金融危机。然而，我们看到，HDB 转售指数在 2008 年房地产危机期间并没有受到太大影响。相反，由于我不知道的原因，它在 2013 年至 2014 年期间出现了更严重的下降。我最初的假设是 SIBOR 利率上升，或者在 2013 年和 2014 年有更好的资产类别可以投资
镇(区)
我认为，我认识的大多数新加坡人对财产有一种近乎不健康的痴迷(包括我在内)。新加坡人几乎愿意花大价钱住在黄金地段/成熟地区，其中一个原因可能是因为连通性和良好的小学。不用进一步滚动，你大概就能猜出哪些区域会卖高价。
图 10。转售价格与城镇
令人惊讶的是，远离市中心的丰戈尔获得了最高的转售价格，然而，从数据来看，丰戈尔是最近 10 年才开发的，这意味着平均零售价格已经很高了。
另一方面，虽然璧山、中央区和武吉提马的平均价格较低，但与其他地区相比，它们的 75%销售价格都较高。
对于成熟的房地产，如皇后镇，芽笼，昂莫桥和大巴窑。虽然它们靠近中心区，被认为是成熟的住宅区，但它们的收入中位数较低。这可能是因为这些地区的房屋一般都比较旧，可能已经过了最佳剩余租期。
图 11。剩余租约与城镇
剩余租约
正如我们可以看到的密度图显示，大多数人在两点出售他们的房子。94 年(紧接着 5 年的 MOP)2。70-85 岁(又名升级者/降级者)。这很直观，因为新加坡的组屋租期为 99 年，最短入住期为 5 年。因此，一旦最低入住期结束，人们可能会升级。第二个最常见的销售期在 70-85 左右是因为房地产价格大幅上涨。
图 12。转售价格与剩余租赁的热图。
小学
对于大多数新加坡人来说，我们听说过父母搬家只是为了更接近一所更好的(主观的)小学。这是因为在投票时，居住在更靠近首选学校的新加坡公民和永久居民会得到优先考虑*。*
学校首先根据其受欢迎程度排名。最后，学校被地理标记到各自的城镇(区),然后进行平均。形成新的特色(学校排名)。
如你所见，武吉提马、璧山和中部地区的平均排名非常高。林珠港没有学校，事实上，快速搜索谷歌会发现林珠港几乎没有任何组屋。
序列神经网络
部署了三层神经网络，因为当与 SGDregressor 和 XGBoost 比较时，它产生最佳结果。基于 MAE 对结果进行评估。
我们利用了一个简单的 3 层序列神经网络，所有特征都被归一化。使用 minmax scaler 将分类特征拆分到单独的列中。
重复 K 倍
使用重复的 K-fold 来评估模型的性能。重复的 K-fold 随机地将数据分成大小大致相等的 5 个块。每个块依次被排除，而另外 4 个块用于训练模型。预测保持的块，并且这些预测被总结为 MAE(平均绝对误差)和 RMSLE(均方根误差)。然后重复该过程 3 次，并对各自的性能指标进行平均。更多关于重复 k 倍的信息，我们可以参考这个。
图 14。k 倍交叉验证，混洗数据集，重复该过程 3 次。
结果
图 15。历元数与 MAE(左)和 MSE(右)
随着时代数量的增加，我们的损失(MAE，MSE，MAPE)减少，除了 MAPE，我相信这是因为它是原始房价的函数。我相信我们还没有找到“最佳点”,因为我们看到，随着我们的训练和验证数据中的时期增加，我们的损失仍在继续减少。
评估指标
重复 k 折叠序列神经网络模型结果是，
平均汇率:22767.46 (+/- 571.83)
MAPE: 9.07% (+/- 0.35%)
均方根误差:0.017(+/- 0.006)
结论
总的来说，模型的表现比我预期的要好。我试图寻找 HDB 预测模型的基准，它徘徊在 0.09 到 0.03 RMSLE 之间。然而，神经网络缺乏可解释性。与线性回归模型相比，我们可以看到不同系数的具体权重。
沙普利值
由于神经网络在很大程度上是黑盒模型，利用 Shapley 值来解释哪些特征在确定价格方面贡献最大。Shapley 值显示了哪些特征对价格预测贡献最大。
图 16。平均 SHAP 值
特征年销售量对价格预测贡献最大，其次是户型(3 房和 4 房)和建筑面积平方米。这意味着销售年份在很大程度上影响了预测，其次是公寓类型是三房还是四房。
未来的工作
最初的数据分析表明，确实有某些特征有助于更高的转售价格。然而，我相信我们可以更多地使用街区号码和街道名称。如果我们能够地理定位每个地址的纬度和经度，我们就能够创建新的距离特征。
此外，我还将尝试不同的算法，为模型提供更好的解释能力。目前，我们无法看到不同的功能如何影响最终的预测价格。
预测新冠肺炎疫情所需的医院容量
预测疾病传播和住院需求，确定应对激增所需的床位数量。
在我的上一篇文章中，我使用时间序列预测了确诊病例的指数增长&应用了修正以考虑社会距离和其他类似的措施。在这篇文章中，我将使用一个简单的流行病学模型来关注传染病传播的建模。在文章的开始，我将简要介绍该模型如何工作的基础知识，并使用该模型预测传播和住院需求
请注意，本文仅用于数据分析目的，本文观点不应被解释为专业建议。
流行病模型是基于分室的模型，将人群分成不同的组。目标是模拟这些疾病如何从一个群体的一个成员传播到另一个群体。一个最简单的房室模型叫做 SIR ，它将种群分为 3 组:
- 可能被感染的易感人群
- 受感染的个人(一)
- 康复个体®。
请注意，术语“康复”指的是可能已经痊愈并产生免疫力的患者，或者可能已经去世的患者。他们不会再被感染，这在现实中还没有被证实。
SIR 模型基于当前观察到的局部参数，如患病率和倍增时间，提供了对未来疾病在人群中传播的洞察。请注意，传播对隔离政策、检测可用性和其他当地因素非常敏感。这个模型的目标是理解
- 疾病传播的速度有多快
- 有百分之几的人口会被感染？
- 有百分之多少的人康复了？
- 我们的医院容量应该是多少？
SIR 模型预测将被感染的潜在个体总数。该数字可用于进一步估计需要医院资源(如住院、ICU、通气等)的患者比例。请注意，该模型对所有人都一视同仁，而实际上，患有更多慢性病的老年人群似乎风险更高。
SIR 模型使用初始(当前)疾病流行参数预测未来感染，如
- 群体大小
- 恢复时间
- 传输速率
方案
我们举个简单的例子。假设我们有一个拥有 100 万人口的县/州，在那个镇上，有五个人被感染。目标是确定它能传播多少和多快。这取决于许多因素。
- 我们是否在跟踪家中的避难所及其对感染者每天可感染人数的影响:每个感染者每天可联系或多或少的人，这取决于是否实行社交距离。假设每个人每天接触两个人。不是每个他/她接触的人都会被感染。存在传播的可能性(易感者和受感染者接触后的感染可能性)。比方说传输概率是 15%。
- 这意味着，如果这个人每天遇到 6 个人，传播概率为 15%，他每天会感染 6 * 15 % = 0.9 人。这个数字被称为 beta (𝛃)又名传播率，代表一个感染者每天感染的人数。出于建模的目的，每天接触的人数可能低至 1 人，高至 10 或 15 人，传播概率在 5%至 15%之间。较高的传播概率可能与具有较高风险的人群(患有其他慢性疾病的老年人群)相关
- **疾病可以从感染者身上传播的天数:**我们已经知道这个感染者每天可以传染 0.9 人。这个人多少天会被感染，传播疾病？假设这个病人与人接触了 5 天。如果他/她每天感染 0.9 人，5 天就感染 0.9 * 5 = 4.5 人。这个数字被称为基本繁殖数，用来估计疾病在人群中传播的速度。这是一个感染者感染的总人数。这是一个需要理解的重要数字，因此让我们在下面对此进行更多的讨论。
- 再现号 Ro 读作“R 零”这是一个数学术语，表示疾病的传染性。告诉你平均有多少人会从一个有传染性的人身上染上疾病。如果一种疾病的 Ro 为 6，只要没有人接种过疫苗或在他们的社区中已经对其免疫，感染者就会将该疾病传播给平均 6 个其他人。2009 年的猪流感或 H1N1 病毒的 Ro 值约为 1.5。由于疫苗和抗病毒药物的作用，影响并不大。在冠状病毒的情况下， R o 值估计在 1.5 到 3.5 之间，其中 1.5 到 2.5 用于良好的社交距离。
- **痊愈率(ɣ)**是单位时间内痊愈的感染者的比率。在这种情况下，感染持续了 5 天，⅕是伽马值，是目前每天恢复的感染人数。
这些参数可用于通过识别感染和康复人群的变化来定义传播。
随着时间的推移，易感人群的数量下降，被感染的人数上升。康复的人数也在上升。(𝛃)传染率是指一个人感染多少人。既然我们有我感染的人，那就是𝛃 我就是感染的力量。当我们将其乘以易感人群，即新感染人群。由于一些种群可能已经恢复，ɣ I 给出了恢复的种群数量。必须减去这个数字才能得到感染人数。这可以用一个简单的等式来描述，如下所示。我不打算在本文中解释这个公式的数学原理，但会解释它的实际用法。
使用 SIR 模型进行预测
预测一个州/县所需的医院床位
现在，我们将上述知识应用于拥有 100 万人口的县/州(如前所述)。假设这个县/州在 3 月 10 日左右宣布了 5 个病例。
这里要解决的问题是
- 疾病传播的速度(到达高峰的时间)
- 有百分之几的人口会被感染？
- 社交距离的影响？
- 对我们医院容量的影响？
截至 2010 年 4 月 20 日，美国主要一般统计数据的简要总结。这些是我们可以在计算中使用/参考的一些公开可用的统计数据。这些数据来源于世卫组织和疾病预防控制中心。
- 被测试人口的百分比相当低(< 1%).
- Out of the tests conducted, on an average anywhere from 10 to 20% of the people tested are becoming positive.
- Out of the individuals who have confirmed positive, 10 to 16% of those are hospitalized.
- Most states in the USA the number of beds range from 1.6 to 4.4 beds for every 1,000 people. The average for the USA is 2.4.
用于模拟的关键数字
- 假设一个人感染了 5 天。如前所述，这意味着γ(回收率)是⅕ = 0.2。为了限制建模的范围，我保持伽玛为常数。在后续文章中，我可能会进一步阐述这种影响。
- 从本文前面的讨论来看，对于新冠肺炎来说基本的再生数 R o 的范围是从 1.25 到 3.5。我们将使用 1.25 的 Ro ，1.5 的& 2 并比较结果。告诉你平均有多少人会从一个有传染性的人身上染上疾病。例如， R o 为 2 表示一个感染者会将疾病传染给另外两个人。如果 R o 值大于 1，则感染率大于恢复率，因此感染将在整个人群中增长。如果 Ro 小于 1，感染将很快消失，因为人们愈合的速度比传播的速度快。
- 相同的β值将分别为( R o * gamma) 0.25，0.3 & 0.4。β代表一个感染者每天感染的人数。
- 对于 SIR 模型，截至 3 月 10 日，最初的易感人群约为 100 万，其中有 5 人被感染。
下一节总结了这个模拟的结果。
Ro 1.5 的结果汇总
对于 Ro = 1.5，从下图可以看出，峰值出现在 3 月 10 日之后的 112 天。这将是 7 月初左右，高峰病例为 63，000 例。注意 y 轴表示人数(%)，x 轴表示天数。
比方说，我们县每个病人有 2.1 张病床& 16%的检测呈阳性的人住院了(根据上面的统计)。让我们用那个计算来估计医院的容量。
从 3 月 10 日起的 112 天中，峰值病例为 63，000，相应的峰值住院需求为 10，080。
添加可用的医院容量，见下图。在现有床位为 2，100 张的情况下，要为 10，080 张床位的激增做好准备，需要通过增加医院容量以及练习社交距离来减少激增。在下一节中，我们将讨论将 Ro 从 1.5 降至 1.25 的影响。
社交距离的影响或缺乏社交距离
下图显示了范围从 1.25 到 2.0 的各种 Ro 值的影响。y 轴显示人数(%)，x 轴显示天数。如果一个人做了很好的社交距离，我们可以看到 Ro = 1.25 时的行为。如果我们开始放松社交距离，你会看到更接近于 2.0(或者更高)的行为。这就是为什么社交距离、测试、追踪和隔离极其重要。
对于 Ro = 1.25，假设继续做好社交距离等隔离方式，197 天左右可以看到 21500 例。这是从 3 月 10 日开始的大约 7 个月，也就是 9 月底/10 月初。
让我们加上每 1000 个可用床位 2.1 张的计算结果(来自统计数据)& 16%(来自统计数据)测试呈阳性的人住院了。考虑到这些因素，所需床位计算如下。
相应的高峰住院需求为 3440 人，而容量仅为 2100 人。
预测不同 Ro 值的医院床位容量
下表总结了 3 个不同 Ro 值的建模结果。峰值表示自模拟开始(2020 年 3 月 10 日)以来可以找到的最大案例数，天数为天。
让我们加上每 1000 个可用床位 2.1 张的计算结果(来自统计数据)& 16%(来自统计数据)测试呈阳性的人住院了。考虑到这些因素，所需床位计算如下。
对于 Ro = 1.25，该模型预测高峰病例将在 197 天后发生，大约高峰日期为 9 月 23 日，需要医院床位 3440 张。假设每 1000 人 2.1 张床位，该县将短缺 1340 张床位。
类似地，对于 Ro = 1.5，该模型预测高峰病例将在 112 天后发生，大约高峰日期为 6 月 30 日，所需医院床位为 10，080 张。假设每 1000 人拥有 2.1 张床位，该县将缺少 7980 张床位…
如果我们假设需要住院治疗的检测呈阳性的人口百分比为保守估计的 10%,则以下是相应的所需医院床位容量。
对于 Ro = 1.25，该模型预测该县不会出现短缺。就 Ro 1.5 而言，预计短缺 4，200 个床位。
通过适当的管理，将 Ro 保持在 1.25 以下，并采取适当的谨慎措施，人们可以处理这一危机。
总结
使用 SIR 模型，我们预测了一个拥有 100 万人口的县的传播范围，该县截至 3 月 10 日有 5 例确诊病例，具有不同的 Ro 值。如果该县能够将 Ro 保持在 1.25 或以下，我们应该能够很好地控制病毒爆发。如果 Ro 达到峰值，那就是我们开始不知所措的时候。这表明了社会距离、接触者追踪和隔离的重要性。更重要的是，知道大概的高峰时间和潜在的住院需求，县可以通过增加适当的临时床位来帮助减少影响。
扫一扫
728
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
