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

admin管理员组

文章数量:1645531

原文：TowardsDataScience

协议：CC BY-NC-SA 4.0

基因组学——终极数据科学

原文：https://towardsdatascience/genomics-the-ultimate-data-science-part-i-ec91683d1e18

思想和理论

基因组学——终极数据科学

从数据科学的角度概述基因组学，强调数据和数据解释在理解生物学中的重要作用

首先，让我们解决所有权要求。二十多年前，当我开始上大学并上第一堂生物课时，我第一次翻开坎贝尔的《生物学》(第四版)，生物学就成了我的主要兴趣。显然，我的观点是主观的，所以其他人可能会觉得金融、运筹学或推荐系统是数据科学的终极主题。

然而对我来说，生物学是一个完整的包裹。我对它了解得越多，它就越让我吃惊。它渗透了我对一切的想法和观点，从生孩子意味着什么，到对政治的看法，到理解不同社会如何互动，到理解当前的问题，如新型冠状病毒病毒的传播和气候变化。

这篇文章是介绍数据科学在基因组学中的作用的第一部分。它证明了我的标题声明，并介绍了基因组学中的一些基本概念，使其成为最终的数据科学。随后的帖子将对此进行扩展，每个帖子都讨论了一个特定的主题，并在相关时提供了代码示例。我的总体目标是与他人分享我对生物学和数据科学的热爱，尤其是它们的交集。希望这将帮助生物学家和数据科学家理解这些领域是如何交叉的。

我花了几年时间才真正理解，但是为什么基因组学是终极数据科学的关键来自生物学的三个不同特征，每一个都是‘啊哈！’属于我自己的时刻。

首先，我最终意识到的最重要的事情之一是，你必须抛弃所有关于生活会有多复杂的先入为主的观念。生活的复杂程度令人震惊。人们不得不怀疑它是否是无限的。这可以在所有尺度上看到，快速瞥一眼的生物数量可能会有助于说服你。

查尔斯·达尔文 1837 年绘制的进化树草图。图片来自维基共享资源。

第二个关键的洞见是理解为什么会这样。生物学中最重要的主题是进化。进化是生命如此复杂的原因。虽然达尔文早期关于进化的思想，包括他著名的草图，是基于合理的推理和一些逻辑跳跃，进化的研究现在是一门定量科学，也是基因组学是终极数据科学的另一个原因。

第三，由于进化是生物复杂性的驱动力，而且这种复杂性似乎是无限的，人们最终必须问自己，人类是否有能力理解它。人类没有任何理由能够理解生命是如何运作的。对我们来说，假设我们能够或者应该能够做到这一点是幼稚的。然而，人类擅长简化论，事实证明——非常幸运——地球上的生命有一个关键特征，它允许我们在生物学上站稳脚跟，形成基因组学领域，并创建最终的数据科学。

解开线——线是通往基因组学的大门

斯蒂夫·约翰森在 Unsplash 上拍照

进化的要素之一是信息存储和复制。生命不是靠魔法运行的，所以生物学家知道这些信息必须在细胞内的某个地方找到。

DNA 的结构，显示了组成其结构的 4 个核苷酸的例子。用户 Zephyris 的图片来自维基共享

每个人现在都知道生命以某种方式编码在 DNA 中，尽管如何做到这一点并不广为人知。众所周知，DNA 是由一系列被称为核苷酸的化学单体组成的。它们是 4 种相关的小化学物质，每种都用一个字母表示:A(腺嘌呤)、C(胞嘧啶)、G(鸟嘌呤)和 T(胸腺嘧啶)。现在这被认为是理所当然的，但这是 1953 年的一个重大发现，它确实令人惊讶，信息存储机制简单到足以被人类理解。这种存储机制打开了基因组学领域的大门，几乎所有生物学的复杂性都以某种方式植根于这 4 个字母的字符串中。

尽管 DNA 占据了媒体的大部分版面，并渗透到了流行文化中，但它实际上并不是第一个被认为是由可以用字符串表示的不同化学物质的线性排序组成的生命成分。这个奖项颁给了蛋白质。

细胞的大部分特性是由其组成蛋白决定的。蛋白质是细胞的主力，负责许多细胞活动。许多人可能只在包装背面的营养标签上听说过与饮食有关的蛋白质，但这掩盖了蛋白质是什么以及它们是如何运作的。

大卫·s·古德塞尔 对大肠杆菌横切面中蛋白质的艺术再现。 RCSB 蛋白质数据库。doi:10.2210/rcsb _ pdb/good sell-gallery-028。

蛋白质真的是原始的乐高积木，尽管更灵活和迷人。大肠杆菌蛋白质的手绘图像是 T2 的大卫·古德塞尔的许多美丽画作中的一幅，这些画作表现了蛋白质的活动。我推荐看看他的其他画作，他的月度分子系列，当然还有他令人惊叹的书生命的机器。这些肯定会改变你对蛋白质是什么以及它们能做什么的印象。该图像代表了大肠杆菌细胞的一部分，并成功地传达了细胞不是无定形的斑点，而是由许多称为蛋白质的精致小机器人组成，每个机器人都有一个特定的三维结构来决定其功能。它们一起工作来协调细胞的活动。一些执行代谢过程并产生能量为细胞提供动力，一些聚集形成驱动细胞运动的复杂马达，一些发挥结构作用以使细胞保持结构完整性，等等。

Digizyme 公司的 Evan Ingersoll & Gael McGill 创作了一个更近的细胞中蛋白质的艺术渲染图(如下所示)。这比 David S. Goodsell 拍摄的大肠杆菌截面图要复杂得多，但是这两张图片使用的方法是相同的。在这两种情况下，数十年来由无数其他人确定的蛋白质结构的 3D 坐标被获取、渲染并放置在场景中，以近似细胞中这些蛋白质角色的样子。这两个图像之间的差异类似于电影《白雪公主和七个小矮人》(1937)中的早期手绘动画与《玩具总动员》(1995)中的计算机生成动画之间的差异。也就是说，大卫·古德塞尔的艺术是手绘的，而现代计算机技术被用来生成下面的细胞景观。他们都以自己的方式引人注目。

Digizyme 公司的 Evan Ingersoll & Gael McGill 对细胞横截面中蛋白质的艺术再现。作者授予的许可。

对于大多数人来说，很容易在生活中没有真正思考过蛋白质的性质和功能，以及它们如何相互作用来执行细胞的功能，或者仅仅是它们的存在是多么多样和惊人。我希望上面的图片可以帮助说服你，但我必须首先提到上面的一点:你必须放弃所有关于生活可以有多复杂的先入为主的观念。上面的图像是细胞的美丽艺术渲染，但一个细胞中可能有数万种不同的蛋白质，它们不是静态的，它们经常在物理上相互作用以形成更大的功能结构，并且它们根据化学梯度、相邻的蛋白质、温度等不断反应和调整它们的行为。事实上，Digizyme 艺术家必须将场景中的蛋白质数量减少 80%以上，才能真正描绘出供人类消费的景观。

牛胰岛素的氨基酸组成。胰岛素由两条不同的氨基酸链组成，A 链和 B 链。每种氨基酸都有特定的化学结构，并用字母表中的一个字符来表示。图片作者。

当然，这些关于蛋白质的细节在 20 世纪上半叶都不为人知。众所周知，蛋白质执行特定的功能，但是确切地说，它们是如何执行这些功能的，以及它们采用了什么样的结构(如果有的话)却是未知的。正是在 20 世纪 40 年代末到 50 年代，蛋白质的大部分性质才被研究出来。由弗雷德·桑格领导的工作表明，这些神秘的蛋白质由一系列被称为氨基酸的化学物质组成，每个氨基酸都有一个单字母符号。这些氨基酸各有不同的结构，就像不同的 DNA 核苷酸有不同的结构一样。Sanger 指出牛胰岛素由两条独立的链组成，每条链都有特定的线性氨基酸序列。尽管在 1953 年发现 DNA 的双螺旋性质后胰岛素的全序列被发表，但是 A 链和 B 链的序列分别在 1952 年和 1951 年被确定。因此，在确定 DNA 可以表示为字符串之前不久，就确定了蛋白质可以表示为字符串的事实。当时，氨基酸排序的重要性还不为人知，因为不知道蛋白质是什么样子，也不知道它们是如何运作的。

在弗雷德·桑格表明蛋白质由线性氨基酸序列组成后不久，第一个蛋白质结构被确定，并于 1958 年发表。约翰·肯德鲁和他的同事确定了肌红蛋白(来自抹香鲸)的三维结构，提供了组成蛋白质的氨基酸线性序列决定其结构和功能的早期证据。

肌红蛋白的带状图，第一个发现的蛋白质结构。图片来自维基共享资源。

字符串和哈希——从 DNA 到蛋白质

伊曼纽尔·埃克斯特伦在 Unsplash 上拍摄的照片

到 20 世纪 50 年代中期，人们知道 DNA 和蛋白质分别是由 4 个不同的核苷酸(字符)和氨基酸(20 个字符)组成的聚合物(字符串)，但它们之间关系的性质实际上在 20 世纪 40 年代就已经确定了。1949 年，奥斯瓦尔德·艾弗里和他的同事们决定性地证明了 DNA——而不是蛋白质——是将一种无害细菌“转化”成一种可能导致疾病的细菌所必需的，并且这是一种可遗传的变化。蛋白质是不够的。这表明 DNA 负责可遗传的变化，并且首次证明细胞中不同的蛋白质集合可能在 DNA 中的某个地方以某种方式编码。当然，细节决定成败。

到了 20 世纪 60 年代，人们意识到 DNA 编码蛋白质，但研究人员还没有真正测序基因组的技术。然而，他们确实有能力生成基因图谱，这使他们能够大致定位基因内的一些突变，即使他们不知道基因的实际序列甚至长度。此外，他们能够突变 DNA，但由于他们实际上不能读取产生的基因组来识别突变，他们不得不通过检测它的影响，或者用生物学的说法，它的表型来推断发生了突变。因此，当时的策略是突变基因和基因组(当然不是人类的)，并确定基因或基因组中突变的大致位置。我们将在 python 中而不是在实验室中近似这种早期工作，以“重新发现”遗传密码的各个方面。

我们从随机突变一个基因开始，插入或删除一个碱基。这在基因的某处可能是这样的:ATCCG → AT G CCG 或 CCGAGT → CCGA T GT 用于插入，或 CG T CCG → CGCCG 或 ACG T AGT → ACGAGT 用于删除。这些只是假设的例子；实际的插入或删除操作(统称为插入)可以在基因的任何地方。在下面的图片中，每条水平线代表一个生物体的 20 个不同个体中长度未知的相同基因，每个圆圈代表随机插入一个碱基(红色带“+”)或随机删除原始基因中的一个碱基(灰色带“-”)。只显示保留基因功能的突变。在这种情况下，“功能性”意味着没有无义突变，即极有可能破坏基因功能的突变。这是一个过于简化的问题，但仍然是这个练习的一个很好的代理。为了节省空间，我省略了生成这些图像的大部分代码，但是大部分工作是由 MutatableSeqRecord 类完成的，该类继承了 biopython 的 SeqRecord 类，增加了向原始序列添加突变的能力:

这导致了下面的遗传图谱:

在 20 种不同的生物体中，沿着基因长度随机产生的可容忍的插入(+)或缺失(-)的位置，显示了沿着基因长度的突变的大致位置(0-开始，1-结束)。图片作者。

以上大致描述了 20 世纪中期遗传学是如何发展的。也就是说，已知基因由一长串有序的核苷酸组成(例如 ACTCGAGAGCCGATTA…或 ATGACACAGTTAGA)，可以延伸数百或数千个核苷酸(字符)，突变的大致位置可以沿其长度测量。除此之外，没有任何细节是已知的。然而，即使在上面的例子中，我们也可以获得一些启示。特别是，这些功能突变体中的所有突变都是在基因的右侧发现的，但不管是突变插入了还是缺失了一个碱基。这本身就很能说明问题，但可以有多种解释。让我们注意这一点:

• 靠近基因左侧的单碱基突变比靠近基因右侧的单碱基突变更有可能破坏基因的功能。

一个基因中也可以产生多个单碱基突变。我们将会看到我们可以从保留具有两个独立 indels 的功能基因中学到什么，使用与上面类似的代码，但为了节省空间，将它省去:

在 20 种不同的生物体中，沿着基因长度随机产生的 2 个容许的插入(+)或缺失(-)的位置，显示了沿着基因长度的突变的大致位置(0-开始，1-结束)。图片作者。

这是黄金。一个类似的模式出现在单一的 indels 中，但是有一些非常有趣的新信息。这些双重突变的功能基因中的 indels 通常仍然位于基因的右侧。然而，少数突变体(在这种情况下为 3/20)在基因左侧具有 indel，这些 indel 非常接近并且是相反的类型:一个插入和一个缺失。注意:突变体#2 具有相同的模式，但是 indel 位置如此接近以至于点是重叠的。我们在列表中添加了一个新的总体趋势:

• 如果靠近相反类型的 indel(插入与缺失相反),基因左侧的单碱基插入更可能被容忍

最后，让我们尝试 3 个独立的 indels，看看是否有更多我们可以提取的见解。根据上述趋势，我们已经知道这些可能会是什么样子。大多数功能突变体可能在基因的右侧有三个任何类型的突变位点，但少数在基因的左侧有两个相邻的相对突变位点，第三个突变位点在基因的右侧。让我们检查一下。这里有 20 个三重突变的功能基因:

在 20 种不同的生物体中，沿着基因长度随机产生的 3 个耐受的插入(+)或缺失(-)的位置，显示了沿着基因长度的突变的大致位置(0-开始，1-结束)。图片作者。

我们的假设被证实了。大多数突变是在基因的右侧发现的，但是在基因左侧突变的突变体与附近相反类型的 indel 配对。这并不能提供更多的见解，所以我们将通过生成三个突变体来缩小我们的关注范围，其中所有三个 indels 都位于左侧。我将继续使用蛮力方法来生成这些，因为它不会花很长时间来计算只有 20 个突变体。这将导致以下结果:

在 20 种不同的生物体中，沿着基因长度随机产生的 3 个容许的插入(+)或缺失(-)的位置，将突变体限制在所有 indels 都在基因前半部分的那些。x 轴显示沿着基因长度的突变的大致位置(0-开始，1-结束)。图片作者。

如你所见，功能性三重突变体中的 indel 位置有一个引人注目的模式，这是这个谜题的最后一块。在所有情况下，三个独立体紧密地聚集在一起，并且总是属于同一类型。我们将在列表中添加另一个观察结果:

• 如果在一个功能性三重突变基因的左边发现三个 indels，它们紧密地聚集在一起，属于同一类型

总之，我们有以下观察结果:

1. 靠近基因左侧的单碱基突变比靠近基因右侧的单碱基突变更有可能破坏基因的功能。
2. 如果靠近相反类型的 indel(插入与缺失相反),基因左侧的单碱基插入更可能被容忍
3. 如果在一个功能性三重突变基因的左侧发现三个 indels，它们紧密聚集在一起，属于同一类型。

为了了解这里发生了什么，我们将使用一个类比。假设你的朋友给你发了一条短信，我们会认为这些短信在通信网络上很容易被破坏。消息存储为内存中连续的字节数组，这些数组通常会因内存连续块中某个位置的随机添加或删除而被破坏。在这种情况下，我们不使用前导位，因为在 ASCII 字符编码中它总是零，所以我们使用七位来表示每个字节。考虑到这一点，你看着你的手机阅读来电信息，内容如下:

Meet me on!)kKgICrACe_k]HAgKmK\AChAiQJA[CYXA]KCdAiQJAM__HAG_ke

这可不好。这句话几乎完全没有意义。你应该在某个时候去见你的朋友，但是你不知道时间和地点。很明显，该消息在靠近消息开头的某处被插入或删除的位所破坏。如果这发生在消息的结尾就好了。那样的话，你也许就能破译它了。照目前的情况来看，这条消息基本上没有什么作用。

假设您真的很不幸，消息被七个不同的位插入或删除破坏了。你现在绝对不走运。信息内容如下:

[KPA[JA_9SWOeK?W;OmK\AChAiQJA[CYXA]KAr the food g_ke

但是，如果这七种突变都是同一类型，并且发生在很近的地方，会怎么样呢？

Meet me on Tuesday a|U.{ seven at the mall near the food court

我们现在可以读了。虽然插入了七位，但消息仍然清晰可辨，因为只有插入附近的位受到影响，有效地将所有下游位向下游移动了一个字节的位。你可能还会想，在购物中心见面不再是任何人真正会做的事情，但不管怎样，希望这个例子能达到它的目的。

这与上述基因突变体中发生的现象完全相同，但它们的“字节”在计算机中相当于 3 个碱基，而不是 8 位。通过在附近插入(或删除)3 个碱基的倍数，只有蛋白质的局部受到影响，而蛋白质的其余部分保持完整。让我们停下来惊叹，地球上所有形式的生命都使用与人类数字存储信息完全相同的信息存储方法，这是多么令人惊讶。进化和人类——他们自己也是进化的产物——都发现了同样的储存信息的方法。

因此，概括地说，我们都知道数字信息可以存储为字节，每个字节由 2 字符字母表(位)中的 8 个字符组成，用于存储 ASCII 字符。这些 ASCII 字符被连续解码以拼写单词和组成句子。这些句子的作用是传达意思。类似地，生物信息可以存储为密码子而不是字节，每个密码子由 4 字符字母表(而不是 2 字符 0/1 字母表)中的 3 个字符(而不是 8 个)组成，以编码氨基酸。这些氨基酸被连续解码，以确定氨基酸的线性序列，该序列决定蛋白质的结构，从而决定其功能。这些蛋白质可能有许多功能。

地球生命中英文的数字位编码(ASCII)和通过 DNA 的蛋白质化学编码的比较。注意，尽管在 ASCII 编码中存在标点符号和其他非英语字符，但为了简单起见，这里只考虑了英语字母表的 26 个字母。图片作者。

这里唯一缺少的一块拼图是，我们知道将字节映射到字符的 ASCII 字符编码，因为人类开发了编码方案，但科学家在 19 世纪中期不知道遗传编码方案。他们只知道它的存在。不幸的是，他们没有谷歌或我们目前的生物学知识库。他们必须在实验室里解决。

在本文的第二部分，我将讨论遗传密码是如何被破译的，并继续详述由此产生的信息到底是做什么的*。*我会继续将它与英语进行比较，因为这样更直观，而且有许多相似之处。这样做，我将为理解基因组学作为一门数据科学奠定基础，并阐述为什么我相信它是终极数据科学。

请继续关注第二部分…

张量流概率简介:分布对象

原文：https://towardsdatascience/gentle-introduction-to-tensorflow-probability-distribution-objects-1bb6165abee1

概率深度学习

介绍

本文属于“概率深度学习”系列。这个每周系列涵盖了深度学习的概率方法。主要目标是扩展深度学习模型，以量化不确定性，即知道他们不知道的东西。

我们使用张量流和张量流概率(TFP)开发我们的模型。TFP 是构建在 TensorFlow 之上的 Python 库。我们将从可以在 TFP 中找到的基本对象开始，并了解如何操作它们。我们将在接下来的几周内逐步增加复杂性，并将我们的概率模型与现代硬件(如 GPU)上的深度学习相结合。

迄今发表的文章:

1. 张量流概率简介:分布对象
2. 张量流概率简介:可训练参数
3. 张量流概率中从零开始的最大似然估计
4. tensor flow 中从头开始的概率线性回归
5. 使用 Tensorflow 进行概率回归与确定性回归
6. Frequentist 与 Tensorflow 的贝叶斯统计

图 1:我们在这个系列中结合了两个世界:概率模型和深度学习(来源)

像往常一样，代码可以在我的 GitHub 中找到。

分发对象

单变量分布

分布对象捕捉概率分布的基本操作。让我们从最简单的形式开始——单变量分布。顾名思义，这些分布只有一个随机变量。高斯分布是完全由其均值和标准差定义的连续概率分布。它的标准型是μ = 0，σ = 1 的特例。

normal = tfd.Normal(loc=0, scale=1)
normal

<tfp.distributions.Normal 'Normal' batch_shape=[] event_shape=[] dtype=float32>

注意属性batch_shape和event_shape。event_shape属性捕捉随机变量的维度。因为我们定义了一个单变量分布，所以event_shape是空的。batch_shape表示我们存储在对象中的不同独立分布对象的数量。在这种情况下，我们只存储一个发行版，因此它是空的。

我们可以用这个分布对象做什么？例如，我们可以从中取样。

# Draw one sample from the normal distribution
normal.sample()

<tf.Tensor: shape=(), dtype=float32, numpy=1.5117241>

# Draw 3 samples from the normal distribution
normal.sample(3)

<tf.Tensor: shape=(3,), dtype=float32, numpy=array([-1.2282027 , -0.01802123,  0.2748567 ], dtype=float32)>

在连续随机变量的情况下，我们还可以评估给定输入的概率密度函数(PDF)。

normal.prob(0.2)

<tf.Tensor: shape=(), dtype=float32, numpy=0.3910427>

当我们开始实施机器学习和深度学习算法时，我们经常会发现自己面临一个共同的问题:乘以概率并不有趣。整个产品开始迅速接近零，你将用尽精度来存储这么小的数字。

解决这个问题的常见方法是用对数概率代替概率。首先，通过应用对数变换，我们将域从[0，1]更改为(-∞，0)。这是相关的，因为我们可以在对数标度中存储大的负数，而不是必须存储非常小的数字(我们会很快用完精度)。第二，对数函数是一个单调递增的函数，这确保我们可以像比较原始对数概率一样比较两个转换后的对数概率。最后，对数转换将所有乘法运算转换为加法运算。因此，概率不再相乘，而是将对数概率相加。

在 TFP 中，处理对数概率是很简单的，我们只需要调用一个不同的方法。

normal.log_prob(0.2)

<tf.Tensor: shape=(), dtype=float32, numpy=-0.9389385>

注意，上面只是由prob方法返回的值的自然对数。

np.log(normal.prob(0.2))

-0.9389385

最后，我们可以绘制一个直方图，它近似于分布的密度。

sns.histplot(data=normal.sample(10000).numpy(), kde=True);

图 2:标准高斯分布的直方图和核密度估计。

现在，让我们探索如何在单个对象中存储一批发行版。让我们创建两个单变量正态分布，一个μ=0，σ=1，另一个μ=1，σ=1。

normal_batch = tfd.Normal([0,1],[1,2])
normal_batch

<tfp.distributions.Normal 'Normal' batch_shape=[2] event_shape=[] dtype=float32>

请注意，这两个分布需要是相同的类型(在这种情况下，它们都是高斯分布),并且它们被认为是独立的，即这不是创建多元分布的方法。在上面的输出中，您可以看到batch_shape现在是预期的 2。

我们可以很容易地从两种分布中抽取样本。

normal_batch.sample(3)

<tf.Tensor: shape=(3, 2), dtype=float32, numpy=
array([[-1.5531085 ,  1.4462973 ],
       [-0.5346463 ,  0.63747466],
       [-2.2433918 , -0.4113649 ]], dtype=float32)>

我们样品的形状不同。现在是(3，2)，因为我们从每个正态分布中抽取 3 个样本。我们可以绘制这两种分布，只是为了确保它们没有任何关联。

samples = normal_batch.sample(10000).numpy()
ax = sns.scatterplot(x = samples[:,0], y = samples[:,1])
ax.set(ylim=(-10, 10), xlim=(-10, 10), xlabel='N(0,1)', ylabel='N(1,2)');

图 3:两个独立高斯分布的样本。

同样，我们可以得到两种分布的 pdf 值。

normal_batch.prob([1,2])

<tf.Tensor: shape=(2,), dtype=float32, numpy=array([0.24197073, 0.17603266], dtype=float32)>

我们得到一个返回的形状为(2)的张量，因为我们得到了两个不同点的分布的 PDF 值。

为了更好地理解形状的行为，让我们增加batch_shape的等级。

normal_batch_D = tfd.Normal([[[0., 2],
                             [0.5, 1],
                             [5, 2]]], scale=1)
normal_batch_D

<tfp.distributions.Normal 'Normal' batch_shape=[1, 3, 2] event_shape=[] dtype=float32>

我们现在有 6 个不同的高斯分布存储在同一个对象中。不运行代码能说出normal_batch_D[0,2,1] 的均值和标准差是多少吗？

像往常一样，我们可以从中取样。

normal_batch_D.sample(10)

<tf.Tensor: shape=(10, 1, 3, 2), dtype=float32, numpy=
array([[[[-0.94896364,  2.1813042 ],
         [ 0.14763275,  0.22235268],
         [ 4.8377185 , -0.19643283]]],
       [[[ 0.6483533 ,  2.3491006 ],
         [-0.11504221,  0.13614637],
         [ 5.2141023 ,  1.9243499 ]]],
       [[[ 0.14039962,  1.5450974 ],
         [ 1.134828  ,  2.2807612 ],
         [ 5.892858  ,  0.6847892 ]]],
       [[[ 1.3779826 ,  2.0819554 ],
         [ 1.0930698 ,  0.5755873 ],
         [ 4.71762   ,  2.248595  ]]],
       [[[ 0.21968068,  1.2137487 ],
         [ 1.3834007 , -0.17327452],
         [ 5.6132197 ,  2.4669297 ]]],
       [[[-0.7178315 ,  1.1999301 ],
         [-0.19561946,  0.14512819],
         [ 3.7195773 ,  1.3497605 ]]],
       [[[ 0.03890136,  2.9174664 ],
         [ 0.37707615, -1.6925879 ],
         [ 4.0377812 ,  2.6682882 ]]],
       [[[ 1.4869312 ,  2.2919848 ],
         [ 1.1833754 ,  0.78834504],
         [ 4.746928  ,  2.398845  ]]],
       [[[-2.2711177 ,  1.9751831 ],
         [ 2.855303  , -0.51687765],
         [ 5.6308627 ,  0.96069396]]],
       [[[-0.5072157 ,  1.7689023 ],
         [ 0.67927694,  0.30969065],
         [ 3.8056169 ,  3.4717598 ]]]], dtype=float32)>

我们得到一个形状为(10，1，3，2)的张量，这意味着我们有 10 个样本(第一维)用于 6 个高斯分布中的每一个。

多元分布

现在，是时候探索我们的分配对象的event_shape属性了。创建多元分布有多种方法，让我们从最简单的一种开始。我们可以定义一个二维高斯分布，并且不包括两个维度之间的任何相关性，这意味着协方差矩阵的非对角项为 0。

mv_normal = tfd.MultivariateNormalDiag(loc=[0, 1], scale_diag=[1., 2])
mv_normal

<tfp.distributions.MultivariateNormalDiag 'MultivariateNormalDiag' batch_shape=[] event_shape=[2] dtype=float32>

最后，我们看到填充了一个event_shape，在本例中值为 2。正如我们所说的，这是我们上面定义的随机变量的维数。

是时候从中抽取样本，并查看与多个批处理分布相比的差异了。

mv_normal.sample(3)

<tf.Tensor: shape=(3, 2), dtype=float32, numpy=
array([[ 1.4743712 , -0.77387524],
       [ 1.7382311 ,  2.747313  ],
       [-1.3609515 ,  4.8874683 ]], dtype=float32)>

形状与我们从两个分批分布中采样 3 次的示例相同。由于我们已经定义了一个具有非对角协方差矩阵的多元分布(维度是独立的),我们得到了类似的结果——比较图 4 和图 3。

图 4:具有独立维度的多元高斯分布的样本。

是时候将这两个概念结合起来，定义一批多元分布了。

normal_diag_batch = tfd.MultivariateNormalDiag(loc=[[0,0], [1,1], [0,0]],
                                              scale_diag=[[1,2], [1,1], [2,10]])
normal_diag_batch

<tfp.distributions.MultivariateNormalDiag 'MultivariateNormalDiag' batch_shape=[3] event_shape=[2] dtype=float32>

我们现在可以看到，batch_shape和event_shape都被填充了，值分别为 3 和 2。这意味着我们已经创建了 3 个二维高斯分布。

samples = normal_diag_batch.sample(10000).numpy()
samples.shape

(10000, 3, 2)

当从上面的对象中采样时，我们得到一个具有以下维度的输出:(number_samples, batch_shape, event_shape)。

让我们探索计算这个对象的对数概率。您希望看到的输出是什么样的？

normal_diag_batch.log_prob(samples)

<tf.Tensor: shape=(10000, 3), dtype=float32, numpy=
array([[-2.5595174, -1.8528149, -5.7963095],
       [-3.219818 , -2.9775417, -9.757326 ],
       [-5.3475537, -1.8487425, -6.300317 ],
       ...,
       [-3.3621025, -3.5017567, -6.5536766],
       [-2.7969153, -2.8572762, -5.0501986],
       [-3.2498784, -2.2819252, -7.9784765]], dtype=float32)>

因为我们有 3 个分布(尽管它们是二维的)，所以我们得到一个输出(number_samples, batch_shape)。为什么？

我们可以绘制 3 个分布的样本，以比较协方差矩阵中不同对角线值的影响程度。

plt_sample_batch = normal_diag_batch.sample(10000).numpy()

fig, axs = (plt.subplots(1, 3, sharex=True, sharey=True, figsize=(10, 3)))
titles = ['cov_diag=[1,2]', 'cov_diag=[1,1]', f'cov_diag=[2,10]']

for i, (ax, title) in enumerate(zip(axs, titles)):
    samples = plt_sample_batch[:,i,:]
    sns.scatterplot(x=samples[:,0], y=samples[:,1], ax=ax)
    ax.set_title(title)
    axs[i].set_ylim(-25, 25)
    axs[i].set_xlim(-25, 25)
plt.show()

图 5:具有不同对角协方差矩阵的 3 个多元高斯分布的样本。

请注意，在此图中，任何分布都没有出现相关性。即便如此，第三分布在 y 轴方向上被极度拉长。这是因为对角线上的第二个值比第一个值大，第二个维度的标准偏差比第一个维度的标准偏差大得多，因此我们看到了 y 方向上的分布。

结论

本文介绍了理解 TFP 中的核心对象——分布对象的第一步。我们首先定义单变量分布，操作对象的batch_shape属性，并测试不同的方法:sample、prob和log_prob。接下来，我们增加了随机变量的维数，并引入了多元分布。在这种情况下，我们探索了处理对象的batch_shape和event_shape属性的方法，允许我们在单个对象中存储多个多元分布。

下周，我们将探讨如何训练这些分布对象的参数。到时候见！

保持联系: LinkedIn

参考资料和材料

[1] — Coursera:深度学习专业化

[2] — Coursera:深度学习的 tensor flow 2专业化

[3] — 张量流概率指南和教程

[4] — TensorFlow 博客中的 TensorFlow 概率帖子

张量流概率简介——可训练参数

原文：https://towardsdatascience/gentle-introduction-to-tensorflow-probability-trainable-parameters-5098ea4fed15

概率深度学习

介绍

本文属于“概率深度学习”系列。这个每周系列涵盖了深度学习的概率方法。主要目标是扩展深度学习模型，以量化不确定性，即知道他们不知道的东西。

我们使用张量流和张量流概率(TFP)开发我们的模型。TFP 是构建在 TensorFlow 之上的 Python 库。我们将从可以在 TFP 中找到的基本对象开始，并了解如何操作它们。我们将在接下来的几周内逐步增加复杂性，并将我们的概率模型与现代硬件(如 GPU)上的深度学习相结合。

迄今发表的文章:

1. 张量流概率简介:分布对象
2. 张量流概率简介:可训练参数
3. 张量流概率中从零开始的最大似然估计
4. tensor flow 中从头开始的概率线性回归
5. 使用 Tensorflow 进行概率回归与确定性回归
6. Frequentist 与 Tensorflow 的贝叶斯统计

图 1:我们在这个系列中结合了两个世界:概率模型和深度学习(来源)

像往常一样，代码可以在我的 GitHub 上找到。

分发对象

在上一篇文章中，我们看到了如何操作 TFP 分布对象。记住分布对象捕捉概率分布的基本操作。我们从单变量分布开始，即只有一个随机变量的分布。然后，我们将我们的理解扩展到如何用分布对象属性来表示多元分布。我们保持简单，因为我们定义了一个 2 维高斯分布，并且没有包括两个维度之间的任何相关性。要回忆的最重要的属性是batch_shape和event_shape。如果您对它们还不满意，请查看我以前的文章。我们将在本系列中广泛使用它们。

在介绍可训练的分布参数之前，我们将再介绍一个关于分布对象的概念。

独立分布

有些情况下，我们希望将事件空间上的一批独立分布解释为事件空间乘积上的单个联合分布。这影响了我们处理batch_shape和event_shape属性的方式。当我们开始构建一些众所周知的算法(如朴素贝叶斯分类器)时，独立分布将非常有用。原因是在朴素贝叶斯的情况下，给定一个类标签，特征是独立的。

为了说明，让我们定义两个正态分布。

第一个是以下形式的多元常态:

为了定义第一个，我们将像以前一样使用MultivariateNormalDiag，因为，再一次，它们之间的维度是不相关的。

mv_normal = tfd.MultivariateNormalDiag(loc=[0, 1], scale_diag=[1,2])

<tfp.distributions.MultivariateNormalDiag 'MultivariateNormalDiag' batch_shape=[] event_shape=[2] dtype=float32>

我们对形状属性越来越满意，因此毫不奇怪我们的event_shape为 2。

像往常一样，我们可以计算对数概率:

mv_normal.log_prob([0.2, 1.5])

<tf.Tensor: shape=(), dtype=float32, numpy=-2.5822742>

我们得到一个单一的值，因为我们有一个单一的分布，即使它是多维的。

让我们从独立的多元高斯分布中取样，并绘制联合分布图。我们以前做过类似的事情。

samples = mv_normal.sample(10000).numpy()
x1 = samples[:,0]
x2 = samples[:,1]
sns.jointplot(x = x1, y = x2, kind='kde', xlim=[-6, 7], ylim=[-6, 7]);

图 2:上面定义的多元高斯分布的近似联合分布图。还显示了单个维度的单变量图。

正如所料，多元高斯分布的维数之间没有相关性。

表示批量高斯分布对象的时间。

locs = [0,1]
scales = [1, 2]

batched_normal = tfd.Normal(loc=locs, scale=scales)
batched_normal

<tfp.distributions.Normal 'Normal' batch_shape=[2] event_shape=[] dtype=float32>

注意batch_shape等于 2。

batched_normal.log_prob([0.2, 1.5])

<tf.Tensor: shape=(2,), dtype=float32, numpy=array([-0.9389385, -1.6433357], dtype=float32)>

由于我们在同一个对象中存储了两个独立的分布，因此计算对数概率会产生两个值。

我们可以绘制两个单变量分布的概率密度函数(PDF)。

x = np.linspace(-4, 4, 10000)
densities = batched_normal.prob(np.repeat(x[:, np.newaxis], 2, axis=1))

sns.lineplot(x=x, y=densities[:, 0], label=f'loc={locs[0]}, scale={scales[0]}')
sns.lineplot(x=x, y=densities[:, 1], label=f'loc={locs[1]}, scale={scales[1]}')
plt.ylabel('Probability Density')
plt.xlabel('Value')
plt.legend()
plt.show()

图 3:作为单个分布对象批处理的两个单变量高斯分布的 PDF。

让我们总结一下上面的内容，这样我们就可以介绍独立分布对象了。很明显，当第一个分布对象返回一个对数概率时，第二个返回 2。不同之处在于，我们传递给第一个的数组被解释为一个二维随机变量的单一实现。在第二种情况下，数组被解释为每个随机变量(批次)的不同输入。

为了帮助我们理解什么是独立分布以及它是如何有帮助的，我们来玩一些概率术语:

• 独立分布是我们从单变量分布转移到单个多变量分布的一种简化方法；
• 独立分布允许我们从单个随机变量的几个分布转移到一组随机变量的联合分布；
• 独立分布提供了从几个批量分布转移到一个多维分布的能力；
• 独立分布是一个接口，将我们想要吸收的任何维度吸收到事件维度；
• 最后，更实用和 TFP 的描述方式——独立分布是一种将分布的batch_shape 维度移动到新分布对象的event_shape的方式。

希望用这么多不同的方式来描述它能使所有这些概率概念以及它们被转化为 TFP 抽象的方式更加清晰。

是时候应用理论概念并查看实际实现了。

independent_normal = tfd.Independent(batched_normal, reinterpreted_batch_ndims=1)
independent_normal

<tfp.distributions.Independent 'IndependentNormal' batch_shape=[] event_shape=[2] dtype=float32>

批量高斯分布现在是一个IndependentNormal分布对象，它是一个独立的多元高斯分布，如我们上面定义的。我们可以通过 2 的event_shape看出来。类似地，对数概率现在应该产生一个值。

<tf.Tensor: shape=(), dtype=float32, numpy=-2.5822742>

最后，让我们将独立高斯分布图与上面显示的图进行比较。

samples = independent_normal.sample(10000).numpy()
x1 = samples[:,0]
x2 = samples[:,1]
sns.jointplot(x = x1, y = x2, kind='kde', space=0, color='b', xlim=[-6, 7], ylim=[-6, 7]);

图 4:独立高斯分布对象的近似联合分布图。还显示了单个维度的单变量图。

可训练参数

变量

现在我们知道了什么是张量流概率对象，是时候了解如何为这些分布训练参数了。这是我们开始应用我们所学的知识和构建算法时所缺少的联系。

在 TensorFlow 中，Variable对象是我们用来捕捉深度学习模型的参数值的东西。这些对象在训练期间通过例如应用从损失函数和数据获得的梯度来更新。

我们来定义一个。注意，要创建一个新变量，我们必须提供一个初始值。

init_vals = tf.constant([[1.0, 2.0], [3.0, 4.0]])
new_variable = tf.Variable(init_vals)
new_variable

<tf.Variable 'Variable:0' shape=(2, 2) dtype=float32, numpy=
array([[1., 2.],
       [3., 4.]], dtype=float32)>

A Variable非常类似于张量。它们具有相似的属性，例如形状和数据类型，以及方法/操作，例如导出到 NumPy。虽然他们有一些不同，例如，他们不能被重塑。

print("shape: ", new_variable.shape)
print("dType: ", new_variable.dtype)
print("as NumPy: ", new_variable.numpy())
print("Index of highest value:", tf.math.argmax(new_variable))

shape:  (2, 2)
dType:  <dtype: 'float32'>
as NumPy:  [[1\. 2.]
 [3\. 4.]]
Index of highest value: tf.Tensor([1 1], shape=(2,), dtype=int64)

请注意，如果出于某种原因，您不希望在训练过程中对某个变量进行微分，您可以使用参数trainable来定义它。

variable_not_diff = tf.Variable(1, trainable=False)

<tf.Variable 'Variable:0' shape=() dtype=int32, numpy=1>

总之，通常，我们希望我们的变量是可微的。张量流允许自动微分，这是用于训练神经网络的反向传播算法的基础。

我们将使用一个 API 来完成自动微分— tf.GradientTape。连接回Variable对象，这个 API 让我们能够计算操作相对于输入的梯度，即一个或多个Variable对象。

让我们用tf.GradientTape API 和Variable对象做一个简单的例子。

x = tf.Variable(3.0)

with tf.GradientTape() as tape:
    y = x**2

一旦我们在tf.GradientTape上下文中定义了一个操作，我们就可以调用gradient方法并传递损失和输入变量。

dy_dx = tape.gradient(y, x)
dy_dx.numpy()

6.0

是时候把这些概念应用到我们的问题上了。回想一下，我们对学习分布的参数感兴趣。

normal = tfd.Normal(loc=tf.Variable(0., name='loc'), scale=5)
normal.trainable_variables

(<tf.Variable 'loc:0' shape=() dtype=float32, numpy=0.0>,)

在这种情况下，上面定义的高斯分布的均值不再是一个简单的值，而是一个可以学习的Variable对象。

对于训练过程，最大似然是深度学习模型中通常的疑点。简而言之，我们正在寻找使数据概率最大化的模型参数。

连续随机变量的概率密度函数粗略地表示样本取特定值的概率。我们将把这个函数表示为𝑃(𝑥|𝜃)，其中 𝑥 是样本值， 𝜃 是描述概率分布的参数:

tfd.Normal(0, 1).prob(0)

<tf.Tensor: shape=(), dtype=float32, numpy=0.3989423>

这可能看起来很奇怪，但事实上，我们已经计算高斯分布的 PDF 有一段时间了，所以这里没有什么特别新的东西。

为了完成对训练参数的介绍，让我们将这个概念与我们在上面分享的独立分布对象联系起来。当从同一个分布中独立地抽取多个样本时(我们通常假设)，样本值 𝑥 1、…、 𝑥𝑛 的 PDF 是每个个体 𝑥𝑖 的 pdf 的乘积。我们可以把它写成:

希望您能在上面的定义中看到这两个概念是如何重叠的。

结论

本文继续探讨 TFP 中的分布对象，但这次我们更进一步，将它们与 TensforFlow 中的Variable对象联系起来。我们从定义什么是独立分布以及它如何帮助我们定义独立联合概率开始。它允许我们从单变量分布转移到独立的多变量分布，从而将我们想要的任何维度吸收到事件维度中。接下来，我们介绍了Variable对象，以及如何利用自动微分来区分它们。有了这些知识，我们将它们与 TFP 中的一个分发对象结合使用。最后，我们讨论了最大似然法，以及当我们从同一个分布中独立抽样时，它与独立联合分布的关系。

下周，我们将更详细地探讨发行版的培训过程。到时候见！

保持联系: LinkedIn

参考资料和材料

[1] — Coursera:深度学习专业化

[2] — Coursera:深度学习的 tensor flow 2专业化

[3] — 张量流概率指南和教程

[4] — TensorFlow 博客中的 TensorFlow 概率帖子

地球同步轨道提升试验 I:理解差异中的差异

原文：https://towardsdatascience/geo-lift-experiments-i-understanding-difference-in-differences-5f35e41a2fdb

这是 2 部分系列的第 1 部分，着眼于产品环境中的 geo lift 实验:

1. 理解差异中的差异
2. Spotify Blend 案例研究

在这一部分，我们将探讨为什么以及何时进行 geo lift 实验。我们还将探索差异中的潜在差异因果推断方法。

由凯尔·格伦在 Unsplash 上拍摄

最近 AB 测试已经成为产品开发团队的黄金标准。在产品上运行实验有助于团队理解每个变化对他们的关键度量的增量效应，并逐渐改进他们的产品。然而，有时随机地向一些用户推出一个特性，而不向另一些用户推出，会导致不准确的结果。这可能是由于网络效应、库存溢出，或者更普遍的情况，即稳定单位处理值假设 (SUTVA)不成立。这是假设一个随机单元参与产品变更的方式独立于任何其他单元参与相同产品变更的方式。随机化单元是我们对每个变量进行随机分配的层次。在大多数情况下，这是个人用户，但也可以是设备、会话或社交群等。让我们看一些例子来更好地理解这个概念。

特别是对于社交产品，许多功能需要双向参与才能发挥作用。随机地向用户显示一个特征变化可能导致那些用户中的一些没有适当的机会参与其中，因为不是他们所有的朋友都可以访问相同的特征变化。因此，测试中衡量的成功标准可能无法真正代表将变更推广到所有用户的效果。

同样，违反 SUTVA 是电子商务和旅游产品在试验决定向用户显示哪些项目的算法时面临的问题之一，如个性化和搜索算法。例如，在 Airbnb 上，当一种算法向用户显示他们继续预订的房间时，该房间不再向任何被分配了不同算法的用户显示。因此，这两个算法的参与度指标是不一致的，因为其中一个算法的参与度并不独立于另一个算法。

通过改变实验的设计，我们有多种方法可以控制这些。例如，在社交网络的情况下，我们可以将网络分成彼此交互频繁的用户群，并在群集级别执行随机分配。然而，实现更复杂的实验设计会有很大的开发成本。因此，有时我们可能没有资源来进行实验。

这个问题的一个潜在解决方案是以 geo lift 实验的形式运行一个准实验。这个想法是向特定地区的所有用户发布特性更改。然后，将该地区的关键成功指标与未发布功能变更的控制地区的成功指标进行比较。然后，我们使用这种比较来理解该特性对我们的关键成功指标的因果影响。值得注意的是，这只能作为一种解决方案，前提是假设将实验限制在特定的地理区域代表了向您的整个用户群发布该功能。以社交产品为例，这只能作为一种解决方案，前提是假设对于大多数用户来说，他们在产品上与之交互的所有其他用户都来自相同的地理位置。

这种方法是基于差异的因果推断方法。使用这种方法测量的关键成功指标的提升不像运行一个设计良好的实验那样是因果影响的有力指标，但它更足智多谋，并且(在大多数情况下)设置和推出更快。

这种方法也是控制基于时间的混杂因素的好方法。例如，在进行实验的同时开展普遍获得活动或销售。或者，对于样本量较小的公司，基于时间的混杂因素可能只是自然变异。然而，该方法的进一步假设是，任何基于时间的混杂因素都适用于两个地理区域。例如，新的普遍获得营销活动必须在控制和测试地区同时进行。为了理解为什么需要这样，让我们更深入地研究一下这种方法背后的因果推理技术。

差异中的差异

差异中的差异方法使用两个相似的用户群，并在功能变更发布前后对关键成功指标进行建模。我们的一个用户群被作为对照，在准实验期间没有接收到特征变化。另一个用户部分进行处理，并接收特征变化。对于每项成功指标，我们建立了以下模型:

• 控制预发布(a)
• 控制后释放(x)
• 治疗前释放(b
• 治疗后释放(y

两个地区之间的发布前差异由b — a给出，两个地区之间的发布后差异由y — x给出。然后，我们可以将由于特征变化导致的整体抬升建模为这些差异之间的差异，因此该方法得名。

uplift = post_release_diff - pre_release_diff
       = (y - x) - (b - a)

此图显示了不同方法的不同之处。geo 之间的发布后差异等于发布前差异+由于特征变化而产生的隆起。(图片由作者提供)

现在我们对 geo lift 实验的工作原理有了更好的了解，让我们来看看在第 2 部分中使用这种方法的真实产品示例。我们将使用一个名为 PyMC 的 python 包，使用我们观察到的数据对我们的成功指标进行建模。它使用贝叶斯统计和蒙特卡罗马尔可夫链 (MCMC)方法来建模指标。我们将把这个应用到 Spotify Blend 功能中。

感谢阅读这篇文章！我希望它能帮助你更好地理解什么时候使用这种因果推理方法是有用的。

如果你喜欢阅读我的文章，愿意支持我的写作，并且正在考虑订阅一个媒体，请随时使用我下面的推荐链接。我会从你的订阅费中提成。

https://medium/@kaushsk12/membership

Geo Lift 实验 II: Spotify Blend 案例研究

原文：https://towardsdatascience/geo-lift-experiments-ii-spotify-blend-case-study-476a81099744

这是两部分系列的第 2 部分，介绍产品环境中的 geo lift 实验:

1. 理解差异中的差异
2. Spotify Blend 案例研究

在这一部分，我们将把我们在第一部分中探索的内容应用到案例研究中。让我们深入案例研究。

照片由迪帕克·乔杜里在 Unsplash 拍摄

产品背景:Spotify Blend

假设我们是 Spotify 播放列表体验团队的项目经理。对于任何不熟悉这个功能的人来说，这是一个允许用户与他们选择的朋友创建播放列表的功能。然后，播放列表会自动生成用户及其朋友都喜欢的歌曲。它还向用户显示他们的音乐品味和朋友的音乐品味之间的匹配分数。当用户收听新音乐时，该品味匹配分数改变，并且播放列表提供的推荐变得更好。

作为这个特性的项目经理，让我们假设我们发现用户在设置混合的体验上有一些摩擦。邀请朋友加入 blend 时，会向用户显示一个本机共享模式，以便与朋友共享邀请链接(如下所示)。当朋友们点击与他们分享的链接时，他们就会被添加到混合中。让我们假设，查看数据，我们看到在用户点击邀请按钮和他们的朋友成功接受他们的邀请之间有相当大的落差。我们没有太多的数据来说明这到底是为什么，因为这可能是由于用户没有通过共享模式发送邀请。也可能是朋友没有回复邀请信息。

点击 Spotify Blend 功能上的邀请按钮时的原生共享模式。(图片由作者提供)

假设

我们假设，通过在 Spotify 应用程序中保持这种体验，用户更有可能成功邀请他们的朋友，并创建一个混合。我们希望更改这部分功能，允许用户通过用户名(或从关注者下拉列表中选择)搜索用户，以便将他们添加到混合中。作为一个数据驱动的项目经理，我们希望将此作为一个实验来验证我们的假设和假说。

问题

由于多个用户之间的交互对该功能至关重要，因此一个用户使用该功能的方式并不独立于另一个用户使用该功能的方式。因此，如果我们将此作为传统的 AB 测试来运行，将每个用户作为随机单元，那么我们将违反 SUTVA。从用户的角度来看，如果用户可以在 Spotify 中搜索并邀请他们的朋友，但他们的朋友在控制组中，看不到邀请，这将中断他们的体验，并影响他们使用该功能的方式。

解决这个问题的一个方法是将网络分成互相跟随的用户群，并在群集级别执行随机分配。然而，由于这是 Spotify 的首批多用户体验之一，让我们假设没有任何基础设施来运行这种类型的实验。就开发资源而言，建立基础设施的成本很高，我们希望更有资源来更快地迭代。

让我们假设，通过一些探索性的数据分析，我们发现 98%的用户只创建与来自他们自己地理位置的其他用户的混合。通过将我们的特征变化限制在特定的地理区域，我们非常接近于重建 SUTVA。因此，这似乎是一个地球同步轨道提升实验的绝佳候选。我们承认，地理提升实验的结果并不像我们要进行的社交聚类实验或如果我们 100%的用户只与来自他们自己地理的用户创建混合那样是一个强有力的因果关系指标。但是我们认为，如果把这个实验作为一个地理提升实验来运行，我们会得到这个变化的影响的一个很好的指示。

我们的控制地理位置将为geo_1，我们的治疗地理位置将为geo_2。我们将只在geo_2发布改变，我们相信这将增加点击邀请按钮的用户被至少一个朋友接受邀请的可能性。因此，我们的成功度量将是表示用户是否成功邀请 1 个朋友加入混合的转换度量。

混杂因素

让我们假设我们已经仔细地选择了geo_1和geo_2，这样我们可以控制的任何混杂因素都可以在两个地理区域一致地应用。任何可能对此功能的指标产生影响的普遍获得活动都适用于两个地区。

现在我们已经清楚了上下文，让我们假设我们已经运行了这个实验一段时间，并且已经收集了足够的数据。所以让我们深入分析一下。

分析

在开始分析之前，值得注意的是，这里使用的数据是人工生成的，因此已经是本练习的理想格式。在现实世界中，我们可能需要执行一些 ETL 操作来获得这种格式的数据。然而，这超出了本文的范围。这些数据也不代表 Spotify 的实际用户参与度数据。

让我们导入依赖项并研究实验数据。

import pandas as pd
import numpy as np
import pymc3 as pm
import theano.tensor as tt
import matplotlib.pyplot as plt
import arviz as azexperiment_data = pd.read_csv('experiment_data.csv')experiment_data.head()

我们看到数据包含一个散列来表示用户 Id，用户在哪个组中，转换后的列表示用户的邀请是否成功。让我们继续总结不同组的转化率。

我们将把每组的转换率建模为一个具有 Beta 分布的变量，就像我们之前在做转换率的贝叶斯建模时所做的那样。我们不会通过检查来估计每个 Beta 分布的参数，而是将它们建模为 0 和相应组样本大小之间的均匀分布变量。所以让我们开始这样做:

diff_model = pm.Model()with diff_model:
    geo_1_pre_a = pm.Uniform('geo_1_pre_a', lower=0, upper=1367)
    geo_1_pre_b = pm.Uniform('geo_1_pre_b', lower=0, upper=1367)    

    geo_1_post_a = pm.Uniform('geo_1_post_a', lower=0, upper=1893)
    geo_1_post_b = pm.Uniform('geo_1_post_b', lower=0, upper=1893)

    geo_2_pre_a = pm.Uniform('geo_2_pre_a', lower=0, upper=1522)
    geo_2_pre_b = pm.Uniform('geo_2_pre_b', lower=0, upper=1522)    

    geo_2_post_a = pm.Uniform('geo_2_post_a', lower=0, upper=1408)
    geo_2_post_b = pm.Uniform('geo_2_post_b', lower=0, upper=1408)

既然我们已经对分布参数进行了建模，我们就可以对每个转换率的分布进行建模了。

with diff_model:
    geo_1_pre_cr = pm.Beta('geo_1_pre_cr', alpha=geo_1_pre_a, beta=geo_1_pre_b)

    geo_1_post_cr = pm.Beta('geo_1_post_cr', alpha=geo_1_post_a, beta=geo_1_post_b)    

    geo_2_pre_cr = pm.Beta('geo_2_pre_cr', alpha=geo_2_pre_a, beta=geo_2_pre_b)    

    geo_2_post_cr = pm.Beta('geo_2_post_cr', alpha=geo_2_post_a, beta=geo_2_post_b)

在对每个转换率进行建模后，我们现在可以对特性发布前后的差异进行建模，同时对差异中的差异进行建模，这就是我们的提升。

with diff_model:
    diff_pre = pm.Deterministic('diff_pre', (geo_2_pre_cr - geo_1_pre_cr))
    diff_post = pm.Deterministic('diff_post', (geo_2_post_cr - geo_1_post_cr))
    diff_in_diff = pm.Deterministic('diff_in_diff', diff_post - diff_pre)

既然我们已经对差异进行了分层建模，我们可以添加我们的观察值。但在此之前，我们不妨也对每个 geo 中的提升进行建模，以帮助我们更好地理解不同混杂因素导致的提升。

with diff_model:
    diff_geo_1 = pm.Deterministic('diff_geo_1', (geo_1_post_cr - geo_1_pre_cr))
    diff_geo_2 = pm.Deterministic('diff_geo_2', (geo_2_post_cr - geo_2_pre_cr))

最后，我们将观察到的转化率添加到模型中，然后对其进行采样。我们将这些转化率建模为伯努利变量，该变量使用之前建模的转化率作为转化率的概率。

with diff_model:
    geo_1_pre_conversions = pm.Bernoulli('geo_1_pre_conversions', p=geo_1_pre_cr, observed=conversion_values['geo_1_pre'])

    geo_1_post_conversions = pm.Bernoulli('geo_1_post_conversions', p=geo_1_post_cr, observed=conversion_values['geo_1_post'])    

    geo_2_pre_conversions = pm.Bernoulli('geo_2_pre_conversions', p=geo_2_pre_cr, observed=conversion_values['geo_2_pre'])    

    geo_2_post_conversions = pm.Bernoulli('geo_2_post_conversions', p=geo_2_post_cr, observed=conversion_values['geo_2_post'])trace = pm.sample()

一旦模型被采样，我们就可以绘制采样分布图。

with diff_model:
    az.plot_trace(trace, compact=False)

查看diff_in_diff图，隆起很可能大于 0，隆起的最高概率约为 5%。

我们现在可以检查所有模型的统计数据，并了解它们的 95%可信区间在哪里。

with diff_model:
    display(az.summary(trace, kind='stats', round_to=2, hdi_prob=0.95))

特别是，我们想看看最近 3 个模型的统计数据— diff_in_diff、diff_geo_1和diff_geo_2。

看起来我们的diff_in_diff平均值为 5 %, 95%的可信区间在 0-9%之间。因此，虽然这一变化很可能确实影响了用户邀请的成功程度，但也很可能只是很小的影响。

查看我们的控制地理，diff_geo_1表明在geo_2发布变更前后平均有 5%的差异，95%可信区间位于 2–9%之间。由于这个原因，我们没有发布特性变化，这种变化很可能是由于混杂因素造成的。

结论

虽然看起来我们想要测试的变化可能对用户的邀请有多成功产生了影响，但这种影响可能非常小。如果我们对此次升级感到满意，我们可以将该功能推广到多个地区，并持续监控成功指标。然而，如果我们觉得影响太小，我们可以重新评估我们的假设，想出其他方法来改善用户体验。

感谢阅读这篇文章！我希望它能帮助你更好地理解当传统的实验不起作用时，如何用这种方法来测试真实的产品特性。

如果你喜欢阅读我的文章，愿意支持我的写作，并且正在考虑订阅一个媒体，请随时使用我下面的推荐链接。我会从你的订阅费中提成。

https://medium/@kaushsk12/membership

自行车共享系统日常运行的地理可视化:赫尔辛基和塔尔图

原文：https://towardsdatascience/geo-visualization-of-daily-movements-of-bike-sharing-system-helsinki-and-tartu-c38b4aed6b06

与伊莱亚斯·威尔伯格合作

图片由 Elias Wilberg 提供。从 2021 年的一个夏季工作日开始，赫尔辛基和埃斯波的自行车共享自行车

最近，我在 GIS4 Wildlife 开始了与大规模运动数据集地理可视化相关的项目。我刚刚想起，在 30 日地图挑战赛期间，我的一位同事 Elias Wilberg 创作了一幅精彩的可视化作品，展示了赫尔辛基自行车共享系统的日常活动。可视化渲染的环境是 PostgreSQL/PostGIS+Qgis+time manager 之间的连接。因此，我决定向 Ellias 询问渲染过程的体验，并将他的回答包含在这篇文章中。

Elias 专注于研究城市和主动移动性，尤其是自行车共享系统数据的使用。我在一篇关于赫尔辛基地区自行车共享运动聚合的旧文章中提到过他的研究。

但是我发现 Elias 在 30DayMapChallenge 中的可视化很吸引人，它代表了自行车共享系统用户在他们的轨迹中从起点到终点的移动。我记得当他在处理它的时候，由于数据集的大小，他在晚上离开了计算机渲染。运行脚本 12 个小时，希望在处理过程中没有错误，这种体验对于大型数据集来说是常见的。因此，在开始我未来与大规模运动数据集可视化相关的项目之前，Elias 的经历对我来说非常有价值，很高兴与地理空间社区分享它。

总的来说，我询问了地理可视化产品的使用体验，Elias 回答道:

“我对自行车共享数据产生了兴趣，因为这是为数不多的关于自行车的详尽且开放可用的数据源之一。通过帮助揭示城市地区人们的时空骑行模式，自行车共享数据可以帮助理解如何让更多的人骑自行车；尤其是由于气候危机，这是城市迫切需要的东西。赫尔辛基地区的当地自行车共享系统已经存在多年，在市民中非常受欢迎，每年约有 350 万人次。这意味着该系统每年都提供易于获取和使用的有趣数据

在我看来，成功分析和可视化数据的关键是使用正确的工具。当你有合适的工具时，即使使用非常简单的硬件，也有可能从大型数据集进行动态可视化和高级分析，就我而言，使用的是一台 5 年前的笔记本电脑。

以自行车共享可视化为例，我使用 Python 进行路由，使用 PostGIS 进行几何运算，使用 QGIS 进行可视化，使用 GIMP 将地图框组合成视频格式。从路线开始，路线规划工具的可用性和效率在最近几年有了很大的提高。有了 Pandana 图书馆，赫尔辛基所有自行车共享站(大约有 450 个)之间的最短路线搜索只需不到 5 分钟。虽然通常是空间操作需要最大的处理能力，但 PostGIS 是一个很好的工具。我能够将最短的路线分成超过 1000 万个节点，并使用 PostGIS 为每个节点插入时间戳，同样在不到 5 分钟的处理时间内。虽然 QGIS 本身可能会塞满那么多点，但根据时间戳将点以更小的块从 PostGIS 输入到时态控制器使可视化成为可能。在这一阶段，QGIS 在 24 小时内每秒生成一个地图框，但正如 Bryan 所提到的，这需要更多的时间。最后，使用 GIMP，可以将帧组合成动态可视化，并优化视频文件的大小。最棒的是所有的事情都是用开源工具完成的。开放万岁！

最后，非常感谢布莱恩。如果没有你的兴趣和热情，我就不会花时间制作这个可视化的东西。很高兴您能和我们一起来到数字地理实验室！"

灵感和伟大的想法不仅仅来自埃利亚斯。由蒂格兰T2【Khachatryan撰写的关于赫尔辛基地区自行车网络的网络分析非常有价值。中型简介包含非常有趣的文章，有许多见解。如果你在 Medium 上，查看他的个人资料，然后关注蒂格兰·哈恰特良。尤其是他的上一篇文章。

https://geometrein.medium/wolt-delivery-network-analysis-cccdc4cb50e3

关于自行车共享系统的许多灵感及其数据集分析的许多可能性。所以，我决定加入我自己的。我可以更像 Elias 那样描述我的可视化，因为我想创建一个地图动画，重点关注自行车共享系统用户在他们的轨迹中的运动。查看我这篇关于爱沙尼亚塔尔图自行车数据地理可视化的文章。

地理可视化

图片由作者提供。自行车街道网络上的自行车运动

塔尔图智能自行车系统由 750 辆自行车和 80 多个站点组成。该网站包含有关自行车使用情况的更新指标，例如当年行驶的距离等。此外，考虑到三分之二的自行车都配备了电动辅助马达，以在踩踏板时提供额外的动力，它为自行车的使用提供了可接受的价格。

我将从 OSM 数据中提取的自行车街道网络添加到地理可视化中。获取数据集的代码可以在本文的知识库中找到:

关于数据集和许可

爱沙尼亚提供了一个开放数据门户，在那里你可以找到自行车共享系统的数据。对于这个地理可视化，我使用了 2019 年 7 月 18 日的每日数据集。

• 自行车共享系统数据:知识共享署名许可-共享 3.0 未发布。
• OSM 自行车网络:获得开放数据共享开放数据库许可(ODbl) 或归属许可。用户可以自由地复制、分发、传输和改编这些数据，只要这些数据是像 OpenStreetMap 贡献者这样的作者所拥有的。

代码和 Web 地图

按下播放并检查移动。建议在“火箭”按钮上减速。

在这里找到代码:储存库
可视化在这里: Tartu 智能自行车

图片由作者提供。塔尔图智能自行车路线

一旦运行了代码并保存了文件。
将它们拖至 KeplerGl 演示。

一条重要的建议是，您可能需要准备时间戳:

data['route_code'] **=** [f'B-{code}' **for** code **in** data['route_code']]
data['timestamp'] **=** data['coord_date'] **+** ' ' **+** data['coord_time']

检查整个笔记本，它有完整的工作流程和添加 OSM 的自行车路线层的过程。

结论

本文旨在强调大规模运动数据集需要支持可视化的环境。目前最好的组合是 PostgreSQL/PostGIS+Qgis+time manager。自行车共享系统数据集是运动数据集的很好的例子，Elias 的经历对未来的工作是有价值的建议。如果你知道更好的环境或集成(当然有)，请让我知道或联系我合作。祝愿你未来的可视化。

带四键的地理围栏

原文：https://towardsdatascience/geofencing-with-quadkeys-7c5b9866ff98

本文介绍了如何创建带有方形分区的地理围栏

上图显示了使用此处描述的算法对葡萄牙大陆进行的等级方形离散化。图片由作者用叶子包和 OpenStreetMap 图像制作。

地理围栏是划定任意地理区域的虚拟边界。我们使用这些构造来分类某个对象是否位于感兴趣的区域内，例如当我们想要确定一辆汽车是否进入了特定的停车场时。地理围栏不必与特定的物理位置相对应，可以传达更微妙的概念，如龙卷风走廊或禁飞区。

为了让地理围栏发挥作用，我们需要一些地理定位设备来报告感兴趣的对象在世界上的位置，并需要一种方法来定义围栏本身。我们通常使用地理空间多边形来处理这最后一点，其顶点是在纬度和经度空间中定义的位置。简单地说，当物体的跟踪设备报告一个我们可以用数学方法确定的位于多边形内部的位置时，物体就在围栏内部。在某种意义上，geofence 是一个二元分类器，用于确定对象相对于 geofence 内部的位置。

地理围栏系统传统上使用特定的存储和查询工具，这些工具允许从地理多边形到透明地映射地理空间顶点的内部表示的直接转换。这些系统是精确的，但是需要定制的存储和索引格式，这些额外的需求通常会降低存储空间和查询性能。例如，想想 PostgreSQL 扩展 PostGIS ，它扩展了基本功能以添加地理空间处理特性，比如函数、数据类型和索引。

通常情况下，应用不需要高精度的地理围栏检测。想一想需要什么来确定你的一辆车是否进入了加油站。没有必要编码位置边界的所有细节，一个粗略的草图可能就是你所需要的，特别是如果它将使位置的编码更简单和检测更快。此外，如果您可以使用更轻便的地理围栏系统，您可能希望在您的 edge 物联网设备上实现它。

哈希备选方案

那么我们有什么选择呢？有趣的是，一些选项使用固定的分层地理空间网格，六边形或正方形。这些替代方法可以将 geofence 区域缩小到 64 位整数索引的列表，您可以轻松地存储和快速查询这些索引。您只需要一个常规的数据库或内存存储以及基本的集合包含逻辑就可以让这些替代方案发挥作用。

这些算法根据它们选择的几何形状将纬度和经度空间中的位置散列成唯一的整数。两者都是通过建立一个固定的层次地理空间网格来工作的，所以你必须接受一些近似误差。幸运的是，您可以使用这两种方法来离散大约 5 米或 16 英尺。

优步·H3 的选择

在其他文章中，我探索了优步的 H3(T1)的潜力，这是一种分层的六边形空间索引。这个空间索引系统是强大而轻便的，因为它是用 C 语言实现的，允许在边缘设备上使用。

由于其哈希特性，我们还可以使用 H3 来实现高速地理围栏查询。正如我在另一篇文章中所展示的，我们可以将地理围栏设想为一组 H3 指数，每个指数都有一个相关的空间六边形。下图显示了这个概念，其中我们用一组相连的 H3 指数来定义地理围栏。包含测试只是将位置转换为相关的 H3 指数，然后运行集合包含测试。

上图显示了如何从一组 H3 六边形设计地理围栏。地理围栏包含所有六边形和所有红点。为了测试任意点的包含性，我们需要将其转换为适当的 H3 指数，并执行集合包含测试。图片由作者用叶子包和 OpenStreetMap 图片制作。

请注意，存储这样的地理围栏是非常便宜的，每个六边形作为一个 64 位整数进行加权。上面描述的地理围栏只有 45 个这样的整数。但是当地理围栏变大时会发生什么呢？我们可以不断将六边形添加到不断增长的 geofence 中，但这种方法有一定的局限性，因为 H3 表示可能会比其他 geo-polygon 表示更加重要。

我们可能会选择使用更大的六边形尺寸来替代日益严重的地理围栏问题。正如你可以从优步的文档中看到的，H3 六边形等级意味着当增长到最接近的级别时，每个六边形代表较低级别的七个六边形。但是这种选择是以丢失细节为代价的。

然而，优步·H3 API 确实提供了一个解决方案，允许我们在层次结构中混合和匹配不同级别的六边形。下面两张葡萄牙大陆的图片展示了这一过程是如何进行的。

上图显示了离散为 2，769 个 6 级 H3 六边形的葡萄牙大陆地图。图片由作者用叶子包和 OpenStreetMap 图像制作。

上图显示了将葡萄牙大陆地图离散为 2，769 个六级六边形。优步提供了将这些六边形合并到更高级别的六边形以填充整个地理围栏形状的选项。我们可以在下图中看到这个过程的结果。

通过压缩前一组 H3 六边形，我们将计数减少到 501。图片由作者用叶子包和 OpenStreetMap 图像制作。

这种方法的显著优点是，它将六边形的数量减少到 501，或整整一个数量级。缺点是，由于其固有的几何形状，六边形不会完全重叠，并且会有一些小区域未被覆盖。正如优步的文献所解释的那样，这种效应对于许多统计应用来说并不重要。但是当涉及到地理围栏时，上面的简化杀死了它。

下图更详细地说明了这种效果。

上图显示了葡萄牙大陆中央区域的 H3 压缩离散化。不同层次的六边形之间的差距非常明显。图片由作者用叶子包和 OpenStreetMap 图片制作。

这是地理围栏的相关问题吗？正如我在上面和之前的一篇文章中所展示的，如果你能忍受这种参差不齐的边缘，这种方法会非常有效。geofence 仅包含无缝连接的相同大小的六边形，因此不会出现包含检测错误。另一方面，对于较大的形状，压缩后的案例显示的间隙可能大到足以导致 geofence 包含检测错误。

四键选择

四键的几何结构确实解决了间隙问题。正如你从开始的图像中看到的，正如四键数学所暗示的，正方形的层次在不同层次的四键之间提供了紧密的配合。你可以从它的来源、微软在线文档和下面的文章中读到更多关于这个主题的内容。

如下图所示，与 H3 不同，四键层次结构允许你紧密地适应不同层次的方块，没有重叠或死区。请注意，方形图块的固定位置意味着我们只能在层次结构中从较低到较高的级别(较低的缩放比例)合并特定的图块集。

上图展示了 Bing 地图切片系统的层次结构。图片由作者用叶子包和 OpenStreetMap 图片制作。

上图揭示了四键的两个优秀特性。首先，四键代码是以四进制表示的整数，这意味着每个“缩放”级别需要两个额外的位来表示。第二，每个四键对应于一个更大的地图上的一个图块，我们可以把它想象成一个位图，每个图块作为一个“像素”我们将在这里使用这两个特性。

数字四键编码

让我们注意，它需要一个数字(0、1、2 或 3)来表示缩放级别。我们在将四键描述为字符串时使用这种方法，因为它方便地打包了图块位置信息和缩放级别(字符串长度)。

通过认识到每个数字只需要两位来编码，我们也可以将这种表示转换成二进制整数。如果我们保持在线地图通常可用的最大 23 个缩放级别，我们需要 48 位来存储完整的四键哈希。有足够的空间来保存缩放级别，只需要五位。这是一个可能的编码的公开提议:

https://github/joekarl/binary-quadkey

就个人而言，我会颠倒编码设计，在低位保留散列，在高位保留缩放级别。这种方法会使下面需要的一些数学运算变得更容易。

为了简单起见，我们将使用基本编码，仅使用低 48 位来存储散列，并在其他地方保持缩放级别。这种方法使得代码更容易实现和遵循。

隐式位图

如前所述，我们可以认为四键仅仅是地理空间散列，地图上的方块，或者是位图上的像素。下图显示了 EaEarth 的 sap 是一组 8x8 的瓦片。为什么不是 8x8 的位图，我们可以在上面画线和多边形？

上图显示了缩放级别为三的 Bing 地图切片。请注意，我们可以将这些看作数组或位图。图片由作者用叶子包和 OpenStreetMap 图像制作。

位图概念在高缩放级别时变得非常有用，在这种情况下，每个图块都像我们喜欢的那样小，我们可以使用它们来“绘制”对应于地理围栏区域的填充多边形。当然，这幅画是一个虚拟的概念。通过“绘制”一个像素，代码将向给定缩放级别的现有集合添加一个四键代码。

将附近的瓷砖合并成更大的瓷砖是一个简单的过程。请参考前面说明缩放操作的图像。请注意，当绘制索引值以数字 3 结束的四键(最后两位设置)时，我们可以检查相邻三个图块的存在。当绘制图块“213”时，我们可以检查是否存在“212”、“211”和“210”如果这些图块存在，代码可以用图块“21”替换它们，并将其放置在适当的缩放级别设置上。如果适用，重复该过程。

现在我们可以看看在位图上绘制填充多边形的过程。

在位图上绘制填充多边形

我们首先为多边形顶点选择一个合适的表示。您可以使用可用的 GitHub 库来跟踪实现。我们从通用多边形填充脚本的代码开始。

上面的类存储多边形顶点的每个实例。(图片来源:作者)

该算法的第一步包括从多边形顶点到边的转换。下面的代码显示了如何存储多边形的边。

多边形边需要特定的顺序，填充算法才能发挥最佳效果。(图片来源:作者)

暂时不要担心“活动边缘”。我们将在后面讨论这个概念。

对于每条边，我们收集 y 坐标的最小值和最大值、与最小值 y 相关的 x值以及边斜率的倒数。存储多边形边信息的类对填充算法所需的特殊排序方法进行编码。此功能与选择接下来要加工的边相关。

为了创建多边形边，我们使用一个专用函数来计算所需的属性，并将它们存储在目标存储对象上。所有这些值将帮助多边形填充算法决定下一条水平线的起点和终点。

上面的代码使用连续的顶点对创建边。(图片来源:作者)

每对连续的顶点都转换成一条多边形边，如以下代码所示:

上面的函数将顶点序列转换成相应的边序列。注意代码是如何复制第一个点来闭合多边形的。(图片来源:作者)

请注意，该函数会过滤掉任何水平边缘(参见第 9 行)。该测试明确地避免了生成边的代码被零除的可能性，并反映了多边形填充算法是如何工作的:通过顺序填充水平线。增加这样的边是多余的。

我们现在可以看看绘制填充多边形的代码。这里我们使用一个 NumPy 数组作为位图。我们首先确定多边形的边，创建支持数组和活动边列表。这个列表是算法工作的基础。

上面的代码包含了位图的多边形填充算法，这里实现为 NumPy 数组。(图片来源:作者)

在第 12 行，代码进入一个循环，遍历所有扫描行，并为每一行设置一个活动边列表。活动边缘与当前绘制的扫描线相交。我们使用一个不同的数据结构来保存这些，因为算法强加了可变性(见第 30 行)。

上面的类定义了活动边的属性。(图片来源:作者)

该算法使用两个函数来管理活动边缘的列表，一个函数插入其最小值 y 等于当前扫描线的边缘，另一个函数移除不再与其相交的活动边缘。

请注意，第一个函数在插入时将边转换为活动边。它对边列表进行操作，并返回活动边列表。第二个函数仅适用于活动边。(图片来源:作者)

在更新和排序活动多边形边的列表后，算法可以绘制所需的水平线。内部循环处理成对的活动多边形边，并在它们之间绘制一条水平线。在外循环的最后，代码更新所有有效边沿的 x 位置。

下面是调用多边形填充函数的方法:

上面的代码定义了一个任意多边形，将其转换成 NumPy 数组，最后显示出来。(图片来源:作者)

这是结果:

上图显示了运行前一幅图像的代码的结果。(图片来源:作者)

现在，我们已经了解了如何在位图上绘制填充多边形，让我们将该过程扩展到使用不同缩放级别的四键生成 geofence。

用四键绘制填充多边形

可以想象，将地理多边形转换为一组四键类似于上面的算法。这个过程在概念上非常简单。我们从要编码的地理多边形开始，使用提供的比例表确定目标细节级别，将所有地理多边形顶点转换为图块坐标，然后绘制填充的多边形。在此过程中，我们递归地将图块聚集到较高的细节级别，从而减少最终的图块数量。请遵循地理多边形四键填充脚本上的代码。这个脚本类似于前面的脚本，但是它不是显示位图，而是将生成的四键保存到数据库中。关于这个问题的更多信息，请参见库文档。

我们首先将一个多边形定义为一组有序的空间坐标，并对源多边形使用纬度和经度空间坐标。这个空间是连续的，纬度从-90 到 90 度，经度从-180 到 180 度。

下一步需要将这些坐标离散到一个等距空间，其粒度取决于我们想要达到的细节水平。根据经验，我使用微软提供的地面分辨率表来粗略估计每个“像素”的大小。选择名为“*地面分辨率(米/像素)*的第三列的值，并将其值乘以 256，这是“滑动地图”图块的标准大小。得到的值是赤道上“像素”宽度的合理估计，单位为米。在 20 级和零纬度，每块瓷砖的尺寸为 38 米(约 125 英尺)。随着远离这些纬度(向北或向南)，像素大小将会减小。

该算法首先将纬度和经度上的所有源数据点转换到目标离散空间。每个坐标现在是一个整数对，而不是实数对。如前所述，离散化取决于“缩放级别”

坐标之间的转换是一个两步过程。首先，我们将纬度和经度转换为所需级别的四键。接下来，我们提取四键瓦片坐标。

上面的函数实现了从纬度和经度到指定细节级别的图块坐标的两步转换过程。(图片来源:作者)

在按顺序通过所有地理多边形顶点运行此转换后，我们就可以使用栅格多边形填充算法了，该算法与上一节中的算法非常相似。它通过从最小(顶部)到最大(底部)循环通过所有离散扫描线( y 坐标)来工作。该函数绘制的“画布”是一个字典，它将一个细节级别映射到同一级别的整数编码四键的集合。让我们通过下面的代码来遵循算法逻辑。

上面的代码显示了四键多多边形填充功能。主循环从上到下遍历扫描线，同时管理活动边。(图片来源:作者)

该函数首先使用下面的代码从多边形中收集边的列表。注意除了边缘列表，该函数如何返回扫描线 y 坐标范围。

(图片来源:作者)

接下来，我们进入主循环，遍历离散化的 y 坐标范围。除了像素绘制(第 20 到 24 行)之外，代码与我们之前看到的类似。如你所见，我们使用了两种不同的绘图功能，一种用于偶数扫描线，另一种用于奇数扫描线。原因是偶数编号的扫描线没有提供机会来合并附近的图块以生成较低分辨率的图块(请参见上面的讨论)。

上面的代码显示了偶数扫描线生成函数。请注意，它只是用四键范围更新缩放级别设置。(图片来源:作者)

下面的函数将图块坐标转换为四键索引，如上所述。

上面的代码显示了如何在给定的缩放级别将图块坐标转换为四键索引。注意，我们生成的位序列没有添加电平数据。(图片来源:作者)

奇数扫描线绘图功能更有趣，因为它必须处理合并图块的可能性。

上图包含了在奇数扫描线上绘图的函数。这些扫描线是图块合并的候选。(图片来源:作者)

它不仅仅更新缩放级别集，而是遍历生成的四键并调用一个专门的函数。对于每个要插入的四键，该函数测试一个终止数字“3”，这是潜在扩展的签名。如果是，它还会测试相邻四键的缩放级别设置。请注意简化的四键代码二进制编码是如何有用的:

上面的代码执行隐式递归图块扩展和插入。(图片来源:作者)

扩展四键包括移除相邻键，并在上面的缩放级别中插入一个更大的键。代码遍历缩放级别，直到无法再找到要展开的四键。

结果是一个 Python 字典，其中缩放级别作为键，缩放级别集作为值。下面的函数显示了如何将新生成的围栏插入数据库。

上面的代码在数据库中以相同的名称插入了一个区域列表。这种方法对于具有多个多边形的国家边界非常有用。(图片来源:作者)

询问

我们现在可以开始查询生成的数据，在这种情况下，通过填充的数据库。这个想法很简单:给定纬度和经度空间中的任意位置，我们想要确定它属于哪个地理围栏(如果有的话)。但是首先要克服一个障碍。我们已经生成了各种缩放级别的四键，那么我们如何知道对任意位置使用哪一个呢？

简单的回答是我们不知道，所以我们查询所有的缩放级别。这个解决方案并不像看起来那么复杂，而且，如果有合适的索引，它也很快。让我们看看如何进行查询。

我们的第一步是确定我们的数据库包含哪些缩放级别。这里，我们在实例化数据库对象时这样做。

数据库查询示例代码。(图片来源:作者)

查询过程从使用存储的最大缩放级别将纬度和经度坐标转换为四键开始。接下来，代码将四键转换为编码整数，包含值和缩放级别。我们去掉最后一个，然后把散列转移到适当的位置。接下来，对于所有当前缩放级别，代码生成一个相应四键值的列表。我们确信所有这些四键都包含输入位置。最终的数据库查询非常简单。只需使用一个集合包含子句。

请注意，为了让这个查询在大规模运行时表现良好，我们需要向表中添加一个适当的索引。对于大规模数据集，我们甚至可以考虑发出并行查询，通过四键前缀分割表以限制索引大小，等等。

结论

地理围栏是地理空间分析行业的工具之一。我们已经看到了它们的常规实现，并讨论了它们的缺点。我们还提出了基于散列或网格的替代方案:基于六边形的优步 H3 和基于正方形的四键。我们看到了 H3 如何合并不同缩放级别的六边形，以及该解决方案如何将一些区域排除在地理围栏之外，从而导致潜在的检测错误。最后，我们看到了基于四键分层正方形的替代方案如何通过生成无缝瓦片集来工作，以及如何高效地查询这些瓦片集。

我希望你喜欢这篇文章！

参考

Git 储存库

https://github/joaofig/quadkey-geofence

地理空间哈希

H3:优步城市等级空间索引

必应地图磁贴系统—必应地图|微软文档

🌐微软(github)提出的使用四键的地理平铺的 Python 实现

多边形填充

多边形填充教学工具

高效扫描线多边形填充算法(PDF)的实现

扫描线填充算法(PDF)

扫描线填充算法—加州大学戴维斯分校(PDF)

国界

https://github/datasets/geo-countries

国家多边形作为 GeoJSON —数据集—数据中心—无摩擦数据

世界行政边界—国家和领土— Opendatasoft

joo Paulo Figueira 在葡萄牙里斯本的TB . LX by Daimler Trucks and bus担任数据科学家。

球形数据的几何深度学习

原文：https://towardsdatascience/geometric-deep-learning-for-spherical-data-55612742d05f

球形 CNN

通过对物理世界平移对称性的理解进行编码，卷积神经网络(CNN)彻底改变了计算机视觉。在这篇博文中，我们研究了 CNN 成功背后的原理是如何转移到一系列问题上的，这些问题的数据表现出复杂的几何形状，例如球体。

这篇博文由来自 Kagenova 的 Oliver Cobb 和 Augustine Mavor-Parker 共同撰写。

球形数据的一个例子。[图片由 NASA 在 Unsplash 拍摄]

在过去的十年中，CNN 从传统(平面)图像和视频中提取语义的能力已经得到了快速的提高。如果有足够的数据，通常可以达到人的水平。然而，对具有空间结构的数据进行分析远不是一个已解决的问题。有许多问题的数据表现出空间的但非平面的结构。示例包括虚拟现实中的 360°图像、大爆炸的宇宙微波背景(CMB)辐射、医学成像中的 3D 扫描以及计算机图形中的网格表面，仅举几例。

对于这些问题中的每一个，我们都希望利用我们对数据结构的了解，特别是它们所涉及的对称变换。正如在 之前的一篇博文 中所讨论的，将对对称性的理解编码到机器学习模型中可以是一种限制所考虑的模型空间的强大方法，允许更有效地学习模型。

对于平面上的图像，平移对称可以通过应用在图像上平移的一堆卷积滤波器来容易且有效地编码。因为相同的卷积滤波器被应用于所有位置，所以结果操作在平移上是等变的，即它尊重平移对称性。这意味着无论特征位于图像中的哪个位置，它都会以相同的方式刺激相应位置的激活神经元。

不幸的是，对于非平面结构的问题，通常不存在这种简单的程序来编码对对称性的理解。然而，对于这些问题，几何深度学习新兴领域的研究人员正在制定新的方法，这些方法利用数据的几何形式的属性，并尊重对称性。最近取得重大进展的一组问题是那些在球面上定义数据的问题。

球形数据的对称性

许多字段涉及到固有地存在于球体上的数据。

当在球面上的每个点进行观测时，会产生球面数据，例如地球的地形图。然而，当在多个方向上进行观测时，也会出现这种情况，例如宇宙学中的宇宙微波背景(CMB)或虚拟现实和计算机视觉中的 360°图像(见下图)。在 Kagenova，我们正在努力为这些问题和其他涉及复杂几何数据的问题(如球体)解锁深度学习的非凡成功。

球形数据的例子。【原创人物由作者创作。]

对于平面图像，CNN 规定，定义特定特征如何变换的规则不应取决于该特征恰好位于平面中的位置。对于在球体上定义的数据，我们希望规定规则不应该依赖于特征如何以及在哪里碰巧在球体上定向。变换特征然后旋转其变换后的形式应该等同于旋转特征并变换其旋转后的形式。关于此属性的操作称为旋转等变(见下图)。

旋转等方差图。给定球形数据(左上)，应用变换(𝒜)以获得特征图(右上)，然后旋转(ℛᵨ)特征图(右下)，相当于首先旋转数据(左下)，然后应用变换(右下)。【原创人物由作者创作。 ]

在物理学中，规定控制系统行为的物理定律不应该取决于系统的取向，这就产生了角动量守恒定律。因此，毫不奇怪，在量子物理中用于研究角动量的一些相同的机器对于在深度学习中定义旋转等变层是有用的(正如我们将在后面看到的)。

标准(平面)CNN 的局限性

在深入球形深度学习的想法之前，也许很自然地会想为什么平面 CNN 的有效性不能被直接利用。我们能不能不把我们的球形数据投射到某个平面表示上，然后简单地用通常的方式应用 CNN？毕竟，我们对球形世界的平面投影(地图)很熟悉。

投影的问题是不存在从球面到平面的既保持形状又保持面积的投影。换句话说，扭曲是不可避免的。

这就是为什么格陵兰岛在地球地图上看起来和非洲差不多大，而实际上还不到非洲的十分之一(见下图)。

球体到平面的投影引入了不可避免的扭曲，这与所用的投影方法是不可避免的。因此，在地球地图上，格陵兰岛看起来和非洲差不多大，但实际上还不到非洲的十分之一。[图片来源于维基共享资源。]

这些失真意味着，当将传统的 CNN 应用于球形图像的平面投影时，特征根据它们所处的位置而出现不同。当应用于球面图像的平面投影时，平面 CNN 的平移等方差不编码旋转等方差。对旋转等变进行编码需要卷积的概念，这是专门为球体的几何形状设计的。

卷积复杂性

不幸的是，由平面 CNN 实现的简单卷积过程不能应用于球形设置。

要了解为什么会出现这种情况，首先考虑平面数据的形式。平面数据被表示为像素值的 2D 阵列。对于定义在平面上的数据，我们可以在水平和垂直方向上均匀地间隔像素位置。平面的这种均匀采样意味着每个像素具有相关的邻居，并且所有像素在相同的相对位置(北、东北、东等)具有邻居。这意味着在相同样本位置定义的任何过滤器都可以通过平移以输入中的任何像素为中心，从而使样本精确对齐。

不幸的是，没有办法对球体进行采样，使得所有像素在相同的相对位置具有邻居。球体上的位置通常使用球面坐标来描述，θ测量极角，ϕ测量方位角。相对于θ和ϕ均匀间隔采样会产生下图左侧所示的球体采样。如果我们使用这些样本位置来定义过滤器，然后旋转过滤器，我们会发现样本位置并不对齐(见下图)。不管我们选择如何对球体进行采样，这都是正确的。

假设我们使用与球形数据相同的样本位置定义一个过滤器。由于样本不对齐，因此无法评估滤波器在各种旋转下与数据的匹配程度。这适用于球体的所有样本。【原创人物由作者创作。]

众所周知，不可能以旋转不变的方式离散化球体。因此，不可能构造出严格旋转等变的纯离散球面卷积。

为了构造一个能捕捉旋转等变性质的卷积的球面概念，我们必须考虑一个连续的表示。值得庆幸的是，存在这样一种表示，对于这种表示，可以执行卷积的自然概念。

考虑球面上连续信号的表示。这些是函数 f: 𝕊 → ℝ将一个值与球面𝕊上的每一点(θ，ϕ)相关联，而不仅仅是在选择的样本位置。正如圆上的连续信号(即周期函数)可以分解为正弦和余弦函数的加权和，球体上的连续信号同样可以分解为谐波基信号的加权和(见下图)。在这两种情况下，权重(系数)可用于表示信号，产生圆上信号的傅立叶级数表示和球上信号的球谐表示。

球面调和函数。[图片来源于维基共享资源。]

虽然这种表示是无限的，但是通过适当地截断系数向量，可以非常精确地近似真实世界的信号。从上图可以看出，低次球谐函数只能捕捉低频变化，而高次球谐函数可以捕捉较高频率的变化。我们截断的点决定了我们的数据表示的分辨率。

球形卷积

回想一下，我们想要执行满足旋转等方差属性的球形数据的转换。

有一个非常自然的球形卷积的概念，在连续设置中类似于在平面情况下执行。

这是取球面信号 f : 𝕊 → ℝ，定义一个球面滤波器 g : 𝕊 → ℝ，计算卷积信号 f * g 定义为

这里我们使用了旋转算子ℛᵨ ，定义为(ℛᵨg)(ω)=g(ρ⁻ω)。换句话说，它具有将相应的反向旋转应用于函数的定义域的效果，类似于我们如何考虑平移定义在 1D 线或 2D 平面上的函数。

对上述等式的解释是，卷积信号 f * g 捕捉到信号 f 在任何给定旋转ρ下与滤波器 g 的匹配程度(见下图)。这类似于平面情况，我们考虑滤波器在各种平移下与输入的匹配程度。

数据(左)与过滤器(中)的球形卷积的可视化，以产生特征图(右)。【原创人物由作者创作。]

这里的主要区别在于，定义卷积信号的空间，即旋转空间(3D 空间)，不同于定义被卷积的信号和滤波器的空间(2D 球)。在上面所示的说明性例子中，滤波器对于方位角旋转是不变的，因此输出保持在球体上。

然而，将输入信号从球体提升到 3D 旋转空间并不是特别有问题。随后可以在信号和定义在旋转空间上的滤波器之间执行卷积的类似概念。因此，给定球形输入，为了分层学习特征，我们可以执行一次球形卷积，从而在 3D 旋转组上产生激活图，然后执行我们想要的任意多的旋转组卷积。

要了解为什么上述卷积概念是旋转等变的，请注意，将输入旋转ρ相当于对积分内的滤波器应用额外的ρ⁻旋转。反过来，这具有旋转由ρ⁻定义的卷积信号的域的效果。换句话说，在执行卷积之前旋转输入相当于先执行卷积，然后旋转输出。

两个球形信号的卷积似乎需要对三维空间中的每个值计算二维积分。不过好在 f * g 和 f 和 g 的调和表示之间的关系很简单。通过在 f 和 g 的调和系数之间执行矩阵乘法，可以在调和空间中计算球形卷积。鉴于深度学习实践者非常习惯于利用 GPU 来高效地执行矩阵乘法，这尤其方便。

球形回旋是不够的

配备了可以有效实现的旋转等变线性运算，看起来我们已经拥有了重复应用该运算和分级学习特征所需的一切。

然而，有一个重要的组成部分，我们迄今为止忽略了提到——非线性的引入。

在平面网络中，非线性是由逐点激活函数引入的，即通过将选择的非线性函数分别应用于每个样本位置的值。由于平面采样方案的均匀性，这实际上是平移等变操作。然而，我们已经过渡到使用没有相关样本位置或值的谐波表示。虽然可以获得基于样本的表示，但是我们不能对球体进行均匀采样(如上所述),这意味着对每个样本应用相同的非线性函数不是严格的旋转等变操作。

然而，以这种方式引入非线性是可能的，并且如科恩等人(2018 年)和埃斯特韦斯等人(2018 年)所示，通常相当有效。然而，为了执行卷积和非线性运算，在谐波和基于样本的表示之间重复转换是麻烦的。此外，人们很自然会怀疑损失的等方差在多大程度上阻碍了性能。

在我们的下一篇文章中，我们将看到如何利用量子物理学的思想在调和空间中直接引入非线性，而不损害我们对旋转对称的尊重程度。

参考

[1]科恩，盖革，克勒，韦林，球形 CNN，ICLR (2018)， arxiv:1801.10130 。

[2]埃斯特韦斯，艾伦-布兰切特，马卡迪亚，达尼利迪斯，学习 SO(3)用球面 CNN 的等变表示，ECCV (2018)， arXiv:1711.06721 。

简单解释了几何分布

原文：https://towardsdatascience/geometric-distribution-simply-explained-9177c816794f

几何分布的简单描述和用途

莫里茨·金德勒在 Unsplash 上拍摄的照片

介绍

在这篇文章中我想讨论一个统计学中常见且容易理解的分布， 几何分布 。这种分布被用于许多行业，如金融、体育和商业。因此，如果您是一名数据科学家，了解这一点非常重要。

在本帖中，我们将详细介绍它的定义、直觉、一点数学知识，最后在一个例题中使用它。

定义、属性和类型

几何分布是一种离散的概率分布，它推断出我们在获得成功之前需要进行的 伯努利试验 的次数的概率。伯努利试验是指单个事件只有两种结果:以一定的固定概率成功或失败。

几何分布通常被称为离散版的 指数分布 。如果你想了解指数分布，我之前写过一篇关于它的短文，你可以在这里找到:

实际上有两种不同类型的几何分布:

• 首先是获得成功所需的试验次数**。**
• 二是成功前的失败次数。

第一种被称为移位几何分布。

实际上，这两种方法都可以使用，但需要从一开始就明确区别，以确保结果的一致性。

在本文中，我们将使用移位版本，因为我觉得它更容易用数学和直觉来处理。

几何分布的一个关键性质是它是https://en.wikipedia/wiki/Memorylessness无记忆的。这意味着，仅仅因为上一个结果是失败的，下一个就不太可能成功。这是因为每个伯努利轨迹都是独立的https://en.wikipedia/wiki/Independence_(probability_theory)。****

数学细节

由于几何分布与伯努利分布和二项式分布 密切相关，其 概率质量函数(PMF) 呈现类似的形式:

作者在 LaTeX 中生成的方程。

其中 p 是成功的概率， n 是获得成功所需的事件数。

【均值】 简单来说就是:

作者在 LaTeX 中生成的方程。

其中 p 是再次试验成功的概率。

让我们通过一个例子让这个理论更具体。

示例和图表

在常规六面模具上第五辊出一个 4 的概率有多大？

在这种情况下我们有:n = 5p = 1/6。** 把这个代入上面的 PMF，我们发现概率为:****

作者在 LaTeX 中生成的方程。

这种情况下的期望值是:

作者在 LaTeX 中生成的方程。

这意味着我们希望使用模具在第六个辊上轧制一个 4 。

我们还可以使用 Python 为一系列试验绘制这个场景，【n】***😗****

***# Import packages
from scipy.stats import geom
import matplotlib.pyplot as plt# Probability and Number of Trials
n = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
p = 1/6# Generate the PMF
dist = geom.pmf(X, p)# Plot the distribution
plt.figure(figsize=(12, 7))
plt.scatter(n, dist, linewidth=2, color='black')
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.ylabel('Probability', fontsize=20)
plt.xlabel('Number of Rolls', fontsize=20)
plt.vlines(n, 0, dist, colors='black', linewidth=2, alpha=0.7)
plt.savefig('plot.png')
plt.show()***

作者用 Python 生成的图。

我们观察到掷出 4 的概率随着掷骰数的增加而呈指数下降。这是有道理的，因为我们的第一个 4 不太可能发生在第 100 次滚动中。

这种概率相对于成功所需试验次数的指数下降是几何分布 PMF 的一般形式。

结论

在本文中，我们讨论、解释并绘制了几何分布图。要记住的关键点是，几何分布计算连续伯努利试验失败指定次数后的成功概率。

和我联系！

• 要在媒体上阅读无限的故事，请务必在此注册！T3💜
• 当我在这里发布注册邮件通知时，可以获得更新！ 😀
• 领英 👔
• 推特 🖊
• github🖥
• https://www.kaggle/egorphysics🏅

(所有表情符号由 OpenMoji 设计——开源表情符号和图标项目。许可证: CC BY-SA 4.0

几何先验 I

原文：https://towardsdatascience/geometric-priors-i-cc9dc748f08

几何深度学习

群、表示、不变性和等方差

一系列博文，总结了 AMMI 计划的 几何深度学习(GDL)课程； 非洲机器智能硕士 ，授课老师 迈克尔·布朗斯坦 ， 琼·布鲁纳 ， 塔科·科恩 ，以及佩塔尔·韦利奇科维奇

维数灾难是高维学习中的一个重大挑战。这篇帖子讨论了几何深度学习(GDL)中的一个主要概念，叫做几何先验 ( 对称性和尺度分离)。几何先验让我们能够解决维数灾难。我们将看到域和信号在这些域上的概念。我们将讨论对称以及它们是如何出现在 ML 和 GDL 中的。我们还将看看一些潜在的数学概念，如群表示、不变量和等价变量。最后，我们将看到如何在深度学习中使用不变&等变网络的概念。

立方体的旋转对称性(Oₕ).群图片来自 GDL 课程第三讲。

参见我上一篇文章 s 上一篇 Erlangen 程序的 ML 和 高维学习 。这些博文是根据四位导师的 GDL 原书 和 GDL 课程atAMMI改编的。

在上一篇关于高维度学习的文章 中，我们看到，由于维数灾难，在没有假设的情况下，高维度的学习是不可能的，也就是说，我们的学习问题所需的样本数量随着维数呈指数增长。我们还介绍了主要的几何函数空间，其中我们在高维空间中的点可以认为是低维几何域上的信号。从这个假设出发，并且为了使学习易于理解，我将呈现对称性(在这篇文章中)和标度分离(在下一篇文章中)。

维度的诅咒。图片来自 GDL 课程第三讲。

此外，我们还讨论了我们需要注意的三种误差，即*近似误差、统计误差和优化误差。*如果我们的函数类减少(我们试图估计的真实函数远在这个类之外)，近似误差会增加，这表明有一个大的函数类。相比之下，统计误差意味着我们不可能基于有限数量的数据点找到真正的函数。随着函数类的增长，这种错误会增加。最后，我们有优化误差，它与我们在假设类中找到最优假设的程度有关(当经验风险最小化可以有效地解决时，误差很小)。

我们将会看到，通过使用几何先验和等变网络:尊重我们问题的对称性的网络，我们可以减少我们假设类的大小或复杂性。通过这样做，我们可以在不放弃有用假设的情况下减少统计误差，这意味着近似误差不会更差。(普适性结果:如果我们的函数类是等变函数的普适逼近器)。

几何域

在几何深度学习(GDL)中，数据存在于由ω表示的域中。域包括网格、组、图和流形(量规&测地线)。定义域也是一个集合(在典型的数学中)，但是它可能有一种不同的结构。我们可以有一个邻域结构(在网格和图形中)。我们也可以有一个出现在流形中的度量结构(点之间的距离),因为我们可以测量表面上点之间的距离。

几何领域:5 Gs。图片来自 GDL 课程第三讲。

GDL 的关键信息是，我们希望设计处理几何数据的神经网络，这些网络应该尊重域结构，无论它可能有许多选择。

几何域上的信号空间

在大多数情况下数据不是域本身，而是域上的信号。ω域上的信号是一个函数𝓍，它输入ω域上的一个元素，并将其作为向量空间𝒞中的一个向量输出，即𝓍∶ω→𝒞(𝒞的维数称为通道)。我们也可以把ω上的𝒞-valued 信号空间定义为所有函数(该空间上的各种信号)的集合，即𝒳(ω，𝒞)={𝓍∶ω→𝒞}.

例如，如果我们的域是一个网格，数学上ω=ℤₙ×ℤₙ，其中ℤₙ是从 0 到 n-1 的整数集合。该信号是采用诸如(0，1)的网格点并将其映射到 rgb 值的空间，即信号𝓍∶ω→ℝ)的函数。类似地，在图中，我们的域是从 1 到 n 的节点集。对于分子图，信号将是从该图到𝑚维度的向量空间的函数，其中𝑚是我们关心的原子的数量(对于每个节点，我们分配一个热向量，该向量告诉该节点的原子类型)。

网格和图形域中的信号。图片来自 GDL 课程第三讲。

希尔伯特空间结构

正如所见，域可以是一个向量空间 n ，比如用于建模图像的连续平面ℝ。所以我们可以在这个域里加点。但是在某些领域，我们没有这个向量空间(例如，在图中，我们没有任何节点添加的概念)。尽管如此，信号空间总是有一个向量空间结构，在这个结构中，我们可以将两个信号相加并乘以标量。从数学上讲，这由以下等式定义:

信号的加法和标量乘法。来自 GDL 原型书第三章的方程式。

换句话说，如果我们有两个信号𝓍和𝓎(两个图像)，并且我们有两个标量 α & β ，我们可以形成信号 α 𝓍 + β 𝓎，并在点𝑢.对其求值这种加法是逐点定义的(在𝑢评估第一个函数，在𝑢评估第二个函数，将这些数与 α 和 β 线性组合，并将它们相加)。

两个信号的相加。图片来自 GDL 课程第三讲。

如果信号空间有无穷多个点，那么它就是一个向量空间。那么，这将是一个无限维的向量空间(就计算而言，我们必须处理有限的东西)。而且，在这个向量空间中定义一个内积会给我们希尔伯特空间。

给定一个内积〈 𝓋，𝓌〉on 𝒞 (𝓋，𝓌 ∈ 𝒞)和ω上的一个测度𝜇(关于它我们将定义积分)，我们通过下式得到𝒞𝒳(ω上的一个内积:

𝒞𝒳(ω上的内积定义)(关于测度𝜇).的积分)来自 GDL 原型书第三章的方程式。

测量𝜇允许测量ω的集合和子集的大小。在大多数情况下，这种数学设置归结起来相当简单。具体来说，如果ω是一个有限集，那么 计数度量 可以用来计算积分，这意味着我们对ω上的点求和，所以当我们做内积时，我们只对向量空间的维数求和。

我们需要内积的原因是我们想做模式匹配当我们使用线性/conv 层时，我们需要以某种方式比较我们的信号𝓍和滤波器𝓎.

几何特征领域

到目前为止，我们已经看到信号𝓍是𝒞).域ω上的函数(取ω的一个元素𝑢并将其映射到向量𝑥(𝑢)稍微一般一点的东西，在物理学上叫做场，在数学上叫做(“丛”)。字段是一个函数的略微一般化，其中不是将ω中的𝑢映射到固定空间𝒞中的向量𝓍(𝑢，而是将其映射到向量空间𝒞的元素，该元素由𝑢进行子索引，即𝓍(𝑢) ∈ 𝒞 ᵤ 。例如，网格上的向量场为流形中的每个点分配一个向量，该向量位于该点的切平面上。所以流形上的每个点都有一个切面。

流形的向量场和切空间。图片来自 GDL 课程第三讲。

为了简单起见，我们现在只讨论𝒞𝒳(ω的函数空间。

域为数据

我们有时候域上没有信号，但是域本身就是数据。这出现在网格、没有节点或边特征的图形以及点云中。有一些方法可以处理这个问题，比如我们可以用图邻接矩阵作为ω×ω上的信号。在流形或网格的情况下，我们可以使用*度量张量，*网格固有的量，它可以被视为信号。

域本身就是数据。图片来自 GDL 课程第三讲。

对称性

这里是我们岗位的核心。为了解释物体的对称性，让我们以三角形为例。我们可以对三角形做几个不改变它的变换。我们可以将它旋转 120 度，这将只是用另一个点替换三角形中的每个点，而不会改变整体形状。我们也可以在垂直线翻转它，它不改变三角形。这些变换叫做对称。

三角形的对称性(由 120 度旋转 R 和翻转 F 产生)。图片来自 GDL 课程第三讲。

一个对象的对称性是该对象保持不变的变换。

我们将看看 ML 中出现的两种对称，参数化的对称和标签函数的对称。在参数化的对称性中，设输入空间为𝒳，输出空间为𝒴(一个标签空间)，权重空间为𝒲(一个参数空间)，模型为𝑓，一个从输入空间𝒳和权重空间𝒲到标签空间𝒴的映射，即f*:𝒳×𝒲→1。我们可以说，将权重𝓌映射到其自身的变换𝔤，即𝔤: 𝒲 → 𝒲，是参数化的对称，如果:*

参数化的对称性。来自 GDL 课程的方程式，第三讲。

比如交换神经网络同一层的两个神经元的进出连接，并不会改变输入输出映射 f (。, 𝑤).

参数化的 S 对称性示例(s 交换两个神经元的输入和输出连接)。图片来自 GDL 课程第三讲。

ML 中出现的第二种对称，也是最会讲的，是标签函数的对称。再次假设输入空间𝒳、输出空间𝒴和标签函数 L 是地面实况标签函数，即从输入空间𝒳到输出空间𝒴.的映射变换𝔤是对称的，如果:

标签函数的对称性。来自 GDL 课程第三讲的方程。

换句话说，如果我们进行转换，然后计算标签，这与只计算标签是一样的。例如，在图像分类任务中，𝔤可以是图像的旋转、平移或缩放。

标签函数的 S 对称性示例(如果我们旋转图像，该类不会改变)。GDL 课程，第三讲。

学习类≅学习对称性

我们可以把所有的学习问题，或者至少是所有的监督学习，看作是关于对称性的学习。假设我们有两个类的输入空间𝒳。我们可以看到，𝒳中任何尊重类边界的可逆映射(它将某类中的一个点映射到该类中的另一个点)都是标号函数 L 的对称。这意味着如果我们知道了所有的对称，我们只需要为每一个类标注一个标签！因为我们可以从一个点开始，应用一个对称很多次，到达这个类中的每一个点。

具有两个类的输入空间𝒳(如果我们知道所有l 的对称性，从某个类中的一点，我们到达该类中的每一个其他点)。图片来自 GDL 课程第三讲。

如果学习问题是非平凡的(如果我们相信我们需要学习来解决我们的问题)，我们不应该期望能够先验地找到完全对称群。

结构域的对称性

在一大类问题中，我们称之为几何深度学习问题，对称性的来源来自于数据赖以生存的域。设ω表示一个几何域。变换𝔤∶ω→ω是对称的，如果它保持ω的结构。这些对称性的一些例子可能是:

• 集合元素的排列保持了集合的成员关系。如果我们的定义域是一个集合，那么这个集合的结构就是一个集合成员关系(无论元素是不是集合的一部分，都不考虑集合中元素之间的顺序或关系)。然后，这个集合的元素的排列保持了这里的结构。
• 欧氏等距(旋转、平移、反射)保持了欧氏空间中的距离和角度ω*=ℝᵈ*。如果我们有欧几里得空间，比方说一个平面， ℝᵈ是一个度量空间(点与点之间的距离是空间结构的一部分)。然后，对称可以是欧几里得等距，即，诸如旋转、平移和反射的距离保持映射。
• *一般微分同胚保持流形ω的光滑结构。*如果我们把 ℝ 或者任何流形仅仅看作一个光滑流形，我们不会以一个度量结束，所以域的任何微分同胚或者光滑翘曲都可以认为是一个对称性。

几何深度学习中对称性的例子。GDL 课程，第三讲。

对称群

给定物体的所有对称的集合称为对称群。

现在我们可以观察到一些关于对称性的东西:

• 恒等变换永远是对称的。
• 给定两个对称变换，它们的合成(一前一后做)也是对称的。
• 给定任何对称，它的逆也是对称。

抽象群体

群是一个集合 𝔊 连同一个二元运算 ⚬:𝔊×𝔊→𝔊 称为组成，表示为 𝔤⚬𝔥 = 𝔤𝔥(为简单起见)，并满足以下公理:

• 关联性 : (𝔤𝔥) 𝔨 = 𝔤 (𝔥𝔨)适用于所有𝔤、𝔥、𝔨 ∈ 𝔊.
• 恒等式:存在唯一的𝔢 ∈ 𝔊满足𝔤𝔢 = 𝔢𝔤 = 𝔤.
• 逆:对于每个𝔤 ∈ 𝔊都有一个唯一的逆𝔤⁻ ∈ 𝔊，这样。𝔤𝔤⁻ = 𝔤⁻ 𝔤 = 𝔢.
• 结束:对于每一个𝔤，𝔥 ∈ 𝔊，我们有𝔤𝔥 ∈ 𝔊.

为了简单起见，通常使用𝔊来指代组，而不仅仅是其元素的集合。注意，一个群是而不是必然可换的，即𝔤𝔥≠𝔥𝔤一般。具有交换运算的群称为阿贝尔群。

立方体的旋转对称性(Oₕ).群图片来自 GDL 课程第三讲。

对称组，抽象组&组动作

在这一节，我们定义群论中的几个概念:

• 对称群是其元素为变换𝔤∶ω→ω的群。
• 分组运算是贴图(函数)的合成。
• 抽象组是一组元素以及复合规则，满足组公理(如上定义)。
• 组动作是一个映射或函数，它获取组元素和域ω的元素，并将其映射到域ω中的新元素。，即𝔊×ω→ω(表示为(𝔤，𝑢) ↦ 𝔤𝑢).群作用应满足:𝔤(𝔥𝑢) = (𝔤𝔥) 𝑢 和 𝔢𝑢 = 𝑢，对于所有𝔤，𝔥 ∈ 𝔊和∈ω。

群体作用的例子:作用于ℝ的欧几里得平面运动。 ℝ 中的欧几里得群具有带角度𝜃的旋转分量以及 x 和 y 方向的平移分量(tₓ，tᵧ)，可应用于二维点(𝓍，𝓎).

作用于ℝ的欧几里得平面运动。 GDL 课程第三讲。

凯莱图表&表格

如果我们的群是有限且离散的，我们有两种方式来表示这个群的结构:凯莱图&表。让我们以我们以前讨论过的三角形的对称群为例。在 Cayley Diagrams 中，每个节点都是对象(三角形)的一个配置，但是我们也可以用一个组元素来标识它。我们可以对任何节点应用旋转或翻转，并移动到其他节点。相比之下，我们可以把群结构表示成一个乘法表。组中的每个元素都列在一行和一列中，表中的条目将是组元素的乘积。

凯莱图和乘法表(三角形的对称群)。GDL 课程，第三讲。

组的种类

这里有一些不是有限的也不是离散的群。大体上，我们可以将群定义为两种类型，离散群和连续群。 离散组可以是:

• 有限的，如旋转对称的三角形。
• 无限，或可数无限，如整数集合的平移和无限扩展网格。

另一方面，我们有连续组，它们可以是:

• 紧凑组，如二维旋转。
• *局部紧群，*如连续平移和旋转平移。
• *非局部紧群，*如微分同胚的无限维群(作为数学句柄最难得到)。

所有这些类型的群都可以是可交换的，例如顺序无关紧要的二维旋转群，或者是不可交换的，例如顺序有关的三维旋转群。

各种团体。GDL 课程，第三讲。

ω作用于信号𝒳(ω的对称性，𝒞)

作用在域ω上的对称性也作用在ω上的信号上。对于群元素𝔤和信号𝓍，我们将𝔤对𝓍的新群作用定义为:

𝔊对𝒞).𝒳(ω信号的行动 GDL 原型书，第 3.1 节。

简而言之，在域中的点𝑢评估的变换信号等于在𝔤⁻ 𝑢.评估的原始信号𝓍

这方面的一个例子可以是翻译图像。𝔊是一个翻译小组，由𝑥和 y 两部分组成。在下图中，我们看到 bug(在右边的原图中)从左边移到了右边( t 是一个正数，沿 x 方向平移)。𝔤.翻译形象的价值在点𝑢的𝓍将是在点𝑢 - t，即 𝔤.的原始图像中的相同值𝓍(𝑢) = 𝓍(𝑢 - t 。注意，如果没有逆(减)，我们就无法满足群体行动的公理。

翻译图像(右边:原始图像，左边:翻译图像)。GDL 课程，第三讲。

如前所述，域可能缺少向量空间结构，但信号的空间总是具有向量空间结构。此外，我们刚刚在信号上定义的组动作是线性(如果我们有信号的线性组合，那么我们应用变换𝔤，这与对单个信号应用𝔤并线性组合它们是一样的)。

信号群作用的线性。GDL 原型书，第 3.1 节。

信号上的群组动作的线性示例(两个图像)。GDL 课程，第三讲。

小组陈述

群𝔊的 n 维实表示是一个映射𝜌 ∶ 𝔊 → ℝⁿ*ⁿ，给每个𝔤 ∈ 𝔊指定一个可逆矩阵𝜌(𝔤，并且满足𝜌

例如，设群𝔊是 1 维中整数平移的集合，即𝔊 = (ℤ,+)，定义域ω=ℤ₅= { 0，1，2，3，4}(如短音频信号)，𝔤 = 𝑛对 u∈ω的作用定义为:𝑛，𝑢↦𝑛+5(mod 5)。那么，在𝒳(ω(5 维空间)上的表示就是移位算子:

团体代表的例子。GDL 课程，第三讲。

要验证这是群表示，我们得证明矩阵𝜌( n 可逆，𝜌 (n+m) = 𝜌 *(n)P(m)(群表示的条件)。*很明显，𝜌(n 的行列式)，即 det(𝜌( n )，不等于 0；(det(𝜌(n)=det(𝜌(1))ⁿ，det(𝜌(1))不能为零，因为单位矩阵是𝜌(1)).的子矩阵因此，𝜌( n 是可逆的。𝜌(n+m)= 𝜌(n)𝜌(m)直截了当。

重要的是要注意，如果一个组满足条件，我们可以对它有多个表示。例如，如果我们有旋转群 SO(2)，我们可以有一个作用于 ℝ 中的点、位于圆上的点、平面上的函数/信号等的表示。在 GDL 蓝图中，每个特征空间都有一个组表示，每个组表示都可能不同。

在下面几节中，我们将讨论不同域中对称群的一些例子。

对称性:集合&图形

S 从集合和图开始，对称群的一个例子是置换 s 的群，即𝔊 = 𝕊ₙ(“对称群”)，定义域是顶点的集合，即ω=𝑉，(可能还有边)。图形的特征可以是标量、矢量或张量(𝕊ₙ).的三组表示

三种图形特征(𝕊ₙ).的表示 GDL 课程，第三讲。

此外，我们有两种作用于图和集合的对称；描述的对称性和对象本身的对称性。首先，请注意:

• 一个图(甚至一个集合)是一个抽象对象(之前定义过)。
• 实际上，我们所能访问的只是计算机内存中的图/集合的描述(它有一些无关的属性，比如图/集合元素的节点顺序)。
• 描述具有属性(例如顺序),这些属性不是对象的固有属性。

通常，我们感兴趣的是描述的对称性，而不是物体本身的对称性。为了解释这一点，让我们讨论下图。对于中间的图，如果我们为每个节点选择一个标签，我们可以将图写成邻接矩阵。如果我们应用一个置换，我们可以得到一个不同的邻接矩阵(右边的那个)。这里是描述的对称性；我们对同一个图有不同的描述(我们希望我们的网络是不变的或等变的)。我们也可以看到这个图形本身有一种对称性。如果我们做另一个置换，例如交换节点 3 和 4，邻接矩阵将是相同的(左边的那个)。

D 描述的对称性和 O 对象的对称性的例子。GDL 课程，第三讲。

对称性:网格

在网格中，定义域是网格点的集合，和对称群𝔊可以是离散 平移、旋转、翻转。底层的连续空间可能有更多的对称性(例如缩放)。群体表征是有规律的；我们可以移动和旋转图像。

常规组代表。GDL 课程，第三讲。

图像的旋转和平移。GDL 课程，第三讲。

对称性:群&齐次空间

我们可以将对称性推广到一般群和齐次空间的集合，其中域ω≅𝔊/ℌ，对于空间中的每两个点，至少有一个对称将一个点映射到另一个点(将在后面的帖子中涉及)。群𝔊是一个局部紧群。这个群的作用是传递的:对于任何𝓍，𝓎∈ω，存在𝔤 ∈ 𝔊使得𝔤.𝓍 = 𝓎.小组代表是正规的(𝜌(𝔤) 𝑥(𝑢) = 𝑥(𝔤⁻ 𝑢)).

一般群和齐次空间。GDL 课程，第三讲。

对称性:流形(测地线&量规)

最后，我们有流形(测地线&量规)，这也将在后面的帖子中详细介绍。群𝔊是群规范变换，它指的是特征空间的参考系中的变化。为了解释这一点，让我们看看下图。我们有一个向量场，它给流形上的每一个点分配一个切平面上的向量。如果我们想用数字表示矢量，我们必须在切平面上选择一个坐标系。改变这个框架将改变矢量的数值表示(这被称为规范变换)。制图表达𝜌取决于要素类型。

流形和网格。GDL 课程，第三讲。

到目前为止的总结

到目前为止，我们已经讨论了各种几何域上的信号空间及其希尔伯特空间结构。我们还看了几何先验中的几个概念，如对称性、对称群和群表示。有趣的问题将是我们如何在深度学习中应用这些概念，以及我们如何建立包含它们的模型？

不变表示法

A 假设我们想对字母“A”、“B”和“C”进行分类，它们可能以不同的字体和不同的方向出现(见下图)。然后，如果我们知道标签函数的对称性(在这种情况下，数字的旋转)，我们可以形成一个不变表示，其中同一字母的每个旋转版本由一种特征向量表示。因此，我们的学习问题变得更容易。

数字分类问题中的不变表示。GDL 课程，第三讲。

然而，在深度学习中，过早地构建不变表示是不明智的。原因是:为了识别整个对象，我们需要首先识别这个对象的部分(这就是为什么神经网络应该是深度的)。但是如果我们让的中间表示不变，我们将丢失关键信息(看下图)。

零件的不变表示将丢失关键信息(见右图)。GDL 课程，第三讲。

等变网络

T 上述问题的解决方案是使用等变网络，其具有以下成分:

• 特征向量空间𝒳ᵢ。
• 网络图层(地图)𝑓ᵢ。
• 对称群𝔊。
• 𝔊的𝜌ᵢ代表每个𝒳ᵢ。

最简单的例子可能是平面图像、RGB 输入( n × n 像素和 3 个通道)。特征空间𝒳₀将是 n × n × 3 。所以我们有一个 n × n × 3 维的𝔊群表示，它可以是平移。𝒳₁可能是具有 64 个通道的卷积层𝑓₁的输出(我们的特征空间将是 n × n × 64 )。这将是同一组𝔊的不同表现(翻译)。直到层数为止。如果满足下面的等式，我们说网络是等变的

等方差方程。GDL 课程，第三讲。

我们可以从下图了解一下。如果我们从𝒳₀开始，使用𝜌₀进行变换(例如，移动图像)，然后应用第一层。如果我们首先应用图层，然后应用相同的变换𝜌₀，我们应该会得到相同的结果，但是现在通过特征空间中的表示来产生𝒳₁.如果所有的网络层都满足这个性质，我们说这个网络是等变。我们还可以证明，如果每一层都满足等方差，那么它们的合成也满足等方差。

等变网络图。GDL 课程，第三讲。

等方差作为对称一致的推广

L et 重复一下我们之前讨论过的分类问题(对数字“A”、“B”、“C”进行分类)。输入空间𝓍和𝓎(在下图中)是位于轨道上的两个信号(如果我们考虑字母“a”的所有旋转副本，这将形成称为“a”的轨道的流形)。网络将该输入映射到由通用字母“A”符号化的特征向量。在非等变网络中，旋转版本可以被发送到特征空间中的不同点(这不可能发生，等变网络必须在整个轨道上一致地概括)。此外，等变网络已经以一种方式被一般化，即符合对称性*。换句话说，如果𝓍和𝓎应该被映射到同一点，那么它们的变换版本也应该被映射到同一点。*

非等变网络(这不可能发生，等变网络必须在整个轨道上一致地推广)。GDL 课程，第三讲。

当我们深入研究对称和等变网络时，你可能想知道为什么不使用数据增强，它可以取代对称？

当然，数据增强在 ML 中应用广泛，简单易行，非常流行，但也有各种优缺点。支持等方差的最突出的论点，在某类问题中，我们有对称性，比如医学图像，真的有平移和旋转。在这些情况下，我们看到等变网有更好的性能。另一个相关的例子是，如果我们有一个大的对称群(如置换群)，实际上不可能使用数据扩充。

T 他的帖子讨论了几何域、信号，以及几何先验的第一要素(对称性)。我们还看到了对称背后的几个数学思想，以及几何领域的例子，以及它们如何用于深度学习。下一篇文章将介绍几何先验的第二个要素(尺度分离)。然后，我们将看到完整的几何深度学习 ( GDL)蓝图。

参考

• GDL 球场，( AMMI ，2021 年夏季)。
• T. Cohen 的第 3 讲[ 视频 | 幻灯片。
• 米（meter 的缩写））m .布朗斯坦、j .布鲁纳、t .科恩和 p .韦利奇科维奇，几何深度学习:网格、组、图形、测地线和量规 (2021)。

我很感谢拉米·阿迈德、汉尼斯·斯特尔克和 M·埃尔法提赫·萨拉赫审阅了这篇文章并提供了有益的评论。

几何先验 II

原文：https://towardsdatascience/geometric-priors-ii-6cd359dfbcf4

几何深度学习

GDL 蓝图

一系列博文，总结了 AMMI 计划的 几何深度学习(GDL)课程； 非洲机器智能硕士 ，授课老师 迈克尔·布朗斯坦 ， 琼·布鲁纳 ， 塔科·科恩 ，以及佩塔尔·韦利奇科维奇

维数灾难是高维学习中最具挑战性的问题之一。在之前的一篇帖子 ( 几何先验 I ) 中，我们讨论了机器学习和几何深度学习中的各种概念，包括对称、不变和等变网络。在这篇文章中，我们将这些数学概念与维数灾难联系起来。然后，我们引入第二个几何先验；尺度分离，代表了战胜维数灾难的有效工具。最后，我们将看到如何将这两个几何先验(对称性和尺度分离)结合起来，构建几何深度学习(GDL)蓝图，它可以作为不同几何域的通用架构。

尺度分离先验(图像分类)。图片来自 GDL 原型书，第 3.4 节。

本帖与M . Elfatih Salah合著。另见上一篇 s 关于 Erlangen 程序的 ML ， 高维学习 ， 几何先验 I 。这些博文是根据四位导师的 GDL 原书 和 GDL 课程atAMMI改编的。

R 调用以前帖子中的学习设置，在统计学习中，我们有一个目标函数 f *，一个从输入空间𝒳到输出空间 ℝ 的映射； f * : 𝒳 *→ ℝ.*我们也把输入空间𝒳写成一组信号𝒳(ω，𝒞)={𝓍∶ω→𝒞}，其中ω是几何域，𝒞是定义在ω上的向量空间。在 几何先验 I 中，我们已经看到了几何域的各种例子，以及如何使用域ω来捕捉与该域中信号的性质密切相关的各种变换。(本帖与 几何先验 I ) 。

这样，从前面介绍的内容来看，目标函数 f 定义为f***:*𝒳(ω，𝒞) → ℝ 。我们的目标是在假设类𝓕中学习这个函数，我们可以按照我们想要的方式参数化它，例如使用神经网络。我们现在的承诺是目标 f *是𝔊-invariant ( f *保持不变，如果考虑变换群𝔊的任何元素和任何信号‘输入’的话)。数学上，

𝔊-Invariant 函数。GDL 课程，第四讲。

因此，自然的问题是:我们如何在假设课上利用这个承诺？这足以打破维度的诅咒吗？

为了回答第一个问题，让我们定义一个组平滑算子，它是一个采用任何假设 f 并用组𝔊.中所有变换的平均值替换它的算子在一个数学设置中，如果𝔊是一个离散的有限群并且 f 是我们的假设，我们将𝑓的𝔊-smoothing 算子定义为，

函数𝑓.的𝔊-Smoothing 算子来自 GDL 课程的方程式，第四讲。

如果我们仔细观察，可以发现函数 f 的这个平滑算子有一个有趣的性质；它也是𝔊-invariant(如果𝑓已经是𝔊-invariant，那么应用平滑运算符不会改变它)。因此，给定假设类𝓕，我们可以通过应用平滑算子使其成为𝔊-invariant，并且更正式地如下:

假设类𝓕.的𝔊-Smoothing 算子来自 GDL 课程的方程式，第四讲。

(同样，如果我们扩展不变性群𝔊，不变性类将变得更小，这似乎是直观的)。

不变性下的学习

通过将平滑算子应用于我们的假设类，我们可以证明近似误差不受该平滑操作的影响，并且我们可以这样正式表述:

近似误差不受平滑操作的影响。GDL 课程，第四讲。

为了证明，我们可以注意到平滑算子是在 L (𝒳).)中的正交投影算子所以预测中产生的 L 误差可以分解为该算子图像中的误差加上正交补中的误差，如下所示，

近似误差不受平滑操作的影响(证明)。GDL 课程，第四讲。

然后，用 f *替换𝔤，知道 f *不受平滑操作影响，我们得到，

近似误差不受平滑操作的影响(证明)。GDL 课程，第四讲。

另一方面，由于平滑后假设类更小，统计误差减少，但问题是:多少？

看一看假设类之一，我们假设它是 Lipschitz 类，我们知道它具有以下属性:

李普希茨函数类。GDL 课程，第四讲。

然后，其平滑版本具有以下属性:

应用𝔊-Smoothing 算子时的李普希茨函数类。GDL 课程，第四讲。

平滑版本使得 Lipschitz 假设更强，因为距离更弱，并且主要结果:

定理[ Bietti，Venturi，B.'21 ]:利用一个𝔊-invariant 核岭回归，学习一个李普希兹，𝔊-invariant 函数 f 的泛化误差满足𝔼𝓡(f̂)≲*θ(|𝔊|n)⁻/ᵈ，其中 n 是样本数，d 是维数。

该定理导致样本复杂度的量化增益，该增益与组的大小成比例，因为泛化误差由样本的数量乘以组的大小来限定，都是 -1/d. 的幂

此外，组大小|𝔊|在维度上可以是指数的，例如，如果我们考虑局部转换(可能的局部转换的数量在域的大小上是指数的)。然而，正如我们所看到的，速率仍然受到维数的诅咒，这表明引入群不变性不太可能打破维数的诅咒。

到目前为止的结论:

到目前为止，我们已经看到使用全局对称性或不变性先验减少了统计误差并保持近似误差不变，但是仍然不足以打破维数灾难。所以问题来了:打破维度诅咒所必须的对称性之外的假设有哪些，我们遗漏了什么？

还有，我们一直在讨论理论问题，没有讨论如何高效地构建不变类。我们所描述的平滑操作符并不是我们想要在计算机中直接实现的东西(因为计算所有转换组的平均值是非常不切实际的)。从基本原则出发，我们如何高效地做到这一点？

深度学习“归纳偏差”:组合性

深度学习研究领域有影响力的人可能会在不同的演讲中提到，比如 Yann LeCun 和 Y oshua Bengio ，深度学习之所以有效，是因为 组合性 *原理。*此外，有来自大脑的证据表明，大脑也是分层组织的，这种想法并不是深度学习研究人员开发的，但它已经存在很长时间了。

深度学习为什么有效？Y 奥舒亚·本吉奥(上)和扬·勒村(下)。GDL 课程，第四讲。

既然我们似乎都同意组合性对深度学习的成功至关重要，我们应该问我们如何才能形式化这种直觉？

多尺度结构

形式化这种直觉的一个方法是观察物理学和计算科学的不同领域。在不同的领域，我们可以看到多尺度结构是一个基本原则和许多学科的基础。

在计算生物学中，为了知道生物学是如何工作的，我们必须理解不同尺度下的过程。有些事情我们可以用分子动力学或者功能图来理解。 比如湍流 ，展现了数学物理中不同尺度之间的交流。 [逾渗](https://en.wikipedia/wiki/Percolation#:~:text=Percolation%20(from%20Latin%20percolare%2C%20%22,is%20described%20by%20Darcy’s%20law.) 另一种模型，描述流体通过多孔材料的运动。因此，我们必须考虑不同的尺度。你可能也有自己的例子，在不同的尺度下，系统有不同的行为。

逾渗模型。来自[维基百科](https://en.wikipedia/wiki/Percolation#:~:text=Percolation%20(from%20Latin%20percolare,%20%22,is%20described%20by%20Darcy’s%20law.)的动画。

多分辨率分析基础

在我们的上下文中，我们希望将多尺度结构的概念与几何域(网格、组、图形和流形)联系起来。另一个非常相关的概念是多分辨率分析(MRA)。(为了简单起见，我们将假设一个网格域，即ω是一个 2D 网格，但大多数东西可以更一般地描述，然后应用于所有的几何域)。

多分辨率分析是一种方案，该方案试图根据存在于较小域 ( 较粗域 域)中的信号来分解存在于该域中的信号。例如，在网格域的情况下，我们想要更少的像素。在多分辨率分析中，我们试图理解如何从一个分辨率(特定数量的像素)到下一个分辨率(减去数量的像素)，以允许我们保留信息。一般来说，我们可以将一幅图像分解成不同分辨率(分辨率更低)的同一幅图像，加上我们需要返回到第一幅图像的细节(见下图)。

图像分解。GDL 课程，第四讲。

在这里，我们必须声明，通过降低网格的分辨率，我们正在攻击维度的诅咒。这是由于在这种情况下像素数量的指数依赖性。我们可以保证，在某一分辨率下，由于样本数量与新分辨率成正比，因此将避免维数灾难。因此，我们的问题将是如何将多尺度结构的想法与学习中的先验结合起来？

秤分离之前

L 让我们讨论一些如何在学习中结合多尺度结构和先验知识的例子。想象一下，目标函数 f *可以在图像的粗糙版本上定义。在下图所示的分类问题中，我们只想知道图像是“海滩还是“山”而不是它的分辨率。

两幅图像的原始版本(左)和粗糙版本(右)(分类问题)。GDL 课程，第四讲。

在这个问题中，我们可以将图像的分辨率降低到能够给出正确分类的某个水平。如前所述，这将减少网格的大小(像素数)。因此，通过粗化，我们可以避免维数灾难。然而，这通常是不正确的，因为我们需要非常小心粗化的程度。如下图所示，太粗化会丢失太多信息，我们无法分辨图像中到底有什么。

过度粗化会丢失图像信息的示例。GDL 课程，第四讲。

让我们看看多尺度结构也能有所帮助的另一个例子。假设目标函数 f* 可以通过局部项的和来近似，

用局部项的和来逼近函数 f* 。GDL 课程，第四讲。

局部术语 x(u) 是指提取以像素 u 为中心的补丁图像。比如在识别一个纹理(下图所示的分类问题)时，如果我们在一堵砖墙或者一堆草里。我们当地的一个术语就属于这一类。因此，如果我们想要一个好的分类，我们可以只取局部描述符的平均值，因为在粗尺度上没有很多可变性(因为纹理是一个具有空间同质性的结构)。

本地术语的示例:黄色块显示的补丁图像。GDL 课程，第四讲。

我们可以看到，维数灾难也被避免了，因为这个问题中的相关维数减少到了面片维数。

我们还必须声明，这是一个强有力的假设，通常是不充分的，如下图所示；所有的本地术语都可能给出错误的分类。

当本地术语可能给出错误分类时的例子。GDL 课程，第四讲。

重要的是要说明这样一个事实，我们可能有原始尺度和粗略尺度的信息，但不知何故它们以一种不太极端的方式相互作用。因此，在某种意义上，我们必须将这两种尺度结合起来。

让我们考虑一个通用的组合模型，

具有两个算子(非线性局部粗化和非线性假设)的一般成分模型。GDL 课程，第四讲。

函数 f 由两个运算符组成；非线性局部粗化结合非线性假设在粗尺度上提取一些信息(比单独平均小块强得多)，两者都是可以学习的。你可能会问为什么这个组合和什么时候可以论证这个组合模型更有效率*。当我们考虑一些具体的例子时，比如动态编程和分治算法，我们是从更小的问题的角度来看待这个问题的。在这些情况下，我们可以认为这种组合模型更有效。然而，这仍然是一个未解决的问题，并且从理论的角度来看，这个组成模型还没有被完全理解。

接下来，我们结合尺度分离和群不变量这两个先验，从第一性原理给出一个强有力的模型。然后，我们介绍完整的几何深度学习(GDL)蓝图。

不变性与尺度分离的结合

S 由于我们有一个函数类是对群体行为不变的(𝔊-invariant 函数)，我们想把它和上面描述的多尺度结构(尺度分离之前)合并成一个架构。因为我们讨论过的这些几何先验给了我们模型中应该有的必要条件，而不是一个特定架构的所有要素。

如果我们从一个对群作用不变的线性算符(线性𝔊-invariant 函数)开始，使用前面介绍的𝔊-smoothing 算符，我们可以写成，

带有𝔊-Smoothing 算子的线性𝔊-Invariant 函数。GDL 课程，第四讲。

其中 x̅ 是群轨道上的平均值(我们将 x 的群平均值称为 A x )。(第一个不等式是因为 f 是𝔊-invariant，我们在应用𝔊-smoothing 算子，第二个不等式是因为 f) 的线性。

然而，如果我们仅仅依靠这个假设，我们会丢失很多信息。因为如果我们有一个平移组，例如，一个图像的轨道将由该图像所有可能的平移组成，然后我们只需将所有这些图像平均为一个。

作为补充，我们可以使用线性𝔊-equivariant 函数，我们在 几何先验 I 中讨论过。(如果𝔊同时作用于𝒳和𝒴，那么映射𝐵 : 𝒳 → 𝒴就是𝔊-equivariant。𝑥) = 𝔤.𝐵(𝑥)).为了实现一个更强大的模型，我们可以用一个元素式非线性函数 𝜎合成这个等变函数𝐵；𝜎 : 𝒳 → 𝒳，𝑤𝑖𝑡ℎ𝜎𝑥(u)= 𝜎(𝑥(u)，以给出一个非线性𝔊-equivariant 函数 u，即 U := 𝜎 ⚬ 𝐵(由于非线性是逐元素应用的， U 将保持𝔊-equivariant).

然后，我们可以通过组合组平均值 A 和非线性𝔊-equivariant 函数u来提供强大的𝔊-invariant 函数

一个迫切的问题是，这种组合是否会产生一个富裕的阶层？换句话说，我们能很好地逼近任何一个目标函数，𝔊-invariant，并且有足够的 b 和𝜎选择吗？

参考 通用逼近定理 ，这种组平均与非线性等变映射的组合产生了通用逼近器。然而，为了更稳定的表示，我们需要使用一个局部等变图。另外，如果我们有一个具有长程相互作用的函数，我们不能用单层的局部𝔊-equivariant 映射来近似它；相反，我们必须组合几个局部等变函数，使区域变粗。(为证明， GDL 原典 ，第 3.5 节)。

最后，为了构建一个丰富的𝔊-invariant 函数类，我们需要一个局部等变映射、一个全局不变映射和一个粗化操作符。

几何深度学习(GDL)蓝图

L 很快，我们可以在所谓的 GDL 蓝图中，为学习目标函数 f 所需的模型勾勒出上面讨论的所有成分。GDL 蓝图是一个建设性的方法，我们可以用它来定义一个通用架构*，它可以应用于各种几何领域(网格、组、图形和流形)。

设ω是一个定义域，𝔊是ω上的对称群，ω′是ω的一个更粗的定义域，即ω′⊆ω。几何深度学习(GDL)蓝图的构建模块是:

• 线性𝔊-equivariant 层𝐵∶𝒞𝒳(ω)→𝒞′𝒳(ω′)，满足𝐵(𝔤. x ) = 𝔤.𝐵(𝑥)适用于所有𝔤 ∈ 𝔊和𝑥∈𝒳(ω、𝒞).
• 非线性𝜎∶𝒞→𝒞′按元素应用 as(𝝈(x)(u)= 𝜎(𝑥(u))。
• *局部汇集(粗化)*𝑃∶𝒳(ω，𝒞)→𝒳(ω′，𝒞)，使得ω′⊆ω。
• *𝔊-invariant 层(全球统筹)*𝐴∶𝒳(ω，𝒞) → 𝒴，满足𝐴(𝔤.𝑥) = 𝐴(𝑥)适用于所有𝔤 ∈ 𝔊和𝑥∈𝒳(ω、𝒞).

几何深度学习(GDL)蓝图(图形输入)。图片来自 GDL 原型书，第 3.5 节。

把这些块拼在一起，我们就能创造出一个丰富的 *𝔊-invariant 函数 f:*𝒳(ω，𝒞) → 𝒴这种采取形式的

丰富的 𝔊-Invariant 函数类。 GDL 原型书，第 3.5 节。

其中选择每个块，使得其输出尺寸与下一个块的输入尺寸相匹配。另外，每个块可能具有不同的对称群𝔊.

我们在这篇文章中已经看到，仅仅对称可能不足以打破维度的诅咒。然后我们提出了第二个几何先验，称为尺度分离。基于多尺度结构的假设，我们讨论了尺度分离如何抵消维数灾难。我们以 GDL 蓝图作为几何领域的**框架来结束，它是由两个几何先验构建的。我们将在随后的文章中详细讨论这个 GDL 蓝图，以及它如何应用于各种几何领域。我们将在这个框架中看到我们最喜欢的架构，如 CNN、transformers 和 GNNs。

参考

• GDL 球场，( AMMI ，2021 年夏季)。
• j .布鲁纳第四讲[ 视频 | 幻灯片。
• 米（meter 的缩写））m .布朗斯坦、j .布鲁纳、t .科恩和 p .韦利奇科维奇，几何深度学习:网格、组、图形、测地线和量规 (2021)。

我们感谢汉尼斯·斯特尔克对草案提出的有益意见。

带有谷歌地球引擎和 Greppo 的地理空间应用程序

原文：https://towardsdatascience/geospatial-app-with-google-earth-engine-and-greppo-2c166b373382

如果没有丰富的 JavaScript 经验，使用 Google Earth 引擎是很困难的。Greppo 让您在 Python 中克服了这个问题。

使用 Greppo 和 GEE 的最终 web 应用程序。图片作者。

谷歌地球引擎是数据科学家工具箱中处理地理空间数据的一个神奇工具。然而，用 GEE 代码编辑器构建 web 应用程序需要很高的学习曲线。基于 JavaScript 的应用创建器需要专门从事 Python 开发的数据科学家投入大量时间

Greppo 是弥合这一差距的完美工具。

在这篇博文中，我将使用一个流行的 GEE 用例 DEM(数字高程模型)来构建一个 web 应用程序。我将带你了解 GEE 的基础知识、客户机-服务器模型、API 如何工作以及 GEE 数据模型。在这种背景下，这篇文章将使用 Greppo 创建一个使用 GEE 的 Python 接口的应用程序，并强调 Greppo 的心智模型和易于使用的界面。

注:这里所有代码都是用 Python 写的。它们是从 文档 中移植来的 GEE 的样本 JavaScript 代码。

入门指南

在我们开始之前，你需要访问谷歌地球引擎。按照此处的说明注册并获得访问权限。

以下是关于 Greppo 及其使用方法的快速教程:

接下来，让我们设置 Python 环境来安装依赖项。要理解什么是 Python 环境以及如何设置它，请阅读这篇。将以下包安装到 Python 环境中。

pip install earthengine-api greppo

web-app 的代码将放入 app.py、 中，通过命令行使用命令greppo serve app.py服务并运行 app。

注意:要在命令行中运行greppo命令，需要激活安装了 greppo 的 python 环境。app.py 文件可以被重命名为任何名称，但是在运行命令greppo serve app.py时一定要在这个文件夹中，或者在一个相对的文件夹结构greppo serve /demo/folder/app.py中。

Greppo 的 GitHub 库:https://github/greppo-io/greppo

如有任何问题，请使用“问题”在 GitHub 上联系我们，或在 Discord 频道中点击。

GEE 认证和初始化

为了能够使用谷歌地球引擎，你需要创建一个服务帐户，并获得与该帐户相关的访问密钥文件。这只需要几分钟的时间，但是请确保按照说明正确操作。遵循此处的说明。要使用服务帐户和密钥文件，请使用以下代码进行初始化。

注意:确保将 key-file.json 保存在另一个位置，最好安全地保存在您计算机的根文件夹中，不要提交到公共存储库中。

了解 GEE 的客户机-服务器模型

正如 GEE 的开发者文档所说，Earth Engine 不像你以前用过的任何 GIS 或地理空间工具。GEE 主要是一个云平台，所有的处理都在云中完成，而不是在你的机器上。您将与 GEE 进行的交互仅仅是翻译并发送到 GEE 云平台的指令。为了更好地理解这一点，我们必须通过 GEE 的客户端与服务器及其懒惰计算模型。

客户端与服务器

先说我之前提到的， GEE 主要是一个云平台 。它让你在云中完成所有的处理。那么，你如何访问这个处理功能呢？

这里是earthengine-api库派上用场的地方。Python 包earthengine-api向客户机(也就是您)提供对象，作为在云中传递和处理的服务器对象的代理。

为了更好地理解客户机-服务器模型，让我们以客户机中的一个字符串变量和服务器中的一个字符串变量为例。在客户端创建一个字符串并打印它的类型时，我们得到 Python class str对象来表示一个字符串。如果我们想将一个字符串发送到服务器，以便在云中使用或操作，我们可以使用ee.String将数据包装在一个代理容器中，该容器可以在服务器中读取。更具体地说，ee. objects是一个eeputedObject，它是代理对象的父类。

>> # Client side string
>> client_string = 'I am a Python String object'
>> print(type(client_string))
<class 'str'>>> # Server side string
>> server_string = ee.String('I am proxy ee String object!');
>> print(type(server_string))
<class 'ee.ee_string.String'>

代理对象不包含任何实际数据或处理功能/算法。它们只是服务器(云平台)上对象的句柄，仅仅是传达要在服务器上执行的指令。可以把它想象成一种使用代码与服务器通信的方式，要做到这一点，您需要将数据和指令包装在特定类型的eeputedObject容器中。

当对数据执行循环或使用条件语句时，这种理解变得更加重要。为了执行这些，指令将被发送到服务器来执行它们。要了解这些是如何实现的查看本页了解更多细节。

惰性计算模型(延迟执行)

因此，从上面我们知道earthengine-api包仅仅是向服务器发送指令。那么，死刑是如何以及何时执行的呢？

客户端库earthengine-api将所有指令编译成一个 JSON 对象并发送给服务器。但是，这不会立即执行。执行被推迟直到有对结果的请求。对结果的请求可以是一个print语句，或者是要显示的image对象。

这种随需应变的计算会影响返回给客户端(您是用户)的内容。来自 earthengine-api 的结果是一个指向要获取数据的 GEE tile 服务器的 url。因此，在提到的感兴趣区域内的图像被选择性地处理。感兴趣区域由客户端显示中地图的缩放级别和中心位置决定。而且，当您移动和缩放时，图像会被处理并发送到客户端进行查看。因此，图像是延迟计算的。

使用 Greppo 和 GEE

使用 Greppo 显示和可视化地球引擎图像对象是相当直接的，你需要使用的是:app.ee_layer()。在 GEE 中存储地理空间数据的基本数据类型是，

• Image:地球引擎中的基础栅格数据类型。
• ImageCollection:图像的堆叠或时间序列。
• Geometry:地球引擎中的基本矢量数据类型。
• Feature:带属性的Geometry。
• FeatureCollection:一套Feature。

在理解了 GEE 的客户机-服务器和惰性计算模型之后，我们可以推断，这些数据类型是根据对其可视化的请求按需处理的。

那么，如何配合 GEE 使用 Greppo 呢？

最好用一个例子来解释。先说 app 的脚手架。你必须首先从Greppo导入app对象，因为这将是你与前端通信的入口点。然后你将不得不import ee，向地球引擎认证你自己，并用你的上述服务账户的凭证初始化你的会话。

接下来，让我们从从目录中选择数据集开始。这里，我们使用USGS/SRTMGL1_003来获取数字高程图。我们需要首先为 DEM 图像数据中所有大于 0 的值获取一个土地掩膜，为此我们使用dem.get(0)。接下来，我们需要在 DEM 上应用蒙版，只显示陆地，为此我们使用dem.updateMask(dem.gt(0))，并将结果指定为我们要显示的ee_dem。由于所有数据都存储为 int 16(32767 和-32768 之间的值的矩阵)，我们必须使用调色板来显示矩阵。

要添加调色板，我们创建一个可视化参数对象，其中包含生成 RGM 或灰度图像的指令。这里我们使用包含Hex values:[‘006633’, ‘E5FFCC’, ‘662A00’, ‘D8D8D8’, ‘F5F5F5’]的调色板，并将其线性映射到与指定的min -> #006633和max -> #F5F5F5相对应的值。

注意:存储在 DEM 中的数据是栅格，表示为矩阵，每个像元包含表示像元的点的高程(米)。

然后使用 Greppo 在 web 应用程序中可视化该地图，您只需使用app.ee_layer()。ee_object是地球引擎图像对象，vis_param是可视化参数字典，name对应于将在 web-app 前端使用的唯一标识符，description是可选的，为 app 用户提供额外的指导。关于这一点的更多内容可以在文档 这里 中找到。

上述步骤中的 web-app 视图。图片作者。

端到端通信:完整的网络应用

到目前为止，我们已经看到了如何只在 Greppo 可视化地球引擎对象。然而，Greppo 能够在前端和后端之间进行复杂的交互。让我们用一个例子来寻找用户指定的一个点的高程。我们将使用 Greppo 的三个 API 特性。

• app.display():在前端显示文本或减价。
• app.number():前端的数字输入功能，供用户输入数值。它在后端绑定到的变量将被更新为用户指定的值。
• app.text():前端的文本输入功能，供用户输入数值。它在后端绑定到的变量将被更新为用户指定的值。

更多详情请参考 文档 。

让我们从使用app.display ( name是惟一的标识符，值是显示的文本，可以是多行字符串)显示一些文本来指导 web 应用程序用户开始。之后，让我们使用app.number()创建两个数字输入，分别代表该点的经度和纬度。

app.number()接受 name、显示在前端的标识符和 value，value 是这个元素的默认值。接下来，我们也创建一个文本输入，使用app.text()和name和value获得点的名称，如前面提到的app.number()一样。

使用该点的纬度和经度，我们现在可以用可视化参数color: ‘red’为该点创建一个地球引擎几何图形对象。我们现在可以使用上面提到的app.ee_layer()来显示。

为了找到该点的高程，我们在 DEM 对象上使用地球引擎方法sample。我们对 DEM 中的点进行采样，以从 DEM 中获取属性。我们从输出中取出第一个点，并使用.get方法找到与属性elevation相关联的值。最后，我们编写了一个多行字符串来显示输出。

注意:要将地图居中到一个点，并在初始加载时缩放，使用app.map(center=[lat, lon], zoom=level)。

具有交互功能的网络应用视图。图片作者。

结论

我们的目标是使用 google earth engine 的数据和计算功能以及 Greppo 的 web 应用程序开发库，完全用 Python 创建 web 应用程序。我们了解了 GEE 的工作原理，了解了如何将 Greppo 与 GEE 集成。学会使用app.ee_layer()、app.display()、app.number()和app.text()创建一个完整的 web 应用程序，与前端和后端进行端到端的通信。

所有的演示文件都可以在这里找到:【https://github/greppo-io/greppo-demo/tree/main/ee-demo

查看一下 GitHub 资源库:此处为 Greppo 上最新更新。如果您的用例有错误、问题或功能请求，请联系 Discord channel 或在 GitHub 上提出问题。用 Greppo 造了什么东西？贴 GitHub。

• GitHub 库:https://github/greppo-io/greppo
• 文献:https://docs.greppo.io/
• 网址:https://greppo.io/

注:文章原载 此处 。

熊猫专家的地理空间数据争论

原文：https://towardsdatascience/geospatial-data-wrangling-for-pandas-experts-96c130c78bd8

熟悉 Python 知识的地理空间数据分析

安德鲁·斯图特斯曼在 Unsplash 上的照片

这个笑话在我高中时代(很多年前)就流传了。一个男孩只记住了一篇考试作文，题目是“牛”。然而，在考试中，他失望地发现作文题目变成了“河流”。不知道该做什么，他想出了一个绝妙的主意。他写了这样一篇文章——从前有一条河。。然后迅速切换到 ……有一头奶牛坐在河边。那是一头中年的黑白条纹奶牛，长着长长的尾巴。男孩继续写文章，就像他熟悉的领域“牛”那样，最后又回到“河”上。

我们很快就会谈到这个故事的寓意。

地理空间数据产品同时具有丰富的信息和华丽的外观。如果你只是给某人看一张地图，什么都不说也不写，它仍然传递了一个故事。然而，对于数据科学家来说，学习地理空间数据分析的前景可能是可怕的。至少发生在我身上。我没有受过这方面的训练，我做过的最令人兴奋的“地理空间”工作是一张我学习地点的地图——用 Microsoft Paint 创建的。我总是被其他人的地理空间工作所吸引，尽管从来没有想过我自己会尝试这样做。我没有时间投入大量精力从头开始学习另一种工具。第二个障碍是我必须购买专有的 GIS 软件许可证(我还不知道 QGIS，它是免费的)。

当我了解到地理空间数据可以表示为 dataframe 对象时，事情很快发生了变化。当我了解到这一点时，我知道我不必从头开始，我可以在 Python 基础上构建我的地理空间能力。

想法很简单:1)使用合适的 Python 库如geopandas、GDAL将地理空间数据导入到您的笔记本环境中；2)然后将其转换成一个pandas dataframe 对象；3)继续分析和操作pandas中的数据；4)最后，使用matplotlib可视化地图。

地理空间数据有多种形式，如多边形、线、点和栅格，这种方法适用于所有这些形式。今天我将讨论多边形，这样你就可以知道如何使用其他的了。作为参考，下面是不同形式的空间数据的可视化表示:

左图表示多边形(湖)、线段(河)和点(井位)。右图代表光栅图像(图像来源:维基百科)

您可以将多边形视为一个州的县/区边界。同样，河流可以表示为线段；和所有的杂货店一样分。另一方面，在栅格数据集中，区域被划分为正方形而不是多边形，每个正方形包含与特定位置相关联的值/变量/要素(例如，气温、人口密度)。

好了，让我们深入研究 dataframe 对象中的地理空间数据(多边形)。

载入库和数据

您只需要两个库就可以开始了。geopandas用于数据争论，而matplotlib用于数据可视化。顾名思义，geopandas将pandas功能的能力用于地理空间数据。

您可以使用您喜欢的软件包管理器(pip、conda、conda-forge)安装 geopandas:

pip install geopandas

安装完成后，让我们导入这个库。

import geopandas as gpd

该库带有内置数据集，因此您可以立即开始使用。稍后您可以随意试验自己的数据，但是现在，让我们使用内置数据集。我们现在将加载数据集，它包含世界上每个国家的多边形。

# load in dataset
dataSource = gpd.datasets.get_path('naturalearth_lowres')
gdf = gpd.read_file(dataSource)

我们现在将检查刚刚创建的对象的数据类型:

type(gdf)

>> geopandas.geodataframe.GeoDataFrame

这是一个地理数据框架，我们很快会看到它只是一个常规的数据框架，但有一个额外的“几何”列。

你可以用matplotlib的本地命令.plot()快速可视化多边形

gdf.plot()

可视化地理数据框架-世界各国的多边形。x 轴和 Y 轴值分别代表经度和纬度(图片由作者生成)

使用 pandas 功能进行数据操作

在上面，我们已经可视化了地理空间数据，其中每个多边形都是一个国家。

每个面(国家)都带有一些以地理数据框架格式存储的属性。这意味着您可以立即开始使用pandas 功能。让我们来看看数据框的前几行:

gdf.head()

(图片由作者生成)

这就是地理数据框架的样子。它就像一个常规的数据帧，但有一个特殊的“几何”列存储地理空间信息(这个几何列有助于绘制多边形)。

通过像对待数据帧一样对待这个表，您现在可以应用许多pandas功能。让我们尝试一些我们通常在数据科学项目中作为探索性数据分析的一部分使用的熟悉方法:

# getting information about the data
gdf.info()

的输出。应用于地理数据框(由作者生成的图像)的 info()方法

使用上面的.info()方法，我们得到了看起来熟悉的输出。它显示有 177 行(每一行代表一个国家)和 6 列(即每个国家的属性)。我们可以用熊猫.shape进一步证实这一点。

# number of rows and columns
gdf.shape

>> (177, 6)

现在让我们再次使用pandas方法，通过调用unique()方法来检查数据集中有多少个洲。

# unique values of a columns 
gdf['continent'].unique()

>>array(['Oceania', 'Africa', 'North America', 'Asia', 'South America',
'Europe', 'Seven seas (open ocean)', 'Antarctica'], dtype=object)

您也可以对行进行条件过滤。让我们只选择位于非洲大陆的国家。

# filtering rows
gdf[gdf['continent']=='Africa'].head()

(图片由作者生成)

不出所料，您还可以操作列，比如创建新的计算字段。让我们基于两个现有的列创建一个名为GDP _ per _ head的新列: gdp_md_est 和 pop_est。

# create calculated column
gdf['gdp_per_capita'] = gdf['gdp_md_est']/gdf['pop_est']

gdf.head()

(图片由作者生成)

我们现在为数据集中的每个国家增加了一个属性列。

这些只是数据操作的几个例子，您可以尝试一些您感兴趣的其他例子。除了这些数据操作技术之外，您还可以生成汇总统计数据，进行高级统计分析等等。让我们生成一些汇总统计数据:

# generate summary statistics
gdf.describe().T

(图片由作者生成)

总结这一部分，首先，我们使用geopandas库导入地理空间数据(多边形，或者用更专业的术语来说是“shapefile”)，然后使用pandas功能来操作和分析地理数据框架。在下一节中，我们将使用另一个熟悉的 Python 库— matplotlib来研究可视化数据。

形象化

地理空间数据的真正威力在于它能够可视化地理数据框架中包含的不同属性。与数据操作的pandas类似，我们将使用matplotlib来可视化地图中的属性。让我们从一个基本的开始——仅仅可视化形状。

# visualizing the polygons
gdf.plot()

用 Python 绘制多边形。plot()方法。它使用几何列来可视化多边形(由作者生成的图像)

上面的地图显示了多边形。在引擎盖下，这些多边形是从数据框的几何尺寸栏创建的。然而，它还没有显示任何数据，但是我们可以通过指定我们感兴趣的数据列来实现:

# visualize a data column
gdf.plot(column = 'pop_est')

显示每个国家估计人口的世界地图(图片由作者生成)

这张地图现在变得有趣且信息丰富，它用颜色渐变显示了世界上每个国家的估计人口。

但是如果你只想放大非洲呢？这很简单，只需在数据框中过滤非洲大陆，然后以类似的方式创建绘图。

# filter Africa data from the dataframe
africa = gdf[gdf['continent']=='Africa']

# plot
africa.plot(column = 'pop_est')

显示非洲国家估计人口的地图(图片由作者生成)

您还可以访问额外的matplotlib功能来定制地图——例如，删除 x 和 y 轴，在右侧添加图形标题和颜色条。让我们做所有这些。

import matplotlib.pyplot as plot

# use matplotlib functionalities to customize maps
africa.plot(column='pop_est', legend=True)
plt.axis('off')
plt.title("Population in the continent of Africa");

使用 pandas(用于过滤)和 matplotlib(用于绘图)的组合来可视化地理空间数据。(图片由作者生成)

这就是了。您刚刚使用两个库:pandas和matplotlib从地理空间数据中创建了一幅漂亮的地图，这完全符合您的 Python 知识。

这只是起点，从这里开始，天空才是极限！

前进

从这里去哪里？您可以尝试我们刚刚使用的库的一些其他功能，创建更多的地图并进一步定制它们。然后尝试使用外部数据集，网上有很多可用的(如果你需要指点，请在评论中告诉我)。接下来，运用你的想象力和创造性思维，提出一个问题，试着回答。举个例子，你可以问——人口密度排名前 10 的国家是哪些？你应该能够使用pandas方法过滤这些信息，然后使用matplotlib在地图上显示它们。

我在 Python 环境中实现了这个练习，但是如果您对其他编程语言感兴趣，它将以类似的方式工作。如果你喜欢 R 编程，试着用tidyverse和ggplot重现我们在这里所做的。

感谢阅读。请随意订阅以获得我即将在媒体上发表的文章的通知，或者通过 LinkedIn 或 Twitter 与我联系。下次见！

地理空间索引 101

原文：https://towardsdatascience/geospatial-index-101-df2c011da04b

如何有效管理地理空间数据

介绍

由于移动设备、物联网和地球观测技术的进步，如今位置数据变得更加可用。大量的数据对开发定位系统和应用程序的人提出了挑战。本文将简要介绍管理地理空间数据最常用的技术地理空间索引。

什么是索引？

索引是一种优化数据检索和搜索的数据管理方法。使用索引的好处是可以在不扫描所有现有数据的情况下定位特定的数据。每个人都熟悉的索引的真实例子是图书馆。

在图书馆里，书是按照一些标准如书名来分组和排序的。所有以 A 开头的书都在第一书架上，而以 B 开头的书在第二书架上，依此类推。当人们去图书馆时，他们只是去与他们找到的书名相对应的书架，而不是浏览所有可用的书籍，这将花费大量的时间。

照片由龙之介菊野在 Unsplash 上拍摄

索引对于当今所有可用的数据库技术都是至关重要的。有了索引，数据库可以加快查询速度，因为要处理的数据量大大减少了，这类似于图书馆帮助用户找到他们想要的书。

地理空间数据

地理空间数据是记录其地理位置的数据。可以把它想象成一个普通的数据库记录，它有一个附加的列来指定这个数据关联到地球上的什么地方。以下是地理空间数据的一个示例:

地理空间数据示例(图片由作者提供)。

除了看起来像数学函数的位置之外，表上的所有记录看起来都像正常的结构化数据。这种数据类型称为几何，其中用于存储地理数据。下面是 4 个最常见的几何对象:

• 点 :地球上单个位置
• 线 :连接的一组点
• 多边形 :表示特定区域边界的一组点
• 多多边形 :具有相同属性的多边形的集合

(注:T 这里有两种地理空间数据:矢量和栅格。点、线和面是表示位置或区域的矢量数据，而栅格数据是特定场景的图像，其中还包含其位置信息，如卫星和无人机图像。本文只关注矢量数据。)

地理空间索引

与可以应用传统索引算法 的原始数据类型(例如整数)不同，索引地理空间数据要复杂得多，因为它有两个维度: x 和 y，更不用说来自一组点的复杂性，例如线和多边形。因此，它需要不同的索引方法。

一种可能的方法是将整个空间分成一组子空间，如下所示。每个子空间都有一个指定的名称，以便我们以后可以正确地引用它。

图片由作者提供。

假设您有 12 条数据记录，在地理上看起来如下:

图片由作者提供。

通过将数据与预定义的区域相结合，这就是您所得到的结果。

图片由作者提供。

现在，如果您按记录所属的区域对记录进行分组，结果如下表所示:

图片作者。

假设您想要检索所有记录，例如，与澳大利亚相关联的记录。第一种方法是查找与国家边界重叠的所有像元。在这种情况下， L 和 P 。

图片由作者提供。

然后，查找位于两个单元格中的所有记录，得到 X ₄、 X ₅和 X ₁₂.

图片由作者提供。

最后，浏览剩余的记录，我们发现只有 X ₄和 X ₅符合搜索条件。

图片由作者提供。

请注意，尽管总共有 12 条记录，但我们只扫描了其中的 3 条。这个例子中减少的运算数量听起来可能很小，但实际上，对于大量的数据，它可以产生巨大的影响。

地理哈希

上面的例子是对 Geohash 的简化解释。Geohash 是一种算法，它将任何地理位置编码为文本，方法是将世界分成大小相等的单元，称为单元(或网格 ) ，其中每个单元都分配有特定的哈希(但是，它使用的不是字母，而是以 32 为基数的数字)。上一个关于 geohash 的例子中没有提到的是它的层次结构。

假设我们像前面的例子一样将整个空间分成网格。

图片由作者提供。

我们不再停留在这里，而是进一步将每个细胞分割成更小的子空间。例如，在下图中，单元格 A 被分割成 16 个更小的单元格。

图片由作者提供。

然后，我们通过在 A 的一个子节点上重复它，比如说 AF ，更深入一层。

图片由作者提供。

应用此过程越多，每个单元格代表的区域就越小，哈希字符串就越长。名称的长度表示哈希的精度级别。

geohash 的一个有趣特性是，分配给每个像元的名称都有其父名称作为前缀。例如，细胞 AFA 到 AFP 都是 AF 的子细胞。

Geohash 的层次结构(图片由作者提供)。

当您想要验证一个区域是否是另一个区域的子区域时，prefix 属性非常方便。它可以被视为关系数据库中的ST_Contains ( spatial contain )函数的近似版本，具有更低的计算成本，因为它只处理文本，而不处理点、线或多边形。

您可能还会注意到，每个单元格的名称从左到右按字母顺序排序。这种机制允许我们使用小区名称来估计它们之间的地理距离。例如，假设有 3 个哈希: AA 、 AB 、 AK 。由于按字母顺序 AA 比 AK 更靠近 AB (字母‘A’比‘K’更靠近‘B’)，所以 AA 在地理上也比 AK 更靠近 AB 。尽管该属性在新行开始的地方有一个例外(例如 AE 比 AC 更接近于 AA ，但是所有单元格的间距相等的事实使得这种边缘情况易于管理。

注意 G eohash 只是可用的地理空间索引之一。其他的索引系统，如 四叉树 和 H3 ，即使基于相同的基础，应用不同的划分策略，所以也有不同的优缺点。

地理空间索引使用

除了作为索引，地理空间索引还可以应用于各种用例。

数据整理

Geospatial data is very different from primitive data like Integer , Float , and String therefore relational databases need a special data type for them: Geography . The geography data type, however, is only a string in the format of Well-Known Text (WKT) that looks like the following:

POINT(103.14411713327522 15.553626958655387)
LINESTRING(103.14261509622688 15.5550326423843,103.14400984491462 15.552738065460904,103.14585520471687 15.552944785150965,103.14594103540534 15.554040396042701)
POLYGON((103.1434304877674 15.553296208147678,103.14435316766853 15.5547639094724,103.14585520471687 15.553482255373645,103.1448896094715 15.552138577185888,103.1434304877674 15.553296208147678))

Notice that the text is very long due to the decimal precision of the coordinate reference system. On numerous occasions, nevertheless, this precision is unnecessary. For instance, in a ride-sharing application where millions of drivers pick up and drop off millions of passengers a day, it would be impossible(and absurd) to analyze all the points these activities occur. Aggregating the data per zone is much more sensible.

Using an index to aggregate geospatial data(Image by author).

A geohash can be an identifier of each zone that can be used to not only reference the area but also infer its spatial properties. Moreover, since there are finite numbers of cells, using geohash can keep the amount of data in the system constant and more manageable.

Data Partitioning

The vast amount of data available nowadays makes the trend of using distributed data stores grows rapidly. Unlike traditional databases, distributed databases split data into chunks and spread them across nodes in the cluster. Designing a data partitioning strategy is necessary (or at l east recommended) for the system owners to utilize distributed data stores effectively.

In developing geographic-oriented systems and applications, partitioning data by its geographical location is very sensible. The data associated with a similar neighborhood should also be in the same chunk and, therefore, stay in the same node in the database cluster. This way, when users retrieve the data, it will come from a minimal number of nodes and files which saves the computational cost significantly.

Geospatial indexes provide a convenient way to split your data geographically. As you can see from the above examples, the closer the coordinates are, the higher chance they will be in the same geohash(or any other index system) cell. Therefore, you can use geohash as a partition key for the distributed database of your choice.

Using geospatial index for data partitioning(Image by author).

Conclusion

Managing geolocation data is critical for developing location-based applications. A geospatial index provides a systematic and elegant approach to inde x geospatial data. Therefore, it’s a key component that any robust geographic software system must have.

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