TowardsDataScience 博客中文翻译 2020（六百四十三）

原文：TowardsDataScience Blog

协议：CC BY-NC-SA 4.0

监控你的机器学习模型

原文：https://towardsdatascience.com/monitoring-your-machine-learning-model-6cf98c106e99?source=collection_archive---------22-----------------------

作者图片

MLOps

部署您的 ML 模型时需要注意什么

在过去的几年里，机器学习和人工智能已经越来越成为利用数据的组织的主要内容。随着这种成熟，出现了需要克服的新挑战，例如部署和监控机器学习模型。

尽管部署和监控软件一直是一种久经考验的实践，但对于机器学习模型来说，这样做却有很大不同。监控您的模型可以帮助您了解随着时间的推移您的预测有多准确。

通过监控模型来防止错误的预测

从个人经验来看，似乎有许多组织，主要是中小企业，现在正面临着机器学习的这个生产阶段的挑战。

本文将重点关注在生产中监控你的机器学习模型。然而，如果您不熟悉部署您的模型，我建议您看看下面的文章，让您快速上手:

[## 如何部署机器学习模型

用 FastAPI+uvicon 创建生产就绪的 API

towardsdatascience.com](/how-to-deploy-a-machine-learning-model-dc51200fe8cf)

**注意:**没有一种正确的方法可以监控你的模型。本文仅作为设计监控框架时的参考。您的里数可能因您正在使用的解决方案而异。

1.为什么要监控你的模型？

所有生产软件都容易出故障，每个软件公司都知道监控其性能以防止问题的发生是很重要的。通常，我们监控 软件 本身的质量，而在机器学习的环境中，人们可能更关注监控 预测 的质量。

并非所有数据都是平等的

您希望监控您的模型有几个原因:

• 模型和输入数据之间的关系会发生变化
• 数据的分布发生了变化，使得模型的代表性降低
• 衡量标准和/或用户群的变化改变了变量的潜在含义

漂移

根据以上原因，通常可以在漂移中找到原因。从本质上来说，漂移是数据的统计属性发生变化的现象，这种现象会导致您的预测随着时间的推移而下降。换句话说，由于数据总是在变化，漂移自然会发生。

由 Unsplash 上的罗兰兹·瓦斯伯格拍摄的照片

以网上商店的顾客数据为例。预测模型可以使用诸如他们的个人信息、购买行为和花在广告上的钱等特征。随着时间的推移，这些特征可能不能代表它们最初被训练时的特征。

漂移通常被称为概念漂移、模型漂移或数据漂移

监控您的模型以确保输入的数据与训练时使用的数据相似是很重要的。

为什么不重新训练你的模特呢？

虽然重新训练一个模型听起来总是一个好主意，但它可能比这更微妙！

如果你没有及时的标签怎么办？如果你的预测是几个星期后的事，那么在实践中验证你的模型将会很困难，因为你需要等待几个星期才能得到事实真相。

或者…如果您自动重新训练您的模型，只是让它运行 10 个小时，然后在第二天重新训练它，会怎么样？这可能会花费组织大量的计算时间，但却有一点积极的影响。

或者…如果您以在线方式重新训练一个模型，但由于只关注新添加的数据而导致其性能下降，该怎么办？

盲目地重新训练一个模型会导致更多的成本，时间的浪费，甚至是一个更差的模型。通过监控你的模型，你可以更精确地决定再培训的最佳方法和时间。

2.监视

监控您部署的模型的方法有很多种。实际上，这主要取决于您的应用程序、模型类型、性能度量和数据分布。但是，在监控模型时，有几件事是您最常看到的。

监控预测

从这里检索。

如果你幸运的话，在你的预测和你的目标预测之间几乎没有时间。换句话说，在做了一个预测之后，你很快就会有真正的标签。在这种情况下，只需监控标准绩效指标，如准确性、F1 评分和 ROC AUC。

但是，当您没有及时的标签时，就不能使用性能测量(例如，疾病预测)，您必须寻找其他方法:

• 预测分布*
  您可以使用回归和分类任务中预测的分布和频率来跟踪新的预测集是否与训练数据具有相似的分布。显著的偏差可能表明性能下降。
• 预测概率 很多机器学习模型都可以输出分类概率。这些表明一个模型有多“确定”这是正确的预测。如果这些概率相对较低，那么该模型可能会在部署中遇到困难。

*注意:如果你想比较分布，那么我们通常会关注统计测试，比如学生 t 检验或非参数 Kolmogorov Smirnov 。这个可能会帮助你为你的数据选择正确的统计检验。

监控输入

从这里的检索。

如果您有及时的标签或者没有，那么监视您的模型的输入通常是一个好的实践。这有助于您了解通常输入到模型中的内容，并帮助您跟踪预测的恶化情况。

方法包括:

特征分布 与预测一样，您可以使用分布输入数据来跟踪它是否具有与训练数据相似的分布。如果您发现显著差异，这可能表明您的培训数据并不代表您在生产中发现的数据。

异常值检测 根据您的预处理步骤，您通常不允许在预测模型中使用某些值。变量中类别的数量可能会随着时间的推移而增加，这是您的模型所没有预料到的。

人工监控

此处检索到。

在自动化时代，我为什么要建议人工监控呢？尤其是当人类的行为对错误极其敏感的时候？

嗯……虽然人工监控应该是最后的手段，但在某些情况下，让人类看一看模型生成的预测是有益的。如果一个被监控的预测即使有超过 90%的概率也是不准确的，那么如果这些异常能被人看到就好了。

当你在实践中没有预测的基础事实标签时，人类的监控甚至可能是必要的。

3.阴影模式

每当您想要部署一个重新训练的模型时，您可能想要推迟将其部署到生产环境中。相反，以影子模式部署它可能更好。影子模式是一种通过新训练的模型运行生产数据的技术，无需将预测反馈给用户。

由林奈创作

阴影模式允许您同时运行两个模型，同时在生产环境中测试新模型的性能。在实际部署到生产环境之前，您可以存储其预测以监控其行为。

此外，您可以使用此模式来检查模型是否按预期工作，而无需模拟生产环境，因为它在技术上是生产环境的副本。

感谢您的阅读！

如果你和我一样，对人工智能、数据科学或心理学充满热情，请随时在 LinkedIn 上添加我，或者在 Twitter 上关注我。

有关生产中的机器学习的更多信息，请单击以下帖子之一:

[## 如何部署机器学习模型

用 FastAPI+uvicon 创建生产就绪的 API

towardsdatascience.com](/how-to-deploy-a-machine-learning-model-dc51200fe8cf) [## 借助 Streamlit 快速构建和部署应用

将您的 Streamlit 应用程序部署到 Heroku，展示您的数据解决方案

towardsdatascience.com](/quickly-build-and-deploy-an-application-with-streamlit-988ca08c7e83) [## 数据科学家的单元测试

使用 Pytest 提高管道的稳定性

towardsdatascience.com](/unit-testing-for-data-scientists-dc5e0cd397fb)

面向自动驾驶的单目鸟瞰语义分割

原文：https://towardsdatascience.com/monocular-birds-eye-view-semantic-segmentation-for-autonomous-driving-ee2f771afb59?source=collection_archive---------2-----------------------

2020 年 BEV 语义分割综述

更新:

• 添加 BEV feat 缝线，2021/01/31
• 添加 PYVA，2021/10/01
• 添加全景 BEV，2021/10/04
• TODO:添加 BEV-Seg ， CaDDN ， FIERY ，HDPE pnet。

我还写了一篇关于 BEV 物体检测的更新博文，尤其是关于变形金刚。

[## 自动驾驶中使用变压器的单目 BEV 感知

截至 2021 年末的学术文献和行业实践综述

towardsdatascience.com](/monocular-bev-perception-with-transformers-in-autonomous-driving-c41e4a893944)

自动驾驶需要自我车辆周围环境的准确表示。环境包括静态元素，如道路布局和车道结构，也包括动态元素，如其他汽车、行人和其他类型的道路使用者。静态元素可由包含车道等级信息的高清地图捕捉。

有两种类型的映射方法，离线和在线。离线映射以及深度学习在离线映射中的应用，请参考我之前的帖子。在没有地图支持或无人驾驶汽车从未去过的地方，在线地图会很有用。对于在线绘图，一种传统的方法是 SLAM(同步定位和绘图),它依赖于对一系列图像的几何特征的检测和匹配，或者利用物体的附加概念的变形。

[## 自动驾驶地图中的深度学习

截至 2020 年的文献综述

towardsdatascience.com](/deep-learning-in-mapping-for-autonomous-driving-9e33ee951a44)

这篇文章将关注另一种在线制图的方式——鸟瞰(BEV)语义分割。与 SLAM 相比，SLAM 需要来自同一移动摄像机的一系列图像，BEV 语义分割基于多个摄像机同时观察车辆的不同方向捕获的图像。因此，与 SLAM 相比，它能够从一次性收集的数据中生成更多有用的信息。此外，当 ego 汽车静止或缓慢移动时，BEV 语义分割仍将工作，而 SLAM 将表现不佳或失败。

BEV 语义分段与 SLAM 的输入差异(图片由本文作者提供)

为什么选择 BEV 语义地图？

在一个典型的自动驾驶堆栈中，行为预测和规划通常以这种自上而下的视图(或鸟瞰图，BEV) 完成，因为高度信息不太重要，自动驾驶汽车需要的大部分信息可以方便地用 BEV 表示。这个 BEV 空间可以被不严格地称为 3D 空间。(例如，BEV 空间中的对象检测通常被称为 3D 定位，以区别于成熟的 3D 对象检测。)

因此，标准做法是将高清地图光栅化为 BEV 图像，并在行为预测规划中与动态对象检测相结合。最近探索这一策略的研究包括 IntentNet (优步·ATG，2018 年)司机网 (Waymo，2019 年)交通规则 (Zoox，2019 年) Lyft 预测数据集 (Lyft，2020 年)，等等。

最近的研究将高清地图渲染为 BEV 语义地图(*图片由本文作者编辑，*来源于参考出版物)

诸如对象检测和语义分割的传统计算机视觉任务涉及在与输入图像相同的坐标框架中进行估计。因此，自动驾驶的感知堆栈通常发生在与车载摄像头图像相同的空间——透视视图空间**。**

感知发生在图像空间(左: SegNet )而规划发生在 BEV 空间(右: NMP ) ( 来源)

感知中使用的表示与预测和规划等下游任务之间的差距通常在传感器融合堆栈中弥合，该堆栈通常在雷达或激光雷达等主动传感器的帮助下，将透视空间中的 2D 观测提升到 3D 或 BEV。也就是说，使用 BEV 表示对跨模态感知是有益的。首先，它是可解释的，并有助于调试每种传感模态的固有故障模式。它也很容易扩展到其他新的形式，并简化了后期融合的任务。此外，如上所述，这种表示中的感知结果可以很容易地被预测和计划堆栈使用。

将透视 RGB 图像提升到 BEV

来自主动传感器(如雷达或激光雷达)的数据有助于 BEV 表示，因为测量在 3D 中是固有的度量。然而，由于全景相机传感器的普遍存在和低成本，具有语义意义的 BEV 图像的生成最近吸引了很多关注。

在这篇文章的标题中，“单目”是指管道的输入是从单目 RGB 相机获得的图像，没有显式的深度信息。自主车辆上捕获的单目 RGB 图像是 3D 空间的透视投影，将 2D 透视观察提升到 3D 的逆问题是固有的不适定问题。

挑战，IPM 及其他

BEV 语义分割的一个明显挑战是视图转换。为了正确地恢复 3D 空间的 BEV 表示，该算法必须利用硬的(但可能有噪声的)几何先验，例如相机的内在和外在先验，以及软的先验，例如道路布局的知识集和常识(车辆在 BEV 中不重叠，等等)。传统上，逆透视映射(IPM) 一直是该任务的常用方法，假设平地假设和固定摄像机外部。但是，当相机外部条件变化时，这项任务不适用于不平坦的表面或崎岖不平的道路。

IPM 会在汽车等 3D 物体周围产生伪像(来源)

另一个挑战在于为这样的任务收集数据和注释。一种方法是让无人机一直跟随自动驾驶汽车(类似于 MobileEye 的 CES 2020 talk )，然后要求人类标注语义分段。这种方法显然不具有实用性和可扩展性。许多研究依赖于合成数据或不成对的地图数据来训练提升算法。

在接下来的会议中，我将回顾该领域的最新进展，并强调其共性。根据使用的监控信号，这些研究可以大致分为两种类型。第一类研究借助模拟进行间接监管，第二类研究直接利用最近发布的多模态数据集进行直接监管。

模拟和语义分割

该领域的开创性研究使用模拟来生成必要的数据和注释，以将透视图像提升到 BEV 中。为了弥合模拟到现实(sim2real)领域的差距，许多人使用语义分割作为中间表示。

VPN(查看解析器网络，RAL 2020)

VPN(用于感知周围环境的跨视角语义分割)是最早探索 BEV 语义分割的作品之一，并将其称为“跨视角语义分割”。视图解析网络(VPN)使用视图转换器模块来模拟从透视图到 BEV 的转换。这个模块被实现为一个多层感知器(MLP ),它将 2D 物理范围扩展到一个 1D 向量，然后在其上执行完全连接的操作。换句话说，它忽略了强几何先验，而是纯粹采用数据驱动的方法来学习透视到 BEV 扭曲。这种扭曲是特定于摄像机的，每个摄像机必须学习一个网络。

视图解析网络的整体管道(来源

VPN 使用合成数据(用 CARLA 生成)和对抗性损失进行训练期间的域适应。此外，它使用语义遮罩作为中间表示，没有照片级的纹理间隙。

视图转换器模块的输入和输出具有相同的大小。该论文提到，这使得它很容易插入到其他架构中。在我看来，这实际上是完全没有必要的，因为透视图和 BEV 本质上是不同的空间，因此不需要强制相同的像素格式，甚至输入和输出之间的纵横比也不需要。代码可在 github 上获得。

渔网(CVPR 2020)

渔网将激光雷达、雷达和相机融合转换成 BEV 空间中的单一统一表示。这种表示使得跨不同模态执行后期融合更加容易。视图转换模块(视觉路径中的紫色块)类似于基于 MLP 的 VPN。视图变换网络的输入是一系列图像，但它们只是在通道维度上连接起来并馈入网络，而不是利用 RNN 结构。

渔网使用 BEV 表示来促进跨传感器模态的后期融合(来源

地面实况生成使用激光雷达中的 3D 注记，主要关注动态对象，如车辆和 VRU(易受影响的道路使用者，如行人和骑自行车的人)。其余的都由一个背景类表示。

BEV 语义网格的分辨率为 10 厘米和 20 厘米/像素。这比离线映射中使用的典型值 4 或 5 厘米/像素要粗糙得多。遵循 VPN 的惯例，图像的尺寸与 192 x 320 的输出分辨率相匹配。CVPR 2020 的演讲可以在 Youtube 上找到。

ICRA 2019 年奥运会

VED(具有卷积变分编码器-解码器网络的单目语义占据网格映射)开发了用于语义占据网格映射预测的变分编码器-解码器(VED)架构。它对驾驶场景的前视视觉信息进行编码，并随后将其解码为 BEV 语义占用网格。

用于 VED 中视图变换的可变编码器-解码器网络(来源

此地面实况是使用城市景观数据集中立体匹配的视差图生成的。这个过程可能会有噪声，这实际上促使使用 VED 和从潜在空间采样，以使模型对不完美的 GT 具有鲁棒性。然而，由于是 VAE，它通常不会产生尖锐的边缘，这可能是由于高斯先验和均方误差。

输入图像和输出图像分别为 256×512 和 64×64。 VED 利用了 vanilla SegNet(传统语义分割的相对强大的基线)的架构，并引入了一个 1 x2 池层，以适应输入和输出的不同纵横比。

学习观察周围的物体(ECCV 2018)

学习环视物体寻找室外场景的俯视图表示在 BEV 中产生遮挡区域的幻觉，并利用模拟和地图数据来帮助。

我个人认为，这是 BEV 语义分割领域的一篇开创性论文，但它似乎没有受到太多关注。也许它需要一个朗朗上口的名字？

学习环视物体的整体流程(来源

通过逐像素深度预测和投影到 BEV 来完成视图变换。这部分克服了 BEV 空间中缺乏训练数据的问题。这也在后面的工作中完成，如下面审查的提升、拍打、射击(ECCV 2020) 。

论文用来学习幻觉(预测遮挡部分)的技巧相当惊人。对于 GT 深度很难找到的动态对象，我们过滤掉损失。随机屏蔽图像块，并要求模型产生幻觉。用亏损作为监督信号。

忽略动态区域，但明确删除静态区域以直接监控幻觉(源)

由于在 BEV 空间中很难获得显式成对的监督，因此本文使用对抗性损失来指导具有模拟和 OpenStreetMap 数据的学习，以确保生成的道路布局看起来像真实的道路布局。这个技巧也用在后面的作品中，如monola layout(WACV 2020)。

它在图像空间中使用一个 CNN 进行深度和语义预测，将预测提升到 3D 空间并在 BEV 中渲染，最后在 BEV 空间中使用另一个 CNN 进行细化。BEV 中的这个细化模块还被用在很多其他作品中，如 Cam2BEV (ITSC 2020)和 Lift，Splat，Shoot (ECCV 2020) 。

Cam2BEV (ITSC 2020)

在 Cam2BEV 中输入语义分割图像和预测的 BEV 图(来源

Cam2BEV (一种 Sim2Real 深度学习方法，用于将来自多个车载摄像头的图像转换为鸟瞰图中语义分割的图像)使用具有 IPM 的空间转换器模块将透视特征转换到 BEV 空间。神经网络架构接收由不同摄像机捕获的四幅图像，并在将它们连接在一起之前对每幅图像应用 IPM 变换。

Cam2BEV 使用确定性 IPM 来转换中间特征图(来源)

Cam2BEV 使用从 VTD (虚拟试驾)模拟环境生成的合成数据。它采用了四个语义分割图像，重点放在提升过程中，避免了处理 sim2real 域的空白。

Cam2BEV 具有相当集中的范围和许多设计选择，这使得它非常实用。首先，它只在语义空间起作用，因此避免了 sim2real domain gap 的问题。它有一个预处理阶段，故意掩盖遮挡区域，以避免对遮挡进行推理，并可以使问题更容易处理。为了简化提升过程，它还将“单应图像”作为输入，该单应图像由语义分割结果的 IPM 生成并连接成 360 deg 图像。因此，Cam2BEV 的主要目标是推理 BEV 中 3D 对象的物理范围，其在单应图像中可能被拉长。

预处理后的单应图像和目标 GT 供神经网络在 Cam2BEV 中预测(来源)

Cam2BEV 的目标是校正 IPM，但在某种意义上，IPM 会扭曲 3D 对象，如不在路面上的汽车。然而，它仍然不能处理不平坦的路面或行驶过程中的倾斜度变化。Cam2BEV 的输入输出都是 256x512 像素。代码可在 github 中获得。它还提供了一个很好的 IPM 基线实现。

你所需要的是(多模态)数据集

最近发布的许多多模态数据集( Lyft 、 Nuscenes 、 Argoverse 等)使得直接监督单目 BEV 语义分割任务成为可能。这些数据集不仅提供 3D 对象检测信息，还提供高清地图以及定位信息，以在高清地图上的每个时间戳精确定位 ego 车辆。

BEV 分割任务有两个部分，即(动态)对象分割任务和(静态)道路布局分割任务。对于对象分割，3D 边界框被光栅化到 BEV 图像中以生成注释。对于静态道路布局，基于所提供的定位结果将地图转换到 ego 车辆框架中，并光栅化为 BEV 注释。

单层布局(WACV 2020)

monola layout:单幅图像的 Amodal 场景布局重点是将单个摄像机提升到语义 BEV 空间。本文的重点是研究造成遮挡区域的模型完成情况。似乎受到了学习观察周围物体的严重影响(ECCV 2018) 。

MonoLayout 在 KITTI(左)和 Argoverse (右)的闭塞区域产生幻觉(来源

视图变换通过编码器-解码器结构来执行，并且潜在特征被称为“共享上下文”。两个解码器用于分别解码静态和动态类。作者还报告了在消融研究中使用组合解码器处理静态和动态对象的负面结果。

MonoLayout 的编解码设计(来源

虽然高清地图 groundtruth 在 Argoverse 数据集中可用，但 MonoLayout 选择仅将其用于评估，而不是用于训练(事后诸葛亮还是有意的设计选择？).对于训练，MonoLayout 使用时间传感器融合过程，通过将整个视频中的 2D 语义分割结果与定位信息聚合来生成弱背景真相。它使用 monodepth2 将 RGB 像素提升为点云。它还会丢弃距离自我车 5 米远的任何东西，因为它们可能会很吵。为了鼓励网络输出可以想象的场景布局，MonoLayout 使用了对抗性特征学习(类似于在 学习环视物体 中使用的)。先验数据分布从 OpenStreetMap 获得。

MonoLayout 的空间分辨率为 30 厘米/像素，因此 128 x 128 的输出相当于 BEV 空间中的 40 米 x 40 米。代码可在 github 中获得。

PyrOccNet (CVPR 2020)

PyrOccNet :使用金字塔占据网络从图像预测语义图表示从单目图像预测 BEV 语义图，并使用贝叶斯过滤将它们融合成连贯视图。

PyrOccNet 在 NuScenes 数据集上的演示结果(来源)

PyrOccNet 中视图转换的核心组件是通过密集转换器模块执行的(注意，该转换器不是基于注意力的！)。它似乎受到了来自同一作者的 OFT (BMVC 2019) 的极大启发。OFT 沿着投射回 3D 空间的光线在像素位置均匀地涂抹特征，这非常类似于计算断层摄影中使用的背投算法。PyrOccNet 中的密集变压器模块通过使用 FC 层沿深度轴扩展而更进一步。实际上，在 BEV 空间中，有多个密集的变压器以不同的比例工作，集中在不同的距离范围。

密集的变压器层的灵感来自于这样的观察:虽然网络需要大量的垂直上下文来将特征映射到鸟瞰图(由于遮挡、缺乏深度信息和未知的地面拓扑)，但是在水平方向上，BEV 位置和图像位置之间的关系可以使用简单的相机几何结构来建立。—来自 PyrOccNet paper

密集变压器的 PyrOccNet 架构图(来源

训练数据来自于 Argoverse 数据集和 nuScenes 数据集的多模态数据集，其中既有地图数据，也有 3D 物体检测地面真实

PyrOccNet 使用贝叶斯过滤以连贯的方式融合跨多个摄像机和跨时间的信息。它从二元贝叶斯占用网格的旧思想中汲取灵感，并增强网络输出的可解释性。时间融合结果非常类似于映射过程，并且非常类似于monola layout中用于生成弱 GT 的“时间传感器融合”过程。

仅静态类的 PyrOccNet 的时间融合结果(来源)

PyrOccNet 的空间分辨率为 25 厘米/像素，因此 200 x 200 的输出相当于 BEV 空间中的 50 米 x 50 米。代码将在 github 中提供。

举起，拍打，射击(ECCV 2020)

nuScenes 数据集上的演示(来源)

提升、拍打、拍摄:通过隐式反投影到 3D 对来自任意相机装备的图像进行编码为视图变换执行密集的逐像素深度估计。它首先使用每个摄像机的 CNN 来执行概率性的像素深度预测，以将每个透视图像提升到 3D 点云中，然后使用摄像机 extrinsics 来绘制 BEV。最后，使用 BEV CNN 来改进预测。“拍摄”部分意味着路径规划，将被跳过，因为它超出了本文的范围。

提升和飞溅管道(左)和概率深度预测(右)(来源

通过预测 RGB 图像中像素的深度分布，提出了概率 3D 提升。某种程度上统一了伪激光雷达 (CVPR 2019)的一热提升和 OFT (BMVC 2019)的统一提升。其实这种“软”预测是 可微分渲染 中常用的一招。伪激光雷达 v3 (CVPR 2020)的并发工作也使用了这种软光栅化技巧，使深度提升和投影变得可区分。

训练数据来自于 Lyft 数据集和 nuScenes 数据集的多模态数据集，其中既有地图数据，也有 3D 物体检测地面真实。

Lift-Splat-Shoot 输入分辨率为 128x352，BEV 网格为 200x200，分辨率为 0.5 m/pixel = 100m x 100m。代码可在 github 中找到。

BEV 特征拼接

PyrOccNet 和Lift-Splat-shot都专注于将来自多个相机的同步图像拼接成一个连贯的 360 度视图。 BEV 特征拼接 ( 使用机载单目摄像机理解鸟瞰语义高清地图)将单目视频(具有估计的自我姿态)融合成连贯的前视。

BEV 特征拼接使用 BEV 时间聚合将单目视频转换成 BEV 图

BEV 特征拼接取单目视频(一系列图像)作为模型的输入。为了融合来自多个帧的信息，它引入了 BEV 时间聚合模块。该模块首先投影中间特征图，然后将具有自我姿态(从里程计管道估计)的特征图聚集到连贯和扩展的 BEV 空间中。此外，每个图像帧的中间特征在摄像机空间中用重新投影的 BEV 地面真实来监督。

BEV 特征拼接具有 200x200 像素的 BEV 网格，分辨率为 0.25 m/pixel = 50m x 50m。

PYVA :用心投射你的观点(CVPR 2021)

由于卷积层的局部受限感受野，CNN 很难直接拟合视图投影模型。因此，大多数基于 CNN 的方法明确地将视图变换编码到 CNN 架构中(具有相机外切的特征地图的单应变换等)。

另一方面，由于全局注意机制，变形金刚更适合做这项工作。 PYVA (

在 VPN 之后，PYVA 使用 MLP 将投影图像特征提升到 BEV 中。然后，它使用交叉注意力转换器来进一步细化提升。提升后的 BEV 特征 X’用作查询，图像特征 X(或循环重投影的图像特征 X”)用作键和值。

使用交叉注意力转换器模块背后的思想是直观的:对于查询中的每个像素(BEV 特征)，网络应该关注键(图像特征)中的哪个像素？不幸的是，本文没有展示 transformer 模块内部的一些注意力地图的例子，这些例子可以很好地说明这种直觉。

PYVA 使用交叉视图转换器模块将前视图提升到 BEV 中

PYVA 遵循monola layout的思想，利用对抗性训练损失，促使解码器输出更合理的道路布局和车辆位置。PYVA 更进一步，不仅鼓励网络输出合理的车辆位置，还输出道路和车辆之间合理的关系，因为道路布局在预测车辆位置和方向之前提供了有用的信息(例如，汽车更有可能停在路边)。

全景 BEV

全景 BEV ( 使用单目正视图图像的鸟瞰全景分割)正确地指出实例的概念对下游是至关重要的。 FIERY 将语义分割思想扩展到实例分割，全景 BEV 更进一步，在 BEV 空间进行全景分割。

Panoptic BEV 使用带有 BiFPN(如 EfficientDet)颈部的 ResNet 主干。在每个级别(P2-P5)，图像特征通过密集转换器模块被投射到 BEV 中(注意，这个转换器不是基于注意力的！)。

全景 BEV 的网络架构(transformer her 不是基于注意力的)

在每个密集变压器中，有一个明显的垂直变压器和一个平面变压器。首先，输入图像被二进制语义分割成垂直掩码和水平掩码。被垂直屏蔽屏蔽的图像被馈送到垂直变换器。垂直转换器使用体积晶格来模拟中间 3D 空间，然后将其展平以生成垂直 BEV 特征。这种立式变压器与升降拍击的变压器非常相似。扁平变压器使用 IPM，后接误差校正模块(ECM ),以产生扁平 BEV 特征。这很大程度上遵循了 PyrOccNet 的思路。

密集变压器提升图像到 BEV，并有一个垂直变压器和平面变压器

通过将 IPM 的先验信息注入到总升降管道中，整个设计似乎是建立在升降拍击和 PyrOccNet 之上。模块设计看起来确实很合理，但它们感觉太手工制作和复杂，一旦模型可以访问足够的数据，可能不会导致最佳性能。

从将视图转换模块命名为“密集转换器”和引用“垂直特征”来看， PyrOccNet 的影响也很明显。从演示视频看，结果很有希望。将在 Github 上发布的代码。

局限性和未来方向

尽管在 BEV 语义分割方面已经取得了很大进展，但是在它被广泛部署到生产系统之前，我认为还存在几个关键的差距。

首先，对于动态参与者来说，还没有实例的概念。这使得在行为预测中很难利用动态对象的先验知识。例如，汽车遵循某种运动模型(如自行车模型)，未来轨迹的模式有限，而行人往往具有更多随机运动。许多现有的方法倾向于在语义分割结果中将多个汽车连接成一个连续的区域。

动态语义类不能被重用，并且在很大程度上是“一次性的”信息，而 BEV 图像中的静态语义类(例如道路布局和道路上的标记)可以被视为“在线地图，并且应该被收获和回收。如何聚合多个时间戳上的 BEV 语义分段来估计一个更好的映射是另一个需要回答的关键问题。monola layout和 PyrOccNet 中的时间传感器融合方法可能有用，但需要针对 SLAM 等传统方法进行基准测试。

如何将线上的逐像素语义图转换成轻量级结构化图以备将来重用。为了不浪费宝贵的车载地图周期，在线地图必须转换成某种格式，以便 ego 汽车或其他汽车在未来可以有效地利用。

外卖食品

• 视图变换:现有的许多工作忽略了相机外部的强几何先验信息。应该避免这种情况。 PyrOccNet 和 Lift-Splat-shot 看起来方向正确。
• **数据和监管:**2020 年之前的大部分研究都是基于仿真数据，并使用语义分割作为中间表示，以弥合 sim2real 领域的差距。最近的工作利用多模态数据集对任务进行直接监督，取得了非常有希望的结果。
• 我确实觉得 BEV 空间中的感知是感知的未来，尤其是借助可微分渲染，可以将视图变换实现为可微分模块，并插入端到端模型，直接将透视图像提升到 BEV 空间中。

参考

• VPN :感知周围环境的跨视角语义分割，RAL 2020

• VED :卷积变分编解码网络的单目语义占据网格映射，ICRA 2019

• Cam2BEV :一种模拟现实深度学习方法，用于将来自多个车载摄像头的图像转换为鸟瞰图中的语义分割图像，ITSC 2020

• BEV-Seg :利用几何和语义点云进行鸟瞰语义分割，CVPR 2020 研讨会

• 学习环视物体寻找户外场景的俯视图，ECCV 2018

• 渔网:网格中语义热图的未来推断，CVPR 2020

• monola layout:单幅图像的模拟场景布局，WACV 2020

• https://www . Reddit . com/r/self driving cars/comments/gxq56t/recreating _ the _ Tesla _ autopilots _ birds _ eye _ view/和 github 上的代码

• https://medium . com/asap-report/from-semantic-segmentation-to-semantic-birds-eye-view-in-the-the-Carla-simulator-1e 636741 af 3f

• 一千零一个小时:自驾运动预测数据集，Lyft 预测数据集

• 司机网:模仿最好综合最差学车，RSS 2019

• 道路规则:用语义交互卷积模型预测驾驶行为，CVPR 2019

• IntentNet :学习从原始传感器数据预测意图，CoRL 2018

• OFT :单目三维物体检测的正交特征变换，BMVC 2019

• 来自视觉深度估计的伪激光雷达:填补自动驾驶 3D 物体检测的空白，CVPR 2019

• 伪激光雷达 v3 :基于图像的三维目标探测端到端伪激光雷达，CVPR 2020

• BEV-Feat-Stitching:使用机载单目摄像机理解鸟瞰语义高清地图，Arxiv 2021/01

• PYVA :用心投射你的视角:通过交叉视角变换的单目道路场景布局估计，CVPR 2021

• :使用单目正视图图像的鸟瞰全景分割，Arxiv 2021/08****

• CaDDN :单目 3D 物体检测的分类深度分布网络，CVPR 2021 口述

• 火热:从环绕单目摄像机鸟瞰未来实例预测，ICCV 2021

• BEV-Seg :利用几何和语义点云进行鸟瞰语义分割，CVPR 2020 研讨会

• 高清地图网:在线高清地图构建与评估框架，CVPR 2021 研讨会

自主驾驶中的单目动态物体 SLAM

原文：https://towardsdatascience.com/monocular-dynamic-object-slam-in-autonomous-driving-f12249052bf1?source=collection_archive---------18-----------------------

截至 2020 年的 monoDOS 回顾

更新:

• 2020 年 7 月 11 日，添加来自 S3DOT 作者的 ST3D (CVPR 2020)。

传统的 SLAM(同步定位和绘图)算法通常有一个静态世界假设。即使对于能够在动态环境中运行的实际 SLAM 系统，它们通常也将动态对象视为异常值，并试图在应用常规 SLAM 流水线之前将其过滤掉以获得静态环境。这严重限制了其在自动驾驶中的在线应用，其中动态对象的显式处理至关重要。

单目动态物体 SLAM (MonoDOS) 以两种方式扩展了传统的 SLAM 方法。它是对象感知的，因为它不仅检测和跟踪关键点，还检测和跟踪具有更高级语义含义的对象。它也是动态的，因为它可以处理带有动态对象的场景，并跟踪这些对象的运动。

最好记住，不是所有的对象 SLAM 系统都是动态的，也不是所有的动态 SLAM 系统都是对象感知的。物体 SLAM 的开创性工作是SLAM++**(CVPR 2013)但是它仍然需要一个有静态物体的静态场景。一些动态 SLAM 系统改进了基于刚体和恒定速度约束的姿态估计，但是没有对象的明确概念。**

这篇文章回顾了动态对象 SLAM 领域的几篇最新论文。它主要关注单目方法，以及一些可以修改为单目设置的立体方法。这绝不是一个详尽的综述，如果你推荐其他相关的研究，请告诉我。

动态对象 SLAM 的要素

动态对象 SLAM 系统引入了对象的概念，这有几个含义。首先，它需要有一个来自单个帧的对象提议阶段，就像传统 SLAM 系统中的关键点提议阶段(如 ORB-SLAM 中的 ORB)。该阶段将给出 2D 或 3D 物体检测结果。单目 3D 物体检测的最新进展将在这里大放异彩。第二，它有一个更复杂的数据关联。静态 SLAM 只关心关键点，数据关联只是指用特征向量进行跨帧的关键点匹配。现在我们引入了对象的概念，我们还必须在每一帧中的关键点和对象之间以及跨帧的对象之间执行数据关联。第三，作为对传统 SLAM 中的束调整的自然扩展，现在我们必须添加被跟踪的对象(轨迹线)和这些对象上的动态关键点，可选地具有来自假定运动模型的速度约束。

我制作了下面的图表来捕捉动态对象 SLAM 的三个基本元素。绿色块表示数据关联过程，蓝色块表示光束法平差过程，红色方块表示光束法平差的因子图表示中需要优化的因子。

典型动态对象 SLAM 的流水线和元素

为什么对象？

记住这三个基本元素，我们可以问这样一个问题:为什么在 SLAM 中引入动态对象会有效？首先，假设物体是刚体，随着时间的推移，形状和大小固定。还假设对象具有简单且物理上合理的运动模型，该模型鼓励随着时间的推移而平滑地改变姿势。因此，物体的大小和姿态变化都可以用少量的参数来描述。另一方面，对象的引入可以在束平差阶段制定更大数量的约束。在 N 帧中引入 M 个对象很大程度上增加了 O(M)个未知数和 O(MN)个约束。这使得 SLAM 问题更加健壮。

在接下来的每篇论文的回顾环节中，我将对每篇作品中的上述三个基本要素进行梳理。

立方体 SLAM:单目 3D 物体 SLAM (TRO 2019)

这可能是动态对象 SLAM 中最全面的作品。cubeSLAM 的一个主要贡献是它以优雅的方式将长方体的大小和位置集成到因子图优化中。它还使用一个运动模型来约束长方体可能的运动，从而优化物体的速度。在这个公式下，3D 对象检测和 SLAM 可以相互受益。对象为 BA 和深度初始化提供几何和比例约束。它还可以增加泛化能力，让 orb slam 在低纹理环境中工作。mono3D 结果使用 BA 和运动模型的约束进行优化。

目标提议

利用 2D 物体检测和低层图像特征进行三维长方体提案的生成和评分。这个看似天真的方法对椅子和汽车都非常有效。但是我确信基于深度学习的方法会在这个领域产生更好的结果。

长方体提案生成与评分(来源

数据关联

正如概述中提到的，数据关联必须在多个级别上执行，包括点-点、点-对象和对象-对象。

• 点对点匹配:首先，按照标准 orb-SLAM 的方式，基于 ORB 特征进行关键点匹配。对于不满足极线约束的匹配关键点，它们属于动态对象。对于动态关键点，按特征匹配可能会失败，直接用稀疏光流(KLT 光流)跟踪。在跟踪之后，动态关键点的 3D 位置将被三角测量。****
• 点-对象匹配:对于静态关键点，如果它们属于同一个包围盒，那么它与那个对象相关联。这里使用了许多试探法来确保稳健的性能。例如，点必须在 2 帧的同一个 bbox 内，并且距离长方体中心小于 1 米。框之间重叠区域中的特征点被忽略。****
• 物体-物体匹配:物体匹配通过匹配的关键点间接完成。如果连续帧中的两个对象共享最多的特征点(并且超过 10 个)，则它们作为同一对象被跟踪。在基于特征的匹配或 KLT 跟踪失败的情况下，在边界框级别上使用视觉对象跟踪方法来完成动态对象跟踪。

目标感知光束法平差

静态关键点与相机姿态联合优化，具有与 ORB-SLAM 相同的相机点误差或重新投影误差。在从每个帧获得最佳 3D 建议后，我们可以将其视为 9DoF 3D“测量”,并公式化光束法平差问题。对于静态对象，我们有以下错误术语。

• ****3D 摄像机-物体误差:被跟踪的地标物体具有 6 自由度姿态+3 自由度长方体尺寸，这可以与 9 自由度 3D 测量进行比较，并形成第一误差项。
• 2D 相机-物体误差:从 3D 测量我们可以投射 8 个角到相机图像中。8 个点的最小边界框应该与每帧的 2d 检测 bbox 一致。

请注意，这种 2D-3D 一致性假设并不总是正确的，但对于自动驾驶中车载摄像头的典型视点(水平或稍微向下看)来说，大多数时间都足够接近。

• 物点误差:对于与 bbox 关联的点，根据长方体的中心和大小，它应该位于长方体内。

对于动态对象，cube SLAM 假设动态对象是遵循物理上可行的运动模型的刚体。这引入了两个额外的误差项。

• 运动误差:从 t 帧的一个物体开始，根据分段匀速运动模型，演化状态为 t+1。那么应该和 t+1 时的观测一致。请注意，这涉及到每个动态对象的新状态速度。
• 动态点误差:如果一个点在动态对象上，那么它与动态对象的相对位置是固定的。换句话说，动态点被锚定到相关联的对象。

立方体的因子图。绿色方块是与动态对象相关的因素。为了清楚起见，博客作者添加了绿色/黑色的小文本。(来源)

裁决

CubeSLAM 可以生成一些令人印象深刻的演示，在 KITTI 中使用时间一致的长方体检测和跟踪动态场景中的 3D 对象。然而 SLAM 的结果可能并不总是优于单眼提议。它还可以估计动态对象的速度轮廓，大约在 1 米/秒以内，即使是在分段恒定运动的假设下。请注意，计算里程计(相机姿态)时会考虑对象约束。

SLAM 的 3D 对象检测及其跟踪轨迹、速度估计(来源)

S3DOT :基于立体视觉的语义三维物体和自主驾驶的自我运动跟踪(ECCV 2018)

这项研究是基于立体视频流，但该框架可以扩展到单目 SLAM。本文的主要贡献是展示了使用视频检测和跟踪截断的 3D 对象的能力，否则很难在单个图像上检测到这些对象。

S3DOT 能够预测极端截断汽车的姿态，预测一致的轨迹，并估计被跟踪汽车的速度(来源

目标提议

该管道的灵感来自 Deep3DBox (CVPR 2017，详细配方见我的评论帖)。它使用 2D 包围盒从对象检测和 8 视点分类。它在推断物距之前使用形状尺寸。这使得管道足够通用，可以用在单目设置中。

2D 盒子+视点= 3D 包围盒(来源)

数据关联

• ****对象-对象匹配:通过相似性得分投票来匹配跨帧的 2D 包围盒。在补偿相机旋转之后，相似性分数考虑了边界框形状相似性和中心距离。
• 点-点匹配:与匹配的边界框或物体轮廓(2D 图像中 8 个投影角对着的凸包)内的圆形特征进行关键点匹配，用于静态背景。
• ****点-对象匹配:没有明确说明，但是当一个点在该对象的对象轮廓内时，该点应该与该对象相关联。这种简单的关联策略是由于高效的 3D 提议导致了多边形 2D 对象轮廓，其非常类似于对象实例遮罩。

目标感知光束法平差

对象轮廓外的关键点被计为静态关键点。静态关键点和相机姿态通过相机点误差以与 ORB-SLAM 相同的方式联合优化。在相机姿态(或自我运动)被求解之后，物体姿态被求解。动态对象 BA 有四个误差项。

• 动态点误差(稀疏特征观察):刚性物体上的特征点在物体坐标系中有固定的坐标。
• 2D 相机-物体误差(语义 3D 物体测量):3D 被跟踪物体的投影应适合 2D 测量。
• 对象尺寸一致性:对象形状在帧之间保持一致。这是 cubeSLAM 中 3D 相机-物体错误的一部分。
• 运动误差:时间预测姿态应与单帧 3D 建议一致。确保一致的方向和运动估计的运动学运动模型。它包括车辆亮度、速度、转向角。
• 动态点云对齐:在最小化上述误差之后，我们获得基于尺寸先验的物体姿态的 MPA 估计。为了纠正这个先验中的任何偏差，S3DOT 执行一个名为的动态点云对齐的步骤，以将 3D 长方体与跟踪的点云对齐。这实质上是 cubeSLAM 中的物点误差。**

裁决

S3DOT 展示了在 KITTI 中检测和跟踪动态场景中的 3D 对象的令人印象深刻的结果。请注意，在计算里程计(相机姿态)时，不考虑对象约束。

S3DOT 的作者在 CVPR 2020 发布了另一个作品， ST3D: 立体 3D 物体跟踪的时空联合优化。这项工作主要利用深度学习进行 3D 检测，然后执行联合时空优化。

ClusterVO :对移动实例进行聚类，并为自己和周围环境估计视觉里程表(CVPR 2020)

ClusterVO 提出了一种更通用的动态 SLAM 方法，通过将对象表示为被跟踪的关键点(或本文中的点界标)的集群。**

ClusterVO 的整体管道(来源)

目标提议

ClusterVO 使用 YOLOv3 作为 2D 对象检测器，为每一帧中的对象提出语义 2D 边界框。它并没有假设物体可以用立方体来描述，比如 cubeSLAM 或者 S3DOT。因此不存在长方体生成阶段。

数据关联

ClusterVO 执行一个非常复杂的数据关联方案。它可以被视为两步过程，将观察到的关键点与被跟踪的界标相关联并将边界框与被跟踪的聚类相关联的多级概率关联，随后是将被跟踪的界标与被跟踪的聚类相关联的异构 CRF。

• 点对点匹配:这属于多级概率关联。最近邻和描述符匹配对于静态关键点可能工作得很好，但是对于动态点则不行。因此，首先预测每个集群的位置及其速度(只有线性，没有旋转)。如果观察值 k 在运动预测标志点 i 的投影附近，则将关键点观察值 k 关联到标志点 i 的概率与描述符相似度成比例。
• 物物匹配:这也属于多级概率关联。如果边界框 m 包含来自一个簇的最多数量的投影点，那么它与该簇 q (对象)相关联。
• 点目标匹配:这是最复杂的部分，使用了一个异质条件随机场(CRF)。它确定界标 i 是否与群集 q 相关联。它有多个能量项。一元能量项包括 2D 能量(如果一个点在与一个簇相关联的边界框内，那么它很有可能与该簇相关联。如果它在多个边界框内，则它可以被分配给多个簇)、一个 3D 能量(如果一个点靠近簇的中心，则它具有与簇相关联的更高概率，由簇的大小来调制)和一个运动能量(地标的投影可以由簇的运动来解释)。成对的标签平滑能量项惩罚邻近的界标，如果它们与不同的聚类相关联的话。****

目标感知光束法平差

在概率数据关联之后，我们可以为静态场景和动态集群制定 BA。它使用滑动窗口和专门设计的双轨帧管理来管理关键帧。

• 相机点误差:对于静态场景，clusterVO 联合优化相机姿态和静态关键点的位置，类似于 ORB-SLAM。当 clusterVO 选择滑动窗口状态估计方法时，它还通过附加的边缘化项来增强。这个边缘化术语捕获了由于滑动窗口的有限宽度而被移除的来自观察的贡献。****
• 运动误差:时间预测的姿态应该与从单个帧推断的 3D 测量一致。采用具有从高斯过程采样的加速度的运动模型。ClusterVO 只考虑平移部分。****
• 动态点误差 : clusterVO 也有这个动态点误差，类似于 cubeSLAM 和 S3DOT。如果一个点在动态对象上，那么它与动态对象的相对位置是固定的。

裁决

ClusterVO 是一种更通用的动态对象 SLAM 方法。从 KITTI 动态场景的结果来看，长方体的质量并不像在地标聚类的后处理步骤中估计的那样好。对于自动驾驶，CubeSLAM 和 S3DOT 似乎更实用。请注意，优化里程计(相机姿态)时会考虑对象约束。

MoMoSLAM: 用于动态环境的多目标单目 SLAM(IV 2020)

“多体单声道 SLAM”的概念似乎来自于“多体 SfM ”，但它本质上与动态对象 SLAM 具有相同的含义。

目标提议

MoMoSLAM 使用了一个相当笨重但精确的 3D 对象提议管道。它使用形状先验和关键点将 2D 检测提升到 3D 形状。它首先检测车辆可区分特征上的 k=36 个有序关键点，以及一系列基本形状的变形系数。然后通过最小化重投影误差将 2D 检测提升到三维，得到 6 自由度的姿态和形状参数。这是作者在以前的出版物(IROS 2018)中使用的，与RoI-10D(CVPR 2019)非常相似。

从 2D 关键点到 3D 形状(来源:地球不是平的，IROS 2018)

数据关联

• 点对点匹配:基于描述符特征的关键点匹配，类似于 ORB-SLAM。
• 对象-对象匹配:这在论文中没有明确提到，但是这是肯定需要的。任何 2D 物体追踪方法都可以。
• 点物匹配:未使用。这是通过检测每一帧中每个对象的语义关键点来隐式地和部分地完成的。

目标感知光束法平差

摄像机-物体姿态图和循环一致性(来源)

MoMoSLAM 使用不同的优化公式。如上所示，MoMoSLAM 不会指定每个错误项并将其最小化，而是会强制在姿势图中创建的每个循环保持一致。但本质上，这应该等同于最小二乘误差的最小化。

• 相机点错误:与 ORB-SLAM 相同。由于单目图像中的比例模糊，这种里程计是按比例的。然后，使用逆透视映射(IPM)将地面区域的语义分割和该区域中的点匹配用于估计 3D 深度。这固定了比例因子，并导致公制比例的里程表。

MoMoSLAM 的公制里程计估算(来源

• 多物体姿态图中的循环一致性:姿态图中的节点是估计，姿态图中的边是测量。 摄像机-摄像机边缘通过公制比例里程计进行约束。摄像机-车辆边缘通过单帧 2D 限制为 3D 提升。车辆-车辆边缘通过两种不同的 3D 深度估计方法进行约束(IPM vs 2D-3D 提升)。注意，没有明确的运动模型。****

我觉得循环一致性有点做作，尤其是车-车边缘。为 IPM 和 2D 到三维提升之间的距离估计的一致性增加一个足够的误差项会更直接。

判决

MoMoSLAM 用于固定单目设置的度量比例的方法非常有用。请注意，在计算里程计(相机姿态)时，不考虑对象约束。

https://www.bilibili.com/video/av90800325/

外卖食品

• 动态对象 SLAM 在 SLAM 的 3D 中加入了对象检测和提升，在后端优化中加入了对象的姿态和大小。
• CubeSLAM 和 ClusterVO 联合优化相机姿态和对象姿态。恒定对象大小和刚体运动的附加约束可用于约束因子图优化。这将有助于在提取的关键点稀疏的挑战性场景中计算相机姿态。相比之下，S3DOT 和 MoMoSLAM 将在 ORB-SLAM 失败时失败，因为它们依赖于 ORB-SLAM 进行相机姿态计算。
• CubeSLAM 和 S3DOT 将物体视为长方体，在自动驾驶中更加实用。虽然 ClusterVO 非常通用，但它没有在自动驾驶中纳入许多有用的先验知识，因此缺乏在 3D 对象检测中实现 SOTA 性能的潜力。
• 当无法获得全球位置信息时，CubeSLAM 似乎是在自动驾驶中执行 VIO 的一个很好的候选框架。

参考

• S3DOT : 基于立体视觉的语义三维物体和自主驾驶的自我运动跟踪，ECCV 2018
• 立方体猛击:单目 3D 物体猛击，TRO 2019
• ClusterVO :对移动实例进行聚类并为自己和周围环境估计视觉里程表，CVPR 2020
• MoMoSLAM : 用于动态环境的多目标单目 SLAM，IV 2020
• SLAM++: 在对象层次上同时定位和绘图，CVPR 2013
• M3D-RPN :用于物体检测的单目 3D 区域提议网络，ICCV 2019
• 利用刚体运动实现动态环境下的 SLAM，ICRA 2018
• ST3D: 立体三维目标跟踪的时空联合优化，CVPR 2020

蒙特卡洛辍学

原文：https://towardsdatascience.com/monte-carlo-dropout-7fd52f8b6571?source=collection_archive---------1-----------------------

用一个小技巧免费改善你的神经网络，获得模型不确定性估计作为奖励。

世上没有免费的午餐，至少根据流行的谚语是这样的。嗯，不再是了！也就是说，涉及到神经网络时就不是了。请继续阅读，看看如何用一个简单而巧妙的叫做蒙特卡洛辍学的技巧来提高网络性能。

拒绝传统社会的人

我们将要介绍的魔术只有在你的神经网络有漏层时才有效，所以让我们先简单介绍一下这些。辍学可以归结为在每个训练步骤中简单地关闭一些神经元。在每一步，一组不同的神经元被关闭。从数学上来说，每个神经元都有某种被忽略的概率 p ，称为*辍学率。*退出率通常设置在 0(无退出)和 0.5(大约 50%的神经元将被关闭)之间。确切的值取决于网络类型、图层大小以及网络过拟合训练数据的程度。

一个完整的网络(左)和在特定训练步骤中有两个神经元退出的同一个网络(右)。

但是为什么要这样做呢？放弃是一种正则化技术，也就是说，它有助于防止过度拟合。在数据很少和/或网络很复杂的情况下，模型可能会记住训练数据，因此，在训练期间看到的数据上效果很好，但在新的、看不见的数据上效果很差。这被称为过度适应，辍学寻求缓解。

怎么会？有两种方法可以理解为什么关闭模型的某些部分可能是有益的。首先，信息在网络中传播得更加均匀。想想网络中某处的单个神经元。有几个其他的神经元为它提供输入。使用 dropout，这些输入源中的每一个都可能在训练期间的任何时候消失。因此，我们的神经元不能只依赖一个或两个输入，它必须分散其权重，并注意所有的输入。因此，它对输入变化变得不太敏感，这导致模型更好地泛化。

从我们蒙特卡罗方法的角度来看，辍学效果的另一个解释甚至更重要。因为在每一次训练迭代中，你随机地对每一层中要被丢弃的神经元进行采样(根据该层的丢弃率)，所以每次都有一组不同的神经元被丢弃。因此，每次模型的架构都略有不同，您可以将结果视为许多不同神经网络的平均集合，每个神经网络仅针对一批数据进行训练。

最后一个细节:辍学仅用于培训期间。在推理时，也就是当我们用我们的网络进行预测时，我们通常不会应用任何丢弃——我们希望使用所有训练过的神经元和连接。

蒙特卡洛

现在我们已经解决了辍学问题，蒙特卡洛是什么？如果你正在考虑摩纳哥的一个街区，你是对的！但是还有更多。

摩纳哥蒙特卡洛。Geoff Brooks 在 Unsplash 拍摄的照片

在统计学中，蒙特卡罗指的是一类计算算法，它依靠重复的随机抽样来获得某个数字量的分布。

蒙特卡洛退出:模型准确性

Monte Carlo Dropout 是由 Gal & Ghahramani (2016) 提出的，它巧妙的实现了常规 Dropout 的使用可以被解释为一个众所周知的概率模型:高斯过程的贝叶斯近似。我们可以将许多不同的网络(剔除了不同的神经元)视为来自所有可用模型空间的蒙特卡罗样本。这为推理模型的不确定性提供了数学基础，事实证明，这通常会提高模型的性能。

它是如何工作的？我们只是在测试时应用 dropout，仅此而已！然后，不是一个预测，我们得到许多，每个模型一个。然后我们可以对它们进行平均或分析它们的分布。最好的部分是:它不需要模型架构的任何改变。我们甚至可以在已经训练好的模型上使用这个技巧！为了看到它在实践中的工作，让我们训练一个简单的网络来识别 MNIST 数据集中的数字。

经过 30 个时期的训练，该模型在测试集上取得了 96.7%的准确率。要在预测时打开 dropout，我们只需设置training=True来确保类似训练的行为，即丢弃一些神经元。这样，每个预测都会略有不同，我们可以生成尽可能多的预测。

让我们创建两个有用的函数:predict_proba()生成所需数量的num_samples预测，并对 MNIST 数据集中 10 位数字中的每一位的预测类概率进行平均，而predict_class()只是选择最高的预测概率来挑选最可能的类。这个和下面的一些代码片段是受 Geron (2019)的启发。该书附有一套优秀的 jupyter 笔记本。

现在，让我们在测试集上进行 100 次预测并评估准确性。

这产生了 97.2%的准确度。与之前的结果相比，我们将错误率从 3.3%降低到了 2.8%，降低了 15%，完全没有改变或重新训练模型！

蒙特卡洛退出:预测不确定性

让我们来看看预测不确定性。在分类任务中，从 softmax 输出中获得的类概率经常被错误地解释为模型置信度。然而，Gal&Ghahramani(2016)表明，即使 softmax 输出很高，模型的预测也可能不确定。我们也可以在 MNIST 预测中看到这一点。让我们比较一下 softmax 输出和蒙特卡洛预测的单个测试示例的退出概率。

soft max _ output:[0\. 0\. 1\. 0\. 0\. 0\. 0\. 0\. 0\. 0.]
MC _ pred _ proba:[0\. 0\. 0.989 0.008 0.001 0\. 0\. 0.001 0.001 0\. ]

双方都同意测试示例最有可能来自第三类。然而，soft max 100%确定是这种情况，这应该已经提醒您有些事情不对劲了。0%或 100%的概率估计通常是危险的。蒙特卡洛辍学为我们提供了更多关于预测不确定性的信息:最有可能是第 3 类，但也有小概率是第 4 类，例如，第 5 类虽然不太可能，但仍然比第 1 类更有可能。

蒙特卡洛辍学:回归问题

到目前为止，我们已经讨论了一个分类任务。现在让我们来看一个回归问题，看看蒙特卡洛退出是如何为我们提供预测不确定性的。让我们使用波士顿住房数据集拟合一个回归模型来预测房价。

对于分类任务，我们定义了函数来预测类别概率和最可能的类别。同样，对于回归问题，我们需要函数来获得预测分布和点估计(让我们使用平均值)。

让我们再次对一个测试示例进行 100 次预测，并绘制它们的分布，标记其平均值，这是我们的点估计，或最佳猜测。

来自波士顿住房数据的一个测试示例的预测价格分布。红线表示平均值。

对于这个特定的测试示例，预测分布的平均值为 18，但是我们可以看到，其他值也不是不可能的—该模型对其预测不是很确定。

蒙特卡洛辍学:一个实施细节

最后一点:在本文中，我们一直通过将模型的training模式设置为true来实现蒙特卡洛脱离。这很有效，但是它可能会影响模型的其他部分，这些部分在训练和推断时表现不同，例如批处理规范化。为了确保我们只打开 Dropout 而不影响其他任何东西，我们应该创建一个自定义的 MonteCarloDropout 层，该层从常规 dropout 继承，并在默认情况下将其training参数设置为true(以下代码改编自 Geron (2019))。

结论

• 蒙特卡洛退出归结为用常规退出训练神经网络，并在推理时保持其开启。这样，我们可以为每个实例生成多个不同的预测。
• 对于分类任务，我们可以平均每个类的 softmax 输出。这往往会导致更准确的预测，从而更恰当地表达模型的不确定性。
• 对于回归任务，我们可以分析预测分布，以检查哪些值是可能的，或者使用其平均值或中值对其进行汇总。
• Monte Carlo Dropout 很容易在 TensorFlow 中实现:它只需要在进行预测之前将模型的training模式设置为true。最安全的方法是编写一个自定义的三行类，继承常规 Dropout。

来源

• Gal Y. & Ghahramani Z .，2016 年，辍学作为贝叶斯近似:表示深度学习中的模型不确定性，第 33 届机器学习国际会议论文集
• Geron A .，2019，第二版，使用 Scikit-Learn 和 TensorFlow 进行机器学习:构建智能系统的概念、工具和技术

感谢阅读！

如果你喜欢这篇文章，为什么不订阅电子邮件更新我的新文章呢？并且通过 成为媒介会员 ，可以支持我的写作，获得其他作者和我自己的所有故事的无限访问权限。

需要咨询？你可以问我任何事情，也可以在这里 为我预约 1:1 。

你也可以试试我的另外一篇文章。不能选择？从这些中选择一个:

** [## SVM 内核:他们实际上做什么？

直观的视觉解释

towardsdatascience.com](/svm-kernels-what-do-they-actually-do-56ce36f4f7b8) [## 校准分类器

你确定你的模型返回概率吗？🎲

towardsdatascience.com](/calibrating-classifiers-559abc30711a) [## 线性回归中收缩法和选择法的比较

详细介绍 7 种流行的收缩和选择方法。

towardsdatascience.com](/a-comparison-of-shrinkage-and-selection-methods-for-linear-regression-ee4dd3a71f16)**

蒙特卡罗积分和抽样方法

原文：https://towardsdatascience.com/monte-carlo-integration-and-sampling-methods-25d5af53e1?source=collection_archive---------15-----------------------

积分近似和随机数发生器

积分是解决问题时经常使用的关键计算。对于概率任务，连续随机变量 x 的期望值由以下积分定义，其中 p(x)是 x 的概率密度函数。

然而，用计算机计算这个值并不容易。为了有效地定义这个积分，使用了数值近似方法。在这里，我将介绍一个简单的近似方法，蒙特卡罗积分

蒙特卡罗积分

蒙特卡罗积分是一种数值积分计算方法，使用随机数来近似积分值。考虑下面 f(x)的期望值的计算。这里，p(x)是 x 的概率密度函数。

在这种方法中，我们从概率恒等函数 p(x)中选择 n 个独立同分布的样本{x_i} (i=1，2，…，n)。对于所有样本 x，积分值近似为 f(x)的平均值。

根据大数定律，当 n →∞收敛于 f(x)的期望值时，蒙特卡罗积分的一个极限。

让我们用蒙特卡洛积分法计算一个***【π】***的近似值。考虑一个边长为 1 的正方形和一个单位圆，如下图所示。

圆的面积是 π。现在假设 f(x，y)是一个函数，当(x，y)位于圆内时输出 1，否则输出 0。设 p(x，y)是[-1，1]上的均匀分布。在这种情况下，单位圆 π 的面积可以写成

用蒙特卡罗积分，这个积分( π )可以用 p(x，y)的 i.i.d 样本来近似。n(in circle)是位于单位圆内的样本数。

下图显示了 n 足够大时蒙特卡罗积分的收敛性。

有了蒙特卡罗积分，我们可以简单地用随机数计算积分值。但是，问题是“*如何生成随机数？”*当然，如果概率分布是众所周知的，例如均匀分布、高斯分布，我们可以很容易地用一些库实现随机数生成器。如果我们必须从任何库中都没有实现的分布密度函数中生成随机数，那该怎么办？为了解决这个问题，使用了抽样方法。

重要性抽样

一种常用的方法是重要性抽样。在这种方法中，引入了一种代理分布来从任意分布中抽取随机数。在大多数情况下，我们选择众所周知的分布如高斯分布、均匀分布作为代理分布。这个方法的主要概念可以简单地写成下面的形式，其中 q(x)是一个代理分布。

利用重要性采样，我们从代理分布 q(x)中选择 i.i.d .样本{x’_i} (i=1，2，…，n ),而不是从 p(x)中生成随机数，并通过以下计算来近似积分值。这里，p(x)/q(x)被称为抽样的*。*

现在，让我们以拉普拉斯分布的方差计算为例。
考虑 f(x)=x 和一个概率密度函数 p(x)=1/2 Exp(-|x|)。具有类似 p(x)的密度函数的分布称为拉普拉斯分布。
如果我们选择均匀分布作为代理分布，拉普拉斯分布的方差可以近似计算为

**

有了纸和铅笔，我们可以很容易地计算 Var[x]。这个计算的值是 2。现在，让我们确认重要性抽样方法的结果。

在重要性抽样中，我们使用从代理分布生成的随机数(样本)来近似积分值。为了从具有任何概率密度函数的分布中直接生成样本，我介绍了以下采样方法。

逆变换采样

在逆变换采样法中，我们用一个由 1 维均匀分布生成的随机数 u 来生成任意 1 维概率密度函数 p(x)的随机数 x。在这种情况下，我们使用 p(x)的累积分布函数的反函数。如果 P(x)的累积分布函数是 p(x)，那么 u=P(x)的反函数是 x =P^-1 (u)。现在，x = P^-1 (u)有 p(x)作为[0，1]上的概率密度函数。因此，对于[0，1] {u_i} (i=1，2，…，n)上均匀分布的 n 个样本，我们可以通过计算 x_i = P^-1(u_i).来生成 p(x)分布{x_i} (i=1，2，…，n)的 n 个样本

再次，让我们考虑拉普拉斯分布的方差计算作为一个例子。这一次，我们使用逆变换采样方法，从拉普拉斯分布的概率密度函数直接生成随机数(样本)。有了这些随机数，我们将再次重新计算 Var[x]的近似值。

让我们检查这个方法的结果。

下图显示了使用逆变换采样方法的拉普拉斯分布方差 V[x]的近似值。

从这个结果可以直接看出，与代理分布的情况相比，使用从分布中直接生成的随机数提供了更好的近似。然而，利用逆变换采样方法，不可能仅从 2 维或更高维分布直接生成随机数。为此，使用了拒绝抽样方法。

拒绝抽样

拒绝抽样的概念

拒绝抽样的思想是使用代理分布(高斯或均匀分布等。)调用 q(x)生成一个随机数并使用另一个均匀分布来评估生成的样本是否接受它为从 p(x)生成的样本。通过这种方法，我们还可以从更高维度的分布中生成随机数。

作为用这种方法生成随机数的准备，我们需要知道 L 的一个有限值，其中 max[p(x)/q(x)] < L. Here, q(x) is a proxy distribution.

• First, we generate a random number x’ from a proxy distribution q(x). This x’ is called a 建议点 。
• 接下来，从[0，1]上的均匀分布生成随机数 v。这个 v 将用于评估建议点，考虑从 p(x)中生成是否是好的。
• 如果 v ≤ p(x’)/q(x ')，则 x '被接受为 p(x)生成的随机数，否则，x '被拒绝。

用拒绝抽样产生 n 随机数的算法是

随机数生成和拉普拉斯分布方差计算的结果(代理分布:高斯分布)

**

马尔可夫链蒙特卡罗方法

在拒绝采样方法中，当不知道上边界 L 时，不可能产生随机数。MCMC 方法是解决这一问题的有效方法。MCMC 方法使用随机过程的概念(在这种情况下是马尔可夫链)。在这种情况下，第 I 个样本 x_i 的生成依赖于前一个样本 x_(i-1)。x1，x2，…，xn 用这个概念叫做马尔可夫链。这里，我介绍一种 MCMC 方法，Metropolis-Hastings 方法。

这种方法的过程类似于拒绝抽样。但是在这里，代理分布密度函数由条件概率 q(x|x_i)表示，并且评价指标 v 由[0，1]上的均匀分布生成。

• 首先，我们从代理分布 q(x|x_i)生成一个随机数 x’。这个 x '叫做一个 建议点 。
• 接下来，从*【0，1】*上的均匀分布生成随机数 v。这个 v 将用于评估建议点，考虑从 p(x)中生成是否是好的。
• 如果 v≤p(x ')q(x _ I | x ')/(p(x _ I)q(x ’ | x _ I))，则 x '被接受为 p(x)生成的随机数，否则，x '被拒绝。

用拒绝抽样产生 n 个随机数的算法是

随机数生成和拉普拉斯分布方差计算的结果(代理分布:高斯分布)

**

拉普拉斯分布的方差近似

示例代码

*import random
import numpy as np
from scipy.stats import norm*

用蒙特卡罗积分计算 Pi(模拟参数 n:样本数)

*def f(x,y):
    if (x*x+y*y<=1):
        return 1
    else:
        return 0n = 10
N = []
P = []
while n <= 100000:
    s = np.random.uniform(-1,1,(n,2))
    nn = 0
    for x in s:
        nn+=f(x[0],x[1])

    pi = 4*nn/n
    N.append(n)
    P.append(pi)
    n = n*2*

拉普拉斯分布方差的重要抽样计算(模拟参数 n:样本数)

 *def f(x):
    return x*x
def p_pdf(x):
    return 0.5*np.exp(-abs(x))
def imp_pdf(x):
    return norm.pdf(x)n = 10
N = []
P = []
while n <= 100000:
    s = np.random.normal(0,1,n)
    nn = 0
    for x in s:
        nn+=f(x)*p_pdf(x)/imp_pdf(x)

    p = nn/n
    N.append(n)
    P.append(p)
    n = n*2*

逆变换采样(随机数生成器和积分近似计算)

*def sgn(x):
    if x>0:
        return 1
    elif x==0:
        return 0
    else:
        return -1

def p_laplace(x):
    return 0.5*np.exp(-abs(x))#def P_laplace_cum(x):
    #return 0.5*((1+sgn(x))(1-np.exp(-abs(x))))def P_laplace_inv(u):
    k = 1-2*abs(u-0.5)
    return -1*sgn(u-0.5)*np.log(k)# random number generatordef reverse_generator(X):
    res = []
    for u in X:
        res.append(P_laplace_inv(u))
    return res n_sample = 50
s = np.random.uniform(0,1,n_sample)
theta = reverse_generator(s)

print("n_sample = ",n_sample)
gen_xx =random.sample(theta,n_sample)
plt.scatter(gen_xx,[0]*n_sample,marker='x')
xx=np.linspace(-15,15,1000)
plt.plot(xx,[p_laplace(i) for i in xx],c='red')
plt.legend(['Laplace Distribution','generated number'])
plt.title("Random Numbers directly generated from p(x)")
plt.show()# Laplace distribution variance calculationN = []
P = []
n= 10
while n <= 100000:
    s = np.random.uniform(0,1,n)
    theta = reverse_generator(s)
    p = 0
    for x in theta:
        p += x*x
    N.append(n)
    P.append(p/n)
    n = n*2

P_r = [2]*len(N)
plt.plot(N,P)
plt.plot(N,P_r,c='red')
plt.legend(["Inverse Transform Sampling value","true value(2)"])
plt.xlabel("number of samples (n) ")
plt.show()*

拒绝采样(随机数生成器和积分近似计算)

*def p_laplace(x):
    return 0.5*np.exp(-abs(x))def p_proxy(x):
    return norm.pdf(x)# random number generatordef rejection_generator(n):
    xx = np.linspace(0,1,100)
    K = max(p_laplace(xx)/p_proxy(xx))+1
    random_samples = []
    while len(random_samples) < n:
        proposal = np.random.randn()
        v = np.random.uniform(0,K)
        if v <= (p_laplace(proposal)/p_proxy(proposal)):
            random_samples.append(proposal)
    return random_samplesprint("generating...")
random_samples = rejection_generator(n=50)
print("n_sample = ",len(random_samples))
plt.scatter(random_samples,[0]*len(random_samples),marker='x')
xx=np.linspace(-15,15,1000)
plt.plot(xx,[p_laplace(i) for i in xx],c='red')
plt.legend(['true','gen'])
plt.show()# Laplace distribution variance calculationN = []
P = []
n= 10
while n <= 100000:
    theta = rejection_generator(n)
    p = 0
    for x in theta:
        p += x*x
    N.append(n)
    P.append(p/n)
    n = n*2

P_r = [2]*len(N)
plt.plot(N,P)
plt.plot(N,P_r,c='red')
plt.legend(["Rejection Sampling value","true value(2)"])
plt.xlabel("number of samples (n) ")
plt.show()*

MCMC 方法(随机数生成器和积分近似计算)

*def p_laplace(x):
    return 0.5*np.exp(-abs(x))def p_proxy(x,theta,sigma):
    return norm.pdf(x,theta,sigma)def MCMC_generator(n):
    random_samples = []
    theta_i = 0
    while len(random_samples) < n:
        proposal = np.random.normal(loc=theta_i,scale=1)
        v = np.random.uniform(0,1)
        a = p_laplace(proposal)*p_proxy(theta_i,proposal,1)
        b = p_laplace(theta_i)*p_proxy(proposal,theta_i,1)
        if v <= (a/b):
            random_samples.append(proposal)
            theta_i = proposal
    return random_samplesn = 50
random_samples = MCMC_generator(n)
print("n_sample = ",n)
plt.scatter(random_samples,[0]*n,marker='x')
xx=np.linspace(-15,15,1000)
plt.plot(xx,[p_laplace(i) for i in xx],c='red')
plt.legend(['true','gen'])
plt.show()N = []
P = []
n= 10
while n <= 10000:
    theta = MCMC_generator(n=n)
    p = 0
    for x in theta:
        p += x*x
    N.append(n)
    P.append(p/n)
    n = n*2

P_r = [2]*len(N)
plt.plot(N,P)
plt.plot(N,P_r,c='red')
plt.legend(["Rejection Sampling value","true value(2)"])
plt.xlabel("number of samples (n) ")
plt.show()*

参考

1. *M.Yamasugi. *统计机器学习-基于生成模型的模式识别-。2019 年，日本奥姆沙

蒙特卡罗积分

原文：https://towardsdatascience.com/monte-carlo-integration-db86b8d7beb3?source=collection_archive---------25-----------------------

Pixabay 上由 Pixapopz 拍摄的照片

很多时候，我们不能解析地求解积分，必须求助于数值方法。其中包括蒙特卡罗积分。你可能还记得，函数的积分可以解释为函数曲线下的面积。蒙特卡罗积分的工作原理是在 a 和 b 之间的不同随机点计算一个函数，将矩形的面积相加，取和的平均值。随着点数的增加，结果越来越接近积分的实际解。

用代数表示的蒙特卡罗积分:

与其他数值方法相比，蒙特卡罗积分特别适合计算奇数形状的面积。

在下一节中，我们将看到当我们知道先验和似然性，但缺少归一化常数时，如何使用蒙特卡罗积分来确定后验概率。

简而言之贝叶斯统计

后验概率指的是贝叶公式中的一个特定项。

Baye 定理可以用来计算一个在特定疾病的筛选测试中测试呈阳性的人实际上患有该疾病的概率。

如果我们已经知道 P(A)，P(B)和 P(B|A ),但还想知道 P(A|B ),我们可以使用这个公式。例如，假设我们正在测试一种罕见的疾病，这种疾病会感染 1%的人口。医学专家已经开发了一种高度敏感和特异的测试，但还不十分完善。

• 99%的病人测试呈阳性
• 99%的健康患者测试呈阴性

贝叶斯定理告诉我们:

假设我们有 10，000 人，100 人生病，9，900 人健康。此外，在对所有人进行测试后，我们会让 99 名病人进行测试，但也让 99 名健康人进行测试。因此，我们会以下面的概率结束。

p(生病)= 0.01

p(未患病)= 1–0.01 = 0.99

p(+|有病)= 0.99

p(+|未患病)= 0.01

贝叶斯定理在概率分布中的应用

在前面的例子中，我们假设一个人生病的先验概率是一个精确到. 001 的已知量。然而，在现实世界中，认为. 001 的概率实际上如此精确是不合理的。一个特定的人生病的概率会根据他们的年龄、性别、体重等有很大的不同。一般来说，我们对给定先验概率的了解远非完美，因为它是从以前的样本中获得的(这意味着不同的总体可能给出不同的先验概率估计)。在贝叶斯统计中，我们可以用先验概率的分布来代替这个值 0.001，该分布捕获了我们对其真实值的先验不确定性。包含一个先验概率分布最终产生一个后验概率，它也不再是一个单一的量；相反，后验概率也变成了概率分布。这与传统观点相反，传统观点认为参数是固定的量。

归一化常数

正如我们在关于吉布斯采样和 Metropolis-Hasting 的文章中所看到的，当归一化常数未知时，蒙特卡罗方法可用于计算后验概率分布。

让我们首先探究一下为什么我们需要一个归一化常数。在概率论中，归一化常数是一个函数必须乘以的常数，因此其图形下的面积为 1。还不清楚？让我们看一个例子。

回想一下，正态分布的函数可以写成:

两个圆周率的平方根是归一化常数。

让我们检查一下我们是如何决定它的。我们从以下函数开始(假设均值为 0，方差为 1):

如果我们把它画出来，它会形成一个钟形曲线。

问题在于，如果我们取曲线下的面积，它就不等于 1，这是概率密度函数所要求的。因此，我们将函数除以积分结果(归一化常数)。

回到手头的问题，即如何在没有归一化常数的情况下计算后验概率…事实证明，对于连续的样本空间，归一化常数可以改写为:

在这一点上，你应该考虑蒙特卡罗积分！

Python 代码

让我们看看如何通过在 Python 中执行蒙特卡罗积分来确定后验概率。我们首先导入所需的库，并设置随机种子以确保结果是可重复的。

import os
import sys
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy.stats as stnp.random.seed(42)

然后，我们设置贝塔分布和二项式分布的参数值。

a, b = 10, 10
n = 100
h = 59
thetas = np.linspace(0, 1, 200)

概率密度函数的范围从 0 到 1。因此，我们可以简化方程。

在代码中，前面的等式如下所示:

prior = st.beta(a, b).pdf(thetas)
likelihood = st.binom(n, thetas).pmf(h)
post = prior * likelihood
post /= (post.sum() / len(thetas))

最后，我们可视化的概率密度函数的先验，后验，和可能性。

plt.figure(figsize=(12, 9))
plt.plot(thetas, prior, label='Prior', c='blue')
plt.plot(thetas, n*likelihood, label='Likelihood', c='green')
plt.plot(thetas, post, label='Posterior', c='red')
plt.xlim([0, 1])
plt.xlabel(r'$\theta$', fontsize=14)
plt.ylabel('PDF', fontsize=16)
plt.legend();

结论

蒙特卡罗积分是一种求解积分的数值方法。它的工作原理是在随机点评估一个函数，将这些值相加，然后计算它们的平均值。

Python 中的蒙特卡罗集成

原文：https://towardsdatascience.com/monte-carlo-integration-in-python-a71a209d277e?source=collection_archive---------6-----------------------

一个著名的赌场启发的数据科学，统计和所有科学的把戏。用 Python 怎么做？

图片来源:维基百科(免费)和作者制作的拼贴画

免责声明:本文灵感来源于 佐治亚理工学院在线分析硕士(OMSA) 项目学习资料。我很自豪能继续这个优秀的在线硕士项目。你也可以在这里查看详情。

什么是蒙特卡罗积分？

蒙特卡洛，实际上是世界著名的赌场的名字，位于摩纳哥城邦(也叫公国)的同名区，在世界著名的法国里维埃拉。

事实证明，赌场启发了著名科学家的思维，设计出一种有趣的数学技术，用于解决统计、数值计算和系统模拟中的复杂问题。

图片来源:维基百科

这项技术最早也是最著名的应用之一是在曼哈顿计划期间，当时高浓缩铀的链式反应动力学向科学家们展示了一个难以想象的复杂理论计算。即使是像约翰·冯·诺依曼、斯坦尼斯劳·乌拉姆、尼古拉斯·大都会这样的天才头脑也无法用传统的方式解决它。因此，他们转向了奇妙的随机数世界，让这些概率量驯服原本难以处理的计算。

三位一体测试:(图片来源:维基百科)

令人惊讶的是，这些随机变量可以解决计算问题，这阻碍了稳健的确定性方法。不确定因素实际上赢了。

就像蒙特卡洛游戏世界中的不确定性和随机性规则一样。这就是这个特殊名字的灵感来源。

图片来源: Pixabay

今天，这是一种广泛应用于各个领域的技术——

• 金融工程风险分析，
• 供应链物流，
• 医疗保健研究，药物开发
• 统计学习和建模，
• 计算机图形学，图像处理，游戏设计，
• 大型系统模拟，
• 计算物理，天文学等。

尽管它获得了成功和名声，但它的基本思想看似简单，却很容易演示。在本文中，我们用一组简单的 Python 代码来演示它。

这项技术最早也是最著名的应用之一是在曼哈顿计划中

这个想法

复杂的积分

虽然一般的蒙特卡罗模拟技术在范围上要广泛得多，但我们在这里特别关注蒙特卡罗积分技术。

它只不过是一种计算复杂定积分的数值方法，缺乏封闭形式的解析解。

比方说，我们想计算，

要得到这个不定形式积分的封闭解并不容易，或者完全不可能。但是数值逼近总能给我们定积分作为和。

这是函数的曲线。

黎曼和

在 黎曼求和 的大类下有很多这样的技巧。这个想法只是将曲线下的区域分成小的矩形或梯形块，通过简单的几何计算来近似它们，然后将这些分量相加。

为了简单说明，我展示了这样一个只有 5 个等间距间隔的方案。

对于程序员朋友来说，其实 Scipy 包里有现成的函数可以快速准确的做这个计算。

如果我随机选择呢？

如果我告诉你，我不需要如此一致地选择区间，事实上，我可以完全随机地选择 100%的随机区间来计算相同的积分，会怎么样？

疯话？我选择的样本可能看起来像这样…

或者，这个…

我们没有时间或范围来证明其背后的理论，但可以证明的是通过合理的大量随机采样，我们实际上可以以足够高的精度计算积分！

我们只需选择随机数(在两个极限之间)，在这些点上评估函数，将它们相加，然后用一个已知的因子缩放。我们完了。

好的。我们还在等什么？让我们用一些简单的 Python 代码来演示这种说法。

尽管它获得了成功和名声，但它的基本思想看似简单，却很容易演示。

Python 代码

用简单的平均法代替复杂的数学运算

如果我们试图计算下面形式的积分——任何积分，

我们只是用下面的平均值代替积分的“估计值”,

其中表示 0 到 1 之间的均匀随机数。请注意，我们如何通过简单地将一堆数字相加并取其平均值来代替复杂的积分过程。

在任何现代计算系统、编程语言，甚至像 Excel 这样的商业软件包中，您都可以使用这个统一随机数生成器。查看我关于这个主题的文章，

* [## 如何从头开始生成随机变量(不使用库)

我们通过一个简单的伪随机生成器算法，并显示如何使用它来生成重要的随机…

towardsdatascience.com](/how-to-generate-random-variables-from-scratch-no-library-used-4b71eb3c8dc7)

我们只需选择随机数(在两个极限之间)，在这些点上评估函数，将它们相加，然后用一个已知的因子缩放。我们完了。

该功能

下面是一个 Python 函数，它接受另一个函数作为第一个参数、两个积分限制和一个可选整数来计算参数函数所表示的定积分。

代码看起来可能与上面的等式略有不同(或者您可能在教科书中看到的另一个版本)。这是因为我通过在 10 个时间间隔内分配随机样本来使计算更加精确。

对于我们的具体例子，自变量函数看起来像，

我们可以通过简单地将它传递给monte_carlo_uniform()函数来计算积分，

在这里，如你所见，我们在积分限值 a = 0 和 b = 4 之间随机抽取了 100 个样本。

计算到底有多好？

这个积分不能用解析方法计算。因此，无论如何，我们需要将蒙特卡罗方法的精度与另一种数值积分技术进行对比。我们为此选择了 Scipy integrate.quad()函数。

现在，您可能也在想— 随着采样密度的变化，精度会发生什么变化。这种选择显然会影响计算速度——如果我们选择降低采样密度，我们需要增加更少的数量。

因此，我们针对一系列采样密度模拟了相同的积分，并将结果绘制在黄金标准之上——下图中以水平线表示的 Scipy 函数。

因此，我们在低样品密度阶段观察到一些小的扰动，但随着样品密度的增加，这些扰动会很好地消除。无论如何，与 Scipy 函数返回的值相比，绝对误差非常小，大约为 0.02%。

蒙特卡洛魔术非常有效！

图片来源: Pixabay

速度呢？

但是它和 Scipy 方法一样快吗？好些了吗？更糟？

我们试图通过运行 100 次循环(总共 10，000 次运行)并获得汇总统计数据来找出答案。

在这个特殊的例子中，蒙特卡罗计算的运行速度是 Scipy 积分方法的两倍！

虽然这种速度优势取决于许多因素，但我们可以肯定的是，在计算效率方面，蒙特卡罗技术并不逊色。

我们在低样品密度阶段观察到一些小的扰动，但随着样品密度的增加，这些扰动会很好地消除

冲洗，重复，冲洗，重复…

对于像蒙特卡罗积分这样的概率技术，不用说，数学家和科学家几乎从来不会只运行一次就停下来，而是多次重复计算并取平均值。

这是一个 10，000 次实验的分布图。

如您所见，该图几乎类似于高斯正态分布，这一事实不仅可以用来获得平均值，还可以围绕该结果构建 置信区间 。

[## 置信区间

4±2 的区间置信区间是一个我们相当确定自己真正价值所在的数值范围…

www.mathsisfun.com](https://www.mathsisfun.com/data/confidence-interval.html)

特别适用于高维积分

虽然对于我们的简单说明(以及教学目的)，我们坚持单变量积分，同样的想法可以很容易地扩展到具有多个变量的高维积分。

与基于黎曼和的方法相比，蒙特卡罗方法在这个更高的维度上表现尤为突出。可以以对蒙特卡罗方法更有利的方式优化样本密度，以使其更快，而不损害精度。

在数学上，该方法的收敛速度与维数无关。用机器学习的话来说，蒙特卡罗方法是你在复杂的积分计算中战胜维数灾难的最好朋友。

阅读这篇文章可以获得很好的介绍，

[## 实践中的蒙特卡罗方法(蒙特卡罗积分)

蒙特卡洛方法在实践中，如果你了解和知道最重要的概率和…

www.scratchapixel.com](https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods-in-practice/monte-carlo-integration)

与基于黎曼和的方法相比，蒙特卡罗方法在这个更高的维度上表现尤为突出。

摘要

我们介绍了蒙特卡罗积分的概念，并说明了它与传统的数值积分方法的不同。我们还展示了一组简单的 Python 代码来评估一维函数，并评估这些技术的准确性和速度。

更广泛的一类蒙特卡罗模拟技术更令人兴奋，并且在与人工智能、数据科学和统计建模相关的领域中以无处不在的方式使用。

例如，DeepMind 的著名 Alpha Go 程序使用了蒙特卡罗搜索技术，在围棋的高维空间中计算效率高。实践中可以找到许多这样的例子。

[## 蒙特卡洛树搜索简介:deep mind alpha go 背后的游戏规则改变算法

一场五局三胜的系列赛，100 万美元的奖金——一场高赌注的枪战。2016 年 3 月 9 日至 15 日之间…

www.analyticsvidhya.com](https://www.analyticsvidhya.com/blog/2019/01/monte-carlo-tree-search-introduction-algorithm-deepmind-alphago/)

如果你喜欢它…

如果您喜欢这篇文章，您可能也会喜欢我关于类似主题的其他文章，

[## 如何从头开始生成随机变量(不使用库)

我们通过一个简单的伪随机生成器算法，并显示如何使用它来生成重要的随机…

towardsdatascience.com](/how-to-generate-random-variables-from-scratch-no-library-used-4b71eb3c8dc7) [## 数学编程——数据科学进步的关键习惯

我们展示了建立数学编程习惯的一小步，这是一个崭露头角的人的关键技能

towardsdatascience.com](/mathematical-programming-a-key-habit-to-built-up-for-advancing-in-data-science-c6d5c29533be) [## 用 Python 实现布朗运动

我们展示了如何模拟布朗运动，在广泛的应用中使用的最著名的随机过程，使用…

towardsdatascience.com](/brownian-motion-with-python-9083ebc46ff0)

喜欢这篇文章吗？成为 中等成员 继续 无限制学习 。如果你使用下面的链接，我会收到你的一部分会员费， 而不需要你额外付费 。

[## 通过我的推荐链接加入媒体

作为一个媒体会员，你的会员费的一部分会给你阅读的作家，你可以完全接触到每一个故事…

medium.com](https://medium.com/@tirthajyoti/membership)*

蒙特卡罗马尔可夫链

原文：https://towardsdatascience.com/monte-carlo-markov-chain-89cb7e844c75?source=collection_archive---------16-----------------------

希瑟·吉尔在 Unsplash 上拍摄的照片

一个蒙特卡罗马尔可夫链 ( MCMC )是一个描述一系列可能事件的模型，其中每个事件的概率只取决于前一个事件达到的状态。MCMC 有广泛的应用，其中最常见的是概率分布的近似。

让我们看一个蒙特卡罗马尔可夫链的例子。假设我们想确定晴天和雨天的概率。

我们得到了以下条件概率:

• 如果今天下雨，明天有 50%的可能是晴天。
• 如果今天下雨，明天有 50%的可能会下雨。
• 如果今天是晴天，明天有 90%的可能是晴天。
• 如果今天是晴天，明天有 10%的可能会下雨。

假设我们从阳光充足的州开始。然后我们进行蒙特卡罗模拟。也就是说，我们生成一个 0 到 1 之间的随机数，如果刚好在 0.9 以下，明天就是晴天，否则就是阴雨。我们又做了一次蒙特卡洛模拟，这一次，明天会下雨。我们重复该过程进行 n 次迭代。

下面的序列被称为马尔可夫链。

我们数出晴天的天数，除以总天数，来确定天气晴朗的概率。如果马尔可夫链足够长，即使初始状态可能不同，我们也会得到相同的边际概率。这种情况下，天晴的概率是 83.3%，下雨的概率是 16.7%。

让我们看看如何用 Python 实现蒙特卡罗马尔可夫链。

我们从导入以下库开始:

import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
sns.set()

我们可以使用状态机和相应的矩阵来表达上一个例子中的条件概率。

我们做一些线性代数。

T = np.array([[0.9, 0.1],[0.5, 0.5]])p = np.random.uniform(low=0, high=1, size=2)
p = p/np.sum(p)
q=np.zeros((100,2))for i in np.arange(0,100):
    q[i, :] = np.dot(p,np.linalg.matrix_power(T, i))

最后，我们绘制结果。

plt.plot(q)
plt.xlabel('i')
plt.legend(('S', 'R'))

正如我们所见，天气晴朗的概率约为 0.833。同样，多雨的概率也趋向于 0.167。

贝叶斯公式

很多时候，我们想知道某个事件发生的概率，假设另一个事件已经发生。这可以象征性地表示为 p(B|A) 。如果两个事件不是独立的，那么这两个事件发生的概率由下面的公式表示。

例如，假设我们从一副标准的 52 张牌中抽取两张牌。这副牌中一半是红色的，一半是黑色的。这些事件不是独立的，因为第二次抽签的概率取决于第一次。

P(A) = P(第一次抽黑牌)= 25/52 = 0.5

P(B|A) = P(第二次抽黑牌|第一次抽黑牌)= 25/51 = 0.49

利用这些信息，我们可以计算连续抽两张黑牌的概率，如下所示:

现在，让我们假设，我们想要开发一个垃圾邮件过滤器，它将根据某个单词的出现与否来将一封电子邮件分类为垃圾邮件。例如，如果一封电子邮件包含单词 【伟哥】 ，我们将其归类为垃圾邮件。另一方面，如果一封电子邮件包含单词 money，那么它有 80%的可能是垃圾邮件。

根据贝叶斯定理，给定包含给定单词的电子邮件是垃圾邮件的概率为:

我们可能知道一封电子邮件是垃圾邮件的概率，以及一个单词包含在被分类为垃圾邮件的电子邮件中的概率。然而，我们不知道给定的单词在电子邮件中出现的概率。这就是 Metropolis-Hastings 算法发挥作用的地方。

大都会黑斯廷斯算法

Metropolis-Hasting 算法使我们能够在不知道归一化常数的情况下确定后验概率。在高层次上，Metropolis-Hasting 算法的工作方式如下:

验收标准只考虑目标分布的比率，因此分母抵消。

对于视觉学习者来说，让我们用一个例子来说明这个算法是如何工作的。

我们首先为θ选择一个随机的初始值。

然后，我们提出一个新的θ值。

我们计算当前θ值的 PDF 和建议θ值的 PDF 之比。

如果 rho 小于 1，那么我们以概率 p 将 theta 设置为新值。我们通过将 rho 与从均匀分布中抽取的样本进行比较。如果 rho 大于 u 则我们接受建议值，否则，我们拒绝它。

我们尝试不同的θ值。

如果ρ大于 1，它将总是大于或等于从均匀分布中抽取的样本。因此，我们接受θ新值的建议。

**

我们重复这个过程 n 次 次迭代。

因为当目标分布大于当前位置时，我们自动接受提议的移动，所以θ将倾向于在目标分布更密集的地方。然而，如果我们只接受比当前位置大的值，我们会卡在其中一个峰值上。因此，我们偶尔会接受向低密度区域转移。这样，θ将被期望以这样的方式反弹，以接近后验分布的密度。

这些步骤实际上是一个马尔可夫链。

让我们看一下如何用 Python 实现 Metropolis-Hasting 算法，但是首先，这里有一个不同类型的发行版的快速复习。

正态分布

在自然界中，随机现象(如智商、身高)往往遵循正态分布。正态分布有两个参数μ和适马。改变μ会移动钟形曲线，而改变适马会改变钟形曲线的宽度。

https://commons . wikimedia . org/wiki/File:Normal _ Distribution _ pdf . SVG

贝塔分布

像正态分布一样，贝塔分布有两个参数。然而，与正态分布不同，β分布的形状会根据其参数α和β的值而显著变化。

https://commons . wikimedia . org/wiki/File:Beta _ distribution _ pdf . SVG

二项分布

与以高度为定义域的正态分布不同，二项式分布的定义域始终是离散事件的数量。

https://commons . wikimedia . org/wiki/File:Binomial _ distribution _ PMF _ sl . SVG

现在我们已经熟悉了这些概念，我们准备深入研究代码。我们从初始化超参数开始。

*n = 100
h = 59
a = 10
b = 10
sigma = 0.3
theta = 0.1
niters = 10000
thetas = np.linspace(0, 1, 200)
samples = np.zeros(niters+1)
samples[0] = theta*

接下来，我们定义一个函数，对于给定的θ值，该函数将返回似然和先验的乘积。

*def prob(theta):
    if theta < 0 or theta > 1:
        return 0
    else:
        prior = st.beta(a, b).pdf(theta)
        likelihood = st.binom(n, theta).pmf(h)
        return likelihood * prior*

我们逐步执行算法，根据前面描述的条件更新θ值。

*for i in range(niters):
    theta_p = theta + st.norm(0, sigma).rvs()
    rho = min(1, prob(theta_p) / prob(theta))
    u = np.random.uniform()
    if u < rho:
        # Accept proposal
        theta = theta_p
    else:
        # Reject proposal
        pass
    samples[i+1] = theta*

我们定义了可能性，以及先验和后验概率分布。

*prior = st.beta(a, b).pdf(thetas)
post = st.beta(h+a, n-h+b).pdf(thetas)
likelihood = st.binom(n, thetas).pmf(h)*

我们将使用 Metropolis-Hastings 算法获得的后验分布可视化。

*plt.figure(figsize=(12, 9))
plt.hist(samples[len(samples)//2:], 40, histtype='step', normed=True, linewidth=1, label='Predicted Posterior');
plt.plot(thetas, n*likelihood, label='Likelihood', c='green')
plt.plot(thetas, prior, label='Prior', c='blue')
plt.plot(thetas, post, c='red', linestyle='--', alpha=0.5, label='True Posterior')
plt.xlim([0,1]);
plt.legend(loc='best');*

正如我们所看到的，Metropolis Hasting 方法在逼近实际后验分布方面做得很好。

结论

蒙特卡洛马尔可夫链是从一组概率分布中提取的一系列事件，这些概率分布可以用来近似另一个分布。Metropolis-Hasting 算法在我们知道似然和先验，但不知道归一化常数的情况下，利用蒙特卡罗马尔可夫链来近似后验分布。

基于 Python 语言伊辛模型的蒙特卡罗方法在 2D 二元合金中的应用

原文：https://towardsdatascience.com/monte-carlo-method-applied-on-a-2d-binary-alloy-using-an-ising-model-on-python-70afa03b172b?source=collection_archive---------6-----------------------

来源:2018 年影响因子，高级科学新闻

使用 Python 中的蒙特卡罗方法，使用 Ising 模型评估材料的属性

诸如一片金属或玻璃样品的材料具有诸如导电性、磁化率、导热性和许多其他的物理性质。对于工业应用，必须了解这些特性，以便选择满足称重传感器要求的最佳材料或合金。例如，汽车的车身必须重量轻、承重大、耐腐蚀。为了满足所有这些需求，使用的三种主要材料是:金属、聚合物和铝合金。但是你如何确定一个物体的属性呢？一个由几十亿几十亿个微观状态组成的物质，怎么可能计算出它的性质？

严格来说，我们需要知道每个微观状态的能量，以便计算一种材料的性质。在本文中，我将重点介绍四个物理量:能量、磁化率、磁化强度和热容量。这些量是相互联系的，并且很好地定义了一个系统，因为它们是基本的物理属性，尽管能量是我们下面定义的所有关系的基础。但是还有一个问题:如何模拟一种材料的行为来计算它的性质？

由于我们受到计算机能力的限制，我们无法模拟硬件的完整行为，因为这将花费我们太多的时间，并导致我们提出一个新的问题:系统的大小应该如何接近现实？

为了回答这些问题，我将考虑 2D 二元合金作为我要研究的材料。对于我感兴趣的所有属性，我将比较不同系统大小的理论值和实验值(由模拟给出的值)。

但首先，我们将看看数学和物理学，它们为我们想要评估的不同属性建立了关系。

发现具有能量 E 的处于特定微观状态的系统的概率由正则分布给出:

𝛃 = 1/kT，k 定义为玻尔兹曼常数，t 定义为温度。z 被称为再分配函数，由下式给出:

有了这两个方程，我们可以测量任何实验参数，但一般来说，求和是不能进行的。幸运的是，发现处于能量 E 的微观状态的系统的概率 P(E)是能量 E 的精确函数。因此，我们可以得到在计算物理性质时要考虑的微观状态数的近似值。在本文中，我将展示系统的大小对物理性质计算的准确性的影响。

我之前说过，我们只对能量，磁化率，磁化强度，热容量感兴趣。有多种方法可以模拟这些性质，如密度泛函理论(最准确的一种)、机器学习或蒙特卡罗方法。

在这篇文章中，我决定建立一个 H=0 的伊辛 2D 模型的蒙特卡罗模拟。使用这个模型，我能够计算 L=4、8、16 和 32 的 L×L 自旋系统的自旋磁化强度绝对值的期望值，作为温度的函数(伊辛模型是自旋在图上的表示)。蒙特卡洛方法基于随机采样的重复(将旋转从-1 改变为 1，反之亦然)来获得新的能量值。如果总能量的值减少，我们接受新的情况，因为这意味着我们越来越接近能量最小的平衡状态。
因为用无量纲单位工作更方便，所以我用 J 单位计算了所有的能量。所有温度均以 J/k 为单位测量。

总系统能量的表达式为:

其中𝛔是自旋，可以等于-1 或 1。此外，我计算了每个温度下每个𝜒站点的磁化率。

这里 M 是系统的总磁化强度，T 是温度。每个位置的磁化强度由下式给出:

最后，每次旋转的热容 C 如下:

用 Tc 定义的临界温度 ln(1 + √2) = (2J)/(k*Tc)。可以看出，热容量在临界温度下呈对数发散。这种发散在有限大小的系统中是不可见的。此外，我没有磁化率的解析解，所以我没有将模拟结果与精确结果进行比较。

现在，数学已经完成，我们可以进入计算部分。

使用 Python，我建立了 H=0 的 2D 伊辛模型的蒙特卡罗模拟，并包含了每个物理属性的不确定性。

**def** init(L):
        state = 2 * np.random.randint(2, size=(L,L)) - 1
        **return** state

**def** E_dimensionless(config,L):
    total_energy = 0
    **for** i **in** range(len(config)):
        **for** j **in** range(len(config)):
            S = config[i,j]
            nb = config[(i+1)%L, j] + config[i, (j+1)%L] + config[(i-1)%L, j] + config[i, (j-1)%L]
            total_energy += -nb * S
    **return** (total_energy/4)

**def** magnetization(config):
    Mag = np.sum(config)

    **return** Mag

我首先初始化了我的系统，它是 LxL 自旋的平方，可以取值-1 或 1。我创建了一个计算能量的函数和另一个计算磁化强度的函数。

**def** MC_step(config, beta):
    *'''Monte Carlo move using Metropolis algorithm '''*
    L = len(config)
    **for** i **in** range(L):
        **for** j **in** range(L):
            a = np.random.randint(0, L) *# looping over i & j therefore use a & b*
            b = np.random.randint(0, L)
            sigma =  config[a, b]
            neighbors = config[(a+1)%L, b] + config[a, (b+1)%L] + config[(a-1)%L, b] + config[a, (b-1)%L]
            del_E = 2*sigma*neighbors
            **if** del_E < 0:
                sigma *= -1
            **elif** rand() < np.exp(-del_E*beta):
                sigma *= -1
            config[a, b] = sigma
    **return** config

然后我用 Metropolis 算法进行蒙特卡罗模拟。基于一个随机的初始配置，我执行了一个随机的自旋变化，并重新计算了能量。如果系统的能量减少，当我们接近平衡时，新的构型被接受。否则，如果新系统具有更高的能量，则配置以 exp()βδE 给出的概率被接受。

**def** calcul_energy_mag_C_X(config, L, eqSteps, err_runs):

    *# L is the length of the lattice*

    nt      = 100         *#  number of temperature points*
    mcSteps = 1000

    T_c = 2/math.log(1 + math.sqrt(2))

    *# the number of MC sweeps for equilibrium should be at least equal to the number of MC sweeps for equilibrium*

    *# initialization of all variables*
    T = np.linspace(1., 7., nt); 
    E,M,C,X = np.zeros(nt), np.zeros(nt), np.zeros(nt), np.zeros(nt)
    C_theoric, M_theoric = np.zeros(nt), np.zeros(nt)
    delta_E,delta_M, delta_C, delta_X = np.zeros(nt), np.zeros(nt), np.zeros(nt), np.zeros(nt) n1 = 1.0/(mcSteps*L*L)
    n2 = 1.0/(mcSteps*mcSteps*L*L)    *# n1 and n2 will be use to compute the mean value and the # by sites*
        *# of E and E^2*

    Energies = []
    Magnetizations = []
    SpecificHeats = []
    Susceptibilities = []
    delEnergies = []
    delMagnetizations = []
    delSpecificHeats = []
    delSusceptibilities = []

    **for** t **in** range(nt):
        *# initialize total energy and mag*
        beta = 1./T[t]
        *# evolve the system to equilibrium*
        **for** i **in** range(eqSteps):
            MC_step(config, beta)
        *# list of ten macroscopic properties*
        Ez = [], Cz = [], Mz = [], Xz = [] 

        **for** j **in** range(err_runs):
            E = E_squared = M = M_squared = 0
            **for** i **in** range(mcSteps):
                MC_step(config, beta)           
                energy = E_dimensionless(config,L) *# calculate the energy at time stamp*
                mag = abs(magnetization(config)) *# calculate the abs total mag. at time stamp*

                *# sum up total energy and mag after each time steps*

                E += energy
                E_squared += energy**2
                M += mag
                M_squared += mag**2

            *# mean (divide by total time steps)*

            E_mean = E/mcSteps
            E_squared_mean = E_squared/mcSteps
            M_mean = M/mcSteps
            M_squared_mean = M_squared/mcSteps

            *# calculate macroscopic properties (divide by # sites) and append*

            Energy = E_mean/(L**2)
            SpecificHeat = beta**2 * (E_squared_mean - E_mean**2)/L**2
            Magnetization = M_mean/L**2
            Susceptibility = beta * (M_squared_mean - M_mean**2)/(L**2)

            Ez.append(Energy); Cz.append(SpecificHeat); Mz.append(Magnetization); Xz.append(Susceptibility);

        Energy = np.mean(Ez)
        Energies.append(Energy)
        delEnergy = np.std(Ez)
        delEnergies.append(float(delEnergy))

        Magnetization = np.mean(Mz)
        Magnetizations.append(Magnetization)
        delMagnetization = np.std(Mz)
        delMagnetizations.append(delMagnetization)

        SpecificHeat = np.mean(Cz)
        SpecificHeats.append(SpecificHeat)
        delSpecificHeat = np.std(Cz)
        delSpecificHeats.append(delSpecificHeat)

        Susceptibility = np.mean(Xz)
        delSusceptibility = np.std(Xz)        
        Susceptibilities.append(Susceptibility)
        delSusceptibilities.append(delSusceptibility)

        **if** T[t] - T_c >= 0:
            C_theoric[t] = 0
        **else**:
            M_theoric[t] = pow(1 - pow(np.sinh(2*beta), -4),1/8)

        coeff = math.log(1 + math.sqrt(2))
        **if** T[t] - T_c >= 0:
            C_theoric[t] = 0
        **else**: 
            C_theoric[t] = (2.0/np.pi) * (coeff**2) * (-math.log(1-T[t]/T_c) + math.log(1.0/coeff) - (1 + np.pi/4)) 

    **return** (T,Energies,Magnetizations,SpecificHeats,Susceptibilities, delEnergies, delMagnetizations,M_theoric, 
            C_theoric, delSpecificHeats, delSusceptibilities)

对于每个系统尺寸(4x4、8x8、16x16、32x32)，我计算了能量、磁化强度、比热和磁化率。每个图形对应一个属性。磁化强度和比热的黑线是精确解，并且对于每种性质，已经考虑了误差估计。

执行上面的脚本后，我得到了以下结果:

评论

模拟说明了有限规模的概念。使用已知的𝑇𝑐和临界指数的精确结果，我们可以检查大系统的模拟结果是否实际上向热力学极限行为收敛。临界现象理论最基本的概念是关联长度，它是一个典型系统长度尺度的量度。相关长度可以根据相关函数来定义，在伊辛模型的情况下，相关函数由下式给出。相关函数在长距离上呈指数下降，由所谓的 Ornstein-Zernicke 形式给出:

𝐶(𝑟⃗)=𝑒xp(−𝑟/𝜉)/r^(d-2/2)

其中𝜉是相关长度。这个表达式只有当𝑟>t13】𝜉时才有效。相关长度𝜉对应于系统中有序畴的典型尺寸。在𝑇𝑐上方的伊辛模型中，自旋向上和自旋向下的有序畴的数量平均相等(材料失去了所有，它们的典型尺寸对应于关联长度)。在 Tc 以下，相关长度对应于相反方向自旋的有序背景中自旋畴的典型尺寸。

当接近临界点(连续相变)时，相关长度在热力学极限内根据幂律发散:**【𝜉∼𝑡^(−𝜈】**其中 t 是测量离临界点距离的降低温度:𝑡=|𝑇−𝑇𝑐|和𝜈是临界指数的一个例子。另一个重要的临界指数是决定临界点处有序开始的指数，即磁化强度。序参量在𝑇𝑐以上为零，在𝑇𝑐以下出现为: < 𝑚 > =(𝑇𝑐−𝑇)^𝛽 。

对于热容量，关系由下面的公式给出 𝐶=𝑡^(−𝛼).当 t 接近𝑇𝑐.时，热容量出现一个奇点

对于易感性，关系如下: 𝜒=𝑡^(−𝛾).随着 t 接近𝑇𝑐.，𝜒似乎出现了一个奇点因此，当 t 接近𝑇𝑐.时，系统对磁场 h 变得无限敏感

不同系统的𝜈,𝛼,𝛽,𝛾指数集属于普适类。属于同一普适类的系统具有相同的指数。普遍性类别不依赖于与系统成分及其相互作用相关的微观细节(只要相互作用是短程的；长程相互作用可以改变普适类)，只是在系统的维数和序参量的性质上。在我们的例子中，系统的大小因此是 l。这就是为什么磁化曲线、热容和温度磁化率随着系统大小的变化而彼此不同。有限尺寸标度的基础是，当相关长度𝜒变得与系统长度 l 相当时，就会出现与无限尺寸临界行为的偏差。这些有限尺寸偏差影响其他量的方式可以通过使用相关长度作为变量来表达它们的温度依赖性来进行研究。由此， 𝑡∼𝜉^(−1/𝜈)

例如，在磁化强度的渐近幂律形式中**<𝑚>=𝑇c^𝛽𝑡𝛽=𝑇c^𝛽𝜉(−𝛽/𝜈).**对于给定的系统尺寸 l，磁化强度的最大值应为 **< 𝑚 > =𝑇c^𝛽 L^(−𝛽/𝜈).**磁化强度向零的收敛取决于系统 l 的维数大小。在 2D 伊辛模型中，我们还具有以下关系 𝛽/𝜈=1/8 。

基于下面的表达式𝜒=𝑡^(−𝛾)，在磁化率的渐近幂律形式中，t 的替换给出了 **𝜒=𝜉^(𝛾/𝜈).**给定系统尺寸 l 的磁化率最大值应为 **𝜒=L^(𝛾/𝜈).**通过该表达式，我们可以看到𝜒的峰值随 l 快速增长，并解释了不同系统尺寸的磁化率曲线的差异。对于二维伊辛模型，我们有 𝛾/𝜈=7/4

关于热容量，关系非常相似，它是下面的一个 𝐶=𝜉^(𝛼/𝜈) 并且对于给定的系统尺寸 l，热容量的最大值应该是 **𝐶=𝜉^(𝛼/𝜈).**因此，就磁化率而言，C 的峰值随 L 快速增长，这解释了不同系统尺寸的热容图中剖面的差异。

从更定性的角度来看，我们注意到能量分布的𝑇=𝑇𝑐有一个拐点，但是磁化强度、热容量和磁化率从这个温度下降到零。这意味着系统行为的改变，这对应于临界温度的定义。事实上，临界温度是我们情况下的居里温度，居里温度是材料失去所有永久磁性的温度，这预示着材料行为的变化。永久磁性是由磁矩的排列引起的。因此，在𝑇𝑐上方，平均有相等数量的有序的向上和向下自旋的畴，这可以通过磁矩、比热和磁化率向零的收敛来说明。

因此，结果的准确性和系统的大小之间的关系取决于我们正在检查的相关性。但是系统规模越大，计算成本越高。因此，由于时间成本非常重要，我们需要考虑权衡:时间/成本/精度。这就是为什么我想知道每个属性对系统大小的依赖性。我所理解的是，对于能量，我们可以用最小的系统规模来执行我们的模拟，并且结果将是可用的，因为能量对系统规模的依赖是指数的。对于系统规模越大、结果越精确的其他属性，这一点不再适用。尽管如此，我们仍然可以通过增加 32×32 系统大小达到平衡所需的步骤数来改进我们的结果。

如果你喜欢阅读这样的故事，并想支持我成为一名作家，考虑注册成为一名灵媒成员。每月 5 美元，你可以无限制地阅读媒体上的故事。如果你使用[我的链接](http://If you enjoy reading stories like these and want to support me as a writer consider signing up to become a Medium member. It’s $5 a month, giving you unlimited access to stories on Medium. If you sign up using my link, I will earn a small commission and you will still pay $5. Thank you !! https://medium.com/@jonathan_leban/membership)注册，我将赚取一小笔佣金，你仍需支付 5 美元。谢谢大家！！

[## 通过我的推荐链接加入媒体-乔纳森·莱班

阅读乔纳森·莱班的每一个故事(以及媒体上成千上万的其他作家)。您的会员费直接支持…

medium.com](https://medium.com/@jonathan_leban/membership)

PS:我现在是柏克莱大学的工程硕士，如果你想讨论这个话题，请随时联系我。 这里的 是我的邮箱。

蒙特卡罗方法

原文：https://towardsdatascience.com/monte-carlo-methods-9b289f030c2e?source=collection_archive---------13-----------------------

深度强化学习讲解— 13

探索-解释困境

在*深度强化学习讲解系列的这篇新帖中，我们将介绍另一种经典的强化学习方法来估计一个策略π的值。最直接的方法是运行几集，收集数百条轨迹，然后计算每个州的平均值。这种估计价值函数的方法叫做蒙特卡罗预测 (MC)。*

在本帖中，我们还将介绍如何评估最优策略和勘探-开采困境**。**

本出版物的西班牙语版本

** [## 5.蒙特卡洛政治评估

请访问第 5 页的自由介绍

medium.com](https://medium.com/aprendizaje-por-refuerzo/5-evaluaci%C3%B3n-de-pol%C3%ADticas-con-monte-carlo-a6d70d1db7d4)

蒙特卡洛与动态规划

在本系列的第 1 部分中，我们介绍了一种 MDP 的解决方案，称为动态编程，由理查德·贝尔曼首创。记住，贝尔曼方程允许我们递归地定义值函数，并且可以用值迭代算法求解。总而言之，动态编程为强化学习提供了基础，但我们需要在每次迭代中循环遍历所有状态(它们的大小可以呈指数级增长，状态空间可以非常大，也可以无限大)。动态编程还需要环境的模型，特别是知道状态转移概率p(s′，r|s，a) 。

相比之下，蒙特卡罗方法都是从经验中学习。任何期望值都可以通过样本均值来近似——换句话说，我们需要做的就是播放一堆剧集，收集回报，然后取平均值。蒙特卡罗方法实际上是基本算法的一组替代方案。这些仅适用于偶发任务，当代理遇到终止状态时，交互停止。也就是说，我们假设体验被分成几集，并且不管选择什么动作，所有的集最终都会终止。

重要的是要注意，蒙特卡罗方法只给我们遇到的状态和动作一个值，如果我们从未遇到一个状态，它的值是未知的。

蒙特卡罗方法

这篇文章将提供一个用于强化学习的蒙特卡罗的实用方法。关于这些方法的更正式的解释，我邀请读者阅读理查德·萨顿和安德鲁·巴尔托的教科书 强化学习:简介 的第五章。

回想一下，最优策略**【π∫规定了，对于每个环境状态*，代理应该如何选择一个行动来实现其最大化报酬的目标。我们还了解到，代理可以通过首先估计最优行动值函数q∑来构建对最优策略的搜索；那么一旦q∫已知，就很快得到π∫。***

代理开始采取一个基本策略，像等概率随机策略，一个随机策略，代理从每个状态中随机选择一组可用的动作，每个动作以等概率被选择。代理使用当前策略π将与环境交互以收集一些情节，然后合并结果以得出更好的策略。

方法是用一个我们称之为 Q 表的表来估算动作值函数。蒙特卡罗方法中的这个核心表为每个状态提供了一行，为每个动作提供了一列。对应于状态 s 和动作 a 的条目表示为 Q ( s ， a )。

我们称之为预测问题:给定一个策略，代理如何估计该策略的价值函数？。我们将预测问题的蒙特卡罗(MC)方法称为 MC 预测方法。

我们将把我们的解释集中于动作值函数 Q，但是 MC 也可以用于估计状态值函数 v。

在用于 MC 预测的算法中，我们通过收集策略的许多片段开始。然后，我们注意到 Q 表中的每个条目对应于一个特定的状态和动作。为了填充 Q-table 的一个条目，我们使用当代理处于那个状态并选择动作时跟随的返回。

我们将一集中状态的每次出现定义为对该状态-动作对的访问**。在一集里，一个状态-动作对可能被访问不止一次。这导致我们有两个版本的 MC 预测算法:**

• 每次访问 MC 预测:在所有事件中，对每个状态-动作对的所有访问后的回报进行平均。
• 初访 MC 预测:每一集，我们只考虑初访状态-动作对。

首次访问方法和每次访问都被认为是保证收敛到真正的行动价值函数。

在本帖中，我们将对 OpenAI 健身房的工作示例进行首次访问:21 点环境。但是实际上，在这个例子中，首次访问和每次访问 MC 返回相同的结果。注意，同一个状态不会在一个情节中重复出现，所以在首次访问和每次访问 MC 方法之间没有区别。首次就诊 MC 预测的伪代码如下:

在这个伪代码中，变量 num _ 剧集表示代理收集的剧集数，有三个相关的表格:

• Q :一个 Q 表，每个状态一行，每个动作一列。
• N :记录我们对每个状态-动作对的第一次访问次数的表格。
• returns _ sum:一个表，记录第一次访问每个状态-动作对后获得的奖励的总和。

**在每一集之后， N 和返回 _ sum 表被更新以存储该集包含的信息。在收集了所有剧集之后，完成对Q(Q 表)的最终估计。

案例研究:21 点

为了更好地理解蒙特卡罗在实践中是如何工作的，让我们使用 OpenAI 的健身房环境Blackjack-v0对 21 点游戏执行一步一步的编码。

二十一点规则

首先，让我们定义游戏的规则和条件:

• 流行的赌场纸牌游戏 21 点的目标是获得数值总和尽可能大而不超过 21 的纸牌。超过 21 导致破产，而双方都有 21 导致平局。
• 我们考虑一个简化的版本，其中每个玩家独立地与庄家竞争。这意味着我们将只和庄家对打，没有其他玩家参与。
• 数字卡的价值是按面值计算的。牌 J、K 和 Q 的值是 10。ace 的值可以是 1 或 11，取决于玩家的选择。
• 双方，玩家和庄家，都有两张牌。玩家的两张牌面朝上，而庄家的一张牌面朝上。在玩家看到他们的牌和庄家的第一张牌后，玩家可以选择击打或站立，直到他对自己的总和感到满意，之后他将站立。
• 然后，发牌者展示他们的第二张牌——如果总数小于 17，他们将继续抽牌，直到达到 17，之后他们将站着。

开放式健身房:21 点环境

我们将使用 OpenAI 的健身房环境[Blackjack-v0](https://github.com/openai/gym/blob/master/gym/envs/toy_text/blackjack.py)。每个状态是一个三元组:

• 玩家目前的总和∈{0，1，…，31}
• 庄家的面朝上的牌∈{1，…，10}
• 无论玩家是否有可用的 a。

代理有两个潜在的操作:

• 棍(动作0):不再取牌(也称“立”或“留”)。
• 击(动作1):从庄家那里再拿一张牌。

21 点的每一局都是一集。赢、输和抽牌的奖励分别为+1、1 和 0。一个游戏内所有奖励都是 0，我们不打折(gamma = 1)；因此，这些最终的回报也是回报。

这篇文章的全部代码可以在 GitHub 上找到，而可以通过这个链接作为一个谷歌笔记本来运行。

我们从导入必要的包开始:

**import sys
import gym
import numpy as np
from collections import defaultdictenv = gym.make('Blackjack-v0')
print(env.observation_space)
print(env.action_space)**

我们可以看到，在环境中有 704 种不同的状态，对应于 32 乘以 11 乘以 2，并且有两种可能的动作对应于选择坚持或击打:

**Tuple(Discrete(32), Discrete(11), Discrete(2)) 
Discrete(2)**

为了在每集的每个时间步骤中查看示例游戏，代理与环境进行交互，我们可以运行以下代码(多次):

**state = env.reset()
while True:
       print(state)
       action = env.action_space.sample()
       state, reward, done, info = env.step(action)
       if done:
          if reward > 0: 
             print('Reward: ', reward, '(Player won)\n')
          else: 
             print('Reward: ', reward, '(Player lost)\n')
          break**

示例游戏的一个例子是:

**(21, 4, True) 
(13, 4, False) 
(14, 4, False) 
Reward:  -1.0**

经纪人输掉比赛的地方。我建议读者运行这段代码几次，看看不同的剧本。

一种简单的 MC 预测方法实现

考虑这样一个策略，玩家只根据她/他当前的分数来决定一个动作。例如，如果我们拥有的卡的总数是 18 或更少，我们认为如果我们要求一辆新的汽车可能是没问题的。我们有 75%的可能性做到这一点。如果卡的总数大于 18，我们认为接受一张新卡太危险了，我们不会以 75%的概率这样做。

总之，如果总和大于 18，我们以 75%的概率选择动作棒;并且，如果总和等于或小于 18，我们选择动作以 75%的概率命中。

请注意，在第一种策略中，我们忽略了状态中的一些信息，例如，庄家的面朝上的牌或我们是否有可用的 a。这是为了简化示例中一段重点代码的解释。以下函数generate_episode使用该策略对情节进行采样:

**def generate_episode(env):
    episode = []
    state = env.reset() 
    while True:
        probs = [0.75, 0.25] if state[0] > 18 else [0.25, 0.75]
        action = np.random.choice(np.arange(2), p=probs)
        next_state, reward, done, info = env.step(action)
        episode.append((state, action, reward))
        state = next_state
        if done:
           break
        return episode**

这段代码返回一个episode(使用策略 π 决定)作为元组的(状态、动作、回报)元组列表。episode[i][0]、episode[i][1]、episode[i][2]分别对应时步𝑖的状态、时步 i 的动作、时步 𝑖+1 的奖励。

为了玩(10 个游戏)，我们可以如下执行前面的功能:

**for i in range(10):
    print(generate_episode(env))**

输出将是(在我的运行中):

**[((19, 3, True), 1, 0.0), ((19, 3, False), 0, 1.0)] 
[((14, 5, False), 1, -1.0)] 
[((15, 10, False), 0, -1.0)] 
[((15, 10, False), 1, 0.0), ((17, 10, False), 1, -1.0)] 
[((13, 2, False), 1, -1.0)] 
[((13, 4, False), 1, 0.0), ((14, 4, False), 0, -1.0)] 
[((11, 10, False), 1, 0.0), ((21, 10, False), 0, 0.0)] 
[((15, 10, False), 1, -1.0)] 
[((9, 9, False), 1, 0.0), ((16, 9, False), 0, 1.0)] 
[((13, 7, False), 1, 0.0), ((18, 7, False), 1, 0.0), ((21, 7, False), 1, -1.0)]**

奖励只有在游戏结束时才会收到，如果我们赢了就是1.0，如果我们输了就是-1.0。我们看到有时我们赢了，有时我们输了。

让我们开始根据前面的伪代码编写一个代码。首先，我们需要初始化字典N、return_sum和Q:

**N = defaultdict(lambda: np.zeros(env.action_space.n))
returns_sum = defaultdict(lambda: np.zeros(env.action_space.n))
Q = defaultdict(lambda: np.zeros(env.action_space.n))**

接下来，该算法在使用所提供的使用该策略的函数generate_episode生成的剧集上循环。每一集都将是状态、动作和奖励元组的列表。然后，我们使用 zip 命令将状态、动作和奖励分成不同的量:

**episode = generate_episode(env)
states, actions, rewards = zip(*episode)**

让我们回到 de 伪代码，看看我们在时间步长上循环，查看每个时间步长对应的状态动作对。如果这是我们第一次访问该对，我们将表N的相应位置增加 1，并将此时步的返回添加到表return_sum的相应条目中。

请记住，在这个 21 点的例子中，第一次访问和每次访问 MC 预测是等效的，因此我们将对每个时间步长进行更新。然后，一旦我们有了return_sum和N的相应更新值，我们就可以用它们来更新我们的Q估算表。这部分的代码如下:

**discounts = np.array([gamma**i for i in range(len(rewards)+1)])for i, state in enumerate(states):
    returns_sum[state][actions[i]] += 
               sum(rewards[i:]*discounts[:-(1+i)])
    N[state][actions[i]] += 1.0
    Q[state][actions[i]] = returns_sum[state][actions[i]] 
               / N[state][actions[i]]**

我们的 MC 预测方法的完整代码如下:

**def mc_prediction(env, num_episodes, generate_episode, gamma=1.0):

   returns_sum = defaultdict(lambda: np.zeros(env.action_space.n))
   N = defaultdict(lambda: np.zeros(env.action_space.n))
   Q = defaultdict(lambda: np.zeros(env.action_space.n))for episode in range(1, num_episodes+1):
      episode = generate_episode(env)
      states, actions, rewards = zip(*episode)discounts = np.array([gamma**i for i in 
                           range(len(rewards)+1)])for i, state in enumerate(states):
             returns_sum[state][actions[i]] += 
                        sum(rewards[i:]*discounts[:-(1+i)])
             N[state][actions[i]] += 1.0
             Q[state][actions[i]] = returns_sum[state][actions[i]] 
                      / N[state][actions[i]]
   return Q**

该算法将 OpenAI 健身房环境的实例、生成的剧集数量以及折扣率(默认值1)作为参数。该算法返回作为输出的 Q 表(动作值函数的估计)，一个字典(一维数组)。

绘制相应的状态值函数，看看哪些状态值更大，会很有意思。我们可以使用这个简单的代码从 Q-table 中完成，它通过根据我们如何定义问题来加权每个动作的值来获得状态的值:

**num_episodes=1000000Q = mc_prediction(env, num_episodes, generate_episode)State_Value_table={}
for state, actions in Q.items():
           State_Value_table[state]= 
                 (state[0]>18)*(np.dot([0.75, 0.25],actions)) +
                 (state[0]<=18)*(np.dot([0.75, 0.25],actions))**

我们可以用从 Udacity 借来的代码来绘制这个图。有两个图对应于我们是否有可用的 a 或者我们没有可用的 a。但是在这两种情况下，我们看到最高的状态值对应于当玩家和是 20 或 21 的时候，这看起来很明显，因为在这种情况下我们最有可能赢得游戏。

(来源:作者)

探索与开发

到目前为止，我们已经了解了一个代理如何采用一个像等概率随机策略这样的策略，使用它来与环境交互，然后使用该经验来填充相应的 Q 表，该 Q 表成为该策略的动作值函数的估计。所以，现在的问题是，我们如何利用这一点来寻找最优策略？

贪婪的政策

为了得到一个更好的政策，不一定是最优的，我们只需要为每个状态选择最大化 Q 表的行动。让我们称这个新政策为【π】。当我们拿一个 Q 表，并使用最大化每一行 的动作来提出策略时，我们说我们正在构建关于 Q 表的贪婪的策略。****

通常将所选动作称为贪婪动作**。在有限 MDP 的情况下，作用值函数估计用 Q 表表示。然后，为了得到贪婪动作，对于表中的每一行，我们只需要选择对应于最大化该行的列的动作。**

ε-贪婪政策

然而，不是总是构造一个贪婪策略(总是选择贪婪动作),而是构造一个所谓的ε贪婪策略**,它最有可能选择贪婪动作，但也有很小但非零的概率选择其他动作。在这种情况下，使用一定在零和一之间的某个小正数εϵ。这是有动机的，正如我们将在后面更详细地解释的那样，由于代理必须找到一种方法来平衡基于他们当前知识的最佳行为的驱动和获取知识以获得更好的未来行为的需要。**

蒙特卡洛控制

我们已经了解了代理如何采用策略 π ，使用它与环境进行多次交互，然后使用结果通过 Q 表来估计动作值函数。一旦 q 表非常接近动作值函数，代理就可以构造策略 π ，即ϵ——相对于 q 表是贪婪的，这将产生比原始策略 π 更好的策略。然后我们可以改进这个策略，把它变成ϵ的贪婪策略。

因此，如果代理人在这两个步骤之间反复交替，直到我们得到越来越好的策略，希望我们收敛到最优策略，我们最终会得到最优策略π∫:

• 步骤 1 :使用策略 π 构建 Q 表，以及
• 第二步:将策略修改为ϵ——贪婪关于 q 表(标注为*)。***

只要我们运行足够长的时间，这个被提议的算法是如此接近给我们最优策略。

我们把这种算法称为蒙特卡罗控制方法*，用来估计最优策略。通常将步骤 1 称为策略评估，因为它用于确定策略的动作值函数。同样，由于步骤 2 用于改进策略，我们也将其称为策略改进步骤。示意性地，我们可以将其表示为:***

因此，使用这个新术语，我们可以总结出，我们的蒙特卡罗控制方法在策略评估和策略改进步骤之间交替进行，以找到最优策略π∫:**

其中箭头中的“E”表示完整的政策评估,“I”表示完整的政策改进。在下一篇文章中，我们将展示蒙特卡罗控制算法的实现。

然而，一些问题出现了，在构造 ϵ 贪婪策略时，如何设置 ϵ 的值？在下一节中，我们将看到如何实现。

勘探开发困境

回想一下，我们的代理最初不知道环境的动态，由于目标是最大化回报，代理必须通过交互来了解环境。然后，在每一个时间步，当代理选择一个动作时，它根据过去对环境的经验作出决定。而且，本能可能是选择基于过去经验的的行动，将获得最大回报。正如我们在上一节中所讨论的，这种对动作值函数估计贪婪的策略很容易导致收敛到次优策略。****

平衡开发与勘探

这可能会发生，因为在早期，代理人的知识非常有限，并拒绝考虑就未来回报而言比已知行为更好的非贪婪行为。这意味着一个成功的代理人不可能在每个时间点都贪婪地行动；相反，为了发现最优策略，它必须继续改进所有状态-动作对的估计回报。但同时保持一定程度的贪婪行动，以尽快维持回报最大化的目标。这激发了我们之前提出的ϵ贪婪政策的想法。

我们将平衡这两个竞争需求的需求称为探索-开发困境*，其中代理必须找到一种方法来平衡基于其当前知识的最佳行为驱动(开发)和获取知识以获得更好判断的需求(探索)。***

设置ε的值

这种困境的一个潜在解决方案是通过在构建 ϵ 贪婪策略时逐渐修改 ϵ 的值来实现的。对于代理人来说，通过选择探索而不是开发尝试各种策略来最大化回报，开始与环境的互动是有意义的。考虑到这一点，最好的开始策略是等概率随机策略，因为它同样有可能探索每个状态的所有可能行为。设置 ϵ =1 产生一个ϵ-贪婪策略，它等价于等概率随机策略。**

在随后的时间步骤中，促进开发而不是勘探是有意义的，在这种情况下，政策相对于行动价值函数估计逐渐变得更加贪婪。设置 ϵ =0 产生贪婪策略。已经表明，最初倾向于通过开发进行勘探，并逐渐倾向于开发而不是勘探是一种最佳策略。

为了保证 MC 控制收敛到最优策略π∑，我们需要保证满足两个条件:每个状态-动作对被访问无限多次，并且该策略收敛到关于动作-值函数估计 Q 的贪婪策略我们将这些条件称为无限探索极限中的贪婪以确保代理在所有时间步继续探索，并且代理逐渐探索更多和探索更少。

满足这些条件的一种方法是修改 ϵ 的值，当指定 ϵ 贪婪策略时，使其逐渐衰减。然而，在设定 ϵ.的衰减率时必须非常小心确定最佳衰变不是小事，需要一点炼金术，也就是经验。我们将在下一篇文章中看到实现的例子。

下一步是什么？

到目前为止，在我们当前的蒙特卡罗控制算法中，我们收集了大量的片段来构建 Q 表。然后，在 Q 表中的值收敛之后，我们使用该表来提出改进的策略。然而，蒙特卡罗预测方法可以在逐集的基础上逐步实施**。在下一篇文章中，我们将介绍如何按照这个想法构建更好的 MC 控制算法。下期帖子再见！****

深度强化学习讲解系列

由 UPC 巴塞罗那理工大学 和 巴塞罗那超级计算中心

一个轻松的介绍性系列以一种实用的方式逐渐向读者介绍这项令人兴奋的技术，它是人工智能领域最新突破性进展的真正推动者。

** [## 深度强化学习解释-乔迪托雷斯。人工智能

本系列的内容](https://torres.ai/deep-reinforcement-learning-explained-series/)

关于这个系列

我在五月开始写这个系列，在巴塞罗那的**封锁期。**老实说，由于封锁，在业余时间写这些帖子帮助了我 #StayAtHome 。感谢您当年阅读这份刊物；它证明了我所做的努力。

免责声明 —这些帖子是在巴塞罗纳封锁期间写的，目的是分散个人注意力和传播科学知识，以防对某人有所帮助，但不是为了成为 DRL 地区的学术参考文献。如果读者需要更严谨的文档，本系列的最后一篇文章提供了大量的学术资源和书籍供读者参考。作者意识到这一系列的帖子可能包含一些错误，如果目的是一个学术文件，则需要对英文文本进行修订以改进它。但是，尽管作者想提高内容的数量和质量，他的职业承诺并没有留给他这样做的自由时间。然而，作者同意提炼所有那些读者可以尽快报告的错误。**

蒙特卡罗方法，变得简单

原文：https://towardsdatascience.com/monte-carlo-methods-made-simple-91758ba58dde?source=collection_archive---------7-----------------------

来源: Unsplash

用混乱来寻找清晰

想象一个坐标网格上的 10 乘 10 的正方形。在那个网格上画了一些形状，但是你不知道它看起来像什么。但是可以查询一个函数 f ( x ， y )，其中( x ， y )是坐标，输出不是 1(在形状中)就是 0(不在形状中)。你如何找到这个形状的面积？

这种形状的一些可能性。由作者创建

答案很简单。它来自于统计学的一个原理，大数定律，即一个函数被随机抽样的次数越多，它的近似值就越精确。然后，我们的解决方案是简单地在 10×10 网格中随机选择点，计算形状中有多少陆地，然后除以采样点的总数。

虽然这种关于随机抽样的本能想法很简单，但它在许多领域都有应用，从法律到气候预测，也许与本文最相关的是机器学习和统计学。它有一个正式的名字:蒙特卡洛方法。

当遇到由确定性原则组成的问题时，如形状的面积、函数的分布或游戏中玩家下一步应该选择哪一步，蒙特卡罗方法基本上假设它可以通过概率和随机性(随机性)来建模。

蒙特卡罗方法依赖于从分布中重复随机取样来获得数值结果。它是一种方法，而不是一种算法。

对于另一个基于形状的例子，请查看 在不使用任何数学 (使用蒙特卡洛采样和多项式回归)的情况下找到圆面积的公式。随机抽样可以作为一种寻找函数积分(曲线下面积)的廉价方法。众所周知，圆周率是一个重要的常数，圆周率的值也可以通过蒙特卡洛采样来近似计算:

来源:维基媒体。图片免费分享。

一般来说，蒙特卡罗抽样有三类/三种应用:

• 直接取样。在没有先验信息的情况下，天真而直接地从分布中抽样。这就是我们如何逼近一个形状的未知区域。
• 重要性抽样。在分布太昂贵而不能取样的情况下，从更简单的近似函数中取样。这是贝叶斯优化和代理优化的核心组成部分。
• 拒绝采样。在分布未知的情况下，提出新的点，如果它们满足某个标准，就接受它们。

蒙特卡罗抽样通常用于两种情况:

• 优化。自然，找到最佳点需要探索和开发之间的健康平衡。当蒙特卡罗抽样(探索)与另一种控制开采的机制相配合时，它可以成为寻找最优解的有力工具。
• 近似概率&函数。当用另一种方法很难间接评估某些概率或函数(通常是概率函数)时，蒙特卡罗抽样是一种很好的方法。

因为蒙特卡罗方法背后的思想很简单，但可以以相当复杂和创造性的方式使用，所以最好通过几个例子来探索这种思维模式。

首先，考虑马尔可夫链蒙特卡罗(MCMC)方法，它试图在不知道分布是什么的情况下从目标分布中生成随机样本。马尔可夫链——图中每个节点是一个状态，有一定的概率 p 移动到另一个状态——被用来表示这些分布。

考虑一个城市中的天气马尔可夫链(天气恶劣),其中唯一可能的天气状态是风、冰雹/雪、雷暴或雨。每天，第二天的天气可以根据当天天气的概率进行预测。例如，如果今天下雪，有 80%的可能性刮风，20%的可能性明天下雨。

图标: Pixabay 。图示:作者。

为了在马尔可夫链中漫游，我们从一个位置 s 开始，以指定的概率移动到另一个位置 s ’ ，。然后， s 成为新的 s 并重复该过程。虽然这个例子的特点是一个小的马尔可夫链，但是具有数千个节点和数十万个连接的巨大图形可以用来模拟错综复杂的概率关系。

如果我们在大量时间(例如 10，000 个模拟日)后运行马尔可夫链，我们开始达到“概率平衡”。这仅仅意味着，我们可以简单地根据我们旅行过的州有多少在下雨来估计下雨的静态概率(基于大数定律)。

例如，如果我们通过运行超过 10，000 个状态的马尔可夫链得到以下(假设的)结果:

• 2754 个州的州风
• 1034 个州有雷暴
• 4301 个州出现冰雹/降雪
• 为 1911 个州下雨

我们将能够产生以下概率:

• p(风)= 0.2754
• p(雷暴)= 0.1034
• p(冰雹/雪)= 0.4031
• p(下雨)= 0.1911

然后，我们可以简单地从分布中进行相应的采样——从马尔可夫链中随机抽取一个状态，而不需要遍历它。通过重复迭代和马尔可夫链的随机遍历(该过程的“蒙特卡罗”部分)，系统能够被折叠并表示为概率分布。

可以构建马尔可夫链来模拟无法直接采样的复杂关系，然后将其简化以找到其潜在的概率分布。

有许多成熟的 MCMC 算法，如 Metropolis-Hastings 算法或 Gibbs 抽样。所有人都试图做一些类似的事情——模拟系统的潜在概率。

要从理论迈出一步，进入应用，看看蒙特卡罗树搜索(MCTS)，智能游戏强化学习系统的一个关键部分。早期的游戏系统，比如 IBM DeepBlue 在 1997 年击败国际象棋冠军加里·卡斯帕罗夫(Gary Kasparov)的最初胜利，都是基于类似极大极小算法的，这些算法从当前的走法中找出所有可能的游戏，并确定哪一个游戏能带来确定的胜利(或最高的机会)。

随着测试游戏变得越来越复杂——最著名的是围棋游戏，甚至是像 DOTA 这样的高级图形射击游戏——计算上不可行，也很愚蠢，去玩完所有可能的场景。相反，对于一个玩家来说，探索潜在的有利棋步，同时放弃失败概率高的棋步，并利用已知的棋步，是更明智的做法。

假设我们开始构建一个这样的树，其中每个节点代表一些游戏状态。我们可以通过采取某些行动在节点之间转换。从游戏最初的开始状态，我们可以取几个可能的节点。

由作者创建

价值神经网络试图为每一步棋赋予一个价值:给定这一步棋，代理人获胜的概率是多少？注意，刚开始的时候，由于价值网没有经过训练，所以会给出随机的猜测。然而，随着系统玩足够多的游戏，它将更善于智能地分析某些移动的效果。

由作者创建

我们不是简单地选择概率最高的节点，而是执行“加权蒙特卡罗采样”。如果网络认为节点 2 将有 90%的机会获胜，那么它有更大的机会被选中，但也有可能节点 3 和 1 被选中。策略——或如何选择潜在节点——也是可以学习的。

这背后的逻辑是，某一步棋的可预见价值是有限的，是以偶然性为界的(也许对手会出其不意)。此外，只开发已知的移动，而不是冒险探索其他领域并潜在地发现更有利可图的移动，将导致无竞争力和懒惰的模型，这在复杂的强化学习环境中是不可取的。

然后，选择的节点被进一步扩展，并且该过程重复，直到游戏终止。将概率纳入决策中，而不是确定性地选择“最佳”值，有助于强化学习平衡开发/探索权衡。

如果我们从统计学的角度来看这个问题，蒙特卡罗抽样被用来模拟最佳概率分布 p ( v )来选择一个节点，假定它具有一些可预见的值 v 。

蒙特卡洛树搜索是 alpha Go(deep mind 的强化学习系统)击败(并将继续击败)顶级围棋选手的框架。MCTS 在其他地方有各种各样的应用。

模拟退火是蒙特卡罗抽样的另一个应用，在寻找全局最优解的任务中用作非梯度函数优化的有效方法。像应用蒙特卡罗方法的其他问题一样，搜索空间是离散的(例如，我们不是试图优化连续值，像神经网络的参数)。

在一些问题中，在固定的时间内找到一个近似的全局最优值比它的精确值更有价值，模拟退火实际上可能比梯度下降等算法更好。

模拟退火的想法来自冶金学，或金属的操作。在冶金学中，退火是对材料进行有控制的加热或冷却，以增加其尺寸并消除缺陷。类似地，模拟退火控制系统中的能量，这决定了它在探索新的可能性时愿意冒多大的风险。

温度从某个初始量开始，它代表某种计时器。在每个状态 s ，模拟退火根据两个状态的值移动到某个相邻状态s’(移动到s’是收益还是损失？)，以及当前温度 T ，以概率 P ( s ，s’， T )。

随着温度随着每个时间步长降低(剩余时间减少)，该模型变得更少探索性，而更多利用性。模拟退火允许在时间上从探索到开发的渐进过渡，这对于寻找全局最优非常有益。

随着温度的降低，寻找复杂函数的全局最大值的模拟退火。来源:维基媒体。图片免费分享。

采用蒙特卡罗抽样和贝叶斯方法对概率函数 P ( s ，s’， T )进行建模。事实上，您可能还记得，Metropolis-Hastings 算法通常是一种马尔可夫链蒙特卡罗方法(或以它为模型的方法)，用于寻找转换阈值(应该转换的概率)。

这很自然，因为模拟退火将解决方案视为状态，并试图找到最佳转移概率——这是马尔可夫链建模的完美场景。

通常，蒙特卡罗方法——或者类似蒙特卡罗的思维——会出现在我们最意想不到的地方。虽然这是一个简单的机制，但它在无数的应用中有着深刻而复杂的根源。

总结/要点

• 蒙特卡罗方法是基于这样一种思想，即在系统中注入随机性通常可以有效地解决这个问题。
• 一般来说，蒙特卡罗抽样有三类:直接抽样、重要性抽样和拒绝抽样。
• 蒙特卡罗的两个常见应用包括优化和复杂概率和函数的评估。
• 在离散(非连续)和确定性问题中，蒙特卡罗方法利用随机性+概率、大数定律和有效的框架来有效地解决它们。

感谢您的阅读，请在回复中告诉我您的想法！

如果你对最新的文章感兴趣，可以考虑订阅。如果你想支持我的写作，通过我的推荐链接加入 Medium 是一个很好的方式。干杯！

[## 贝叶斯优化的美妙之处，用简单的术语解释

巧妙算法背后的直觉

towardsdatascience.com](/the-beauty-of-bayesian-optimization-explained-in-simple-terms-81f3ee13b10f)

采用美元成本平均法的蒙特卡罗模拟

原文：https://towardsdatascience.com/monte-carlo-simulation-with-dollar-cost-averaging-653ae47ec7d5?source=collection_archive---------26-----------------------

艾萨克·史密斯在 Unsplash 上拍摄的照片

利用 Python 中的蒙特卡罗模拟分析标准普尔 500 ROI

声明:我不是金融专业人士。此处提供的内容仅用于信息、教育和娱乐目的，不应被解释为专业的财务建议。

从我完成研究生学业开始我在分析领域的新事业已经有几个月了。在我拿到最初的几份工资后，我开始思考如何更明智地使用每月支付费用后剩余的现金，因为把它留在利息为 0.01%的支票账户上似乎有点乏味。经过一番研究后，我决定为什么不把它投入股市，至少可以跟上通货膨胀，这样我就可以将美元成本平均到标准普尔 500 交易所交易基金(ETF)中了。

美元成本平均法(DCA)是一种常见的投资策略，其中固定金额的资本定期投资于某项资产，以减少市场波动的影响。这种策略与 ETF 完美契合了我的需求，因为它是自动多样化的，不需要多少市场知识。

在我把辛苦赚来的钱投入股市之前，我想我应该好好利用我的分析学位所获得的新知识。我决定在 ETF 上用 Python 实现一个简单的蒙特卡罗(MC)模拟，看看 DCA 的表现如何。通过 MC 模拟，我可以在对历史数据进行重复随机采样后，获得预期投资回报(ROI)的分布。

数据

为了做到这一点，我研究了两个最受欢迎的基金，SPDR 标准普尔 500 ETF 信托基金(SPY)和先锋标准普尔 500 ETF (VOO)。我从雅虎上下载了这些数据。财务，从开始日期到当前日期的所有日期。模拟中使用的两个关键变量是交易日和每个交易日结束时的调整收盘价(因为收盘价不包括股息)。没有考虑成交量和日内价格。

间谍数据帧的前五行

通过快速绘制数据，我注意到曲线的形状几乎相同(这是意料之中的，因为它们跟踪的是相同的股票)。间谍存在的时间也是 VOO 的两倍多，在过去的三月里，大概是由于新冠肺炎·疫情的原因，它的大量减少非常明显。

幸运的是，我不必为这个简单的时间序列做任何数据清理，所以我们可以直接进行模拟。

仿真法

构成 MC 模拟的主要组件是过程、常量变量、随机变量和样本大小。在这种情况下:

• 流程: DCA 在 12 个月内每个日历月进行一次间谍和 VOO
• **恒定变量:**每月 1000 美元投资
• **随机变量:**选择一个连续 12 个月的时间段和每个日历月一次的购买日期
• 样本量: 20，000 次(间谍)和 10，000 次(VOO)迭代

模拟开始时，银行里有 1000 美元作为第一个月的购买力，在股票存在的范围内随机选择 12 个月的时间，间谍从 1993 年到 2020 年，VOO 从 2010 年到 2020 年。在随机过程生成器中，允许 12 个月的周期从一年开始，到下一年结束。所有随机过程发生器都运行在均匀分布上。

对于所选日期范围内的每个月，从历史数据中随机选择该月的一个日期作为购买日期。日期随股价而定，用可获得的购买力购买最大数量的股票。除了每月 1000 美元的新收入外，剩余的购买力会滚动到下个月。

在 12 个月的期间内，每个月都重复该过程的这一阶段，以创建一个样本。整个过程对于 SPY 重复 20，000 次，对于 VOO 重复 10，000 次。

每个数据点都在一个循环函数中计算并存储在 Python 中。请查看我在 GitHub 上的代码了解更多细节。

注意:SPY 中的额外迭代是为了说明其更长的存在时间(需要更多的样本来收敛到总体均值)。

结果

在对间谍和 VOO 进行模拟后，可以计算出每个样本从第 1 个月到第 12 个月的年投资回报率，从而得出样本分布。当我们可视化结果时，我们看到一致的分布，稍微偏向负回报，中心在 6–7%的年回报。SPY 中的分布更加广泛，偏离平均值更多，这可能来自 VOO 没有的 1993 年至 2010 年间的额外数据。

间谍和 12 个月投资回报样本分布

12 个月回报统计

无论某人投资哪只 ETF，都有 75%的机会获得至少 1%的回报，这已经是 0.01%利息支票账户的 100 倍。更妙的是，SPY 的最高回报是 27%(虽然概率很低)！

出于好奇，我还对输入进行了一些调整，并模拟了 60 个月的 ROI，以观察它的表现。间谍变得不那么一致，变得更加多变，而 VOO 没有太大变化。

间谍和 VOO 的 60 个月投资回报样本分布

60 个月退货统计

在 5 年的时间里，投资者可以预期获得 30%左右的平均回报。最好的情况是，一项投资可能会因为间谍而翻倍，或者在最坏的情况下损失 40%。

结论

根据这一分析，每月随机选择一天用 DCA 将钱投入标准普尔 500 ETF，几乎肯定比把钱留在支票账户要好。与任何投资一样，风险越高，回报也越高。正如我们在这里看到的，平均 12 个月的投资回报率为 6–7%也有可能出现负投资回报率。SPY 提供了一个略有不同(也许更好？)由于起步较早，对市场的估计比 VOO 大。周期性将是一个有趣的分析点，可以看出每月分期付款与每季度、每两个月分期付款相比如何。最后，这只是一个有趣的项目，给了我一个练习分析技能的借口，并希望获得一些对我个人财务有用的见解。

如果你有兴趣查看我的代码，请查看我的 Github ！