用 Qiskit 模拟量子计算机
Qiskit 量子模拟器-如何使用它们以及它们的用途
量子机器学习要不要入门？看看 动手用 Python 学习量子机器。
在量子计算中，我们使用量子物理的基本性质来执行计算:叠加和纠缠。
叠加是指量子系统存在于量子态|0⟩和|1⟩.的复杂线性组合中的量子现象纠缠是量子粒子之间极强的关联。即使相隔很远，纠缠的粒子仍然保持完美的关联。
仅仅是我们能够利用这些现象进行计算这一事实就令人震惊。但是，这种类型的计算使我们能够解决看似无法解决的问题的事实是不可思议的。
令人惊讶的是，我们甚至不需要量子计算机。我们可以用一个日常笔记本来模拟它们。
假设我们有下面的量子电路。
作者图片
我们在|0⟩基态的量子位上应用哈达玛门，把它放入|+⟩.态在这种状态下，我们有 50%的几率将量子位测量为 0，有 50%的几率将其测量为 1。
下面的清单在 Qiskit——IBM 的量子计算 SDK——中描述了这个电路。
Qiskit 提供了Aer包。它为模拟量子电路提供了不同的后端。先说第一个，那个qasm_simulator。
一旦我们用qasm_simulator后端(或任何其他后端)执行了我们的量子电路(qc，我们就可以使用job.result()方法获得结果。我们可以将这个结果转换成数字(get_counts())，然后输入到最终分布的直方图中。
尽管我们的量子电路产生了一个量子位，我们以 50%的几率将其测量为 0，但结果显示了略有不同的分布。
作者图片
这是因为qasm_simulator根据经验检索测量计数。shots参数告诉 Qiskit 运行电路多少次并获得测量结果。因此，结果并不完全准确。但是拍摄次数越多，结果就越准确。
qasm_simulator只考虑我们测量量子位时使用的经典位。按照量子力学，我们看不到一个量子位所处的真实状态。我们只能反复测量才能得出它的状态。
当然，知道量子位的确切状态会有所帮助——尤其是在量子电路的开发过程中。
幸运的是，Qiskit 提供了另一个后端，即statevector_simulator。这个模拟器计算一个量子位的确切状态。当我们使用这个模拟器时，去除所有的测量值是很重要的，因为测量会破坏量子叠加，不可避免地导致系统可能处于一个确定的状态。但是我们对坍缩前的量子叠加感兴趣。
下面的代码描述了更新的量子电路。注意，statevector_simulator不接受 shots 参数。
类似于qasm_simulator，我们可以从执行结果中得到计数。
作者图片
结果显示了精确的测量概率。
但是它是怎么做到的呢？
为了理解经典计算机如何计算量子电路，我们需要考虑量子态(这里是|0⟩)和量子算符(这里是 T9 h T10)到底是什么。
在量子力学中，我们使用看起来像|0⟩.的狄拉克符号这不是什么奇特的东西，只是一个简单的向量。
量子算符通常用大写字母表示，比如 H 。但它们也没什么特别的。量子算符是一个矩阵。
最后，量子电路一点也不花哨。就是矩阵乘法。是的。整个电路将矩阵 h 与向量|0⟩相乘，并产生另一个向量，如下式所示。
结果向量表示测量幅度，其平方是测量概率。1/sqrt(2)的平方是 1/2。
从结果对象中，我们也可以使用get_statevector()方法获得结果向量。输出是状态向量的数组。在我们的例子中，它是:
array([0.70710678+0.j, 0.70710678+0.j])
注:1/2 的平方根约为 0.707 左右。
statevector_simulator后端计算给定量子系统的状态。这是检查叠加态量子位的理想方法。
我们的量子电路由一个单一的算子组成——哈达玛门。比方说，我们用一系列的门来代替。例如，我们可以应用三个门，H-Z-H(在这篇文章中了解更多关于 HZH 序列的信息)。
作者图片
当我们运行电路时，我们看到结果为 1(如果量子位处于|0⟩).态如果你想从数学上验证这个结果，我们需要将三个矩阵相乘。
即使对于小矩阵，它们的乘法也很麻烦。你需要做相当多的乘法，你必须确保不要混淆行和列。
所以，我一直用电脑做矩阵乘法。有无数的计算器可以做到这一点。然而，最方便的是在 Qiskit 内部。就是这个unitary_simulator。
这个模拟器执行电路一次，返回电路本身的最终变换矩阵。不管我们的量子电路看起来像什么，最终，它是一个矩阵。在输入状态下运行电路就是将转换矩阵与状态向量相乘。同样，当我们使用这个模拟器时，我们的电路不应该包含任何测量。
以下代码输出 H-Z-H 电路的矩阵。正如我们所看到的，它类似于非门的矩阵。因此，毫不奇怪，当应用于|0⟩态的量子位时，结果是|1⟩态的量子位，我们总是测量为 1。
array([
[ 2.22044605e-16+6.1232340e-17j, 1.00000000e+00-1.8369702e-16j],
[ 1.00000000e+00-6.1232340e-17j, -2.22044605e-16+6.1232340e-17j]
])
如果量子计算就是矩阵乘法，有什么了不起？
矩阵乘法没什么花哨的。但问题是，这需要大量的计算。矩阵越大，所需乘法运算的次数(指数)就越高。因此，矩阵乘法花费的时间越多。我们可以用传统的计算机乘任意大小的矩阵。唯一的问题是，对于大型(假设是巨大的)矩阵来说，它(几乎)需要花费很长时间。大多数经典计算机无法(有效地)处理大于 30x30 的矩阵。
这就是量子计算机发挥作用的地方。他们在一个步骤中乘以许多矩阵。矩阵的大小无关紧要。因此，它们能够执行经典计算机无法执行的计算(在合理的时间内)。
结论
量子计算机是令人着迷的设备，尤其是当我们想到它们的硬件以及它们如何利用量子力学来执行计算时。这远远超出了大多数人(包括我)所能理解的范围。
但是量子计算并没有那么复杂。就是矩阵乘法，本质上就是。在小范围内，经典计算机在矩阵乘法方面做得相当好(而且准确)。因此，经典计算机可以模拟量子计算机做的事情。
但是随着矩阵规模的增加，经典计算机已经达到了极限。然后，没有模拟器会帮助你了。
量子机器学习要不要入门？看看 动手用 Python 学习量子机器 。
在这里免费获得前三章。
用 Python 中的 SimPy 模拟酿造操作
python 中的离散事件模拟教程
罗伯特·北原惠子·桑塔纳在 Unsplash 上拍摄的照片
工业工程师使用离散事件模拟来评估生产车间的替代方案，但模拟还有更多用途。有了 Python 的 SimPy 包——也有其他语言版本——我们可以模拟个人项目，而无需昂贵的软件许可。SimPy 的 ReadtheDocs 文档很好；我建议从头到尾多看一遍。但是这篇文章是为那些喜欢额外的学习例子的人准备的。
SimPy 简介及示例
SimPy 使用带有 yield 而不是 return 的生成器函数，以便存储模拟中的最新步骤。关于生成器函数的介绍，我推荐凯文·康拉德的 youtube 教程。
这个初始模型利用 SimPy 框架来启动酿造过程，然后清洗酿造罐。后续的模型会更复杂。
在导入 simpy 之前，如果您还没有安装 SimPy，您可能需要安装它。要安装 SimPy，您可以使用 pip 或 SimPy 文档中编写的替代方法。
pip install simpy
输出:
Starting to brew a lager at time 0
Finished brewing the lager at time 30
Starting to brew a lager at time 35
Finished brewing the lager at time 65
Starting to brew a lager at time 70
Finished brewing the lager at time 100
在这种情况下，模型在时间单位达到 101 时终止。一旦处理所有酿造过程的循环终止，下一个模型将终止。
在 SimPy 中使用资源和建模流程
现在我有 3 个同样大小的啤酒缸。我仍然在酿造淡啤酒，但是我想在第 0、5、11 和 20 天开始一个新罐，如果这些天有可用的罐的话。我在第 12-18 天休假，所以我不想在那几天开始酝酿过程。
在这个模型中，贮藏啤酒的酿造过程时间仍然假定为 30 天，尽管我意识到酿造过程时间实际上是可变的。我们将在下一个模型中添加处理时间的概率分布。
目前，这个模型的焦点是资源的类型，称为资源(而不是容器资源或商店资源)。
- 这个 brewery_tanks 资源将作为这个模型中的流程。
- 资源被请求，然后在它们不再被使用后释放。
- 在使用中，资源的计数将显示有多少资源的容量单位在使用中。
输出:
Decide to start brewing a lager at time 0
Number of tanks available: 3
Started brewing at time 0
Number of tanks now in use: 1
Decide to start brewing a lager at time 5
Number of tanks available: 2
Started brewing at time 5
Number of tanks now in use: 2
Decide to start brewing a lager at time 11
Number of tanks available: 1
Started brewing at time 11
Number of tanks now in use: 3
Decide to start brewing a lager at time 20
Number of tanks available: 0
Finished brewing a lager at time 30
Finished brewing a lager at time 35
There are 1 tanks available at time 35
Started brewing at time 35
Number of tanks now in use: 3
There are 1 tanks available at time 40
Finished brewing a lager at time 41
There are 2 tanks available at time 46
Finished brewing a lager at time 65
There are 3 tanks available at time 70
我选择了一个很好的时间去度假，因为当时没有任何可用的坦克。事实上，直到时间 35 都没有可用的坦克。
如您所见，仿真模型的输出值得思考。文档中的一个章节叫做 monitoring 讨论了保存和监控模型输出的更高级的方法。但是对于这个介绍和许多其他场景，print 语句就足够了。
随机模拟
虽然可以用数学方法计算确定性模拟结果，而不是利用模拟模型，但模拟对于讲述故事和自动化数据处理非常有用。随机模型将为更接近真实情况的建模场景提供更多的见解。
我阅读了多种描述酿造啤酒的资源，并分享了典型的酿造时间表是 4-8 周。酵母会做酵母做的事情，但不总是按照你的计划。让我们假设一个均值=6 周，标准差=1 周的正态分布。一旦啤酒冷却(我们假设这是过程的一部分)，然后我们将罐或瓶啤酒并储存在库存中，比率为每个啤酒厂罐 100 箱。
为了跟踪库存，这个模型中引入了另一种称为容器的资源。虽然容器类似于资源，因为它们都具有容量参数，但是容器也具有初始的级别参数，并且容器在给定时间的级别将揭示容器存储的量。
容器的级别可以随着一个 put 而增加，或者随着一个 get 而减少，这实现了与资源不同的建模能力。
输出:
Number of tanks available at time 0: 3
Started brewing at time 0 and 1 tanks now in use
Number of tanks available at time 5: 2
Started brewing at time 5 and 2 tanks now in use
Number of tanks available at time 11: 1
Started brewing at time 11 and 3 tanks now in use
Number of tanks available at time 20: 0
Finished brewing a lager at time 38.30356066389249
At time 43.30356066389249 there are now:
1 tanks available
100 cases of cans available for sale
Started brewing at time 43.30356066389249 and 3 tanks now in use
Finished brewing a lager at time 44.85515309121678
At time 49.85515309121678 there are now:
1 tanks available
200 cases of cans available for sale
Finished brewing a lager at time 53.146254345432425
At time 58.146254345432425 there are now:
2 tanks available
300 cases of cans available for sale
Finished brewing a lager at time 85.11565243779629
At time 90.11565243779629 there are now:
3 tanks available
400 cases of cans available for sale
处理订单
现在您已经了解了如何使用两种类型的资源，即资源和容器，下一个模型将分享一种处理订单和销售库存的方法。虽然这是一个处理订单的简单方法，但是更复杂的方法可能需要使用任何简单事件的 AllOf。
输出:
Before order processing, cases available: 100
After order of 1, cases available: 99
After order of 10, cases available: 89
After order of 5, cases available: 84
After order processing, cases available: 84
我希望本教程对你有所帮助。如果你对第二部分啤酒厂模拟有什么建议，或者想在下一个 SimPy 模型上合作，我很乐意阅读你的想法！
用气流模拟连续学习模型
配置工作流管理工具，在 13 分钟内定期重新创建模型
C 持续学习是系统自适应学习外部世界的过程。在机器学习的背景下，模型通过自主和增量开发不断地学习更多关于它被训练的主题领域的知识。这使模型能够学习更复杂的逻辑、知识和功能，从而有助于做出更好的预测。
持续学习对于更新模型尤其重要。假设你被一家生产哑铃的制造公司聘为数据科学家，比如鹦鹉螺公司。你的目标是预测下个季度的收入和销售额。从 2017 年 3 月到 2020 年 3 月，您的模型在收入预测上的准确率为 95%。但是，从 2020 年 3 月到 2021 年 1 月，你的模型的准确率是 40%。发生了什么事？
2020 年 3 月—2021 年 1 月发生的重大事件是冠状病毒。检疫条例迫使人们呆在室内，不要频繁外出。不能再去健身房的健身爱好者不得不建立他们的内部健身房，以跟上他们的锻炼养生。鹦鹉螺公司是哑铃和家用健身设备的领先公司，因此在疫情期间销售额飙升。在这种情况下，你在 2017 年创建的模型将严重低估 2020 年的预测销售额。这是一个很好的例子，说明为什么一个模型需要根据最近的和新出现的当前数据进行重新训练。世界的发展速度令人难以置信，3 年前的数据见解可能不适用于不断发展的数据。
您可能已经意识到需要创建 3 个独立的模型:一个在疫情之前(2020 年 3 月之前)接受数据培训，一个在疫情期间(2020 年 3 月-2021 年 5 月)接受数据培训，一个在疫情之后(2021 年 5 月至今)接受数据培训。但是，您可能没有足够的时间按照特定的时间段清晰地分割数据集。因为您可以假设疫情的开始和结束日期(例如，假设疫情的结束时间是疫苗推出的时间)，所以您可以准确地知道在哪里分割数据集。
随着时间的推移，各种数据集可能会发生微妙的变化，您需要应用一种持续学习的方法，以便您的模型能够随着新信息的出现而学习。这些数据既反映了收集数据的动态环境，也反映了提供数据的形式和格式。作为后一种发展的例子，考虑一个图像分类器，它确定 1000 页 pdf 中的一页是文本还是表单。随着时间的推移，表单的格式会发生变化，原始图像分类器会将其识别为文本。因为您不知道表单更改的确切时间，所以您需要根据所有当前数据训练一个模型。此外，图像分类器必须在较新的表单数据上重新训练，同时跟踪较旧的表单数据。
一个模型需要一个持续的学习方法来基于更新的数据进行再训练和更新。也就是说，用数学方法实现深度学习模型和递归神经网络的这种方法可能非常棘手。许多数据科学家不具备在公司期限内构建这种定制模型的知识。但是，有一种方法可以使用工作流管理工具(如 Airflow )来模拟持续学习。在预定的基础上，Airflow 可以获取新的训练数据，并使用以前的和新的数据重建模型。然后，数据科学家可以构建简单的模型(无论是他们自己构建的简单 PyTorch 神经网络，还是他们利用现有的 Python 机器学习库，如 scikit-learn)，并专注于持续更新训练数据。
本文将通过一个例子来说明如何对时间序列模型应用连续学习。我们的模型将使用 ARIMA 分析两只股票的历史股票市场数据:Fastly 和 Nautilus。
**免责声明:**这不是投资建议。话虽如此，我确实持有这两只股票。所以我倾向于提升他们，这样我就能从中受益。在花你的血汗钱之前做你的研究。
另外，在你决定用 ARIMA 预测股市之前，请多读一些时间序列方面的书。当数据集中存在可识别的模式时，时间序列模型是有效的。我没有花时间去分析 ARIMA 是否适合预测 Fastly 和 Nautilus 的股票历史数据。本文从数据工程的角度讨论构建持续学习架构。ARIMA 和时间序列被用来表明，这种架构可以用于各种类型的模型(时间序列，分类，神经网络)。分析这种时间序列模型的预测指标超出了本文的范围。
为什么是气流？
Airflow 是一个用 Python 编写的开源工作流管理工具。它已经成为行业标准，并且很容易建立。在处理复杂和动态的工作流时，气流也有缺点，最好使用不同的工作流工具，如 Prefect。对于本教程，我们将创建一个简单的工作流。因此，气流是完美的。
我们将在 AWS EC2 上设置气流。如果你愿意，你也可以利用 AWS MWAA(亚马逊管理的 Apache Airflow 工作流)。但是，请记住，对于一些开发人员的预算来说，AWS MWAA 可能很昂贵。例如，我每天只需支付 27 美元就可以让它正常运行，其他 AWS 弹性云服务(例如 NAT 网关)每月支付 80 美元。
术语
- AWS MWAA-亚马逊管理的 Apache Airflow 工作流。处理气流的所有设置和配置，但维护费用昂贵。
- AWS EC2 —运行应用程序的虚拟服务器。
- AWS S3——存放腌制模型的桶。
- ARIMA——自回归综合移动平均。它是一类模型，捕捉时间序列数据中一组不同的标准时间结构。
推荐阅读
本教程假设您对气流和 Dag 有基本的了解。如果没有，请阅读气流:如何以及何时使用。
要理解气流的架构，请阅读Apache 气流架构概述。
本教程使用 ARIMA 进行时间序列分析。这个主题超出了本文的范围。如果你想更多的了解时间序列和 ARIMA，我推荐Arima 简介:非季节模型。
体系结构
我使用 Cloudcraft 创建的气流架构图
辅导的
在本教程中，我们将在 Linux Ubuntu Server 20.04 LTS 上安装 Airflow。我们可以使用类型 t2.medium 创建一个 AWS EC2 映像。要阅读更多关于 AWS EC2 的内容，我推荐亚马逊的教程这里。
Airflow 需要运行一个 web 服务器和一个调度程序，因此它需要比自由层提供的内存更多的内存。t2.medium 有足够的资源来安装和运行 Airflow。
我们还使用 Ubuntu，因为它是为数不多的适合机器学习的 Linux 操作系统之一。其他 Linux 操作系统(包括 AWS Linux)将不会有一些机器学习库所需的确切包和依赖项，包括脸书先知。他们有这些包的替代品，但如果脸书先知不能正确编译它们也没关系。脸书预言家需要一个被称为 C++14 的 C++标准来编译它所有的依赖项。Ubuntu(和 Mac OS)支持 C++14，但 CentOS 或 AWS Linux 不支持。Ubuntu 和 MacOS 对于下载未来的机器学习包并正确编译非常有用。
创建 EC2 实例
登录 AWS 控制台并导航到 EC2 >实例>启动实例。配置您的实例，使其与以下内容匹配。
**注意:**我删除了安全组中 Source 下的 IP 地址。但是为了保护您的实例以确保没有其他人可以访问它，请确保您选择了我的 IP 而不是 Anywhere)
**编辑:**除了这些安全组，还应该添加 HTTP，TCP，端口 8080。默认情况下，该端口由 airflow webserver 使用。如果您想使用不同的 airflow 服务器端口，您应该更改端口。
启动实例时，它会询问您从哪里下载. pem 文件。将此文件命名为 airflow_stock.pem，并将其存储在您选择的文件夹中。
等待 EC2 实例运行，然后在本地终端上调用这个命令。
ssh -i "airflow_stock.pem" ubuntu@ec2-3-131-153-136.us-east-2.compute.amazonaws.com
SSH 是一个安全的 Shell 协议，允许您在本地命令行上登录到云实例。它需要一个. pem 文件来验证您是否可以登录到服务器。
**注意:**如果你得到一个未受保护的私钥错误，这意味着你的。pem 是大家看得见的。您只需要调整. pem 的权限。在执行之前运行此命令
chmod 400 airflow_stock.pem
有关此错误的更多信息，请参见本文这里的。
安装 Pip 和气流
假设你有一个干净的 Ubuntu 实例，你需要安装 pip 和 apache-airflow。要安装 pip，请运行以下命令
sudo apt update
sudo apt install python3-pip
我们将在教程中使用 Python 3。要安装 apache-airflow，请运行以下命令
pip3 install apache-airflow
下一步是创建文件作为气流项目的一部分。
创建文件夹结构
您可以创建以下文件，或者从 git repo 中获取。
该文件夹有 4 个组件:一个用于添加气流环境变量的 bash 脚本(add_airflow_env_variables.sh)，一个用于安装 pip 包的需求文本(requirements.txt)，一个启动气流脚本(start_airflow.sh)，一个停止气流脚本(stop_airflow.sh)，以及一个 dags 文件夹，该文件夹提取股票数据并从所有数据中重建 ARIMA 时间序列模型。
Visual Studio 代码中气流文件夹结构的屏幕截图
add _ airflow _ env _ variables . sh包含了我们在运行 air flow 之前需要设置的所有环境变量。现在，它只列出了 AIRFLOW_HOME，我们存储 AIRFLOW 的日志/输出/数据库的路径。我在 ec2 实例中克隆了这个 repo，所以我使用了 github repo 的文件路径。
这个 requirements.txt 文件包含了我们需要的所有 pip 包。
statsmodels
pandas
pandas_datareader
boto3
awscli
要安装 requirements.txt 中的所有包，只需运行
pip3 install -r requirements.txt
start_airflow.sh 包含在 EC2 服务器上启动 airflow 的所有命令。
stop_airflow.sh 包含停止当前在 EC2 服务器上运行的 air flow web 服务器和调度程序的所有命令。
stock_data_dag.py 将包含触发工作流的逻辑。工作流分为两个功能:提取和加载。
Extract(get _ stock _ data)将从 TIINGO 获取股票从 2015 年 1 月 1 日到触发日的所有历史价格。对于本教程来说，这个平台已经足够了。然而，在生产环境中，平台将要求只获取当天的数据，并追加到存储在 S3 中的以前的历史数据。
Load(store _ arima _ model _ in _ s3)从历史股票价格中创建 ARIMA 模型，并将腌制文件存储在 S3 桶中。我们可以获取那些根据新的训练数据重新训练过的新模型。
**注意:**您需要创建一个 TIINGO 帐户来获取 get_stock_data 中使用的 API 键(第 25 行)。
将 python 函数 get_stock_data 和store _ ARIMA _ model _ in _ S3视为任务。为了在 DAG 上调用这些任务，我们需要使用操作符。操作符是 DAG 的主要构件，用于执行某些功能/脚本/API/查询等。在本例中，我们创建了 4 个不同的 Python 运算符:
- 获取 NLS 股票数据
- 商店 _ NLS _ ARIMA _ 模型 _in_s3
- 获取 _ FSLY _ 股票 _ 数据
- 存储 _FSLY_arima_model_in_s3
工作流时间表的格式为 crontab:0 17 * * 1–5。这意味着该工作流将在每个工作日的下午 5 点被触发。
stock_data_dag 为每只股票列出了两个不同的任务:get_stock_data 和 store_stock_arima_model 到 s3。 get_stock_data 获取自 2015 年 1 月 1 日以来的 TIINGO 历史股票数据，并将其存储在 JSON 中。 store_stock_arima_model 获取股票 JSON，获取所有日期的调整后收盘价，根据这些价格构建 arima 模型，并将 ARIMA 模型存储在名为stock-model的 S3 存储桶中。然后我们可以从 S3 下载模型，用它们来预测第二天的价格。
**注意:**这种架构不仅仅适用于 ARIMA 车型。各种其他 skikit-learn 模型和算法，如 XGBoost、RandomForest、NaiveBayes、LinearRegression、LogisticRegression 和 SVM，都可以从这种持续学习方法中受益。您可以简单地重用 stock_data_dag.py 文件的store _ ARIMA _ model _ in _ S3函数，或者直接编辑 python DAG 并将第 44 行更改为您想要的模型。您甚至可以添加自己的神经网络逻辑来代替 skikit-learn 模型，也就是说，您必须重新考虑您的回归/分类模型是基于哪些数据进行训练的。您可能需要更改 get_stock_data 来为您的模型获取不同的数据集。无论您试图解决什么样的业务问题，这种气流持续学习架构都可以互换，以适用于各种用例。
我们现在有了本教程计算所需的所有文件。下一步是执行某些脚本，并为运行配置气流。
创建 AWS S3 时段
我们想创建一个 S3 桶来存储我们的模型。导航到 S3 并创建一个名为stock-model的新存储桶。然后，我们将配置 EC2 aws 凭证，以便我们可以使用 boto3 来存储进入 S3 的气流的模型。
配置 AWSCLI
我们想把我们的 ARIMA 模型存放在一个 S3 桶里。为此，我们需要在 EC2 实例上配置本地 AWS 设置，以便它能够识别所需的 AWS 帐户。
首先，导航到 IAM 管理控制台的“安全身份证明”页面: IAM 管理控制台(Amazon.com)。
IAM 的安全身份证明页
如果您还没有这样做，向下滚动访问键，并点击创建新的访问键。这将生成一个 excel 文件，您可以下载该文件来获得您的访问密钥和秘密访问密钥。不要与任何人共享此文件。如果您这样做，您的 AWS 帐户将面临黑客攻击和可能的超额收费的风险。
在 EC2 上，键入
aws configure
AWS 将提示您输入访问密钥和秘密访问密钥。将 excel 文件中的值放在那里。您可以在默认区域名称和默认输出上按 enter 键。
AWS 凭证现在设置在 EC2 上，因此您可以创建 boto3 会话。
配置气流
首先，添加环境变量。因为我们有一个脚本，您可以简单地调用以下命令。
source add_airflow_env_variables.sh
bash add_airflow_env_variables.sh
如果您想检查是否添加了环境变量 AIRFLOW_HOME，请在命令行中运行env。
下一步是创建用于存储 dag 运行的气流数据库。如果这是您第一次创建 airflow 数据库，请在 AIRFLOW_HOME 路径中运行下面的命令。
airflow db init
您还需要添加一个用户，因为 airflow 会提示用户登录。以下命令摘自本地运行气流文档中的示例。
airflow users create \
--username admin \
--firstname Peter \
--lastname Parker \
--role Admin \
--email spiderman@superhero.org
它会提示您输入自己选择的密码。
接下来，使用以下命令安装 requirements.txt 中的所有包
pip3 install -r requirements.txt
最后，运行以下命令来启动 daemon 中的 airflow webserver 和调度程序。
bash start_airflow.sh
这将运行气流调度和网络服务器没有任何问题。要检查 UI 是否工作，请导航到此 url。
http://ec2–11–111–111–11 . us-east-1 . compute . Amazon AWS . com:8080
ec2–11–111–111–11.us-east-1.compute.amazonaws.com 是一个虚构的 Ec2 公共 IPv4 地址。您应该能够在 EC2 -> Instances 控制台上的 Details 选项卡中为您的 EC2 实例获取它。确保您使用 http (不是 https)和端口 8080(如果您设置气流使用该端口，否则使用您配置的任何端口)。
您应该会看到一个登录页面弹出。输入用于在 Airflow 中创建用户帐户的用户凭据(用户名:admin，密码:您键入的任何内容)。现在，您应该可以看到列出的所有气流 DAG，包括我们之前创建的 stock_data DAG。
气流中所有 Dag 的列表，启用了 stock_data。
单击股票数据 dag。然后，您可以在图形视图中看到这个 dag 的 UI。
股票数据的 DAG 视图
在类型的右上角(在“时间表”下)，您会看到一个“播放”按钮。点击该按钮手动触发 DAG。如果一切正常，您应该会看到所有操作员都成功通过绿色认证。如果有错误(失败，等待重试)，单击操作符检查日志。
成功的 DAG 运行
现在，让我们导航到 S3 桶stock-model，看看我们的模型是否已经被添加。
添加了模型的库存模型桶
如果让 EC2 整夜运行，并且气流调度程序和 web 服务器不中断，您将看到这两个模型最后一次修改是在 7 月 23 日下午 5:00 UTC-05:00 左右。
注意:确保完成后终止 EC2
结论
本教程解释了如何使用工作流工具构建持续学习架构，以在当前和新的训练数据上重建 ARIMA 模型。我们学习了如何利用气流来创建这种持续学习架构，以预测两种不同股票的收盘价:Nautilus 和 Fastly。我们学习了如何在亚马逊 EC2 上部署气流，并在 S3 上存储新建立的模型。
持续学习不需要从头开始创建循环神经网络。这是通过使用流行的数据工程工具来创建持续学习模型的一种更简单的替代方法。仍在学习错综复杂的神经网络的入门级数据科学家可以依靠这一解决方案来创建一个健壮、可持续的架构，用于模型“学习”更新的数据。此外，这种架构可以在各种不同的算法上复制:时间序列、线性回归和分类算法。虽然读取和训练数据模型的逻辑需要改变，但总体架构应该适应各种业务问题。
Github Repo:【github.com】Data-Science-Projects/Medium/air flow _ CL/air flow at master hd2zm/Data-Science-Projects(
参考文献:
1.Matthew Kish 鹦鹉螺公司报告称，由于疫情推动了家庭锻炼趋势，销售额增长了 94%
感谢阅读！如果你想阅读更多我的作品，请查看我的目录。
如果你不是一个中等收入的会员，但对订阅《走向数据科学》感兴趣，只是为了阅读类似的教程和文章，单击此处注册成为会员。注册这个链接意味着我可以通过向你推荐 Medium 来获得报酬。
模拟曼卡拉:当我把这个游戏推到极限时会发生什么？
我用 MATLAB 模拟了 10 万个曼卡拉游戏，看看会发生什么。结果如下:
图片由作者提供
简介:
由于曼卡拉是世界上最古老的游戏之一，所以选择这款游戏作为我的下一个棋盘游戏 MATLAB 项目是有意义的。曼卡拉是一个零和游戏，意思是一个人的收获直接导致另一个人的损失。这使得一次模拟几十万个游戏变得非常容易。
曼卡拉规则:
这是一个非常简单的游戏。每个玩家有 6 个坑和一个商店，棋盘上总共有 12 个坑和 2 个商店。每个坑开始时有 4 粒种子。目标是在您的商店中获得最多的种子。每回合，一名玩家从他们的 6 个坑中选择一个，取出所有的种子。逆时针方向移动，玩家从手中丢下一粒种子到下一个坑或仓库，直到手中的种子用完。
图片由pexels.com提供
还有一些特殊的规则需要注意:
- 如果你的最后一粒种子在你的店里用完了，你可以再去一次。
- 如果你的最后一粒种子落在你这边的一个完全空的坑里，你可以直接收集棋盘上的所有种子。
当一方完全空了，游戏结束，另一方从他们的坑里收集剩余的种子，放入他们的仓库。商店被计算出来，种子最多的人获胜。
编码 Mancala:
这种模拟的工作原理是让两个玩家互相随机移动，直到一个人获胜。随机选取一个有种子的坑，这些种子沿着棋盘逆时针分布(见下文)。
while gameover == 0 % while nobody has wonwhile p1_turn == 0 % while still player 1 turnif sum(layout(1:6)) ~= 0 % while player 1 still has seeds in pitsmoves = moves + 1; % plus 1 for total moveswhile good_pick1 == 0 % loop to find a pit with seedsp1 = randi(6); % picks random pitamount = layout(p1); % finds amount of seeds in that pitif amount > 0 % if pit has seedsbreak % breaks out of loop (found a good pit)endend
在每场比赛结束时，感兴趣的数据点被记录下来，并一直进行到所有 100，000 次模拟完成。
move_data(repeats) = moves; % adds moves to move datasetp1_data(repeats) = mean(p1_moves); % adds average move to player 1 move datasetp2_data(repeats) = mean(p2_moves) - 7; % adds average move to player 1 move datasetscore1_data(repeats) = layout(7); % adds player 1 score to player 1 score datasetscore2_data(repeats) = layout(14); % adds player 2 score to player 2 score dataset
结果是:
一个标准的曼卡拉游戏是什么样子的？完成一局游戏平均需要多少回合？最短的游戏是什么？所有这些问题都可以通过 10 万次曼卡拉模拟来回答。
曼卡拉(每坑 4 粒种子)的平均游戏需要 41.42 步才能完成一局。找到的最长的游戏是 77 步，最短的是 11 步(见下面的柱状图)。
图片由作者提供
最短可能的游戏发生在一个玩家很快用完他们所有的种子，然后游戏结束，另一个玩家收集剩下的种子。
每个玩家的获胜概率是多少？先走对你的胜算有影响吗？嗯，开始时每个坑有 4 粒种子，先走的玩家有 48.624%的机会获胜，而另一个玩家有 44.604%的机会获胜(剩下的 6.772%的机会是平局)。
开始游戏时有没有一个最好的第一坑可以选择，以增加自己的胜算？如果你是第一个出手的玩家，选择 6 号坑(最右边的坑)将使你平均每局获得 24.9033 粒种子，比第二高的坑平均多 0.3 粒种子。第二个玩家也一样。第 6 号坑将让他们平均每场比赛获得 23.9039 粒种子。
将曼卡拉推向极限:
开始时每个坑有 100 个种子的游戏是什么样子的？1000 呢？这对胜负概率有什么影响？
我对下面显示的每个初始种子数量(1-10、15、20、50 和 100)运行了 100，000 次模拟，以查看获胜百分比如何变化。
图片由作者提供
随着首发种子数量的增加，比赛平局的几率下降。这是有道理的，因为随着种子数量的增加，商店最终拥有准确数量种子的几率越来越低。
这里值得注意的是获胜概率。第一次行动的玩家和第二次行动的玩家之间的差异似乎变小了，这意味着随着更多种子的加入，游戏变得更加公平。当每个坑有 10 粒种子时，这种差异又跳回来重新开始。这意味着，如果你想要最公平的比赛，每个坑 9 个种子是理想的设置(100 个种子也有一点点差别，但这种游戏将永远持续下去)。
我还想看看开始种子的数量如何影响完成游戏的总平均移动。
图片由作者提供
每场比赛的平均移动不是线性增加的，而是以对数的方式增加的。我将其扩展到 1000 个种子和 5000 个种子，找到了一个 R 值为 0.9715 的最佳拟合方程(如下所示)。
图片由作者提供
假设你从 3 个或更多的种子开始，最佳拟合方程可以用来找出一个游戏在任何初始种子数量下大概要走多少步。
结论:
当您以不同的方式调整初始设置时，Mancala 会出现一些有趣的模式。先走一步可以大大增加你获胜的机会，如果你想降低比赛以平局结束的几率，就多加些种子。观察棋盘游戏的极限以及游戏在这些时刻的反应是非常有趣的。
如果您对这类内容感兴趣，我推荐您阅读我的其他一些关于用 MATLAB 模拟棋盘游戏的文章:
模拟垄断:使用 MATLAB 寻找最佳属性
我用 MATLAB 模拟了 1000 万次大富翁游戏来寻找最好的属性
图片由作者提供
简介:
“大富翁”有 40 个空间和 28 种不同的属性，起初看起来完全是随机的，但在游戏过程中，某些属性可能对你更有益。我想通过 MATLAB 模拟 10，000，000 次掷骰子来找出这些属性。
工作原理:
我用随机数发生器模拟两个骰子滚动，然后将两个值相加。然后，我将该金额添加到玩家的当前位置，以获得他们的新位置。我还有一个变量来跟踪一行中的 doubles 数量。如果玩家连续掷出 3 个双打，他们会因为超速被送进监狱。
die1 = randi(6); % roll first diedie2 = randi(6); % roll second dieif die1 == die2 % if same valuedoubles = doubles + 1; % increase number of doubles in a rowelsedoubles = 0; % set doubles to 0 if different dice valuesend
接下来的代码是所有的特殊情况。例如，如果位置加起来大于 40，它将返回到起点作为参考，因为棋盘是一个循环。位置值 43 将变为 3。此外，如果玩家降落在位置 31(“去监狱”空间)，他们将被发送到位置 11(监狱)。
if doubles < 3 % if less than 3 doubles in a rowposition = position + die1 + die2; % finds new positionif position > 40 % if above 40position = position - 40; % reverts position back to starting position as a referencearound = around + 1; % indicates the player went around the boardendif position == 31 % if player lands on the "go to jail" spacelands(position) = lands(position) + 1; % indicates player landed on position 31position = 11; % sends player to jailend
其他特殊情况的例子是公益金空间和机会空间。16 张社区公益金卡中，只有 1 张改变了玩家的位置(坐牢)。在 16 张机会卡中，有 9 张会改变玩家的位置:
- 进监狱
- 前进到木板路
- 前进到伊利诺斯大街
- 前进到圣查尔斯广场
- 前进到雷丁铁路
- 前进到最近的铁路
- 前进到最近的设施
- 继续前进
- 后退 3 格
在每一轮结束后，玩家结束的空间被记录在所有 40 个点的数组中，然后每个点被模拟的总量(1000 万)除。这给了我们所有 40 个空间着陆的百分比可能性。
结果:
掷骰子 10，000，000 次后，每个空格落在上面的可能性如下所示:
图片由作者提供
下面显示的是上面条形图中的前 20 个电路板空间:
图片由作者提供
虽然一个属性被登陆的次数确实会影响它给玩家带来的利益，但我们也需要考虑购买该属性的成本，以及玩家每次登陆时会付给你多少钱。
上面显示的概率也可以描述为玩家在每次掷骰子后落在该地产上的机会。我们可以利用它，以及每处房产的租金，来计算如果你拥有该房产，你在每个玩家的掷骰子中将会赚到的预期金额。
例如，Ventnor Ave(如上所示)在每一次掷骰中有 2.691%的几率被击中。Ventnor Ave 没有房子的租金是 22 美元，这意味着每一次掷骰子的结果是所有者预期的 0.59202 美元。
每一卷的收入并不是比较房产的最佳指标，因为一个比 Ventnor 更贵的房产应该产生更多的收入，但事实并非总是如此。如果我们将购买房产的初始成本除以每卷的收入，我们可以找到一笔房产需要多少卷才能达到收支平衡。例如，Ventnor 的购买成本为 260 美元，因此所有者应该预计在 439.18 卷后可以收回所有的钱。不需要太多投资就能达到收支平衡的房产会被认为是更好的投资，因为你可以在其他玩家之前开始盈利。
我们还需要考虑在地产上建造房屋和酒店会如何影响盈亏平衡点。下面是一个垄断委员会上的每个属性的细分，给出了关于房屋和酒店的情况(不包括铁路和公用事业)。每个单元格中的值表示平均收回在该资产上花费的资金所需的滚动量。
图片由作者提供
正如你所看到的，从长远来看，有些房产几乎肯定会让你赔钱。购买地中海大道而不在上面建造任何房屋，在整个游戏过程中肯定会让你损失金钱。这不会是一个显著的数额(只需 60 美元购买)，但它仍然不是最好的购买。
所有没有房子的房产都需要非常长的时间才能达到收支平衡，这就是为什么买房子和酒店如此重要。
图片由作者提供
上面的图显示了如何为你拥有的每一个垄断权购买至少 3 个房子，可以大大减少回本所需的滚动次数。即使是开始时是可怕投资的棕色房产(上面的棕色线),在有 3 个或更多房子的游戏中也可以开始盈利。
值得注意的另一点是，随着房产变得越来越贵，3 套房子、4 套房子和一家酒店之间的变化越来越小。以浅蓝色和黄色属性为例:
作者提供的图像(每个单元格中的值代表平均收回在该属性颜色上花费的钱所需的卷数)
淡蓝色的房产在 3 至 4 套房子之间，以及 4 套房子和一家酒店之间的价格有相当大的下降。黄色属性显示最右边的 3 列之间几乎没有变化。这意味着花在 3 套房子上的钱和花在 4 套房子上的钱一样多，但是收支平衡后的收益会更大(对于更贵的房产)。
但是铁路和公用事业呢？
下面显示了铁路，以及您拥有的铁路的所有不同排列。每个单元格中的值表示平均收回在一处或多处房产上花费的钱所需的滚动量。
图片由作者提供
对于仅仅一条铁路的所有权，最好拥有 B&O，因为它会登陆更多。如果你想拥有 2 个，那就去雷丁和 B&O 组合吧。对于 3，寻找阅读，宾夕法尼亚州和 B&O。总的来说，铁路的最佳投资要么获得 3 或全部 4，只是因为人们开始在他们的财产上建造房屋，你的 1 或 2 铁路将不会比游戏早期赚得更多。
至于公用设施:
图片由作者提供
在游戏早期拥有公用事业可能是最好的举措之一。与其他没有房子的房产相比，这两处房产的盈亏平衡点最低。然而，就像铁路一样，你不能在上面建房子，所以随着游戏的进行，这些变得越来越没用。考虑到其他 22 处房产中有 21 处房产的盈亏平衡点低于 3 处房产，即使同时拥有这两处房产，它们的排名仍然垫底。
结论:
橙色的属性是游戏中最好的点。橘子园位于最好的 5 处房产内，有 3 栋房子，是酒店的前 3 名。orange properties New York Ave 和 Tennessee Ave 在最受欢迎的景点中排名第三和第五。
图片由作者提供
木板路和公园广场(深蓝色)是棘手的。Boardwalk 在盈亏平衡的滚动次数方面一直名列前茅，而 Park Place 总是垫底。与橙子相比，这种垄断是一种高风险、高回报的选择。
绝对最差的酒店群体肯定是格林一家。这三家酒店的盈亏平衡营业额都排在后 5 名。就登陆的可能性而言，它们也低于平均水平。
图片由作者提供
只有当你拥有 3 个或者全部 4 个的时候，铁路才会在游戏接近尾声的时候变得有利。如果你能早点得到它们，并开始收支平衡，那么这是一个伟大的举动。公用事业是最好的财产，直到人们开始建造房屋。仅仅依靠公用事业或铁路是无法取胜的，所以不要把所有的精力和金钱都投入其中。
总的来说，这个模拟表明，有一些策略可以应用到你的下一个垄断游戏中。如果你对这种类型的内容感兴趣，我推荐阅读我的其他棋盘游戏文章，我教了一台机器如何玩 Connect 4 ，以及最好和最差的卡坦棋盘设置。
用 Python 模拟多智能体群集事件
以鸟类迁徙为例模拟多主体强化学习分析中的群集现象
成群结队在自然界频繁发生:从一群候鸟到一群执行特定任务的无人机。因此，能够模拟群集如何工作将使您能够构建自己的多智能体强化学习算法，以密切模仿现实世界的场景。为了让你对我们今天正在构建的东西有所了解，下面的 GIF 将是我们使用候鸟为例的群集模拟模型的输出。事不宜迟，我们开始吧！
作者 GIF
**注:**可以在帖子末尾找到完整代码。
动机
为什么模拟在多智能体强化学习(RL)问题中很重要？这主要是因为模拟是 RL 中的探索性数据分析(EDA)方法之一，它允许我们可视化系统中不同代理的相互关系。在这篇文章中，我们将讨论群集事件，因为它们是多智能体强化学习的许多例子中的一些，我在下面编译了它们在现实生活中是如何行动的！
群集算法:Viscek 模型
我们将实现 Viscek 模型(1995 ),作为我们的多智能体群集基线算法，在此基础上已经构建了许多更高级的实现。
想象你正在试图模拟 300 只鸟的群集( N =300)，其中每只鸟的索引为 i=0…( N -1)。一只鸟只有在半径为 R 的圆内与另一只鸟相互作用时，才会改变方向和速度矢量。
交互发生的半径(来源: Viscek )
让我们从导入所需的 Python 库并初始化上述常量开始。
import matplotlib.pyplot as plt
import numpy as npN = 300 # number of birds
R = 0.5 # radius of effect# Preparing the drawing
L = 5 # size of drawing boxfig = plt.figure(figsize=(6,6), dpi=96)
ax = plt.gca()
对于每个时间步， n ，每只鸟沿 x 轴和 y 轴( r )的位置更新如下:
每个时间步长的位置更新(来源: Viscek
这是用 Python 实现的。
Nt = 200 # number of time steps
dt = 0.1 # time step# bird positions (initial)
x = np.random.rand(N,1)*L
y = np.random.rand(N,1)*L
v 为速度恒定的速度矢量， v0 ，角度θᵢ(读作:theta) *。*我们还将引入一个角度扰动(读作: eta )来模拟鸟类运动的自然随机性。
时间步长 n 时 x 和 y 方向的速度矢量计算(来源: Viscek
速度向量在下面的 Python 中初始化。
v0 = 0.5 # constant speed
eta = 0.6 # random fluctuation in angle (in radians)
theta = 2 * np.pi * np.random.rand(N,1)# bird velocities (initial)
vx = v0 * np.cos(theta)
vy = v0 * np.sin(theta)
现在我们将进入主 for 循环。对于每个时间步长，位置r更新为:
# Update position
x += vx*dt
y += vy*dt
并且通过找到 R 内的(1) 平均角度，(2) 加上随机扰动以产生新的θ，以及(3)更新新的速度向量**、 v 来更新v😗*
(1)R 内的平均角度
mean_theta = theta# Mean angle within R
for b in range(N):
neighbors = (x-x[b])**2+(y-y[b])**2 < R**2
sx = np.sum(np.cos(theta[neighbors]))
sy = np.sum(np.sin(theta[neighbors]))
mean_theta[b] = np.arctan2(sy, sx)
(2)给速度角增加一个随机扰动。随机项是通过从[-0.5，0.5]内的正态分布中采样，将其乘以预定义的η，并将该扰动项与计算的平均θ相加而生成的。
正态分布(作者图片)
# Adding a random angle perturbation
theta = mean_theta + eta*(np.random.rand(N,1)-0.5)
(3)用新的θ更新速度矢量:
# Updating velocity
vx = v0 * np.cos(theta)
vy = v0 * np.sin(theta)
总的来说， r 和 v 的更新将在 dt * Nt 个时间步长内完成。
最后，现在一切都准备好了，让我们画出更新的样子来结束吧！
plt.cla()
plt.quiver(x,y,vx,vy,color='r')
ax.set(xlim=(0, L), ylim=(0, L))
ax.set_aspect('equal')
plt.pause(0.001)
如果一切正常，你会在屏幕上看到这个消息。恭喜你！您可以试验并更改一些参数，看看您的模拟模型如何变化。
蜂拥事件(图片由作者提供)
如果你已经到达终点，那就做得很好。你可以在这里找到完整的代码:
离别的思绪
模拟是强化学习中探索性数据分析的一个步骤。一旦您对现实世界的数学行为有了更清晰的认识，您将能够设计更鲁棒的 RL 算法。如果你喜欢这篇文章，我会做更多真实世界的模拟，作为我的多代理 RL EDA 系列的一部分。敬请关注！
【https://tinyurl.com/2npw2fnz】订阅我的电子邮件简讯:https://tinyurl.com/2npw2fnz我定期在那里用通俗易懂的英语和漂亮的可视化总结 AI 研究论文。
最后，您可以在这里找到更多关于 RL 和仿真主题的资源:
** https://github.com/AndreasKuhnle/SimRLFab https://deepsense.ai/what-is-reinforcement-learning-the-complete-guide/ https://www.guru99.com/reinforcement-learning-tutorial.html **
用 SimPy 在 Python 中模拟和可视化真实事件
我们将介绍活动行业完整模型的开发过程，并展示三种不同的可视化结果的方法(包括 AR/VR)
离散事件仿真(DES)已成为专门产品的领域，如 Simulink 8[1]和 MatLab/Simulink [2]。然而，在用 Python 进行分析的时候(过去我会使用 MatLab ),我很想测试 Python 是否也有 DES 的答案。
DES 是一种使用统计函数对现实生活事件建模的方法，通常用于医疗保健、制造、物流和其他领域的队列和资源使用[3]。最终目标是获得关键的运营指标，如资源使用和平均等待时间，以便评估和优化各种实际配置。SIMUL8 有一个描述如何模拟急诊室等待时间的视频[4]，MathWorks 有许多教育视频来提供该主题的概述[5]，此外还有一个关于汽车制造的案例研究[6]。SimPy [7]库支持在 Python 中描述和运行 des 模型。与 SIMUL8 等软件包不同，SimPy 不是一个用于构建、执行和报告模拟的完整图形环境；但是，它确实提供了执行模拟和输出数据以进行可视化和分析的基本组件。
本文将首先介绍一个场景，并展示如何在 SimPy 中实现它。然后，它将查看可视化结果的三种不同方法:Python 原生解决方案(使用 Matplotlib [8]和 Tkinter [9])、基于 HTML5 canvas 的方法和交互式 AR/VR 可视化。最后，我们将使用 SimPy 模型来评估替代配置。
场景
在我们的演示中，我将使用我以前工作中的一个例子:一个活动的入场队列。然而，遵循类似模式的其他示例可能是在杂货店或接受在线订单的餐馆、电影院、药店或火车站排队。
我们将模拟一个完全由公共交通服务的入口:定期会有一辆公共汽车让几个顾客下车，然后他们需要在进入活动前扫描他们的门票。一些参观者会有他们预先购买的徽章或门票，而其他人则需要先去售票亭购买门票。更复杂的是，当参观者走近售票亭时，他们会成群结队(模拟家庭/团体购票)；但是，每个人都需要单独扫描门票。
下面描述了该场景的高级布局。
图片作者和 Matheus Ximenes
为了模拟这种情况，我们需要决定如何使用概率分布来表示这些不同的事件。在我们的实施中所做的假设包括:
- 平均每三分钟就有一辆公共汽车到达。我们将使用λ为 1/3 的指数分布来表示这一点
- 每辆巴士将包含 100 +/- 30 名游客，使用正态分布确定(μ = 100，σ = 30)
- 游客将使用正态分布(μ = 2.25，σ = 0.5)组成 2.25+/–0.5 人的群体。我们将把它四舍五入到最接近的整数
- 我们假设 40%的访问者需要在销售亭购买门票，另外 40%的访问者将带着已经在线购买的门票到达，20%的访问者将带着员工凭证到达
- 参观者平均需要一分钟走出巴士，走到卖家摊位(正常，μ = 1，σ = 0.25)，再花半分钟从卖家走到扫描仪(正常，μ = 0.5，σ = 0.1)。对于那些跳过卖家的人(预先购买的门票或佩戴徽章的员工)，我们假设平均步行时间为 1.5 分钟(正常，μ = 1.5，σ = 0.35)
- 游客到达时会选择最短的队伍，每一队都有一个卖家或扫描仪
- 销售需要 1 +/- 0.2 分钟完成(正常，μ = 1，σ = 0.2)
- 完成一次扫描需要 0.05 +/- 0.01 分钟(正常，μ = 0.05，σ = 0.01)
记住这一点，让我们从输出开始，并从那里向后工作。
作者图片
左手边的图表代表每分钟到达的游客数量，右手边的图表代表在该时刻离开队列的游客在被接待之前需要等待的平均时间。
简单模拟设置
包含完整可运行源代码的存储库可以在https://github.com/dattivo/gate-simulation找到，下面的片段摘自 simpy example.py 文件。在本节中，我们将逐步完成简单的特定设置；但是，请注意，为了将重点放在 SimPy 的 DES 特性上，省略了连接到 Tkinter 进行可视化的部分。
首先，让我们从模拟的参数开始。分析起来最有趣的变量是销售者行数( SELLER_LINES )和每行的销售者数( SELLERS_PER_LINE )以及它们对于扫描器的等价物( SCANNER_LINES 和 SCANNERS_PER_LINE )。此外，请注意两种可能的队列/卖家配置之间的区别:尽管最普遍的配置是拥有多个不同的队列，访问者将选择并停留，直到他们得到服务，但在零售业中看到多个卖家排队也变得更加主流(例如，在一般商品大盒子零售商的快速结账队列)。
配置完成后，让我们开始 SimPy 过程，首先创建一个“环境”、所有队列(资源)，并运行模拟(在本例中，直到第 60 分钟)。
请注意，我们正在创建一个 RealtimeEnvironment ，它旨在近乎实时地运行一个模拟，特别是为了在它运行时将它可视化。随着环境的建立，我们生成我们的销售者和扫描器线路资源(队列),然后我们将依次传递给我们的总线到达的“主事件”。 env.process() 命令将开始流程，如下图的 bus_arrival() 函数所述。此函数是顶级事件，所有其他事件都从该事件调度。它模拟每隔 BUS_ARRIVAL_MEAN 分钟到达一辆公交车，车上有 BUS_OCCUPANCY_MEAN 人，然后相应地触发销售和扫描过程。
由于这是顶级事件函数，我们看到该函数中的所有工作都发生在一个无限的 while 循环中。在循环中，我们用 env.timeout() 来“放弃”我们的等待时间。SimPy 广泛使用了生成器函数，这些函数将返回生成值的迭代器。关于 Python 生成器的更多信息可以在[10]中找到。
在循环的最后，我们将分派两个事件中的一个，这取决于我们是直接进入扫描仪，还是随机决定该组需要先购买门票。请注意，我们不会屈服于这些过程，因为这将指示 SimPy 按顺序完成这些操作中的每一个；取而代之的是，所有下车的游客将同时排队。
请注意，使用了 people_ids 列表，以便为每个人分配一个唯一的 ID 用于可视化目的。我们使用 people_ids 列表作为待处理人员的队列；随着游客被分派到他们的目的地，他们从 people_ids 队列中被移除。
purchasing_customer() 函数模拟三个关键事件:走向队列、排队等候，然后将控制传递给 scanning_customer() 事件(与 bus_arrival() 调用的函数相同，用于那些绕过销售者直接进入扫描器的事件)。这个函数根据选择时最短的线来选择它的线。
最后，我们需要实现 scanning_customer() 的行为。这与 purchasing_customer() 函数非常相似，但有一个关键区别:尽管游客可能会成群结队地到达和步行，但每个人都必须单独扫描他们的门票。因此，您将看到对每个被扫描的客户重复扫描超时。
我们将步行持续时间和标准偏差传递给 scanning_customer() 函数，因为这些值会根据访问者是直接走到扫描仪前还是先在卖家处停留而变化。
使用 Tkinter(原生 Python UI)可视化数据
为了可视化数据，我们添加了一些全局列表和字典来跟踪关键指标。例如，到达字典按分钟跟踪到达人数，而 seller_waits 和 scan_waits 字典将模拟的分钟映射到在这些分钟内退出队列的人的等待时间列表。还有一个 event_log 列表，我们将在下一节的 HTML5 Canvas 动画中使用。当关键事件发生时(例如，一个访问者退出队列)，调用 simpy example.py 文件中 ANALYTICAL_GLOBALS 标题下的函数来保持这些字典和列表最新。
我们使用一个辅助的 SimPy 事件向 UI 发送一个 tick 事件，以便更新时钟、更新当前的平均等待时间并重新绘制 Matplotlib 图表。完整的代码可以在 GitHub 储存库中找到(https://GitHub . com/dattivo/gate-simulation/blob/master/simpy % 20 example . py)；但是，下面的代码片段提供了如何从 SimPy 发送这些更新的框架视图。
使用标准 Tkinter 逻辑来表示用户进出销售者和扫描者队列的可视化。我们创建了 QueueGraphics 类来抽象销售者和扫描者队列的公共部分。这个类中的方法被编码到上一节中描述的 SimPy 事件函数中，以更新画布(例如， sellers.add_to_line(1) ，其中 1 是卖方编号，以及sellers . remove _ from _ line(1))。作为未来的工作，我们可以在流程的关键点使用事件处理程序，这样 SimPy 模拟逻辑就不会与特定于该分析的 UI 逻辑紧密耦合。
使用 HTML5 画布制作数据动画
作为另一种可视化方式，我们希望从 SimPy 模拟中导出事件，并将它们放入一个简单的 HTML5 web 应用程序中，以便在 2D 画布上可视化场景。我们通过在简单事件发生时添加一个 event_log 列表来实现这一点。特别是，公交车到达、步行到销售者、在销售者队伍中等待、购买门票、步行到扫描仪、在扫描仪队伍中等待和扫描门票事件都被记录为单独的字典，然后在模拟结束时导出到 JSON。你可以在这里看到一些示例输出:https://github . com/dattivo/gate-simulation/tree/master/output
我们开发了一个快速的概念验证来展示如何将这些事件转化成 2D 动画，你可以在 https://dattivo.github.io/gate-simulation/visualize.html 进行实验。你可以在https://github . com/dattivo/gate-simulation/blob/master/visualize . ts中看到动画逻辑的源代码。
作者图片
这种可视化得益于动画，然而，出于实际目的，基于 Python 的 Tkinter 接口组装起来更快，Matplotlib 图形(可以说是该模拟最重要的部分)在 Python 中设置起来也更平滑和更熟悉。也就是说，让行为生动起来是有价值的，尤其是在向非技术利益相关者传达结果的时候。
使用虚拟现实制作数据动画
让画布动画更进一步， Matheus Ximenes 和我一起工作，使用 HTML5 画布也在使用的相同 JSON 模拟数据，将以下 AR/VR 3-D 可视化放在一起。我们使用我们已经熟悉的 React [11]和 A-FRAME [12]实现了这一点，A-FRAME 令人惊讶地易于使用和学习。
你可以在 https://www.dattivo.com/3d-visualization/的亲自测试模拟
分析卖方/扫描仪队列配置备选方案
虽然这个例子是为了演示如何创建和可视化 SimPy 模拟，但是我们仍然可以展示几个例子来说明平均等待时间如何依赖于队列的配置。
让我们从上面动画中演示的案例开始:六个销售人员和四个扫描仪，每行一个销售人员和一个扫描仪(6/4)。60 分钟后，我们看到卖家的平均等待时间为 1.8 分钟，而扫描仪的平均等待时间为 0.1 分钟。从下面的图表中，我们看到销售者的时间在大约 6 分钟的等待时达到峰值。
我们可以看到，卖家们都是一致备份的(虽然 3.3 分钟可能不会太不合理)；所以，让我们看看，如果我们再增加四个卖家，总数达到 10 个，会发生什么。
正如预期的那样，卖家的平均等待时间减少到了 0.7 分钟，最长等待时间减少到了 3 分钟多一点。
现在，让我们说，通过降低在线门票的价格，我们能够将持票到达的人数提高 35%。最初，我们假设所有访客中有 40%需要购买门票，40%已经在线预购，20%是凭凭证进入的工作人员和供应商。因此，随着 35%以上的人带着票到达，我们将需要购买的人数减少到 26%。让我们用初始的 6/4 配置来模拟一下。
在这种情况下，卖家的平均等待时间减少到 1.0 分钟，最长等待时间超过 4 分钟。在这种情况下，增加 35%的在线销售额与增加更多的卖家队列有相似的效果；如果等待时间是我们最想减少的指标，那么在这一点上，我们可以考虑这两个选项中哪一个具有更强的业务案例。
结论和未来工作
Python 可用的数学和分析工具的广度令人生畏，SimPy 完善了这些功能，还包括离散事件模拟。与 SIMUL8 等商业打包工具相比，Python 方法确实留给编程更多的空间。对于快速分析来说，组装模拟逻辑和从头构建 UI 和度量支持可能是笨拙的；然而，它确实提供了很大的灵活性，对于已经熟悉 Python 的人来说应该相对简单。如上所述，SimPy 提供的 DES 逻辑产生了清晰易读的代码。
如前所述，Tkinter 可视化是三种演示方法中最直接的一种，尤其是在包含 Matplotlib 支持的情况下。HTML5 canvas 和 AR/VR 方法可以方便地将可共享的交互式可视化放在一起；然而，他们的发展并不平凡。
比较队列配置时需要考虑的一个重要改进是销售者/扫描者利用率。减少排队时间只是分析的一个组成部分，因为在得出最佳解决方案时，还应该考虑销售者和扫描仪空闲时间的百分比。此外，如果有人看到队列太长而选择不进入，添加一个概率也会很有趣。
参考
[2]https://www . mathworks . com/solutions/discrete-event-simulation . html
[3]https://en.wikipedia.org/wiki/Discrete-event_simulation
https://www.simul8.com/videos/
[5]https://www . mathworks . com/videos/series/understanding-discrete-event-simulation . html
https://simpy.readthedocs.io/en/latest/
[9]https://docs.python.org/3/library/tkinter.html
【https://wiki.python.org/moin/Generators
[11]https://reactjs.org/
(本文改编自之前发表的两篇文章:https://dattivo . com/simulating-real-life-events-in-python-with-simpy/和https://dattivo . com/3d-visualization-of-a-simulated-real-life-events-in-virtual-reality/)
在线性建模中模拟多重共线性的影响
使用 R、purrr 和并行计算
在数据科学中，使用线性回归有时会被低估。作为 t 检验和方差分析等非常基本的统计概念的概括，它与统计学领域有着很深的联系，可以作为解释数据方差的强大工具。由于线性模型只有参数是线性的，所以也可以描述多项式甚至一些乘法关系。
但是，由于参数的性质，线性回归也更容易受到极值和数据中多重共线性的影响，我们希望通过模拟对后者进行更详细的分析。在许多其他问题中，多重共线性可能产生不稳健的模型(交错系数),并可能破坏统计显著性。
在运行模拟之前，我们将研究多重共线性对简单线性模型的影响。我们将使用来自 Stat Trek 的数据集，其中包含学生的考试成绩和一些其他表现指标(如智商或平均绩点):
我们观察到所有变量之间的高度相关性，尤其是iq和gpa相关性非常强:
GGally::ggpairs(data)
成对测量图—图片由作者完成
如果我们使用iq和study_hours拟合一个简单的模型来解释测试分数，我们可以看到它们的影响在统计上是显著的(p 值约为 0.03，总体显著性):
simple_model <- lm(test_score ~ iq + study_hours, data = data)
sqrt(mean((simple_model$residuals)^2))
## [1] 3.242157
summary(simple_model)
##
## Call:
## lm(formula = test_score ~ iq + study_hours, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.341 -1.130 -0.191 1.450 5.542
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.1561 15.9672 1.450 0.1903
## iq 0.5094 0.1808 2.818 0.0259 *
## study_hours 0.4671 0.1720 2.717 0.0299 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.875 on 7 degrees of freedom
## Multiple R-squared: 0.9053, Adjusted R-squared: 0.8782
## F-statistic: 33.45 on 2 and 7 DF, p-value: 0.0002617
如果我们将gpa加入到模型中，我们预计会得到一个更小的(样本内)误差，但截距和斜率的估计值会有很大不同。另外，iq和gpa不再重要。这是由于它们的共线性，在这里有很好的描述。
# estimates and p values very different
full_model <- lm(test_score ~ iq + study_hours + gpa, data = data)sqrt(mean((full_model$residuals)^2))
## [1] 3.061232summary(full_model)
##
## Call:
## lm(formula = test_score ~ iq + study_hours + gpa, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.3146 -1.2184 -0.4266 1.5516 5.6358
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 50.30669 35.70317 1.409 0.2085
## iq 0.05875 0.55872 0.105 0.9197
## study_hours 0.48876 0.17719 2.758 0.0329 *
## gpa 7.37578 8.63161 0.855 0.4256
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.952 on 6 degrees of freedom
## Multiple R-squared: 0.9155, Adjusted R-squared: 0.8733
## F-statistic: 21.68 on 3 and 6 DF, p-value: 0.001275BIC(full_model, simple_model)
## df BIC
## full_model 5 62.26804
## simple_model 4 61.11389
在正常的预测建模场景中，我们将通过使用正则化来自动进行特征选择，从而在保持良好的可预测性的同时，提高我们模型的泛化能力和鲁棒性，从而尝试避免这种情况。在经典统计学中，我们不仅关心可预测性(低误差)，还关心变量的统计显著性(p 值)和因果关系。
这里的传统方法是使用方差膨胀因子从等式中移除变量。我们不仅仅着眼于预测因子的成对相关性，而是试图着眼于预测因子的线性相关性:我们本质上试图通过所有其他预测因子来估计一个预测因子，然后方差膨胀因子作为未解释方差的倒数给出(即1 / (1 - R.sq.))。如果预测值可以由其他预测值线性表示，则该值较高，否则该值较低。根据经验，如果它高于 4，则该预测值可能会被丢弃。
vif_gpa_model <- lm(gpa ~ iq + study_hours, data = data)
vif_gpa_factor <- 1 / (1 - summary(vif_gpa_model)$r.squared)
vif_gpa_factor
## [1] 19.65826
我们可以很容易地对所有的预测值都这样做，并决定去掉高膨胀系数的特征。
car::vif(full_model)
## iq study_hours gpa
## 22.643553 2.517786 19.658264
领域知识可能会告诉我们，GPA 是智商的结果，而不是智商的结果，因此将智商作为一个独立变量可能比使用代理变量 GPA 更合理。
模拟 1:用 purrr 运行“许多模型”
为了了解数据中多重共线性的有害影响，让我们模拟一下当我们重复拟合具有不同特征的线性回归时，模型系数可能会发生什么情况。
为了做到这一点，我们将使用purrr，它帮助我们在数据帧列表的顶部以矢量化的方式运行我们自己的函数，并因此在数据和/或模型上运行(推荐阅读:R 中的用于数据科学以及官方文档和purrr的)。
为了理解我们的模拟，让我们首先模拟从简单的线性关系y = 5 + 3 * x中提取一些虚拟数据，并多次拟合线性模型。我们创建了一些辅助函数来绘制可重复的样本并拟合模型:
generate_dummy_data <- function(data, seed = NULL) {
set.seed(seed = seed)
data <- dplyr::tibble(x = rnorm(100, 50, 10),
y = 5 + 3 * x + rnorm(100, 0, 0.1))
return(data)
}generate_dummy_model <- function(data) {
model <- lm(y ~ x, data = data)
return(model)
}
我们现在可以多次运行该模拟，并使用purrr::map提取模型系数(截距和斜率):
我们看到估计值非常稳定，接近真实值。
dummy_simulation %>%
tidyr::unnest(coefs) %>%
ggplot(aes(x = estimate)) +
geom_histogram(bins = 30) +
facet_wrap(~ term, scales = "free") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
非常简单的线性关系的估计分布—图片由作者完成
现在，让我们通过使用自举样本，将此应用于我们的测试分数数据:
我们不仅要看一个模型，还要比较两个模型:一个是使用所有特征的完整模型(即存在多重共线性),另一个是使用简化特征集的模型:
让我们绘制两个模型中截距和斜率的估计值(使用ggplot2和theme_set(theme_linedraw()))。我们看到截距的巨大差异，还有iq-系数:
model_coefs <-
simulations %>%
dplyr::select(name, coefs) %>%
tidyr::unnest(coefs)model_coefs %>%
ggplot(aes(x = estimate, color = name)) +
geom_freqpoly(bins = 30) +
facet_wrap( ~ term, scales = "free")
完整和简单模型的系数分布——作者完成的图像
p 值看起来也不是很有希望，模型很难获得关于iq和gpa相关性的信心，p 值没有真正集中在低端:
model_coefs %>%
ggplot(aes(x = p_value, color = name)) +
geom_freqpoly(bins = 30) +
facet_wrap( ~ term) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
完整和简单模型的 p 值分布—图片由作者制作
模拟 2:并行模型评估
让我们在一个稍大的数据集上运行另一个模拟，在这里我们额外计算一些测试数据误差度量。
我们将以稍微不同的方式执行模拟:我们将并行评估不同的模型，并且由于我们只对系数估计和误差度量感兴趣，我们将避免从工人传输太多数据，并且不存储训练数据和模型对象本身。
让我们来看一下鱼市场数据集，它由不同种类的不同鱼类尺寸(长、宽、高)及其重量组成。测量值再次高度相关。我们将尝试预测权重，并在此过程中使用一些简单的领域知识。让我们去除一些明显的异常值，即测量误差。
url <- "[https://tinyurl.com/fish-market-dataset](https://tinyurl.com/fish-market-dataset)"
raw_data <- readr::read_csv(url)
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## Species = col_character(),
## Weight = col_double(),
## Length1 = col_double(),
## Length2 = col_double(),
## Length3 = col_double(),
## Height = col_double(),
## Width = col_double()
## )data <-
raw_data %>%
dplyr::rename_all(tolower) %>%
dplyr::mutate(species = tolower(species)) %>%
dplyr::filter(weight > 0)dplyr::sample_n(data, 5)
## # A tibble: 5 x 7
## species weight length1 length2 length3 height width
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 perch 700 34.5 37 39.4 10.8 6.26
## 2 smelt 8.7 10.8 11.3 12.6 1.98 1.29
## 3 perch 820 37.1 40 42.5 11.1 6.63
## 4 bream 700 30.4 33 38.3 14.9 5.29
## 5 perch 78 16.8 18.7 19.4 5.20 3.12GGally::ggpairs(data) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
所有测量的混沌对图—图片由作者完成
data %>%
ggplot(aes(x = length1, y = weight)) +
geom_point() +
geom_smooth(method = "loess")
## `geom_smooth()` using formula 'y ~ x'
长度 1 与重量—图片由作者制作
由于重量与体积成比例，我们可能会找到类似于weight ~ length^3的关系，因此我们可以拟合log(weight) ~ log(length)，这将导致log(weight) ~ a + b * log(length)的估计值a, b。因此，在求幂运算weight ~ exp(a) * length^b之后，我们期望b大约为 3，而a为某个比例因子:
simple_log_log_model <- lm(log(weight) ~ log(length1), data = data)summary(simple_log_log_model)
##
## Call:
## lm(formula = log(weight) ~ log(length1), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.90870 -0.07280 0.07773 0.26639 0.50636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.62769 0.23481 -19.71 <2e-16 ***
## log(length1) 3.14394 0.07296 43.09 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3704 on 156 degrees of freedom
## Multiple R-squared: 0.9225, Adjusted R-squared: 0.922
## F-statistic: 1857 on 1 and 156 DF, p-value: < 2.2e-16
很公平，b在 3.14 左右。添加height和width可能会再次提高模型性能。在我们开始之前，让我们快速构建另一个没有对数变换的简单模型，通过查看残差图来了解由于非线性关系和不断增加的方差而导致的不良拟合。
naive_model <- lm(weight ~ length1 + width + height + species,
data = data)par(mfrow=c(2,2))
plot(naive_model)
朴素模型的残差诊断—由作者完成的图像
# plot residuals against all features
data %>%
dplyr::mutate(residual = naive_model$residuals) %>%
dplyr::select(-species) %>%
# convert to wide to long for plotting
tidyr::gather(key, value, -residual) %>%
ggplot(aes(x = value, y = residual)) +
geom_point() +
facet_wrap(~ key, scales = "free")
朴素模型的剩余图——作者完成的图像
现在回到双对数模型，我们假设它会表现得更好。
log_log_model <- lm(log(weight) ~ log(length1) + log(width) + log(height), data = data)par(mfrow=c(2,2))
plot(log_log_model)
对数模型的残差诊断—图片由作者完成
data %>%
dplyr::mutate(residual = log_log_model$residuals) %>%
dplyr::select(-species) %>%
tidyr::gather(key, value, -residual) %>%
ggplot(aes(x = value, y = residual)) +
geom_point() +
facet_wrap(~ key, scales = "free")
对数-对数模型的残差图—图片由作者完成
我们在残差中没有看到趋势或严重的漏斗效应。一般来说，当我们预期 x%的增加将导致 y%的响应增加时，双对数模型是有意义的。
我们使用caret::dummyVars对我们的species特性进行一键编码(在lm中，我们也可以使用因子，但是那些不能与glmnet一起工作):
one_hot_encoder <- caret::dummyVars( ~ ., data = data)
data <- dplyr::as_tibble(predict(one_hot_encoder, newdata = data))
我们将实际对我们的特征进行对数变换，而不是使用 R 公式:
numeric_features <- c("length1", "length2", "length3", "height", "width")data <-
data %>%
dplyr::mutate_at(numeric_features, log) %>%
dplyr::rename_at(numeric_features, function(x) paste0(x, "_log"))
x <- data %>% dplyr::select(-weight)
y <- data %>% dplyr::pull(weight)generate_data <- function(seed=NULL, test_size=1/4) {
n <- nrow(x)
set.seed(seed)
train_idx <- sample(1:n, size = (1-test_size) * n)
return(list(x_train = x[train_idx,],
x_test = x[-train_idx,],
y_train = y[train_idx],
y_test = y[-train_idx]))
}
在正常的线性模型之上，我们还将使用glmnet包来拟合正则化的线性模型(山脊、套索和弹性网)。我们创建了一个包装器函数，可以用来传递种子和带有模型参数的命名列表。然后，函数本身将使用上面定义的函数根据种子生成数据。
让我们用一些参数来测试我们的函数:
params <- list(features = c("length1_log", "height_log", "width_log"))fit_glmnet(seed = 1, params = params)
## $coefs
## # A tibble: 4 x 2
## term estimate
## <chr> <dbl>
## 1 (Intercept) -1.95
## 2 length1_log 1.50
## 3 height_log 0.652
## 4 width_log 0.876
##
## $rmse
## [1] 26.5608
用相同的种子和参数再次拟合将产生相同的估计:
fit_glmnet(seed = 1, params = params)
## $coefs
## # A tibble: 4 x 2
## term estimate
## <chr> <dbl>
## 1 (Intercept) -1.95
## 2 length1_log 1.50
## 3 height_log 0.652
## 4 width_log 0.876
##
## $rmse
## [1] 26.5608
我们现在可以使用purrr::map轻松运行多次模拟，并使用purrr::map_dbl提取错误。
simulations <- purrr::map(1:1000, fit_glmnet)
rmse <- purrr::map_dbl(simulations, "rmse")
hist(rmse)
模拟的误差分布—图片由作者完成
现在让我们用不同的参数组合来运行模拟。我们将把我们的参数添加到数据帧中，以便能够在模型和种子上按行并行运行模拟:
params <- list(ridge = list(lambda = 0.1, alpha = 0),
lasso = list(lambda = 0.01, alpha = 1),
elastic_net = list(lambda = 0.01, alpha = 0.5),
full = list(lambda = 0),
simple = list(features = c("length1_log",
"width_log",
"height_log"),
lambda = 0))params_df <- dplyr::tibble(name=factor(names(params)),
params=params)params_df
## # A tibble: 5 x 2
## name params
## <fct> <named list>
## 1 ridge <named list [2]>
## 2 lasso <named list [2]>
## 3 elastic_net <named list [2]>
## 4 full <named list [1]>
## 5 simple <named list [2]>
我们现在将添加和取消种子嵌套:
n_simulations <- 1000
n_workers <- 10simulations <-
params_df %>%
dplyr::mutate(seed = list(1:n_simulations)) %>%
tidyr::unnest(seed)head(simulations)
## # A tibble: 6 x 3
## name params seed
## <fct> <named list> <int>
## 1 ridge <named list [2]> 1
## 2 ridge <named list [2]> 2
## 3 ridge <named list [2]> 3
## 4 ridge <named list [2]> 4
## 5 ridge <named list [2]> 5
## 6 ridge <named list [2]> 6
使用该数据帧，我们现在可以使用furrr::future_map2作为purrr::map2的并行版本，按行并行运行模拟。我们通过将furrr::furrr_options(seed = NULL)传递给函数调用来通知furrr我们已经设置了种子。
tictoc::tic()future::plan(future::multisession, workers = n_workers)simulations <-
simulations %>%
dplyr::mutate(simulation = furrr::future_map2(.x = seed,
.y = params,
.f = fit_glmnet,
.options = furrr::furrr_options(seed = NULL))) %>%
dplyr::mutate(rmse = purrr::map_dbl(simulation, "rmse"),
coefs = purrr::map(simulation, "coefs"))tictoc::toc()
## 6.596 sec elapsed
让我们看看不同模型的误差分布:
simulations %>%
ggplot(aes(x = rmse, color = name)) +
geom_boxplot()
不同模型的误差—图片由作者制作
我们看到，完整模型表现最好，但正则化模型的表现不会差太多。
simulations %>%
dplyr::group_by(name) %>%
dplyr::summarise(rmse_q90 = quantile(rmse, probs = 0.9),
rmse = median(rmse)) %>%
dplyr::ungroup() %>%
dplyr::mutate(perc = rmse / min(rmse)) %>%
dplyr::arrange(perc)
## # A tibble: 5 x 4
## name rmse_q90 rmse perc
## <fct> <dbl> <dbl> <dbl>
## 1 full 72.5 53.4 1
## 2 ridge 74.4 53.6 1.00
## 3 lasso 76.5 53.9 1.01
## 4 elastic_net 74.1 54.2 1.02
## 5 simple 83.9 61.0 1.14
尽管不同模型的系数变化很大，但不同模拟的系数也不同，尤其是完整模型:
model_coefs <-
simulations %>%
dplyr::select(name, coefs) %>%
tidyr::unnest(coefs)model_coefs %>%
ggplot(aes(x = estimate, color = name)) +
geom_freqpoly(bins = 30) +
facet_wrap(~ term, scales = "free", ncol = 2) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
不同模型的系数估计—图片由作者完成。
结论
我们看到了多重共线性如何给线性回归中的参数估计带来麻烦，例如在参数的高方差中寻找表达式或缺乏统计显著性。与此同时，我们看到，这并不一定会导致更差的模型性能，如均方根误差等测试误差指标。
我们还学习了一些处理这些问题的经典技术，例如查看方差膨胀因子、进行简单的特征选择或应用正则化。还有许多其他有用的技术，如降维、逐步或基于树的特征选择。
此外，许多其他更复杂的机器学习模型对这些问题更稳健，因为它们在拟合过程中结合了信息增益和正则化，或者使用不同的数据表示，这些数据表示不太容易受到单次观察和数据内相关性的影响。例如，像树模型这样的非参数候选模型，有时还有像神经网络这样的“极端参数”模型。
但是，由于线性回归允许我们从数据中得出强有力的见解和结论，并且通常作为一个简单的基准模型是有用的，所以理解它的统计主干和它的一些缺点是非常有帮助的。
最初发表于【https://blog.telsemeyer.com】。
模拟 NBA:2021–22 赛季结果
我用蒙特卡罗方法模拟了 1230 场比赛和季后赛。你的队会赢几次？
图片由作者提供
简介:
随着 NBA 2021-22 季前赛的开始，我开始模拟每支球队的 82 场比赛，加上季后赛，以找出联盟最好和最差的表现。
NBA 编码:
类似于我的 NFL 赛季模拟，所有 30 支球队都有一个评分。今年的范围从 1304(克利夫兰骑士队)开始，到 1752(密尔沃基雄鹿队)，其他球队在两者之间。
每个队赢得比赛的机会是由这个等级决定的。以赛季揭幕战为例:主场作战的雄鹿(1752)对篮网(1717)。
图片由作者提供
利用这些等式，雄鹿有 55.20%的胜算，而篮网有 44.80%的胜算。
一个随机数生成器决定获胜者，然后根据他们是否处于劣势来调整两个评级。如果雄鹿击败了网队，那么他们的排名会比网队击败雄鹿时增加的少，这纯粹是因为雄鹿有望获胜。这使得球队在整个赛季中获得动力。这方面的代码如下所示:
for i = 1:length(home_schedule)Ha = 1/(1+10^((team_data(home_schedule(i),1) - team_data(away_schedule(i),1) + 2.5) * -1/400)); % approx home team chance of winningAa = 1/(1+10^((-1 * team_data(home_schedule(i),1) + team_data(away_schedule(i),1)) * -1/400)); % approx away team chance of winningT = Ha + Aa - 1; % chance of a tieH = Ha - T/2; % adjusted home team winning chancesA = Aa - T/2; % adjusted away team winning chancesteam_data(home_schedule(i),1) = team_data(home_schedule(i),1) + 20 * (HW - H); % adjusts home team's ranking based on resultteam_data(away_schedule(i),1) = team_data(away_schedule(i),1) + 20 * (AW - A); % adjusts away team's ranking based on resultend
这种模式适用于所有 1230 场比赛，以决定季后赛的赛程。NBA 季后赛由来自每个联盟的 10 支球队组成。每个小组的前 6 名将自动进入第一轮，而 7-10 号种子将在一个小范围内决出第 7 和第 8 号种子。一旦两个联盟的所有 8 支球队都确定下来，季后赛就像一个标准的括号一样开始，两个联盟总决赛的获胜者将在冠军赛上相遇。每场比赛都是 7 场系列赛中的一场，最高分的种子将获得第 1、2、5 和 7 场比赛的主场比赛(如果必要的话)。
常规赛记录、季后赛记录和季后赛机会的所有数据都被记录下来，因为赛季重复了 50，000 次。
结果是:
下面显示了所有 30 支球队最有可能的常规赛记录，以及 50，000 个模拟赛季后的统计数据。最小和最大列代表该队在一个赛季中最少和最多的胜利。10%、50%和 90%是所有 50，000 次模拟的百分位数。
图片由作者提供
迈阿密热火队预计将取得 43 胜 39 负的战绩，有 10%的几率他们会赢不到 33 场比赛，有 10%的几率他们会赢超过 55 场比赛。这使得他们有 80%的机会赢得 33 到 55 场比赛。
每支球队的季后赛机会如下所示:
图片由作者提供
第一列代表那支球队进入前 10 名的机会，这意味着他们进入了季后赛。下一列是一个队进入第一轮比赛的机会(前 8)。
只有一支球队在所有 50，000 场比赛中没有赢得一个总冠军，那就是克利夫兰骑士队。奥兰多魔术队和俄克拉荷马雷霆队都只获得过一次总冠军。
NBA 总决赛夺冠热门前 3 名分别是密尔沃基雄鹿队(16.150%)、洛杉矶快船队(12.158%)、菲尼克斯太阳队(11.802%)。
我也对一支球队在 NBA 历史上拥有最好或最差赛季的几率感兴趣。至少有一支球队打破赛季胜利纪录的几率为 6.44%，其中密尔沃基雄鹿队占这些赛季的 34.70%。接下来两个可能的球队是洛杉矶快船(14.20%)和菲尼克斯太阳(13.77%)。
至少有一支球队打破赛季亏损纪录的几率为 12.08%，其中克利夫兰骑士队占这些赛季的 40.31%。接下来两个可能的球队是奥兰多魔术队(19.20%)和俄克拉荷马雷霆队(16.55%)。
在所有 50，000 个模拟赛季中，只有一支球队有机会打破赛季输赢记录。金州勇士队有 0.026%的几率输掉 NBA 历史上最多的比赛，但有 0.031%的几率赢得 NBA 历史上最多的比赛。
结论:
这个 NBA 赛季模拟显示，密尔沃基雄鹿队显然是今年 NBA 总决赛的热门球队(16.15%)，即使西部联盟占据了大多数冠军(54.05%)。
一支球队经历 NBA 历史上最糟糕赛季的几率几乎是一支球队经历 NBA 历史上最好赛季的两倍。
在所有 50，000 次模拟中，只有一支球队赢得了零冠军，那就是克利夫兰骑士队，预计战绩为 17 胜 65 负。
如果你对这类内容感兴趣，我推荐你阅读我的另一篇关于职业体育模拟的文章，NFL 模拟。
用 Python 模拟交通流
实践教程
实现微观交通模型
约翰·马特丘克在 Unsplash 上的照片
虽然交通并不总是顺畅，但汽车无缝穿越十字路口，在交通信号灯前转弯和停车，看起来会非常壮观。这种沉思让我想到了交通流量对人类文明的重要性。
在这之后，我内心的书呆子忍不住想了一个模拟交通流量的方法。我花了几周时间做一个关于交通流量的本科项目。我深入研究了不同的模拟技术，最终选择了一种。
在这篇文章中，我将解释为什么流量模拟是重要的，比较不同的方法可能建模流量，并提出我的模拟(连同源代码)。
为什么要模拟交通流？
模拟交通的主要原因是在没有真实世界的情况下生成数据。您可以使用软件上运行的模型来预测交通流量，而不是测试如何在现实世界中管理交通系统或使用传感器收集数据的新想法。这有助于加速交通系统的优化和数据收集。模拟是真实世界测试的一种更便宜、更快速的替代方式。
训练机器学习模型需要庞大的数据集，收集和处理这些数据集可能很困难，成本也很高。通过模拟交通按程序生成数据可以很容易地适应所需数据的确切类型。
建模
为了分析和优化交通系统，我们首先要对交通系统进行数学建模。这种模型应根据输入参数(路网几何形状、每分钟车辆数、车速等)真实地表现交通流量。
微观模型(左)宏观模型(右)。图片作者。
交通系统模型通常分为三类，这取决于它们运行的级别:
- **微观模型:**分别代表每一辆车，试图复制驾驶员的行为。
- **宏观模型:**用交通密度(每公里车辆数)和交通流量(每分钟车辆数)来描述车辆整体的运动。它们通常类似于流体流动。
- **介观模型:**是结合了微观和宏观模型特点的混合模型；他们将流量建模为车辆的“数据包”。
在本文中，我将使用一个微观模型。
微观模型
微观驾驶员模型描述了单个驾驶员/车辆的行为。因此，它必须是一个多智能体系统，也就是说，每辆车都使用来自其环境的输入来独立运行。
图片作者。
在微观模型中，每辆车都有编号 i 。第 i 辆跟随第( i-1 )辆。对于第 i 辆车，我们将用 xᵢ 表示其在道路上的位置， vᵢ 表示其速度， lᵢ 表示其长度。每辆车都是如此。
我们将用 sᵢ 表示保险杠到保险杠的距离，用δvᵢ表示第 i 辆车与其前车(车号 i-1 )之间的速度差。
智能驱动模型(IDM)
2000 年，Treiber、Hennecke 和 Helbing 开发了一种被称为智能驾驶模型的模型。它将第 i 辆车的加速度描述为其变量和其前方车辆变量的函数。动力学方程定义为:
在我解释这个模型背后的直觉之前，我应该解释一些符号代表什么。
我们已经讲过 sᵢ 、 vᵢ、和t31】δvᵢ**。其他参数是:**
- s₀ᵢ :车辆 i 和 i-1 之间的最小期望距离。
- v₀ᵢ : 是车辆的最大期望速度 i.
- δ : 是加速指数，它控制加速的“平滑度”。
- Tᵢ : 是第辆车驾驶员的反应时间。
- aᵢ : 为车辆的最大加速度一、
- bᵢ : 是车辆 i. 的舒适减速度
- s* :车辆 i 和 i-1 之间的实际所需距离。
首先，我们将查看 s* ，它是一个距离，由三项组成。
图片作者。
- s₀ᵢ : 前面说过，是最小期望距离。
- vᵢTᵢ : 是反应时间安全距离。它是驾驶员做出反应(刹车)之前车辆行驶的距离。
既然速度是距离随时间的变化，那么距离就是速度乘以时间。
- (vᵢδvᵢ)/√(2aᵢbᵢ):这是一个有点复杂的术语。这是基于速度差的安全距离。它代表车辆减速(不撞上前方车辆)，不过度制动(减速度应小于 bᵢ )所需的距离。
智能驱动模型如何工作
假设车辆沿直线行驶，并遵循以下等式:
为了更好地理解这个方程，我们可以把它分成两部分。我们有一个自由道路加速度和一个交互加速度。
自由道路加速度是自由道路上的加速度，也就是前方没有车辆的空旷道路。如果我们将加速度绘制成速度 vᵢ 的函数，我们得到 :
加速度是速度的函数。图片作者。
我们注意到，当车辆静止时( vᵢ=0 )，加速度最大。当车速接近最大速度 v₀ᵢ 时，加速度变为 0。这表明自由道路加速将使车辆加速至最大速度。
如果我们为不同的 **δ、**值绘制 v-a 图，我们会注意到，它控制着驾驶员在接近最大速度时减速的速度。这又控制加速/减速的平滑度
加速度是速度的函数。图片作者。
互动加速度与前方车辆互动相关。为了更好地理解它是如何工作的，让我们考虑以下情况:
- 在自由路面上(sᵢ > > s):* 当前方车辆很远，也就是距离 sᵢ 是主导期望距离 s* 时，交互加速度几乎为 0。
这意味着车辆将受到自由道路加速度的控制。
- 在高接近速率(δvᵢ): )下，当速度差很大时，交互加速度试图通过使用分子中的**(vᵢδvᵢ)**项制动或减速来补偿，但太难了。这是通过分母 4bᵢsᵢ 实现的。(老实说，我不知道它如何将减速精确地限制在 bᵢ 。
- 在小距离差时(sᵢ < < 1 和δvᵢ≈0): )加速度变成简单的排斥力。
交通道路网络模型
有向图的例子。图表(左)集合(右)
我们需要建立一个道路网的模型。为此，我们将使用一个有向图 G=(V，E) 。其中:
- V 是顶点(或节点)的集合。
- E 是代表道路的边的集合。
每辆车都有一条由多条道路(边)组成的路径。我们将对同一条道路(同一条边)上的车辆应用智能驾驶员模型。当一辆车到达路的尽头时，我们把它从那条路上移走，并把它附加到下一条路上。
在模拟中，我们不会保留一组节点(数组),相反，每条道路都将由其起始和结束节点的值明确定义。
随机车辆生成器
为了将车辆添加到我们的模拟中，我们有两个选项:
- 通过创建一个新的
Vehicle类实例并将其添加到车辆列表中，将每辆车手动添加到模拟中。 - 根据预定义的概率随机添加车辆。
对于第二个选项，我们必须定义一个随机车辆生成器。
随机车辆生成器由两个约束条件定义:
- 车辆生成率(τ): (每分钟车辆数)描述平均每分钟应添加多少车辆到模拟中。
- 车辆配置列表(L): 包含车辆配置和概率的元组列表。L = [(p₁，V₁)，(p₂，V₂)，(p₃，V₃)，…]
随机车辆生成器以概率 pᵢ 生成车辆 Vᵢ 。
交通灯
图片作者。
交通灯位于顶点，由两个区域表征:
- **减速区:**以一个减速距离 和一个减速系数为特征，是一个车辆使用减速系数降低其最大速度的区域。
- **停车区:**以停车距离为特征，是车辆停车的区域。这是通过以下动力学方程使用阻尼力实现的:
模拟
我们将采用面向对象的方法。每辆车和每条路都将被定义为一个类别。
我们将在许多即将到来的课程中重复使用下面的__init__函数。它通过一个函数set_default_config设置当前类的默认配置。需要一个字典，并将字典中的每个属性设置为当前类实例的属性。这样我们就不用担心更新不同类的__init__功能，也不用担心以后的变化。
路
我们将创建一个Road类:
在屏幕上绘制道路时，我们需要道路的length及其角度的余弦和正弦值。
模拟
还有一个Simulation班。我添加了一些方法来给模拟添加道路。
我们必须在屏幕上实时显示我们的模拟。为此，我们将使用pygame。我将创建一个期望一个Simulation类作为参数的Window类。
我定义了多个绘图函数来帮助绘制基本形状。
loop方法创建一个pygame窗口，并在每一帧调用draw方法和loop参数。当我们的模拟需要每帧更新时，这将变得有用。
我将名为trafficSimulator的文件夹中的每个文件与一个导入所有类名的__init__.py文件组合在一起。
每当一个新的类被定义时，它应该被导入到这个文件中
将trafficSimulator文件夹放在我们的项目文件夹中，我们就可以使用这个模块了。
模拟测试。图片作者。
车辆
现在，我们必须增加车辆。
我们将使用泰勒级数来近似本文建模部分讨论的动力学方程的解。
无限微分函数f的泰勒级数展开式为:
用▲x代替a，用x+▲x代替x，我们得到:
用位置x替换f:
作为近似，我们将在位置的 2 阶停止，因为加速度是最高阶导数。我们得到方程 (2) :
对于速度，我们将用v代替x:
我们将停在 1 阶，因为我们的最高阶导数是加速度(1 阶是速度)。方程式 (2) :
在每次迭代(或帧)中，使用 IDM 公式计算加速度后，我们将使用以下两个公式更新位置和速度:
在代码中是这样的:
因为这只是一个近似值，所以速度有时会变成负值(但是模型不允许这样)。当速度为负时，不稳定性出现，位置和速度发散到负无穷大。
为了克服这个问题，每当我们预测一个负速度时，我们将它设置为零，并从那里计算出:
在代码中，这是按如下方式实现的:
为了计算 IDM 加速度，我们将前导车辆表示为lead，并计算当lead不是None时的交互项(表示为alpha)。
如果车辆停止(例如在红绿灯处)，我们将使用阻尼方程:
然后，我们在一个Vehicle类中的一个update方法中将所有东西组合在一起:
在Road类中，我们将添加一个deque(双端队列)来跟踪车辆。队列是存储车辆的更好的数据结构，因为队列中的第一辆车是道路上最远的一辆车，它是可以从队列中移除的第一辆车。要从deque中移除第一个项目，我们可以使用self.vehicles.popleft()。
我们将在Road类中添加一个update方法:
和Simulation类中的update方法:
回到Window类，我添加了一个run方法来实时更新模拟:
现在，我们将手动添加车辆:
车辆在移动！图片作者。
车辆发电机
一个VehicleGenerator有一个包含(odds, vehicle)的元组数组。
元组的第一个元素是在同一个元组中生成车辆的权重(不是概率)。我使用权重是因为它们更容易处理，因为我们可以只使用整数。
例如，如果我们有 3 辆车，重量分别为1、3、2。这与1/6、3/6、2/6和6=1+3+2.相对应
为了实现这一点，我们使用以下算法
- 生成一个介于 1 和所有权重之和之间的数字。
- 当 r 为非负时:
遍历所有可能的车辆，并在每次迭代中减去其权重。 - 归还最后使用的车辆。
如果我们有权重: W₁ 、 W₂ 、 W₃ 。该算法将第一辆车分配到 1 和 W₁ 之间的号码，第二辆车分配到 W₁ 和 W₁+W₂ 之间的号码，第三辆车分配到 W₁+W₂+W₃ 之间的号码。
至于何时添加车辆，每当生成器添加车辆时，名为last_added_time的属性被更新为当前时间。当当前时间与last_added_time之间的持续时间大于车辆生成周期时，增加一辆车辆。
添加车辆的周期是60/vehicle_rate，因为vehicle_rate在每分钟车辆，60是 1 分钟或 60 秒。
我们还必须检查道路是否还有空间来容纳即将到来的车辆。我们通过检查道路上最后一辆车之间的距离以及即将到来的车辆的长度和安全距离之和来做到这一点。
最后，我们应该从Simulation的update方法中调用update方法来生成update车辆。
车辆正在产卵！图片作者。
红绿灯
交通信号的默认属性是:
self.cycle是一个元组数组，包含self.roads中设置的每条道路的状态(True表示绿色，False表示红色)。
在默认配置中，(False, True)表示第一组道路是红色，第二组是绿色。(True, False)则相反。
使用这种方法是因为它易于扩展。我们创造的交通灯包括两条以上的道路，交通灯有单独的右转和左转信号，甚至有多个交叉路口的同步交通信号。
交通信号的update功能应该是可定制的。其默认行为是对称的固定时间循环。
我们需要将这些方法添加到Road类中:
而这个，在Road的update方法中。
并用Simulation的update方法检查交通灯状态:
停下来。图片作者。
曲线
在现实世界中，道路是有曲线的。虽然从技术上讲，我们可以通过手写大量道路的坐标来近似一条曲线，从而在这个模拟中创建曲线，但我们也可以在程序上做同样的事情。
为此，我将使用贝塞尔曲线。
我创建了一个curve.py文件，其中包含了帮助创建曲线并通过道路索引引用它们的函数。
测试:
优美的曲线。图片作者。
例子
这些示例的代码可以在本文底部链接的 Github 资源库中找到。
高速公路上
高速公路上。图片作者。
双向交叉路口
双向交叉路口。图片作者。
迂回的
迂回。图片作者。
分叉菱形立交
分叉菱形立交。图片作者。
限制
虽然我们可以修改Simulation类来存储关于我们稍后可以使用的模拟的数据，但是如果数据收集过程更加精简的话会更好。
这个模拟还欠缺很多。弯道的实施是糟糕和低效的，并导致车辆和交通信号之间的相互作用的问题。
虽然有些人可能认为智能司机模型有点矫枉过正，但拥有一个能够复制真实世界现象的模型很重要，如交通波(又名幽灵交通蛇)和司机反应时间的影响。出于这个原因，我选择使用智能驱动程序模型。但是对于准确性和极端真实性不重要的模拟，如视频游戏，IDM 可以由更简单的基于逻辑的模型代替。
完全依赖基于模拟的数据会增加过度拟合的风险。您的 ML 模型可以针对只存在于模拟中而不存在于真实世界中的零食进行优化。
结论
仿真是数据科学和机器学习的重要组成部分。有时，从现实世界中收集数据要么是不可能的，要么是昂贵的。生成数据有助于以更好的价格建立庞大的数据集。模拟也有助于填补真实世界数据的空白。在某些情况下，真实世界的数据集缺乏可能对开发的模型至关重要的边缘案例。
这个模拟是我参与的一个大学项目的一部分。目的是优化城市十字路口的交通信号。我做这个模拟是为了测试和验证我的优化方法。
我从来没有想过要发表这篇文章，直到我在看特斯拉的 AI day，在这篇文章中，他们谈到了他们如何使用模拟来为边缘案例生成数据。
源代码和贡献
这里有一个到 Github 库 的链接，包含本文中的所有代码，包括例子。
如果您对代码有任何疑问或问题，请随时联系我，或者在 GitHub 上提交请求或问题。
数据科学模拟
有三个代码示例
在本文中，我们将看三个模拟示例。
- 统计理论 101
- 蒙蒂·霍尔问题
- 航线优化
介绍
在高中和研究生院的时候，我总是纠结于与以下相关的课程:
- 证明
- 统计理论
通过实践，我最终理解了这些课程中的信息。现在，我作为一名实践数据科学家进入了现实世界，显然这些课程会为我提供日常工作和生活所需的丰富知识…我不会走那么远！
说实话，我确实在这些课程中学到了很多。但实际上，我有了一种新的方法来回答这些课上提出的很多问题。其实这几天几乎每个问题我都用同样的方法回答！解决方法是 模拟 。
让我们看一个统计理论导论课上的例子。
统计理论 101
假设你的任务是找出将一枚硬币抛 50 次，恰好得到 20 条反面的概率？
你知道这是一个二项分布吗？如果你这样做了，很可能你很快就计算出这些几率是 4.2%。我个人在statisticshelper.com用在线计算器找到了答案，上面显示的是概率密度函数(以下所有公式/解都是从上述网站截取的图片):
这张图片显示了我们的问题设置中每个字母代表的含义。
如果硬币被操纵了，它 80%的时候都是正面朝上。现在，20 条尾巴的几率是多少？如果你对上面的公式感到满意，你会发现答案是 0.06%。
但是如果你不知道这是二项分布呢？或者如果那个公式对你没有意义呢？或者如果问题不同呢？
如果问题是第一条尾巴在第五次翻转时翻转并且前四条都是正面的几率有多大呢？第一条尾巴出现在第 100 次翻转时，前 99 条是正面的可能性有多大？
现在，你需要知道几何分布。🤦♂️
你在第 100 次空翻中完成第 5 次空翻的几率有多大？而另外 4 次空翻可以分散在前 99 次空翻中的任何地方。几率有多大？
现在你需要了解负二项分布。️️🤦♀️.这已经失控了！
模拟示例 1:二项式分布
如果你不理解上面的任何解决方案，没关系。这就是这个帖子的力量。你不需要！
除了了解所有这些统计分布，您还可以模拟结果。让我们重温一下第一个问题:假设你的任务是找出将一枚硬币抛 50 次，恰好得到 20 条反面的概率？
我们可以用上面的方法模拟一次翻转。random.random()以相等的概率返回 0 和 1 之间的值。所以默认情况下，coin_flip()模拟公平的硬币投掷，但是通过更新默认的probability值也允许不公平的硬币。
现在我们可以模拟 50 次抛硬币，看看结果。
作者图片
您的结果可能会有所不同，因为这是随机的，我们没有设置种子。你可以看到在第一组 50 次抛硬币中，我们得到了 26 个正面和 24 个反面。问题是我们多久才能得到 20 条尾巴。我们现在可以模拟这个过程 10，000 次，看看 20 条尾巴出现的频率。
运行这段代码，我得到了大约 4.2%的结果！🎉 🎊
您可以通过将probability更改为 0.2 来轻松进行修改。我的模拟结果提供了 0.0005，这与理论结果并不完全相同。然而，在现实世界中，这是一个非常合理的值。
作者图片
模拟示例 2:蒙蒂霍尔问题
如果你对三门蒙蒂霍尔问题不熟悉，大卫伯顿的这篇文章可以帮助你加快速度。你也可以在 UCSD 网站这里自己试题。总的来说，有三扇门。一扇门后面有辆车，另外两扇门是哑弹。你选择了一扇门，然后给你看了一个哑弹。此时，您可以将您的选择更新到另一个门，或者继续使用您最初选择的门。
对我来说，在看到一扇坏门后换门以增加选择汽车的几率从来都不是直觉。然而，我们可以模拟这个问题，以证明更新我们的选择是正确的决定。
下面的代码将迭代一个游戏，其中我们选择用门的变化来更新我们最初的选择。
我们现在可以运行这个场景的 10，000 次迭代，并跟踪我们赢得汽车的频率。这可能会有一些变化，但我的模拟运行提供了几乎准确的理论解决方案，其中 66%的情况下我们会选择更换门，而只有 33%的情况下我们会选择保留原来的门。
模拟示例 3:航线问题
作为最后一个模拟示例，假设您为一家航空公司工作。想象一下，我们知道一架飞机有 100 个座位，每个买了票的乘客有 10%的几率不会出现。另外，每张票的利润是 75 美元。如果你超售了飞机，你会为每个因没有座位而需要转到新航班的乘客赔钱。这些超售席位的损失是 25 美元。
想象一下，无论乘客是否出现，你都会获得利润。
如果你卖出 100 张票，你会赚 7500 美元。
如果你卖出 101 张票，每个人都来了，你就赚了
75 美元 X100-25 美元 x1 = 7，475 美元
但是，如果至少有 1 人缺席，您将获得:
75 美元 x101 = 7 575 美元
为了实现利润最大化，应该出售的最佳座位数是多少？
我们可以使用下面的函数计算收入，默认情况下销售 100 张门票，这应该总是提供 7，500 美元的收入。
然后，我们可以模拟多次运行calculated_revenue函数，但是销售不同数量的门票，以了解最佳数量。假设我们在每一个可能售出的门票数量上运行 10，000 次迭代。然后我们看看每张售出的票的可能收入。下面的代码执行模拟。
然后我们可以画出可能的收入值。
作者图片
如果我们想要最大可能的收入，我们可能会选择增加风险，直到卖出 130 张以上的门票。然而，潜在的不利因素导致售出的座位数非常高时的预期利润低于我们在座位数低得多时的预期利润。取决于我们的风险承受能力，没有正确的答案。该图提供了做出明智决策所需的信息。
使用平均值(在大多数传统课程中建议的最优值)提供了以下图表和售出 115 个座位的解决方案。
作者图片
结论
我们的模拟解决方案之旅到此结束。下次你遇到问题时，如果你不确定正确答案，记得你可以模拟各种可能性！然后利用结果做出正确的决定！
基于数据块的单节点分布式深度学习
理解大数据
卢克·范·德·维尔登& 里克·琼格里乌斯
Databricks 是一个分析生态系统，现在可以在大多数主要的云提供商 Google、AWS 和 Azure 上使用。数据块集群计算使用分布式 Spark 引擎。最近推出的单节点 Spark 集群不支持分布式计算，为什么？
单节点数据块集群
在多节点集群上，带有 PySpark 的 Python 解释器在驱动节点上运行以收集结果，而工作节点执行 JVM jar 文件或Python UDF。运行单节点集群的选项于 2020 年 10 月日推出，其动机如下。
标准数据块 Spark 集群由一个驱动节点和一个或多个工作节点组成。这些集群最少需要两个节点——一个驱动程序和一个工作程序——来运行 Spark SQL 查询、从增量表中读取数据或执行其他 Spark 操作。然而，对于许多机器学习模型训练或轻量级数据工作负载，多节点集群是不必要的。
单节点集群可以运行 Databricks 机器学习运行时，该运行时经过精心策划，开箱即用。您可以选择 GPU 或 CPU 运行时，并将其部署在各种规模的虚拟机上。非常适合训练我们的机器学习模型。
熊猫数据帧
在单节点集群上，Spark 命令以本地模式运行。单节点集群仍然可以使用 Spark 数据源和相关的云认证/授权支持，比如 Azure 上的 ADLS 认证直通。除了数据提取和加载，Spark 功能在单节点集群上并不那么有用。为了训练机器学习模型，我们经常从 Spark 转换到 Pandas 数据帧。自从 Spark 2.3-2.4 以来，熊猫之间的转换已经通过箭头加速。有时，我们甚至在数据砖块上立即将拼花地板(数据湖)转换成熊猫，见下文。
# azure data lake parquetdf = spark.read.parquet("abfss://<storage>/path").toPandas()
Azure 用例
让我们讨论一下在计算管道中使用单节点集群。我们希望训练多个 DeepAR 模型，每个模型在我们的数据集中有不同的时间序列。这是令人尴尬的并行，因为每个模型都可以在内存中有一个时间序列的单个 CPU 虚拟机上训练。每个培训工作相当于运行一个 Databricks 笔记本，它导入一个特定的时间序列和一个具有单元测试功能的 Python 包。我们有数百个时间序列，并希望旋转数十个单节点来执行我们的训练工作和注册我们的模型。我们使用 Azure Data Factory 来安排我们的 Databricks 笔记本和启动集群。
DataFactory-Databricks 架构，图片由作者提供
Azure 数据工厂的并行性
DataFactory 管道可以并行运行 Databricks 笔记本，并等待它们完成，然后继续进行管道的下一个活动。ForEach 操作符为序列中的元素(例如数据湖拼花路径)启动一个记事本。它可以按批量顺序或并行运行。请参见下面数据工厂 JSON 导出中的“isSequential”和“batchCount”。我们启动了 20 个运行相同数据块的单节点笔记本活动来处理我们的时间序列。数据工厂通过数据块链接服务管理这些集群资源。
{
"name": "train_each_time_series",
"type": "ForEach",
"dependsOn": [
{
"activity": "generate_time_series_paths",
"dependencyConditions": [
"Succeeded"
]
}
],
"typeProperties": {
"items": {
"value": "[@json](http://twitter.com/json)('[{\"time_series\": \"some/path\"}]')",
"type": "Expression"
},
"isSequential": false,
"batchCount": 20,
"activities": [
{
"name": "train_a_time_series",
"description": "Train a time series DeepAR",
"type": "DatabricksNotebook",
"dependsOn": [],
"userProperties": [],
"typeProperties": {
"notebookPath": "/our_model/train_model"
},
"linkedServiceName": {
"referenceName": "[parameters('db_pool')]",
"type": "LinkedServiceReference"
}
}
]
}
}
数据块实例池
要从 DataFactory 运行 Databricks 笔记本，我们需要在 DataFactory 中创建一个链接服务。启动笔记本活动时，有三种方法可以配置链接服务来选择或创建 Databricks 集群。
- 互动集群是一个预先存在的运行集群，我们可以选择我们的笔记本活动。
- 为启动的每个笔记本活动创建作业群。作业集群在笔记本活动期间存在。这会导致集群创建和终止的等待时间很长(两者都需要几十分钟)。
- 数据块实例池包含可配置数量的就绪虚拟机实例，这些实例会等待被包含到集群中。空闲和运行之间的延迟是几分钟。从 2019 年末开始支持，实例池允许作业之间的快速周转时间。
待续
关于我们使用 Python SDK 的 Databricks 笔记本和数据工厂管道的 CICD 工作流的后续博客将进一步深入这个用例。在下一节中，我们想知道如果我们不使用 Spark 进行模型训练，多节点数据块集群有什么用。
分布式深度学习
当我们在具有足够内存和 GPU/CPU 资源的单节点上运行 Databricks 笔记本时，我们已经看到了单节点 Databricks 集群对于机器学习的价值。让我们回溯并考虑分布式计算。多节点数据块集群只对 Spark 计算有用吗？我们能否在多节点数据块集群之上搭载一个不同的分布式计算框架？
Databricks Horovod 转轮
Horovod ( 开源优步)是一个使用 MPI 和 NCCL 的分布式深度学习框架，支持 TensorFlow、Keras、PyTorch 和 Apache MXNet。Databricks 的 Spark-Deep-Learning 通过机器学习运行时支持 data bricks 集群上的 Horovod。它提供了一个 HorovodRunner 来对 Spark 任务中的多个工人运行 Python 深度学习。Spark 任务在 worker 节点(spark executor)上的现有 SparkContext 中运行，并涉及 worker 节点上的 Python 运行时，它与共存的 JVM 通信。
使用分布式数据
通常，分布式 Horovod 优化器使用梯度下降计算，通过 MPI 在节点之间进行通信。数据提取和预处理在 worker 节点上的 Python 解释器中完成。坚持培训成果；DBFS 上的一个检查点目录被一个工作人员用来存储模型检查点，参见下面的简明代码。.)
import horovod.torch as hvdfrom sparkdl import HorovodRunnerlog_dir = "/dbfs/ml/horovod_pytorch"def train_hvd(learning_rate): hvd.init() train_dataset = get_data_for_worker(rank=hvd.rank())
train_loader = torch.utils.data.DataLoader(
train_dataset,
batch_size=batch_size,
sampler=train_sampler,
) optimizer = optim.SGD(
model.parameters(),
lr=learning_rate,
momentum=momentum,
) # hvd.DistributedOptimizer handles the distributed optimization
optimizer = hvd.DistributedOptimizer(
optimizer, named_parameters=model.named_parameters()
) # all workers start with the same initial condition
hvd.broadcast_parameters(
model.state_dict(), root_rank=0
) for epoch in range(1, num_epochs + 1):
train_epoch(
model, device, train_loader, optimizer, epoch
) # save model checkpoints only from one worker
if hvd.rank() == 0:
save_checkpoint(
log_dir, model, optimizer, epoch
)# initialize a runner with two workers
hr = HorovodRunner(np=2)
# launch our training function across the workers
hr.run(train_hvd, learning_rate=0.001)
在为 Horovod runner 编写 Python 代码时，我们不能使用 Spark SQL 功能从类似 HDFS 的文件系统中查询数据。只有 Spark driver 节点可以协调该功能并收集其结果。但是，也有从类似 HDFS 的数据湖存储中读取的选项。
佩塔斯托姆
Petastorm(也是优步的开源)符合我们的需求。该库允许在深度学习模型的单节点和分布式训练期间使用数据湖存储的拼花文件。它还可以使用内存中的 Spark 数据帧。Petastorm 将物化的 Spark 数据帧转换为 PyTorch 和 Tensorflow 的不同数据加载器。预处理转换可以作为 TransformSpec 添加到 Petastorm 转换器中，它作用于 Pandas 数据帧。
我们可以使用 Petastorm 进行分布式深度学习，如果我们将 Petastorm 转换器中的’ shard_count '配置为 Horovod runners 的数量，并让 Horovod 根据工人的等级选择合适的 shard(py torch 示例)。我们还会根据 Horovod 工作线程的数量小心地划分底层的 parquet 文件或内存中的数据帧。这太棒了！HDFS 式文件系统上的分布式数据可用于深度学习模型的分布式训练。
通过流媒体实现未来协同效应
我们已经讨论了 Horovod(深度学习)与 Databricks spark cluster 和 Petastorm 的集成，用于数据加载。如何改进 Spark 和 Horovod 集成？我们希望每个 worker 节点上的 worker JVM 和 Python 解释器之间的数据传输更加容易。Arrow 只是一项技术，并且已经得到支持。一个解决方案方向可能是火花结构化流。
结构化流抽象是基于一个连续的数据框架，随着事件从流源进入而增长(f.i. Kafka)。操作以微批处理或连续方式执行，有三种输出格式(追加、更新、完成)。PySpark API 落后于 Java/Scala API ，不支持有状态或无状态的任意操作。性能限制可能是不支持任意 Python 操作的原因，也是 API 关注 Spark SQL 聚合的原因。
愿望列表:任意 PySpark 流操作
在训练期间，我们的 DL 模型表示我们希望在每个训练步骤中更新的状态。传递给 HorovodRunner 的 Python UDF 是对一批数据的有状态操作。例如，在数据集中的每个图像上训练 DNN。我们可以通过 Spark 结构化流将 HDFS 或卡夫卡的每幅图像传输到 Python 进程，正好赶上培训，就像在在线学习一样。下面的非功能性代码让我们对它的外观有了一些了解(我们忽略了事件窗口逻辑)。我们将训练数据流与模型权重流相结合，并对模型进行增量训练，将更新后的模型权重添加到模型权重流中，用于下一个处理步骤。
# Pseudo code not functional
schema = ...
train_data = spark.readStream.format("kafka").option(
"subscribe", "train_data"
)model_weights = spark.readStream.format("kafka").option(
"subscribe", "model_weights"
)train_data.join(
model_weights,
expr("train_data.model_id == model_weights.model_id"),
).filter(col("train_data.rank") == hvd.rank).groupByKey(
model_id
).mapGroupsWithState(
train_model
).writeStream.format(
"kafka"
).option(
"topic", "model_weights"
).outputMode(
"append"
).start()
包裹
这个博客围绕着数据块集群的非 Spark 中心的使用。我们从在单节点集群上并行运行深度学习训练转移到在多节点集群上运行分布式训练。我们可以使用分布式数据进行模型训练。此外，我们还提出了缩小 Spark/Databricks 集群的数据密集型和计算密集型用途之间差距的方法。在后续博客中，我们将更详细地讨论我们如何通过 Databricks 笔记本和数据工厂管道来接近 CICD。
最初发布于https://codebeez . nl。
奇异值分解
Python 中的解释、推导和应用
来自 Pexels 的 Nicole Avagliano 的照片
介绍
对数据科学、机器学习和人工智能至关重要的线性代数经常被忽视，因为大多数入门课程未能展示整体情况。从从业者的角度来看，像特征分解和奇异值分解(SVD)这样的概念非常重要；它们是包括主成分分析(PCA)和潜在语义分析(LSA)在内的降维技术的核心。本文旨在展示 SVD，方法是结合有形的 Python 代码，温和地介绍所需的数学知识。
奇异值分解
矩阵乘法
首先，让我们考虑下面的向量， x ，作为两个基向量 i 和 j 的和。
由作者生成的图像
我们可以很容易地使用 matplotlib 和 Python 来可视化这个向量…
由作者生成的图像
矢量相对简单，它们有方向和大小。在这个例子中，我们从原点画出向量 x 。然而，这个相同的向量可以在ℝ的任何地方等价地绘制。
现在让我们考虑一个矩阵和矩阵乘法。数学上…
由作者生成的图像
使用标准矩阵乘法…
由作者生成的图像
在 Python 中…
此外，我们可以画出矩阵和向量的乘积…
由作者生成的图像
我们可以看到矩阵乘法( Ax )改变了我们原始向量 x 的方向和大小。
E 本质上，矩阵旋转并且拉伸我们的原始向量。
这可以通过在同一张图上绘制两个向量来更有效地看出…
由作者生成的图像
更具体地说，我们可以看到原始向量 x 被旋转了一个角度 θ ，并被拉伸了一个标量 α 。分解成旋转和拉伸的概念可以推广。
然而，在我们继续之前，代替我们的原始向量， x ，现在让我们考虑向量的集合作为矩阵 V.
由作者生成的图像
在这种情况下，我们的向量集合包含原始基向量( i 和 j )，它们相加产生 x 。
维度一般可以扩展到mxn—使用2x2辅助可视化。
这种旋转和拉伸的概念也可以在 AV 中看到…
由作者生成的图像
旋转矩阵
这是负责旋转原矢量的矩阵 θ 。
由作者生成的图像
重要的是要注意旋转矩阵是酉变换。酉变换的一个特殊性质是它们的逆是(平凡的)复共轭转置。
旋转很容易可视化，尤其是在二维空间中。然而，该矩阵可以被修改以实现更高维度的旋转，同时指定特定的旋转轴。
拉伸矩阵
这是负责拉伸原向量由 α 的矩阵。
由作者生成的图像
向更高维度的扩展是非常明显的——通过简单地对角化相关维度空间中的拉伸参数。
简化奇异值分解
这个难题的所有部分现在对我们来说都是可用的。我们知道，乘积 AV ，包含一个旋转和一个拉伸。为了进一步推广这一观点，我们还将考虑…
- dim(A)=mxn
- dim(U _ hat)=mxn
- dim(σ_ hat)=nxn
- dim(V )=nxn*
数学上…
由作者生成的图像
如果我们将右边两边乘以 V 的倒数…
在这种情况下，我们将把 V 的倒数称为 V*
由作者生成的图像
由作者生成的图像
由作者生成的图像
这就是约化奇异值分解。
奇异值分解
我们可以通过向旋转矩阵添加 m-n 列，向拉伸矩阵添加 m-n 行来实现这一点。利用这个想法，我们将 U_hat 和σ_ hat重新定义为 U 和σ…
- dim(A)=mxn
- dim(U)=mxm
- dim(σ)=mxn
- dim(V )=nxn*
由作者生成的图像
这被正式称为奇异值分解。其中σ包含拉伸元素，奇异值，降序排列。这种分解的主要好处是它适用于任何矩形或正方形矩阵。
解析奇异值分解
如前所述，特征分解与奇异值分解密切相关。我们可以把上面已经得到的转化为初等线性代数中的经典特征值问题来寻找奇异值分解。
由作者生成的图像
由作者生成的图像
回想一下 U 是酉变换…
由作者生成的图像
由作者生成的图像
将两边乘以矢量集合V…
由作者生成的图像
由作者生成的图像
由作者生成的图像
回想一下特征值问题…
由作者生成的图像
我们可以在这个问题的上下文中匹配以下元素…
- 甲→ A^TA
- x → V
- λ→σ
通过从 AA^T.开始，可以采用相同的程序找到 U
Python 中的奇异值分解
Python 使得使用 numpy 找到矩阵的奇异值分解变得非常容易。
*array([[ 2., 3.],
[-2., 4.]])*
在上面的代码片段中，我们找到了矩阵 A 的奇异值分解，也展示了原始矩阵的 SVD 重构。
在其分解形式中，我们还可以通过线性变换原始基向量集合 V 来可视化奇异值分解的元素。
由作者生成的图像
由作者生成的图像
图像压缩中的应用
也许奇异值分解最直观的例子之一来自图像压缩。首先，我们将读入一幅图像，并找到奇异值分解。接下来，我们将把秩降低到包含奇异值的矩阵的任意三层(σ)。最后，我们将重建降秩后的图像。
在图像中阅读
注意:确保你的 cwd 有一个图像，这样你就可以跟随
照片由 Aaron Burden 在 Unsplash 上拍摄——情节由作者生成
移除颜色通道
在奇异值分解之前，我们要移除与这张图片相关的颜色通道…
Aaron Burden 在 Unsplash 上的照片——作者生成的情节
奇异值分解
现在我们可以计算代表这个图像的矩阵的奇异值分解。我们将像以前一样使用 numpy，并使用 matplotlib 可视化这个分解的组件…
Aaron Burden 在 Unsplash 上的照片——作者生成的情节
为了更好地理解奇异值，我们可以按降序绘制前 50 个值…
由作者生成的图像
降级σ
现在让我们计算对角矩阵 S 的秩(σ)…
*480*
我们可以选择保留哪些奇异值，并将其余的设置为零以降低秩——因此，我们可以通过取分解和绘图的乘积来查看对原始图像的影响…
10 个奇异值
Aaron Burden 在 Unsplash 上的照片——作者生成的情节
25 个奇异值
图片由 Aaron Burden 在 Unsplash 上拍摄——由作者生成的情节
50 个奇异值
Aaron Burden 在 Unsplash 上的照片——作者生成的情节
结论
本文介绍了奇异值分解及其应用。首先，我们对矩阵乘法有了一个基本的了解，并将运算分解为一个旋转和一个原始向量(或向量集合)的拉伸。然后我们通过解析地和计算地寻找奇异值分解来推广旋转和拉伸的概念；同时示出了分解对基向量集合的影响。最后，我们通过降低含有奇异值的矩阵的秩(σ)给出了奇异值分解在图像压缩中的应用。
扫一扫
745
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
