TowardsDataScience 博客中文翻译 2022（一百一十四）

admin管理员组

文章数量:1645531

原文：TowardsDataScience

协议：CC BY-NC-SA 4.0

使用 Python 创建漂亮的独立交互式 D3 图表

原文：https://towardsdatascience/creating-beautiful-stand-alone-interactive-d3-charts-with-python-804117cb95a7

应用于 D3 力有向网络图

截图来自:【https://d3js/

可视化数据可能是项目成功的关键，因为它可以揭示数据中隐藏的见解，并增进理解。说服人们的最好方法是让他们看到自己的数据并与之互动。尽管 Python 中提供了许多可视化包，但制作漂亮的独立交互式图表并不总是那么简单，这些图表也可以在您自己的机器之外工作。D3 的主要优势在于它符合 web 标准，因此除了浏览器之外，你不需要任何其他技术来使用 D3。重要的是，互动图表不仅有助于告诉读者一些东西，还能让读者看到、参与和提问。在这篇博客中，我将概述如何使用 Python 构建你自己的独立的、交互式的力导向 D3 网络。 注意步骤与任何其他 D3 图表相似。如果你需要工作版本， d3graph 库就是给你的！

如果你觉得这篇文章很有帮助，可以使用我的 推荐链接 继续无限制学习，并注册成为中级会员。另外， 关注我 关注我的最新内容！

使用 D3 的动机

D3 是*数据驱动文档的缩写，*是一个 JavaScript 库，用于在 web 浏览器中生成动态、交互式数据可视化。它利用了可缩放矢量图形(SVG)、HTML5 和级联样式表(CSS)标准。 D3 又称 D3.js 或 d3js。 我会互换使用这些名字。在交互式 Python 解决方案上面使用 D3 的优势有很多，我来列举几个[1]；

• D3 是一个 Javascript 库。因此，它可以用于您选择的任何 JS 框架，如 Angular.js、React.js 或 Ember.js。
• D3 专注于数据，所以如果你在数据科学领域，数万个数据点保持交互而不需要付出任何额外的努力，这是一种解脱。
• D3 是轻量级的，直接与 web 标准一起工作，因此速度非常快，并且可以很好地处理大型数据集。
• D3 是开源的。因此您可以使用源代码并添加自己的特性。
• D3 支持网络标准，所以除了浏览器，你不需要任何其他技术或插件来使用 D3。
• D3 支持 HTML、CSS 和 SVG 等 web 标准，在 D3 上工作不需要新的学习或调试工具。
• D3 没有提供任何特定的特性，所以它可以让你完全控制你的可视化，以你想要的方式定制它。

D3 图表。

D3 是被设计成一起工作的模块的集合；您可以单独使用这些模块，也可以将它们作为默认构建的一部分一起使用。D3 网站提供了 168 个工作图表 ，允许交互过程中的性能增量更新，支持拖动、刷、缩放等流行交互。图表可用于多种目的，例如*定量分析、可视化层次结构、创建网络图、以及条形图、线图、散点图、辐射图、地理投影、*和各种其他交互式可视化，用于探索性解释。一些精选的图表如图 1 所示。这些 d3js 图表的各种各样在 D3Blocks 库 *中很容易被 python 化。*阅读 D3Blocks 中型文章了解更多详情。

图 1:D3 图表的各种例子。截图来自https://d3js/

入门！

让我们开始我们的 D3-Python 项目吧！在我们迈向 python 化 D3 图表之前，我们需要理解 D3 是如何工作的。我将通过在 D3 中创建一个非常小的力导向网络图(没有 Python)来演示它，之后，我将演示 Python 中的集成。解释 D3 图表发展的最好方法是将发展分成四个不同的部分；

1. 层叠样式表(CSS)。
2. D3 模块包含所有需要的库。
3. Javascript 构建图表。
4. 数据为 JSON 文件。

图表中的每个部分都有自己的角色，为了得到一个工作版本，我们需要连接所有部分，这可以在一个 HTML 文件中完成，如图 2 所示。在接下来的部分中，我将描述每个部分的作用和实现。

图 2:在一个最终的 HTML 文件中，四个部分被连接起来以构建 D3 图。图片来自作者。

1.级联样式表(CSS)

CSS 文件是一种简单的机制，用于以结构化的方式向 HTML 文件添加样式。例如，CSS 文件可以定义 HTML 元素的大小、颜色、字体、行距、缩进、边框和位置。我们将创建强制定向网络图，为此我们将定义整体字体大小、字体系列、样式、颜色以及特定于节点和边的属性。您可以相应地更改 CSS，并设置您认为最好的属性。我创建的文件可以在这里下载，看起来如下:

级联样式表(CSS)

2.D3 模块

D3 模块是最简单的部分，因为您只需要导入它或将整个内容嵌入到最终的 HTML 文件中。D3 模块包含了创建任何图表的所有函数。不需要对该文件进行编辑。最新版本是 v7 ，可以从本地资源或网站导入。

<script src="https://d3js/d3.v7.min.js"></script>

我将使用稍微老一点的版本( v3 )，因为我已经用这个版本创建了很多脚本。可以在这里下载，它包含在我下面几个步骤的示例代码中。如果您构建自己的图形，我强烈建议使用最新版本。

3.构建图表的 Javascript。

构建图表不需要您从头开始编码。我们可以在 D3 网站上找到感兴趣的图表，并使用开源代码。让我们转到强制导向图页面并阅读“ 我如何使用这段代码？ 。描述如下:

图 3:截图来自https://observablehq/@d3/force-directed-graph。

它告诉我们复制粘贴整个" ForceGraph "函数，如下图所示。我创建了另一个截图(图 4)，你可以看到“函数 Forcegraph({ )的一部分，这就是 ForceGraph 函数。您需要复制整个函数并将其粘贴到一个新文件中。我将我的文件命名为D3 graph script . js，可以在这里下载。请注意，我的函数与这里显示的最新函数略有不同。该网站还告诉我们导入 D3 模块，我们在前面的部分中已经介绍过了。

图 4:部分 ForceGraph 代码。截图摘自:https://observablehq/@d3/force-directed-graph。

4.数据

图表的燃料是我们需要得到正确形状的数据。尽管 D3 有导入本地文件的功能。使用d3.json()的 json-file，它可能无法工作，因为导入带有 D3 的本地 csv 或 json 文件被认为不安全。一种解决方案是将数据直接嵌入到最终的 HTML 文件中。但是，根据数据量的不同，它可能会产生大量的 HTML 文件。然而，将数据直接嵌入到 HTML 中会将所有脚本和数据集成到一个文件中，这非常实用。下面的方框显示了一个 json 数据文件的例子，其中描述了一些节点和边。我们可以在 HTML 中嵌入这样的数据块。完整的 json 数据文件可以从这里下载。

graph = {"links":[{"node_weight":5,"edge_weight":5,"edge_width":20,"source_label":"node_A","target_label":"node_F","source":0,"target":2},{"node_weight":3,"edge_weight":3,"edge_width":12,"source_label":"node_A","target_label":"node_M","source":0,"target":4}"]}

5.连接零件

我们现在有了构建图表所需的四个部分( CSS、D3、javascript 图形和数据)。我创建了一个 HTML 样板文件，其中的四个部分将被连接起来；下载了这里。这些零件包括如下:

1. 层叠样式表为Style . CSS
2. D3 模块为D3 . v3 . js
3. 将 javascript 构建为D3 graphscript . js
4. 数据为JSON _ Data

要创建最终的 HTML，您可以替换每个文件中包含它的行上的内容。或者换句话说，用style . CSS中的内容替换 {% include “style.css” %} 以此类推。你会得到一张看起来像这样的工作图。检查页面源代码或下面的脚本，查看四个部分的组合。就是这样！你在 D3 中创建了一个有向力网络图！因为我们创建了四个构建模块，所以在 Python 中集成这个图是一小步。

制作力定向图的最终 HTML。在这里下载，粘贴到一个纯文本文件中，去掉标签<！- TEXT - >，并重命名(如 forcedirected.html)。双击它。

把 D3 变成 Python

Pythonize 化 D3 脚本的步骤是将静态值改为变量。我们可以为最终 HTML 文件和强制定向 javascript 文件中的属性执行此操作。想法如下:

1. 使用 Python 导入/读取最终的 HTML 和/或 javascript 文件。
2. 使用 Python 将变量替换为任何所需的值。
3. 将最终调整后的 HTML 文件写入磁盘。

但是首先，我们需要在最终的 HTML 中将静态值手动更改为变量名。让我们把的宽度、高度、电荷、距离、方向和碰撞转换成变量。创建唯一的变量名以避免意外的错误替换是很重要的。我的解决方案是创建变量名，如 {{ width }} ，现在可以很容易地在 Python 中找到并替换为真实值。除了最终 HTML 文件和 javascript 图形文件之外， json 数据文件也包含可以相应更改的变量。注意 D3 模块将保持不变。在开始实现所有这些之前，请阅读下一节！

值被转换成变量的 HTML 文件。

D3 图形库

d3graph 库的设计方式与以上章节所述类似。这是一个 Python 库，构建于 D3 之上，创建了一个独立的、交互式的力导向网络图。输入数据是邻接矩阵，其中列和索引是节点，值为 1 或更大的元素被视为边。输出是一个包含交互式力定向图的 HTML 文件。这里的是泰坦尼克号使用案例中的网络演示。我的 D3 版本有一些扩展，其中一个滑块可以根据边值断开网络的边，双击一个节点会高亮显示该节点及其连接的边。要制作自己的网络图， pip 安装 d3graph 包，其中包含:

pip install d3graph

现在，您可以使用各种参数创建交互式网络:

图 5:使用 d3graph 创建的各种图表。图片由作者提供。

最后的话

恭喜你！你刚刚学习了如何创建漂亮的 D3 图表，特别是使用 D3 的力定向网络，以及如何将它集成到 Python 中。我希望这个博客能给你创建任何你想要的 D3 图表所需要的知识，不管有没有 Python。D3 图形库 将帮助你使用 Python 创建自己的 D3 力定向图。输出是一个单独的 HTML 文件，可以独立工作，可以共享或发布在网站上，除了浏览器，你不需要任何其他技术。您还可以将 D3 网络图作为可视化集成到您现有的库中。阅读这篇关于我如何为 HNet 库 学习联想的博客。如果您想创建其他 d3js 图表，如散点图、小提琴图、移动气泡图、桑基图、热图、和弦、时间序列、Imageslider 或粒子图，我们创建了 D3Blocks 库 ，其中包含 10 个可以用 Python 创建的漂亮图表。阅读 D3Blocks 中帖【2】了解更多详情。随意摆弄库，叉源代码。

注意安全。保持冷静。

干杯，E.

如果你觉得这篇文章很有帮助，可以使用我的 推荐链接 继续无限制学习，并注册成为中级会员。另外， 关注我 关注我的最新内容！

软件

• D3 图形库
• HNet 媒体博客
• HNet 库
• D3Blocks 库

我们连线吧！

• 让我们在 LinkedIn 上联系
• 在 Github 上关注我
• 跟随我在媒体上

参考

1. https://www.tutorialsteacher/d3js
2. E.Taskesen， D3Blocks:创建交互式和独立 D3js 图表的 Python 库 ， Medium，迈向数据科学，2022 年 9 月

用 Python 制作漂亮的地形图

原文：https://towardsdatascience/creating-beautiful-topography-maps-with-python-efced5507aa3

实践教程

用 Python 制作漂亮的地形图

当您可以用 Python 构建引人注目的 3D 地形图时，谁还需要 GIS 呢？

制图 twitter 目前充斥着一些引人注目的单个国家或地区的地形图。这主要是由于可用数据源的扩展和许多无聊的地理书呆子在疫情期间让自己忙碌。有很多工具和方法可以用来生成上面分享的美丽的意大利地图，但是在这篇文章中，我将使用 Python 向您展示一种非常规的方法，希望您会像我一样确信，如果某件事情值得做，就值得用 Python 来做。

有许多可用的数据源，都有不同程度的分辨率。在纯粹的美学层面上，图像的分辨率越高，最终的图像看起来就越好。从精度的角度来看，较低分辨率的数据集也可能会夸大要素，使相对平坦的区域看起来相当多山。

使用高分辨率(7.5 弧秒)和低分辨率(1 弧分)数据源绘制的意大利地形图。作者图片

虽然它们看起来没有什么不同，但左图使用的数据(1)比右图(2)的分辨率高得多，因此相对平坦的区域，例如意大利北部的波河流域，看起来比右图中的实际丘陵要多。因此，今后我们将使用来自美国地质调查局(1)的高分辨率数据。数据可以免费使用，关于数据使用的更多信息可以在下面的链接中找到。数据被分割成 30 x 20 度的 tif 文件，覆盖了全球的不同地区，幸运的是意大利位于其中一个文件中，所以我们不需要担心合并不同的 tif 文件，尽管我将在以后的文章中探讨这个问题。下载数据非常简单，只需点击您想要的图块，然后按照下载说明进行操作。然后把数据复制到你想编程的任何地方。

美国地质勘探局运营的数据交付平台截图。作者图片

首先要做的是使用rasterio打开并读取数据，这是一个相对简单的过程，但我已经包含了代码。然后，重要的是绘制数据，以了解实际情况和您正在处理的问题。有各种各样的函数可以汇总数据，但是我认为它们很少能很好地代替图表。需要注意的是，rasterio.open创建了一个rasterio 数据集对象，它包含高程值的 2D 数组(纬度 x 经度)以及关于所使用的投影和图像范围的信息，这在后面会很重要。read 方法将 rasterio 数据集对象中的数据读入高程值的 2D numpy 数组(在本例中)。2D 数组用于绘图，而rasterio 数据集对象将在以后使用。

原始数据，使用光谱色彩图绘制。作者图片

有一个迫在眉睫的问题需要解决，数据集包含南欧大部分地区和北非部分地区的信息，而实际上，我们想要的只是意大利。幸运的是，rasterio 提供了一种基于地理配准形状(如多边形)裁剪栅格的有用方法。在本例中，我使用了 NaturalEarth (3)数据集的一部分意大利多边形，通过rasterio.mask.mask 函数来裁剪栅格。我加载了 NaturalEarth shapefile 并使用geopandas 提取了意大利的几何图形，然后使用了蒙版功能。mask 函数采用rasterio 数据集对象，并返回所提供的多边形内的栅格部分。

作者图片

现在我们有了仅与意大利相对应的海拔值，但是仍然有一个问题。默认情况下,rasterio.mask.mask将用 0 填充所有不在意大利多边形内的值。虽然这是合理的，但这些 0 将使绘图变得棘手，因为它们在色彩图的底部充当锚，如果 0 和意大利数据中的最小高程值之间有很大的差距，那么您会看到上面看到的图，其中一半色彩图中有真实数据，另一半色彩图中没有真实数据，因为 0 和真实数据的最小值之间有差距。例如，如果 0 和最小值之间的差距等于最小和最大高度之间的差距，那么实际上只有一半的颜色图将用于真实数据。

为了解决这个问题，mask 函数允许您使用关键字nodata明确设置什么值将被应用到不在意大利多边形内的值。我在下面加入了一个函数来解决这个问题。在下面的函数中，我传递了意大利GeoDataFrame 和rasterio 数据集对象。您会注意到掩膜函数被调用了两次，第一次是如上所述调用，不在意大利面内的值被返回为 0。第二次，使用了nodata 参数，不属于意大利面的值被设置为比意大利地形数据集中的最大值大 1(使用第一个掩膜计算)。结果是，我们现在有了一个在值中没有自然间隙的数据集，并且图表开始成形。这个函数返回的还有value_range 变量。这对应于数组中最小和最大值之间的差距，并且在以后构造颜色图时需要。

作者图片

现在我们需要构建一个合适的颜色图。绝对明确地说，本文的目的是向您展示如何生成有趣但至关重要的美丽地形图。如果你计划创造一些军事单位将在战斗中使用的东西，那么我建议你使用更定量的颜色图。在这个例子中，我们将基于意大利国旗构建一个颜色图。下面我用意大利国旗的颜色做了一个。仅仅使用这三种颜色会产生一个中间有太多白色的颜色图，所以我在这个颜色图的 2 号和 4 号位置添加了额外的绿色和红色，以最小化白色的优势。

作者图片

我们仍然需要处理不属于之前创建的意大利高程数据集的值。我们将这些值设置为比高程数据中的最大值大 1。解决方案是用足够的颜色构建一个色彩图，这样意大利高程数据中的每个唯一值都有自己的颜色，然后用我们的背景色替换色彩图中的最后一种颜色(红边)。例如，考虑一个场景，其中最小高程值为 10，最大高程值为 100，我们将非意大利值设置为 101。如果我们用 91 种颜色创建一个色彩映射表，用我们的背景色替换第 91 种颜色，那些非意大利值将被映射到第 91 种颜色，这是我们的背景色。

现在，有了新的颜色图和剪辑数据，我们就可以绘制数据了。

作者图片

虽然我认为这仍然看起来很好，它仍然是 2D 和地形是三维的。所以最后要做的是添加山体阴影来模拟光照在地形上。山体阴影是表面的 3D 表示，通常以灰度渲染。较暗和较亮的颜色表示您在地形模型中看到的阴影和高光。山体阴影通常用作地图中的参考底图，以使数据看起来更加三维，从而在视觉上更加有趣。我们将使用earthpy 山体阴影函数来生成山体阴影数据。有两个参数可以调整，根据数据集的不同，它们会给出明显不同的结果。第一个是azimuth 值，范围为 0-360 度，与光源的发光位置有关。0 度对应于指向正北的光源。第二个是光源所在的altitude ，这些值的范围是 1-90。下面是几个例子，强调了改变这两个值会产生非常不同的结果。

不同的方位角值及其对山体阴影渲染方式的影响。所有的高度都设置为 1。作者图片

不同的高度值及其对山体阴影渲染方式的影响。所有都设置了 180°的方位角。作者图片

最后我决定方位角值为 180，这样我在阿尔卑斯山的南边得到了一个漂亮的阴影，高度值为 1。不过，我会鼓励你自己尝试这些价值观。最后是绘制成品的时候了。首先绘制主要的意大利地形数据，在顶部绘制带有小 alpha 值的山体阴影。

瞧啊。作者图片

放大可以让你真正看到图像的细节！

作者图片

这就是我们的意大利地形图，一幅非常赏心悦目的意大利地形图，随时可以裱起来挂在你的墙上。这种方法可以应用于任何国家或地区，但是如果一个国家的数据跨越多个 tif 文件，则需要额外的步骤。我将在以后的文章中讨论这个问题，所以请订阅以确保你能看到它。

(1)-高分辨率数据源- Danielson，J.J .和 Gesch，D.B .，2011，全球多分辨率地形高程数据 2010 (GMTED2010):美国地质调查局公开文件报告 2011–1073，26 页 doi: 10.5066/F7J38R2N 数据免费使用，关于使用的进一步信息可在此处找到- https://topotools.cr.usgs

(2)—低分辨率数据源- doi:10.7289/V5C8276M

(3)—natural earth—https://www.naturalearthdata/

使用 Seaborn Python 库创建箱线图

原文：https://towardsdatascience/creating-boxplots-with-the-seaborn-python-library-f0c20f09bd57

Seaborn 箱线图快速入门指南

图片来自 Pixabay

箱线图是可视化数据的伟大统计工具，通常用于数据科学项目的探索性数据分析(EDA)阶段。它们为我们提供了数据的快速统计摘要，帮助我们了解数据是如何分布的，并帮助我们识别异常数据点(异常值)。

在这个简短的教程中，我们将看到如何使用流行的 Seaborn Python 库生成箱线图。

什么是箱线图？

一个箱线图是一种基于五个关键数字显示数据分布的图形化和标准化方法:

• “最低”
• 第一个四分位数(第 25 个百分位数)
• 中位数(第二个四分位数/第五十个百分位数)
• 第三个四分位数(第 75 个百分位数)
• “最大值”

最小值和最大值分别定义为 Q1-1.5 * IQR 和 Q3 + 1.5 * IQR。任何超出这些限制的点都被称为异常值。

箱线图的图形描述，突出显示关键组成部分，包括中位数、四分位数、异常值和四分位数间距。作者创建的图像。

箱线图可用于:

• 识别异常值或异常数据点
• 来确定我们的数据是否有偏差
• 了解数据的分布/范围

为了构建箱线图，我们首先从中间值(第 50 百分位)开始。这代表了我们数据中的中间值。

然后在第 25 和第 75 个百分点之间形成一个方框(分别为 Q1 和 Q3)。这个方框表示的范围称为四分位间距(IQR)。

从这个盒子延伸出两条线，也就是众所周知的胡须。这些延伸到 Q1-1.5 * IQR 和 Q3 + 1.5 * IQR，或者延伸到小于该值的最后一个数据点。

任何超出晶须极限的点都称为异常值。

资料组

我们在本教程中使用的数据集是作为 Xeek 和 FORCE 2020 (Bormann et al .，2020) 举办的机器学习竞赛的一部分使用的训练数据集的子集。

完整的数据集可以通过以下链接获得:https://doi/10.5281/zenodo.4351155。

竞赛的目的是利用测井测量从现有的标记数据预测岩性。完整的数据集包括来自挪威海的 118 口井。

此外，您可以从 GitHub 资源库下载本教程中使用的数据子集:

https://github/andymcdgeo/Petrophysics-Python-Series

锡伯恩图书馆

Seaborn 是一个建立在 matplotlib 之上的高级数据可视化库。它为创建更高级的绘图提供了更容易使用的语法。与 matplotib 相比，默认数字也更具视觉吸引力

使用 Seaborn 构建箱式地块

导入库和数据

首先，我们首先需要导入将要使用的库: pandas 用于加载和存储我们的数据，以及 Seaborn 用于可视化我们的数据。

import seaborn as sns
import pandas as pd

导入库后，我们可以从 CSV 文件导入数据并查看文件头。

df = pd.read_csv('Data/Xeek_train_subset_clean.csv')
df.head()

在数据集内，我们有关于油井、地质分组和地层的详细信息，以及我们的测井测量。如果您不熟悉这些数据，请不要担心，因为下面的技术可以应用于任何数据集。

创建简单的箱线图

我们可以生成第一个箱线图，如下所示。在括号内，我们传入想要从数据帧中访问的列。

sns.boxplot(x=df['GR']);

由 Seaborn 生成的简单箱线图。图片由作者提供。

我们也可以旋转我们的绘图，使方块是垂直的。为了做到这一点，我们为y而不是x提供一个值。

sns.boxplot(y=df['GR']);

Seaborn 生成的垂直箱线图。图片由作者提供。

我们可以结合x和y参数来创建多个盒状图。在本例中，我们将 y 轴设置为 GR(伽马射线),它将被 LITH(岩性)列分割成单独的箱线图。

sns.boxplot( x=df['LITH'], y=df['GR']);

根据岩性划分的伽马射线数据的海底生成的垂直箱线图。图片由作者提供。

从表面上看，我们现在有一个按岩性划分的多个箱线图。不过，有点乱。我们可以整理一下，用几行额外的代码使它变得更好。

整理默认的 Seaborn 箱线图

更改图形大小和旋转 x 轴标签

由于 Seaborn 构建在 matplotlib 之上，我们可以使用 matplotlib 的功能来提高我们的绘图质量。

使用 matplotlibs .subplots函数，我们可以使用figsize定义图形的大小，还可以调用图形的元素，如 xticks。在下面的例子中，我们将图形大小设置为 10 乘 10，并将xtick标签的旋转角度设置为 90 度。

import matplotlib.pyplot as pltfig, ax = plt.subplots(1, figsize=(10, 10))sns.boxplot(x=df['LITH'], y=df['GR']);
plt.xticks(rotation = 90)
plt.show()

当我们运行这段代码时，我们得到了一个更容易阅读的图。

在定义图形大小和旋转 x 轴标签后，由按岩性划分的伽马射线数据的海底生成的垂直箱线图。图片由作者提供。

更改 Seaborn 箱线图的图形大小方法 2

改变 Seaborn 地块大小的另一种方法是调用sns.set(rc={“figure.figsize”:(10, 10)})。使用这个命令，我们可以很容易地改变绘图的大小。

然而，当我们使用这条线时，它会将所有后续的图设置为这个大小，这可能并不理想。

sns.set(rc={"figure.figsize":(10, 10)})
sns.boxplot( x=df['LITH'], y=df['GR']);

Seaborn 在使用 sns.set()更改图形大小后生成的垂直箱线图。图片由作者提供。

设计 Seaborn 箱线图

Seaborn 提供了五种预设样式(深色网格、白色网格、深色、白色和刻度),可以快速轻松地改变整个地块的外观。

要使用其中一种样式，我们调用sns.set_style()，并传入其中一种样式作为参数。在本例中，我们将使用白色网格。

sns.set_style('whitegrid')
sns.boxplot( y=df['LITH'], x=df['GR']);

当我们运行代码时，我们得到了下面的图。请注意，我还交换了 x 轴和 y 轴，这样方框就可以水平绘制了。

Seaborn 箱线图显示应用 Seaborn 主题后不同岩性的伽马射线值。图片由作者提供。

如果我们想改变方框图的颜色，我们只需使用color参数，并传入我们选择的颜色。

sns.set_style('whitegrid')
sns.boxplot( y=df['LITH'], x=df['GR'], color='red');

这将返回以下带有红框的图。

Seaborn 盒状图显示了设置盒的颜色后不同岩性的伽马射线值。图片由作者提供。

除了固定的颜色，我们还可以对箱线图应用调色板。这将使每个盒子有不同的颜色。在这个例子中，我们将调用蓝调调色板。你可以在这里了解更多关于 Seaborn 调色板的细节。

sns.set_style('whitegrid')
sns.boxplot( y=df['LITH'], x=df['GR'], palette='Blues');

Seaborn 盒状图显示了应用调色板给盒子着色后不同岩性的伽马射线值。图片由作者提供。

设定 Seaborn 图的 X 轴和 Y 轴标签的样式

默认情况下，Seaborn 将使用轴标签的列名。

首先，我们必须将我们的箱线图赋给一个变量，然后访问所需的函数:set_xlabel、set_y_label和set_title。当我们调用这些方法时，我们还可以设置字体大小和字体粗细。

p = sns.boxplot(y=df['LITH'], x=df['GR'])
p.set_xlabel('Gamma Ray', fontsize= 14, fontweight='bold')
p.set_ylabel('Lithology', fontsize= 14, fontweight='bold')
p.set_title('Gamma Ray Distribution by Lithology', fontsize= 16, fontweight='bold');

当我们运行这段代码时，我们得到了一个更好看的图，带有易于阅读的标签。

对标题、x 轴和 y 轴标签应用格式后的 Seaborn Boxplot。图片由作者提供。

设计 Seaborn 箱线图的异常值

除了能够设计盒子的样式，我们还可以设计离群值的样式。为了做到这一点，我们需要创建一个变量字典。在下面的例子中，我们将改变标记的形状(marker)、标记的大小(markersize)、异常值的边缘颜色(markeredgecolor)和填充颜色(markerfacecolor)以及异常值透明度(alpha)。

flierprops = dict(marker='o', markersize=5, markeredgecolor='black', markerfacecolor='green', alpha=0.5)p = sns.boxplot(y=df['LITH'], x=df['GR'], flierprops=flierprops)
p.set_xlabel('Gamma Ray', fontsize= 14, fontweight='bold')
p.set_ylabel('Lithology', fontsize= 14, fontweight='bold')
p.set_title('Gamma Ray Distribution by Lithology', fontsize= 16, fontweight='bold');

更改默认异常值(flier)属性后的 Seaborn 箱线图。图片由作者提供。

摘要

在这个简短的教程中，我们看到了如何使用 Python Seaborn 库来生成测井数据的基本箱线图，并按岩性对其进行拆分。与 matplotlib 相比，Seaborn 直接提供了更好的绘图。

我们可以使用箱线图来可视化我们的数据，并了解数据的范围和分布。然而，它们是识别数据异常值的优秀工具。

感谢阅读。在你走之前，你一定要订阅我的内容，把我的文章放到你的收件箱里。 你可以在这里做！或者，您可以 注册我的简讯 免费获取更多内容直接发送到您的收件箱。

其次，通过注册会员，你可以获得完整的媒介体验，并支持我自己和成千上万的其他作家。它每个月只花你 5 美元，你可以完全接触到所有令人惊叹的媒体文章，也有机会用你的写作赚钱。如果你用 我的链接报名，你直接用你的一部分费用支持我，不会多花你多少钱。如果你这样做了，非常感谢你的支持！

参考

博尔曼，彼得，奥桑德，彼得，迪里布，法哈德，曼拉尔，投降，&迪辛顿，彼得。(2020).机器学习竞赛 FORCE 2020 井测井和岩相数据集[数据集]。芝诺多。http://doi/10.5281/zenodo.4351156

使用 Python 的叶库创建 Choropleth 地图

原文：https://towardsdatascience/creating-choropleth-maps-with-pythons-folium-library-cfacfb40f56a

如何用 Python 制作不同数据结构的 choropleths

【2019 年 4 月纽约市可出租公寓一览表(GitHub)

C horopleth 地图用于显示地理区域的数据变化( 人口教育 )。我使用 choropleths 显示了纽约市不同邮政编码的可用出租公寓数量，并显示了给定时期每个邮政编码的抵押贷款交易数量。Python 的叶子库使用户能够构建多种定制地图，包括 choropleths，你可以将它们作为.html文件与不知道如何编码的外部用户共享。

加载和查看地理数据

美国政府网站通常有创建地图所需的地理数据文件。纽约市的 OpenData 网站和美国人口普查局的网站有多种数据类型的地理边界文件。Python 允许你加载多种文件类型，包括 geo JSON*(*.geojson*)文件和 shape file(*.shp*)*。这些文件包含给定位置的空间边界。

关于folium.Choropleth()方法的 folio 文档说明，geo_data参数接受 GeoJSON 几何作为字符串来创建映射，“URL、文件路径或数据 (json、dict、geopandas 等)到您的 GeoJSON
几何“()folio 文档 ) 。无论我们如何加载文件，我们都必须用这种方法转换几何数据才能正常工作。该方法的key_on参数将每个特定位置 *(GeoJSON 数据)的数据与该位置(即人口)*的数据绑定。

杰奥森

GeoJSON 文件存储几何形状，在这种情况下是位置的边界及其相关属性。例如，加载带有纽约市邮政编码边界的 GeoJSON 文件的代码如下:

# Code to open a .geojson file and store its contents in a variable**with** open ('nyczipcodetabulationareas.geojson', 'r') **as** jsonFile:
    nycmapdata **=** json**.**load(jsonFile)

变量nycmapdata包含一个至少有两个键的字典，其中一个键叫做features，这个键包含一个字典列表，每个字典代表一个位置。第一个位置的主要 GeoJSON 结构摘录如下:

{'type': 'FeatureCollection',
 'features': [{'type': 'Feature',
   'properties': {'OBJECTID': 1,
    'postalCode': '11372',
    'PO_NAME': 'Jackson Heights',
    'STATE': 'NY',
    'borough': 'Queens',
    'ST_FIPS': '36',
    'CTY_FIPS': '081',
    'BLDGpostal': 0,
    '@id': 'http://nyc.pediacities/Resource/PostalCode/11372',
    'longitude': -73.883573184,
    'latitude': 40.751662187},
   'geometry': {'type': 'Polygon',
    'coordinates': [[[-73.86942457284177, 40.74915687096788],
      [-73.89143129977276, 40.74684466041932],
      [-73.89507143240859, 40.746465470812154],
      [-73.8961873786782, 40.74850942518088],
      [-73.8958395418514, 40.74854687570604],
      [-73.89525242774397, 40.748306609450246],
      [-73.89654041085562, 40.75054199814359],
      [-73.89579868613829, 40.75061972133262],
      [-73.89652230661434, 40.75438879610903],
      [-73.88164812188481, 40.75595161704187],
      [-73.87221855882478, 40.75694324806748],
      [-73.87167992356792, 40.75398717439604],
      [-73.8720704651389, 40.753862007052064],
      [-73.86942457284177, 40.74915687096788]]]}}, ... ]}

folium.Choropleth()方法的key_on参数要求用户以字符串形式引用 GeoJSON 文件中位置字典的唯一索引键:

key_on ( 字符串，默认无)—geo _ data geo JSON 文件中的变量，数据绑定到该变量。必须以“feature”开头，并且采用 JavaScript 异议表示法。例如:“feature.id”或“feature.properties.statename”。

在上面的例子中，索引键是邮政编码，与每个位置相关的数据也必须有一个邮政编码索引键或列。上例中的key_on参数是以下字符串:

‘feature.properties.postalCode’

注意:字符串的第一部分必须始终是单数单词 *feature* ，它不像父字典那样包含每个单独位置字典的列表。

key_on参数是访问每个特定位置的properties键。properties键本身包含一个有 11 个键的字典，在这种情况下,postalCode键是索引值，它将几何形状链接到我们想要绘制的任何值。

地质公园

另一种加载地理数据的方法是使用 Python 的 GeoPandas 库 ( 链接 ) 。这个库在加载 Shapefile 时很有用，shape file 在美国人口普查网站 ( 制图边界文件— Shapefile ) 上提供。GeoPandas 的工作方式类似于 Pandas，只是它可以存储和执行几何数据的功能。例如，用美国所有州的边界加载 shapefile 的代码如下:

# Using GeoPandasimport geopandas as gpd
usmap_gdf = gpd.read_file('cb_2018_us_state_500k/cb_2018_us_state_500k.shp')

us map _ GDF 数据帧的头部

如果您调用 Jupyter 笔记本中第一行的 (Mississippi) 几何列，您将看到以下内容:

usmap_gdf[“geometry”].iloc[0]

当调用特定的几何值时，您会看到一个几何图像，而不是表示形状边界的字符串，上面是第一行(密西西比)的几何值

与 GeoJSON 字典的内容不同，没有用于访问内部字典的features键，也没有properties列。folium.Choropleth()方法的key_on参数仍然要求字符串的第一部分为feature，但是该方法将引用 GeoPandas 数据帧中的列，而不是引用 GeoJSON 的位置字典。在这种情况下，key_on参数将等于“feature.properties.GEOID”，其中GEOID是包含将我们的数据绑定到地理边界的唯一州代码的列。GEOID列有前导零，加州GEOID是06。您也可以使用STATEFP列作为索引，确保使用的列、格式和数据类型都是一致的。

查看 Choropleth 的群体数据

地理数据和要绘制的相关数据可以作为两个单独的变量存储，也可以一起存储。重要的是跟踪列的数据类型，并确保索引(*key_on*)列对于地理数据和位置的关联数据是相同的。

我访问了美国人口普查 API 的美国社区调查 ( 链接 ) 和人口估计与预测 ( 链接 ) 表格，得到 2019 年到 2021 年的人口和人口统计数据。数据帧的标题如下:

美国人口普查数据框架的头

我将数据保存为一个.csv文件，在某些情况下，这将改变列的数据类型；例如，字符串可以变成数值。调用.info() 时的数据类型如下:

将数据框保存和加载为 CSV 文件之前和之后的人口普查数据的数据类型

另一个需要注意的重要事情是，在加载数据帧后，state列中的所有前导零都不会出现。这必须得到纠正；id 必须匹配并且是相同的数据类型*(即它不能在一个数据框中是整数而在另一个数据框中是字符串)*。

基本的 Choropleth 映射有五种不同的方式

如上所述，follow 允许您使用地理数据类型创建地图，包括 GeoJSON 和 GeoPandas。这些数据类型需要被格式化，以便与 least 库一起使用，为什么会出现某些错误并不总是直观的(对我来说是*，至少是*)。以下示例描述了如何准备地理数据*(在本例中为美国州界)和相关的绘图数据(各州的人口)*以用于folium.Choropleth()方法。

方法 1:使用 Pandas 和 GeoJSON，不指定 ID 列

这种方法非常类似于文档中的 choropleth 地图的例子。该方法使用包含州边界数据的 GeoJSON 文件和 Pandas 数据帧来创建地图。

当我开始使用 GeoPandas 文件时，我需要使用 GeoPandas 的to_json() 方法将其转换为 GeoJSON 文件。提醒一下 usmap_gdf GeoPandas 数据帧看起来像:

us map _ GDF 数据帧的头

然后我应用.to_json()方法，并指定我们从数据帧中删除id，如果它存在的话:

usmap_json_no_id = usmap_gdf.to_json(drop_id=True)

注意: *usmap_json_no_id* 是本场景中存放 json 字符串的变量

这个方法返回一个字符串，我对它进行了格式化，以便于阅读，并显示在下面的第一组坐标上:

'{"type": "FeatureCollection",
  "features": [{"type": "Feature",
  "properties": {"AFFGEOID": "0400000US28", 
                 "ALAND": 121533519481,
                 "AWATER": 3926919758,
                 "GEOID": 28,
                 "LSAD": "00",
                 "NAME": "Mississippi",
                 "STATEFP": "28",
                 "STATENS": "01779790",
                 "STUSPS": "MS"},
  "geometry": {"type": "MultiPolygon",
               "coordinates": [[[[-88.502966, 30.215235]'

注意:“属性”字典没有名为“id”的键

现在，我们准备将新创建的 JSON 变量与上一节中获得的美国人口普查数据框架连接起来，其标题如下:

美国人口普查数据框架的头部，下面称为 all_states_census_df

使用 folium 的Choropleth()方法，我们创建地图对象:

创建一个Choropleth with a GeoJSON variable which does not specify an id的代码

geo_data参数被设置为新创建的*usmap_json_no_id*变量，而data参数被设置为 all_states_census_df 数据帧。由于在创建 GeoJSON 变量时未指定 id，key_on参数必须引用地理数据中的特定关键字，并且其工作方式类似于字典( GEOID 是“属性”关键字的值)。在这种情况下，GEOID键保存州代码，该代码将州几何边界数据连接到 all_states_census_df 数据帧中相应的美国人口普查数据。下面是 choropleth:

由上述方法得到的 choropleth

方法 2:使用 Pandas 和 GeoJSON，并指定 ID 列

除了在调用.to_json()方法之前使用一个索引之外，这个过程与上面的几乎完全相同。

在上面的例子中，usmap_gdf数据帧没有索引，为了纠正这个问题，我将把索引设置为GEOID列，然后立即调用.to_json()方法:

usmap_json_with_id = usmap_gdf.set_index(keys = “GEOID”).to_json()

直到第一个州的数据的第一对坐标，结果字符串如下:

'{"type": "FeatureCollection",
  "features": [{"id": "28",
                "type": "Feature",
                "properties": {"AFFGEOID": "0400000US28",
                               "ALAND": 121533519481, 
                               "AWATER": 3926919758,
                               "LSAD": "00",
                               "NAME": "Mississippi",
                               "STATEFP": "28",
                               "STATENS": "01779790",
                               "STUSPS": "MS"},
                "geometry": {"type": "MultiPolygon",
                             "coordinates": [[[[-88.502966, 30.215235],'

“properties”字典不再有GEOID键，因为它现在作为一个名为id的新键存储在外层字典中。您还应该注意到，id值现在是一个字符串，而不是一个整数。如前所述，您必须确保连接数据的数据类型是一致的。如果涉及前导零和尾随零，这会变得很乏味。为了解决这个问题，我从all_states_census_df中的state列创建了一个名为state_str的新列:

all_states_census_df[“state_str”]=all_states_census_df[“state”].astype(“str”)

现在我们可以创建 choropleth:

创建 c 的代码horopleth with a GeoJSON variable which specifies an id

该代码与之前使用的代码的区别在于key_on参数引用的是id而不是properties.GEOID。生成的地图与方法 1 完全相同:

由上述方法得到的 Choropleth

方法 3:使用熊猫和 GeoPandas 的 Python 要素集合

该方法从具有__geo_interface__ 属性的原始 GeoPandas 数据帧创建 GeoJSON like 对象( python 要素集合)。

我将usmap_gdfdata frame(US geographic data)的索引设置为STATEFP列，该列以字符串形式存储州 id，并以零开头:

usmap_gdf.set_index(“STATEFP”, inplace = True)

然后，我通过添加一个前导零，在all_states_census_df数据框(美国人口普查数据)中创建了一个匹配列:

all_states_census_df[“state_str”] = all_states_census_df[“state”].astype(“str”).apply(lambda x: x.zfill(2))

最后，我使用us_data_gdf GeoPandas dataframe 的__geo_interface__ 属性来获取几何状态边界的 python 特性集合，存储为一个字典，类似于前两种方法:

us_geo_json = gpd.GeoSeries(data = usmap_gdf[“geometry”]).__geo_interface__

下面是us_geo_json变量的摘录:

{'type': 'FeatureCollection',
 'features': [{'id': '28',
   'type': 'Feature',
   'properties': {},
   'geometry': {'type': 'MultiPolygon',
    'coordinates': [(((-88.502966, 30.215235), ...))]

最后，我们创建了 choropleth:

创建一个 c horopleth with a GeoPanda's **__geo_interface__**属性的代码

这张地图看起来和上面的一样，所以我把它排除了。

方法 4:使用 Geopandas 的几何类型列

在这里，我们坚持 GeoPandas。我创建了一个名为us_data_gdf的 GeoPandas 数据框架，它将几何数据和普查数据组合在一个变量中:

us_data_gdf = pd.merge(left = usmap_gdf,
                       right = all_states_census_df,
                       how = "left", 
                       left_on = ["GEOID", "NAME"],
                       right_on = ["state", "NAME"]
                       )

注:all _ States _ Census _ df是美国人口普查数据的熊猫数据帧，US map _ GDF是存储州几何边界数据的 GeoPandas 数据帧。

使用 GeoPandas 数据帧创建 choropleth 的代码如下:

使用 GeoPandas 数据帧创建 choropleth 的代码

在上面的例子中，geo_data参数和data参数都引用相同的 GeoPandas 数据帧，因为信息存储在一个地方。因为我没有设置索引，所以key_on参数等于“feature.properties.GEOID”。即使使用 GeoPandas，folium 也需要使用key_on参数，就像引用一个类似字典的对象一样。

和以前一样，这张地图看起来和上面的一样，所以我把它排除了。

方法 5:使用 Geopandas 几何类型和 Branca

在这里，我们使用 Branca 库和 leav 的示例创建了一个更加时尚的地图。除了安装之外，Branca 的第一步是创建一个ColorMap对象:

colormap = branca.colormap.LinearColormap(
    vmin=us_data_gdf["Total_Pop_2021"].quantile(0.0),
    vmax=us_data_gdf["Total_Pop_2021"].quantile(1),
    colors=["red", "orange", "lightblue", "green", "darkgreen"],
    caption="Total Population By State",
)

在上面的代码中，我们访问了branca.colormap.LinearColormap 类。在这里，我们可以设置我们使用的颜色以及色阶的值。对于这个 choropleth，我希望颜色与美国人口普查数据中最低和最高的人口值成比例。为了设置这些值，我使用了上述的vmin和vmax参数。如果我忽略了这一点，那么色标中将考虑没有值的区域，没有这些设置参数的结果如下:

没有设置**vmin**和***vmax***参数的 Branca choropleth

一旦创建了ColorMap对象，我们就可以创建一个 choropleth ( 完整代码在下面):

使用 GeoPandas 数据框架和 Branca 库创建 choropleth

我修改了 leav 网站上的例子，使用了us_data_gdf GeoPandas 数据框架。该示例允许我们排除地理数据中不具有相关人口普查数据的部分(显示为透明的)(如果一个州的人口为空，那么除非被排除，否则 choropleth 上的颜色将为黑色。生成的 choropleth 如下:

用麦麸和地瓜做成的鸡蛋饼

布兰卡是可定制的，但如何使用它的解释是少之又少。其存储库的自述文件指出:

没有文档，但是你可以浏览范例库。

你必须练习使用它来制作你想要的那种地图。

摘要

对那些有或没有编码知识的人来说，叶子可以用来制作信息丰富的地图，如 choropleths。政府网站通常拥有为您的数据创建位置边界所需的地理数据，这些数据也可以从政府网站获得。理解你的数据类型和文件类型是很重要的，因为这会导致不必要的麻烦。这些地图是高度可定制的，例如你可以添加工具提示来注释你的地图。充分利用这个图书馆的潜力需要实践。

我的这篇文章可以在 这里找到 。快乐编码。

使用优步的 H3 六边形创建一致的空间计算

原文：https://towardsdatascience/creating-consistent-spatial-calculations-with-ubers-h3-hexagons-1af032802a77

将人口普查数据添加到 H3 六边形以进行空间分析

唐尼·姜在 Unsplash 上的照片

所有的分析技术、观点和代码都是独立开发和我自己的工作。

获取一个城市或邮政编码的人口统计数据非常简单——有几个 API 和数据下载工具可以做到这一点。但是，由于许多原因，这些地理位置有局限性(举几个例子，邮政编码是 USPS 递送的功能，城市边界是在几个世纪前建立的——并不完全是为数据分析设计的)。

我测试了几种不同的方法来克服这些边界的任意性——从创建统一的盒子到开发驾驶(或行走)时间多边形，这取决于用例。

我最近开始意识到 H3 的好处，这是优步在 2018 年推出的一套全球六边形。主要思想是六边形是最有利于空间分析的形状，可用于创建一致的计算。

示例 H3 六边形；作者图片

每个六边形都有一系列更小的六边形，这些六边形(大部分)位于另一个六边形的内部，这创建了一个可用于一致参考和分析的层次结构，一直到边长为 2 英尺。

作者图片

这种方法适用于点数据——骑手是在哪里上车的？离开了？任何有经纬度的地址或兴趣点都可以映射成六边形。

如果我们想利用六边形分析的优势，但使用人口统计数据，会发生什么？我发现了很多关于使用六边形的信息，但是没有引入人口统计数据，所以我开始将两者结合起来。

处理人口普查数据时，稳健数据报告提供的最详细视图是人口普查区块组(在此阅读更多关于区块组的信息)。块组是通常有 600 到 3，000 人的边界，这使得它们可以同质学习和使用(这么小的区域往往有相似的人居住在其中)。

然而，因为它们受人口限制，所以大小缺乏一致性。例如，看看下面北卡罗来纳州罗利市中心的人口普查区块组，根据人口分布，有些较小，有些较大。

相比之下，六边形(为了有效分析，我随机选择了一个尺寸，但我们可以根据需要选择更小或更大的尺寸)，大小一致，边界一致。

北卡罗来纳州罗利市中心的人口普查区块组(左)和 H3 六边形(右)

因此，问题变成了，我们如何将人口普查区块组中可用的人口统计信息映射到六边形的统一结构？我们如何从人口普查区块组的粒度和六边形的一致性中获益？

总体情况——我们如何解决这个问题？

假设我们想要找到每个六边形内的人口，更具体地说，让我们看一个六边形。我们可以用相应的块组边界覆盖六边形边界，得到类似下面这样的结果…

覆盖在交叉六边形上的块群总体；作者图片

完全位于一个六边形内的块组是容易的，我们考虑将全部人口分配到六边形(人口 708、969 和 1476 在六边形的中间)。然而，由于形状不匹配，许多有部分重叠，这只会更常见于较小的六边形。

对于与六边形相交的任何块组，我们可以计算两个形状之间的面积重叠，并根据此计算按比例分配人口。如果一半的块组重叠，我们认为一半的人口是在六边形。

这种假设并不完美——人口并不均匀地分布在一个街区组中——但出于实用目的，它为我们提供了一种“足够好”的方法来组合这两种地理形状。

实施

我将用 R 编程语言的代码来演示，但这种方法在 Python 或其他语言中在概念上是相同的。

我们需要两条主要信息来进行分析——带有相应人口统计数据的街区组形状和六边形形状。

在 R 中，我们可以利用tidysensus 包从美国人口普查 API 中获取数据(您将需要一个免费的开发者密钥)。[get_acs](https://www.rdocumentation/packages/tidycensus/versions/1.2.1/topics/get_acs) https://www.rdocumentation/packages/tidycensus/versions/1.2.1/topics/get_acs功能是获取块组数据的最简单方法，其典型语法如下:

get_acs(geography = 'block group', year = 2020, variables = c(total_population = 'B01001_001'), state=37, county=183, output='wide')

可选参数定义了我们需要输出的块组、2020 年 5 年 ACS 的数据、总人口变量标识符、北卡罗来纳州(FIPS 37)、威克县(代码 183)和宽输出格式。

get_acs 函数的输出；作者图片

从那里，我们需要为每个块组匹配一个形状(边界)。我们可以使用来自 tigris 包的block_groups函数，并提供州和县标识符，这将返回一个包含几何多边形的形状文件。

获取北卡罗来纳州威克县每个街区组的形状文件的示例语法；作者图片

通过大地水准面的简单连接可以将 shapefile 与人口统计信息放在一起，以创建单个块组文件。然后使用 sf 包将输出转换回 shapefile(稍后将详细介绍),并转换到坐标参考 4326。

acs_block_data %>%
  left_join(wake_block_groups, by='GEOID') %>%
  st_as_sf() %>%
  st_transform(4326)

为了获得威克郡的所有六边形，我获取了威克郡的 shapefile(也使用了 tigris 包，但是使用了counties函数)，然后利用 h3jsr 包来访问 H3 的六边形。

polyfill函数查找给定边界内的所有六边形(这里是威克郡，罗利的家乡)，而h3_to_polygon函数给出每个六边形的形状。res = 7参数定义了六边形的粒度级别(15 是最精确的，0 是最大的)。

wake_hex <- polyfill(wake_boundary, res=7, simple=F)
wake_hex_all <- h3_to_polygon(unlist(wake_hex$h3_polyfillers), simple=F)

输出是一系列六边形参考 id 和多边形形状。

作者图片

现在，对于每个六边形，我们想要:

• 查找与六边形的边界相交(重叠)的任何块组
• 找出块组的重叠百分比
• 根据重叠百分比分配块组中的人口

R 中的 sf 包是空间分析的首选包。st_intersection函数查找两个形状的重叠区域，下面的 mutate 函数计算重叠百分比。此代码的输出将是所有重叠的块组，以及它们的重叠区域。

st_intersection(wake_hex_all[1,1], wake_blocks) %>%
mutate(intersect_area = st_area(.),
       pct_intersect = intersect_area / (ALAND + AWATER) )

带有附加数据操作的 st_intersection 函数的输出示例；作者图片

最后一步是将总体转换为数字，将总体乘以百分比交集以获得调整后的总体，然后将调整后的总体相加得到六边形中的总数。

st_intersection(wake_hex_all[1,1], wake_blocks) %>%
mutate(intersect_area = st_area(.),
       pct_intersect = intersect_area / (ALAND + AWATER),
       pop_total = as.numeric(total_populationE),
       pop_adj = pop_total * pct_intersect ) %>%
summarize(total_population = sum(pop_adj, na.rm = T)

例如，在上面的图像中，具有约 11%重叠的 2，911 人口成为约 317 的调整人口。

回到我们之前的例子，这种方法给出了估计的六边形人口为 9，980。再次观察每个区块组中的人口，视觉比较表明这可能不会太远。

底层块组群体的群体计算示例

当六边形/地块组与受保护的土地(如公园或水体)有较大重叠时，这种方法有局限性，这对于更偏远地区的小六边形最有影响。考虑到大多数城市和近郊的人口密度，我认为这种方法在大多数情况下会有方向性的帮助。

扩展分析

人口只是用作示例，但是可以使用在块组级别报告的任何内容。六边形的家庭收入中位数是多少？有多少户人家？有孩子的家庭比例是多少？添加到在get_acs函数中收集的变量中，可以启用如何分析和汇总数据的其他选项。

这也可以与可以聚合到六边形中的任何兴趣点数据相结合，例如，通过查找给定兴趣点的高密度区域，然后描述该六边形的人口构成。

虽然并不完美，但这提供了一种为空间分析和趋势创建一致计算的方法。总的来说，我们希望数据在方向上是有用的——如果我们计算出一个六边形的人口比另一个六边形多 20%,我们可以合理地确定第一个六边形的人口更多，即使不是正好 20%。

有疑问或者对这些话题感兴趣？ 在 LinkedIn 上和我联系，或者在 jordan@jordanbean 给我发邮件。

在 Django 中创建动态长度表单

原文：https://towardsdatascience/creating-dynamic-length-forms-in-django-53709c23464e

费萨尔在 Unsplash 上拍摄的照片

讨论如何创建一个 Django 表单，它可以有动态数量的输入字段

最近，我正在解决一个问题，创建一个应用程序，可以使用搜索词收集维基百科数据，并收集文章简介。该应用程序进一步拆分所有句子，并允许用户将文章中每个句子的翻译保存在数据库中。这个应用程序将帮助任何人通过使用搜索词来注释或翻译维基百科的内容。

该应用程序主要是使用 Django 和 HTML 模板构建的。我使用了维基百科 API 来收集数据。最终输出可以在这里找到和源代码在我的 GitHub 上。

https://github/prakharrathi25/wikipedia-translator

下面，我还展示了应用程序的工作原理。输入搜索词后，API 抓取结果并显示在页面上。

输入搜索词后的应用程序输出(图片由作者提供)

搜索之后，应用程序的第二部分对文本进行标记(将文本分解成句子)并为翻译做准备。第二列中的每一行都是一个表单输入，用户可以在其中输入他们的翻译并保存。

准备翻译的标记化句子(图片由作者提供)

我面临的问题——动态表单长度

在构建项目时，我遇到了一个有趣的新问题— **创建动态表单。**这是因为在对介绍段落进行标记后，我们需要创建一个表单，它可以根据句子的数量拥有不同数量的输入字段。然后，我们需要在后端处理所有的输入字段，并保存翻译。一个人也可以只输入一些句子的翻译，而把其余的留空。

我能够使用我最喜欢的 web 开发框架之一来解决这个问题——Django！

解决方案

现在，我不会深入研究 Django 的基础知识，因为这篇文章的目标读者是那些有编写 Django 代码经验的人。如果你想入门，这里有个好教程。

我将从定义我在这个应用程序中使用的模型开始。我使用的主要模式有两种——项目和**句子。**每个新的文章标题都保存在项目中，标题的介绍被标记成句子。然后，每个句子都保存在句子模型中，并带有返回到项目的外键。

显示句子和项目模型及其字段的模型文件

接下来，我将展示我用来在同一个表单中动态显示多个输入字段的 HTML 文件。这意味着输入字段的数量将根据文章引言段落中的句子数量而变化。

显示动态条目的表格和表单(代码片段链接

这是 Django 的模板格式，你可以在这里阅读更多关于它的内容。这个表格循环了我通过视图上下文传递的句子的数量，稍后我会再看一遍。因此，根据句子的数量，我可以在表单中添加更多的字段，并从中获取输入。我做的另一件事是使输入字段 name 属性动态化，这样我就能够从单个输入中收集数据。

用户在文本区输入的数据处理如下。

还有一个 GET 请求处理程序，我已经省略了。完整的功能可以在这里找到。POST 请求处理器处理每个输入字段，并使用文本区域的 name 属性收集翻译。这是一个非常简单而优雅的解决方案，我花了一些时间才想出来。除了这个线程之外，这个没有很多堆栈溢出的答案。

结论

这不应该是一篇精心设计的文章。目的是分享一种解决 Django 论坛上经常提出的问题的新方法。希望这在将来能派上用场。

[## 使用 Streamlit 创建多页面应用程序(高效！)

towardsdatascience](/creating-multipage-applications-using-streamlit-efficiently-b58a58134030) https://levelup.gitconnected/5-new-features-in-python-3-11-that-makes-it-the-coolest-new-release-in-2022-c9df658ef813

用强化学习和虚幻引擎创建涌现行为

原文：https://towardsdatascience/creating-emergent-behaviors-with-reinforcement-learning-and-unreal-engine-4cd89c923b7f

使用虚幻引擎和免费的 MindMaker 机器学习插件生成人工智能角色的紧急行为

在下面的文章中，我将讨论如何使用虚幻引擎、强化学习和免费的机器学习插件 MindMaker 在人工智能角色中生成紧急行为。目的是让感兴趣的读者可以将此作为在他们自己的游戏项目或具体的人工智能角色中创建紧急行为的指南。

什么是突现行为，为什么要使用它？

首先是关于突现行为的初级读本。突现行为是指没有预先编程，而是响应某些环境刺激而有机发展的行为。突现行为是许多生命形式的共性，是进化本身的一种功能。这也是最近具体化人工代理的一个特征。当一个人采用涌现行为方法时，他不会为人工智能严格地编程特定的动作，而是允许它们通过一些适应性算法“进化”，如遗传编程、强化学习或蒙特卡罗方法。在这样的设置中——行为在开始时没有被预言，而是被允许基于级联的一系列事件“出现”,在某种程度上依赖于机会。

为什么人们会选择使用突现行为？一个主要原因是，涌现行为方法可以创造出行为更类似于碳基生命形式的人工智能代理，更不可预测，并显示出更大的战略行为多样性。

为娱乐或人类互动而设计的各种各样的具体化人工智能代理可以从紧急行为中受益，以便对它们的用户来说看起来不那么静态，更有吸引力。《自然》杂志上最近的一篇期刊文章证明了机器人在反应时间和运动模式方面表现出更多的可变性，被认为更像人类。

紧急行为也可以克服用传统编程方法创造的其他形式的人工智能所不能克服的障碍。通过允许人工智能发现开发人员没有预见到的解决方案，代理可以探索比传统编程方法更广阔的解决方案空间。甚至有可能出现的行为方法可能导致人工一般智能的产生，这是一种可以与人类拥有的技能多样性相媲美的人工智能。

使用紧急行为方法有一个警告，但是，如果你想要创建一个非常特定的行为，例如精确地重复模仿特定动物的动作，紧急行为技术可能不是最好的选择。突现行为不太适合不需要变化的重复性高保真任务。这是为什么人类和其他有机生命形式可能在高度重复的任务中挣扎的原因之一，我们天生是可变的。此外，当试图复制特定的行为时，可能很难重现导致该行为模式在现实世界中出现的相同偶然事件。例如，某种生物适应或行为可能会进化出各种各样的解决方案来实现，其结果受偶然事件的影响。因此，并不是所有的新方法都一定会得到相同的解决方案。一个例子是趋同进化在几个场合下重建眼睛的方式，但不一定是完全相同的形式。

涌现行为方法有利于创造多种多样的行为，其中一些行为可能与真实动物的行为相似，但不太可能完全相同。另一方面，人们可以期待用这些方法产生大量有趣和独特的行为。例如，人们可以看到由强化学习代理发现的穿越风景的各种各样的移动方法。

创造突发行为的第一步是决定你想让哪种行为成为突发行为，以及该行为的目标是什么。从某种意义上说，这是容易的部分。下一部分稍微复杂一点，需要概述一些涌现理论。

不同类型的紧急行为

当我们讨论突现的话题时，我们需要区分开放式突现和静态或“固定点”突现。在定点涌现中，一开始可能会出现各种各样的行为，但它们会稳定地收敛到一个单一的解决方案或策略，涌现的数量会随着时间的推移而减少。本质上，这发生在当问题有一个静态的全局解决方案时，而这个行为是用来解决这个问题的。考虑井字游戏:虽然采用强化学习或遗传算法等新兴技术的人工智能最初可能会在井字游戏中显示各种各样的策略，但它会相当快地收敛到游戏的单一主导解决方案。在这一点上，不会出现进一步的涌现或策略。

当我谈到突现行为时，我相信当人们接近这个话题时，他们通常想到的是“开放式突现”。在这种形式的涌现中，没有固定的解决方案，一个主体会不断产生新的行为。创造这种开放式涌现比创造固定点涌现更复杂，你必须考虑创造它所必需的结构。令人欣慰的是，最近在理解和编纂开放式涌现的要求方面取得了一些进展，特别是 Joel Liebo 和其他人在自动课程方面所做的工作。在一篇开创性的论文“自主课程和来自社会互动的创新的出现”中，作者列出了一些可以预期开放式涌现发生的条件。下面的表格将有助于理解什么时候会出现无止境的涌现。

开放式涌现的秘诀

1.其中，行为或目标的解决方案非常大，因此，对于所有的意图和目的，涌现的数量在检查的时间周期内看不到顶端。这不是真正的开放式涌现，因为可能存在全局解决方案，如果给予足够的时间，代理将发现并停止适应。然而在许多情况下，这可能超出了人类生命的时间跨度，所以从观察者的角度来看，这种出现将是无限的。一个例子可能是国际象棋游戏，其中一个全局解决方案被认为存在于井字游戏中，但由于游戏的复杂性，这还没有被发现，而且在我们的有生之年也不太可能。

2.创造开放式涌现的另一个方法是采用依赖于不断变化的环境的行为或目标状态。想想地球上的碳基生命形式——由于行星气候条件，环境在不断变化，这确保了动物必须始终适应才能生存，并且在出现水平上没有上限。

3.开放式涌现的第三种方法是多智能体场景，其中涉及合作或竞争，智能体采用一些适应性学习策略。在这种情况下，一个主体必须适应其同伴或竞争对手的策略，这反过来确保其他主体也必须适应，从而形成一个不断适应的反馈循环，这是一种进化军备竞赛。虽然这可能导致涌现，但也可能导致重复行为的循环，无限循环，直到环境中的某些因素将它们推出循环。这不是一个真正的均衡，因为行为或策略不是固定的，但循环本身成为一种均衡。确保环境足够复杂和动态，这是避免这些周期性行为或策略的一种方式。

展示开放式涌现

在下面用虚幻引擎创建的例子中，我选择使用上面列出的第三个配方，一个涉及两个虚拟物种的多代理场景。我们称其中一只为犀牛，另一只为老虎。两组人都在探索他们的环境，寻找“浆果”，在虚幻引擎游戏环境中，这些浆果以大型圆形物体的形式出现。让事情变得更复杂的是，这些浆果丛中随机分布着野生鸟类，它们会尖叫，吓跑老虎和犀牛。然而，如果一只老虎或犀牛与另一只老虎或犀牛一起接近浆果，它们会吓得鸟儿们不敢吱声。他们现在可以吃浆果了，尽管他们必须在彼此之间把浆果分开。

除了这种合作行为，它们还有一种竞争行为，在这种行为中，它们可以发出吼声(老虎)或开始跺脚(犀牛)，吓跑对手物种的代理人，但也消耗自己宝贵的能量。在许多方面，这种情况复制了自然界中动物和人类企业的一些合作和竞争权衡。

竞争或合作是最古老的问题之一，通过在虚拟绅士中模拟这一点，我们可以获得对导致各种竞争或合作结果的条件的宝贵见解。使用这样的虚拟代理，我们也可以创造出在动物王国甚至人类身上可以观察到的行为复杂性。OpenAI group 最近发表的一篇论文能够证明工具使用的出现是玩捉迷藏游戏的 AI 角色之间多智能体竞争的函数。

类似的事情也可以在任何一群适应性虚拟代理人身上发生，这些代理人会在彼此之间玩一个游戏，在这个游戏中，收益取决于其他玩家的策略，而环境会受到他们行为的影响。目标可以是任何东西，只是必须有一些方法让代理接收关于他们的环境的反馈，包括竞争者采取的行动和作为响应改变他们自己的行动的能力。这个目标是由玩家设定的还是由程序员设定的，从复杂程度上来说没有太大的区别。选择由玩家来设置它可以带来关于游戏叙事的有趣选项。需要注意的重要一点是——谁选择行为的目标并不重要，重要的是这符合上述竞争或合作的要求。

但是涌现实际上是如何“涌现”的呢？一种方法是使用强化学习或遗传算法。在我们的例子中，我们有一群虚拟代理人，他们随机改变行为来观察目标状态是否受到影响。代理优先考虑那些带来好结果的行为。在遗传算法中，这是通过一种在群体中发生的交叉函数来实现的，在这种情况下，随机发生的良好适应被保留，而相对不那么有用的适应被丢弃。这是由适应度函数决定的。这不需要看起来像有性生殖或任何类似的东西，因为现有的代理可以无缝地被它们适应的对应物取代，而用户看不到任何正在发生的交叉。在我们的例子中，我们将使用强化学习，这是一种算法技术，我已经在另一系列文章中广泛介绍过。

与遗传算法类似，强化学习依赖于大数定律，即给定足够多的随机行为，好的行动将会出现，并可以在未来进行优先排序。最初，代理人选择一系列随机的行为，但是如果这些行为导致代理人获得奖励，那么奖励的价值就归因于代理人采取的行为，因此在未来更有可能重复这些行为。当这个随机行为的过程被充分重复时，那些与接受奖励的代理人有因果关系的行为将与那些仅仅是偶然事件的行为区分开来。这就是为什么强化学习是核心，一种因果算法并可用于检测因果关系。

人们可以这样想象这个过程。

作者图片

为了这个项目，我使用了免费的思维制造插件和稳定基线强化学习算法套件。通过 MindMaker，我们可以使用 OpenAI Gym 格式在虚幻引擎游戏环境中轻松部署各种 RL 算法。这为在各种模拟环境中部署强化学习算法提供了一个通用的结构。

总结和未来发展

上面描述的简单设置导致一种情况，其中差异数字物种的行为不可能达到稳定的平衡。鸟类在浆果斑块中的随机分布确保了环境总是在变化，而这反过来又不断改变着代理商合作或竞争的潜在利益。代理人本身的活动也会影响有鸟类和没有鸟类的浆果斑块的比例——随着更多的代理人相互合作，合作的动力就会减少，因为没有鸟类的浆果斑块的比例会增加。

这个想法是，这些相互作用的力量创造了一个总是波动和不稳定的环境，推动着战略适应。鉴于我们的代理人与他们的环境互动的方式很少——他们不能利用环境中的物体诱捕竞争对手或隐藏浆果，我们的范式不太可能导致像紧急工具使用这样有趣的事情。然而，它确实避免了固定的均衡，并创造了不断变化的合作和竞争策略。这可以提供比开放世界视频游戏中通常出现的更有趣的一组行为。因此，我相信，像这里概述的这种新兴行为技术很可能会主导下一代视频游戏人工智能，创造出更有趣、更吸引人的角色和行为剧目。

Aaron Krumins 是“智胜——强化学习的承诺和危险”的作者，拥有 web 应用和机器学习的背景。他目前是一名自由职业的机器学习顾问。

使用 Python 创建 Excel 报表

原文：https://towardsdatascience/creating-excel-reports-with-python-e8c94072f09d

如何在不打开 Excel 的情况下创建 Excel 文件

照片由米卡·鲍梅斯特在 Unsplash 上拍摄

Excel 非常受欢迎。

从 Investopedia 来看，“在商业中，从字面上看，任何行业的任何职能都可以受益于那些拥有强大 Excel 知识的人。Excel 是一个强大的工具，已经成为全球业务流程中不可或缺的一部分”

业务用户喜欢 Excel，因此大量的报告需要在 Excel 中。

但这并不意味着我们必须使用 Excel 来处理数据。

在这里，我们将学习使用 Python 和 Pandas 创建 Excel 报表，而无需打开 Excel。

背景

作为一名数据分析师，我的第一个任务是在 Excel 中运行生产力报告。运行这个报告需要我手动下载两个大型数据集，将数据导入 Excel，合并数据，并执行计算。

因为数据集很大——至少对 Excel 来说——我经历了持续的滞后和可怕的死亡之轮。

它花了整整 8 个小时的工作日才完成。

然后，在我可以(当然是手动的)分发之前，它必须被发送给一个主管进行双重检查。

该报告必须每两周用新数据重新创建一次。在一年的时间里，仅复制一份报告就要花费超过 5 周的时间！

我知道一定有更好的方法。

任务

在这里，我们将使用 Python 为我们想象的小部件工厂创建一个 Excel 报表。

widget factory 委托我们创建一份报告，详细说明工人和团队创建 widget 的平均时间。他们希望报告是一个 Excel 文档，每个团队都有一张表格。

资料组

我们将继续使用我们在本文的中创建的小部件工厂数据集，在那里我们通过 Python 生成假数据。

创建报告

我们所有的数据操作都是使用 Python 和 Pandas 完成的。

首先，我们使用 Pandas 读入数据集，并检查每个数据集的形状。

接下来，我们使用.agg()通过worker_id查找汇总统计数据。对数据进行分组后，我们就有了一个需要扁平化的具有多索引的分组数据集。

我们还想处理一些列标题的格式。当我们写到 Excel 时，我们希望标题是漂亮的，所以我们现在可以处理它了。

最后，我们将分组的数据集合并回工人数据集，以获得工人姓名和团队。

按员工分类的汇总统计数据—按作者分类的图片

我们将数据集形状打印为双重检查，以确保我们对数据集的操作和合并是正确的。我们可以看到我们有 5000 行，并且确信这是正确的，因为 worker 数据集也有 5000 行。

下一步是为每个团队将数据分割成单独的数据帧。

首先，我们需要一个独特的团队名称列表。我们可以使用unique()来寻找唯一的名字，而不是硬编码这些名字。这将有助于防止将来数据发生变化——我们的代码不会崩溃！

接下来，我们创建一个空字典，用于存储团队数据帧。

最后，我们遍历各个团队，为每个团队创建一个新的数据帧，并将其存储在字典中，关键字是团队名称。

我们通过打印每个团队数据帧的形状来仔细检查这是否如预期的那样工作。正如所料，这些行加起来有 5000 行！

现在，我们准备将最终报告写入 Excel。

写入 Excel

我们将使用 Pandas 和 XlsxWriter 创建我们的 Excel 报表。

XlsxWriter 是熊猫用来编写 Excel 的 XLSX 格式文件的默认模块。

可以使用以下命令安装它:

pip install xlsxwriter

首先，我们创建一个变量来存储文件名和路径。

接下来，我们使用 Pandas 的ExcelWriter来编写我们的 Excel 文件。

熊猫的文档指示:“作者应该被用作上下文管理器。”非常简单，这告诉我们应该将代码包装在一个with块中，以确保我们的外部资源(这里是我们正在创建的 Excel 文件)在我们完成后关闭。

最后，遍历我们的数据帧字典，将每个数据帧写到它自己的表中。

我们现在可以在 Excel 中打开完成的报告了！

已完成的 Excel 报告—作者图片

虽然我们已经成功地将报告写到 Excel 中，并为每个团队提供了单独的表格，但它看起来并不漂亮。

让我们添加一些格式来改善该报告的外观:

我们为每列添加了下拉箭头，冻结了顶行，并为每列的宽度添加了适当的间距。

带格式的 Excel 报表-按作者排序的图像

这样看起来好多了！我们的报告现在可以分发了！

结论

在这里，我们学习了使用 Python、Pandas 和 XlsxWriter 创建 Excel 报表。

虽然我们只触及了 XlsxWriter 功能的表面，但我们能够使用 Python 创建 Excel 报表，而无需打开 Excel！

为训练机器学习模型创建基础设施

原文：https://towardsdatascience/creating-infrastructure-for-training-machine-learning-models-aec52cf37929

引入用于训练、跟踪和比较机器学习实验的自动管道

Firmbee在 Unsplash 上拍照

让我们想象一下下面的场景:你有一个新项目要做。对于这个项目，你需要开发一个机器学习模型，这将需要运行和训练几个实验。每个实验可能需要几个小时甚至几天，并且需要跟踪。

在开发阶段，你有自己的笔记本电脑，但是用自己的笔记本电脑进行培训和运行所有的实验是不现实的。首先，您的计算机可能没有所需的硬件，例如 GPU。第二，浪费时间。如果你可以并行运行和训练所有的实验，为什么要在一台计算机上顺序训练？

如果你平行训练，每次实验在不同的机器上进行，你如何容易地跟踪和比较他们？你如何访问训练数据？你能实时意识到任何故障以便你能够修复并再次运行它吗？

所有这些问题都将得到解答。在这篇博文中，我想描述一个解决上述所有问题的可能方案。

该解决方案结合了不同的服务，并将它们联合起来，以创建一个用于训练机器学习模型的完整管道解决方案。

在下图中，您可以找到该管道的方案。

作者图片

流水线由三个步骤组成——将代码容器化，以便我们以后能够运行它，在几台独立的机器上执行代码，以及跟踪所有的实验。

(1)集装箱化

大规模运行训练实验的最佳方式是构建一个 docker 映像。创建这种图像的一种方法是使用 詹金斯 。开源自动化服务器，支持构建、测试和部署软件。另一种方法是使用 GitHub 动作 **。**持续集成和持续交付(CI/CD)平台，允许您自动化您的构建、测试和部署管道。这可以直接从你的 GitHub 账户上完成。

让我描述一个可能的场景——每个对我们的 git 分支的推送都会触发一个 Jenkins 管道。这个管道构建 docker 映像(除了它能做的其他事情之外)，标记它，并将其存储在 docker 注册表中。

一旦我们有了这个图像，我们就可以执行它，并使用虚拟云机器大规模地训练我们的模型。

(二)执行

我们已经准备好代码，我们想开始训练。假设我们有 10 个实验要运行，我们希望并行运行它们。一种选择是使用 10 个虚拟机，使用 ssh 连接每个虚拟机，并在每个虚拟机上手动执行我们的代码。我想你明白这不是最优的。

由于我们已经在上一步中容器化了我们的应用程序，我们可以利用最流行的容器编排平台来运行我们的培训— Kubernetes 。Kubernetes 是一个开源系统，用于自动化容器化应用程序的部署、扩展和管理。这意味着，通过使用它，我们可以获得 docker 映像并大规模部署它。我们可以使用一个简单的 bash 脚本在 10 台独立的机器上运行 10 个实验。

你可以在这里阅读更多相关信息。

(3)跟踪

训练数据采集

为了执行培训，我们首先需要访问我们的培训数据。最好的方法是从训练机器直接访问存储器中所需的训练文件( S3 、 GCS 等)以及数据库中所需的表格信息。为了确定这一点，我们可能需要帮助来访问相关的凭证和秘密管理器(一个用于 API 密钥、密码、证书和其他敏感数据的安全而方便的存储系统)。一旦配置完成，我们就可以开始了。

NFS 磁盘—共享磁盘

用于培训的所有机器都连接到一个 NFS 磁盘(在我们的例子中，我们使用 GCP 的托管解决方案，称为文件存储)，这是一个额外的共享磁盘，所有机器都可以对其进行读写访问。我们添加这个共享磁盘有三个主要原因—

1. 当在不同的机器上训练时，我们不想为每台机器反复下载我们的训练数据。相反，我们可以将所有机器连接到一个共享文件存储，这是一个所有机器都可以访问的附加磁盘。这样我们可以只下载一次数据。
2. 培训模型可能需要时间，因此成本可能很高。让训练更便宜的一个选择是使用 现货实例 ，价格低得多。但是，如果需要回收计算容量以分配给其他虚拟机，您的云提供商可能会停止(抢占)这些实例。使用 Kubernetes，一旦实例停止，它会被自动替换并再次执行相同的实验。在这种情况下，我们希望从停止加载最后一个保存的检查点的地方继续训练。只有使用共享磁盘才能做到这一点，没有共享磁盘，我们将无法访问先前机器中保存的检查点。其他解决方案，如将所有检查点上传到云存储并在其上应用逻辑，也是可行的，但开销(和成本)很高。
3. TensorBoard —使用 TensorBoard 比较不同机器的不同实验，所有信息必须写入同一父目录。我们将把所有信息写入同一个父日志目录，TensorBoard UI 将从中读取。

实验跟踪

训练模型中最重要的一点是能够跟踪训练进度并比较不同的实验。对于每个实验，我们希望保存我们使用的参数、我们训练的数据、损失值、结果等等。我们将使用这些数据进行跟踪、追踪和比较。

有许多跟踪工具。当我想要高级可视化时，我使用 MLflow 和 TensorBoard 。它们都是以这样一种方式部署的，即所有团队都可以访问同一个实例，并看到团队的所有实验。

MLflow 是一个管理端到端机器学习生命周期的开源平台。它提供了一个简单的 API 来集成，并提供了一个很好的 UI 来跟踪和比较结果。

TensorBoard 提供了机器学习实验所需的可视化和工具。它提供了一个 API，将我们的信息记录到机器磁盘上的文件中。TensorBoard UI 从传递的输入目录中读取这些文件。

错误警报

由于培训可能需要时间、几个小时甚至几天，因此如果出现错误，我们希望得到实时通知。一种方法是将我们的代码与 Sentry 集成，将 Sentry 与 Slack 集成。这样，代码中出现的每个错误都将作为事件发送给 Sentry，并触发 Slack 通知。

现在，我可以实时收到错误通知。这不仅减少了我监视运行的需要(不管它们是崩溃了还是还在运行)，而且还减少了错误发生和触发下一次运行之间的时间，代码是固定的。

摘要

当我第一次训练模型时，是在我的笔记本电脑上完成的。很快，它就不可伸缩了，所以我使用 ssh 转移到虚拟机上。这种选择也有缺陷。对于每台机器，我不得不登录、提取代码、下载数据、手动运行代码，并希望我没有混淆我传递的参数。数据部分是最容易解决的，我们为我的所有虚拟机添加了一个共享磁盘。但我仍然对参数感到困惑，并在几台机器上传递了相同的参数。**正如我在这篇文章中所描述的，是时候使用可扩展的、更加自动化的基础设施了。**与 DevOps 团队一起，我们定义并构建了我们今天拥有的强大渠道。现在，我有了一个可伸缩、易于执行的管道，甚至在每次出现错误时都发送一个 Slack 通知。这种流水线不仅由于实验的并行性而节省了时间，而且最大限度地减少了参数错误和由于错误而导致的空闲时间，除非主动检查每台机器，否则我们不会知道这些错误。

https://www.docker/ https://www.jenkins.io/ https://kubernetes.io/ https://mlflow/ https://www.tensorflow/tensorboard https://cloud.google/ https://sentry.io/welcome/ https://docs.github/en/actions

使用 Python 创建交互式地理空间可视化

原文：https://towardsdatascience/creating-interactive-geospatial-visualisation-75a2dfaf7e59

从提取几何数据到创建交互式图表—逐步指南

德尔菲·德拉鲁阿在 Unsplash 上拍摄的照片

背景

最近，在澳大利亚的新南威尔士州，许多地理位置被政府宣布为新冠肺炎传播的关注区域，与其他区域相比，记录了相对更多的病例。

当时(大约在 2021 年 7 月)，澳大利亚针对新冠肺炎采取了一种传播最小化策略。因此，对这些受关注地区的居民实施了更严格的公共卫生命令和限制(“封锁措施”)，这意味着更少的人被允许离开家园进行非必要的活动(如外出就餐)，提供非必要服务的当地企业(如咖啡馆和餐馆)被迫暂时关闭。这影响了消费，并对人们的收入产生了流动效应。请注意，这也可以推广到其他有严格封锁措施的地区，如 2021 年 9 月的整个维多利亚州，或 2021 年的一些欧洲国家。

从人寿保险公司(本文作者是该公司的精算师)的角度来看，锁定措施可能会间接导致人寿保险保单的投保人因家庭可支配收入减少而停止向人寿保险公司续保(“终止”)。

为此，我的公司委派我了解封锁措施是否对澳大利亚某些地理位置的中断率有影响。要做到这一点，没有比在地图上可视化地理区域的中断率更好的方法了，这是我意识到对一个人来说，实施提取几何数据的工作流程，用正确的软件包在 Pyhon 中应用可视化，并最终创建一个综合的同时显示所需信息的交互式图表是多么困难。

鉴于上述动机，本文提供了一步一步的指南，以创建一个交互式图表，显示在悉尼按地理区域的中断率。虽然这个用例是用澳大利亚的几何数据呈现的，但它可以很容易地扩展到世界上的其他地理区域。

中止及其重要性

在澳大利亚，零售市场上的寿险保单主要是由我们称之为金融顾问的中介机构销售的。金融顾问会因成功销售而获得初始佣金。

大多数保单都是以高额初始佣金售出的。在 2018 年之前，这可能会远远超过首年保费的 110%。这种佣金结构通常被称为“预付”结构，这意味着在保单有效期内支付的大部分佣金在保单生效时支付。

如果投保人在出售后不久(比如一两年内)终止保单，人寿保险公司将遭受损失，因为它无法收回支付给财务顾问的初始佣金。事实上，人寿保险公司通常需要 5-7 年的时间才能从保单中赚钱。这就是为什么终止保险很重要，而且大多数公司都投入大量资金来留住现有的投保人。

提取几何数据

澳大利亚地图可以通过各种数字边界划分成地理区域。一个例子是将澳大利亚地图按州界分成九个州。同样，澳大利亚的地图可以用其他类型的边界划分成更多的自然区域。下图分别显示了以州界和地方政府区域(“ LGA ”)为界划分的澳洲地图。

图 1:澳大利亚的州界。图片作者。

图 2:LGA 边界旁的澳洲。图片作者。

前面提到的宣布为关注区的地理位置是在 LGA 一级宣布的。对于悉尼的居民来说，你可能熟悉下面显示的一些地方政府。

表 3: LGA 的例子。按作者分类的表格

几何数据通常以 json 或 shapefile 格式提供。这些都是由澳大利亚政府机构发布的。出于本文的目的，我将提取新南威尔士的 LGA 的几何数据，这些数据可通过 LGA 数字边界在这里获得。

可以使用 Python geopandas 包读取几何数据，如下所示:

## path_LGA_geo refers to the path to the json file downloaded 
## from the link abovedf_LGA_geo = geopandas.read_file(path_LGA_geo) ## In practice, we'll only need two fields from the dataframe,  
## namely the locality name, denoted by 'nsw_lga__3', and the 
## geometry arrays, denoted by 'geometry'df_LGA_geo_final = df_LGA_geo[['nsw_lga__3','geometry']].drop_duplicates(subset ="nsw_lga__3",keep = False, inplace = False)

如果你没有太多的计算资源，你可以按照这个链接中的说明选择缩小 json 或 shapefile。

按地理区域提取中断率

在高层次上，中止率可以被校准为中止的保单数量与有效保单总数之比。由于我们跟踪已中止和有效的投保人，中止率可以在邮政区域级别(即通过邮政编码)进行校准。

然而，与 LGA 相比，邮政区域代表了不同类型的数字边界。幸运的是，澳大利亚统计局在这本出版物中提供了邮政区域和 LGA 之间的地图。通过下载 ASGS 通信(2016)，可以在“CA _ POSTCODE _ 2018 _ LGA _ 2018 . xlsx”文件中找到该映射。

我顺便提一下，同一出版物中还提供了数字边界的其他映射，如统计区域到 LGA、统计区域到地区(即郊区)等。这为用户提供了在不同粒度级别切割几何数据的选项。

一旦邮政编码被映射到新南威尔士州的 LGA，就可以通过地方政府机构汇总中止的保单数量和有效的保单数量，并通过地方政府机构校准中止率。然后，这些可以以表格格式存储，如下所示。请注意，出于本文的目的，显示的中断率是合成的(尽管受到实际范围的限制)。

表 4:示例数据帧。按作者分类的表格

然后，可以使用下面的 Python 代码将数据帧与几何数据帧连接起来:

## Joining geometry data with discontinuance rate data
## df_LGA_lapse is a dataframe containing discontinuance rate by 
## LGAsdf_geo_disc_merge = df_LGA_geo_final.merge(df_LGA_lapse,left_on = 'nsw_lga__3', right_on = 'LGA_Name')

创建可视化

我将使用 Python 包来创建可视化。关于这个包的文档使用，请参考这个 Github 页面。

还有其他可视化工作所需的包。如果你正在使用 Google Colab，你可以安装和导入如下的包。

图 5:必需的包。作者图片

然后可以通过以下四(4)个步骤创建交互式可视化。

第一步:获取中心坐标

这是为了确保观想“打开”在被绘制的集体地理区域的中心，在这种情况下，是新南威尔士州的中心，也就是澳大利亚地图右下角附近的州。

# Obtain center cordinates of the mapx_center = df_geo_disc_merge.centroid.x.mean()y_center = df_geo_disc_merge.centroid.y.mean()

第二步:创建地图的基础图层

## Create the base layernsw_disc_map = folium.Map(location = [y_center, x_center], zoom_start = 6, tiles = None)folium.TileLayer('OpenStreetMap', name = "Actl Map", control = True).add_to(nsw_disc_map) nsw_disc_map

上面的 Python 代码创建了地图的基础图层，该图层在中心坐标处打开。这个基础图层有一个“OpenStreetMap”的图块样式，如下图所示。您可以根据 Folium 文档在这里使用样式选项定制图块集(并搜索图块图层类)。

图 6:树叶地图的基础层。作者图片

第三步:将 LGA 的中止率添加到基础层

可以使用下面的 Python 代码将 LGA 的中断率图层添加到地图的基础图层上。

set_quantile = (df_geo_disc_merge['Discontinuance_Rate'].quantile((0.0, 0.2, 0.4, 0.6, 0.8, 1.0))).tolist()nsw_disc_map.choropleth(geo_data = df_LGA_geo_final,name = 'Choropleth',data = df_geo_disc_merge,columns=['LGA_Name','Discontinuance_Rate'],key_on = 'feature.properties.nsw_lga__3',fill_color = 'YlGn',fill_opacity = 0.6,line_opacity = 0.3,threshold_scale = set_quantile,legend_name = 'Discontinuance Rate by LGA',smooth_factor = 0,highlight = True) nsw_disc_map

如下图所示，在添加的图层中引入了一些不错的功能。即:

• LGA 的数字边界已经在新南威尔士州绘制出来，每个 LGA 都由一块以绿色突出显示的土地表示。
• 中断率以 5 个分位数的色标显示(代码中有 set_quantile 自定义)，由地图右上角显示的不同绿色色标表示。

图 7:地图中 LGA 图层的中断率。作者图片

步骤 4:创建工具提示交互

让我们通过添加工具提示层来使地图具有交互性。想法是，当您将鼠标光标悬停在地图上的特定 LGA 上时，将会出现一个工具提示，并显示有关 LGA 的某些信息，在本例中为中断率。

这可以通过以下 Python 代码来实现。

## Add a tooltip layer add_tooltip = folium.features.GeoJson(df_geo_disc_merge, tooltip = folium.features.GeoJsonTooltip(fields=['LGA_Name', 'Discontinuance_Rate'],aliases = ['LGA_Name: ' ,'Disc_Rate: '],))nsw_disc_map.add_child(add_tooltip)nsw_disc_map.keep_in_front(add_tooltip)nsw_disc_map

下面的短片展示了本练习的最终结果。

视频 8:交互式地图演示。作者提供的视频

您可以通过自定义样式和高亮功能来进一步设置地图的样式，这些功能可以嵌入到工具提示层中。互动 Choropleth 地图教程是一个学习的好地方。

工具提示上有更多！

实际上，您可以在工具提示上显示更多信息，只要您可以在与地图相同的位置级别上收集信息(在本例中为 LGA)。例如，我在工具提示中添加了收入数据，如 LGA 的收入者的计数、中值年龄、平均值和平均收入，如下图所示。收入数据可以从这里获得。

此外，我设置的不透明度让我可以看到地图顶层下面的郊区，这是另一个很好的功能。

图 9:带有输入数据工具提示。作者图片

作为一名精算师，可视化之后可以进行一些按地理区域划分的收入与中止率的分析，这可能证明是有见地的。

我能想到的其他使用案例包括:

• 将此与自然风险数据(如干旱、森林火灾、地震或与洪水相关的地区)叠加，以确定公司的投资组合受影响的程度，并告知公司是否要向受影响的投保人提供保费减免。
• 通过地理区域(可能还有财务顾问)可视化销售数据来发现机会。

结论

我将几何数据视为半非结构化数据的一种形式，与文本和图像等其他非结构化数据一样，对其进行导航并最终将其可视化在地图上可能很困难。

在本文中，我提供了一个逐步指南，介绍如何使用 Python follow 包创建交互式地图。请注意，虽然本文中讨论的用例是在澳大利亚的环境中，但是您可以很容易地将其扩展到其他管辖区域，就像替换澳大利亚的几何数据一样简单。

对于对保险行业的其他数据科学应用感兴趣的读者，请关注我，并愉快地阅读我在媒体上的文章。

参考

[1]澳大利亚统计局(2021) 新南威尔士州地方政府区域—地理景观行政边界 ， ABS 网站，2022 年 1 月 14 日访问，(在https://creativecommons/licenses/by/4.0/获得许可)

[2]澳大利亚统计局(2021)【ASGS 地理通讯(2016) ， ABS 网站，2022 年 1 月 14 日访问，(在 https://creativecommons/licenses/by/4.0/获得许可)

[3]澳大利亚统计局(2021) 澳大利亚个人收入 ， ABS 网站，2022 年 1 月 15 日访问，(https://creativecommons/licenses/by/4.0/许可)

创建电影分级模型第 3 部分:测试候选模型算法

原文：https://towardsdatascience/creating-movie-rating-model-part-3-testing-out-model-algorithm-candidates-612c7cf480ce

尝试五种二元分类和回归算法来支持我们的电影分级模型项目

朋友们好！是的，我们已经连续第二天在我们的电影分级模型系列中有了新的条目。正如我在昨天的帖子中所说的，我仍然在以博客的形式玩着我迄今为止作为 YouTube 直播系列所做的事情。在那篇文章中，我们介绍了从原始系列文章中收集的原始数据集，并执行了一个特征工程来筛选一些干净的特征。如果你想进一步了解，你可以在这个 GitHub 库中找到所有这些工作和支持数据。GitHub 项目的自述文件解释了这个项目的所有内容，所以我建议你去看看，而不是在这里重复。(我想我刚刚创下了一个段落中超链接数量的个人记录！)

既然我们的数据已经进行了适当的特征工程，我们准备开始测试一些不同的机器学习模型算法候选，以查看哪一个将在本系列的下一篇文章中最适合创建正式训练的模型。请记住，我们实际上将为这个项目创建两个模型。这是因为电影评论家卡兰·比恩(Caelan Biehn)为他的电影提供了两个评分。第一个分数更像是一个“是或否”的二元评级，第二个分数是一个介于 0 和 10 之间的数字分数。后者被称为“比恩等级”，可以用一个小数位来表示，例如，电影可以得到 4.2 或 6.9 的评级。因为有两个不同的分数，我们将为二进制“是或否”分数训练一个二进制分类 算法，为 Biehn 标度分数训练一个回归算法。

因为我们希望创建尽可能精确的模型，所以我们将为二元分类和回归模型尝试五种不同的算法。在我们开始测试各种算法之前，让我们为如何继续这篇文章制定建模策略。

建模策略

虽然我们已经注意到，我们将测试每种算法中的五种，但是在执行建模时，我们还需要做一些特定的活动。这些事情包括以下内容:

• 超参数调整:为了确保每个算法都以最佳状态运行，我们将执行超参数调整，为每个模型寻找理想的超参数。具体来说，我们将使用 Scikit-Learn 的 GridSearchCV 工具来帮助我们。
• K 倍验证:因为我们将要训练的数据集相对较小(大约 125 条记录)，所以我们不能像对待通常较大的数据集那样进行典型的训练测试分割。因为我们想最有效地利用我们的数据集，所以我们将使用 Scikit-Learn 的 k 倍验证。这个过程会将数据集分成几个小的训练和验证批次，这将发生多次。该过程的输出将允许我们最大程度地评估数据集。(如果你想了解更多关于 k 倍验证的知识，我最近发表了一篇不同的文章正是关于这个话题！)
• 指标验证:对于经过训练的模型，我们希望通过将它们与适当的验证指标进行比较，来确保它们能够有效地执行。我们将在下一节详细讨论验证指标。
• 特征缩放(可选):根据我们使用的算法，我们可能需要也可能不需要对数据集执行特征缩放。为了执行这种缩放，我们将使用 Scikit-Learn 的 StandardScaler 对数据执行基本的标准化。

现在，我们已经在一个高层次上制定了我们的策略，让我们更多地讨论验证指标，然后继续创建一些帮助函数，使测试五个算法中的每一个都变得非常容易。

模型验证指标

确保任何模型按预期执行的基础是，我们必须为每种不同类型的算法定义一些模型验证度量。对于二元分类和回归模型，我们将使用三种不同的度量来验证每个模型的性能。所有这些验证指标都将在 Scikit-Learn 内置函数的帮助下进行计算。

对于我们的二元分类算法，下面是我们将测试的模型验证指标:

• 准确性:在这个列表的所有指标中，这个指标显然是最基本的。归根结底，这是我最关心的指标，但是这个验证指标本身并不能说明全部情况。为了全面了解模型是否正常运行，我们需要其他验证指标的帮助。
• ROC AUC :这个指标的全称其实是“受试者工作特征曲线/曲线下面积”自然，它经常被缩写为 ROC AUC，因为全称很拗口。这个指标的作用是计算算法在一条曲线上多个不同阈值下的性能，我们可以从技术上把这条曲线画在图上。一旦计算出这条曲线，我们就可以用曲线下面积(AUC)来为这个模型的表现给出一个分数。
• F1 :如上所述，准确性本身并没有多大帮助。我可以做的是计算模型的精确度和召回率，或者我可以通过计算 F1 分数让生活变得简单一点。F1 评分发现了精确性和回忆性评分之间的和谐，从而产生单一评分。

对于回归算法，我们将计算以下三个验证指标:

• 平均绝对误差(MAE) :对于每个数据点，此验证度量计算观察点离模型生成的回归线有多远。回归线和观察点之间的距离通常被称为误差，因为一些观察点可能低于回归线，这意味着我们可能会有负误差。为了抑制这种情况，我们取误差的绝对值，然后取所有这些绝对值的平均值来产生我们的单个 MAE 分数。
• 均方根误差(RMSE) :这个验证指标与上面的非常相似，除了我们不是取误差的绝对值，而是取误差的平方。这些误差的平均值也将被平方，如果我们愿意，我们可以停下来计算均方误差(MSE)。因为我不是特别喜欢 MSE 不代表实际误差，所以我要取 MSE 的平方根，这自然会产生 RMSE。
• R 平方:也称为“决定系数”，这个验证度量计算模型输出相对于模型输入的可解释性。值为 1 的 R 平方得分意味着模型输入完美地解释了模型输出，而接近 0 的数字意味着模型输入不能很好地解释输出。是的，R 平方度量可以是负的，这意味着当试图解释模型输出时，模型输入是一种倒退。请记住这一点，因为这个特定的指标将在这篇文章的结尾非常有说服力。

助手功能

因为我们将测试许多不同的算法，所以我将创建两个可重用的辅助函数:一个用于二进制分类算法，另一个用于回归算法。这些函数将做一些事情，如特征缩放(根据需要)，执行超参数调整，用理想的超参数训练模型，用 k-fold 验证来验证训练的模型，并打印出来自每个 k-fold 训练的平均验证度量。

(因为这些是大块的代码，这里有一个友好的提醒，你可以在这个 Jupyter 笔记本中找到这些相同的功能。)

以下是二进制分类算法的辅助函数:

作者创建的图像和代码

这是回归算法的辅助函数:

作者创建的图像和代码

唷！这些是一些很大的函数，但是好消息是，你将会看到这两个函数将会使测试每一个算法变得多么容易。说到这里，让我们继续讨论我们将要测试的不同算法。

候选模型算法

如上所述，我们将为每种不同类型的模型测试五种不同的算法。让我们简单地谈谈每个候选人将要参加的选拔。

(M1 Mac 用户请注意:对于这两种类型的模型，我想尝试一下 Catboost 库中的一些算法，因为我的同事已经发现了这些算法的许多成功之处。不幸的是，Catboost 库还不能在基于 M1 的 Mac 电脑上运行。因为我在 M1 Mac mini 和微软 Surface Pro 8 之间工作，所以我必须在我的 Surface 上测试 Catboost 算法。对不起，M1 用户！)

二元分类候选

• Scikit-Learn 的逻辑回归算法:虽然名字中的“回归”可能具有欺骗性，但逻辑回归是一种非常简单但功能强大的二元分类算法。因为我们想要测试各种算法类型，所以我们选择 Scikit-Learn 的逻辑回归算法作为更简单的变体。
• Scikit-Learn 的高斯朴素贝叶斯(GaussianNB)算法:朴素贝叶斯算法最流行的实现，我们将测试 Scikit-Learn 的 GaussianNB 实现，看看它在我们的数据集上表现如何。
• Scikit-Learn 的支持向量机(SVM)算法:虽然不像逻辑回归算法那么简单，但 SVM 是一种更简单的算法。这种算法往往在更高维度(即具有更多特征的数据集)中表现更好，尽管我们的数据集维数更少，但我仍然认为它值得一试。
• Scikit-Learn 的随机森林分类器算法:这是 ML 行业中最流行的二进制分类算法之一。这是因为它经常产生非常准确的结果，并且具有更容易的算法解释能力。随机森林分类器也是被称为集合模型的经典例子。
• CatBoost 的 CatBoostClassifier 算法:你可能以前没有听说过这个算法，但它在我的财富 50 强公司的同事中非常流行。这是因为它经常被证明可以提供最佳的性能结果。因此，很自然地，我想看看它与我的数据相比如何！

回归算法候选

• Scikit-Learn 的线性回归:就像我们用二元分类器分析的逻辑回归算法一样，这可能是我们可以测试的回归算法的最简单实现。鉴于它的简单性，我显然不抱太大期望，但无论如何它总是值得一试！
• Scikit-Learn 的 Lasso 回归:这个算法和上面的算法是一个家族，Lasso 实际上是一个缩写，代表最小绝对选择收缩算子。坦白地说，我并不精通这个算法背后的数学，所以我甚至不打算解释它。😂
• Scikit-Learn 的支持向量回归机:类似于我们如何尝试二进制分类模型的支持向量分类器，我们将在这里看到支持向量回归机的表现。因为支持向量机通常计算特征之间的距离，所以这里的特征将需要被适当地缩放。
• Scikit-Learn 的随机森林回归:就像我们对二元分类模型使用随机森林分类器一样，我想在这里尝试一下它的回归变体，特别是因为这被认为是一个集合模型。
• CatBoost 的 CatBoostRegressor :最后，就像我们对前面提到的一些变体所做的那样，我想尝试这个最终算法作为第二个集合选项。

当我们使用上面的特殊助手函数时，我们可以非常简单地用几行简短的代码就能完成每一个候选函数。例如，下面是我需要为逻辑回归算法运行的所有附加代码:

*# Setting the hyperparameter grid for the Logistic Regression algorithm*
logistic_reg_params **=** {
    'penalty': ['l1', 'l2'],
    'C': np**.**logspace(**-**4, 4, 20),
    'solver': ['lbfgs', 'liblinear']
}*# Instantiating the Logistic Regression algorithm object*
logistic_reg_algorithm **=** LogisticRegression()*# Feeding the algorithm into the reusable binary classification function*
generate_binary_classification_model(X **=** X,
                                     y **=** y,
                                     model_algorithm **=** logistic_reg_algorithm,
                                     hyperparameters **=** logistic_reg_params)

看看助手函数是如何让测试所有这些算法变得如此容易的？我停在了五种不同的算法上，但是因为助手函数使它变得如此简单，如果我愿意，我可以非常容易地添加额外的候选算法。但是正如您将在我们的下一节中看到的，这并不是真正必需的。

模型验证结果

因为这个帖子已经有点过时了，所以我打算鼓励你去看看这个 Jupyter 笔记本看看具体的结果。此外，我不打算在这里进入每个算法的性能的本质细节，因为每个算法都或多或少地表现相同。当然，有些人比其他人过得好，但差别几乎可以忽略不计。让我们更多地讨论每种模型类型的验证结果。

二元分类验证结果

就准确性而言，每一个算法都有 78%或 79%的准确率。同样，每一个算法的 F1 分数也相当不错，要么是 87%，要么是 88%。ROC AUC 是事情变得有趣的地方，原因有二。首先，ROC AUC 得分在这里变化最大，与准确性和 F1 相反。算法之间的差异不是 1 点，而是 8 点。但是第二点和坏消息 ROC AUC 分数对每一个算法来说都很糟糕。它们在 50%和 58%之间。我并不特别惊讶，我会在即将到来的结论中解释我的不惊讶。

回归验证结果

哦，天啊…我们回归候选人的情况很糟糕。虽然与二元分类算法相比，我们有更多的结果分布，但结果并不太好。对于 MAE，我们看到的平均误差约为 1.6，对于 RMSE，我们看到的平均误差约为 2.0。请记住，我们讨论的是 0 到 10 分之间的分数范围，所以这个误差相当大。甚至谈论 R 平方分数都让我感到痛苦…记住，零或负的验证分数真的很糟糕，而且…嗯，我们看到几个候选人的分数是负的。在最好的情况下，一些算法几乎没有超过零。呀。

最终候选人选择和结论

好吧，伙计们…这并不像我希望的那样好！诚然，这只是一个帮助我们更好地了解数据科学概念的有趣项目，老实说，结果并没有让我感到惊讶。影评人卡兰·比恩(Caelan Biehn)并不是一个严肃的影评人，所以他给电影打出一些非常怪异的分数并不罕见。例如，他可能会给一部电影二进制“是的，去看这部电影”，但随后会给它一个像 1.3 这样极低的比恩评分。请记住，卡兰给出的评论是喜剧播客的一部分，所以荒谬有助于喜剧！此外，我正在处理的这个数据集非常小，目前只有大约 125 条记录(评论)。如果我有一个更大的数据集，也许我们会看到更强的信号。

那么，未来我们将采用什么样的候选算法呢？对于二进制分类模型，我们将使用 Scikit-Learn 的随机森林分类器算法。这是在这个和逻辑回归算法之间的一次掷硬币，但是在一天结束时，每次我重新运行网格搜索和模型训练/验证时，随机森林分类器仍然表现得更好一些。我们也可以使用 Catboost 分类器，但正如你在帖子前面看到的，显然 CatBoost 还不能在 M1 的 Mac 电脑上使用。鉴于我计划继续在我的 M1 Mac mini 上进行正式的模型训练和推理的直播，我不打算使用 Catboost，因为它根本无法在那台计算机上工作。

对于回归模型，我们将使用 Scikit-Learn 的 Lasso 回归算法。同样，我们看到所有回归模型的结果非常相似，但我们看到验证指标之间最大的差异是 R 平方得分。当然，所有的回归算法在 R 平方得分方面都表现得非常糟糕，但 Lasso 回归算法在我每次重新运行模型训练和验证时都表现得更好一些。(如果我完全诚实的话，我有点高兴 Lasso Regression 挤掉了他们，因为我有点希望那个算法“赢”，因为奇妙的电视节目 Ted Lasso 。😂)

这篇文章到此结束！即使我们的模型在以后的推理中不会很棒，我仍然认为这个项目为了学习的目的是值得继续的。在本系列的第 4 部分中，我们将创建一个正式的模型训练管道，它还捆绑了所有适当的特性工程，以将所有这些转换序列化为两个 pickle 文件。即使这里的结果很糟糕，我仍然希望你学到了一些东西，并且一路上玩得开心。下期帖子再见！

创建电影分级模型第 4 部分:创建完整的模型培训渠道

原文：https://towardsdatascience/creating-movie-rating-model-part-4-creating-a-full-model-training-pipeline-d52f54fd03c4

将我们的特征工程和模型训练代码正式化为单一管道

作者创作的标题卡

你好，朋友们！我们回到创建电影分级模型系列的第 4 部分。如果你想看这个系列之前的帖子，请查看这个链接，你也可以在我的个人 GitHub 库找到整个项目的代码。简单回顾一下我们停下来的地方，我们之前的帖子测试了许多不同的算法候选项，用作我们的最终模型。提醒一下，我们在技术上构建了两个模型:一个预测电影评论者 Caelan Biehn 给出的二元“是/否”电影评级，另一个预测评论者以基于回归的模型形式给出的 0-10 分。

在测试了许多不同的算法后，我们决定使用 Scikit-Learn 的 RandomForestClassifier 用于二元分类模型，使用 Scikit-Learn 的 Lasso 回归用于回归模型。现在，我们已经准备好创建一个正式的模型培训管道。现在，如果到目前为止您已经了解了这个系列，那么我们主要是在 Jupyter 笔记本上完成我们的工作。因为理想情况下，我们希望将这样的训练脚本集成到类似 MLOps 管道的东西中，所以最好将它转换为 Python 文件。

让我们假设，当我们在将来收集用于推理的数据时，在将数据输入模型之前，在运用特征工程来清理数据之前，它将看起来是相同的原始格式。如果是这种情况，我们可以利用 Scikit-Learn 库提供的一个特殊的 管道 对象。这个 Scikit-Learn 管道本质上允许您做的是将一组事物链接在一起，如数据工程和模型推理，在一个单一的无缝管道中。只需插入您的原始数据，并在另一边得到一个推论，而无需在中间忙乱数据工程！

我知道我们这样做已经有一段时间了，但是你还记得我们为执行特征工程而创建的所有自定义函数吗？我们将再次使用所有这些文件，为了方便起见，我将它们都放在了 this helpers.py 文件中。这是我们写的其中一个函数的例子:

作者截图

因为我们需要对几乎所有的原始特征执行特征工程，所以我们需要在实际使用 Scikit-Learn Pipeline 对象之前做一些事情。也就是说，我们需要建立一个column transformer，它也是由 Scikit-Learn 库提供的。让我继续向您展示代码，我们将在下面解释它在做什么，因为它非常忙。

作者截图

好吧，尽管这里有很多东西，但实际上很容易理解。基本上，我们正在定义如何正确地转换每一列。在以ohe_engineering开始的行中，我们使用一个 OneHotEncoder 对象对右边的特征primary_genre和secondary_genre进行一次热编码。每个转换器有 3 个部分:转换的名称(您可以随意称呼它)，需要进行什么转换，以及在哪些列上执行转换。

在我们上面提到的ohe_engineering例子中，我们使用了一个“默认的”OneHotEncoder 对象来执行简单的一键编码。但是如果我们想使用我们自己的自定义函数呢？正如你在其他每个转换器中看到的，我们正是在 Scikit-Learn 的 FunctionTransformer 的帮助下做到这一点的。看看以movie_age_engineering开头的那一行。请注意，我们在中间有这段代码:

FunctionTransformer(generate_movie_age, validate = False)

这正是使用我们在上面分享的第一个截图中定义的函数。同样，我们在每个其他转换器中使用我们各自的定制函数，并针对适当的原始列运行它们。非常酷！在继续之前，还要注意最后一个 transformer 是一个默认的columns_to_drop功能，正如您所猜测的，它删除了不需要的特性。此外，确保将remainder参数设置为passthrough，以确保任何其他已经可以从原始状态进入另一端的列。如果这个remainder参数设置不正确，默认行为是删除任何没有明确提到的特性。对我们来说不理想！

好了，我们的 ColumnTransformer 一切就绪，我们现在准备实例化并利用 Scikit-Learn 管道对象。到目前为止，我们已经完成了繁重的工作，所以从这里开始，一切都是下坡路了！请记住，我们在这里有两个单独的模型要训练，所以您在这里基本上会看到双倍的代码。对于这些模型中的每一个，以下是实例化管道对象并执行训练的代码:

作者截图

正如您所看到的管道实例化，实际上集成模型算法作为最终的管道步骤与我们在不使用管道的情况下实例化算法的方式完全相同。此外，请注意，当我们根据训练数据训练管道时，我们会调用与模型算法对象完全相同的fit()函数，而不使用管道。同样，当我们稍后使用这些模型来执行推理时，我们将对经过训练的管道使用相同的predict()，就像我们对任何独立的经过训练的算法一样。不错！

说到推理，我们将在稍后使用这些经过训练的管道，所以让我们用下面的代码将它们导出到序列化的 pickle 文件中:

作者截图

如果你想知道我为什么使用cloudpickle来执行序列化，这是因为我的新 Mac mini 采用了新的苹果硅芯片，与标准的pickle库和 Scikit-Learn 的joblib都不完全匹配。在做了一些研究之后，我发现cloudpickle库工作得很好，它在我的 Windows 笔记本电脑上也工作得很好。

(如果你想看看完整的、有凝聚力的培训脚本是什么样子，请查看这个 train.py 文件。)

这个帖子到此为止！我们还没有完成这个系列。在接下来的几篇文章中，我们将把这些经过训练的模型转换成一个简单的网页，并带有后端 API，允许人们插入电影标题来生成预测的电影评论分数。敬请关注这些帖子，并感谢您的阅读！🎬

创建回归模型以预测数据响应

原文：https://towardsdatascience/creating-regression-models-to-predict-data-responses-120b0e3f6e90

学习从头开始创建和编写回归模型，以预测结果

凯利·西克玛在 Unsplash 上的照片

回归分析可以描述为一种统计技术，用于在给定一个或多个自变量(预测因子或特征)值的情况下，预测/预报因变量(响应)的值。回归被认为是监督机器学习的一种形式；这是一种建立数学模型以确定输入和期望输出之间关系的算法。可以生成的回归模型数不胜数，所以让我们先介绍一下回归数学，然后跳到三个常见的模型。

回归数学

在我们看任何代码之前，我们应该了解一点回归模型背后的数学。如前所述，回归模型可以有多个输入变量或特征，但是对于本文，为了简单起见，我们将使用单个特征。回归分析包括猜测哪种类型的函数最适合您的数据集，是直线函数、n 次多项式函数还是对数函数等。回归模型假设数据集遵循以下形式:

这里， x 和 y 是我们在观察点 i 的特征和响应， e 是误差项。回归模型的目标是估计函数 f ，使其最接近数据集(忽略误差项)。函数 f 是我们对哪种类型的函数最适合数据集的猜测。当我们看到例子时，这将变得更加清楚。为了估计函数，我们需要估计 β 系数。最常见的方法是普通最小二乘法，它可以最小化数据集和函数之间的误差平方和。

最小化平方误差将导致我们的 β 系数的估计，标记为 β -hat。可以使用由普通最小二乘法得出的正规方程来获得估计值:

反过来，我们可以使用这些估计的 β 系数来创建我们的响应变量 y 的估计值。

这里的 y -hat 值是一个观察值 i 的预测值。通过使用我们估计的 β 系数和特征值 x 来生成预测。

为了确定我们的回归模型表现如何，我们可以计算 R 平方系数。此系数评估回归模型预测周围数据集的分散程度。换句话说，这是回归模型与原始数据集的拟合程度。通常，该值越接近 1.0，拟合度越好。

为了正确评估您的回归模型，您应该比较多个模型(函数， f )的 R 平方值，然后决定使用哪个模型来预测未来值。让我们开始编码，以便更好地理解回归模型的概念。

导入库

对于这些示例，我们将需要以下库和包:

• NumPy 用于创建数值数组(为方便调用，定义为 np
• NumPy 中的 random 用于在函数中创建噪声，以模拟真实世界的数据集
• 来自 matplotlib 的 pyplot 用于显示我们的数据和趋势线

# Importing Libraries and Packages
import numpy as np
from numpy import random
import matplotlib.pyplot as plt

创建和显示测试数据集

为了创建数据集，我们可以创建一个作为我们的特征的x数组和一个作为我们的响应的y数组。y数组是一个任意指数函数。使用随机包，我们将在我们的数据中引入噪声来模拟某种真实世界的数据。

# Creating Data
x = np.linspace(0, 4, 500)  # n = 500 observations
random.seed(10)  # For consistent results run-to-run
noise = np.random.normal(0, 1, x.shape)
y = 0.25 * np.exp(x) + noise

# Displaying Data
fig = plt.figure()
plt.scatter(x, y, s=3)
plt.title("Test Data")
plt.xlabel("x")
plt.ylabel("y")
plt.show()

以下是我们应用回归分析的数据集:

测试数据[由作者创建]

让我们花点时间来看看这个情节。我们知道使用了一个指数方程来创建该数据，但是如果我们不知道该信息(如果您正在执行回归分析，您就不会知道)，我们可能会查看该数据并认为二次多项式最适合数据集。然而，作为最佳实践，您应该评估几个不同的模型，看看哪一个性能最好。然后，您将使用最佳模型来创建您的预测。现在让我们分别来看几个回归模型，看看它们的结果如何比较。

线性回归

正如在数学部分提到的，对于回归分析，你猜测一个模型，或函数。对于线性模型，我们的函数采用以下形式， f :

下一步是为估计我们的 β 系数的正规方程(在数学部分)创建矩阵。它们应该如下所示:

使用我们的数据集，我们估计的 β 系数以及线性回归模型将为:

# Linear Regression
X = np.array([np.ones(x.shape), x]).T
X = np.reshape(X, [500, 2])

# Normal Equation: Beta coefficient estimate
b = np.linalg.inv(X.T @ X) @ X.T @ np.array(y)
print(b)

# Predicted y values and R-squared
y_pred = b[0] + b[1] * x
r2 = 1 - sum((y - y_pred) ** 2)/sum((y - np.mean(y)) ** 2)

# Displaying Data
fig = plt.figure()
plt.scatter(x, y, s=3)
plt.plot(x, y_pred, c='red')
plt.title("Linear Model (R$^2$ = " + str(r2) + ")")
plt.xlabel("x")
plt.ylabel("y")
plt.legend(["Data", "Predicted"])
plt.show()

生成的回归模型如下所示:

测试数据[由作者创建]

正如图的标题所示，线性回归模型的 R 平方值约为 0.75。这意味着线性模型符合数据，但如果我们收集更多的数据，我们可能会看到 R 平方的值大幅下降。

多项式回归

多项式回归的一般形式如下:

本文将函数 f 的形式显示为二阶多项式。您可以以类似的方式添加更多的多项式项来满足您的需求；但是，要警惕过度拟合。

创建正规方程的矩阵给出如下:

使用这些矩阵，估计的 β 系数和多项式回归模型将为:

# Polynomial Regression
X = np.array([np.ones(x.shape), x, x ** 2]).T
X = np.reshape(X, [500, 3])

# Normal Equation: Beta coefficient estimate
b = np.linalg.inv(X.T @ X) @ X.T @ np.array(y)
print(b)

# Predicted y values and R-squared
y_pred = b[0] + b[1] * x + b[2] * x ** 2
r2 = 1 - sum((y - y_pred) ** 2)/sum((y - np.mean(y)) ** 2)

# Displaying Data
fig = plt.figure()
plt.scatter(x, y, s=3)
plt.plot(x, y_pred, c='red')
plt.title("2$^{nd}$ Degree Poly. Model (R$^2$ = " + str(r2) + ")")
plt.xlabel("x")
plt.ylabel("y")
plt.legend(["Data", "Predicted"])
plt.show()

由此产生的情节:

测试数据[由作者创建]

请注意，二次多项式在匹配数据集方面做得更好，其 R 平方值约为 0.91。然而，如果我们看数据集的末尾，我们会注意到模型低估了这些数据点。这应该鼓励我们研究其他模型。

指数回归

如果我们决定不采用多项式模型，我们的下一个选择可能是以下形式的指数模型:

同样，创建法线方程的矩阵给出如下:

现在，有了估计的 β 系数，回归模型将是:

请注意，指数项的领先系数与模拟数据领先系数非常匹配。

# Exponential Regression
X = np.array([np.ones(x.shape), np.exp(x)]).T
X = np.reshape(X, [500, 2])

# Normal Equation: Beta coefficient estimate
b = np.linalg.inv(X.T @ X) @ X.T @ np.array(y)
print(b)

# Predicted y values and R-squared
y_pred = b[0] + b[1]*np.exp(x)
r2 = 1 - sum((y - y_pred) ** 2)/sum((y - np.mean(y)) ** 2)

# Displaying Data
fig = plt.figure()
plt.scatter(x, y, s=3)
plt.plot(x, y_pred, c='red')
plt.title("Exponential Model (R$^2$ = " + str(r2) + ")")
plt.xlabel("x")
plt.ylabel("y")
plt.legend(["Data", "Predicted"])
plt.show()

最终的数字:

测试数据[由作者创建]

请注意，与二次多项式模型相比，R 平方值有所提高。更有希望的是，数据集的末端比二次多项式更适合此模型。给定这些信息，我们应该使用指数模型来预测未来的 y 值，以获得尽可能精确的结果。

文章到此结束。希望这能让您了解如何在 Python 中执行回归分析，以及如何在现实场景中正确使用它。如果你学到了什么，给我一个关注，看看我关于空间、数学和 Python 的其他文章！谢谢大家！

为人工智能平台创建路线图

原文：https://towardsdatascience/creating-roadmap-for-internal-ai-platform-product-in-an-enterprise-ce60b77df25f

AI 平台团队的产品路线图应该考虑哪些？

灰色滚柏油路下多云的天空。照片由皮查拜从派克斯拍摄。

组织的技术战略目标主要通过其技术平台来实现。这些产品的路线图对于获得组织内关键利益相关者的认可至关重要。因此，路线图对于将技术战略引向正确的方向至关重要。

**路线图在高层次上描述了团队实现其产品愿景的方式。**它可以用许多不同的方式来表现，比如看板、幻灯片、excel 文件等。，显示时间线上的详细信息。

人工智能平台产品路线图

一个组织的人工智能愿景需要具有不同能力的平台产品。本组织需要以下平台

• 应用研究
• 代码版本控制、持续集成和持续部署
• 数据摄取、功能管理
• 模型版本控制，模型管理
• 大规模在线/批量培训
• 大规模在线/批量服务
• …

通常情况下，组织无法为每一项功能提供单独的产品团队。相反，一个团队选择在一个人工智能平台产品下解决其中的几个功能，该产品结合了供应商管理/提供的以及自行开发的解决方案。这种产品的路线图定义了如何在可预见的时间内实现产品的最终愿景，并指示团队如何组合漏斗以实现愿景。

打算使用这样一个平台的人工智能团队是客户。参与批准预算和承担 AI 平台开发和维护成本的组织内的团体是其利益相关者。其他团队，如数据平台团队，其路线图受其路线图影响，也是其利益相关方。

路线图的质量

为人工智能平台产品提出一个好的路线图是很棘手的，特别是当多个团队使用的产品可能在功能上没有明确的界限，并且将购买和开发的解决方案结合在一起时。我们举几个例子来说明下面的复杂性。

问题 1:产品路线图与组织 AI 愿景不一致

通常，提供人工智能平台解决方案的技术供应商，例如可扩展的笔记本电脑环境，正在向团队追加销售/推广他们的功能，以确保他们的未来。同时，团队成员也对新的平台范例和展示/提高他们的技能感到兴奋。当平台通过这些漏斗的目标与组织的人工智能愿景相匹配时，这就是天作之合。然而，团队的目标与组织的目标不一致，很难理解为什么事情会在团队的路线图中。这样的团队获得了太多不必要的想法，面临着对特性和里程碑的大量怀疑，在最坏的情况下，不得不将产品束之高阁。

问题 2:产品输出未能交付客户成果

对于一个人工智能平台产品团队来说，当它有太多的产出，很少或根本没有为客户的人工智能雄心提供什么时，这是一个不好的迹象。这种情况经常发生，因为平台团队不清楚他们的交付如何影响他们客户的底线。在这些情况下，团队观察到实现里程碑对他们的人工智能驱动的业务和客户满意度的 KPI 没有什么影响。

问题 3:团队几乎没有时间去学习什么是有效的

通常情况下，客户和利益相关者会强迫或推动平台产品团队完成一个紧凑的特性交付时间表。然而，找出高影响力特性的道路通常并不简单。相反，团队需要尝试一些东西，并从中学习，直到他们发现什么是真正有效的。致力于详细的特性交付，迫使团队专注于单个解决方案轨道。这太冒险了，因为这样的策略对团队来说是一个障碍，而尝试可能有助于团队找出适当的细节。

问题 4:客户和利益相关者对路线图不感兴趣

有时，人工智能团队会推迟使用平台的新功能。此外，尽管源于组织人工智能愿景的人工智能团队稳步增长，但平台中客户的采用率很低。人工智能团队也避免挑战路线图或在适当的时候向平台团队提供反馈，即使这样做，他们也太迟了。另一方面，涉众可能会羞于提升平台团队的能力，即使他们这么做了，他们也是针对旧的特性。这些症状表明，客户和利益相关者都对平台产品的路线图不感兴趣。

确保良好的路线图

为了提供良好的路线图，平台产品团队应该考虑以下几点:

1. 团队应该将平台的路线图与组织的关键人工智能愿景联系起来。这要求团队中的每个人都应该意识到组织的宏观愿景，并尽可能经常地讨论它。
2. 平台团队应该专注于为人工智能团队取得成果，而不仅仅是交付功能。一旦一个平台团队掌握了它的全局，它就不应该进入解决方案的细节，即使这对团队来说是显而易见的。相反，团队应该关注人工智能团队的高层次价值交付。
3. 平台团队不应该为了说服涉众和客户而致力于不切实际的、雄心勃勃的和不必要的特性。它应该强调结果的重要性并坚持下去。利益相关者和客户不应该规定实现结果的方法。
4. 平台团队不应该选择太容易的议程，并且在组织的人工智能愿景中产生很小的影响。同样，它应该避免采取没有长期持续时间的速赢措施。

人工智能平台路线图框架

像大多数软件产品一样，人工智能平台产品路线图的良好框架应该包含以下组件:

• 愿景:一旦平台产品完全实现，人工智能团队将如何从中受益？
• **时间框架:**按时间顺序排列的交付里程碑。每个时间段的持续时间可能不同，时间越早，持续时间越短。
• **结果:**如果组织中的人工智能团队在给定的时间框架内使用平台产品中已交付的功能，他们会有什么可测量的不同？
• **主题:**里程碑期间要实现的成果所涉及的共同主题。
• **风险:**哪些问题会导致无法在规定时间内交付成果？
• **免责声明:**可能改变路线图的警告。

在下图中，我们展示了一个 AI 平台产品的路线图示例。框架的内容并不那么重要。更重要的是把内容作为理解作品设计的素材。

人工智能平台产品路线图示例。图片由作者提供。

评论

路线图有多种形状和形式。你的产品路线图遵循什么样的设计？你同意或同意以上陈述吗？让我知道你的想法。

用微软行星计算机创建哨兵 2 号(真正的)无云镶嵌图

原文：https://towardsdatascience/creating-sentinel-2-truly-cloudless-mosaics-with-microsoft-planetary-computer-7392a2c0d96c

使用 GEE 的 S2 无云图层，在微软行星计算机上有效地遮蔽云(和云的阴影)

美国宇航局在 Unsplash 拍摄的照片

介绍

在地理空间领域，微软行星计算公司是久负盛名的谷歌地球引擎的有力竞争者。无需下载每张图像即可访问数 Pb 卫星信息的可能性以及由 Dask clusters 提供的计算能力是开发区域/全球应用和研究的一大变革。行星计算机仍处于预览阶段(需要访问请求)，但使用其计算中心与开源工具如 XARRAY 和 STAC 的可能性是优于谷歌专有 API 的一个优势(IMHO)。

另一个优点是，我们不必像在 GEE 中那样受大小约束的困扰，大小约束会阻止我们轻松地访问 Numpy 数组格式的值。好吧，我知道 GEE 的概念是不同的，这意味着迫使人们在服务器端运行计算(使用惰性数组的概念)，但有时直接访问数据更容易。在 Planetary Computer 中，甚至可以使用“自制”下载器直接从微软的目录中下载全部资产(更多内容将在以后的文章中介绍)。

不利的一面是，数据集目录不像在谷歌的竞争对手那里找到的那样广泛。在处理光学(哨兵 2 号)影像时，有一件事真的很重要:云层遮罩。

在微软的官方无云镶嵌教程(这里是)中，他们陈述了以下内容:*“在云是短暂的假设下，合成图像不应该包含(许多)云，因为它们不应该是许多图像在该点的中间像素值。”。那是…呃…不完全是我们所期待的。*在 L2A S2 影像中总会有场景分类图层——SCL 图层，但之前用过 Sen2Cor 的人都知道这不是精度最好的。 Baetens et al. (2019) 在这个课题上带来了全面的比较。

这又把我们带回了我们的故事:如果我们运行在微软的行星计算机上，我们如何才能有效地从哨兵 2 的图像中去除云。正如在行星计算机范例 GitHub 知识库中提到的(这里是这里是)这是他们正在做的事情，但是在那之前…

S2 无云算法

S2 无云包是 Sinergise 开发的机器学习云检测算法，在 GitHub 上有售(此处)。欧空局正在使用这种算法，它最近作为一张云概率图被列入了 GEE 的目录。因此，我们可以直接从 GEE 访问概率图，而不是为每个场景运行训练好的算法。这些问题是:

• 我们想使用行星计算机，但我们不在地球环境中！而且，
• 云影呢？

为了解决前面提到的两点，我结合使用了 geeS2Downloader 包(关于这个故事的更多信息这里)和 GEE 的页面上提供的教程 Sentinel-2 Cloud Masking with S2 Cloud less(这里)。该教程显示了如何根据太阳方位角投影云阴影，以找到实际的地面阴影，geeS2Downloader 包使从 GEE 下载资源更容易，克服了其大小限制。如果看起来太复杂，让我们来看看完整的解决方案…

解决方案

让我们首先从微软行星中心的一个空笔记本开始，安装依赖项。在这种情况下，我们将需要两个未预装在 PC(行星计算机)中的软件包和一个用于更快可视化的可选软件包:

• earthengine-api:从谷歌地球引擎访问资产。
• geeS2Downloader:从 GEE 下载资产。
• Geemap:来自吴秋生教授的伟大软件包，它的最新版本也可以在微软的个人电脑上运行。

初始化后，我们将通过 TILE_ID 和日期搜索特定的图块。为此，我们将使用一个名为search_tiles的函数。

<Item id=S2B_MSIL2A_20181214T133219_R081_T22KFV_20201008T100849>

然后，我们需要在 GEE 中得到相应的云概率图。另外，为了投射阴影，我们也需要完整的 S2 图像。为此，我们将创建一个名为get_gee_img 的新函数，它将负责在给定 STAC 商品的情况下定位 GEE 目录中的图像。

现在是时候看看我们用 Geemap 下载的图片了:

代码输出。图片作者。

正如我们所看到的，云被正确识别，但我们仍然需要摆脱阴影。为此，我们将创建一个受 GEE 教程启发的函数create_cloud_mask。它将扩大最终蒙版与 50 米的缓冲区，并重新调整到 20 米的分辨率。

{'type': 'Image',
 'bands': [{'id': 'cloudmask',
   'data_type': {'type': 'PixelType', 'precision': 'int', 'min': 0, 'max': 1},
   'dimensions': [5490, 5490],
   'crs': 'EPSG:32722',
   'crs_transform': [20, 0, 600000, 0, -20, 7500040]}]}

正如我们可以从代码的输出中看到的，掩码似乎是用 20m 正确创建的，正如预期的那样。要在 geemap 上显示它，我们只需使用以下命令添加该层:

Map.addLayer(mask.selfMask(), {'min': 0, 'max': 1, 'palette': ['orange']}, 'Final Mask', True, 0.5)

这是最终的蒙版输出。请注意，我们不仅有云彩，还有阴影(橙色)。

代码输出。图片作者。

将掩模下载到 PC

现在我们已经在 GEE 中处理了遮罩，是时候下载它了，看看它是否正确地符合我们的原始图像。为此，我们将使用 geeS2Downloader 包。这个软件包将根据 GEE 的限制自动分割蒙版并重建最终的矩阵。

代码输出。图片作者。

显示结果

现在，让我们显示最终结果，以与 PC 中的图像进行比较。代替 geemap ，我们将使用普通的 Matplotlib 来完成这个任务。

代码输出。图片作者。

放大输出。图片作者。

结论

正如我们在这个故事中看到的，将谷歌地球引擎和微软行星计算机的资产结合起来，从两个平台中提取精华是可能的。虽然微软的个人电脑不包括可靠的云端层，但这是我为了继续我的项目而创建的一个变通办法。在未来，我预计这种变通办法将不再是必要的。直到那里，希望它能帮助我们中的一些人！

谢谢，下一个故事再见！

保持联系

如果你喜欢这篇文章，并想继续无限制地阅读/学习这些和其他故事，考虑成为 中等会员 *。你也可以在 https://cordmaur.carrd.co/*查看我的作品集。

https://cordmaur.medium/membership

参考

Baetens，l .，Desjardins，c .，Hagolle，o .，2019。使用有监督的主动学习程序生成的参考云掩膜对从 MAJA、Sen2Cor 和 FMask 处理器获得的哥白尼 Sentinel-2 云掩膜的验证。遥感 11433。https://doi/10.3390/rs11040433

创建合成数据

原文：https://towardsdatascience/creating-synthetic-data-3774391c851d

使用模型来提供数据中心性

我以前写过关于合成数据的博客。实际上是两次。我第一次使用荷兰癌症登记处的合成数据进行生存分析。第二次我使用相同的数据集来应用规格曲线分析。

合成数据将在未来几年发挥关键作用，因为模型变得越来越复杂，但受到可用数据的限制。一段时间以来，以模型为中心的运动正慢慢被更以数据为中心的观点所取代，随着时间的推移，模型将再次占据中心舞台。

合成数据实际上是两个世界的结合，还处于起步阶段。它的数据，而是我们习惯于看到和看待数据的方式。对于大多数研究人员来说，数据来自实地观察或传感器。合成数据是模型数据。本质上，它的数据伪装成其他数据，其中敏感组件被删除，但关键组件保留。如果成功，建模者甚至不应该知道数据是合成的。或许一个很好的类比是描绘你站在蒙娜丽莎面前的卢浮宫，但它不是真正的蒙娜丽莎。真正的蒙娜丽莎藏在地下室，你看到的是一个戒备森严的复制品。

合成数据是关于尊重原始的本质，而不是原始的。从这一点出发，合成数据甚至可以更多。可以成为现实生活中很少会出现的模型的训练素材，但是可以作为模拟素材分享给大家。因为数据失去了敏感性，但没有失去价值。

在这个模块中，我要做的是获取 R 中的数据集，即钻石数据集，并从中创建几个模型来建立合成数据。当然，钻石数据集已经开放，因此它没有敏感信息，但是这个练习的本质是在不重新创建数据集的情况下重新创建数据集的本质。

我使用 diamonds 数据集的原因是它有足够的行和列，但没有直接和明确的关系。

让我们开始吧。

rm(list = ls())
options("scipen"=100, "digits"=10)
#### LIBRARIES ####
library(DataExplorer)
library(ggplot2)
library(dplyr)
library(tidyr)
library(ggExtra)
library(skimr)
library(car)
library(GGally)

library(doParallel)
cl <- makePSOCKcluster(5)
registerDoParallel(cl)

加载完库之后，我想立即查看数据本身。如果我要对数据建模，以识别其本质(通过交叉相关)，并从本质上复制它，那么我必须牢牢掌握数据集。为此，我将首先查看协方差和相关矩阵。

skimr::skim(diamonds)
str(diamonds)
diamonds$price<-as.numeric(diamonds$price)
diamonds%>%
  dplyr::select(where(is.numeric))%>%
  cor()
diamonds%>%
  dplyr::select(where(is.numeric))%>%
  corrgram::corrgram()
diamonds%>%
  dplyr::select(where(is.numeric))%>%
  cor()%>%
  ggcorrplot::ggcorrplot(., hc.order = TRUE, type = "lower",
             lab = TRUE)

diamonds 数据集有许多行，但只有几列。这应该会使建立模型变得更容易，但是关系并不明确。图片作者。

我将首先从数值开始，因为这是最容易的开始。这并不意味着我们不能为序数或名词性数据重建合成数据。然而，非连续数据的相关矩阵需要一些额外的变换，我们现在不考虑这些变换。

所有数值变量的相关矩阵。图片作者。

和三种类型的相关矩阵。在左侧，您可以看到一个附有树状图的热图。树状图是一种涉及相关性的聚类技术。因此，将它们结合起来是很简单的。正如你所看到的，有一些严重的相关性和一些部分几乎是分离的。图片作者。

我们甚至可以创建更奇特的图，但是 R 在绘制所有 53k 行时会遇到一点问题。因此，我将首先选择一个 1000 行的随机样本，然后绘制每个数值的分布及其潜在的相关性。

diamonds%>%
  dplyr::select(where(is.numeric))%>%
  sample_n(., 1000)%>%
  scatterplotMatrix()

这种分布应该是多峰值的，这很有趣。此外，变量之间的关系有时是极端的，有时是不存在的(你也可以称之为极端)。图片作者。

而且，有一个很好的方法可以使用 GGally 库来绘制每个特定类别(这里是 cut )的相关性和分布。

diamonds%>%
  dplyr::select(x, y, z, carat, depth, table, price, cut)%>%
  sample_n(., 1000)%>%
  ggpairs(.,
  mapping = ggplot2::aes(color = cut),
  upper = list(continuous = wrap("density", alpha = 0.5), 
               combo = "box_no_facet"),
  lower = list(continuous = wrap("points", alpha = 0.3), 
               combo = wrap("dot_no_facet", alpha = 0.4)),
  title = "Diamonds")

图片作者。

如果我们要基于只包含数值的原始数据重新创建数据集，我们需要立即考虑需要注意的两个主要问题:

1. 每个变量的汇总统计和分布需要尽可能接近原始数据，但不需要完全相同。这意味着需要维护每个变量的结构。如果不是这样，数据集的描述部分就不能使用。造型部分仍可保留。
2. 数据点之间的协方差/相关性需要相同。这意味着需要维护变量之间的潜在关系。这对造型部分至关重要。

从数据中提取协方差/相关矩阵并不困难，但正如我们所说，标准公式只适用于数值。尽管如此，我们可以用数字值来开始我们的创作，让这个想法继续下去。

这里你可以看到协方差矩阵。

diamond_cov<-diamonds%>%dplyr::select_if(., is.numeric)%>%cov()

                 carat           depth           table             price
carat    0.22468665982   0.01916652822    0.1923645201     1742.76536427
depth    0.01916652822   2.05240384318   -0.9468399376      -60.85371214
table    0.19236452006  -0.94683993764    4.9929480753     1133.31806407
price 1742.76536426512 -60.85371213642 1133.3180640679 15915629.42430145
x        0.51848413024  -0.04064129579    0.4896429037     3958.02149078
y        0.51524781641  -0.04800856925    0.4689722778     3943.27081043
z        0.31891683911   0.09596797038    0.2379960448     2424.71261297
                     x                y                z
carat    0.51848413024    0.51524781641    0.31891683911
depth   -0.04064129579   -0.04800856925    0.09596797038
table    0.48964290366    0.46897227781    0.23799604481
price 3958.02149078326 3943.27081043196 2424.71261297033
x        1.25834717304    1.24878933406    0.76848748285
y        1.24878933406    1.30447161384    0.76731957995
z        0.76848748285    0.76731957995    0.49801086259

从这个过程中，我们获得了协方差，但我们不能创建数据。为此，我们还需要可以很容易获得的汇总统计数据。

diamonds%>%dplyr::select_if(., is.numeric)%>%summary()

     carat               depth             table              price         
 Min.   :0.2000000   Min.   :43.0000   Min.   :43.00000   Min.   :  326.00  
 1st Qu.:0.4000000   1st Qu.:61.0000   1st Qu.:56.00000   1st Qu.:  950.00  
 Median :0.7000000   Median :61.8000   Median :57.00000   Median : 2401.00  
 Mean   :0.7979397   Mean   :61.7494   Mean   :57.45718   Mean   : 3932.80  
 3rd Qu.:1.0400000   3rd Qu.:62.5000   3rd Qu.:59.00000   3rd Qu.: 5324.25  
 Max.   :5.0100000   Max.   :79.0000   Max.   :95.00000   Max.   :18823.00  
       x                   y                   z            
 Min.   : 0.000000   Min.   : 0.000000   Min.   : 0.000000  
 1st Qu.: 4.710000   1st Qu.: 4.720000   1st Qu.: 2.910000  
 Median : 5.700000   Median : 5.710000   Median : 3.530000  
 Mean   : 5.731157   Mean   : 5.734526   Mean   : 3.538734  
 3rd Qu.: 6.540000   3rd Qu.: 6.540000   3rd Qu.: 4.040000  
 Max.   :10.740000   Max.   :58.900000   Max.   :31.800000 

从这些组合中，我们应该能够使用各种程序重建数据。也许最直接的方法是使用多元正态分布，它存在于质量包中。我只需要每个变量的平均值和相关矩阵。多元正态分布将完成剩下的工作。为了获得相等的比较，我将创建与原始数据集一样多的观察值。

sigma<-diamonds%>%dplyr::select_if(., is.numeric)%>%cor()%>%as.matrix()
mean<-diamonds%>%dplyr::select_if(., is.numeric)%>%as.matrix()%>%colMeans()
df<-as.data.frame(MASS::mvrnorm(n=dim(diamonds)[1], mu=mean, Sigma=sigma))
> dim(df)
[1] 53940     7
> head(df)
           carat       depth       table       price           x           y
1 -1.34032717822 61.49447797 57.19091527 3931.162855 3.545669821 3.931061364
2 -0.47751648630 61.16241371 57.86509627 3931.306295 4.906696688 5.057863929
3  2.24358594522 63.09062530 56.73718104 3932.957682 7.386140807 7.405936831
4 -0.03108967881 60.99439588 57.58369767 3931.677322 5.194041948 5.431802322
5  1.16179859890 62.39235813 57.96524508 3933.322044 6.213609464 6.592872841
6 -0.16757252197 60.84783867 56.68337288 3932.501268 4.987489939 5.118558015
            z
1 1.138732182
2 2.596622246
3 5.674154202
4 3.089565271
5 4.387662667
6 2.577509666

创建完成后，接下来的任务是应用两种方法检查数据的有效性和可用性:

1. 检查单变量特征。
2. 检查多元特征。

有了 ggpairs 函数，我可以两者兼而有之，并初步了解该过程及其生成的数据的适用性。

diamonds%>%dplyr::select_if(., is.numeric)%>%ggpairs()
ggpairs(df)

左边是原始数据集，右边是从原始数据中提取平均值和相关矩阵的模拟数据。当然，合成数据的分布遵循正态分布。数据之间的相关性得到了维护，但是汇总统计数据肯定没有得到维护(除了平均值)。图片作者。

现在，我们已经说过，我们的目的是构建合成数据，这意味着构建本质上相同的数据，但在前景上不一定相同。我们已经实现了这个目标。我将使用 caret 包在两个数据集上构建一个快速模型，看看合成数据是否会给我与原始数据集完全相同的结果。

diamonds_num<-diamonds%>%dplyr::select_if(., is.numeric)
trainIndex <- caret::createDataPartition(diamonds_num$carat, 
                                         p = .6, 
                                         list = FALSE, 
                                         times = 1)
> wideTrain <- diamonds_num[ trainIndex,];dim(wideTrain)
[1] 32366     7
> wideTest  <- diamonds_num[-trainIndex,];dim(wideTest)
[1] 21574     7

fitControl <- caret::trainControl(
  method = "repeatedcv",
  number = 20,
  repeats = 20)

lmFit1 <- caret::train(carat ~ ., 
                        data = wideTrain, 
                        method = "lm",
                        trControl = fitControl,
                        verbose = FALSE)
> summary(lmFit1)

Call:
lm(formula = .outcome ~ ., data = dat, verbose = FALSE)

Residuals:
        Min          1Q      Median          3Q         Max 
-0.54124488 -0.03836279 -0.00665401  0.03530248  2.71983984 

Coefficients:
                    Estimate       Std. Error   t value               Pr(>|t|)    
(Intercept) -2.5273956761070  0.0294158001479 -85.91966 < 0.000000000000000222 ***
depth        0.0188998891767  0.0003766002370  50.18555 < 0.000000000000000222 ***
table        0.0046678986326  0.0002266381898  20.59626 < 0.000000000000000222 ***
price        0.0000330063643  0.0000002519416 131.00798 < 0.000000000000000222 ***
x            0.2956915579921  0.0032315175272  91.50238 < 0.000000000000000222 ***
y            0.0130670340617  0.0028968529147   4.51077            0.000006482 ***
z           -0.0026197077889  0.0025495576103  -1.02751                0.30419    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.0844222 on 32359 degrees of freedom
Multiple R-squared:  0.9684143, Adjusted R-squared:  0.9684085 
F-statistic:   165354 on 6 and 32359 DF,  p-value: < 0.00000000000000022204

以上是基于原始数据的简单线性回归的总结。下面是基于合成数据的简单线性回归的总结。是的，程序并不完全相同，根据训练和测试数据的选择，以及反复的交叉验证，可以预期会有一些变化，但是程序的本质是看数据的本质是否被保留。

trainIndex <- caret::createDataPartition(df$carat, 
                                         p = .6, 
                                         list = FALSE, 
                                         times = 1)
> wideTrain <- df[ trainIndex,];dim(wideTrain)
[1] 32364     7
> wideTest  <- df[-trainIndex,];dim(wideTest)
[1] 21576     7

fitControl <- caret::trainControl(
  method = "repeatedcv",
  number = 20,
  repeats = 20)
lmFit2 <- caret::train(carat ~ ., 
                        data = wideTrain, 
                        method = "lm",
                        trControl = fitControl,
                        verbose = FALSE)
summary(lmFit2)

Call:
lm(formula = .outcome ~ ., data = dat, verbose = FALSE)

Residuals:
        Min          1Q      Median          3Q         Max 
-0.68929576 -0.11680934  0.00078743  0.11741046  0.74064948 

Coefficients:
                   Estimate      Std. Error    t value               Pr(>|t|)    
(Intercept) -1079.584562077     8.210331042 -131.49099 < 0.000000000000000222 ***
depth           0.052707918     0.001155147   45.62874 < 0.000000000000000222 ***
table           0.019422785     0.001035408   18.75859 < 0.000000000000000222 ***
price           0.272537301     0.002088731  130.47987 < 0.000000000000000222 ***
x               0.708827904     0.006118409  115.85167 < 0.000000000000000222 ***
y               0.015186425     0.004344156    3.49583             0.00047322 ***
z               0.007592425     0.004694642    1.61725             0.10583335    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1736463 on 32357 degrees of freedom
Multiple R-squared:  0.9697885, Adjusted R-squared:  0.9697829 
F-statistic: 173109.8 on 6 and 32357 DF,  p-value: < 0.00000000000000022204

一目了然的是截距的大小发生了变化，这是因为合成数据具有相同的均值，但分布不同。了解数据的另一种方法是查看试图预测钻石克拉的模型中每个预测值的方差重要性。

> lm1Imp;lm2Imp
lm variable importance

         Overall
price 131.007979
x      91.502384
depth  50.185548
table  20.596258
y       4.510769
z       1.027515

         Overall
price 130.479871
x     115.851671
depth  45.628736
table  18.758588
y       3.495829
z       1.617253

似乎这些模型确实具有相同的可变重要性，并且它们之间的距离似乎被保留了下来。当然，它并不完美，但也不必如此。

了解合成数据可用性的第二种方法是查看累积分布图。累积密度图是显示变量分布形式的另一种方式，但在这种方式中，如果两个数据集显示相同的本质，就会立即变得清楚。

重要的是要认识到，没有单一的方法来确定合成数据是否保持了与原始数据相同的本质。

df$Model<-"Synthetic"
diamonds_num$Model<-"Original"
combined<-rbind(df, diamonds_num)
ggplot(combined, aes(carat, colour = Model)) +
  stat_ecdf()+
  theme_bw()

在这里，您可以看到我尝试使用简单的线性回归建模的结果。正如你所看到的，合成数据非常清晰，相比之下，原始数据是凹凸不平的。因此，包含相关矩阵的简单多元正态分布是不够的。图片作者。

下一步是为数据集中的每个数值变量绘制图表。

combined_long<-combined%>%tidyr::pivot_longer(!Model, 
                                              names_to = "Variable",
                                              values_to = "Value",
                                              values_drop_na = TRUE)
ggplot(combined_long, aes(Value, colour = Model)) +
  stat_ecdf()+
  theme_bw()+
  facet_wrap(~Variable, scales="free")

如你所见，累积分布函数(CDF)在原始数据和合成数据之间没有太大差异。合成数据更加原始，对于克拉，我们已经看到了某些偏差，这可能会导致原始数据和合成数据之间的分析差异。然而，这与绝对错误的价格相比根本不算什么。图片作者。

让我们深入挖掘价格并比较使用密度。我将首先绘制严重偏离的价格，然后绘制几乎相同密度结构的深度。

g1<-ggplot()+
  geom_density(data=diamonds_num, aes(x=price, fill="Original"), alpha=0.5)+
  theme_bw()+
  labs(fill="")
g2<-ggplot()+
  geom_density(data=df, aes(x=price, fill="Synthetic"), alpha=0.5)+
  theme_bw()+
  labs(fill="")
gridExtra::grid.arrange(g1,g2,nrow=2)

合成数据集中价格的 cdf 看起来如此奇怪的原因是因为它必须以与原始数据相同的比例绘制。而且原版的最大值超过 15000，不像合成版停在 3950 左右。图片作者。

这种差异是构成多元正态的函数的直接结果，这意味着每个变量都有一个平均值，并且与其他函数相关。在这个特殊的例子中，平均值并不能说明全部情况。

ggplot()+
  geom_density(data=diamonds_num, aes(x=depth, fill="Original"), alpha=0.5)+
  geom_density(data=df, aes(x=depth, fill="Synthetic"), alpha=0.5)+
  theme_bw()+
  labs(fill="")

原始数据和合成数据并不相同，但确实表现出相似的特征。然而，如果从原始数据中提取汇总统计数据，并将其与合成数据进行比较，这种方法充其量也是有限的。这是因为描述性统计需要完全相同的值，这意味着分布应该完全相同。这就是综合数据集的描述部分和建模部分之间的区别。图片作者。

如果我采用质量包的多元正态分布的经验形式会怎样——这意味着样本大小会发挥更大的作用。

df_emp<-as.data.frame(MASS::mvrnorm(n=dim(diamonds)[1], 
                                    mu=mean, Sigma=sigma, empirical = TRUE))

g1<-ggplot()+
  geom_density(data=diamonds_num, aes(x=price, fill="Original"), alpha=0.5)+
  theme_bw()+
  labs(fill="")
g2<-ggplot()+
  geom_density(data=df, aes(x=price, fill="Synthetic"), alpha=0.5)+
  theme_bw()+
  labs(fill="")
g3<-ggplot()+
  geom_density(data=df_emp, aes(x=price, fill="Synthetic Emp"), alpha=0.5)+
  theme_bw()+
  labs(fill="")
gridExtra::grid.arrange(g1,g2,g3,nrow=3)

对于 53k 行数据，两个多变量正态程序之间没有实际差异，这是可以预期的。样本大小在这里不是问题。图片作者。

好的，多元正态分布确实为创建合成数据集提供了一个良好的开端，但只是从建模的角度来看。不是从描述的角度，人们需要自己决定是否有必要。数据合成的部分本质是确保本质得到保持，对于建模者来说，这意味着建模时能够得到相同的结果。

现在，创建合成数据的另一种方法(这意味着模拟相关变量)是深入到连接函数的世界中，我们现在将在某种程度上使用高斯函数。连接函数是理解和建立多元分布的连接概率的一个很好的方法。copula 这个词的意思是“链接”,这正是他们所做的。

根据维基百科，一个 copula 是:一个多元 累积分布函数 其中 边际概率 每个变量的分布是https://en.wikipedia/wiki/Uniform_distribution_(continuous)**均匀分布在区间[0，1]上。如果我们将这些步骤分解开来，看起来会是这样的(这是我从博客上引用的):

1. 多元正态分布的样本相关标准化(N[0，1])分布。
2. 用正态 CDF 将它们转换成相关的均匀(0，1)分布(概率积分转换)。
3. 用概率分布的逆 CDF 将它们转换成你想要的任何相关概率分布(逆变换采样)。

下面是一个函数，显示了一个构建的多元正态分布函数，它等于 MASS 包的 mvnorm 函数。因此，我将从多元正态分布中获得(就像以前一样)数据，但这次是标准化的。

 *mvrnorm <- function(n = 1, mu = 0, Sigma) {
  nvars <- nrow(Sigma)
  # nvars x n matrix of Normal(0, 1)
  nmls <- matrix(rnorm(n * nvars), nrow = nvars)
  # scale and correlate Normal(0, 1), "nmls", to Normal(0, Sigma) by matrix mult
  # with lower triangular of cholesky decomp of covariance matrix
  scaled_correlated_nmls <- t(chol(Sigma)) %*% nmls
  # shift to center around mus to get goal: Normal(mu, Sigma)
  samples <- mu + scaled_correlated_nmls
  # transpose so each variable is a column, not
  # a row, to match what MASS::mvrnorm() returns
  t(samples)
}
df_new <- mvrnorm(dim(diamonds)[1], Sigma = sigma)
mean2<-rep(0,dim(sigma)[2])
names(mean2)<-colnames(sigma)
df_new2 <- MASS::mvrnorm(dim(diamonds)[1], 
                         mu=mean2, 
                         Sigma = sigma)
> cor(df_new)
              carat          depth         table          price              x              y
carat 1.00000000000  0.02548515939  0.1756632819  0.92081606682  0.97483969924  0.95123398113
depth 0.02548515939  1.00000000000 -0.2966914653 -0.01267406029 -0.02968316709 -0.03193223836
table 0.17566328191 -0.29669146532  1.0000000000  0.12037212430  0.18934326066  0.17732479479
price 0.92081606682 -0.01267406029  0.1203721243  1.00000000000  0.88323677536  0.86373468972
x     0.97483969924 -0.02968316709  0.1893432607  0.88323677536  1.00000000000  0.97460644946
y     0.95123398113 -0.03193223836  0.1773247948  0.86373468972  0.97460644946  1.00000000000
z     0.95342001221  0.09193958958  0.1450393656  0.86003163701  0.97049061902  0.95189367284
                  z
carat 0.95342001221
depth 0.09193958958
table 0.14503936558
price 0.86003163701
x     0.97049061902
y     0.95189367284
z     1.00000000000
> cor(df_new2)
              carat          depth         table          price              x              y
carat 1.00000000000  0.02401766053  0.1879023338  0.92211539384  0.97544578111  0.95144071623
depth 0.02401766053  1.00000000000 -0.3013205641 -0.01515267326 -0.02925527573 -0.03155516412
table 0.18790233377 -0.30132056412  1.0000000000  0.13488569935  0.20173295788  0.18904146957
price 0.92211539384 -0.01515267326  0.1348856993  1.00000000000  0.88524849733  0.86534901465
x     0.97544578111 -0.02925527573  0.2017329579  0.88524849733  1.00000000000  0.97451646301
y     0.95144071623 -0.03155516412  0.1890414696  0.86534901465  0.97451646301  1.00000000000
z     0.95428998681  0.09122255285  0.1570967255  0.86180784719  0.97117225153  0.95206791685
                  z
carat 0.95428998681
depth 0.09122255285
table 0.15709672554
price 0.86180784719
x     0.97117225153
y     0.95206791685
z     1.00000000000*

因此，我们首先获得的值是来自多元正态分布的标准化值。这意味着变量都在相同的尺度上，并携带原始的互相关矩阵。我们可以很容易地检查两者(标准化规模和相关性)。

*pairs(df_new)
hist(df_new[,1])*

**

左图:相关矩阵。右图:配送。都来自标准化的多元正态分布。图片作者。

下一步是转换到均匀分布，同时保持基本的互相关矩阵。

*U <- pnorm(df_new, mean = 0, sd = 1)
hist(U[,1])
cor(U)
              carat          depth         table         price              x              y             z
carat 1.00000000000  0.02617867093  0.1682925760  0.9140135553  0.97247036750  0.94681926838 0.94931607318
depth 0.02617867093  1.00000000000 -0.2837151691 -0.0102210807 -0.02657673258 -0.02879604409 0.08933946454
table 0.16829257601 -0.28371516908  1.0000000000  0.1160484259  0.18057212310  0.16815896730 0.13787984049
price 0.91401355528 -0.01022108070  0.1160484259  1.0000000000  0.87409215149  0.85353166518 0.84964860456
x     0.97247036750 -0.02657673258  0.1805721231  0.8740921515  1.00000000000  0.97228684415 0.96785128673
y     0.94681926838 -0.02879604409  0.1681589673  0.8535316652  0.97228684415  1.00000000000 0.94769217826
z     0.94931607318  0.08933946454  0.1378798405  0.8496486046  0.96785128673  0.94769217826 1.00000000000*

均匀分布矩阵中的第一个变量是克拉，现在也是均匀分布的。保持与所有其他数据的互相关。图片作者。

我们可以对不同的变量做同样的处理，比如价格*。下面，你会看到我构建了新的价格和克拉的变量，但这次我是从泊松分布中对它们进行采样。这是三步中的最后一步，我可以通过逆变换采样，从包含均匀分布数据的矩阵 U 中创建任何我想要的分布。这是一个相当酷的技术！*

*price <- qpois(U[, 4], 5)
par(mfrow = c(2, 1))
hist(price)
hist(diamonds_num$price)

carat <- qpois(U[, 1], 30)
par(mfrow = c(2, 1))
hist(carat)
hist(diamonds_num$carat)*

**

原始分布和我用 copula 做的分布。图片作者。

*> cor(diamonds_num$carat, diamonds_num$price)
[1] 0.9215913012
> cor(carat, price)
[1] 0.9097580747*

如您所见，这种相关性在一定程度上得以保持。出现偏差的原因是，离散数据的相关性(泊松分布)与连续数据的相关性(高斯分布)不同。

我们现在可以看看建模部分。我现在不使用 carat，而是选择一种更直接的方法来确保训练和测试集的采样以及交叉验证不会碍事。下面是两个简单的线性回归。

*fit1<-lm(carat~price, data=diamonds_num)
fit2<-lm(carat~price, data=data.frame(cbind(carat,price)))
fit1;fit2

Call:
lm(formula = carat ~ price, data = diamonds_num)

Coefficients:
 (Intercept)         price  
0.3672972042  0.0001095002  

Call:
lm(formula = carat ~ price, data = data.frame(cbind(carat, price)))

Coefficients:
(Intercept)        price  
  18.887067     2.224347* 

很明显，系数不同，但这是意料之中的，因为我建立了不同的描述符。

*par(mfrow = c(2, 4))
plot(fit1);plot(fit2)*

模型拟合当然也有差异——上述四个图的矩阵来自原始数据，这些数据并不原始。后四个图来自我创建泊松分布的连接函数。假设正常数据是错误的，使用线性回归分析离散数据。图片作者。

一个更好的测试，仍然要记住泊松数据的线性回归是错误的，是对交互作用建模。还有什么比使用花键更好的方法呢！

*depth <- qpois(U[, 2], 30)
fit1<-lm(price~ns(carat,3)*ns(depth,3), data=diamonds_num)
fit2<-lm(price~ns(carat,3)*ns(depth,3), data=data.frame(cbind(carat,price, depth)))

> fit1

Call:
lm(formula = price ~ ns(carat, 3) * ns(depth, 3), data = diamonds_num)

Coefficients:
                (Intercept)                ns(carat, 3)1                ns(carat, 3)2  
                  2218.8390                   17936.5330                 -126840.8539  
              ns(carat, 3)3                ns(depth, 3)1                ns(depth, 3)2  
               -231168.4360                   -1091.1322                    1362.9843  
              ns(depth, 3)3  ns(carat, 3)1:ns(depth, 3)1  ns(carat, 3)2:ns(depth, 3)1  
                  7416.9316                    -329.7165                   68951.7911  
ns(carat, 3)3:ns(depth, 3)1  ns(carat, 3)1:ns(depth, 3)2  ns(carat, 3)2:ns(depth, 3)2  
                118173.3027                  -19264.2274                  263279.9813  
ns(carat, 3)3:ns(depth, 3)2  ns(carat, 3)1:ns(depth, 3)3  ns(carat, 3)2:ns(depth, 3)3  
                480346.8166                  -34626.0634                   14745.8991  
ns(carat, 3)3:ns(depth, 3)3  
                 73406.4320  

> fit2

Call:
lm(formula = price ~ ns(carat, 3) * ns(depth, 3), data = data.frame(cbind(carat, 
    price, depth)))

Coefficients:
                (Intercept)                ns(carat, 3)1                ns(carat, 3)2  
               -1.148528255                  8.348102184                 17.258451783  
              ns(carat, 3)3                ns(depth, 3)1                ns(depth, 3)2  
               15.909197993                 -0.242271773                 -1.327453563  
              ns(depth, 3)3  ns(carat, 3)1:ns(depth, 3)1  ns(carat, 3)2:ns(depth, 3)1  
               -0.862625027                 -0.393382870                 -0.090925281  
ns(carat, 3)3:ns(depth, 3)1  ns(carat, 3)1:ns(depth, 3)2  ns(carat, 3)2:ns(depth, 3)2  
               -0.029317928                 -0.006505344                  0.323633739  
ns(carat, 3)3:ns(depth, 3)2  ns(carat, 3)1:ns(depth, 3)3  ns(carat, 3)2:ns(depth, 3)3  
               -1.189341372                  0.188006534                 -0.109592519  
ns(carat, 3)3:ns(depth, 3)3  
               -0.984283018* 

当然，系数是不一样的，因为数据没有标准化，但让我们看看相互作用图。人们会假设，如果数据的本质得到维护，变量之间的关系也会得到维护。

*sjPlot::plot_model(fit1, type="pred")
sjPlot::plot_model(fit2, type="pred")*

**

上面你看到的是来自原始数据的变量之间的关系，下面你看到的是来自合成数据的关系。克拉和价格之间的关系似乎在某种程度上得以维持，但深度和价格肯定不是。我之所以选择样条曲线，是因为它们经常被使用，并且非常依赖于数据。因此，从原始到合成的转换中的错误转折很可能会被样条曲线拾取。图片作者。

检查转换和模型有效性的一个好方法是绘制原始数据并与合成数据进行比较，因为深度和价格根本不相关。两个模型都显示了联系。

我将使用 ggplot 并在原始数据的图形中拟合一条样条线。

*ggplot(diamonds_num, 
       aes(x=depth, y=price))+
  geom_point()+
  geom_smooth()+
  theme_bw()

ggplot(diamonds_num, 
       aes(x=carat, y=price))+
  geom_point()+
  geom_smooth()+
  theme_bw()*

**

清楚地显示了数据样条拟合的问题。如你所见，价格和克拉在这个二维层面上没有相关性，但样条曲线确实倾向于在中间跳动一点。在右边，我们看到克拉和价格被显示出来，在它们的顶部，一条样条线首先画出了一个清晰的关系。然后，它需要一个沉重的曲线来寻找它能找到的任何点。样条的伟大和危险显示在两个图中。图片作者。

以上图为原始数据。让我们也观察一下，如果我在合成数据上绘图会发生什么，现在合成数据具有与原始数据完全不同的分布属性(离散的，而不是原始的连续标度)。

*g1<-ggplot(diamonds_num, 
       aes(x=carat, y=price))+
  geom_point()+
  geom_smooth()+
  theme_bw()
g2<-ggplot(data.frame(cbind(carat,price, depth)), 
       aes(x=carat, y=price))+
  geom_point()+
  geom_smooth()+
  theme_bw()
gridExtra::grid.arrange(g1,g2,nrow=1)

g1<-ggplot(diamonds_num, 
           aes(x=carat, y=depth))+
  geom_point()+
  geom_smooth()+
  theme_bw()
g2<-ggplot(data.frame(cbind(carat,price, depth)), 
           aes(x=carat, y=depth))+
  geom_point()+
  geom_smooth()+
  theme_bw()
gridExtra::grid.arrange(g1,g2,nrow=1)* 

**

这里我们看到了两次观察同一关系的图。我们有克拉和价格，我们有克拉和深度。两者似乎都坚持原来的关系，即使来自不同的分布，但它并不完美。图片作者。

看上面的图，我们可以看到大部分的原始关系(或缺失)被保持。合成数据不会完美也就不足为奇了。此外，原始数据和合成数据的模型显示不同的系数也就不足为奇了。我已经制作了数据，所以它会有不同的描述特征，即使来自不同类型的分布，但仍然能够保持它的本质。

这篇博文只是一个简短的介绍，介绍了一种构建合成数据的方法，而且只是关于数值。使用 copulas，我们可以构建许多不同类型的合成数据。此外，我们还没有冒险进入深度学习，如 GANs，它主要用于建立合成数据。这个例子表明，我们不必走那么远。

如果有什么不对劲，请告诉我！

本文标签： 中文翻译一百一十博客TowardsDataScience