如何优化 BigQuery 的使用和成本
使用谷歌数据仓库时你应该知道的最佳实践
在 Unsplash 上由 Delphine Ducaruge 拍摄的照片
BigQuery 是一种 SaaS (软件即服务)数据仓库技术，将按照付费模式向您收费。在这个小故事中，我想告诉你价格是多少，以及如何优化你的使用。
哪些行动需要成本？
首先是好的方面 BigQuery 中的一些操作是免费的，比如:
- 加载数据(例如从存储桶)
- 提取数据(例如，提取到 Google 工作表)
- 复制数据(例如，复制到新表中)
- 删除表、视图和数据集
- 集群和分区表
但是你大概可以猜到，最常见的活动如:
- 数据存储
- 查询数据
确实花了一些钱。特别是最后一个功能是谷歌如何赚钱的。更多详情请访问谷歌大查询定价【1】。
这些行动的成本是多少？
通用概述将让您初步了解 BigQuery 中最常见任务的成本:
- 查询:按处理的字节收费，每 TB 5 美元(每月第一个 TB 免费)
- 存储:对于 BigQuery 存储 API，每月存储的 GB 数将为 0.02 美元(活动存储)和 1.10 美元/TB
- 流式插入将花费您每 GB 0.05 美元
这些是谷歌显示的美国地区的当前价格，其他地区的价格有所不同。
BigQuery 与其他技术相比如何？
将价格与亚马逊的 Redshift 或 Snowflake 等其他技术进行比较很难，因为这取决于用例以及您如何优化每项技术。同样，价格也非常相似 —这里 fivetran 发布了一个超级好的基准【2】。在这个基准测试中，他们得出了相同的结论。然而，重要的是要强调 BigQuery 提供了最高的 SaaS 感受——因此您不必太担心供应和操作，这同样对成本有积极的影响。根据我自己的经验，我可以肯定 BigQuery 不仅易于使用，而且运行起来比传统的数据仓库更便宜。
最后，您如何优化您的使用？
- 使用 web UI 中的验证器来帮助计算成本——UI 和 CLI 将返回查询中处理的字节数。
已处理的预测字节数—按作者分类的图像
2.使用 SuperQuery——使用 super query(例如可以通过 Chrome 安装),你甚至可以直接计算每次查询的费用。此外，该附件还提供了许多其他有用的功能[3]。
超级查询价格预测-作者图片
3.仅包含您需要的列和行 —尤其是不必要的列会产生额外的成本，这很容易避免。
4.尽可能使用缓存结果和永久表来代替视图。
5.监控您的使用情况 —借助 GCP 的日志和监控服务，您可以创建仪表盘和警报。在标准监控中，您已经有了一个很好的概述。
GCP 标准监控—用户提供的图像
但是要了解更多细节，例如哪个用户查询了哪些数据以及查询的成本，您可以使用 Stackdriver。
6.使用 Limit 仍然会处理所有行——所以要注意那个陷阱。
7.删除重复和不必要的数据 —人们经常犯保存不必要的数据的错误(我以后可能需要这个)。这应该通过严格的数据治理来避免，这样就不会出现数据沼泽。
8.使用集群 —它可以提高某些类型查询的性能，例如使用筛选子句的查询和聚合数据的查询。
9.使用分区表 —表将被分成段，称为分区，这使得管理和查询数据更加容易。此外，它会节省你的钱。
结论
我希望这篇文章在成本和成本优化的最佳实践方面有所帮助。我本人是 BigQuery 的忠实粉丝，并在专业和私人场合使用它。它有很多优点，速度超级快。但是，我也知道，考虑一些事情来对抗不必要的成本是很重要的。尤其是在企业家的意义上，人们真正谈论的是大数据。在这里，即使是小错误和反模式也会很快导致大的成本劣势。如果你也有任何建议，请随时告诉我。
资料来源和进一步阅读
[1]谷歌，定价 (2021)
[2] fivetran， 2020 数据仓库基准测试 (2020)
【3】super Query，写出完美的查询。每次都是。 (2021)
大查询获取+多重处理
多重处理是否提高了 BigQuery API 请求的抓取速度？
作者图片
BigQuery 存储读取 API 可用于从 BigQuery 表中获取数据。然而，在我写这些代码的时候，还没有为了更快地处理读取而使用这个 API 结合多处理进行基准测试。
在本文中，我将展示我所做的一些研究和基准测试，以便找到从 BigQuery 获取数据的最佳方法。
抓取+多重处理
从 BigQuery 获取数据的最常见和最简单的方法是只使用机器/实例上的一个内核来线性处理获取。
然而，如果您正在使用具有多个内核的计算机或 GCP 计算引擎(GCE ),您可能想知道使用它们来并行处理更多数据是否有用，但是有一个关键概念需要理解:
内核越多，读取速度不一定越快。从互联网上获取一些数据所需的时间很大程度上取决于路由器/网络上可用的互联网带宽。
因此，在您的设备上拥有 200 个内核并不意味着您的获取处理速度会快 200 倍，而是可以创建 200 个进程，每个进程将获取一个表块(使用较少的带宽)，以便最终将它们合并到一个表/数据帧中。
参考 GCE 文档，每个 vCPU 可以处理 2GB 的入口，最多 6 个 vCPU。所以我构建了一个工具:**bqfetch**，它给出了要使用的核的数量，可以获取一个或多个核的表。如果nb_cores =-1，那么算法将使用机器上可用的虚拟 CPU 的数量，它将为每个内核创建一个进程，并并行处理它们。该工具适用于multiprocessing、billiard和joblib平行后端。
为了证明 BigQuery 存储在多处理中是否有效，我使用 1 到 6 个内核在 2GB 的表上独立运行了 6 次读取，并测量了所用的时间。这是使用该工具的实际代码:
L 油墨至 bqfetch repo 。
中篇关于 bqfetch 快取库。
注意，为了获取巨大的表(多个 TB ),我们必须将整个表分块/分割成可以放入内存的更小的块。因此，第一步是使用chunks()方法划分表格。我们必须指定要用作索引的列，以便将该列中的所有不同值分成多个块。我们还必须指定我们想要处理的每个块的大小，因此这个值必须小于机器上的可用内存(如果是这样，将会出现一条警告消息)。详细模式打印内存信息以及区块的数量/大小。
在中，下面的例子将整个表分成大小为 2GB 的小块和连续使用 1 到 6 个内核获取它们，仅用于演示目的。注意，在这种情况下，块的大小和表的大小是 2GB，所以只有一个块将被
chunks()方法返回，但是在实际使用情况中，您将指定一个更大的表，并且该方法将返回一个块的列表以进行顺序处理。
在 12 个虚拟内核的 MacBook Pro 上的结果
作者图片:获取 2GB 表作为 pandas df 的时间，取决于核心/进程(PC)的数量
如图所示，在使用标准 Wi-Fi 连接的个人电脑上使用多个内核并不能提高性能，甚至会更差。这是可以预料的，因为网络上可用的互联网带宽保持不变，即使我们使用多个内核来获取数据。由于每个进程必须创建一个新的 ReadSession 以从 BigQuery 读取数据，这需要更多的时间，因为我们必须并行初始化多个连接，并且我们还必须在每个进程完成时合并结果。
在 30 个虚拟内核的 GCP 计算引擎上的结果
作者图片:获取 2GB 表作为 pandas df 的时间，取决于核心/进程的数量(GCE)
使用 GCE 时，会出现同样的问题，当内核数量增加时，性能仍然会更差。根据 GCE 文档，我们应该增加网络带宽，那么问题是什么呢？发生这种情况的原因对我来说很模糊，但我猜这是由于 BigQuery 限制关于每个请求的行过滤器的大小只有 10MB。
结论
- 在我写这几行文字的时候(2021 年 8 月)，Google 还没有提出任何原生实现来高效地使用多处理 BigQuery 存储。因此，由于 API 的限制，我们实现这一目标的唯一方法并没有像我们希望的那样优化。我们必须等待 API 贡献者的进一步改进，以便有可能使用这种方法获得性能。
- 您可以使用我已经实现的工具通过 BigQuery API 使用多处理来做您自己的测试。代码是开源的，所以欢迎你来查看或投稿。
用于排列测试的 BigQuery 存储过程
在 Google BigQuery SQL 中运行置换测试的简单方法
存储的过程允许我们执行打包成一个“函数”的多个Google biqquery SQL操作。了解如何使用存储过程快速有效地对任何数据集应用置换测试。
只需要密码？从这里开始别忘了启动它！⭐️
向您的 BigQuery SQL pantry 添加另一罐技巧|照片由 Luisa Brimble 在 Unsplash 上发布
我还有几个 BigQuery 教程:
要获得所有媒体文章的完整访问权限，包括我的文章，请考虑在此订阅。
什么是存储过程？
不，存储过程不存储在软盘上|照片由 Fredy Jacob 在 Unsplash 上拍摄
在 BigQuery 中，可以编写自定义，用户自定义函数(UDF)。这些函数非常棒，因为您可以减少代码中的重复，使代码对其他人更具可读性。然而，你也可以在 BQ 中编写脚本。这包括循环、if 语句、变量声明等。不允许在函数中使用它，因为函数是作用于表元素的，而不是作用于脚本中的变量。
这就是存储过程的用武之地。把它们看作是函数的等价物，但不是脚本的等价物。有了它们，您可以重复所有复杂的 SQL 逻辑，即使它们包括变量和 for 循环的赋值。
这篇来自 Google Cloud 的文章介绍了存储过程的概念，如果你需要复习 BigQuery 脚本，这里有关于如何编写循环和条件的官方文档。
谈到存储过程，有几个问题。最大的问题是它们可以接受一个或多个输入，但是输出必须是声明的变量。例如，要运行文档中的示例，我们必须首先声明一个适当的变量来存储结果:
这个例子摘自官方公告。
如果你不知道上面的程序在做什么，也不用担心。我们稍后会在这里做一个例子。
一旦您用CREATE PROCEDURE命令创建了您的过程，它将像 BigQuery 中的任何其他表一样显示，并且您将能够用通常的dataset.procedure()符号调用它。
让我们看看这一切是如何工作的。
带有存储过程的斐波那契数
斐波那契无处不在|照片由Yair MeJA在 Unsplash 上拍摄
我总是发现，如果例子中的逻辑很容易理解，那么例子是最好的。为此，我将创建一个存储过程，当调用时，计算第 n 个斐波纳契数。
作为复习，斐波纳契数列是一个自然数序列，其中第 n 个元素由F(n)表示，并由F(n) = F(n-1) + F(n-2)和 F(0) = 0以及F(1) = 1定义。所以你得到:0，1，1，2，3，5，8，13…
为了做到这一点，我们将首先使用脚本来完成，然后将其包装在一个过程中。我们需要做的是:
- 初始化我们的变量——我们需要记录我们已经计算了多少个数字，以及何时完成计算
- 做一个 for 循环，其中每一步我们更新我们当前跟踪的 2 个数字
- 选择我们的答案并打印到控制台
我在上面的脚本中添加了一堆注释，所以要确保你明白是怎么回事。
我总是建议先写脚本，然后把它改成存储过程。您会发现这样调试容易多了。
这些关键字执行以下操作:
DECLARE创建一个具有特定类型的新变量，并将其设置为默认值SET将已经声明的变量更新为新值WHILE <cond> DO和END WHILE缠绕我们的 while 循环。SELECT是你的普通 SQL 语句打印到控制台的答案
如果我们运行上面的代码，我们会将13打印到控制台，这是第 7 个斐波那契数列的正确答案。通过更改代码的第一行并将i设置为其他值来尝试其他数字。
现在剩下要做的就是将上述内容包装在一个过程中，并将其存储在一个数据集中。我们是这样做的:
变化不大，但有几件重要的事情需要注意:
- 我们通过
CREATE PROCEDURE声明我们的过程,后面是数据集和我们想要存储它的过程名。确保您已经创建了数据集*ds*！ - 我们必须声明我们的过程输入—
i INT64在过程的顶部。 - 我们对它的输出做同样的事情:
OUT answer INT64。 - 输入和输出不再需要在过程内部声明，所以我们
SET将答案设为 0，而不是使用默认值。 - 我们将整个脚本包装在
BEGINEND语句中。
如果您运行上面的代码，您应该会得到类似于:该语句创建了一个名为 your-project-name.ds.fibonacci 的新过程。
为了使用这个过程，我们需要声明一个整数变量和CALL过程来填充它。
如果你都做对了，你应该得到:55 (8+13 =21，13+21 = 34，21+34=55)。🎉
仅此而已。您现在知道了如何使用 BigQuery 脚本编写存储过程。给自己一个读到这里的掌声吧！
我想，即使是岩石也赞同 BigQuery 过程… | 来源
什么是排列测试？
由卡洛斯·鲁伊斯·华曼在 Unsplash 上拍摄的照片
我听到你在说什么了。"这些都很酷，但为什么有用呢？”。好吧，如果我告诉你，你可以用 SQL 在一些巨大的数据集上运行一些超级的统计测试，而且很便宜，会怎么样呢？这正是我们在这里要做的😉！
我知道我不会很好地解释什么是排列测试，所以看看 Jared Wilber 的这个很棒的网站,在那里你会读到一个羊驼洗发水的实验。
简而言之，这就是排列测试:
- 无效假设是两组的统计数据相同。
- 你用统计学来衡量差异(在我们的例子中是指)。
- 然后，假设零假设是真的，你在 2 组之间排列观察。这里的逻辑是，如果零假设确实是真的，那么你如何给你的实验对象贴上标签并不重要。
- 然后你数一数有多少次你观察到排列的统计差异比正确标记的设置更极端。
- 这一比率将给出您的 p 值，即我们由于随机运气观察到这样或更极端值的概率。
- 如果 p 值很低(例如< 5%) then you reject the null hypothesis.
But as I said, Jared 的网站做得更好，所以请去看看。
使用 BigQuery SQL 进行排列测试
使用 BigQuery 的强大功能，我们可以快速有效地调整数据集|图片来自Amol ty agionUnsplash
我再次建议，在用程序把事情复杂化之前，先弄清楚你的脚本的逻辑。在我们做任何事情之前，我们需要设定一些要求。我们需要一个至少有两列的表:一列是我们要取平均值的值，另一列是我们的组标签的 0 和 1。我们还假设我们的样本大小相等。这不是问题，因为我们将使用 BigQuery 中 iris 数据集的子集，其中 Virginica 和 Versicolor 物种都有 50 个观察值。
我们准备好实施我们的计划了。这是我们要做的:
- 将我们的表格列重命名为 measurement and treatment(0 和 1 ),以使脚本的其余部分更加一致。
- 添加另一个
row_id列以确保我们能够以一种简单的方式对进行采样。我在这里解释了我们如何以及为什么要这样做。 - 用整数列表进行交叉连接——每个置换样本一个。
- 记录原始表格的观察差异。
- 记录我们排列的平均差异。
- 计算我们得到的差高于或等于我们原始均值差的次数的比率。
下面是 SQL 的样子:
以上只是普通的 SQL，我一直建议你以这种方式开始。当我们想把它变成一个程序时，事情就变得复杂了…
我们希望我们的过程接受下面的输入:
- 我们要对其应用测试的表的名称。
- 包含测量值的列的名称。
- 包含我们的治疗标志的列的名称。
- 我们想要进行的排列的数量。
据我所知，我们不能只是将表名粘贴到脚本中，所以我们将在 BigQuery 中使用可爱且非常强大的EXECUTE IMMEDIATE命令来完成这项工作。这个家伙基本上运行一个字符串作为 SQL 语句，允许我们动态地编写 SQL 字符串，然后执行它。相当整洁！
是的，立即执行就是那种强大的… | 来源
由于我们的查询很长，我决定使用CONCAT函数将多行字符串组合成一行。这对于调试来说非常烦人，所以我建议您运行您的函数的第一遍，将合成的 SQL 复制到一个新窗口中，并使用 BigQuery UI 在那里调试您的查询——相信我，您会看到许多遗漏的逗号和空格弄乱了您的查询😉，我当然做了。
第 15、16 和 19 行是奇迹发生的地方。
这比我们之前的查询更难理解，但基本上是一样的。我只是在列名和表名中使用了字符串串联||到粘贴以及我们在前 2 个表中要求的排列数。其余的不变——因此有了标准化的列名😃。
运行以上程序来存储您的过程！然后我们准备准备我们的桌子。我们将使用经典的免费 Iris 数据集，它可以作为 BigQuery 表在bigquery-public-data:ml_datasets.iris下获得。我们将扔掉setosa品种，因为我们只需要 2 个，我们将询问海滨锦鸡儿花是否有更大的萼片宽度。如果你需要复习虹膜数据集，我推荐这款 Kaggle 笔记本。
我特意将所有的专栏都留在那里，以便您也可以对其他专栏尝试这个过程。
我们有桌子和程序。是时候创建一个浮点变量并调用我们的排列测试了:
运行上面的程序，我得到的 p 值是 0.0008499，非常小。因此，我们得出结论，弗吉尼亚种的萼片宽度明显大于云芝种。
BigQuery 的美妙之处在于，我们为查询的字节付费，而不是为 CPU 时间付费。因此，我们可以将排列增加到 100 万，但仍然只为 1.56KBs 的查询付费。
🙏请在一些较大的数据集上尝试上述方法并评论一下它完成查询的速度，我很想知道你的结果。
结论
照片由 Massimo Sartirana 在 Unsplash 上拍摄
希望您会发现以上内容对您的数据科学工作流程有用。现在，您不必将大型 BigQuery 数据集导出到磁盘，租用半 TB 的内存机器并在其上运行 Python 来做一些统计测试，您可以在 BigQuery 中便宜而有效地完成所有这些工作。
如果没有别的，请记住以下几点:
- 一个存储过程是一个可重用的 BigQuery 脚本(带有循环和变量)，它接受输入并产生一个输出变量。
- 程序以
CREATE PROCEDURE dataset.procedure_name(inputs..., output) BEGIN开始，以END结束。 - 不要
DECLARE输入或输出变量。使用SET改变变量值。 EXECUTE IMMEDIATE可能是最强大的脚本关键字，允许您将任何字符串作为 SQL 语句执行。- 用
CONCAT和||将表名和变量名粘贴到 SQL 字符串中。 - 总是首先在过程之外调试 SQL 脚本！
感谢一路读到最后。我写了一些关于 BigQuery、Julia 和 Python 的文章，所以如果这些教程对你有帮助的话，请随时关注我。
想要更多 BigQuery 采样教程吗？看看这个。
用于细粒度视觉识别和多模态深度学习的双线性池
巴德·赫利松在 Unsplash 上的照片
高级神经网络架构通过学习特征交互来工作
双线性池起源于计算机视觉社区，作为一种细粒度视觉识别的方法。或者用不太花哨的语言来说，一种在识别和分类视觉对象时寻找特定细节的方法。在高层次上，该方法的工作方式如下。给定一个输入图像 I ，我们将 I 送入两个不同的深度卷积神经网络 A 和 B，见图 1。在应用若干池化和非线性变换之后，我们从 A 和 b 输出特征图。这两个网络可以被预训练以解决不同的任务。直觉是，A 和 B 以这种方式从输入图像中学习不同的特征。例如，A 被训练来检测基本的物体形状，而 B 检测纹理特征。然后，来自 A 和 B 的输出特征由所谓的双线性池层组合。这仅仅意味着我们通过取它们的内积，将 A 的每个特征与 B 的每个特征结合起来。读者可能会注意到，这类似于支持向量机中的 2 次多项式核。双线性池层背后的直觉是，这些特征交互允许我们检测图像的更多特定细节。
图一。双线性池的整体结构。来源[1]
原始建筑
我假设读者熟悉卷积神经网络的基础知识。不然我推荐这些入门讲义。
图 1 中网络 A 的输出是一个维数为 n×d 的矩阵 U，或 n×d*，*，网络 B 的输出是一个维数为 m×d 的矩阵 V，在 CNN 的上下文中，n 和 m 分别是 A 和 B 的输出层中滤波器(或内核)的数量。每个滤波器的维数是 *d，*它是通过展平二维特征图获得的，即经过几次核卷积和汇集变换的输出图像。双线性池操作被定义为
等式 1。双线性池。
然后，矩阵 Z 被展平成维数为 m*n 的向量 Z。Z 和可学习向量 w_k 的内积被提供作为 softmax 激活函数的输入，其中 K 个类别中的类别 K 的概率被计算为:
等式 2。使用双线性映射进行分类。
在上面，w_k 是维数为 m*n 的可学习向量，它模拟向量 z 属于第 k 类的概率，1≤k≤K。
双线性池为什么有效？
如前所述，双线性池的核心是明确考虑特征交互的思想，类似于支持向量机中的多项式核。让我们考虑一个具体的例子。作为第一个应用，双线性池用于细粒度的图像分类。考虑图 2，我们想解决一个鸟类分类问题。显然，这些都是彼此相似的鸟类，我们需要具体的细节来正确区分不同的物种。对于图 1 中的两个神经网络 A 和 B，假设 A=B，并对 A 和 B 使用预训练的 VGG16 模型，VGG16 是一个深度 CNN 模型，已经在 ImageNet 2014 挑战赛的大量图像上进行了训练，其中图像属于 1，000 个类别中的一个。很明显，VGG16 可以检测鸟类羽毛、喙、颜色等形状。然而，许多鸟类都有相同的喙形。
图二。即使对一个人类专家来说，正确分类所有的鸟类也是一个挑战。来源[2]
VGG16 足以正确区分鸟类和猫，但要正确区分不同的鸟类是有挑战性的，差异已经非常微妙了。为了区分不同的物种，我们需要明确地考虑这些特征的组合。这样的组合对每个物种来说更有可能是独一无二的。羽毛颜色和鸟喙形状的结合更有可能唯一地识别一种鸟。因此，使用正确识别某些基本形状的预训练模型，我们可以了解这些形状的哪些组合对于手头的分类问题是重要的。
用于多模态学习的双线性池
双线性池的另一个主要应用是结合不同领域的神经网络模型。例如，在视觉问答中，训练数据由自然语言和图像中的问题组成，参见图 3 中的示例。
图 3。集合视觉和语言特征进行视觉问答。来源[3]
双线性池的变体
有几种不同的方法可以使用双线性池。
- **单个神经网络 vs 两个网络。**网络 A 和 B 可以是相同的网络，即 A=B。如在计算机视觉示例中所讨论的，我们然后从输出神经网络的不同通道学习重要的特征对。对于多模态问题，当我们组合来自不同领域的特征时，网络 A 和 B 当然彼此不同。
- 端到端架构与预训练模型。两个网络 A 和 B 可以是像 VGG16 那样的具有固定权重的预训练网络。在这个设置中，新的输入通过 A 和 B 传递，然后我们只学习双线性层的权重。在端到端架构中，来自双线性层的误差反向传播到 A 和 B，并且更新它们各自的参数。双线性层是可微分的，它可以由众所周知的深度学习框架来处理，如 TensorFlow 或 PyTorch。
计算的复杂性
显然，双线性池的一个主要缺点是它的计算复杂性。矩阵乘法
时间复杂度为 O(mnd)。通常，t 尺寸 d 相当小。想象一个 CNN，我们从 256x256 像素的图像开始逐渐粗化维度，更深的层具有 16x16 维度的核图。另一方面，为了获得强大的模型，核的数量需要更大。高级架构使用几千数量级的 m 和 n。对于具有 K 个类别的分类问题(见等式 2)，这意味着我们需要 Kmn 个可学习的参数，其数量级可以是几百万。这可能使双线性池难以训练，并容易过度拟合。
可伸缩双线性池的矩阵分解
已经设计了几种旨在降低原始算法的计算复杂度的方法。我将介绍一种矩阵分解方法。
设 W_k 是一个 m 乘 n 矩阵，具有可学习的双线性权重。方程 2 中的向量 w_k 是平坦矩阵 W_k，网络 A 的输出是 m-次-d 维的矩阵 U，网络 B 的输出是 n-次-d 维的矩阵，设第 I 列 U 为 u_i，V 的第 I 列为 v_i，特殊情况 d=1 时，矩阵 U 和 V 只是单列。
z_i 是 u_i 和 v_i 的平坦的外积，我们可以证明以下等式成立:
上面简单地写出了等式两边的量。
根据定义，它成立
因此，我们可以将内积重写为
为了提高计算复杂度，我们可以用低秩分解来代替双线性矩阵 W_k:
在上面，r 是一个小常数，使得 r < min(m，n)。
因此，我们获得内积的以下表达式:
注意，右边求和的计算复杂度是 d*(mr + nr)。因为 r 是一个小常数，比 m 和 n 小得多，所以我们实现了更好的运行时间。此外，双线性池图层的可训练参数的数量从 Kmn 减少到 K*(m+n)r。为了更好地了解这两个量的比较情况，请考虑一个 CNN，其典型值为 m=n=1，000 个过滤器和 K=10 个类。使用 r = 20，我们从 10⁷到 410⁵减少了可学习参数的数量，即因子分解双线性池中的参数数量仅为原始设置的参数数量的 4%。
双线性池的缺点和限制
最后，我们应该提到双线性池的一些限制:
- 对于两个输出，特征向量 d 的维数必须相同。对于两个 CNN 模型来说，这不是一个大问题，但是对于多模态双线性池来说，这可能是一个限制。例如，单词嵌入可以具有与 CNN 非常不同的最佳维度，调整这两者可能会导致信息的丢失。
- 我们仅限于考虑成对的相互作用。研究人员提出了像三线性池这样的一般化方法，其中我们考虑了三个特征。但这使得模型更难训练，容易过度拟合。
履行
在 FashionMNIST 数据集上用于图像分类的双线性池的简单实现可以在这个 Jupyter 笔记本中找到。请注意，双线性池并不真正适合像 FashionMNIST 这样的小型数据集。该实现作为一个示例，可以根据项目的需要进行调整和改进。
参考
[1]宗-林玉，阿鲁尼·罗伊·乔杜里，苏博兰苏·马吉:
用于细粒度视觉识别的双线性卷积神经网络。IEEE Trans。肛门模式。马赫。智能。40(6)，可用此处
[2]舒孔，查尔斯·c·福尔克斯。用于细粒度分类的低秩双线性池。CVPR 2017。可用此处，项目页面此处
[3]周瑜，，，陶大成。用于视觉问答的多模态分解双线性池联合注意学习。ICCV 2017，此处可用
台球与弗雷德的常客和拜耳的贝叶斯[拜耳赢:)
我在这里讨论贝叶斯自己的一个实验，它强调了贝叶斯方法的美丽、优雅和正确性，即使频率主义方法给出了不正确的结果。
在 1763 年死后发表的论文[1]中，Bayes 展示了一个例子，其中 Frequentist 点估计给出了不正确的结果(与实验不匹配)，而 Bayes 分析给出了与实验匹配的优雅解决方案。本文中的问题确切问题陈述是文献[1][2][3]中的陈述的推广。
在这里，我讨论了这个例子，并用 colab 笔记本[4]链接了推导和蒙特卡罗模拟。
设置
不可否认，这个设置有点做作，但我们可以忽略它。让我们在宽度为 1 的桌子上打一场台球。卡罗尔进来，随机地打了一个台球，当它落在离一面墙 p 的距离后，她把它藏在了窗帘后面。然后爱丽丝和鲍勃，不知道球的位置，射击他们自己的球，如果一个球落在卡罗尔的球的左边，爱丽丝赢了这一轮(所以概率是 p)，如果它在卡罗尔的球的右边结束，那么鲍勃赢了这一轮(所以概率是(1-p))。先得 6 分的人赢得比赛。
图片由 Klara Kulikova 提供
一些见不得人的角色正在观看这场游戏，当鲍勃在 K 杆时赢得了 nB 回合时，他们叫来了常客弗雷德来赌这场游戏。他们请常客弗雷德预测鲍勃赢得比赛的概率。在给了弗雷德频繁主义者多次尝试解决问题，他悲惨地失败了，这些人物听到贝叶斯统计，并邀请你，拜尔贝叶斯显示你的把戏包里有什么。
常客弗雷德如何尝试解决方案
常客 Fred 首先通过查看 Alice 获胜的分数来估计 p 的点数
注意 nA+nB=K。
在这之后，Frequentist 简单地计算 Bob 在 11-K 剩余游戏中赢得至少 6-nB 的概率(游戏数量永远不能超过 11。说服你自己这个论点。)
常客弗雷德对鲍勃获胜概率的回答
我们以后会看到他错得有多离谱。
拜耳贝叶斯如何拯救世界
在这一点上，你被叫来尝试你的贝叶斯手在这个问题上。你采取了一种非常不同的方法。你知道这一点(如果你不知道这一点，那么你应该考虑一下和/或阅读一下条件概率)
这里发生的是，我们对 p 的所有可能值进行积分，因为我们不知道它是什么。模型中的这种参数称为讨厌参数。
现在既然 Carol 均匀采样 P P§= 1。出于可读性的考虑，我在这里不做数学计算，推导过程可以在我下面引用的 github 笔记本中找到。
最后的表述是
鲍勃获胜概率的贝叶斯答案
实验
哪个答案是正确的？这是一个非常容易模拟的问题，已经在下面的笔记本中完成了。让我展示几个例子。
首先注意，如果 K=0，那么 nB =0，虽然直觉上答案是 0.5，Bob 会赢，但常客 Fred 说这是未定义的！
在 K=1，2，3，4，5 的另一个极端，如果 K=nb，Frequentist 的点估计给出 p=0 和 PF(B)=1，所以他肯定 Bob 会赢，但稍加思考就会发现这是不正确的。
甚至连常客弗雷德也认为两个答案都不正确，但他不知道该怎么办，而是继续坚持认为其他答案是正确的。
然而，弗雷德不仅在上面两个例子中是错误的，他还经常错误(双关语！).下面比较不同的 K 和 nb，贝叶斯答案如何总是匹配 MC(人们看到的小的不匹配来自 MC 采样中的有限大小效应),并且频率主义者的答案总是错误的，除了当 nb=K/2 时。
K 和 nb 的其他值的结果可以在上面引用的笔记本[4]中找到。
帮助常客弗雷德认识他的错误[高级]
那么，常客弗雷德做错了什么？为什么他经常得到错误的答案？让我们再来看看 Bayer the Bayesian 是如何处理讨厌的参数 P 的，它可以是从 P§中提取的任何值(在这种情况下恰好是均匀的)，然后对结果进行积分。另一方面，frequent ist Fred 取 p=(K-nB)/K，这是 P(D|p)的最大似然估计量，并将其插入 P(B|p，D)。现在，这将在可能性非常尖锐的情况下工作(这就是我们物理学家所说的狄拉克δ函数)。然而，事实证明情况并非如此。同行的物理学家将会看到这里与路径积分有很大的相似之处，在量子力学中，对所有路径求和，但在经典极限中，只有一条路径对路径积分有贡献。用这个类比，频率主义者弗雷德在一个问题中使用了经典极限，他应该做量子力学。
更快的方法
在经历了所有的困难之后，让我指出，有一种更简单的方法可以实现这一点，那就是我们所做的是对二项数据的后验预测分布求和，即 Bob 的相关获胜次数的贝塔二项分布。
鲍勃以文明的方式获胜的概率的贝叶斯答案。
结论
贝叶斯方法更直观，因此更难出错。虽然常客们显然也能得到正确的答案，但要做到这一点却困难得多。虽然关于贝叶斯和频率主义者的哲学辩论将持续很长时间，但在实践层面上，贝叶斯是你最好的选择。
参考资料:
- T.贝叶斯。一篇解决机会主义中一个问题的论文。伦敦皇家学会哲学汇刊 53(0):370–418，1763
- 艾迪先生。什么是贝叶斯统计？。自然生物技术 22:1177–1178，2004
- 关于使用 Python 的贝叶斯统计的非常好的介绍。我自己的例子是对这个例子的概括，最初的引用来自贝叶斯https://arxiv.org/pdf/1411.5018.pdf
- 我的笔记本https://github . com/borundev/DNN _ 讲座/blob/master/Bayesian _ pool . ipynb
二值化神经网络综述
二值化神经网络是一种有趣的神经网络变体，可以节省内存、时间和能量。了解它们是如何工作的。
亚历山大·辛恩在 Unsplash 上的照片
使用神经网络的一个障碍是运行网络所需的功率、内存和时间。这对于没有强大 CPU 或大内存的移动和物联网设备来说是个问题。二值化神经网络是解决这一问题的方法。通过使用二进制值而不是浮点值，可以更快地计算网络，并且使用更少的存储器和功率。
力学:
从概念上讲，二进制神经网络(BNN)类似于常规前馈神经网络(NN)。一个区别是 BNN 中的权重和激活被限制为只有两个值:1 和-1，因此得名“二进制化”。因此，只使用了一个激活函数:符号函数。我们不能使用像 sigmoid 或 relu 这样的常规激活函数，因为它们有连续的输出，正如我们刚才所说的，激活只能是 1 和-1。
来源:整个神经网络的布尔掩蔽
当我们使用梯度下降来训练网络时，这些限制会导致一些问题。首先，符号函数的梯度为 0。这是不好的，因为它使所有权重和激活的梯度也为 0，这意味着实际上不会发生任何训练。我们可以通过在进行反向传播时忽略符号激活函数来解决这个问题。然而，我们有另一个问题，即梯度更新将导致权重不再是 1 或-1。我们通过保留一组实值权重并对这些权重进行梯度更新来解决这个问题。然后，网络权重是这些实值权重的二进制化。
因此，我们的最终算法是:
- 前向传递:我们有实值权重 W_r 和输入向量 x，首先，我们对所有权重应用 sign 函数，得到 W_b = sign(W_r)。然后我们像往常一样用 W_b 和符号激活函数计算神经网络的输出。
- 反向传递:我们像往常一样进行反向传播，并计算权重 W_b 的梯度，只是我们忽略了符号激活函数。我们通过减去这些梯度来更新权重 W_r(记住，W_b 来自 W_r，所以不直接更新 W_b！).
来源:[ 2 ]
性能:
我们已经看到了二进制神经网络是如何工作的。现在让我们比较一下普通神经网络和 BNN 之间的理论性能差异。普通的神经网络使用 32 位浮点数。如果我们有两个 32 位寄存器，我们可以用一条计算机指令在两个 32 位数之间执行一次乘法。二进制网络使用 1 位数(我们将+1 编码为 1，将-1 编码为 0)。要将两个 1 位数字与我们的表示相乘，我们可以使用 XNOR 指令。
来源:[ 2
我们能做的是将 32 个 1 位数字放入我们的两个 32 位寄存器中的每一个，然后运行一个 XNOR 指令来同时相乘所有的数字。其结果是，我们可以用一条指令做 32 次 1 位乘法，而用一条指令只能做 1 次 32 位乘法。因此，与普通神经网络相比，BNN 理论上有 32 倍的加速比。
实际上，32x 这个数字很难达到。CPU/GPU 的编程方式会影响指令调度，指令本身并不总是需要相同数量的时钟周期来执行。此外，现代的 CPU/GPU 并没有针对运行按位代码进行优化，因此必须注意代码的编写方式。最后，虽然乘法是神经网络中总计算的很大一部分，但也有我们没有考虑的累加/求和。尽管如此，[ 1 ]报告称，在比较优化的 BNN 和未优化的正常神经网络时，加速比提高了 23 倍，这表明明显的加速是完全可以实现的。
因为我们使用的是 1 位数字而不是 32 位数字，我们还希望使用大约 32 倍的内存。当我们使用的内存减少 32 倍时，我们的内存访问也减少了 32 倍。这也将减少大约 32 倍的功耗。我找不到测试这个假设的实验数据，但这个结论对我来说似乎是合理的。
最后，需要注意的是，所有这些 BNN 优势仅适用于运行时，不适用于训练时。这是因为，正如我们前面提到的，我们保留了一组实值权重来进行训练。这意味着梯度是一个 32 位浮点数，并不受制于我们所描述的 1 位优势。因此，BNN 模型不应该在内存/功率受限的设备上训练，但是它们可以在这些设备上运行。
精度:
当然，如果网络由于测试集上的不良准确性而不可用，那么无论速度/内存/功率提高多少都没有关系。二值化激活和权重比常规神经网络的表达能力差，所以我们预计准确性会更差。问题是会恶化到什么程度。我找不到任何回答这个问题的理论论文，所以现在我认为像下面这样的经验数据是我们最好的。我们在这里描述的 BNN 实现是第一部分，常规 NN 变体是最后一部分。
来源:[ 1
我们看到，误差肯定会增加，特别是在 CIFAR-10 数据集上(10.15%对 7.62%)。然而，考虑到 BNNs 如前所述的所有性能优势，这大约 3%的差异并不是世界末日。对于每个精度点都至关重要的任务，例如医疗 x 射线筛查，我们不想使用 BNNs。但是在精度不那么重要或者性能更重要的情况下，bnn 是一个可行的选择。
在这篇文章中，我们已经简要了解了 BNN 的机制。我们还看到了 BNNs 如何显著提高速度、内存使用和功耗，代价是精度稍低。要更深入地了解 BNNs，请看下面的参考资料。
[1] M. Courbariaux 和 I. Hubara，二值化神经网络:训练深度神经网络，权重和激活限制为+1 或-1 (2016)
[2] T. Simons 和 D. Lee，二值化神经网络综述 (2019)
基于自动机器学习的二元分类
使用开源的 MLJAR auto-ML 更快地构建精确的模型
自动化机器学习工具的兴起让开发者能够更快地建立准确的机器学习模型。这些工具通过执行特征工程、算法选择、调整以及记录模型来减少工程师的工作。一个这样的库是开源的 MLJAR 包。在本文中，让我们看看如何使用这个包进行二进制分类。
入门指南
MLJAR 包可以用来建立一个完整的机器学习管道，具有特征工程和超参数调整功能。该软件包还支持流行的机器学习算法，包括:
MLJAR 根据选择的算法建立几个模型，并通过集合或堆叠模型获得最终预测。一旦训练完成，该软件包还提供可视化，可用于解释模型。其中包括功能重要性以及显示各种型号性能的排行榜。
为了开始使用这个包，你需要通过pip安装它。
pip install mljar-supervised
获取数据
对于这个示例，让我们使用 Scikit-learn 来创建一个分类数据集。这可以使用数据集模块中的make_classification函数来完成。
下一步是将这个数据集分成训练集和测试集。
使用 MLJAR 进行二元分类
令人惊讶的是，使用 MLJAR 进行二进制分类只需要几行代码。MLJAR 在幕后负责所有的机器学习魔法。这里的第一步是导入AutoML类。
下一步是创建该类的一个实例，同时指定您希望用于解决该问题的算法。可以设置的其他参数包括:
- 模式 —该套件配有四种内置模式。
Explain模式是解释和理解数据的理想模式。它产生了特征重要性的可视化以及树的可视化。在为生产建立 ML 模型时使用Perform。Compete是为了建立用于机器学习竞赛的模型。Optuna模式用于搜索高度调谐的 ML 模型。 - 算法 —指定您想要使用的算法。它们通常作为列表传入。
- 结果路径— 存储结果的路径
- total_time_limit —训练模型的总时间(秒)
- train_ensemble — 指示是否在训练过程结束时创建一个集合
- stack_models —决定是否创建模型堆栈
- eval_metric —将被优化的度量。如果
auto``logloss用于分类问题，而rmse用于回归问题
下一步是使算法适合训练和测试集。
MLJAR 解释
模型拟合后，您可以使用report功能获得各种模型解释。例如，您可以看到符合此二元分类问题的模型的排行榜。
当您单击单个模型时，您还会看到几个可视化效果和模型摘要。
例如，你会看到混淆矩阵。
MLJAR 还会默认显示训练过程的学习曲线。为验证集和训练集绘制学习曲线。垂直线显示最佳迭代次数，用于进行预测。
你也可以从报告中发现特性的重要性。这是排列重要性的例子。使用 Sckit-learn 的排列重要性进行计算。
MLJAR 还绘制了 SHAP 重要性。
它还显示了 SHAP 重要依赖图。
姆加魔法
当用 MLJAR 训练一个模型时，在幕后会发生一些事情。机器学习模型的训练分几个步骤进行。这些步骤需要:
- 训练简单的算法，如决策树，以获得快速的洞察力。
- 使用默认超参数训练模型。这里每个算法都有一个模型。
- 在给定的超参数上执行随机搜索。
- 从原始数据中构建新的特征，使模型具有更强的预测能力。MLJAR 将这些特性称为黄金特性。这些特性可以在结果路径中的一个名为
golden_features.json的文件中找到。这些特征被添加到数据中，并且使用 Xgboost、CatBoost 和 LightGBM 算法来训练新的模型。 - 为每个算法获得的最佳模型被选择，并且其超参数被用于基于一组选择的特征来训练新的模型。
- 模型被进一步调整，然后被集合。
为什么要使用 MLJAR
人们为什么会考虑在机器学习工作流中使用 MLJAR 有几个原因:
- 它执行特征预处理，例如缺失值的插补和类别的转换
- 执行高级特征工程
- 能够转换时间和日期
- 执行超参数调谐
- 提供广泛的模型解释和可视化
- 为您运行的每个实验创建降价报告
最后的想法
在本文中，您已经看到了如何构建二进制分类模型 MLJAR auto-ML。您已经看到这个过程只需要几个步骤。您还看到了 MLJAR auto-ML 使用多种算法为分类问题构建多种模型。下一步是在您自己的二元分类数据集上使用该包。
https://github.com/mljar/mljar-supervised/
二元分类:通过非穷举网格搜索和交叉验证的 XGBoost 超参数调整场景
平衡模型性能和计算资源限制的实际例子——代码和可视化
简介
XGBoost 或极端梯度提升是当今使用最广泛的机器学习算法之一。众所周知，它在赢得卡格尔比赛方面效率很高。许多文章称赞它，并指出它优于替代算法的优势，因此它是练习机器学习的必备技能。
虽然 XGBoost 相对较快，但在标准笔记本电脑上运行脚本仍然具有挑战性:当拟合机器学习模型时，它通常会附带超参数调整和交叉验证(尽管不是必需的)。如果考虑到超参数的所有组合，调整可能会指数地增加计算资源需求。因此，穷尽搜索最佳参数可能不是一个可行的选择。
本文介绍了超参数调优的替代方法，旨在找到模型性能和计算资源限制之间的良好平衡。我自豪地拥有一个英特尔酷睿 i5–8250 u CPU，根据 CPUBenchmark 的数据，它的分数低于提交的中位数。对于有这些限制的人来说，我的代码可能会很有趣。
技术术语概述
尽管我通过分享该过程的代码从实用的角度出发，但我想对所使用的术语和参数名称做一个概述。如果您熟悉它们，请跳过这一点，进入数据集描述和代码结果。
自相关
决策树学习目标变量可以取一组离散值的树模型称为 分类树*；在这些树结构中，树叶代表类别标签，树枝代表导致这些类别标签的特征的合取*
梯度推进
梯度推进是一种用于回归、分类和其他任务的机器学习技术，它以弱预测模型的集合的形式产生预测模型，通常是决策树。当决策树是弱学习器时，产生的算法称为梯度提升树，它通常优于随机森林。它像其他 boosting 方法一样以分阶段的方式构建模型，并通过允许优化任意可微分损失函数来推广它们【3】注:XGBoost 通过其正则化机制来防止过度拟合，从而区别于其他梯度 boosting 技术。
交叉验证
" 交叉验证通过“将数据样本划分为互补的子集，对一个子集(称为训练集)执行分析，并对另一个子集(称为验证集或测试集)的分析进行验证”来组合(平均)预测中的适合性度量，以获得更准确的模型预测性能估计值……使用不同的分区执行多轮交叉验证，并在各轮中组合(例如平均)验证结果，以给出模型预测性能的估计值
网格搜索
网格搜索是一个过程，它通过手动指定的目标算法超参数空间的子集进行彻底搜索…并基于生成的超参数集评估成本函数【5】
随机搜索
" 随机搜索…使用概率分布为每个超参数独立选择一个值…并基于生成的超参数集评估成本函数【5】
贝叶斯搜索
Eta (学习率，[0，1]，默认值= 0.3)
步长收缩用于更新，以防止过拟合。【7】
“类比于 GBM 中的学习率”[8]
Gamma ([0，∞]，默认值=0 ):
" 在树的一个叶节点上做进一步划分所需的最小损失减少量。【7】
仅当结果分裂给出损失函数的正减少时，节点才被分裂【8】
最大深度 ([0，∞]，默认= 6):
一棵树的最大深度。增加这个值会使模型更加复杂，并且更有可能过度拟合
λ(未提供界限，缺省值= 1):
关于权重的 L2 正则化项。【7】
……类似于岭回归【8】
Alpha (未提供界限，缺省值=0):
" L1 正则化项关于权重。【7】
……类似于套索回归【8】
树方法(auto、exact、approx、hist、gpu_hist、默认= auto ):
、exact——精确贪婪算法。列举所有分割候选
approx — 近似贪婪算法使用分位数草图
hist — 更快的直方图优化近似贪婪算法
gpu_hist — GPU 实现 hist 算法【7】
增长策略 ([ depthwise，lossguide，默认= depthwise ):
" 控制新节点添加到树中的方式。
depthwise — 在最接近根的节点处分割。
lossguide — 在损耗变化最大的节点处分割【7】
准确性
"…正确预测的观测值与总观测值的比率 " [9]
精度
" 精度是正确预测的正观测值与总预测正观测值的比率 " [9]
召回
召回是正确预测的正面观察与实际类中所有观察的比率【9】
F1 评分
F1 评分是准确率和召回率的加权平均值。因此，这个分数将假阳性和假阴性都考虑在内…如果假阳性和假阴性具有相似的成本，准确性效果最好。如果假阳性和假阴性的成本相差很大，那么最好同时考虑精确度和召回率[9]
混淆矩阵
"…一个表格，通常用于描述一个分类模型对一组真实值已知的测试数据的性能 " [9]
ROC AUC
ROC曲线(接收器操作特性曲线)是显示分类模型在所有分类阈值的性能的图表。该曲线绘制了两个参数:真阳性率和假阳性率…AUC代表“ROC 曲线下面积”也就是说，AUC 测量从(0，0)到(1，1)的整个 ROC 曲线下的整个二维面积(考虑积分)…AUC 提供了对所有可能的分类阈值的综合性能测量。解释 AUC 的一种方法是将模型对随机正面例子的排序高于随机负面例子的概率*【10】*
问题陈述阐述
机器学习模型带有默认参数:如果你没有为可选参数分配特定的值或字符串，算法会通过预设值或字符串自动完成。然而，最佳参数可能因数据集而异，因此默认参数通常不是最佳的。需要通过网格搜索找到最优参数，其中网格代表每个参数的实验值(n-维空间)。
如上所述，计算能力需求的指数增长问题通过应用蛮力方法和穷尽搜索每个组合而出现。考虑到损失函数，穷举搜索产生全局最小值，因此每次尝试都是有吸引力的。存在替代的、更快的方法——应用的方法在上面列出——以产生对应于局部最小值的一组参数值的概率为代价。
全局和局部损失函数最小值的简化表示(来源:作者)
供考虑的备选方案:
- 逐步搜索又称坐标下降——一次优化一个参数。可能导致局部最小值。在几个超参数的情况下，可能存在许多局部最小值。它取决于超参数的阶数，达到哪个最小值。也可能有这样的情况，当两个或多个超参数的组合修改会产生更好的分数，但修改其中一个不会[11]
- 随机搜索-从搜索空间中随机选择超参数值。搜索空间中的随机点:可能不是最小值，但可能是全局最小值
- 贝叶斯搜索——当函数的数学形式未知或计算成本很高时，这是一种极其强大的技术。背后的主要思想是根据数据计算目标函数的后验分布(使用著名的贝叶斯定理)，然后选择好的点来尝试这个分布【12】
我的研究目标是提供一种方法，通过检查模型得分和运行时间来比较上述方法，并在资源有限的情况下选择最合适的方法。值得注意的是，每个数据集可能需要使用有限的子样本进行单独调查，然后在考虑最佳结果时，将最合适的方法应用于整个数据集。
数据集
选择用于演示的数据是 CERN 希格斯玻色子挑战的子样本，并已在 Kaggle 上发表为“希格斯玻色子和背景过程”[13]。
CERN ATLAS 团队已经发表了原始出版物，著名的希格斯玻色子已经通过对这些数据的分析得到了证实。
简而言之，该数据包含 ATLAS 测量的 21 个特征(运动学特性)、从低级别特征导出的 7 个(高级别)特征和一个表明该过程是希格斯过程还是背景噪声的结果的二进制特征。
代码和结果
既然一张图抵得上千言万语，请让我来评估一下下面的流程图中的代码，随便用 Paint 做的(开个玩笑，用 Sketchviz 做的)。请在我的 Jovian 简介这里看完整的代码。
流程图(作者用 Sketchviz 制作)
注意:在每种方法的情况下，存储墙壁时间，并且测试用例的数量都相等，以便进行有意义的比较。
导入了几个评估函数以及标准库:
从 CSV 格式加载数据，检查列车数据。id 被丢弃，因为它作为每个过程的唯一值对分类没有贡献:
检查测试数据帧发现，一些列包含不适当的字符，因为它们被格式化为对象。原来一些缺失的值用“？”来表示，因此是对象格式。问号被替换为 0，作为代表缺失值的标准值:
作为一个标准的过程，自相关和偏相关已被评估。
因为自相关不存在于任何特征中，所以可以肯定地说我们不是在处理时间序列数据。因此，随机化数据帧不会造成差异，但是，出于比较目的(与已公布的子样本分析结果)，它在训练/测试分割期间应用，除了在分割(分层)后特别注意保持类别频率比率相同。
在执行网格搜索算法之前，必须拟合一个基准模型。通过调用fit()方法，默认参数被获取并存储以备后用。由于GridSearchCV在列表中接受输入，单个参数值也必须被包装。通过在GridSearchCV实例上调用fit()，执行交叉验证，提取结果，计算分数并存储在字典中。
重要提示:除了搜索最佳网格搜索算法之外，可能还会设置一些速度参数。此类参数为 *tree_method* ，设置为 *hist* ，将连续特征组织在桶(箱)中，读取列车数据变得明显更快【14】。请阅读参考资料，了解 XGBoost 情况下的更多技巧。
迭代整个参数网格需要很多时间，因此将详细度设置为 1 有助于监控该过程。但是，墙壁时间不等于打印的安装时间，因此循环时间也被跟踪和打印。
坐标下降是一种特殊情况，在这种情况下，后续迭代的性能被跟踪，并且可以被可视化以选择最佳参数设置。情节如下:
- 所选超参数网格每次试验的评估分数，以及最佳超参数值的混淆矩阵
- 后续迭代的最佳参数值的评估分数的演变
重要提示:由于交叉验证是在训练数据上进行的，而手动 AUC 分数是使用测试数据计算的，因此测试分数不一定在每次后续迭代中增加。因此，最终的模型参数是手动选取的。
随机搜索和贝叶斯搜索以相似的方式执行，在相同的结果字典中捕获结果。
最后，比较 AUC 分数和搜索时间。所有的方法都做了同样多的工作，最终，贝叶斯搜索胜出。随机和贝叶斯搜索的执行时间都比 CD 短。
重要提示:执行多次搜索，例如嵌入循环中，执行时间和分数可能会有所不同。测试分数差异的显著性超出了我的文章范围。然而，我重新运行笔记本几次，贝叶斯搜索总是获胜，尽管执行时间也有所不同，并在 45 分钟左右达到峰值。
总结
有趣的是，对于所有三种方法，一些超参数值是相似的还是不同的。例如，对于 CD 和贝叶斯搜索，最大树深度被设置在顶部网格值，但是 lambda 参数对于每一个都是完全不同的。在每种情况下，学习率都保持在低水平。当然，在可能的组合中已经进行了 68 次试验(631 800)，但是模型已经得到了改进，同时在预计的 256 478 分钟的强力网格搜索中节省了至少 256 448 分钟。
可以尝试扩展作业的数量，使其与总搜索空间更具可比性，并可能进一步增加 AUC 分数。
值得注意的是，原始论文分别使用超参数调节的 NN(浅层神经网络)和 d NN(深层神经网络)报道了 0.841 和 0.883 的 AUC 分数。表现最差的 CD 算法在看不见的数据上得到了 0.8033/0.7241(AUC/准确度)的分数，而数据集的发布者使用决策树分类器获得了 0.6831 的准确度分数，使用支持向量机(SVM)获得了 0.6429 的准确度分数。考虑到所使用的硬件，这将 XGBoost 算法和结果放在上下文中。
参考文献
[1]
【2】【https://en.wikipedia.org/wiki/Decision_tree_learning】
【3】
【4】https://en . Wikipedia . org/wiki/Cross-validation _(Statistics)
【5】Simon Chan，Philip Treleaven，第 5 章——大规模推荐系统的连续模型选择，主编:Venu Govindaraju，Vijay V. Raghavan，C.R. Rao，Handbook of Statistics，Elsevier，taskId = 3370
【14】https://towards data science . com/do-you-use-xgboost-heres-how-to-make-it-200 x-faster-16cb 6039 a16e
二元神经网络——低成本神经网络的未来？
二元神经网络能代替全精度网络吗？
亚历山大·辛恩在 Unsplash 上的照片
E 每年都有更深入的模型被开发出来执行各种任务，比如对象检测、图像分割等。它总是能打败最先进的模型。但是，人们越来越关注让模型更轻、更高效，以便它们可以在边缘设备和移动设备上运行。这对于弥合机器学习的研究价值和生产价值之间的差距极其重要。
一种减少深度神经网络的存储和计算成本的方法是二进制神经网络的概念。二进制神经网络的概念非常简单，其中权重和激活张量的每个值使用+1 和-1 来表示，使得它们可以以 1 位而不是全精度来存储(-1 在 1 位整数中表示为 0)。浮点值到二进制值的转换使用如下所示的符号函数—
符号阈值函数
现在，使用上述阈值函数的一个主要问题是函数的梯度趋于零。一个解决方案是使用一个直通估计器。直通估计器是这样一种估计器，它在反向传播期间完全按照原样传递梯度，而没有任何变化。这简化了二进制神经网络中阈值函数的反向传播机制，并且显示出非常好的效果。
直通估计器解释(图片由作者提供)
在梯度累积阶段，使用二进制权重和激活来累积每一层的梯度。但是权重更新是在实值原始权重上进行的。为什么会这样呢？我举个例子解释一下。
- 我们有一个神经网络中的节点(原始参数)，值为 0.05。
- 在将该值通过阈值函数时，我们获得值 1 ( 二进制参数)。
- 让我们假设在二进制参数处累积的梯度为 3。我们使用直通估计器的概念将其传递给原始参数。
- 现在，一旦梯度累积，我们就可以更新值。
如果我们使用二进制值来获得新的参数值，那么参数的新值将是 1–0.1 * 3(其中 0.1 是学习率)，即 0.7 。
如果我们使用原始参数值获得新参数，那么新值将是 0.05–0.1 * 3，即 -0.25 。 - 这就是可以注意到差异的地方。当我们执行下一次迭代时，我们将再次通过二进制阈值函数传递新的参数值，这将给出不同的结果。如果我们继续使用参数的二进制值，我们可能永远无法改变该位，因为损失总是从 1 或-1 开始计算。
当参数值被更新时，这些值被裁剪为介于-1 和 1 之间。主要原因是这些值只是增加/减少，否则对网络没有任何影响。需要注意的一点是，最后一个激活层中的值没有被二进制化，而是直接用于分类/回归。这些是二进制神经网络的基本概念。
可以添加到权重和激活层二进制化的另一个改进是使用比例因子来表示权重和激活。这里，比例因子只是权重向量中所有值的平均值的绝对值。如果我们有一个值为[[0.2，0.2]，[-0.6，-0.6]]的 2*2 矩阵，那么比例因子α将为 0.2，二进制矩阵表示为[[1，1]，[-1，-1]]。
重量的α(比例因子)计算
因此，一旦权重和激活以二进制形式表示，以及它们各自的比例因子，它们就可以使用下面的等式来表示。这里，I 和 W 分别代表激活层和权重层，而激活和权重的比例因子分别用κ和α表示。
重量和激活层的 XNOR 操作
现在，作为卷积运算基础的乘累加( MAC )函数是一种非常昂贵的运算。这现在可以用一个 XNOR+popcount 操作来代替。受益？几乎每一个 CPU 都固有地实现了按位操作，并且实现起来明显更快更便宜。popcount 操作只是检查设置的位。下面的例子说明了 MAC 操作可以被 XNOR+popcount 操作所取代。
XNOR + POPCOUNT 运算取代 MAC 运算，其中 n 是原始位数(图片由作者提供)
二进制神经网络的主要缺点是它们不能达到与全精度深层网络一样好的精度。但这种情况一直在慢慢改变，每年(更像是每个月都有论文发表)，随着差距的缩小，都有很多进步。由于人们越来越关注在计算资源有限的设备上实现机器学习模型，未来几年该领域的研究工作将会增加。
参考文献
[1]对 PyTorch 实现的直通估计器的直观解释。(2020 年 9 月 19 日)。检索于 2021 年 1 月 28 日，来自https://www . hassanaskary . com/python/py torch/deep % 20 learning/2020/09/19/intuitive-explain-of-straight-through-estimators . html
[2] Courbariaux，Matthieu 等人，“二值化神经网络:训练深度神经网络，其权重和激活被限制为+ 1 或-1。” arXiv 预印本 arXiv:1602.02830 (2016)。
[3] Rastegari，Mohammad 等，“Xnor-net:使用二进制卷积神经网络的图像网络分类”欧洲计算机视觉会议。施普林格，查姆，2016。
[4] Sush16。(2017 年 10 月 02 日)。理解二元神经网络。检索于 2021 年 1 月 28 日，来自https://sush science . WordPress . com/2017/10/01/understanding-binary-neural-networks/
浮点数的二进制表示形式
照片由米卡鲍梅斯特
你有没有想过计算机如何在内存中存储浮点数，比如3.1415 (𝝿)或9.109 × 10⁻³¹(以千克为单位的电子质量)，而内存是由有限数量的 1 和 0(也就是比特)组成的？
对于整数来说似乎很简单(例如17)。假设我们有 16 位(2 字节)来存储这个数字。在 16 位中，我们可以存储范围为[0, 65535]的整数:
(0000000000000000)₂ = (0)₁₀(0000000000010001)₂ =
(1 × 2⁴) +
(0 × 2³) +
(0 × 2²) +
(0 × 2¹) +
(1 × 2⁰) = (17)₁₀(1111111111111111)₂ =
(1 × 2¹⁵) +
(1 × 2¹⁴) +
(1 × 2¹³) +
(1 × 2¹²) +
(1 × 2¹¹) +
(1 × 2¹⁰) +
(1 × 2⁹) +
(1 × 2⁸) +
(1 × 2⁷) +
(1 × 2⁶) +
(1 × 2⁵) +
(1 × 2⁴) +
(1 × 2³) +
(1 × 2²) +
(1 × 2¹) +
(1 × 2⁰) = (65535)₁₀
如果我们需要一个有符号的整数，我们可以使用二进制补码并将[0, 65535]的范围移向负数。在这种情况下，我们的 16 位将代表[-32768, +32767]范围内的数字。
您可能已经注意到，这种方法不允许您表示像-27.15625这样的数字(小数点后的数字将被忽略)。
虽然我们不是第一个注意到这个问题的人。大约 36 年前，一些聪明人通过引入浮点运算的 IEEE 754 标准克服了这个限制。
IEEE 754 标准描述了使用这些 16 位(或 32 位或 64 位)来存储更大范围的数字的方式(框架)，包括小浮点数(小于 1 且更接近 0)。
为了理解标准背后的思想，我们可以回忆一下科学符号——一种表示过大或过小的数字(通常会产生一长串数字)的方式，以便于用十进制形式书写。
插图由 trekhleb 绘制
正如您在图中看到的，数字表示可能被分成三部分:
- 符号
- 分数(也称为有效数)—数字的有价值的位数(含义、有效载荷)
- 指数-控制分数中小数点移动的距离和方向
我们可以省略基础部分，只要同意它等于什么。在我们的例子中，我们将使用2作为基础。
我们可以共享这些位，同时存储符号、指数和分数，而不是使用所有 16 位(或 32 位或 64 位)来存储数字的分数。根据我们要用来存储数字的位数，我们最终得到以下分割:
使用这种方法，分数的位数已经减少(即，对于 16 位数，它从 16 位减少到 10 位)。这意味着分数现在可能取更窄的值范围(失去一些精度)。然而，由于我们也有一个指数部分，它实际上会增加最终的数字范围，也允许我们描述 0 和 1 之间的数字(如果指数是负的)。
例如，带符号的 32 位整数变量的最大值为 1 = 2，147，483，647，而 IEEE 754 32 位二进制浮点变量的最大值约为 3.4028235 × 1⁰ ⁸.
为了使负指数成为可能，IEEE 754 标准使用了偏置指数。想法很简单——从指数值中减去偏差，使其为负。例如，如果指数有 5 位，它可能从范围[0, 31]中取值(这里所有值都是正的)。但是如果我们从中减去15的值，范围就是[-15, 16]。数字15称为偏差，通过以下公式计算:
exponent_bias = 2 ^ (k−1) − 1k - number of exponent bits
我试图在下图中描述浮点数从二进制格式转换回十进制格式的逻辑。希望它能让你更好地理解 IEEE 754 标准是如何工作的。为了简单起见，这里使用 16 位数字，但是同样的方法也适用于 32 位和 64 位数字。
插图由 trekhleb 绘制
为简单起见，上述示例中省略了几个极限情况(即-0、-∞、+∞和NaN(非数字)值)
以下是不同浮点格式支持的数字范围:
代码示例
在 javascript 算法库中，我添加了一个二进制到十进制转换器的源代码，它在上面的交互示例中使用过。
下面是一个例子，说明如何在 JavaScript 中获得浮点数的二进制表示。JavaScript 是一种相当高级的语言，这个例子可能太冗长，不像低级语言那样简单，但仍然可以在浏览器中直接进行实验:
const singlePrecisionBytesLength = 4; // 32 bits
const doublePrecisionBytesLength = 8; // 64 bits
const bitsInByte = 8;function floatAsBinaryString(floatNumber, byteLength) {
let numberAsBinaryString = ''; const arrayBuffer = new ArrayBuffer(byteLength);
const dataView = new DataView(arrayBuffer); const byteOffset = 0;
const littleEndian = false; if (byteLength === singlePrecisionBytesLength) {
dataView.setFloat32(byteOffset, floatNumber, littleEndian);
} else {
dataView.setFloat64(byteOffset, floatNumber, littleEndian);
} for (let byteIndex = 0; byteIndex < byteLength; byteIndex += 1) {
let bits = dataView.getUint8(byteIndex).toString(2);
if (bits.length < bitsInByte) {
bits = new Array(bitsInByte - bits.length).fill('0').join('') + bits;
}
numberAsBinaryString += bits;
} return numberAsBinaryString;
}function floatAs64BinaryString(floatNumber) {
return floatAsBinaryString(floatNumber, doublePrecisionBytesLength);
}function floatAs32BinaryString(floatNumber) {
return floatAsBinaryString(floatNumber, singlePrecisionBytesLength);
}// Usage example
floatAs32BinaryString(1.875); // -> "00111111111100000000000000000000"
参考
您可能还想查阅以下资源，以更深入地了解浮点数的二进制表示:
Python 中“臭名昭著”的算法:二分搜索法
数据结构和算法
为什么它的性能优于线性迭代？
凯西·马蒂亚斯在 Unsplash 拍摄的照片
介绍
让我们从两个问题开始今天的话题。
你学习算法的第一印象是什么？
到底什么是算法？
如今，有几十个像“比特”、“智能”、“人工智能”、“学习”、“监督学习”和“算法”这样的热门词汇在架子顶上。
但是很少有人理解这些词的真正含义，更不用说它们是如何工作的了。
“这是什么？
它是做什么的？"
你问过了。
对我来说，算法是我们希望计算机如何处理信息的过程，什么先什么后。
就是这样。
就这么简单。
我在两篇文章( Python 和 FAANG )中详细阐述了数据科学面试变得越来越受编程驱动和面向 CS 的主题。在今天的帖子中，我将介绍二分搜索法的基本知识和意识形态，并回答 FAANG 提出的一些真实的面试问题。二分搜索法绝对是经过最严格测试的算法之一，如果不是唯一的话，它非常容易掌握。
二分搜索法:什么和如何
官方将其定义为“在一个排序数组中找到目标值位置的搜索算法”，并将目标值与中间元素进行比较，检查它们是否等价。如果不相等，则移到下一个区间进行比较(改编自 Wiki )。
非官方的，我的直观定义是将搜索范围一分为二，询问计算机所选区间是否包含该值:如果区间小于目标值，则上移搜索空间，重新搜索；否则，向下移动。
一场猜谜游戏即将开始。它叫做“我的号码是多少？”基本上，我从 1 到 100 中挑选一个数字，你的工作就是通过问我问题来找到这个数字。
求数有两种方式:无效率和有效率。第一，你可以问我数字是不是 1；如果不是，检查范围内的其余数字，直到 100。如果我选了 100，你需要 100 次才能选对。这就是所谓的野蛮暴力法，虽然有效，但效率很低。如果数字在 1 到 100 万之间，你会怎么做？
第二，更聪明有效的方法是问这样的问题:
#1 问题:“你的数字高于 50 吗？”
“没有。”
所以，数字必须是从 1 到 49！我们甚至不需要检查任何高于 50 的值。
#2 问题:“你的数字高于 25 吗？”
“是的。”
现在，数字必须在 26 和 49 之间！没有必要检查从 1 到 25 的数字。
两个问题之下，我们缩小了真数的范围！你不会花太长时间最终得到它。
这是玩这个游戏时的算法思维。尽管二分搜索法的飞行名气很大，但它有一种脚踏实地的方法。最终，算法会告诉计算机如何处理信息。谈到算法，这是我的第一条规则。
Python 实现
基础学完了，就该用 Python 写点代码了。在这一部分，我通过科技公司问的三个面试问题来体验代码。每个问题都有几种解决方案，在应用二分搜索之前，我将介绍使用内置函数的方法。
通常对于 Python 编码挑战，你的面试官允许你使用任何内置的库和方法。只要您获得了预期的输出，您就可以开始了。随后，他们会增加额外的限制，比如你不能使用这个库/包，或者提高你的内存和速度。
出于教学目的，我介绍了这两种方法，请在准备技术面试时练习这两种方法。
以下问题按难度从大到小排序 。这是二分搜索法算法的内部笑话。学完今天的内容，你会笑。
问题 1:Sqrt(x ), FAANG
-给定一个非负整数 x，计算并返回 x 的平方根
-由于返回类型为整数，小数位数被截断，只返回结果的整数部分。
——比如 8 的平方根是 2.82842，但是这个函数应该只返回整数部分，2。
-https://leetcode.com/problems/sqrtx/
走过我的思考
该问题要求得到一个数的平方根，如果它不是完美的平方根，则只返回整数部分。两步之内就能完成。
1.求 x 的平方根；
2.只返回整数部分。
下面是 Python 代码。
2
numpy floor()和 int()方法的作用相同，它们只返回整数部分，而忽略小数部分(如果有的话)。但是它们的输出是不同的数据类型:floor()得到一个浮点数，int()返回一个整数。
你可能同意，这很容易被包括在 FAANG 的技术面试中。最有可能的是，你的面试官要求你不要用内置的方法，应用一些算法。
在这里，二分搜索法派上了用场。二分搜索法算法有三个步骤。在我们开始之前，请确保数组已经排序，或者减少或者增加。
#第一步。定义搜索空间:左、右和中间
#第二步。我们大胆猜测，从中间开始搜索算法。
#第三步。使用“if-else”语句检查目标值和搜索空间中的每个元素。
3
前两步不言而喻，让我们将第三步分解成不同的场景。第一个 if 条件(步骤 3.1)很简单，如果 root_squared 和 x 相同，它将返回根值。
elif 条件值得解释。如上所述，二分搜索法就像在玩一个从 1 到 100 的猜谜游戏，每次你都可以问一个问题:它是高于还是低于一个数字。条件说如果平方根小于 x，那么根号不够大，我们需要上移，下次再猜一个更大的数。
else 语句规定了当上述两个条件不满足时该做什么。如果平方根大于我们的数字 x，计算机应该向下移动搜索空间，就像我们之前玩的猜谜游戏一样。
所有这些条件都会在 while 循环中自动结束运行。最后，我们返回搜索空间的上界，aka。右。
如此简洁的解决方案！我总是说一个好的算法应该有直观的意义，二分搜索法说的很有道理。
#问题 2:二分搜索法，亚马逊、微软、贝宝和彭博
-给定一个排序的(按升序)n 个元素的整数数组 nums 和一个目标值，写一个函数在 nums 中搜索目标。如果目标存在，则返回其索引；否则，返回-1。https://leetcode.com/problems/binary-search/
走过我的思考
对于这个问题，我的首选方法不是二分搜索法，而是应用一个 enumerate ()方法来迭代链表，如下图所示:
4
结果显示该数字等于其在位置 4 的索引。可以使用列表理解将 for 循环进一步压缩成一行代码。
现在，让我们改变思路，使用二分搜索法。
解决办法
4
这是二分搜索法算法的标准实现，分别检查了 if-elif-else 语句。
到目前为止，您应该对算法的基本形式有了很好的理解。是时候实践该算法更复杂的应用了。
在 Unsplash 上 engin akyurt 拍摄的照片
#问题 3:检查一个数是否是排序数组中的多数元素，作者脸书
-给定一个按非降序排序的数组 nums 和一个数字目标，当且仅当目标是多数元素时，返回 True。
-多数元素是在长度为 N 的数组中出现 N/2 次以上的元素
-https://leet code . com/problems/check-if-A-number-is-majority-element-in-A-sorted-array/
走过我的思考
相比较而言，这是一个稍微复杂一点的问题。我的第一个方法是线性方法，即找到数组中的多数元素，然后判断目标值是否等于多数元素。
如果允许内置函数，我们可以这样做:
Result for the 1st test case is: True
Result for the 2nd test case is: False
简单地说，我们从集合中导入 Counter，遍历数组，如果条件成立，则返回 True。
如果反制方法不可行，事情会变得很糟糕，但我们不会完蛋。
二分搜索法可以扭转乾坤。
像我一样，如果这是你第一次学习算法，可能会有些不知所措。所以，我创建了这个问题的三个版本(难度递增)。一旦我们理解了如何回答前两个问题，最终(原始)版本就变得更“可解”
基本版本 1:使用二分搜索法检查目标值是否在一个排序数组中。
轻松点。这是算法的标准实现。我们已经涵盖了物流，并在问题 1 和 2 中提供了详细的解释。
基本版本 2:使用二分搜索法查找排序数组中某个元素的第一次和最后一次出现。
接下来，让我们增加一些趣味，找到一个元素的第一次出现(最左边)和最后一次出现(最右边)。
找到元素的第一个和最后一个位置很重要，因为它将多数元素定义为长度为 N//2 的元素。如果我们能够确定最左边/最右边的索引，我们就能很好地找到最终的解决方案。
向左搜索第一个匹配项
The Leftmost Occurrence is at Position: 2
有趣的是。在基础版本 2.1 中，搜索空间内的唯一区别是只有两个条件需要检查。
1.当选择的元素小于目标值时，我们需要向上移动左边界。
2.否则，元素大于或等于目标值，我们决定中间为上限。
向右搜索最后一次出现的内容
The Rightmost Occurrence is at Position: 11
类似地，我们可以使用上面的代码找到该元素的最后一次出现。
完整的 Python 代码在我的 Github 上有。
#高级版本 3.1:使用二分搜索法检查目标是否为多数情况
原问题 3
True
哇！
代码看起来很复杂，但坦白说很简单。
初始部分叫做 一个辅助函数 ，和其他任何函数是一回事。唯一的区别是我们将在外部函数中使用 helper 函数提供的功能。
helper 函数( 从第 8 行到第 16 行 )查找元素最左边的位置，这在计算第 24 行和第 25 行中元素的第一次和最后一次出现时很方便。最后，我们检查最右边和最左边位置之间的差值范围是否大于 N//2。
我在我的 Github 上包含了一个 3.1 版本的改进解决方案。
恭喜你的算法之旅。
外卖食品
- 成长心智( 学会如何学习 )。踏入未知的领域是可怕的。有了成长的头脑，如果我们坚持做下去，我们会做得更好。我过去通常害怕 Python 和编码。我担心如果我不知道如何写一个简单的 for 循环，我会显得多么笨拙。所有的负面情绪都在持续，很难集中精力学习。用 Python 写代码，解决有挑战性的谜题，带给我如此多的快乐和满足感。Python 是我的 禅 。
- 算法。这是我们用来告诉计算机如何处理信息的另一种语言。在其花哨的标题下，算法简单明了，就像我们玩的猜谜游戏一样。
- Practice makes perfect. 熟能生巧 (It’s the Chinese equivalence). Rome wasn’t built overnight. All of these popular expressions tell us one thing: no substitute for hard work. Practice daily, and we will be better at coding, Python, algorithm, and everything and anything.
Medium 最近进化出了它的 作家伙伴计划 ，支持像我这样的普通作家。如果你还不是订户，通过下面的链接注册，我会收到一部分会员费。
https://leihua-ye.medium.com/membership
我的数据科学面试序列
</6-python-questions-you-should-practice-before-coding-interviews-f958af55ad13> </5-python-coding-questions-asked-at-faang-59e6cf5ba2a0>
喜欢读这本书吗？
还有，看看我其他关于人工智能和机器学习的帖子。
宁滨机器学习特征工程奖
使用宁滨作为一种技术来快速轻松地创建用于机器学习的新特征。
蒂姆·莫斯霍尔德在 Unsplash 上的照片
如果您已经训练了您的模型，并且仍然认为准确性可以提高，那么可能是时候进行特征工程了。特征工程是使用现有数据创建新特征的实践。这篇文章将关注一种叫做“宁滨”的特征工程技术。
这篇文章假设读者对 Python、Pandas、NumPy 和 matplotlib 有基本的了解。大多数情况下，提供链接是为了更深入地了解正在使用的内容。如果有什么说不通的地方，请留言评论，我会尽力阐述。
什么是宁滨？
宁滨是一种技术，它实现了听起来的效果。它将采用一列连续的数字，并根据我们确定的范围将这些数字放入“箱”中。这将给我们一个新的分类变量特征。
例如，假设我们有一个汽车的数据框架。
汽车样本数据框架
我们将重点放在 MPG 栏，并制作 3 个箱子，燃油效率，平均燃油效率，和耗油量大。
制作垃圾桶
我们的第一步将是确定这些箱的范围。这可能很棘手，因为它可以通过几种不同的方式来完成。
一种方法是根据值的分布均匀划分容器。基本上，我们将经历相同的过程，就好像我们正在创建一个直方图。
由于 Python 使用 matplotlib 简化了直方图的生成，所以让我们在直方图中可视化这些数据，并将 bin 参数设置为 3。
import import matplotlib.pyplot as plt
mpgs = cars['MPG']
plt.hist(mpgs, bins=3)
plt.show()
mpg 直方图
看起来我们的 3 个箱子大约是 10–25，25–40 和 40–56。
让我们通过返回一个显示箱子宽度的数组来确认我们是正确的。
plt.hist(mpgs, bins=3)[1]
\\array([10\. , 25.33333333, 40.66666667, 56\. ])
我们还可以对燃油效率进行一些研究，看看是否有更好的方法来创建我们的 bin 范围。用谷歌快速搜索，很明显，燃料效率是相对的。这是由许多因素决定的，通常与车龄和车型/级别有关。为了使这个例子简单，我们将坚持使用直方图方法。但是，请记住，理解数据是宁滨的一个重要组成部分，一点点研究可以走很长的路。
将 MPG 放入正确的箱子
现在我们知道了箱子的范围，我们可以创建一个简单的函数，利用熊猫超酷的剪切函数。Cut 将使用我们提供的标签名称和范围拆分我们的列。请注意，要包含 10 个，我们必须在 9 点开始第一次切割。
def make_bins(df):
label_names = ["Gas Guzzler", "Average Fuel Efficiency","Fuel Effecient" ]
cut_points = [9, 25.333, 40.667, 56]
df["MPG_effeciency"] = pd.cut(df["MPG"], cut_points, labels=label_names)
return df
我们现在要做的就是将原始数据帧作为参数传递给函数。
new_cars_df = make_bins(cars)
这会给我们……请敲鼓。带有 MPG 类别列的新数据框！
假装一下
洛伊·赫尔曼先生和他的口技木偶。威廉·詹姆斯·哈丁(1826–1899):旺格努伊区的底片。参考文献:1/4–006817-g .亚历山大·特恩布尔图书馆，新西兰惠灵顿。/记录/22322062
最后一步是得到一些虚拟模型，因为我们的模型很可能只能接受整数。
dummies = pd.get_dummies(new_cars_df["MPG_effeciency"])
new_cars_df = pd.concat([new_cars_df, dummies], axis=1)
瞧，我们已经快速生成了新的特征用于训练，并有望提高我们模型的准确性。
在您使用这些新设计的功能进行训练之前，请确保删除 MPG 列，因为重复数据可能会导致一些不良结果。
仅此而已。几分钟之内，我们已经设计了新的功能！
宁滨用熊猫切和 QCut 记录了一个连续变量
入门
何时、为何以及如何将数字特征转换为分类特征
今天，我将使用“西雅图市工资:按性别比较-工资增长职位”数据集来探索宁滨-也称为分组记录-以及单个数字变量。在西雅图开放数据门户的这里找到数据。
给 binners 一个警告:宁滨降低了粒度，并不总是有用的。宁滨通常不用于机器学习模型。
对分箱数据消费者的一个警告:分箱边缘的选择会有很大的影响，特别是在小样本中。当心有人利用宁滨来欺骗或误导你。就推理和其他宁滨选项的结果提问。
何时以及为什么要装箱
我使用宁滨将连续数据分组进行比较。例如，我可以把人们分成不同的年龄组。这允许统计比较和流线型视觉的创造。
同时，这种划分可能有些武断:如果我们以 10 岁为增量分组，30 岁的人可能比 21 岁的人更接近 31 岁的人。
我喜欢在两种情况下绑定:
- 当存在现实生活或商业门槛时。例如，在航空里程方面，当你达到某个阈值时，福利和激励就会改变，从而形成逻辑分组。对于联合航空公司:MileagePlus、卓越银牌、卓越金牌、卓越白金、卓越 1k。
- 确定模式后，当创建适合我的观众的可视化效果时。在今天的例子中，我绘制了每种工作的女性百分比(在城市报告中，具有该职位的员工中有多少百分比被确定为女性)与该职位的平均时薪(在该职位级别的所有员工中):
西雅图市女性工作百分比与平均时薪散点图。趋势线随着平均时薪的增加而减小。
通过最佳拟合线，我们可以看到，随着平均时薪的增加，女性的比例通常会以平均时薪每增加 1 美元女性减少 0.9%的速度下降。我们有一个 3.6 倍 10^-11 的 p 值，因此这种趋势不太可能是由于抽样。这对一些观众来说太难了。
此外，如此多的职位只有很少的雇员，我们得到的这些水平线分别是 0%、50%和 100%的女性。这些突出的方式分散了外卖的注意力。我想把这个有点抽象的散点图转化成一个快速可消化的可视化。我决定去垃圾箱。
如何装箱
有三个主要的宁滨决策:垃圾箱的数量、垃圾箱的类型和标签。我将逐一介绍这些决策的一些考虑因素:
垃圾箱数量:在网上快速搜索可以找到各种各样的经验法则，帮助你决定垃圾箱的数量。(试试这篇维基百科的文章。)
箱柜类型:
- 您的数据有自然的、有意义的间断吗？例如，您可能想要调查在设定的航空里程范围内不同常旅客类别的行为差异。也许你想用联邦退休年龄和典型的高中毕业来将人们分为“青年”、“工作年龄”、“符合退休条件”。
- 您是否希望根据记录的相对定位来拆分记录？例如，你可能想比较收入最高的 1%的人和收入最低的 10%或 50%的人，而不是考虑特定的收入水平。
标签:还是由你决定——你想给你的小组起什么名字？您也可以给它们加上范围标签或整数标签。
一旦你决定了这三件事(或者把这些决定记在心里)，你就可以使用 cut(如果你想根据你的自变量进行分割的话)或者 qcut(qcut)进行样本百分比分割。
对连续变量进行拆分，然后用 cut 对记录进行分类
要根据小时工资率将职称分成五组，并使用相等的 x 轴大小的框:
df['pay_grp_cut_n'] = pd.cut(df['total_avg_hrly_rate'], 5)
这将向 df 添加一个列“pay_grp_cut_n ”,其中每个值都是一个记录所在的 bin 范围。
在 y 轴上显示职称数量会创建一个直方图:
我们可以看到时薪组分为五个 16 美元/小时的区间，从最低平均时薪(5 美元)到最高平均时薪(86 美元)。这种划分导致绝大多数的职位都属于中心和中心以下的群体。
我可以画出这些群体中女性的平均百分比，来说明不同时薪群体中女性百分比的变化:
# get sums within group - esp for employee counts
grouped = df.groupby(df['pay_grp_cut_n'], as_index=False).sum()# set x and y
x = grouped['pay_grp_cut_n'].astype(str)
y = grouped['no_female_empl'] / grouped['total_no_empl'] *100# create visual
fig, ax = plt.subplots(figsize=(12,8))
plt.xlabel('Average Hourly Rate')
plt.ylabel('Percent Female')
plt.title('Percent Female vs. Average Hourly Rate for Seattle City Jobs')
plt.ylim(0,100)
plt.bar(x, y, alpha=.5)
这种形象化将观众的注意力集中在女性代表人数的下降上，因为职称工资增加了。我可以使用这些分类来确定每个分组的统计显著性和/或根据上下文添加误差线。
根据您的需求优化 cut 的一些其他选项(一如既往，在文档中了解更多):
- 调整标签
df[‘pay_grp_cut_l’] = pd.cut(df[‘total_avg_hrly_rate’]
,5
,labels=[‘Very Low’, ‘Low’, ‘Medium’
,‘High’, ‘Very High’])
- 设置容器边界。你可以这样做，为航空公司奖励水平，或获得整数斌休息。
df['pay_grp_cut_set_boundaries'] = pd.cut(df['total_avg_hrly_rate']
,[5,20,35,50,65,80])
使用 qcut 将记录分组到基于卷的组中
QCut或“分位数切割”不是基于标量值进行分割，而是基于落入每个桶中的样本数进行分割。
如果你想要五个大小相等的组，按小时工资排序，你可以使用类似于cut的qcut 。
df['pay_grp_qcut_n'] = pd.qcut(df['total_avg_hrly_rate'], 5)
使用qcut偶数大小的分组，您会看到现在我们的 bin 大小从大约 5 变化到大约 43。每个框代表大致相同数量的职位。
比较cut 和qcut 在每个时段中的职位数量可以看出这一差异:
使用qcut，您可以像设置cut一样设置标签。df[‘pay_grp_qcut_l’] = pd.qcut(df[‘total_avg_hrly_rate’], 5, labels=[‘Very Low’, ‘Low’, ‘Medium’, ‘High’, ‘Very High’])
您可以设置每个桶中的体积百分比来创建不同大小的分组，而不是显式设置边界。例如，如果你在一条曲线上打分，那么你可能想把学生分成(对我来说，任意的):10% F，20% D，20% C，40% B，10% A。你可以使用 qcut 给学生打分:
df['letter_grade'] = pd.qcut(df,
[.1, .2, .2, .4, .1],
labels=['F', 'D', 'C', 'B', 'A'])
何时使用 Cut 与 QCut
这个选择真的取决于你的用例。你有没有很强的商业逻辑或者现有的分类是你想要匹配的或者是你的受众习惯思考的？这将是一个主要的原因削减！其他时候，考虑类似规模的分组或客户的百分比可能更有意义。
无论您决定在您的环境中哪一个有意义，您都可能想要考虑未选择的路径的影响。使用 cut 时，小组规模可能会有很大不同——强调这种变化的幅度(为您自己和/或您的观众)。使用 qcut，铲斗大小可能会有很大的变化——忘记这一点，后果自负。
确定容器边界的另一种方法
Fisher-Jenks 算法在一维数据中找到自然边界，像 k-means，但是 1D。查看jenkspy获得一个不错的实现。您仍然需要仔细考虑要寻找的条柱/间隔的数量，但是可能会得到更多的逻辑分组，尤其是当您的数据集具有不同的间隔时。例如，在一个 4 岁孩子的生日聚会上，身高可以自然地分为“成年人”、“4 岁儿童”、“6 岁儿童”、“2 岁儿童”，尽管偶数或 qcuts 的划分可能与这些分组相去甚远。
编码快乐！你可以在这里找到回购。
照片由 Luca Upper 在 Unsplash 上拍摄
Python 中的二项式分布和二项式检验—统计学
在本文中，我们将探讨 Python 中的二项分布和二项测试。
由 Unsplash 上的 Edge2Edge 媒体拍摄
目录
- 介绍
- 什么是二项分布
- 在 Python 中创建和绘制二项式分布
- 什么是二项式检验(举例)
- Python 中的二项式测试(示例)
- 结论
介绍
为了继续学习本教程，我们需要以下 Python 库:scipy、numpy 和 matplotlib。
如果您没有安装它，请打开“命令提示符”(在 Windows 上)并使用以下代码安装它:
pip install scipy
pip install numpy
pip install matplotlib
什么是二项分布
二项式分布和正态分布是统计学中最流行的分布之一。二项式分布是一系列独立实验( n )中若干成功的离散概率分布( X )。每个实验都有两种可能的结果:成功和失败。成功结局有概率( p )，失败有概率( 1-p )。
注意:一个单独的实验也叫做伯努利试验，一个只有两种可能结果的实验。并且一个实验的二项式分布( n =1)也是伯努利分布。
换句话说，二项式分布模拟了在重复多次的独立实验中观察到成功或失败结果的概率。
假设成功的概率等于: p
那么故障概率等于: q=1-p
所以获得( k )成功和( n-k )失败的概率等于:
作者图片
获得( k )成功的方法数量计算如下:
作者图片
使用上面的符号，我们可以求解概率质量函数(在( n 个)实验中获得( k 个)成功的总概率):
作者图片
注:概率质量函数(PMF)-给出离散随机变量恰好等于某个值的概率的函数。
二项式累积概率函数的公式为:
作者图片
例如:
你正在掷出一个 6 面骰子 12 次，你想找出 5 次结果为 3 的概率。在这里，得到 3 是一个成功的结果，而得到任何其他(1，2，4，5，6)是一个失败的结果。显然，在每一次掷骰子中，你得到 3 的概率是(1/6)。根据我们这里的数据，掷骰子 12 次，结果应该是 3(12 *(1/6))。
现在，我们实际上如何计算 5 次观察到结果 3 的概率？
使用上面的公式，我们可以很容易地解决它。我们有一个实验，它发生了 12 次( n = 12)，有问题的结果数是 5 ( k = 5)，概率是(1/6)或四舍五入后的 0.17(p= 0.17)。
代入上述等式，我们得到:
作者图片
这个实验的二项式分布看起来像这样:
作者图片
你可以清楚地看到，观察 3 作为结果，在 2 次的概率最高，观察 5 次的概率小于 0.05。
在 Python 中创建和绘制二项式分布
现在让我们探索如何创建二项式分布值并使用 Python 绘制它。在本节中，我们将使用三个 Python 库:numpy、matplotlib 和 scipy。
我们将首先导入所需的模块:
对于数据，我们将继续上一节的示例，我们掷骰子 12 次( n = 12)，其中观察到从 1 到 6 的任何数字的概率是(1/6)或 0.17 ( p = 0.17)。
现在我们将在 Python 中为它们创建值:
其中 x 是一个从 0 到 12 的数字数组，代表任何数字可以被观察到的次数。
使用这些数据，我们现在可以计算二项式概率质量函数。概率质量函数(PMF)是给出二项式离散随机变量恰好等于某个值的概率的函数。
在我们的示例中，它将显示从 12 次掷骰子中，您可以观察到任何概率为 0.17 的数字的次数。
构建 PMF:
你应该得到一个有 13 个值的数组(这是我们的 x 值的概率):
[1.06890008e-01 2.62717609e-01 2.95952970e-01 2.02056244e-01 9.31162813e-02 3.05152151e-02 7.29178834e-03 1.28014184e-03 1.63873579e-04 1.49175414e-05 9.16620011e-07 3.41348087e-08 5.82622237e-10]
现在我们有了二项式概率质量函数，我们可以很容易地将其可视化:
您应该得到:
作者图片
现在，试着解释我们看到的东西怎么样？
该图显示，如果我们从 1 到 6(骰子面)中选择任意一个数字，并掷骰子 12 次，观察到这些数字的最高概率是 2 次。
换句话说，如果我选择数字 1，掷骰子 12 次，很可能 1 会出现 2 次。
如果你问，1 出现 6 次的概率是多少？通过查看上面的图表，你可以看到它略高于 0.02%或 2%。
什么是二项式检验
二项式检验是一种单样本统计检验，用于确定二分分数是否来自二项式概率分布。
利用上一节的例子，让我们以一种可以做一些假设检验的方式来重新表述这个问题。以下是情况:
你怀疑骰子偏向数字 3(三个点)，你决定掷骰子 12 次( n = 12)，观察到数值 3(三个点)5 次( k = 5)。您想了解骰子是否偏向数字 3(回想一下，观察到 3 的预期概率是(1/6)或 0.17)。
制定假设我们会有:
作者图片
现在计算概率:
作者图片
这里的概率是显著性检验的 p 值。由于 0.03 < 0.05，我们拒绝零假设，接受骰子偏向数字 3 的替代假设。
Python 中的二项式测试(示例)
现在让我们使用 Python 对上面的例子进行二项式测试。
这是一个非常简单的 scipy 库函数的几行实现。
第一步:
导入函数。
第二步:
定义成功次数( k )，定义尝试次数( n )，定义预期成功概率( p )。
第三步:
在 Python 中执行二项式测试。
我们应该得到:
0.03926688770369119
这是用于显著性检验的p-值(类似于我们通过求解上一节中的公式得到的数值)。
注意:默认情况下，计算的测试是双尾测试。如果您正在处理单尾测试情况，请参考此函数的 scipy 文档。
结论
在本文中，我们探讨了二项式分布和二项式测试，以及如何在 Python 中创建和绘制二项式分布，并在 Python 中执行二项式测试。
如果你有任何问题或对一些编辑有建议，请随时在下面留下评论，并查看更多我的统计文章。
最初发表于 2021 年 10 月 1 日https://pyshark.com。
实践中的二项式分布
数据科学家的数学复习
用 R 和 Python 语言清晰地解释了理论和方便的函数
你知道概率分布很重要，二项分布是基本分布。但是你仍然对此有些怀疑。我应该什么时候使用它？如何理解概率质量函数和累积密度函数的方程？R 和 Python 中有哪些得心应手的函数？最后，从理论和实践的角度来看，你会觉得使用二项式分布很舒服。让我们开始吧！
我们将在本文中讨论以下主题:
- 伯努利试验
- 二项分布
2.1 二项系数
2.2 二项密度函数(PMF)
2.3 累积密度函数(CDF) - 二项分布的行为
- R 和 Python 中的二项式分布函数
4.1 R
4.2 Python - 实践中的二项分布
1.伯努利试验
伯努利试验是一个随机实验，它有两种可能的结果，通常表示为“成功”(1)和“失败”(0)。“成功”符合给定的标准，例如，大于 7 的数字、女性、年龄小于 10 岁、负回报等。它并不意味着这个结果在这个词的伦理意义上是“好”的。
如果 q 是成功的概率，而 p 是失败的概率，那么:
由于除了 0 或 1 没有其他选择，所以成功和失败的概率之和总是等于 1。
2.二项分布
让我们先看看定义。二项式分布描述了 N 次独立伯努利试验中成功次数的分布，其中成功概率为常数，p【1】。简单地说，它的意思是:
- 我们重复同一个实验很多次。比如抛硬币 N=1000 次。
- 在每次试验中，恰好有两种结果是可能的:头和尾(成功和失败，伯努利试验)。因此，这个定义中的伯努利试验将二项分布缩小为离散分布，其中实验结果只能取 0 或 1。任何介于之间的，比如 0.3，0.01，或者既不成功也不失败，都是不可能的。
- 在每次试验中，获得正面的概率是相同的(p=0.5)。它不依赖于任何先前的结果，所以它是独立的。
2.1.二项式系数
为了表达二项式密度函数，我们首先需要定义二项式系数:
我们可以理解为“N 选 m”。它描述了 X 可以取给定值 m 的方式的数量。当选择的顺序无关紧要时，它从集合中选择 m 个项目。如果集合有 N 个元素，m 个组合的数目等于二项式系数[2]。
让我们看看这个例子。在抽签中，从 20 个球中选出 5 个球。需要购买多少张彩票才能确保中奖？一开始我们有 20 个球可供选择，第二选择有 19 个，第三选择有 18 个，依此类推。所以，我们有
二十个球中选择五个球的不同方法。但是我们并不关心彩票的顺序，重要的是有相同的号码。这意味着{A，B，C，D，E}顺序的球和{B，E，A，C，D}对我们来说是同一套球。
所以我们需要减少 1.860.480 种按重复选择的方式(所选球的排序方式数)。既然我们有五个球，我们可以在五分钟内完成！=120 种不同的方式。最后，我们将 1.860.480 种选择方式除以 120 个订购选项，得到 15.504 张票可以购买！
值得吗？
让我们回到二项式系数。在分子中，我们可以看到我们可以用多少种方法来排列 20 个球！(N！).在分母中，我们通过排序选项来减少结果:
- 订购方式 5！选定的球(m！)，以及
- 对池中剩余的球进行排序的可能性(N-m)！，这里 15！
像被选中的球的顺序一样，我们不关心未被选中的球的顺序。最后，我们得到:
和我们上面得到的完全一样，但是方式更简单。
“！”符号是阶乘。正如你可能看到的，对于非负整数 x ，它被计算为所有数字的乘积，直到 x ，例如:
2.2.二项式密度函数(PMF)
现在，我们准备将二项式密度函数定义为在 N 伯努利试验中获得 m 次成功的概率:
因此，二项式分布是一系列给定实验的 N 次独立重复中成功次数( m )的离散概率分布，该实验询问是-否(成功-失败，1–0)问题，成功的概率是 *p，*而失败的概率是 q=1-p 。
让我们考虑一个例子。篮球训练的时候你在练罚球。从赛季统计中，我们知道你得一分的概率是 75%。你的教练告诉你，如果你在 20 次尝试中得到 17 分，你将开始下一场比赛。你正好得 17 分的概率是多少？
我们需要假设罚球成功的概率与之前的结果无关(这里心理力量不发挥作用)。我们也不关心得分的顺序，所以不管你是第一次、第三次还是最后一次失败都没关系。因此，这是一个二项式分布。我们可以使用上面给出的二项式密度函数，得到:
我们可以为其他分数重复这个练习。结果，我们得到了二项分布图(PDF):
给定 20 次伯努利试验中的分数，获得的概率给定成功试验的概率为 75%。作者的情节。
但是你至少要拿 17 分，不是正好 17 分。那么你在下一场比赛中首发的可能性有多大？
2.3.累积密度函数
这里，一个二项分布的累积分布函数将有助于回答这个问题。如果你想知道为什么，请先看看这篇文章:
累积分布函数(CDF)描述了 X 取值等于或小于 k 的概率(机会)。二项式分布的 CDF 函数如下:
其中*【k】为 k 下的“下限”，即等于或小于k*的最大整数
因此，我们需要将你获得 17 次、18 次、19 次或 20 次罚球的概率相加，在 PMF 图上用红色标出:
在 20 次尝试中获得 a 分数的概率，假设成功概率为 75%。作者的情节。
由于概率总和为 1，我们也可以采取相反的方法——从 1 中减去最高得分 16 分(17 分的剩余部分)的累积概率:
这显示在 CDF 图中:
20 次尝试的累积概率，假设成功概率为 75%。作者的情节。
这两种方法的结果是一致的。你的教练给你的任务很有挑战性，但是考虑到概率还是可行的，所以试一试吧！
3.二项分布的行为
如果你从教练那里得到更多的尝试，机会会有怎样的变化？你仍然需要达到 85%的成功投掷，但是你接受了 20、50 或 100 次尝试。让我们看下一个情节。
考虑不同样本量的二项分布的概率密度函数。作者的情节。
样本量越大，分布越广。获得至少 85%罚球成功率的概率如下:
鉴于成功概率为 75%，获得 85%的试验次数的概率。表作者。
考虑到表格中的结果，你不应该强迫你的教练给你额外的试验。尝试的次数越多，获得比你长期成功概率更好的分数的可能性就越小。
接下来，让我们看看艰苦的训练是否增加了你的机会。我们考虑三种成功的可能性:20%、50%和 90%。轨迹的数量保持不变(100)。
考虑不同成功概率的二项分布的 PDF。作者的情节。
下表总结了 85%成功投掷的概率，给出了 75%的长期成功概率。
努力训练比和教练争着跑更多的路更有意义！
4.R 和 Python 中的二项式分布函数
让我们仔细看看 R 和 Python 中有助于处理二项分布的函数。
4.1.稀有
在 r 中至少那四个函数是值得知道的，下面的例子中， m 是成功尝试的次数， N 是样本的大小(所有尝试的次数)， p 是成功的概率。
- dbinom(m，N，p): 该函数计算在成功概率为 p 的 N 次随机试验中恰好有 m 次成功的概率。例如，20 次投掷中的 17 次得分给出 dbinom(17，20，0.75)。如果我们为每个可用的分数计算 dbinom，我们可以绘制一个 PMF 图(如上例所示)。
图片作者在 @carbon_app
- pbinom(m，N，p): 该函数计算累积概率函数，使得成功次数等于或小于 m 。对于至少有 17 个分数的示例，它是 1-pbinom(16，20，0.75)或 pbinom(16，20，0.75，lower.tail=FALSE)。如果我们加上 lower.tail=FALSE，我们考虑获得至少 m 次成功的概率。默认值较低。tail=TRUE。
图片作者在 @carbon_app
-
qbinom(probability P，N，p) 返回逆累积密度函数的值。它找出二项分布的第 p 个分位数。对于给定的概率，你需要多少成功的分数？qbinom(0.55，size=20，prob = 0.75)给出 15。这意味着有 55%的几率你最高可以得 15 分。
-
rbinom (m，N，p) 是生成一个二项分布随机变量的向量。这对于训练目的是有用的。
4.2.计算机编程语言
Python 中那些来自scipy . stats . binom库的函数值得了解。在下面的例子中， m 是成功尝试的次数， N 是样本的大小(所有尝试的次数)， p 是成功的概率。
- pmf(m，N，p): 这个函数计算概率质量函数，所以在 N 随机试验中有精确 m 成功的概率，有 p 成功的概率。比如 20 投 17 分给 scipy.stats.binom.pmf(17，20，0.75)。如果我们计算每个可用分数的 dbinom，我们就可以绘制 PMF。
图片作者在 @carbon_app
作者图片
- cdf(m，N，p) 允许计算累积分布函数，使得成功次数等于或低于 m 。对于我们这个至少有 17 个分数的例子，就是:1-scipy.stats.binom.cdf(16，20，0.75)。
图片作者在 @carbon_app
作者图片
- ppf(m，N，p) 此函数给出了百分点函数，因此它返回了反累积密度函数的值。它找出二项分布的第 p 个分位数。对于给定的概率，你需要多少成功的分数？scipy.stats.binom.ppf(0.55，size=20，prob = 0.75)给出 15。这意味着有 55%的几率你最高可以得 15 分。我们还可以检查 PFM 和 CDF 的准确性:
图片作者在 @carbon_app
- stats(N，p，moments=‘mvsk’) 允许计算前四个矩，即均值(’ m ‘)、方差(’ v ‘)、偏斜(’ s ‘)、峰度(’ k ')。如果你需要重温统计时刻，这篇文章可能适合你:
图片作者在 @carbon_app
- rvs (n，p，size=) 是生成一个二项分布随机变量的向量。这对于训练目的是有用的。
5.实践中的二项分布
二项式分布是统计显著性二项式检验的基础。它旨在检查只有两种可能结果的实验结果是否与预期的结果有统计学差异。
例如，你期望你的投资回报在 95%的时间里都在-100 美元以上。然后，您检查上个月(22 个交易日)的实际回报，您会得到以下结果:
我们能说结果与我们的期望不同吗？首先，我们需要假设每天的收益是独立的。然后，我们将测试的零假设定义为实验结果与预期结果没有显著差异。我们想检查真实数据是否支持另一种选择，即成功的真实概率(这里的回报低于-100 美元)大于 0.05 。假设每天只有成功(低于-100 美元)和失败(高于或等于-100 美元)选项，我们可以使用上面给出的 CDF 等式:
考虑到 5%的阿尔法水平，我们不能拒绝零假设，即结果没有明显不同于预期。因此，即使例外情况的数量高于预期，保持冷静，不要对市场走势反应过度也是一个好策略。
感谢阅读！
我很高兴你看到了这篇文章的结尾。我们讨论了与二项分布相关的最重要的主题。我希望这对你来说是一次激动人心的旅行！
我很乐意在下面的评论区听到你的想法和问题，你可以直接通过我的 LinkedIn 或者 akujawska@yahoo.com 的 T2 联系我。后会有期！
您可能还喜欢:
参考
[1] C .亚历山大(2008):《市场风险分析》。第一卷:金融中的定量方法”，约翰·威利父子有限公司，ISBN 978-0-470-99800-7
[2]维基百科中的组合:【https://en.wikipedia.org/wiki/Combination】T4
【3】a . b .唐尼:“想统计。Python 中的探索性数据分析"
云上的生物信息学
将 RNAseq 与 Google 云平台上的 STAR 对齐
云中的运行之星。图片作者。
关联
分析一个人的 DNA 及其转录指令(RNA)的价格正在迅速变得便宜。因此，不仅在科学实验室，而且在日常医疗保健中，测序技术的应用正在增长。但是，为单个个人生成的数据的大小很容易超过几个 GB。因此，这确实是一个大数据问题！由于云是专门为解决大规模问题而设计的，我们可以利用它的工具和基础设施来促进我们的生物信息学分析。通过在云上运行您的生物信息学工具，您可以在很短的时间内完成分析。
范围
图 1:校准程序的示意图。STAR 使用参考基因组和基因注释将 FASTQs 转换为 SAM 文件。图片作者。
本教程将带你完成在谷歌云平台(GCP)上比对人类 RNA 测序数据(简称 RNAseq)的步骤。
目标:
- 在云上存储基因组数据。
- 使用 STAR 比对 RNA 测序数据。
背景
从一条链上读出核苷酸碱基的测序机器只能读出片段的一小部分。因此，机器会产生大量的小块碱基，称为读作。通常，读取有两种方式:左读和右读。这些读数由成对末端测序机器生成，是机器开始读取碱基的片段的两边。完成测量后，测序仪通常会生成两个文件(每种口味一个)，其中包含所有这些碱基块以及相应的质量评估，格式为FASTQ。所有这些小拼图，现在都是 FASTQ 格式，必须重新组合在一起。这是与参考基因组比对的过程，其中读数被分配到基因组上的一个位置。当成功时，这将产生一个指示读取的基因组位置的 SAM 文件。我们今天的任务是对来自 RNA 材料的碱基进行比对。或者换句话说，将 FASTQ 文件合并并转换为 SAM 文件。
初步步骤
打开谷歌云平台仪表盘，点击右上角的云壳图标(图 2)启动终端。
图 2:谷歌云平台中的云壳图标。图片作者。
我们将利用 dsub 在 Google 云平台上分配我们的生物信息学工作量。要安装dsub，请执行:
**sudo pip3** install dsub
在您的云 shell 中，或者在您的本地机器上，如果您已经在本地安装并配置了gcloud的话。
图 3:在 Google 云平台的云壳中执行 shell 命令。作者截屏。
接下来，选择存储数据的名称，例如gs://rna-seq-tutorial，并创建一个存储桶(mb):
export BUCKET=gs://rna-seq-tutorial
export REFDATA=${BUCKET}/refdata
export INPUT=${BUCKET}/input
export GENOME_INDEX=${BUCKET}/genome_index
**gsutil mb** -l europe-west4 ${BUCKET}
并将注释文件下载到新创建的存储桶中:
**wget** -O - \
ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_36/gencode.v36.chr_patch_hapl_scaff.annotation.gff3.gz \
| **zcat** \
| **gsutil cp** - ${REFDATA}/gencode.v36.chr_patch_hapl_scaff.annotation.gff3
该命令下载文件(wget)、解包(zcat)并将数据复制到 bucket ( gsutil cp)。接下来，复制要与 bucket 对齐的 FASTQ 文件(原始格式的 RNAseq 数据)。如果您没有使用自己的文件，请使用以下示例文件 [1]:
**wget** [http://genomedata.org/rnaseq-tutorial/practical.tar](http://genomedata.org/rnaseq-tutorial/practical.tar)
**tar** -xf practical.tar*# Sample 1.*
**gsutil cp** \
hcc1395_normal_rep1_r1.fastq.gz \
${INPUT}/sample1_R1.fastq.gz
**gsutil cp**\
hcc1395_normal_rep1_r2.fastq.gz \
${INPUT}/sample1_R2.fastq.gz
*# Sample 2.*
**gsutil cp** \
hcc1395_tumor_rep1_r1.fastq.gz \
${INPUT}/sample2_R1.fastq.gz
**gsutil cp**\
hcc1395_tumor_rep1_r2.fastq.gz \
${INPUT}/sample2_R2.fastq.gz
最后，启用云生命科学 API ，以使用dsub。
将表达式数据与星号对齐
为了将 RNA 转录物与参考基因组进行比对，我们将利用STAR【2】。与 STAR 的比对是一个两步过程:
- 使用基因组参考信息生成基因组索引。
- 使用基因组索引比对测序数据。
生成基因组索引
STAR 使用参考基因组和注释文件来生成索引文件。创建一个名为step1.sh的 bash 脚本，内容如下:
#!/bin/bash
**STAR** --runThreadN 8 \
--runMode genomeGenerate \
--sjdbOverhang 100 \
--genomeFastaFiles ${FASTA_FILE} \
--sjdbGTFfile ${GTF_FILE} \
--genomeDir ${GENOME_INDEX}
这将使用 8 个内核来生成索引文件，并将它们写入$GENOME_INDEX(位置将在下面指定)。我们使用--genomeFastaFiles标志传递参考基因组，而基因注释文件可以用--sjdbGTFfile指定。关于命令行参数的更多细节，请参考 STAR 手册。接下来，使用dsub在一个具有 8 个虚拟内核(--min-cores 8)和 48gb 内存(--min-ram 48)的工作线程上执行工作负载，因为 STAR 在这一步需要大量内存:
**dsub** \
--provider google-cls-v2 \
--project <my-project> \
--location europe-west4 \
--zones europe-west4-a \
--preemptible \
--min-ram 48 \
--min-cores 8 \
--logging "${BUCKET}/logging/" \
--input FASTA_FILE=gs://gcp-public-data--broad-references/hg38/v0/GRCh38.primary_assembly.genome.fa \
--input GTF_FILE=${REFDATA}/gencode.v36.chr_patch_hapl_scaff.annotation.gff3 \
--output-recursive GENOME_INDEX=${GENOME_INDEX} \
--image registry.gitlab.com/hylkedonker/rna-seq \
--script step1.sh
决定你是否要冒着工作终止的风险去寻找更便宜的先占计算(--preemptible)，并用你的谷歌云平台项目 ID 替换<my-project>。您可以根据需要调整位置，但建议您相应地调整铲斗位置。请注意，不需要启动虚拟机，也不需要安装我们的生物信息学工具。相反，我们使用了一个 Docker 映像(--image registry.gitlab.com/hylkedonker/rna-seq)，它包含了我们需要的所有东西，而dsub负责处理剩下的事情。当我们的工作负载在云上执行时，我们可以为下一步做准备。
参考基因组比对
为了进行实际校准，创建一个包含以下内容的文件step2.sh:
#!/bin/bash
**STAR** \
--runThreadN 4 \
--readFilesCommand zcat \
--genomeDir ${GENOME_INDEX} \
--readFilesIn ${R1} ${R2} \
--outFileNamePrefix "${OUTPUT}/"
该脚本调用 STAR，使用上一步(--genomeDir ${GENOME_INDEX})中生成的索引在四个内核(--runThreadN 4)上进行实际对齐。注意，每个样本的对齐(在计算上)是一个独立的问题。这意味着所有样品可以平行对齐！使用dsub执行任务数组很容易。简单地创建一个制表符分隔的文件(一个.tsv文件)，其中的参数将被传递给dsub。对于本教程，创建一个名为job2.tsv的文件，其内容如下:
--input R1 --input R2 --output-recursive OUTPUT
gs://rna-seq-tutorial/input/sample1_R1.fastq.gz gs://rna-seq-tutorial/input/sample1_R2.fastq.gz gs://rna-seq-tutorial/output/step2/sample1
gs://rna-seq-tutorial/input/sample2_R1.fastq.gz gs://rna-seq-tutorial/input/sample2_R2.fastq.gz gs://rna-seq-tutorial/output/step2/sample2
用适当的值替换<my-bucket>。另外，请注意，这些列是用制表符分隔的(没有空格)。在执行比对工作负载之前，确保STAR完成了基因组索引的生成。可以使用dstat查询状态。验证作业已成功执行后，使用以下方法将调整工作负载提交到云:
**dsub** \
--provider google-cls-v2 \
--project <my-project> \
--location europe-west4 \
--zones europe-west4-a \
--preemptible \
--min-ram 32 \
--min-cores 4 \
--logging "${BUCKET}/logging/" \
--input-recursive GENOME_INDEX=${GENOME_INDEX} \
--image registry.gitlab.com/hylkedonker/rna-seq \
--tasks ./job2.tsv \
--script step2.sh
这里，我们使用--tasks ./job2.tsv来指代输入和输出参数，每个样本一行。最后，等待工作人员完成(使用dstat来轮询状态)，并且不要忘记在完成计算后删除您的资源，以防止额外的存储成本。
**gsutil rm** -r ${BUCKET}
结论
就这样，您已经使用 STAR 在云上对齐了您的第一个 RNAseq 数据！请注意扩展我们的分析是多么容易:对齐 1000 个而不是 2 个样本只是扩展制表符分隔的文件。
当然，比对本身通常是下游分析所需的预处理步骤。在第二部分中，我们将继续我们的云冒险，从我们比对的 RNAseq 中计算基因表达。
如果您有有用的STAR技巧或有趣的云技巧要分享，或者对代码有任何改进，请发表评论。
参考文献
[1]:玛拉基·格里菲斯*、杰森·r·沃克、尼古拉斯·c·斯皮斯、本杰明·j·安斯考、奥比·l·格里菲斯*。2015.RNA-seq 信息学:云上分析的网络资源。公共科学图书馆比较生物学。11(8):e1004393。
[2]: Dobin，Alexander，等. STAR:超快通用 RNA-seq 比对仪。《生物信息学》29.1(2013):15–21。
云上的生物信息学第二部分
理解大数据
在 Google 云平台上使用 HTSeq 计算 RNAseq 表达式
使用 STAR 和 HTSeq-count 的云上的 RNAseq 分析管道。图片作者。
关联
仅在医疗保健领域，估计到 2025 年将有 6000 万患者的基因组编码被测序[1]。测序成本的下降和医学的进步推动了这一增长。重要的是，细胞水平上的信息可以用来帮助指导患者治疗。但是如何分析这些海量数据呢？我们需要的是能够根据问题的规模进行扩展的工具和基础设施。
进入云技术！
范围
本教程将带你通过谷歌云平台上的 RNAseq 计算基因表达(即每个基因的 RNA 物质数量)的步骤。
目标
- 了解什么是基因表达。
- 使用
samtools为下游分析预处理 SAM 文件。 - 用
htseq-count计算基因表达。
背景
图 1:电池装配线的示意图。DNA 指令通过 RNA 转化成蛋白质。图片作者。
细胞的大部分艰苦工作都是由蛋白质完成的。为了构建它们，细胞表达指令，这些指令被编码在被称为基因的 DNA 中的特定区域。细胞在级联过程中合成蛋白质，其中基因的 DNA 首先作为前体信使-RNA 被读出(转录)。这种前体随后被转化为信使 RNA — 或简称为mRNA。最后，细胞在转移 RNA* 的帮助下将*mRNA 翻译成蛋白质。通过分析 RNAseq 数据，我们可以推断出样本中表达了哪种类型的 RNA 物质，以及有多少。因此，RNAseq 数据向我们展示了细胞装配线中途的快照。虽然这些快照不能给出细胞状态的完整图像，但它们无疑是细胞下游发生的实际过程的有用替代物。通过这种方式，RNAseq 可以帮助了解细胞机制的基础，例如，特定的疾病。
就我们的目的而言，重要的是要认识到 DNA 中的指令并没有一对一地映射到最终产品(蛋白质)上。在转录阶段，只有 DNA 的编码区(基因的外显子)被保留下来，在一个称为剪接的过程中形成转录本。去除非编码区(内含子*)后，剩余的外显子可以以各种方式组合。如图 1 所示，单个 DNA 片段因此可以产生组合数量的相关转录物，称为同种型。反过来，由于控制蛋白质折叠过程的竞争性偶极相互作用，每个同种型翻译成独特的三维结构。*
概括
在之前的中，我们讨论了如何理解在测序仪上生成的大量小碱基阵列(A、C、T、&、G)，称为读数*。我们讨论了如何在一个称为参照基因组比对的过程中把拼图(参照基因组)的碎片(阅读)放在一起。使用 STAR ，我们能够将测序机的原始数据以 FASTQ 格式转换成一组 SAM 文件。本教程将会更进一步。我们将使用 SAM 文件，使用 HTSeq-count [2]计算每个基因产生了多少 RNA，这是一种使用谷歌云平台的可扩展方式。*
基因表达
量化 RNAseq 从计数开始:我们跟踪所有可能的转录物，并计算出现的次数。在其最基本的形式中，这正是基因表达的含义:每个基因的 RNA 计数表。这听起来很简单，但请记住，我们的转录本在测序前就被粉碎成小片段。因此，为了将一段 RNA 分配给一个转录本，我们需要一个转录本目录——一个转录组—* 来选择。与排列好的 RNA 碎片(SAM 格式)一起，htseq-count将汇总数据，建立一个计数表。*
初步步骤
说够了，该编码了！打开谷歌云平台仪表盘，点击右上角的云壳图标(图 2)启动终端。
图 2:谷歌云平台中的云壳图标。图片作者。
我们将利用 dsub 在 Google 云平台上分配我们的生物信息学工作量。要安装dsub，请执行:
***sudo pip3** install dsub*
在你的云壳里。
图 3:在 Google 云平台上的云壳中创建一个存储桶。图片作者。
接下来，我们将需要一个位置，一个桶*，来在云上存储我们的文件。选择一个唯一的名称，例如gs://rna-seq-tutorial，并创建存储桶(mb)😗
***export** BUCKET=gs://rna-seq-tutorial
**export** REFDATA=${BUCKET}/refdata
**export** INPUT=${BUCKET}/input
**gsutil mb** -l europe-west4 ${BUCKET}*
如果你已经成功地完成了之前的教程，在那里我们使用 STAR 对齐了 FASTQ 文件，你应该已经生成了两个文件:sample1.sam和sample2.sam。如果没有，下载并复制样本文件到输入目录:
***wget** -O - [https://tinyurl.com/jnetj2u3](https://tinyurl.com/jnetj2u3) \
| **gsutil cp** - *"${INPUT}/sample1.sam"*
**wget** -O - [https://tinyurl.com/7exr8jsn](https://tinyurl.com/7exr8jsn) \
| **gsutil cp** - *"${INPUT}/sample2.sam"**
我们还需要人类基因组的基因转录本目录(转录组)，通常是gff3或gtf格式。如果你愿意，你可以直接从头构建一个转录组——这被称为转录组组装——但为了简单起见，我们将依赖于 Gencode 之前策划的人类基因组注释。
将转录组文件下载到新创建的 bucket 的 refdata 目录中，如下所示:
***wget** -O - [https://tinyurl.com/vteb57e2](https://tinyurl.com/vteb57e2) \
| **zcat** \
| **gsutil cp** - *"${REFDATA}/gencode.v36.chr_patch_hapl_scaff.annotation.gff3"**
最后，启用云生命科学 API 开始使用dsub。
预处理校准文件
图 4:使用samtools的预处理步骤概述。图片作者。
在我们可以通过计算 RNA 转录物来量化基因表达之前，我们需要预处理我们的 SAM 文件。HTSeq-count 要求 SAM 文件按照染色体位置分类，并附带索引文件，以便快速访问特定位置的读数。这两者都可以通过一个叫做 samtools 的小工具来完成。用以下内容创建一个名为step1.sh的文件(见图 4):
*#!/bin/bash
*# a) Sort SAM file on genomic position.*
**samtools** **sort** \
--threads 2 \
--output-fmt BAM \
-o ${BAM} \
${SAM}*# b) Create index of sorted BAM file.*
**samtools** **index** \
-@ 2 \
${BAM} \
${BAI}*
顾名思义，第一个命令对 SAM 文件进行排序，并将其压缩成(更小的)二进制 BAM 格式(--output-fmt BAM)。然后在两个线程上使用samtools index对生成的 BAM 文件进行索引(-@ 2)。注意，在这一点上，Bash 变量${SAM}、${BAM}和${BAI}还没有定义。当我们将工作负载(step1.sh)提交给云时，这些将由dsub提供。为此，首先创建一个制表符分隔的值文件 step1.tsv，包含以下字段
*--input SAM --output BAM --output BAI
<my-bucket>/input/sample1.sam <my-bucket>/output/step1/sample1_sorted.bam <my-bucket>/output/step1/sample1_sorted.bam.bai
<my-bucket>/input/sample2.sam <my-bucket>/output/step1/sample2_sorted.bam <my-bucket>/output/step1/sample2_sorted.bam.bai*
用您的 bucket 位置替换<my-bucket>(在我们的例子中是gs://rna-seq-tutorial)，并仔细检查这些值是否由制表符(而不是空格)分隔。正如您可以从标题中推断的那样，这些列指的是脚本中的变量${SAM}、${BAM}和${BAI}。创建这个tsv文件允许我们快速地将我们的计算扩展到任意大小，因为dsub将为step1.tsv文件中的每一行启动单独的工作线程。要启动实际的工作负载，请执行
***dsub** \
--provider google-cls-v2 \
--project <my-project> \
--location europe-west4 \
--zones europe-west4-a \
--preemptible \
--min-ram 8 \
--min-cores 2 \
--logging "${BUCKET}/logging/step1/" \
--image registry.gitlab.com/hylkedonker/rna-seq \
--tasks step1.tsv \
--script step1.sh*
决定你是否要冒着工作终止的风险去寻找更便宜的先占计算(--preemptible)，并用你的谷歌云平台项目 ID 替换<my-project>。这将启动两个工作线程(step1.tsv中的每条线一个)，具有 8gb 的 RAM ( --min-ram 8)和两个虚拟内核(--min-cores 2)。每个工人运行一个 Docker 映像(--image)，其中包含运行step1.sh所需的生物信息学工具。Dsub 自动负责将文件从step1.tsv中的--input SAM列复制到step1.tsv中的--output BAM和--output BAI列。
当我们的工作负载在云上运行时，让我们为下一步做准备。
计算基因表达
图 5:使用 htseq-count 估计 RNA 丰度的概述。图片作者。
计算 RNAseq 听起来似乎很简单。但是要知道，大多数现代测序机器并不能完全测量 RNAseq 片段。相反，他们产生仅测量片段两端的成对末端测序数据。更糟糕的是，如何统计一个跨越多个基因的片段？幸运的是，htseq-count用非常合理的缺省值解决了这些问题。
为此(如图 5 所示)，创建一个名为step2.sh的文件，其内容如下:
*#!/bin/bash
# Compute expression
**htseq-count** \
--stranded=reverse \
--order pos \
--format bam \
--additional-attr=gene_name \
--type gene \
-n 2 \
${BAM} \
${GTF} \
-c ${COUNTS}*
这里，我们使用位置排序(--order pos ) BAM 文件(--format bam)和两个内核(-n 2)上的 Gencode 基因注释，运行htseq-count来合计每个基因的计数(--type gene)。需要注意的一个重要方面是，在样品制备过程中是否使用了链识别文库，该文库可以在处理过程中保留 RNA 的链。(如果您不确定， RSeQC 的infer_experiment.py工具可以帮助您识别样品的滞留情况。)在我们的例子中，我们的数据是使用 TruSeq，一个链特定的库[3]生成的，因此我们使用--stranded=reverse标志。类似于我们的预处理计算作业，我们创建一个制表符分隔的文件step2.tsv，它列出了输入 BAM +索引文件和相应的计数文件(以tsv格式)，每个样本一行:
*--input BAM --input BAI --output COUNTS
gs://<my-bucket>/output/step1/sample1_sorted.bam gs://<my-bucket>/output/step1/sample1_sorted.bam.bai gs://<my-bucket>/output/step2/sample1.tsv
gs://<my-bucket>/output/step1/sample2_sorted.bam gs://<my-bucket>/output/step1/sample2_sorted.bam.bai gs://<my-bucket>/output/step2/sample2.tsv*
不要忘记用您的存储位置替换<my-bucket>,并仔细检查您是否正确地复制了制表符(而不是空格)。在提交下一个工作负载之前，让我们通过查询它的状态来检查上一个作业是否成功完成
***dstat** \
--provider google-cls-v2 \
--project <my-project> \
--location europe-west4 \
--status '*'*
如果一切正常，我们准备提交下一个工作负载。
***dsub** \
--provider google-cls-v2 \
--project <my-project> \
--location europe-west4 \
--zones europe-west4-a \
--preemptible \
--min-ram 8 \
--min-cores 2 \
--logging "${BUCKET}/logging/step2/" \
--input GTF=${REFDATA}/gencode.v36.chr_patch_hapl_scaff.annotation.gff3 \
--image registry.gitlab.com/hylkedonker/rna-seq \
--tasks step2.tsv \
--script step2.sh*
与前面的提交命令唯一的细微差别是我们直接向dsub而不是tsv文件提供了--input GTF标志，因为它对于每个工作负载都是相同的。
成功完成工作量后(使用dstat查询状态)，您应该有一个每个样本的计数表(sample<i>.tsv)。您可以通过执行以下命令轻松验证这一点
***gsutil ls** "${BUCKET}/output/step2/"*
显示每个相应样本的tsv文件。(如果没有，您可能想看看日志gsutil ls "${BUCKET}/logging/step2/看看哪里出错了。).
瞧，您刚刚量化了样品的基因表达！
讨论
基因表达计数
在总结之前，我们应该提到，如果您想要比较单个样本，您的计数需要一个额外的标准化步骤。在每次实验运行中，测序仪将给出不同总数的读数*。此外，因为 RNA 转录物在测序前被粉碎成片段，较长的转录物将产生更多的片段。因此，适当选择的变换有助于协调样本的表达估计。幸运的是，大多数下游分析包，如修边机、 limma 和 DeSeq2 会为您完成这项工作，事实上，需要计数数据作为输入。*
Dsub和朋友
在本教程中，我们转向dsub将我们的计算分派给虚拟工作者*。Dsub 非常适合在云上开始，尤其是当您想要做一次或几次特定的分析时。但它绝不是唯一的工具，尤其是如果你想以可重复、容错的方式一遍又一遍地运行相同的工具序列——所谓的工作流。 Snakemake 和 Nextflow 是帮助您在云上运行生产工作流的工具的很好例子。然而，如果你喜欢工具不可知的方法，你可以考虑通用工作流语言 (CWL)，它是一个以 YAML 格式表达工作流的开放标准。*
结束语
就这样，您已经成功地使用 HTSeq 计算了云上的基因表达！注意扩展我们的分析是多么容易:预处理和计数 1000 个样本的 RNAseq 而不是 2 个样本，只需扩展制表符分隔的文件。
如果你有有用的htseq-count技巧或有趣的云技巧要分享，或者对代码有任何改进，请发表评论。
参考文献:
本教程是主要与 Bram van Es 合作的副产品，我们试图揭示免疫疗法(一种用于抗癌的药物)背后的生物模式。
为什么生物学和代码是天作之合
利用 Python 寻找大肠杆菌基因组的 ori 区域
你知道吗，在美国几乎一半的婚姻都以离婚告终。生物学和代码的结合不是那种婚姻，在本文结束时我会向你证明这一点。
在早期的文章中，我在霍乱弧菌的 ori 区域寻找 DnaA 盒。但是如果我们不知道 ori 在基因组中的位置呢？
秘密在于将生物知识转化为基于代码的解决方案。
(This article assumes you have a comfortable grasp on Python and DNA replication.)
“我们来谈谈 DNA 复制吧，小蜜蜂！”—佩帕盐沼，1990 年。
细菌不对称地复制它们的基因组。
- 一半被连续复制: DNA 聚合酶把核苷酸串在一起，没有任何间断。
- **另一半不连续复制:**合成了许多 DNA 的短片段，然后缝合在一起，形成一条连续的链。
更多信息 此处 。
细菌染色体是圆形的，DNA 复制从或区域开始，因此我们可以将染色体分成 4 个半链:
- 反向半股:正向半股的连续合成模板，从 ori 开始，往 3’到 5’方向走。
- 正向半股:反向半股的不连续合成模板，从 ori 开始，往 5’到 3’方向走。
复制从复制泡的两侧进行，并终止于与 ori 区相对的 ter 区。|作者。
在视频中，你有没有注意到前面的半股通常是单股的？这增加了 C 脱氨基(突变)成 T 的机会。
我们可以从这些知识中得出三个结论:
- 脱氨作用在进化史中反复发生，因此正向半链通常含有较少的胞嘧啶。
- 当反向半链在突变的正向半链模板上合成时，DNA 聚合酶将收集 A (与 T 配对)而不是 G (与 C 配对)。反向半链往往会有较少的鸟嘌呤。
最重要的是…
- 沿着基因组的 5’到 3’方向，反向半链在或区域与正向半链相遇。
“我是歪斜的哒-巴-迪-哒-巴-哒”——埃菲尔 65，1998。
让我们称 Gs 和 Cs 的数量之差为“偏斜”:
skew = #G - #C
- 在前半股上，#G > #C
- 在反向半链上，每个窗口中的# G< #C
If we window across the genome and see that skew is rising, we’re finding Cs 更少—我们在正向半链上。偏斜正在减少？更少的 g——我们在反向半链上。
代码
由于 ori 是在反向半股与正向半股相遇的地方发现的，所以我们可以通过观察偏斜何时从减小切换到增大来发现 ori 。也就是最小偏斜。
作者使用 ray.so
取一个字符串Genome作为参数。我们首先创建一个列表Skew，其中第 0 个索引值为 0。Skew中的每个元素都是基因组中该点的#G — #C的连续计数。在第 0 个核苷酸处(第一个核苷酸之前)，Skew = 0。
接下来，我们搜索整个基因组，并将Genome中的第 i 个字符与“G”或“C”进行比较:
- 如果是前者，我们对偏斜的运行计数加 1，并将其添加到列表
Skew。 - 如果是后者，我们-1 和做同样的事情。
- 如果都不是，偏斜保持不变。
下面是它如何发挥作用的示意图:
按作者
这个函数的另一个版本(归功于 Cristian Alejandro Fazio Belan)更加直观:
作者使用雷.索
这里，我们更新运行总数score并将该值追加到列表中。
我们可以把这个列表画出来:
按作者
我们通过编写一个快速函数来返回偏斜最小的所有值，从而找到偏斜最小的点:
由作者使用 ray.so
>>> positions = MinimumSkew(Genome)
>>> print(positions)[3923620, 3923621, 3923622, 3923623]**Alternatively, you can simply do the following:**>>> skew = SkewArray(Genome)
>>> print(skew.index(min(skew)))3923620
对于大肠杆菌，似乎或在第 3923620 个核苷酸附近。
所以我们知道 ori 在哪里。也许吧。我们如何证实这一点？
在早期的文章中，我们发现了如何在基因组序列中找到 DnaA 盒。它们通常是 9 个碱基的序列(9 聚体)，DnaA 蛋白与之结合以启动 DNA 复制。
如果我们在已经找到的点(平均长度为或)切割 500 个碱基，并应用上一篇文章中的 FrequencyMap 和 FrequentWords *函数…
>>>for i in range(3,10):
print(i, ":", FrequentWords(Text, i))3 : ['GAT']
4 : ['GATC']
5 : ['GGATC']
6 : ['AGGATC']
7 : ['AAGGATC']
8 : ['AAAGGATC', 'AAGGATCC', 'AGGATCCG', 'GGATCCGG', 'CTGGGATC', 'GATCCTGG', 'GATCGCAC', 'TGGATAAC']
9 : ['AAGGATCCG', 'AGGATCCGG']
…我们找不到任何出现 3 次以上的 9 聚体/反向互补体。
*(FrequencyMap scans across a sequence and creates a dictionary that contains all possible k-mers in that sequence as keys, and the number of times they occur as values. FrequentWords returns the most frequents k-mers)
生物知识再一次拯救了我们。 DnaA 蛋白质也可以结合到大致相同的 DnaA 盒上。
两个有少量错配的 9 聚体仍能结合 DnaA。
不匹配的次数可以称为汉明距离。让我们写一个找到这个的函数:
作者使用雷.索
将两个字符串p和q作为参数，并遍历两者的每个字符。如果不相等，增加count。
我们可以从频数图中提取 9 聚体，并将其用作更新的模式计数函数的输入:
由作者使用 ray.so
这里，Pattern是我们的任何一个 9 聚体，Text是我们在 ori 取的 500 碱基切片，d是我们能容忍的最大错配数。
我们在假设的或区域范围内，将window与Pattern进行比较，T7 是我们通过频率图找到的任何一个 9 聚体/它们的反向互补体。如果它们有 1 个或更少的不匹配，我们增加计数。
这个函数不是非常高效，但是可以完成工作。|作者使用 ray.so |注意:ReverseComplement 函数的代码可以在这里找到。
这样做，我们就找到了TTATCCACA。它和它的反向互补TGTGGATAA 就是实验确定的大肠杆菌的 DnaA 盒！
但如你所见，他们不是唯一的。其他九聚体同样频繁出现。那么，为什么这很重要呢？重点是什么？
通过收集一些可能的九聚体及其在基因组上的位置，我们可以给实验者指明方向，并帮助他们节省实验室时间。
代码清理干草，这样生物学家可以更快地找到针。
生物学和计算机科学的结合是不寻常的，因为它通常是成功的结合。
哈！找到了。
理解生物学有助于我们做出高效有用的算法。
使用生物数据代码有助于生物学家完善他们的实验。
如果你能驾驭生物和 CS 的联姻，你就拥有超能力。你可以积累特定的知识，很可能变得非常富有。最重要的是，你可以做重要的工作来帮助人们。
感谢阅读！这篇文章是基于我在 Coursera 上参加的 初学者生物信息学 课程的一课。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
