如何通过操作符融合和 CPU/GPU 代码生成加速深度学习
介绍了 PyTorch 中的深度学习编译器技术，包括图捕获、中间表示、操作符融合以及优化的 C++和 GPU 代码生成
·
关注 发表在 Towards Data Science · 17 分钟阅读 · 2023 年 4 月 20 日
–
插图由作者提供
计算机编程是神奇的。我们用人类可读的语言编写代码，然后仿佛魔法一样，它被转换成通过硅晶体管的电流，使其像开关一样工作，并实现复杂的逻辑——只是为了让我们能在互联网上欣赏猫咪视频。在编程语言和运行它的硬件处理器之间，有一个重要的技术组件——编译器。编译器的工作是将我们人类可读的语言代码翻译并简化为处理器能理解的指令。
编译器在深度学习中扮演着非常重要的角色，提升训练和推理性能，提高能源效率，并针对各种 AI 加速器硬件。在这篇博客文章中，我将讨论支持 PyTorch 2.0 的深度学习编译器技术。我将带你了解编译过程的不同阶段，并通过代码示例和可视化讨论各种基础技术。
什么是深度学习编译器？
深度学习编译器将用深度学习框架编写的高级代码转换为优化的低级硬件特定代码，以加速训练和推理。它通过执行层和操作符融合、更好的内存规划，生成目标特定的优化融合内核来减少函数调用开销，从而在深度学习模型中发现性能优化的机会。
插图由作者提供
与传统软件编译器不同，深度学习编译器必须处理高度并行化的代码，这些代码通常在专用的 AI 加速硬件（如 GPU、TPU、AWS Trainium/Inferentia、Intel Habana Gaudi 等）上加速运行。为了提高性能，深度学习编译器必须利用硬件特定的功能，如混合精度支持、性能优化的内核，并尽量减少主机（CPU）与 AI 加速器之间的通信。
虽然深度学习算法正在以迅猛的速度不断进步，但硬件 AI 加速器也在不断发展，以跟上深度学习算法的性能和效率需求。我在之前的博客文章中讨论了算法与 AI 加速器的共同演进：
在 AI 加速器（如 NVIDIA GPU、Intel Habana Gaudi 和 AWS 等）中，机器学习的高效算法和方法——
towardsdatascience.com
在本文中，我将专注于软件方面的内容，特别是接近硬件的软件子集 —— 深度学习编译器。首先，让我们来看看深度学习编译器中的不同功能。
PyTorch 2.0 中的深度学习编译器
PyTorch 2.0 包括新的编译器技术，以提高模型性能和运行时效率，并通过简单的 API 针对不同的硬件后端进行目标化：torch.compile()。虽然其他博客文章和文章详细讨论了 PyTorch 2.0 的性能优势，但在这里，我将专注于当您调用 PyTorch 2.0 编译器时内部发生的情况。如果您想要量化的性能优势，可以在来自 huggingface、timm 和 torchbench 的不同模型的性能仪表板找到。
在高层次上，PyTorch 2.0 深度学习编译器的默认选项执行以下关键任务：
-
图捕获：为您的模型和函数表示计算图。PyTorch 技术：TorchDynamo，Torch FX，FX IR
-
自动微分：使用自动微分进行后向图追踪并降低到基本运算符。PyTorch 技术：AOTAutograd，Aten IR
-
优化：前向和后向图级优化以及操作融合。PyTorch 技术：TorchInductor（默认）或其他编译器
-
代码生成：生成硬件特定的 C++/GPU 代码。PyTorch 技术：TorchInductor，OpenAI Triton（默认）和其他编译器
通过这些步骤，编译器将转换您的代码并生成逐渐“降低”的中间表示（IR）。“降低”是编译器词汇中的一个术语，它指的是通过自动转换和重新编写将广泛的操作（例如由 PyTorch API 支持的操作）映射到狭窄的操作集（例如由硬件支持的操作）。PyTorch 2.0 编译器流程：
如果您对编译器术语不熟悉，请不要被这一切吓到。我也不是编译器工程师。继续阅读，随着我将使用一个简单的示例和可视化来详细解释，一切将变得清晰起来。
通过 torch.compile() 编译器过程的详细步骤
注意：整个步骤详见此处的 Jupyter Notebook
为了简单起见，我将定义一个非常简单的函数并将其通过 PyTorch 2.0 编译器流程运行。您可以用深度神经网络模型或 nn.Module 子类替换此函数，但这个例子应该能帮助您更好地理解底层发生的事情，而不是复杂的多参数模型。
该函数的 PyTorch 代码：
def f(x):
return torch.sin(x)**2 + torch.cos(x)**2
如果你在高中三角学课程中认真听讲，你会知道我们的函数值对于所有实值 x 始终为 1。这意味着它的导数是一个常数的导数，必须等于零。这在验证函数及其导数的作用时会很有用。
现在，是时候调用 torch.compile()。首先，让我们说服自己编译这个函数不会改变其输出。对于相同的 1x1000 随机向量，函数输出与 1 的向量之间的均方误差应在编译和未编译的函数之间为零（在某些误差容限下）。
作者截图
我们所做的只是添加了一行额外的代码 torch.compile() 来调用我们的编译器。现在让我们看看每个阶段发生了什么。
图形捕获：PyTorch 模型或函数的计算图表示
PyTorch 技术： TorchDynamo, FX 图, FX IR
编译器的第一步是确定编译内容。进入 TorchDynamo。TorchDynamo 截取你的 Python 代码的执行，并将其转换为 FX 中间表示（IR），并将其存储在一个称为 FX 图的特殊数据结构中。你问这看起来是什么样的？很高兴你问了。下面，我们将查看生成这些的代码，但这是转换和输出的结果：
作者截图
重要的是要注意，Torch FX 图只是 IR 的容器，并不真正指定它应该包含哪些操作符。在下一部分，我们将看到 FX 图容器再次出现，但带有不同的 IR 集。如果你比较函数代码和 FX IR，它们之间几乎没有差别。实际上，它是你编写的相同 PyTorch 代码，只是以 FX 图数据结构期望的格式进行布局。当执行时，它们将提供相同的结果。
如果你在没有任何参数的情况下调用 torch.compile()，它将使用默认设置，这会运行整个编译器栈，包括默认的硬件后端编译器 TorchInductor。但如果我们现在讨论 TorchInductor 会有点超前，因此暂时搁置这个话题，等我们准备好时再回来。首先我们需要讨论图形捕获，我们可以通过截取 torch.compile() 的调用来实现这一点。以下是我们将要做的：torch.compile() 允许你提供自己的编译器，但由于我不是编译器工程师，也不知道如何编写编译器，我将提供一个虚拟编译器函数来捕获 TorchDynamo 生成的 FX 图 IR。
以下是我们为 torch.compile() 函数编写的虚拟编译器后端函数 inspect_backend，在这个函数中我做了两件事：
-
打印由 TorchDynamo 捕获的 FX IR 代码
-
保存 FX 图形可视化
上述代码片段的输出是 FX IR 代码和图表，显示了我们的函数 sin²(x)+cos²(x)
作者的屏幕截图
请注意，我们的虚拟编译器 inspect_backend 函数是在我们使用一些数据调用编译函数时才调用的，即当我们调用 compiled_model(x) 时。在上述代码片段中，我们只评估函数，或者用深度学习术语说，“进行前向传播”。在接下来的部分，我们将利用 PyTorch 的自动微分引擎 torch.autograd 来计算导数和“反向传播”图。
自动微分：前向和后向计算图
PyTorch 技术: AOTAutograd, 核心 Aten IR
TorchDynamo 给出了前向传播函数评估作为 FX 图，但是后向传播呢？为了完整起见，我将偏离我们的主要话题，说一点为什么我们需要评估函数相对于其权重的梯度。如果你已经了解数学优化如何工作，请跳过这节。
反向传播和反向图是什么？
深度学习和机器学习中的“学习”部分是一个数学优化问题，简单地说就是：找到一个变量 w 的值，使得某个关于 w 的函数取得最小值。或者更简洁地说：
在机器学习中，f(w) 是由权重参数化的损失函数。 f(w) 可以更清晰地表示为训练标签与模型根据训练数据对训练标签进行预测之间的错误度量：
结果发现，如果我们能够计算损失相对于权重的“减少速率”，我们就可以更新我们的权重，使其向一个越来越小的损失 f(w) 前进一步。换句话说，我们必须朝着更好地拟合训练数据集的模型前进。我们可以通过计算在给定 w 处损失 f(w) 的最陡斜率，来找到下一个权重值，并扰动 w 朝着那个方向前进。函数关于权重的斜率，就是其相对于权重的导数。由于存在多个权重值，导数变成一个矢量量，称为梯度，是一个针对每个权重的分量的偏导数的矢量。权重 w 在每次迭代时通过梯度的某个函数 g() 进行扰动，如下所示：
其中函数 g(.) 取决于优化器（例如 sgd, sgdm, rmsprop, adam 等）。
对于 SGD 权重更新步骤如下：
PyTorch 2.0 如何跟踪反向传播图？
首先，让我们计算我们期望反向传递图应该是什么样的，然后与 PyTorch 生成的结果进行比较。对于我们的简单函数，前向图和反向图应该实现以下函数。如果 sin 和 cos 让你感到困扰，你可以想象 f(x) 是应用于神经网络的损失函数。
PyTorch 使用反向模式的 自动微分 来计算梯度，PyTorch 的自动微分实现称为 Autograd。PyTorch 2.0 引入了 AOTAutograd，它在执行之前预先跟踪前向和反向图，并生成一个联合的前向和反向图。然后，它将前向图和反向图分成两个独立的图。前向图和反向图都存储在 FX 图数据结构中，可以如下面所示进行可视化。
作者截图
你可以通过遍历图上的节点来验证数学计算是否正确。AOTAutograd 生成的反向传递确实计算了我之前分享的方程中的导数，该导数应该等于零，因为原始函数仅生成了单位矩阵。
现在我们将通过扩展我们伪编译器函数 inspect_backend 来运行 AOTAutograd，调用 AOTAutograd 并生成我们的反向图。更新后的 inspect_backend 定义了一个前向 (fw) 和反向 (bw) 编译器捕获函数，该函数读取 AOTAutograd 中的前向和反向图，打印降低的 ATen IR，并保存前向和反向图的 FX 图。
这将生成以下的前向和反向图。请注意，前向图也与我们在图 x 中看到的略有不同。例如，FX 图 IR 中的 torch.sin(x) 和我们原始代码中的已被替换为 torch.ops.aten.sin.default()。你可能会问，这个叫做 aten 的奇怪东西是什么，如果你还不熟悉的话。ATen 代表一个张量库，它是一个非常富有创意命名的低级库，具有 C++ 接口，实现了许多在 CPU 和 GPU 上运行的基本操作。
在急切模式操作中，你的 PyTorch 操作会被路由到这个库，然后调用适当的 CPU 或 GPU 实现。AOTAutograd 自动生成代码，将较高级的 PyTorch API 替换为前向和反向图的 ATen IR，你可以在下面的输出中看到：
作者截图
你还可以看到，除了前向传递的输出，前向图还输出了一些额外的张量 [add, sin, cos, primals_1]。这些张量被保存用于反向传递的梯度计算。你也可以在之前分享的前向和反向传递的计算图中看到这一点。
PyTorch 中不同类型的 IR 有哪些？
ATen IR 是我们在上一节讨论的 ATen 库所支持的操作符列表，你可以在 ATen library here 查看实现的完整操作列表。PyTorch 中还有两个其他的 IR 概念你应该了解：1/ Core Aten IR 2/ Prims IR。Core Aten IR 是更广泛的 Aten IR 和 Prims IR 的一个子集，而 Prims IR 是 Core Aten IR 的一个更小的子集。假设你正在设计一个处理器并希望在你的硬件上支持 PyTorch 代码加速。要在硬件中支持完整的 PyTorch API 列表几乎是不可能的，因此你可以构建一个仅支持 Core Aten IR 或 Prims IR 中定义的较小基本操作符子集的编译器，并让 AOTAutograd 将复合操作符分解为核心操作符，正如我们将在下一节中看到的那样。
ATen IR、Core ATen IR 和 Prims IR 之间有什么区别？
作者截图
Core Aten IR（以前称为标准 Aten IR）是 Aten IR 的一个子集，可用于组合 Aten IR 中的所有其他操作符。针对特定硬件加速器的编译器可以专注于支持 Core Aten IR 并将其映射到其低级硬件 API。这使得为 PyTorch 添加硬件支持变得更容易，因为他们不必实现对完整 PyTorch API 的支持，而 PyTorch API 将继续随着更多抽象的增加而增长。
Prims IR 是 Core Aten IR 的一个更小的子集，它将 Core Aten IR 操作进一步分解为基本操作，使得针对特定硬件的编译器更容易支持 PyTorch。但将操作符分解为更低级别的操作将不可避免地导致性能下降，因为会增加额外的内存写入和函数调用开销。预期是硬件编译器可以将这些操作符融合回去，以支持硬件 API 并恢复性能。
尽管我们不需要进一步将我们的函数分解为 Core Aten IR 和 Prims IR，我将下面演示如何操作。
（可选话题）分解为 Core Aten IR 和 Prims IR
如果你正在设计硬件或硬件编译器，在硬件中支持完整的 PyTorch API 列表几乎是不可能的，尤其是考虑到深度学习和 AI 发展的速度。但硬件设计师的优势在于，大多数深度学习功能可以映射到很少的基本数学操作中，而计算密集型的操作是矩阵-矩阵和矩阵-向量操作。像 PyTorch API 支持的复合操作符可以使用 AOTAutograd 分解为这些基本操作，如我们将在本节中讨论的。如果你不涉及底层硬件，你可以跳过本节。
你可以更新 AOTAutograd 函数以传递一个分解字典，该字典可以将 Aten IR 降级到 Core Aten IR 和 Prims IR。我只会在这里分享相关的代码片段和输出，因为你可以在 GitHub 上找到完整的笔记本。默认情况下，操作符不会被分解为 Core Aten IR 或 Prims IR，但你可以传递一个分解字典。
在下面的代码片段中，我将我们的函数 f 转换为损失函数 f_loss，通过将均方误差 (MSE) 的计算包含到我们的函数中。我这样做是为了演示 AOTAutograd 如何将 MSE 分解为其基本操作符。
作者提供的截图
分解的输出是 mse_loss 被分解为更基本的操作：减法、幂(2)、均值。
作者提供的截图
这是因为两个向量 x 和 y 之间的 MSE 或均方误差定义为以下公式，只需要这 3 个操作：减法，其中幂是逐元素操作。如果你为你的硬件编写编译器，你可能已经支持这 3 个操作，通过分解，你的 PyTorch 代码可以在不做进一步修改的情况下运行。
你还可以在 FX 图形可视化中看到这一点
作者提供的截图
现在让我们进一步将其分解为 Prims IR，这是 ~250 个操作符的一个更小的子集。同样，我只会在这里分享相关的代码片段和输出，因为你可以在 GitHub 上找到完整的笔记本。
作者提供的截图
prim IR 分解的输出如下。所有标记为红色的 aten 操作都被替换或分解为使用绿色的 prim 操作符。
作者提供的截图
图优化：层和操作符融合以及 C++/GPU 代码生成
讨论的 PyTorch 技术： TorchInductor、OpenAI Triton（默认）其他编译器
在博客的最后一部分，我们将讨论使用 TorchInductor 进行操作符融合和自动代码生成。首先是一些基础知识：
深度学习优化编译器是什么？
深度学习的优化编译器擅长发现代码中的性能差距，并通过将代码转换以减少诸如内存访问、内核启动、针对目标后端的数据布局优化等代码属性来解决这些差距。TorchInductor 是 torch.compile() 的默认优化编译器，可以使用 OpenAI Triton 为 GPU 生成优化内核，并使用 OpenMP pragma 指令为 CPU 生成优化内核。
深度学习中的操作符融合是什么？
深度学习由许多基本操作组成，例如矩阵-矩阵和矩阵-向量乘法。在 PyTorch 的急切执行模式中，每个操作都会导致在硬件上进行单独的函数调用或内核启动。这会导致 CPU 启动内核的开销，并导致在内核启动之间进行更多的内存读写。像 TorchInductor 这样的深度学习优化编译器可以将多个操作融合为一个复合运算符，并为其生成低级 GPU 内核或 C++/OpenMP 代码。这样可以由于较少的内核启动和较少的内存读/写而实现更快的计算。
作者提供的截图
前一节中从 AOTAutograd 输出的计算图由许多在 FX 图中表示的 Aten 运算符组成。TorchInductor 优化不会改变图中的底层计算，只是通过运算符和层的融合对其进行重组，并且为其生成 CPU 或 GPU 代码。由于 TorchInductor 可以提前看到完整的前向和反向计算图，因此可以对不依赖彼此的操作进行乱序执行，并最大化硬件资源利用。
在底层，对于 GPU 目标，TorchInductor 使用 OpenAI 的 Triton 生成融合的 GPU 内核。Triton 本身是一个独立的基于 Python 的框架和编译器，用于编写优化的低级 GPU 代码，否则需使用 CUDA C/C++编写。但唯一的区别是，TorchInductor 将生成编译为低级 PTX 代码的 Triton 代码。
对于多核 CPU 目标，TorchInductor 会生成 C++代码，并插入 OpenMP 指示以生成并行内核。从 PyTorch 用户级别的角度来看，这是 IR 转换流程：
当然，这是高级视图，我在这里省略了一些细节，并鼓励您阅读TorchInductor 论坛帖子和OpenAI 的 Triton 博文。
现在我们将放弃在前一节中使用的虚拟编译器，并使用使用 TorchInductor 的完整 PyTorch 编译器栈。
注意，我已经传递了一个可选参数，启用了两个调试功能：
-
trace.enabled：生成中间代码以检查 TorchInductor 生成的代码 -
trace.graph_enabled：生成在运算符融合后的优化计算图可视化
对于我们的简单示例，TorchInductor 能够将函数中的所有中间操作融合到一个自定义运算符中，您可以在下面看到这如何简化了前向和反向计算图。
作者提供的截图
你一定好奇这个融合运算符在代码中是什么样子。融合运算符的代码是由 TorchInductor 自动根据目标设备 — CPU 或 GPU — 生成的，它不需要明确指定目标设备，可以根据数据和模型设备类型进行推断。
要查看生成的代码，你必须启用调试，使用 trace.enabled=True，这将创建一个名为 torch_compile_debug 的目录，包含调试信息。
前向和反向图代码的完整路径为：
-
torch_compile_debug/run_<DATE_TIME_PID>/aot_torchinductor/model__XX_forward_XX/output_code.py -
torch_compile_debug/run_<DATE_TIME_PID>/aot_torchinductor/model__XX_backward_XX/output_code.py
如果你设置 device = ‘cuda’（假设你的计算机有 GPU 设备），那么在 forward 文件夹生成的代码将是 Open AI Triton。
作者截图
如果你设置 device = ‘CPU’，那么生成的代码将是带有 OpenMP pragmas 的 C++ 代码。
作者截图
总结
深度学习编译器的内部机制复杂，有如瑞士表般精致。在本博文中，我希望给你提供一个关于这个主题的简明易懂的入门，并介绍这些技术如何驱动 PyTorch 2.0。
作者截图
-
我从一个简单的 PyTorch 函数开始。
-
我展示了 TorchDynamo 如何捕获图并在 FX IR 中表示它。
-
我展示了 AOTAutograd 如何生成反向传递图，将 PyTorch 运算符降低为 Aten 运算符，并在 FX 图容器中表示它。
-
我讨论了如何进一步将 Aten 运算符分解为 Core Aten IR 和 Prims IR，以减少其他编译器支持的运算符数量，而无需支持完整的 PyTorch API 列表。
-
我展示了 TorchInductor 如何执行运算符融合，并为 CPU 和 GPU 目标生成优化的代码。
如果你跟随我，应该能够对以下问题提供高层次的回答：
-
什么是深度学习编译器？
-
当你调用 torch.compile() 时，PyTorch 2.0 编译器会做什么？
-
为什么我们需要提前进行前向和反向传递图？
-
PyTorch 中有哪些不同的中间表示（IR）？
-
ATen IR、Core ATen IR、Prims IR 之间有何区别？
-
什么是运算符融合，以及它为何重要？
感谢你阅读到最后！
如果你觉得这篇文章有趣，请考虑在 medium 上关注我，以便在我发布新文章时收到通知。也请查看我在 medium 上的其他博文，或者关注我的推特 (@shshnkp)，LinkedIn 或在下方留言。如果你希望我写关于特定主题的文章，我很乐意听取你的建议！
当 1+1≠2：量子物理学如何打破统计学定律
原文：
towardsdatascience.com/how-quantum-physics-broke-the-laws-of-statistics-86fb8941ed2c
揭示 2022 年物理学诺贝尔奖背后的数据科学
·发表于 Towards Data Science ·阅读时间 12 分钟·2023 年 5 月 3 日
–
统计学是数据科学的核心支柱，但其假设并非总是经过充分检验。这在量子计算的崛起中更加明显，因为即使是统计学公理也可能被违反。在这篇文章中，我们将深入探讨量子物理学如何打破统计学，并揭示使用数据科学类比来理解它的方法。
2022 年的物理学诺贝尔奖归结于掷硬币，发现量子物理学违背了统计学的基本定律（背景： Dan Dennis，硬币 3D 模型： hyperionforge，组合由作者提供）
让我们玩一个掷硬币的游戏：掷三枚硬币，并尽量让它们都落在不同的面上。这看似是不可能完成的任务，因为无论硬币如何伪造，它最多只能有两个面。所有三次掷硬币的结果都不同的可能性实在太少。
然而，凭借量子物理学的力量，这样一个看似不可能的成就可以在统计上实现：三次掷硬币可以全部落在不同的面上。而获胜的奖励？2022 年物理学诺贝尔奖，授予了阿兰·阿斯佩、约翰·克劳斯和安东·蔡林格，颁奖日期为 2022 年 10 月 4 日。
根据 nobelprize.org，他们的成就
“对于纠缠光子的实验，确立了贝尔不等式的违反，并开创了量子信息科学。”
这句话充满了术语：纠缠光子、贝尔不等式和量子信息科学。我们需要一个更简单的、通俗的描述来解释这样一个重要的成就。以下是翻译：
科学家们通过展示量子物理学可以违背看似不可能的几率，证明了我们对世界的统计学视角存在缺陷。
这些不可能几率的细节被称为贝尔不等式的数学公式捕捉到。研究人员通过激光（使用纠缠光子束）展示了这些不可能的几率，而不是掷硬币。
这与数据科学有什么关系？由于我们的量子力学世界是数据的终极来源，统计规律中的缺陷可能会破坏数据科学的根基。如果统计学确实不完整，我们将无法信任从中得出的结论。
幸运的是，在我们的宇宙中，这些统计缺陷通常非常微小和可忽略。尽管如此，理解经典统计学需要如何修改仍然很重要，因为在遥远的未来，数据科学可能需要结合这些缺陷（例如，量子计算机）。
在回答量子物理如何违背统计规律之前，我们首先需要理解统计学如何有效地描述我们的世界。
以概率取代不可预测性
一堆硬币中有多少个正面/反面？答案大约是 50%。概率是混沌世界的近似模型（图片来源：Claudio Schwarz于Unsplash）
掷硬币，你得到正面/反面。然而硬币并不完全是随机的：一个控制完美的机器人可以严重操控掷硬币的结果。
50/50 的概率是什么意思？硬币的方向对其周围微小细节非常敏感。这使得预测硬币的着陆方向变得困难。因此，我们选择一个非确定性的结果，而不是解决非常复杂的方程来得出确定的结果。我们观察到典型的硬币在正面/反面方面相当对称。在没有任何特定偏差的情况下，50/50 的几率将是一个很好的近似值（尽管研究表明这些几率可以被改变，例如，Clark MP 等人）。
总结一下，
概率是建模复杂系统细节的近似值。复杂的物理学被不确定性取代，以简化数学运算。
从天气模式到经济和医疗保健，不确定性都可以追溯到复杂的动态。数学家将这些近似值转换为基于公理的严格定理，以帮助我们操作和从不可预测的结果中得出见解。
计算概率
统计公理帮助我们每天做出决策，这包括医学领域（照片由Towfiqu barbhuiya拍摄，来源于Unsplash）
量子物理学如何打破统计规律？它违反了可加性公理。
这个公理是如何工作的？让我们考虑一些常见的场景，其中我们使用统计数据来做出决策：
-
当外面下雨🌧时，我们会带上雨伞☔️。
-
当我们生病时，医生会开药💊帮助我们康复。
在雨天的场景中，尽管可能有数万亿种雨滴落下的方式，但这些可能性中的大多数会使我们变湿变冷，因此我们带上雨伞。
在医生场景中，给定一个诊断，有多个可能性：不同的疾病进展、副作用、恢复率、生活质量，甚至误诊……等等。我们选择能够带来最佳整体结果的治疗方案。
加法公理是一个形式化的声明，我们可以将概率分解成可能性：
这个公理是合理的，因为统计学旨在量化我们对系统的无知。就像我们给硬币翻转分配 50/50 一样，我们使用加法公理通过对其组成部分的所有可能轨迹进行平均来推导系统的属性。
尽管这一切听起来很直观，但这真的就是自然的运作方式吗？通过实验，我们可以确认宏观物体是这样工作的，但当我们放大到微观时会发生什么？它与宏观世界一样，亚原子演员从一个场景移动到下一个场景吗？还是更像是电影屏幕，上面的抽象像素在闪烁开/关，创造出一个故事的幻觉？
结果是，像素类比更为准确。当我们放大时，可能性的不确定路径变得更加模糊。因此，加法公理被违反。
我们的公理的替代是什么？它是量子物理学的定律。
量子可能性的计数
结果是，我们的基本量子世界与 AI 模型非常相似（照片由Steve Johnson拍摄，来源于Unsplash）
虽然量子物理学相当复杂，但我们可以通过数据科学类比来理解其要点。量子物理学基于线性代数，因此可以被看作是一个特殊的机器学习模型。
以下是与机器学习类比相关的关键量子公理：
-
世界由一长串（复杂的）数字描述，称为量子状态——类似于图像的像素值，或更抽象的机器学习嵌入向量。
-
随着时间的推移，这个量子状态会发生变化。这个更新可以通过将我们的量子状态传递通过类似神经网络的函数来计算，这个函数称为算符（技术上是单位矩阵）：
量子状态变化的示意图（作者提供）
继续我们的机器学习类比，我们可以把宇宙看作一个巨大的神经网络。每个算子代表一个（线性）网络层。通过这个网络，每次发生的交互都被印刻在我们宇宙的量子态上。这个计算从时间开始以来就不断运行。这是看待我们世界的一种深刻方式：
我们的一致现实源于量子态中的孤立分组。
我们对物体存在的宏观感知源于我们算子的特定神经网络连接。
听起来有些抽象，我们来考虑一个具体的例子：量子物理如何描述降落在我们头上的雨滴？
-
空气分子和我们在开放环境中的数据被捕捉在一个量子态中。
-
当水分子感受到地球的引力时，量子态会通过相应的算子进行更新。
-
在经历了类似神经网络的多层更新之后，量子态会获取一些特定的数值。
-
物理定律决定了这些数字往往会形成簇。这些簇中的一些转化为这些雨滴的一致存在，这最终与我们的神经元感知到这些雨滴相关联。
从现代观点来看，没有理由认为加法公理必须成立。因为
类似于机器学习黑箱，我们不总是能够跟踪量子态的所有物理属性。因此，物理结果并不总是附带一系列中间可能性。
在雨滴场景中，这意味着我们不能总是找到导致特定水分子下落的量子态中的具体数字。事实上，量子态通常包含多个位置的分子数据（例如，叠加态），我们对其物理位置的感知可能是所有这些数据的复杂总和。
这可能看起来很矛盾，因为我们在日常生活中并没有感觉到奇怪的差异和叠加态！但原因是这些差异非常微小，它们的微小性可以通过技术退相干理论来证明，这超出了我们的范围（虽然这里是我的一篇文章，可能会提供一些启示）。
尽管如此，小并不等于零。量子效应有时可能是显著的，并且可能导致看似不可能的统计数据。
如何做呢？让我们来找出答案。
测试 1+1 = 2
为了使普通统计学定律无效，我们需要考虑一些简单但不可能的情景。其中最简单的涉及 3 枚硬币。
想象有 3 个机器人进行 3 次独立的掷币实验。在经典统计学中，我们可以利用加法公理完全指定统计数据：通过列出所有 8 种结果及其概率（注意：机器人/硬币可能被操控）：
翻转 3 枚硬币的所有 8 种可能性（作者提供）
从实验上来看，我们可以通过重复这些抛币来测量这些概率。
无论概率选择如何，都存在一个理智约束：一个硬币只有 1+1 = 2 面，因此当我们抛 3 枚硬币时，至少有 2 枚会落在相同面上。因此，如果我们随机（均匀地）选择一对硬币进行检查，我们应该预期至少有 1/3 的机会观察到它们是相等的。
让我们尝试一些示例，将三枚硬币标记为A、B、C
-
如果 3 枚硬币都是公平和独立的，那么我们选择一对相等的硬币的机会是 1/2。
-
如果A = B，但A ≠ C。无论A 如何落地，只有一对是相等的。选择这一对的机会是 1/3。
我们看到相同对的概率始终至少为 1/3。这可以总结为贝尔不等式（参考 L. Maccone 的论文）
贝尔不等式的最简单示例（作者提供的图像）
尽管测试如此显而易见的事物可能看起来荒谬，但事实证明，这一不等式实际上可以被*违反——*这证明它们并非如此显而易见。
诺贝尔级别的抛币
量子抛币通常使用类似于上述设置的激光（维基媒体：Adelphi Lab Center，许可证）
为了观察贝尔不等式的违反，物理学家不能仅依赖传统的硬币。相反，他们需要利用由激光制成的量子硬币，这些硬币具备了抛币所需的所有要素：
-
抛硬币：将激光发送到光束中
-
观察正面/反面：在两个探测器中的一个上获取读数*
-
随机性：读数通常是不可预测的，除非进行操控
(* 如果没有探测器观察到任何内容，可能会有故障读数)
现在，我们可以将激光设置为不同的方向，以模拟 3 次不同的抛币。那么量子硬币如何处理不可能的情况呢？如果我们观察三次抛币的字面结果，看到三种不同的结果在逻辑上是不可能的。
这就是贝尔不等式的作用：它将关于 3 枚硬币的逻辑陈述分解为仅涉及 2 枚硬币的概率陈述。因此，如果我们抛 3 枚硬币，但每次仅观察 2 枚，那么在保持逻辑的同时，可以违反统计规律。在量子物理学中，抛硬币与观察硬币遵循两个不同的相互作用：
- 量子：抛硬币和观察硬币由两个不同的算符来支配。尚未被观察的抛币不需要被赋予一个确定的结果*。
这与经典统计学形成对比
- 经典：硬币被掷出时，正面/反面已确定。这是由加法性公理保证的。无论我们是否决定观察它，这一点都不影响结果。
(*这就是“远距幽灵行动”发挥作用的地方，因为任何时刻任何人都可以启动探测器来观察第三个硬币，从而破坏我们的结果。)
那么我们如何进行实验呢？我们需要准备我们的硬币处于特定的量子态。这里，我们构建一个系统，其中三个硬币的量子态可以由平面上的三个向量表示，如下图所示*：
(* 从技术上讲，量子态涉及更复杂的纠缠光子，但为了简洁起见，我们将跳过这些细节)
违背贝尔不等式的三硬币量子态（作者）
两次掷币结果相同的概率是多少？答案来自物理学，并被设计为余弦相似度的平方：
相同掷币概率（作者）
现在，如果我们随机选择一对量子硬币进行检查*，它们相同的概率只有 1/4；这比逻辑上的 1/3 保证还要低！
(*实验需要设置成这样，选择是在硬币被掷出之后进行的，以便排除粒子和设备之间的幽灵勾结)
用我们的贝尔不等式重新表述，我们有
贝尔不等式违背（作者）
我们的理智检查被违背了！如果我们假装经典统计仍然适用，这将意味着至少 1/4 的时间，所有三次掷币结果都不同！
请注意，虽然我们的三硬币实验易于理解，但在实验结果中存在实验难度和潜在漏洞。因此，典型实验往往涉及更多的掷币和更复杂的观察（例如，GHZ 实验由潘建伟等进行）。
新的统计显著性
那么，我们看到量子概率有时会导致意想不到的结果，这有什么大不了的，我们为什么要在意？
首先，让我们从实际应用开始。随着技术不断推进，将更多计算能力压缩到更小的体积中，量子物理学将变得越来越重要。最终，我们的计算范式需要进行彻底的改革，以充分利用量子设备。因此，虽然贝尔不等式的违背可能很微妙，但这表明我们在设计量子算法时需要仔细思考。
其次，这些违背暴露了传统统计推理的基本限制。例如，如果有人赢得了彩票，将其原因归因于彩票球以特定方式出来是完全合理的。然而，我们不能放大并因果关联中奖彩票与房间内所有分子（量子）状态。因此，我们的因果推理统计理论有一个物理极限！
最后，量子效应挑战我们重新思考我们的宇宙。尽管量子物理已经被反复验证，但它仍然可能只是一个近似值。未来，我们可能会发现它被更抽象的基本法则所取代。
作为历史教训，甚至爱因斯坦也被量子物理的怪异性所劝阻，以至于他宣称“上帝不会掷骰子”来拒绝它。然而，量子物理继续取得胜利，并在推进我们现代技术和对世界的理解中发挥了基础性作用（参见我的文章）。
总结来说，量子物理主宰了世界，而 2022 年的物理学诺贝尔奖突显了其与统计学和数据科学的深刻联系。虽然量子物理不常被教授，但我们都应该努力理解并接受它的重要性。
如果你喜欢这篇文章，你可能会喜欢下面这些相关的文章。请留下评论/反馈，因为我花费了无数小时来完善上述许多物理和数据科学的解释，我很想听听你的想法。如果你想支持我，只需分享我的文章，并帮助我传播科学知识。👋
从我们的存在到计算机芯片的工作方式，量子物理对于理解宇宙至关重要。
medium.com](https://medium.com/physicist-musings/quantum-physics-is-everywhere-and-we-should-all-learn-it-26c9b1dedd2a?source=post_page-----86fb8941ed2c--------------------------------) [## 量子物理没有所有怪异性：一致性历史方法
量子物理不一定是怪异的。在这里，我们探讨了一种现代且直观的量子系统解释……
将现实描述为模拟大大低估了我们世界的复杂性。这里是为什么模拟……
我们都听过“相关性不意味着因果性”这句话，但没有人真正谈论因果性到底是什么……
足球进球的随机性如何？
原文：
towardsdatascience.com/how-random-are-goals-in-soccer-8a822c1f3bc
体育分析
通过频率统计理解进球事件
·发表于 Towards Data Science ·7 min read·2023 年 7 月 29 日
–
足球进球——照片由 Chaos Soccer Gear 提供，来自 Unsplash
足球（对于美国读者来说是 soccer）是一项了不起的运动。它不可能凭借巧合成为世界上最受欢迎的运动。
足球将人们聚集在一起，它是从忙碌生活中断开的借口，因为比赛时间是娱乐时间。我们点一些快餐，边吃边看梅西用球创造魔法——我们能享受他的表现真是太幸运了。我们还能观看像 2010 年的巴萨或 2023 年的曼城这样许多了不起的球队。
许多人会说没有两场比赛是相同的。这就是足球，没有其他运动可比。但是我认为这是错误的。
尽管它如此杰出，但它仍然受数学的支配。就像其他一切一样。
生活充满了数学模型。足球也不例外。
我一生都是巴萨的忠实粉丝。再加上我目前的职业状况，结果就是对体育分析产生了真正的兴趣——显然是倾向于足球。
这篇文章是我关于体育分析的第一篇，因此我会保持相对简单。不过，我计划写更多的文章，深入了解数学如何应用于足球（以及其他运动如手球）——并与大家分享这些见解。
数据科学家被聘用从事体育分析工作的数量正在强劲增长，这种趋势似乎不会很快停止。利用数据进行体育分析比以往任何时候都更有意义，尤其是考虑到数据生成的速度也在快速增长。
所以，这篇文章将是所有有志于成为体育分析师或对体育感兴趣的数据相关人员的绝佳入门工具。
在这里，我将使用 StatsBomb 的[1] 开放且免费的数据[2] 来检查我随机选择的 2015–2016 赛季的西甲。我邀请你做同样的分析，看看是否对其他赛季和联赛也适用。
那么，让我们开始深入挖掘吧！
准备数据
现在有一个很棒的 Python 模块可以让我们获取所有需要的数据：statsbombpy[3]。我们首先需要做的是逻辑上安装它：
pip install statsbombpy
然后打开你的 Python 文件或笔记本，开始导入下列模块：
import pandas as pd
from statsbombpy import sb
import seaborn as sns
接下来，我们需要获取 2015–16 赛季的所有西甲比赛：
competition_row = sb.competitions()[
(sb.competitions()['competition_name'] == 'La Liga')
& (sb.competitions()['season_name'] == '2015/2016')
]
competition_id = pd.unique(
competition_row['competition_id']
)[0]
season_id = pd.unique(
competition_row['season_id']
)[0]
matches = sb.matches(competition_id=competition_id, season_id=season_id)
这就是比赛数据框的样子：
数据框的前 5 行和前 10 列——图片由作者提供
现在，我们要检查进球情况，但没有一列显示比赛中的进球数。我们可以通过将 home_score 和 away_score 列相加来创建它：
matches['goals'] = matches['home_score'] + matches['away_score']
我们现在准备开始分析。
之前的进球是否会影响未来的进球？
在进行任何计算之前，我喜欢先可视化数据。这是理解你所处理数据的最佳方式。让我们构建一个直方图：
import seaborn as sns
sns.histplot(
x='goals',
data=matches,
bins=matches['goals'].nunique(),
binwidth=0.9
)
西甲 2015/16 赛季的进球分布——图片由作者提供
正如预期的那样，大多数比赛中进球数在 1 到 4 个之间。我们都知道最常见的结果总是 1–0、2–0、2–1、3–0、3–1（客队同样如此）。有一个进球数为 12 的比赛案例是明显的异常值——皇马对阵拉约·巴列卡诺，最终比分为 10–2。
到目前为止没有新内容，仅仅是有趣的事实。
该赛季的平均进球数为 2.74（四舍五入）。这转化为每分钟平均 0.030497（四舍五入）个进球（或 1/32.79）。换句话说，我们可以说在每个 90 分钟的时间段内有 0.030497 的进球概率。
有趣的部分来了：我们将模拟整个赛季。由于我们只关注比赛进球——而不是比赛胜负或赛季排名——我们将使用的唯一参数是进球的分钟概率。
# Creating a simulation
import random
mean_goals = matches['goals'].mean()
def simulate_match():
goals = 0
for i in range(90):
goals += np.random.choice(
np.arange(0, 2),
p=[1-(mean_goals/90), mean_goals/90]
)
return goals
def simulate_season(n_games):
goals_per_game = []
for i in range(n_games):
goals_per_game.append(simulate_match())
return goals_per_game
现在我们有了 simulate_match() 和 simulate_season() 函数，它们每个都包含一个简单的循环（即在每场比赛的每分钟循环），然后随机计算这一分钟是否有进球。
为此，我们使用 random.choice()，它以指定的概率从 0 或 1 中选择（~0.030497 选择 1，~0.969503 选择 0）。
使用我们之前看到的直方图，让我们绘制模拟分布线叠加在上面：
goals_per_game = simulate_season(len(matches))
mu = np.var(goals_per_game)
pmf = poisson.pmf(goals_per_game, mu)
pmf *= (most_repeated_count/pmf.max())
sns.lineplot(
x=goals_per_game,
y=pmf,
color='red',
label="Simulated"
).set(ylabel='count')
sns.histplot(
x='goals',
data=matches,
bins=matches['goals'].nunique(),
binwidth=0.9
).set(title='Goal vs Simulated Poisson distribution')
西甲 2015/16 赛季的进球和模拟的泊松分布——图片由作者提供
根据模拟数据，实线红色曲线显示了泊松分布。当之前事件的时间对未来事件没有影响时，这种分布适用。
换句话说，使用维基百科的话：“泊松分布表示在固定时间或空间间隔内发生给定数量事件的概率，如果这些事件以已知的恒定平均速率发生，并且独立于上一个事件以来的时间。”[4]
这个图表的要点非常直接：这是一个通过随机模拟生成的模型，它与实际结果非常接近。
所以我们有了答案：进球不受目前进球数或比赛时间的影响。 它们在时间上是不可预测的。
正如 David Sumpter 在《Soccermatics》[5]中所说：“正是足球比赛从一刻到下一刻的不可预测性，在 90 分钟后产生了泊松分布。我们知道比赛中的平均进球数，但它们的时机是不可预测的。因此，一些比分比其他比分更可能出现。这里的悖论是，比分是由随机性解释的。进球时间的高度随机性使得结果中的模式是可预测的。”
结论
目标本身并非随机的，它们发生是有原因的。但它们确实是不可预测的。当某些条件满足时，我们才会有目标：防守者犯错，进攻者传球顺利，梅西得到球，守门员无法阻止球进门……
一次关键的比赛瞬间可以在任何时候完全改变比赛的方向。我们今天已经用数学证明了这一点——进球不受之前进球或发生时间的影响。
但我们不应局限于足球来得出结论。我们所做的这个简单分析也可以应用于任何与足球无关的数据。我不仅仅是在谈论冰球进球、手球进球或篮球得分。
数据科学家可以期望在合理假设事件可以意外发生、在任何时间、独立于之前发生的事件数量时，出现泊松分布。
你看到这有多强大了吗？一个简单却强大的数学建模工具。
如果你不相信我，可以尝试用你喜欢的任何数据进行分析。例如，试着分析某一天有多少人去超市、工厂中发生的事故数量，或跨越两个国家边界的汽车数量。你会发现它们也无法逃脱泊松分布。
再次强调，这对数据科学家来说是一个理解我们所处理数据的强大工具。
回到足球上来，显然它不是一项随机的运动。它远不止于此：它涉及挫折、逆转、阵容、战术、技能……
如果我们不想深入探讨，随机性可能足以研究一个赛季中的进球数。但足球远比仅仅解释进球分布更有趣。
在我接下来的足球分析帖子中，我们将进一步探讨超越简单随机性的内容。
**Thanks for reading the post!**
I really hope you enjoyed it and found it insightful.
Follow me and subscribe to my mailing list for more
content like this one, it helps a lot!
**@polmarin**
如果你想进一步支持我，可以通过下面的链接订阅 Medium 的会员：这不会增加你的额外费用，但会帮助我完成这个过程。
[## 通过我的推荐链接加入 Medium - Pol Marin
阅读 Pol Marin（以及 Medium 上其他数千名作家）的每一个故事。你的会员费用直接支持 Pol…
medium.com](https://medium.com/@polmarin/membership?source=post_page-----8a822c1f3bc--------------------------------)
资源
[1] StatsBomb
[4] 泊松分布 — 维基百科
比率的可靠性如何？
原文：
towardsdatascience.com/how-reliable-is-a-ratio-1467f50e943d
学习如何使用 Python 中的经验贝叶斯分析评估比率的可靠性
·发表于Towards Data Science ·9 分钟阅读·2023 年 12 月 8 日
–
图片来源：rupixen.com 在 Unsplash
介绍
在数据科学领域，我的一位参考对象是Julia Silge。在她的干净星期二视频中，她总是制作代码跟随类型的视频，教授/展示给定的技术，帮助其他分析师提升技能并将其纳入自己的技能库。
上周二的话题是经验贝叶斯（她的博客文章），引起了我的注意。
那是什么呢？
经验贝叶斯
经验贝叶斯是一种统计方法，用于处理像成功/总尝试数这样的比率。当我们处理这些变量时，常常会遇到 1/2 的成功率，这转化为 50%的成功百分比，或者 3/4（75%），0/1（0%）。
这些极端百分比并不代表长期的现实，因为尝试次数太少，很难判断是否存在趋势，大多数情况下这些案例被忽略或删除。需要更多的尝试来确定真实的成功率，比如 30/60，500/100，或其他对业务有意义的数据。
然而，使用经验贝叶斯，我们能够利用当前的数据分布来计算其自身数据在早期或晚期阶段的估计值，正如我们将在本文接下来的部分看到的那样。
我们使用数据分布来估计每个观察值在早期和晚期阶段的比率。
分析
让我们跳到分析部分。需要遵循的步骤是：
-
加载数据
-
定义成功并计算成功比率
-
确定分布的参数
-
计算贝叶斯估计
-
计算可信区间
让我们继续。
导入
# Imports
import pandas as pd
import numpy as np
import scipy.stats as scs
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px
from distfit import distfit
数据集
本文中使用的数据集是客户及其电子商务订单和实体店订单的日志。它是基于真实数字生成的，但所有数据集都经过了缩减、随机化和修改，以使其匿名化以供公开使用。
变量如下：
-
customer: 客户 ID
-
brick_mortar: 实体店发生了多少次交易。
-
ecomm: 发生了多少次电子商务交易。
数据的头部。图像由作者提供。
成功比例
假设我们的公司愿意增加其在线渗透率。因此，在这个分析中，我们的成功定义为当客户在线购买时。 客户在线购买[成功]，实体店购买[不成功]。
这样，我们的指标可以是简单的电子商务交易占总数的*[ecomm]/[ecomm + brick_mortar]*。这将给我们一个比例。我们现在的问题是：这两个客户：808 vs. 672 有多大不同。
808 的比例是 1/11（9%），而 672 的比例是 1/2（50%）。哪个更好？我应该把更多时间投资在哪里以促进增长？好吧，客户 672 的渗透率为 50%，但总共只有 2 次交易。客户 808 与公司有更长的历史，但刚刚开始在线购买。
808 还是 671？投资在哪里？图像由作者提供。
让我们计算所有的比例。在这里，我们将忽略那些从未在线购买过的客户，因为无法从中计算出估计值。
# Creating success ratio
df2c = (
df
.query('ecomm > 0 & brick_mortar > 0')
.assign(total_trx = df['brick_mortar'] + df['ecomm'],
ratio = df['ecomm']/(df['brick_mortar'] + df['ecomm']) )
)
display(df2c)
我们得出了结果。
计算出的比例。图像由作者提供。
接下来，让我们找出分布参数。
找出分布参数
这里的下一步是计算我们比例的分布参数。但在深入之前，让我们建立对为什么需要执行这一步骤的直觉。
好吧，托马斯·贝叶斯定理说概率是条件性的。因此，某事发生的可能性应该基于在类似情况下发生的先前结果。
如果某事发生了 10 次，那么它在第 11 次发生的可能性比某事在第 10 次未发生的情况下发生的可能性要大。
我们将使用贝塔分布来表示我们的先验期望，或到目前为止发生的情况。这就像观察我们的数据并理解模式，以便能够对任何阶段的任何数据点给出更精确的估计，考虑到它会遵循相同的分布。
为了知道哪个α和β使得贝塔分布适合我们的数据，我们将使用 Python 中的distfit模块。
看看分布，我们可以直接绘制直方图。
# Looking at the distribution
px.histogram(df2c,
x='ratio', nbins=20,
template="simple_white", width = 1000,
title='Histogram of the Ratio Brick & Mortar vs e-Comm')
比例分布。图像由作者提供。
现在，让我们看看哪些参数使得贝塔分布适合数据。
# Our distribution
X = df2c['ratio']
# Alternatively limit the search for only a few theoretical distributions.
dfit = distfit(method='parametric', todf=True, distr=['beta'])
# Fit model on input data X.
dfit.fit_transform(X)
# Print the bet model results.
dfit.model
--------------
[OUT]
[distfit] >INFO> fit
[distfit] >INFO> transform
[distfit] >INFO> [beta] [0.11 sec] [RSS: 0.823271] [loc=0.011 scale=0.941]
[distfit] >INFO> Compute confidence intervals [parametric]
{'name': 'beta',
'score': 0.8232713059833795,
'loc': 0.011363636363636362,
'scale': 0.9411483584238516,
>>>'arg': (0.850939343634336, 1.4553354599535102),<<<
'params': (0.850939343634336,
1.4553354599535102,
0.011363636363636362,
0.9411483584238516),
'model': <scipy.stats._distn_infrastructure.rv_continuous_frozen at 0x7838c17ad570>,
'bootstrap_score': 0,
'bootstrap_pass': None,
'color': '#e41a1c',
'CII_min_alpha': 0.030238213140192628,
'CII_max_alpha': 0.8158034848017729}
我标记了这行>>>arg: (0.850939343634336, 1.4553354599535102),<<<，这些是我们需要的α和β。
如果我们想查看拟合情况，可以使用这段代码，其中我们创建了一个包含 3 个子图的图形，并添加了概率密度函数、累积分布函数和 QQ 图。
import matplotlib.pyplot as plt
# Create subplot
fig, ax = plt.subplots(1,3, figsize=(18, 7))
# Plot PDF with histogram
dfit.plot(chart='PDF', ax=ax[0])
# Plot the CDF
dfit.plot(chart='CDF', ax=ax[1]);
dfit.qqplot(X, ax=ax[2])
DistFit 生成的图表。图片由作者提供。
好的。现在来计算我们的估计值。
贝叶斯估计
为了计算估计值，我们需要设置刚刚发现的 alpha（a）和 beta（b）。
# Alpha and Beta values
a,b = 0.850939343634336, 1.4553354599535102
使用 Beta 分布的估计参数，我们将用经验贝叶斯方法对电子商务交易比例进行建模。我们将从总体先验开始，然后根据个体证据进行更新。计算方法很简单，就是将 a 添加到成功次数（电子商务订单）中，再将 a+b 添加到订单总数中。
# Calculating Bayes Estimates
df2c = (df2c
.assign(emp_bayes = (df2c['ecomm'] + a)/(df2c['total_trx'] + a + b) )
)
很好！现在我们已经可以绘制简单比例与估计比例的图表，看看模型的效果如何。
# Scatterplot
fig = sns.scatterplot(data = df2c,
x= 'ratio', y= 'emp_bayes',
s=50)
# Slope line
plt.plot([0,1], [0,1], linestyle='--', color='red');
结果如图所示。
模型估计值与简单比例的对比。图片由作者提供。
这里，点越接近红线，比例就越可靠。这意味着简单比例[成功]/[总数]非常接近贝叶斯估计，贝叶斯估计已经考虑了先验概率。
看一下表中的一些数字。
估计值。图片由作者提供。
查看表格，我们发现电子商务交易数量和总交易数量越高，贝叶斯估计值与实际比例越接近。这意味着那些 1/2 的情况将大大偏离目标，因此作为趋势的可靠性较低。
我们还可以计算估计值的置信区间。这就是我们接下来要做的。
置信区间
置信区间是使我们的分析更为完整的一步。我们可以给决策者提供一个信任范围。
首先，让我们计算每个观察值的 alpha 1 和 beta 1，这些是更新先验后的 Beta 分布的后验形状参数。
# Calculate a1 and b1
df2 = (df2c
.assign(a1 = df2c.ecomm + a,
b1 = df2c.total_trx - df2c.ecomm + b,)
)
接下来，我们可以使用scipy.stats.beta模块，通过interval方法进行计算。95%置信区间的输出是一个包含两个数字的元组，其中索引[0]是下界，[1]是上界。
# Calculating the CI and add to the dataset
df2 = (df2
.assign( low = scs.beta.interval(.95, df2.a1, df2.b1)[0],
high = scs.beta.interval(.95, df2.a1, df2.b1)[1] )
.sort_values('total_trx', ascending=False)
)
最后，通过观察绘制区间图。我们创建一个图形，然后为中点绘制散点图和两个条形图，一个表示低范围，颜色为白色（在白色背景上看不见），另一个表示高范围。
# customer ID to string for better visualization
df2['customer'] = df2.customer.astype('string')
# Number of customers to plot
n_cust = 100
# Plot
plt.figure(figsize=(25,7))
plt.scatter(df2.customer[:n_cust], df2.emp_bayes[:n_cust], color='black', s=100 )
plt.bar(df2.customer[:n_cust], df2.low[:n_cust], label='low', color='white', zorder=1 )
plt.bar(df2.customer[:n_cust], df2.high[:n_cust], label='high', zorder=0 )
plt.legend();
这里是结果。
按电子商务订单数量排序的置信区间。电子商务订单越多，区间越小。图片由作者提供。
这个图表按电子商务订单数量递减排序。从多（左）到少（右），我们注意到置信区间的大小只会变大，这显示出当电子商务交易数量过少时，预测可靠估计真的很困难。
回到我们的问题之一：我应该投入更多时间在客户 672 还是 808 的电子商务增长上，以下是这些区间。
客户 808 和 672 的估计。图像由作者提供。
客户 672 的置信区间要大得多，因此不够可靠。我们应该花更多时间在开发客户 808 上，这与品牌的历史关系一致。
客户 808 和 672 的区间。图像由作者提供。
在你离开之前
好的，这个工具表现出强大的功能。我们看到它在我们的比较案例中带来了良好的结果。很多时候，我们会直接淘汰交易少的客户。但很多时候，这可能不太可能。此外，我们实际上可能希望分析小型客户。或者，它也可以帮助我们比较两个下单量大的客户，基于 CI 查看他们的潜力。
客户 2360 的置信区间可能达到 78%，而客户 361 的置信区间为 75%。图像由作者提供。
经验贝叶斯对于较大的数据集是一个不错的主意，它可能是其他模型的额外变量。
好吧，我敢打赌你可能已经在考虑这在你工作中的实用性了。那么，我们到此为止。
如果你喜欢这些内容，请关注我的博客，订阅我的页面，并在发布后获取帖子。
在 Medium 上阅读 Gustavo Santos 的文章。数据科学家。我从数据中提取洞察，帮助人们和公司…
gustavorsantos.medium.com](https://gustavorsantos.medium.com/?source=post_page-----1467f50e943d--------------------------------)
也可以在 LinkedIn 上找到我。
这是完整的代码，供参考。
[## Studying/Python/statistics/Empirical_Bayes.ipynb at master · gurezende/Studying]
这是一个包含我的测试和新包研究的库 - Studying/Python/statistics/Empirical_Bayes.ipynb 在…
github.com](https://github.com/gurezende/Studying/blob/master/Python/statistics/Empirical_Bayes.ipynb?source=post_page-----1467f50e943d--------------------------------)
参考
[## 经验贝叶斯用于 #TidyTuesday 《神秘博士》剧集 | Julia Silge]
数据科学博客
这两个比例哪个更高：10 中 4，还是 1000 中 300？这听起来像是个愚蠢的问题。显然……
本系列之前
免责声明：我不是统计学家，而是软件工程师。我对统计学的大部分知识来自于……
## scipy.stats.beta - SciPy v1.11.4 手册
一个贝塔连续随机变量。作为 rv_continuous 类的一个实例，beta 对象从中继承了一个集合……
Self-RAG 如何革新工业 LLMs
原文：
towardsdatascience.com/how-self-rag-could-revolutionize-industrial-llms-b33d9f810264
说实话——普通的 RAG 真的很笨。没有保证返回的响应是相关的。了解 Self-RAG 如何显著帮助
·发表在 Towards Data Science ·7 分钟阅读·2023 年 11 月 14 日
–
Self-RAG 演示 | Skanda Vivek
大型语言模型（LLMs）已经准备好革新各个行业。以金融行业为例，LLMs 可以用来分析大量文档，并在短时间内以低于分析师完成同样任务的成本找到趋势。但问题在于——你得到的答案很多时候只是部分的和不完整的。举个例子，假设你有一份包含公司 X 在过去 15 年的年收入的文档，但这些信息分布在不同的部分。在下面所示的标准检索增强生成（RAG）架构中，你通常会检索前 k 个文档，或选择固定上下文长度内的文档。
RAG 原型 | Skanda Vivek
然而，这可能会有几个问题。一个问题是前 k 个文档不包含所有的答案——例如可能只对应于过去的 5 年或 10 年。另一个问题是计算文档片段与提示之间的相似性并不总是能得到相关的上下文。在这种情况下，你可能会得到错误的答案。
一个实际的问题是你已经开发了一个普通的 RAG 应用程序，它在你测试的简单案例中工作良好——但当你将这个原型展示给利益相关者并且他们提出一些超出常规的问题时，它就会失败。
这就是 Self-RAG 大显身手的地方！作者开发了一种巧妙的方法，让经过微调的语言模型（Llama2–7B 和 13B）输出附加在 LM 生成的特殊标记 [Retrieval]、[No Retrieval]、[Relevant]、[Irrelevant]、[No support / Contradictory]、[Partially supported]、[Utility] 等，用于决定上下文是否相关/无关，上下文中的 LM 生成文本是否得到支持，以及生成的实用性。
训练 Self-RAG
Self-RAG 通过2 步层次化过程进行了训练。在第一步中，训练了一个简单的语言模型（LM），用于对生成的输出（无论是仅有提示还是提示+RAG 增强输出）进行分类，并在末尾附加相关的特殊标记。这个“评判模型”由 GPT-4 注释进行训练。具体来说，GPT-4 通过特定类型的指令进行提示（“给定指令，判断从网络上查找一些外部文档是否有助于生成更好的回应。”）
在第二步中，生成模型使用标准的下一个标记预测目标，学习生成延续内容，以及检索/评估生成内容的特殊标记。与其他微调或 RLHF 方法不同，这些方法可能会影响模型输出并使未来生成的内容产生偏差，通过这种简单的方法，模型仅被训练为在适当的时候生成特殊标记，而不会改变基础 LM！这真是太棒了！
评估 Self-RAG
作者对公共健康事实验证、多项选择推理、问答等进行了大量评估。共有 3 种任务。封闭集任务包括事实验证和多项选择推理，准确率被用作评估指标。短文本生成任务包括开放领域的问答数据集。作者评估了模型生成中是否包含黄金答案，而不是严格要求精确匹配。
长文本生成包括传记生成和长篇问答。对于这些任务的评估，作者使用 FactScore 来评估传记——基本上是生成的各种信息及其事实正确性的度量。对于长篇问答，使用了引用精度和召回率。
Self-RAG 在非专有模型中表现最佳，在大多数情况下，13B 参数的模型优于 7B 模型。在某些情况下，它甚至优于 ChatGPT。
推理
对于推理，self-RAG 存储库 建议使用 vllm——一个用于 LLM 推理的库。
在安装 vllm 后，你可以按如下方式加载库并进行查询：
from vllm import LLM, SamplingParams
model = LLM("selfrag/selfrag_llama2_7b", download_dir="/gscratch/h2lab/akari/model_cache", dtype="half")
sampling_params = SamplingParams(temperature=0.0, top_p=1.0, max_tokens=100, skip_special_tokens=False)
def format_prompt(input, paragraph=None):
prompt = "### Instruction:\n{0}\n\n### Response:\n".format(input)
if paragraph is not None:
prompt += "[Retrieval]<paragraph>{0}</paragraph>".format(paragraph)
return prompt
query_1 = "Leave odd one out: twitter, instagram, whatsapp."
query_2 = "Can you tell me the difference between llamas and alpacas?"
queries = [query_1, query_2]
# for a query that doesn't require retrieval
preds = model.generate([format_prompt(query) for query in queries], sampling_params)
for pred in preds:
print("Model prediction: {0}".format(pred.outputs[0].text))
对于需要检索的查询，你可以按下面的示例提供必要的信息作为字符串。
paragraph="""Llamas range from 200 to 350 lbs., while alpacas weigh in at 100 to 175 lbs."""
def format_prompt_p(input, paragraph=paragraph):
prompt = "### Instruction:\n{0}\n\n### Response:\n".format(input)
if paragraph is not None:
prompt += "[Retrieval]<paragraph>{0}</paragraph>".format(paragraph)
return prompt
query_1 = "Leave odd one out: twitter, instagram, whatsapp."
query_2 = "Can you tell me the differences between llamas and alpacas?"
queries = [query_1, query_2]
# for a query that doesn't require retrieval
preds = model.generate([format_prompt_p(query) for query in queries], sampling_params)
for pred in preds:
print("Model prediction: {0}".format(pred.outputs[0].text))
[Irrelevant]Whatsapp is the odd one out.
[No Retrieval]Twitter and Instagram are both social media platforms,
while Whatsapp is a messaging app.[Utility:5]
[Relevant]Llamas are larger than alpacas, with males weighing up to 350 pounds.
[Partially supported][Utility:5]
在上述示例中，对于第一个查询（与社交媒体平台相关），段落上下文是不相关的，如检索开始时的 [Irrelevant] 令牌所示。然而，外部上下文与第二个查询（与美洲驼和羊驼相关）相关。正如你所见，它在生成的上下文中包含了这些信息，并用 [Relevant] 令牌标记。
但在下面的示例中，上下文 “I like Avocado.” 与提示无关。如下面所示，模型预测在两个查询中都开始为 [Irrelevant]，并仅使用内部信息来回答提示。
paragraph="""I like Avocado."""
def format_prompt_p(input, paragraph=paragraph):
prompt = "### Instruction:\n{0}\n\n### Response:\n".format(input)
if paragraph is not None:
prompt += "[Retrieval]<paragraph>{0}</paragraph>".format(paragraph)
return prompt
query_1 = "Leave odd one out: twitter, instagram, whatsapp."
query_2 = "Can you tell me the differences between llamas and alpacas?"
queries = [query_1, query_2]
# for a query that doesn't require retrieval
preds = model.generate([format_prompt_p(query) for query in queries], sampling_params)
for pred in preds:
print("Model prediction: {0}".format(pred.outputs[0].text))
Model prediction: [Irrelevant]Twitter is the odd one out.[Utility:5]
[Irrelevant]Sure![Continue to Use Evidence]
Alpacas are a much smaller than llamas.
They are also bred specifically for their fiber.[Utility:5]
重点总结
Self-RAG 相比于普通 LLM 具有几个优势。
-
自适应段落检索：通过这一点，LLM 可以持续检索上下文直到找到所有相关上下文（当然是在上下文窗口内）。
-
更相关的检索：很多时候，嵌入模型在检索相关上下文方面表现不佳。Self-RAG 通过相关/不相关特殊令牌潜在地解决了这个问题。
-
超越其他类似模型：Self-RAG 超越了其他类似模型，并且在许多任务中意外地超过了 ChatGPT。对比 ChatGPT 未经过训练的数据——即更多的专有工业数据，将会非常有趣。
-
不改变基础 LLM：对我来说，这是一个巨大的卖点——因为我们知道微调和 RLHF 很容易导致模型偏见。Self-RAG 似乎通过添加特殊令牌来解决这个问题，并保持文本生成不变。
不过，处理固定上下文长度仍有改进的空间。这可能通过在 Self-RAG 中添加一个总结组件来实现。实际上，之前已有一些相关的工作（参见：通过压缩和选择性增强改进检索增强型语言模型）。另一个令人兴奋的方向是最近 OpenAI 发布的增加上下文窗口长度——GPT-4 128k 上下文窗口更新。然而，正如在 论坛 中提到的，这个上下文窗口表示输入长度，而输出限制仍然是 4k 令牌。
RAG 代表了工业界将 LLM 融入其数据以产生实际业务影响的最激动人心的方式之一。然而，对于语言模型的 RAG 特定调优还不多。我对这个领域未来的改进感到兴奋。
推断代码在这个 GitHub 仓库中：
[## GitHub - skandavivek/self-RAG: 关于 Self-RAG 的教程。
关于 Self-RAG 的教程。通过在 GitHub 上创建一个账户来为 skandavivek/self-RAG 的开发做出贡献。
如果你喜欢这篇文章，请关注我——我撰写有关生成式 AI 在实际应用中的内容，以及更广泛的数据与社会之间的交集。
欢迎在 LinkedIn上联系我！
如何让简单的数值积分在方程建模任务中让你的生活更轻松
仿真和数值建模
以迈克利斯-门腾方程作为酶催化的例子，并包括一个可运行的网络应用程序及其代码
LucianoSphere (Luciano Abriata, PhD)
·发表于Towards Data Science ·阅读时间 10 分钟·2023 年 10 月 5 日
–
图像由作者从自己的和免费软件的截图中组成。
在自然科学、工程和经济学的数学中，数值积分的艺术作为描述动态系统随时间变化的行为的强大工具。该技术在这些领域中的许多问题上都有大量具体应用，尤其是在处理微分方程时。
从技术上讲，主要挑战出现在当微分方程将一个或多个变量及其导数混合在一起，使得解析积分最终纠缠了变量的不同出现方式，从而无法将其组合在一起并作为其他变量的函数进行求解。在这种情况下，数值积分，而不是常规积分，将会派上用场。
尽管有整套研究和方法可以高效地实现微分方程的数值积分，例如龙格-库塔方法，但我们可以回归基础，以一种非常简单的“类似欧拉”的方式直接解决问题，即通过使用导数的瞬时值并对小变化进行时间积分。这正是我将在这里展示的示例应用程序的情况：对描述酶催化的迈克利斯-门腾方程进行数值积分。全部包括完整的数学、一个可工作的网络应用程序和所有简单的 JavaScript 代码。
本文/教程的风格类似于我最近分享的另一篇文章，该文章介绍了看似复杂的过程实际上是相当简单且极具帮助的：
掌握数据分析和模型拟合中的不确定性，配合实际代码和示例
[towardsdatascience.com
进入微分方程及其简单、直接的数值积分
微分方程常常在建模动态系统时出现，这些系统中的变量彼此变化且随时间变化。这些方程通常涉及时间的导数，代表变化率。然而，它们固有的复杂性可能使得提取直观见解或直接计算特定结果变得困难，例如预测系统何时会达到特定状态。
目前，微分方程的主要问题在于它们未能提供时间和系统变量之间的直接关系。它们描述了变量如何持续变化，但不容易提供在特定时间点达到特定条件或浓度的见解。这一限制可能会妨碍我们进行精确预测和理解所考虑系统的动态。
解决这个问题的显而易见的方法是积分包含导数的项。但在科学领域，当你积分一个已经出现在方程中的变量时，它通常会分离成不同的项，无法合并在一起求解变量。你会在我这里讨论的米氏-门腾酶促反应方程的例子中看到这一点。
另一种方法是使用一些数值方法，其中所谓的 Runge-Kutta 方法特别设计用来解决这些问题。但还有一种更简单的数值方法，它直接且明确地与导数的定义相结合：一个变量因另一个变量的小变化而发生的小变化，例如位置随时间的变化，浓度随时间的变化等。
通过离散化时间和近似变量的连续变化，非常简单的数值积分方法允许我们在小的、可管理的步骤中计算系统如何随时间演变。这些方法使我们能够在离散时间点近似动态系统的行为，从而能够预测何时会满足特定条件或分析系统在特定时刻的行为。
在这里，我将向你展示如何通过数值积分轻松处理像迈克利斯-门腾方程这样的复杂方程，使其在酶学中作为基础模型变得可行和有用。这个方程描述了酶催化反应的动力学，为展示数值积分的优雅和多功能性提供了一个理想的案例研究，它本质上将复杂的动态过程分解为可管理的步骤。此外，该方程的解析积分形式面临的挑战在于其最相关的变量不能互相求解，因此使得数值积分方法更具相关性。
设定示例问题
迈克利斯-门腾方程是一个基础的酶学模型，描述了涉及底物（S）和酶（E）的酶促反应的速率，以生成产物。这个方程在理解酶催化反应的动力学中起着关键作用。然而，由于没有办法将其项移动得到底物浓度作为时间的解析函数，因此不可能得到底物浓度作为时间的显式解。
但正如我将向你展示的，我们可以通过数值积分方程来模拟在给定初始底物和酶浓度以及一组动力学参数下的底物浓度随时间的变化。我还将把这些应用到一个基于网络的工具中，你可以非常轻松地在浏览器中运行你自己的底物浓度模拟，并且你可以查看和复制代码以便进行临时修改。
迈克利斯-门腾方程
迈克利斯-门腾方程表示为：
其中：
V是给定酶和底物浓度下的反应速度，单位为 M/s（摩尔每秒），[E]是酶浓度，通常以纳摩尔到微摩尔表示，但为了统一性在此以摩尔单位表示，[S]是底物浓度，也以 M 表示，Km是描述催化周期中底物结合部分的迈克利斯常数（此处也以 M 表示），而kcat是周转数或催化常数（单位为 1/s）。
该方程描述了反应速度V（即底物浓度[S]对时间t的导数）如何依赖于酶和底物的浓度，给定一组催化参数kcat和Km。
隐式变量，此处为底物浓度
虽然迈克利斯-门腾方程提供了对酶动力学的宝贵见解，但它无法将底物浓度解为时间的显式函数（这仅在对V作为[S]关于t的导数进行积分后才出现）。其核心限制来自方程的非线性和微分性质。
然而，我们可以通过数值方法克服这个限制，特别是通过数值积分反应速度，这可以被看作是当 deltaTime 非常小时的 deltaConcentration/deltaTime，即 d[S]/dt。
数值积分对于许多问题是一个非常强大的数学技术，它通过将动态系统分解为小的时间步长，从而近似系统随时间的行为，在这些小时间步长中会发生一些事情，变量也会变化。
让我们看看它是如何工作的：
-
首先，我们使用已知值初始化系统：酶浓度（[E]）、初始底物浓度（[S]0）、催化常数（kcat）和米氏常数（Km）。此外，我们定义一个时间步长（dt）用于底物浓度变化率的积分，d[S]/dt = V。经过一个时间步长后，我们将得到时间 t = tprevious + dt 的浓度[S]，而追踪底物浓度的总时间长度由迭代次数给出，我们也必须设定这个值。
-
在每次迭代中，我们计算由于酶促反应引起的[S]的瞬时变化率，这由上述的米氏-门腾方程计算出的反应速度 V 表示，代入剩余底物浓度[S]。
-
这个瞬时速度描述了反应在特定时刻的进行速度。我们可以认为在反应的那个点，在一个无穷小的时间内，会消耗掉相应量的底物。由于计算机中无法使用无穷小的数值，我们将其近似为一个非常小的数值（我们可能需要调整这个数值，稍后会有相关说明）。然后我们将步骤 2 中计算的速度乘以我们的 deltat 或时间步长，以确定在这个短时间间隔内底物浓度的变化。这表示在这个微小时间步长内消耗（或生成）的底物量。
-
然后我们通过从当前浓度[S]中减去上一步计算的变化量（V x dt）来更新底物浓度。这一步反映了反应在短时间间隔内的实际进展，并给出了当前时间的新[S]。
-
我们重复这个过程，进行指定次数的迭代，每次根据更新的底物浓度重新计算速度。
-
一旦在增加的时间t下计算了一系列[S]，我们就得到了我们想要的轮廓。
实现准确的积分
时间步长（dt）的选择在数值积分中至关重要。较小的时间步长允许我们捕捉反应动态的更细微的细节，但需要更多的计算资源。相反，较大的时间步长加快了计算速度，但可能导致结果不准确，甚至可能破坏整个过程——在这个例子中，可能会导致底物浓度出现不现实的大变化，从而引入显著的偏差，甚至可能产生比当前底物浓度还大的变化，导致底物浓度变为负值。
因此，在实际操作中，用户应微调时间步长，以在计算效率和准确性之间取得平衡。根据我的经验，最好尝试找到一个小的数值，使曲线形状合理地接近预期（见下文示例），然后慢慢减少或增加dt，同时检查曲线的形状是否受到影响。如果没有，那么你的数值积分就进行得很好。
使用在线网页应用的示例模拟
这种计算瞬时速度、更新底物浓度并重复多个小时间步长的迭代过程，使我们能够模拟反应的动态行为，即使在无法获得明确的数学解时也是如此。而且，制作一个这样的程序非常简单，就像你可以在这里尝试的网页应用一样：
[## 底物与时间的关系，来自米氏-门腾方程
编辑描述
lucianoabriata.altervista.org](https://lucianoabriata.altervista.org/jsinscience/michaelis-menten/test1.html?source=post_page-----7b5a90b97056--------------------------------)
这是一个模拟示例：
作者截图，来自他自己的网页应用。
这样的工具使用户能够非常轻松地探索酶催化反应的动态行为。特别是，了解在不同参数下底物轮廓如何变化有助于设计实验；例如，向反应中添加多少酶、预期等待多久才能消耗掉给定的底物分数等等。这些都是实验室研究酶性质时的相关问题，从基础化学生物学到制药公司。所有这些都由数学和计算机代码辅助！
一些代码——虽然实际上，它非常简单
我觉得很美妙的是，这样一个复杂的数学问题可以通过像这样的策略来解决，虽然使用了蛮力，但可以如此简单地编码。
虽然我邀请你通过在大多数网页浏览器中按 Ctrl-U 来查看我网页应用的完整代码，但这里是它的核心部分：
var lastS = S0 //Set to the starting concentration of substrate
var time = 0
for (var i = 0; i <= iterations; i++) {
var Vinst = kcat * E * lastS / (Km + lastS) //Instantaneous velocity at the current substrate concentration
var deltaS = Vinst * dt //Integration by multiplying the velocity by a sufficiently small timestep
lastS = lastS — deltaS //Update substrate concentration
time = time + dt //Update total time
data.addRow([time, lastS, S0 — lastS]); //Data is then plotted using Google's chart library in my web app
}
与更复杂的数值解微分方程的方法比较
还有一些专门设计来解决这些及其他相关问题的更复杂的方法，其中主要的是 Runge-Kutta 系列方法。虽然我可以在未来讨论 Runge-Kutta 方法，但目前让我与您分享这篇专门的维基百科文章以及与我在这里提出的方法的比较：
在数值分析中，Runge-Kutta 方法（RUUNG-ə-KUUT-tah）是一系列隐式和显式的迭代…
我上面介绍的方法并没有明确实现特定的 Runge-Kutta 方法，但它确实类似于解决常微分方程的数值积分方案。我们看到的数值方法简单得多，并且具有这样的优势：它在概念上与导数的定义直接相关，这里速度作为“在短时间内的微小变化”，并使用基本的算术运算。
我们使用的简单方法的一个重要缺点是它对时间步长非常敏感，不仅是当该参数过长时，而且当涉及的速度过高时也是如此。在 Web 应用中，这种限制通过可以更改时间步长的选项得到了一定程度的弥补，你可以通过重复模拟来轻松检查模拟是否稳健。
相比于设计为随着时间步长趋近于零而收敛于真实解的 Runge-Kutta 方法，这里使用的方法可能需要如此小的时间步长以达到可接受的精度，以至于计算可能变得计算上过于昂贵（因为更小的时间步长需要更多的迭代来实现长时间的模拟）。
要了解 Runge-Kutta 及其他不同方法在迈克利斯-门腾方程问题上的应用，请查看这些论文：
使用经典的亨利-迈克利斯-门腾（HMM）模型（或简称迈克利斯-门腾模型）来研究底物…
Springer [## 迈克利斯-门腾方程的分解法解 - PubMed
我们提出了一种低阶递归解法来解决迈克利斯-门腾方程，使用了分解法。这…
收获
我希望我已经说服你，数值积分是一种简单而强大的数学技术，它能够解决像这样无法获得积分变量显式解的问题。正如你所看到的，这个过程将复杂的动态分解为可管理的步骤，采用类似欧拉的方法，这种方法根植于导数的基本概念，因此从底层数学和编程的角度来看非常简单。
www.lucianoabriata.com 我写作和拍摄关于我广泛兴趣领域的所有内容：自然、科学、技术、编程等。 订阅以获取我的新故事 通过电子邮件。* 要* 咨询小型工作 请查看我的 服务页面。你可以 在这里联系我。
辛普森悖论如何误导统计数据
原文：
towardsdatascience.com/how-simpsons-paradox-can-mislead-statistics-f63d1c6a8e15
以及为什么机器学习的解释并不总是可靠
·发表于 Towards Data Science ·阅读时间 12 分钟·2023 年 1 月 18 日
–
图片由 Jason Leung 提供，来源于 Unsplash
如果有我真正讨厌听到的话，那就是经典的权威论证“统计数据显示插入一个事实”。
随着统计工具和机器学习的民主化，将一些数据转化为洞察比以往任何时候都更加容易。然而，每个人都应该注意不要陷入一些非常违反直觉的陷阱，这些陷阱可能导致研究偏差，尤其是在研究速度过快或没有专家视角的情况下。
在这篇文章中，我们将深入探讨一个经典的统计陷阱，即辛普森悖论，它可能导致由于存在隐藏的相关因素而对模型中特征的解释产生误解。
我设计这篇文章是为了那些具备基本统计和机器学习知识的人，特别是那些对这个悖论不了解的数据科学家和分析师。不过，第二部分中展示的案例可能会引起任何有意发现数据被操控以表达相反观点的人的兴趣。
这个悖论是什么？
你正在对一个数据集进行分析。你有一堆特征 X1, … Xn，你正在使用这些特征来预测目标变量“y”。
辛普森悖论出现的条件是：
-
你至少缺少一个变量（我们称之为 Xs），这个变量解释了你的目标“y”方差的一部分。
-
这个隐藏变量 Xs 和你的一个特征 Xb 之间存在直接的因果关系。
在这种情况下，忽略 Xs 可能会直接导致对 Xb 对目标变量 y 影响的误解。
悖论的插图
线性回归的情况
为了理解这个悖论的背后，我们来研究一个简单的线性回归案例，其中目标“y”与两个变量“X1”和“X2”相关。我们还添加了一个误差项，该误差项独立且分布相同：
在正常情况下，当变量 X1 和 X2 不相关时，忽略 X2 并对 X1 进行回归将会得到：
具有：
让我们用一堆模拟数据来验证这一点，取 a = 1 和 b = 1：
a = 1
b = 1
#Let's draw some random points
X1 = np.random.normal(0, 1, 5000)
X2 = np.random.normal(0, 1, 5000)
eps = np.random.normal(0,1, 5000)
#Create our target value
y = a*X1 + b*X2 + eps
#Get the coefficients of the linear regression of y on X1
lr = LinearRegression().fit(X1.reshape(-1,1),y)
theta_1,theta_0 = lr.coef_[0], lr.intercept_
我们确认，当 X2 和 X1 之间没有相关性时：
-
y 对 X1 的线性回归的斜率等于“a”（第二张图）
-
回归的误差 epsilon1 等于 bX2 + epsilon（第三张图）。
当 X1 和 X2 不是独立的
我们现在假设 X1 和 X2 之间存在直接相关性：
在这种情况下，我们不能直接将 X2 再注入误差项，因为误差将不再是独立的（由于 X1 和 X2 之间的相关性）
由于 X1 中的项，误差不再是独立的
要恢复一个独立的误差，我们需要重新排列 X1 的方程：
这次，不像之前，系数不再等于“a”，而是包含一个考虑 X1 和 X2 之间相关性的项（a + bc）。
这就是辛普森悖论所在：根据 b（X2 对 y 的解释方差）和 c（X2 对 X1 的解释方差）的值，我们现在可以得到 y 和 X1 之间的相关性，这种相关性可能更高，更低，甚至… 相反！
注意，如果 b = 0 或 c = 0，则悖论消失，确认了在这个简化案例中之前暴露的两个条件。
为了直观地展示这个效果，让我们修改我们的模拟：
a = 1
b = -1
c = 2
X1 = np.random.normal(0, 1, 5000)
#Define X2 as correlated to X1
eps_2 = np.random.normal(0,1,5000)
X2 = c * X1 + eps_2
#Define the target as before
eps = np.random.normal(0,1, 5000)
y = a*X1 + b*X2 + eps
lr = LinearRegression().fit(X1.reshape(-1,1),y)
theta_1,theta_0 = lr.coef_[0], lr.intercept_
当 b = -1 和 c = 2 时，我们应该在对 X1 进行回归时发现一个负系数，而不是一个正系数…
… 这正是我们发现的（第二张图）。注意我们还找到了线性回归误差为 epsilon + b epsilon_2（第三张图），这是我们的小计算所预期的。
就是这样，我们揭示了我们线性系统中的辛普森悖论。通过忽略变量 X2 的影响，X2 与 X1 和 y 都相关，我们使自己暴露于对 y 对 X1 回归系数的错误解释中。
注意，根据我们的简单例子，悖论可能发生在 b 和 c 的广泛值范围内。例如，b = 2 和 c = -1 也会创建一个辛普森悖论，其中 X1 和 X2 负相关，而 X2 对 y 具有强影响。
这当然是一个非常简单的例子，但有助于发展对这一特定现象的直觉。
一些著名的辛普森悖论
这看起来像一个有趣的理论问题，但这种悖论在现实情况中有很多表现。让我们深入探讨两个真实案例，以查看统计数据如何可能误导那些不完全了解这些细微差别的人，并看看它如何与我们在上一部分中强调的内容相关。
我参考了David Louapre的优秀文章（法文），这激发了我深入分析悖论。
吸烟与预期寿命
在“忽略一个协变量：辛普森悖论的一个例子”[1]的第一部分中，David R. Appleton、Joyce M. French 和 Mark PJ Vanderpump 研究了吸烟和非吸烟女性（X0）的生存率（y）。
为了方便可视化这种情况中的悖论，我们将模拟一个类似的数据集：
#Proba to use for non-smokers
proba_not_smoke = {
"18-25":0.52,
"26-35":0.55,
"36-44":0.51,
"45-55":0.46,
"56-67":0.53,
"68-75":0.8,
"76+":0.85
}
#Proba to use for people dying naturally
proba_die_naturally = {
"18-25":0.01,
"26-35":0.01,
"36-44":0.05,
"45-55":0.13,
"56-67":0.3,
"68-75":0.6,
"76+":0.95
}
#Proba to use for people dying smoking
proba_die_smoking = {
"18-25":0.001,
"26-35":0.004,
"36-44":0.008,
"45-55":0.01,
"56-67":0.01,
"68-75":0.01,
"76+":0.01
}
def get_group(age):
'''Simple function to transform age into a category'''
if age<=25:
return "18-25"
if age<=35:
return "26-35"
if age<=44:
return "36-44"
if age<=55:
return "45-55"
if age<=67:
return "56-67"
if age<=75:
return "68-75"
else:
return "76+"
POPULATION = 5000
#Create a random population
population_ages = np.random.randint(18,90, POPULATION)
dataset = []
for age in population_ages:
group = get_group(age)
#draw the probas
p_not_smoke = proba_not_smoke[group]
p_die_naturally = proba_die_naturally[group]
p_die_smoking = proba_die_smoking[group]
#calculate the condition of the person based on the proba drawns
smoker = np.random.random()>p_not_smoke
died_naturally = np.random.random()<p_die_naturally
died_smoking = np.random.random()<p_die_naturally
died = died_naturally | ((died_smoking) & (smoker))
dataset.append({"smoker":smoker*1,"age_group":group,"lived_after_20_years":1-died})
df = pd.DataFrame(dataset)
数据集样本（合成的）
下图显示了我们模拟数据集中的预期寿命：
在模拟数据集中，按我们关注的特征（吸烟）拆分数据集后的预期寿命
如我们所见，在对总体人群进行简单分析时，初步结果似乎表明吸烟者的预期寿命较高。
这个非常奇怪的结果来源于一个混杂变量，即所研究人群的年龄（我们的变量 X2，参见第一部分的例子）在图表中被省略了。
现在让我们看看在分析中添加年龄变量后结果如何。
在模拟数据集中，按我们关注的特征（吸烟）和年龄拆分数据集后的预期寿命
这次我们可以看到，每个组的 20 年后的预期寿命在每个年龄组中对非吸烟者更好。
当我们从不同的角度查看数据时，我们可以发现年龄（X2）与吸烟者/非吸烟者类别（X1）之间的相关性。
在我们模拟的数据集中，按年龄组分配吸烟者的情况
我们回到第一部分中暴露的条件集：
-
我们的变量年龄“X2”与目标（y）“20 年后的预期寿命”（b << 0）呈强负相关，因为老年人不太可能比年轻人多活 20 年。
-
同时，年龄（X2）与吸烟者（X1）之间呈负相关（c < 0）
-
如果忽略年龄 X2 的话，结果会给人一种吸烟能提高预期寿命的印象（a + bc > 0），而实际上它会降低预期寿命（a < 0）。
性别和学校录取
另一个非常著名的例子是《研究生招生中的性别偏见：来自伯克利的数据》[2]中突出的辛普森悖论。
原始数据集展示了 1973 年伯克利按部门和性别的招生情况，我们还将准备类似的合成变体以作说明。
#Probabilities of admission in the major
admissions_proba = {
"Major_1":0.74,
"Major_2":0.64,
"Major_3":0.35,
"Major_4":0.34,
"Major_5":0.24,
"Major_6":0.06
}
#Probabilities of selection a major as a man
gender_proba = {
"Major_1":0.08,
"Major_2":0.04,
"Major_3":0.65,
"Major_4":0.47,
"Major_5":0.67,
"Major_6":0.47,
}
def get_major(p):
"""This function is used to select a major randomly based on some
probabilities"""
if p<0.25:
return "Major_1"
if p<0.37:
return "Major_2"
if p<0.54:
return "Major_3"
if p<0.70:
return "Major_4"
if p<0.82:
return "Major_5"
else:
return "Major_6"
POPULATION = 12000
dataset = []
for i in range(POPULATION):
p_major = np.random.random()
major = get_major(p_major)
gender = ['M','F'][np.random.random()<gender_proba[major]]
admission = ['Accepted','Rejected'][np.random.random()>admissions_proba[major]]
dataset.append({"Major":major,"Gender":gender,"Admission":admission})
df = pd.DataFrame(dataset)
合成数据用于说明研究生招生中的性别偏见悖论
当我们简单查看按性别（X1）的统计数据时，我们发现男性性别的招生百分比（y）存在统计学上显著的偏见。
按性别分割的简单招生率，在合成数据上
另一方面，添加主修（X2）进行分析会在大多数部门中颠倒性别偏见的结论，揭示出悖论。
添加主修暴露了悖论，合成数据上的示例
与之前的例子一样，这个悖论可以通过主修（X2）与目标 y（选择性）和原始特征 X1（性别）的关联来解释。
下图显示，在数据集中，女性倾向于申请非常有选择性的专业，而男性则申请选择性较低的专业。
主修（X2）与我们的目标以及原始特征之间的隐藏关联
高级机器学习工具中的悖论
直到现在，我们讨论了相对简单数据集中的辛普森悖论。但也有可能生成更复杂的例子，当高级机器学习和特征解释算法可能会产生误导，如果不包括所有相关数据的话。
在这一部分，我模拟了一个由 5 个特征组成的“房价数据集”：
-
资产表面
-
车库面积
-
花园面积
-
房屋状况
-
区域豪华度（可选参数）
我们将根据这 5 个特征预测资产价格。
作为一个心机深沉的人，我设计了稍微不同的区域豪华度特征：
我构建了这个特征，使得区域的豪华度对我们的目标（销售价格）有不可忽视的影响，同时对车库面积也有影响，这将成为悖论的受害特征（前面例子的“X1”）：例如，我们可以想象在豪华区域，可能没有足够的空间容纳大车库，从而解释负相关性。
这当然是完全虚构的，所以不要太当真，但这将帮助我阐明我的观点。
N = 3000
n_clusters = 10
inc = N//n_clusters
n_samples = inc*n_clusters
#Create our main features
asset_surface = np.random.normal(90, 25, n_samples)
garden_surface = np.random.normal(500, 100, n_samples)
house_condition = np.random.randint(0, 5, n_samples)/10+1
#Build the two last features so they have a negative correlation
garage_surface = np.array([])
fanciness_area = np.array([])
for i in range(n_clusters):
gs = np.random.normal(i, 1, inc)
fa = np.random.normal(-i, 1, inc)
garage_surface = np.hstack([garage_surface, gs])
fanciness_area = np.hstack([fanciness_area, fa])
#Just make it a bit more non-linear...
fanciness_area = np.exp((fanciness_area + 12)/10)
df = pd.DataFrame([asset_surface, garden_surface, house_condition, garage_surface, fanciness_area])
#We build our target so that the formula is non-linear and stay in the range of a house price...
sale_price = (asset_surface*(fanciness_area**0.5)*2 + garden_surface/5 + garage_surface/2)*house_condition*100+50000+np.random.normal(0, 20000, n_samples)
我将跳过特征工程和清理，直接进入结论。我使用了经过微调的xgboost来预测 y，使用数据集中的所有特征（包括和不包括我们的“区域豪华度”特征），并利用shap库检查特征解释的结果，就像我在经典机器学习项目中所做的一样。
注意：对于那些不熟悉这个库的人， shap 是一个出色的工具，用于执行特征解释，特别是对于基于树的算法特别高效。
一些解释图表的关键，对于那些不熟悉 shap 的人：
-
特征按重要性从上到下排序，列表顶部的特征对目标有较大影响，而底部的特征影响较小。
-
每个点代表一个样本。红色表示特征在给定样本中值较高，而蓝色表示值较低。
-
一个点位于左侧越多，该特征对该样本的负面影响越大。越往右，该特征的正面影响越大。
Shap 值根据区域时尚特征（X2）的存在，对车库表面（X1）给出了不同的解释。
在我们的场景中，我们可以清楚地观察到悖论：如果我们忘记包括“时尚区域”特征，即使分数没有那么不同，它也会影响我们对结果的解释。它将带来两个重大后果：
-
高估车库表面的影响：在没有区域时尚的情况下排名更高
-
改变我们对车库表面特征的解释：如果我们不考虑区域时尚，shap 模型告诉我们，高车库表面倾向于对销售价格产生负面影响，而解释是高车库表面如果包括该区域时尚将提升销售价格。
在这种情况下，我们可以简单地理解存在一个问题，并且一个大的车库理论上应该增加资产的价值。
然而，在其他情况下，当特征数量较多且更难以解释时，你可能会在未意识到的情况下陷入悖论，并对你的结果做出错误的解释。
结论
在这篇文章中，我尝试向你解释著名的辛普森悖论背后的直觉，并突显一些可能隐藏在统计分析背后的危险。
我们特别看到，拥有最好的机器学习模型和最先进的解释算法是不足以避免这种情况的，如果你对所研究的数据没有全面的了解。
我个人没有避免这种悖论的解决方案，但作为数据科学家和数据分析师，我们需要记住“天下没有免费的午餐”，在可能的情况下寻求领域专家的帮助，并始终质疑我们的结果。
这篇文章对我特别重要，因为当我第一次了解到这个悖论时，我变得非常渴望学习更多的统计学知识，这也使我走到了今天的道路上。
我希望你阅读它时和我写作时一样享受。
[1] Appleton, David R., Joyce M. French, 和 Mark PJ Vanderpump. « 忽略协变量：辛普森悖论的一个例子。 »《美国统计学家》50.4 (1996): 340–341。
[2] Bickel, Peter J., Eugene A. Hammel 和 J. William O’Connell. « 研究生招生中的性别偏见：来自伯克利的数据。 » 科学 187.4175 (1975): 398–404.
SQL 执行顺序在不同数据库中的差异
为什么在 SQL Server 中不能按顺序位置进行 GROUP BY，但其他数据库可以
·
关注 发表在 Towards Data Science ·4 min 阅读·2023 年 12 月 7 日
–
Transact-SQL 与 MySQL 执行顺序（作者提供的图片）
在与 MySQL 和 PostgreSQL 等开源数据库定期工作之后，我最近有机会参与一个 SQL Server 项目，并发现了 SQL 环境中的一个微妙但重要的区别。我观察到在 SQL Server 中，我无法按顺序位置进行 GROUP BY（GROUP BY 1, 2, 3…），而这是我在其他数据库中经常使用的功能，特别是用于快速测试。
这一发现让我探讨了这两个数据库系统的几个细微差别，特别是SQL 执行顺序，这将是本文的重点。
为什么这很重要？在使用数据库系统时，了解这些细微差别可以极大地影响你的工作流程，提高你的生产力。它可以节省你大量的故障排除时间。此外，通过 理解各种数据库的 SQL 执行顺序，你可以根据你所使用的系统编写更优化的 SQL 查询。
在本文中，我们将探讨一个主要的用例——GROUP BY——并调查其原因。然而，这些见解也可以应用于 HAVING、WHERE 或任何其他 SQL 命令子句。
让我们开始吧
让我们看看下面这个查询示例。虽然在 MySQL 中它可以正常工作，但在 SQL Server 中这将不起作用：
SELECT
DATEPART(year, day) AS order_date,
SUM(cost) as cost
FROM clean
GROUP BY 1;
如果你运行这个，你可能会得到类似这样的错误：
每个 GROUP BY 表达式必须包含至少一个不是外部引用的列。
然而，在将 GROUP BY 序号引用替换为显式表达式后，这个修订的查询可以正常工作。你还会注意到，在 ORDER BY 子句中可以引用序号位置，这让我觉得很奇怪：
SELECT
datepart(year, day),
sum(cost) as cost
from clean
GROUP BY datepart(year, day)
ORDER BY 1;
在 SQL Server 中，我很快了解到我必须在 GROUP BY 子句中使用显式的列名或表达式。这被认为是一种最佳实践，因为它使代码更易于理解。然而，我对为什么不同数据库之间会有这种行为差异感到好奇。此外，我发现 SQL Server 中的ORDER BY子句与序号位置一起使用，这进一步激起了我的好奇心。
探索 SELECT 语句执行顺序
为了了解情况，我们来看一下SELECT语句在 SQL Server 和其他数据库中的执行/处理顺序。需要注意的是，在 SQL 数据库中，每部分查询都是按顺序执行的，这一顺序与书写的顺序不同。
例如，在 SQL Server 中，我们可以从下面的图片和Microsoft 文档中看到，FROM 子句是第一个被评估的命令。此外，SELECT 子句在 GROUP BY 子句之后运行。这就是为什么在第一个示例中我们无法在 GROUP BY 子句中引用列的位置或其别名的原因！
然而，我们可以在 ORDER BY 子句中自由引用序号位置和/或别名，因为这是在 SELECT 子句之后被评估的。SELECT 子句告诉数据库哪些列将被返回，因此此时位置是已知的。很酷，对吧？
SQL Server 执行顺序
SQL Server SELECT 语句处理顺序（图片由作者提供）
MySQL
然而，在 MySQL 中，我发现很难找到明确的文档说明 SQL 查询的执行顺序。执行顺序似乎依赖于查询的内容以及查询优化器定义的最佳路径。
但从我们从 MySQL 文档中看到的情况来看，这些线索显示了执行顺序可能的方式，即在评估 GROUP BY 子句之前，SELECT 子句先被评估：
对于 GROUP BY 或 HAVING 子句，它在搜索select_expr值之前先搜索 FROM 子句。（对于 GROUP BY 和 HAVING，这与 MySQL 5.0 之前的行为不同，后者使用与 ORDER BY 相同的规则。）
GoogleSQL
如果我们再看看 GoogleSQL（前称标准 SQL）的文档，这是 Google BigQuery 中使用的语法，你会看到与 SQL Server 查询执行方式不同的类似偏差：
GROUP BY 和 ORDER BY 还可以引用第三组：整数文字，它们指的是 SELECT 列表中的项。整数 1 指的是 SELECT 列表中的第一项，2 指的是第二项，依此类推。
如你所见，这种行为在 SQL Server 中不受支持。Google 文档还提到，GROUP BY、ORDER BY 和 HAVING 可以引用来自 SELECT 列表的别名。
综上所述，我们可以有很大把握认为其他数据库的执行顺序与下图类似：
MySQL、PostgreSQL 和 BigQuery 可能的执行顺序
MySQL SELECT 语句执行顺序（作者提供的图片）
结论
这篇文章简短地介绍了 MySQL、GoogleSQL 和其他数据库 SQL 语法的执行顺序，与 SQL Server 的区别，基于观察到的行为和文档。SQL Server 强调在 GROUP BY 子句中显式声明以提高代码清晰度，而 MySQL 的执行顺序明确将 SELECT 子句评估在 GROUP BY 子句之前，允许我们引用其中的序号位置。
欢迎分享您对这个主题的想法，并期待在下一篇中与您见面。
您可以成为 Medium 会员来支持我，享受更多类似的故事。
参考资料
-
ROWNUM — Oracle
-
PostgreSQL - EXPLAIN
变量之间的关联强度有多大？
原文：
towardsdatascience.com/how-strongly-associated-are-your-variables-80493127b3a2
使用 Cramer’s V 检验来检查两个分类变量的关联强度
·发表于 Towards Data Science ·阅读时间 7 分钟·2023 年 2 月 28 日
–
图片由 Susan Holt Simpson 提供，刊登于 Unsplash
介绍
特征选择是任何数据科学项目中的一个重要步骤。我可以想象如果你对这个领域不陌生，你已经听到过无数次了，如果你是新手，我相信你也听说过，但我还是要再说一遍：如果你给模型喂垃圾，你将得到垃圾作为结果。
好了，现在我们把这件事放在一边，接下来继续。选择最佳特征有几种好的方法，例如运行随机森林模型，然后检查 feature_importances_ 属性，使用 sklearn 的 SelectKBest，分别进行统计测试，以及其他技术。
sklearn 提供的这些自动化测试非常实用，是我们快速进行特征选择的优秀选择。它们是使用统计测试（如 F 检验、相关性、卡方检验）的自动化方法，可以快速执行假设检验，以根据结果选择变量。
当我们谈论分类变量时，例如，如果我们运行 SelectKBest，我们将不得不使用评分函数 chi2 来找出 p 值是否低于变量依赖的统计显著性阈值。
Ho = 变量独立
Ha = 变量不独立
常用的显著性水平是 0.05。
然而，返回的值仅包括 p 值和检验统计量。这意味着该工具旨在仅提供快速结果，例如 p 值低于阈值，确认我们有证据拒绝零假设，并指出两个分类变量之间存在依赖关系。但它不会告诉你这种关联有多强。
一个简单的解决方案是 执行 Cramer’s V 检验，将在本帖中介绍。
在继续之前，让我介绍一下本帖中示例使用的数据集。这是 diamonds 数据集，一个来自 Seaborn 包的开放样本数据。
import seaborn as sns
# Load the dataset
df = sns.load_dataset('diamonds')
来自 seaborn 的钻石数据集。图片由作者提供。
该数据集包含钻石的切割、颜色、大小、克拉和价格的观察值。我们的目的是检查分类变量之间的关联 **cut, color, clarity** 以及 **price**。
可以自由导入包进行代码编写。
import pandas as pd
import scipy.stats as scs
from sklearn.feature_selection import SelectKBest
from sklearn.feature_selection import chi2
使用 Select K Best
好的，一旦我们介绍了数据集和 SelectKBest 工具要执行的检验，让我们看看它是如何工作的以及结果如何。
首先，由于我们将检查两个分类变量之间的关联，让我们将 price 变量变成分类变量。为此，我将钻石的价格分成几个区间：
-
便宜：从 0 到低于平均值的 20%
-
平均值：从平均值 - 20% 到平均值
-
高价格：从平均值到平均值 + 20%
-
昂贵：从平均值 + 20% 到最大值
# Create a price bin variable
df['price_bins'] = pd.cut(
df['price'],
bins= [0, df.price.mean()*0.8, df.price.mean(), df.price.mean()*1.2, np.inf],
labels= ['cheaper', 'on_average', 'high_price', 'expensive']
)
这就是它的样子。
添加了 price_bins 变量。图片由作者提供。
接下来，我们将数据分成 X（解释变量）和 y（被解释变量）。
# Split X and y
X = df.drop(['price', 'x', 'y', 'z', 'depth',
'table', 'carat', 'price_bins'], axis=1)
y= df.price_bins
在此过程中，我们选择分类变量并对其值进行编码。
# Select categorial variables
categorical_vars = X.select_dtypes(include='category').columns.to_list()
# Encode the categorial variables
X[categorical_vars] = X[categorical_vars].apply(lambda x: x.cat.codes)
有了这些，我们就准备好拟合数据并提取结果了。
# Instance of SelectKBest
fsel= SelectKBest(score_func=chi2, k=3)
# Fit
fsel.fit(X, y)
# Show a dataframe of the results
(
pd.DataFrame({
'variable': X.columns,
'chi2_stat': fsel.scores_,
'p_value': fsel.pvalues_})
.sort_values(by='p_value', ascending=False)
)
结果数据框见下图。
SelectKBest 工具的结果。图片由作者提供。
太好了！现在我们可以看到这 3 个特征在 95% 的统计显著性水平下与价格相关（p 值 < 0.05；拒绝 Ho）。但我们如何确保这些关联是否强烈呢？
让我们进行 Cramer’s V 检验并找出结果。
Cramer’s V 检验
如果我们做个快速研究，会发现：
Cramer’s V 是两个分类变量之间的关联度量，其值介于 0（弱）和 1（强）之间。
要在 Python 中执行此操作，我们需要 Pandas 来创建列联表，scipy 用于运行 Chi² 检验，这将导致 V 值的最终计算。
好的，从我们的数据集中，让我们创建一个 cut 和我们最近创建的 price_bins 变量之间的列联表。
# Creating a contingency table
cont_table = pd.crosstab(index= df['cut'],
columns= df['price_bins'])
接下来是结果列联表。这仅仅是计算每对关联中存在的观察数量。例如，有 14,181 颗钻石的切割为“理想”，价格为“更便宜”。735 颗钻石的切割为“优质”，价格“平均”。
切割 x 价格分箱的列联表。图像由作者提供。
要执行 Chi²检验，我们将使用scipy中的chi2_contingency函数。因此，我们将列联表传递给该函数。它返回很多有趣的信息：（1）Chi²统计量；（2）p 值；（3）自由度；（4）期望值。由于我们只需要 Chi²进行此测试，因此我们取第一个索引chi_stat = X2[0]。
# Chi-square value
X2 = scs.chi2_contingency(cont_table)
chi_stat = X2[0]
# Print X2
X2
[OUT]
(1603.5199669055353,
0.0,
12,
array([[12378.84006303, 1318.4705228 , 1509.44898035, 6344.24043382],
[ 7921.51562848, 843.72080089, 965.93248053, 4059.8310901 ],
[ 6939.87033741, 739.16573971, 846.23277716, 3556.73114572],
[ 2817.9940304 , 300.14460512, 343.6200964 , 1444.24126808],
[ 924.77994067, 98.49833148, 112.76566555, 473.95606229]]))
现在，让我们计算 Cramer 的 V 值。V 的公式是：
V 的公式：X²是卡方统计量；N 是样本大小；k 是行和列中类别数的最小值。
在 Python 中，计算通过以下代码片段完成。
# Size of the sample
N = len(df)
# Minimum dimension
# Minimum between Number of categories in rows-1, # categs columns-1
minimum_dimension = (min(cont_table.shape)-1)
# Calculate Cramer's V
result = np.sqrt((chi_stat/N) / minimum_dimension)
# Print the result
print(result)
[OUT]
0.09954537514956
切割与价格分箱之间的 Cramer’s V 强度值为 0.099，即 9.9%，可以理解为从小到中等效果之间的范围（参见此处的值解释表）对于 3 个自由度（minimum_dimension = 3）。
为了clarity，这是结果：
# Creating a contingency table
cont_table = pd.crosstab(index= df.clarity,
columns= df['price_bins'])
# Chi-square value
X2 = scs.chi2_contingency(cont_table)
chi_stat = X2[0]
# Performing Cramer's V calculation
# Size of the sample
N = len(df)
# Minimum dimension
minimum_dimension = (min(cont_table.shape)-1)
# Calculate Cramer's V
result = np.sqrt((chi_stat/N) / minimum_dimension)
# Print the result
print(result)
[OUT]
0.18476912508901078
为了clarity，数字 0.18 大约是中等强度。对于color，V = 0.115，也在小效果范围内。
在你离开之前
这个主题引起了我的注意，因为确实，我们很多时候只是使用像 sklearn 的SelectKBest这样的自动化解决方案，并且将它们的 p 值视为决定模型变量的唯一标准。
现在你读完了这篇文章，你拥有了另一个统计工具来选择模型的最佳变量。
总结：
-
选择两个分类变量
-
使用
pd.crosstab()创建列联表 -
运行
scs.chi2_contingency(contingency_table)以收集 Chi²统计量。 -
计算 Cramer 的 V：
np.sqrt((chi_stat/N) / minimum_dimension)
代码：
[## Studying/Cramers_V.ipynb at master · gurezende/Studying
你目前不能执行该操作。你在另一个标签或窗口中登录。你在另一个标签或窗口中注销了……
如果你喜欢这个内容，关注我的博客获取更多。在Linkedin上找到我。
参考资料
[## Cramér’s V - 维基百科
来自维基百科，自由百科全书。在统计学中，Cramér’s V（有时称为 Cramér’s phi，记作）…
en.wikipedia.org Statology [## 如何解释 Cramér’s V（带例子） - Statology]
它的范围从 0 到 1，其中：0 表示两个变量之间没有关联。1 表示完全关联…
www.statology.org spss-tutorials.com [## 列联表、卡方检验与 Cramer’s V]
如何轻松检查分类变量之间的关联
towardsdatascience.com SciPy 文档 [## scipy.stats.chi2_contingency - SciPy v1.10.1 手册]
卡方独立性检验用于检验列联表中变量的独立性。此函数计算卡方统计量…
docs.scipy.org Statology [## 如何解释 Cramer’s V（带例子） - Statology]
它的范围从 0 到 1，其中：0 表示两个变量之间没有关联。1 表示完全关联…
www.statology.org spss-tutorials.com [## Cramér’s V - 初学者教程]
Cramér’s V 是一个介于 0 和 1 之间的数字，表示两个分类变量之间的关联强度。如果我们…
不要在没有这 3 步策略的情况下申请数据科学职位
成功的数据科学家如何在今天的劳动危机中获得工作邀请
·发表于数据科学前沿 ·13 分钟阅读·2023 年 12 月 10 日
–
我在 Blind 上进行的调查快照——作者提供的图像
我最近对数据科学家进行了调查，发现一个令人震惊的数字——86%的人盲目地发送求职申请，希望能有好运。希望不是策略，在这样的时期，它可能是获得工作和继续失业之间的区别。
你不是在申请快餐连锁店，所以别再做通用化的求职申请了！数据劳动市场已经变成了一个只有那些脱颖而出的人才能成功的战场。这意味着你不再能随大流。
你需要脱颖而出，为此你需要一个策略。
如果你符合以下条件，这篇文章适合你：
-
你是那些没有策略的自由求职者中的一员，并且已经厌倦了一个接一个的拒绝
-
你想了解成功的数据科学家如何在经济危机中找到工作。
我收集了成功的数据科学家的推荐，并与其他几位成功者交流，以了解他们的获胜求职策略。在这篇文章中，我将他们的智慧和经验汇聚成一个 3 步策略，帮助你度过这些困难时期，并取得成功。
没有策略就是失败的策略。所以让我们带你走向成功的一方！
1. 将期望与现实对齐
图片由Yosep Surahman提供，来源于Unsplash
在求职时，了解哪些因素在你控制范围内，哪些因素超出你的控制范围是至关重要的。这是制定成功策略的关键，因为它帮助你专注于可以优化的方面，以更好地导航你的求职过程。
超出你控制范围的因素及你能做的事
-
职位空缺： 许多职位不会公开招聘。此外，一些职位被填补得很快，通常是由第一个合格的候选人填补。一般来说，简历提交得越晚，越有可能没有被审阅。这意味着时机和运气发挥了重要作用。
该怎么办？尝试尽早申请职位。
-
竞争强度： 高需求的职位吸引成千上万的申请者。这意味着你将与经验更多的竞争者一较高下。这并不意味着他们会得到这个职位，但他们更有可能通过简历筛选。
该怎么办？社交网络！这可以让你的简历迅速送到招聘人员的桌面上。
-
地点限制： 工作地点的接近程度可能是决定性因素，即使是远程职位也是如此。公司可能更倾向于选择靠近他们办公室或在相同时间区的候选人，特别是在一些公司回归办公室（RTO）政策的背景下。所以你已经处于劣势。
该怎么办？在你的简历中明确表示如果你准备好可以搬迁。
你可以控制的因素及如何应对
-
期望： 重要的是针对与你当前技能水平相符的职位。这将节省你申请不符合资格职位的精力。我过去经常这样做，认为不会有害，但这只是浪费时间 + 令人精神疲惫。你还需要了解如何驾驭招聘流程，你会明白为什么。
-
简历： 简历的目的是让你获得面试机会。招聘人员和招聘经理是你的主要受众，因此将简历量身定制以满足他们的需求和期望是至关重要的。
-
社交网络： 单独申请是无济于事的。你必须通过社交网络更好地了解你申请的公司，并尽可能获得推荐。
技术职位的典型招聘流程及如何应对
作者根据TechResume的资源和个人经验绘制的图示
再次强调，如果你的目标是优化你的方法，了解外部情况非常重要。你需要专注于了解典型的招聘流程。你会很快明白，制作一份出色的简历为何可能决定你的成败！
阶段 #1：你的简历首先进入申请追踪系统
大多数大型公司使用系统来管理申请。这些系统不会拒绝申请，这是一个误区；是人类拒绝的。然而，ATS 会从你的简历中提取特定的关键词，寻找与职位相关的关键词。许多 ATS 甚至可以根据这些关键词对你的简历进行评分。
阶段 #2：如果幸运的话，它会被人工审阅
一些公司将这一过程外包给第三方，这些第三方会扫描简历，通常在 5-20 秒内评估你是否合适。
在其他公司，招聘人员通常会在 ATS（申请跟踪系统）中搜索与他们需要填补的职位相关的特定关键词。这就是为什么你必须在简历中使用职位描述中的关键词。你需要制作一个适合 ATS 的简历，以引起招聘人员的注意。
招聘人员随后会花费几秒钟时间审阅你的简历，并决定你是否适合这个职位。在这两种情况下，你的简历必须立即脱颖而出！在这里，通常会有很大的淘汰率，只有几十名候选人能够通过这一筛选。
阶段 #3：招聘人员筛选 在一些公司，招聘人员可能会首先打电话给你，与您进行简短交流，然后将你推进下一阶段。我在育碧公司经历过这种情况，但在 Spotify 没有。
这次面试的目标，通常是一个电话，是为了再次确认你的动机和经验。在这时，你要展示你有多么兴奋！
阶段 #4：技术 + 文化适配筛选 只有有前途的简历才能进入这一阶段。这包括文化适配面试和编码/业务挑战，有时还会有不止一个。
对于技术挑战，候选人通常需要在全职员工面前解决这些挑战。这通常包括 SQL 和 Python 案例，有时还会有数据业务案例。对于其他公司，你可能会遇到一个需要带回家的 ML 挑战。
当涉及到文化适配面试时，招聘经理通常是执行面试并决定你是否适合他们的团队，在软技能、价值观和文化方面。这是你大放异彩的时刻！
最终关卡：恭喜！你收到了工作邀请 现在还不是放松的时候。在这一阶段，你需要谈判，确保你为所有的辛勤工作得到一个好的交易。不要害怕提出你想要的待遇，谈判总是至关重要的，没有人比你更了解自己的价值。只有在这时，你才能给自己一个拥抱，而不是在之前。
2. 一份出色的简历是你的金钥匙
图片来源：作者（Dall.E）
如果你没有收到至少 5-10%的申请回音，而你申请的是你应该有资格的职位，那么更新简历应该是你的首要任务。
我见过很多成功的数据科学家花费数百美元在专业简历服务上，以完善他们的简历。这说明了简历部分的重要性。
以下是制作一份出色简历的逐步指南：
-
停止盲目申请；针对性投递简历。 发送通用简历不起作用——至少在这个被激烈经济压制的丛林中是这样的。挑选 5-6 家公司，并针对它们量身定制你的简历。
-
了解他们的业务并成为他们想聘用的人。 查阅他们的年度报告。了解他们的策略；他们在投资什么？如果你花 5 分钟发送简历，就要期望他们花 5 秒钟浏览并拒绝你。他们为什么要雇佣你？专注于这个问题！
-
量化你的工作影响。 始终记住，你简历中的每一点都应回答你做了什么，如何做到的，最终结果是什么。包括团队人数、使用的数据类型、用户数量、技术栈等。尽量将它们压缩到你描述的每个经历的一两行中。
-
使用主动语言。 始终使用“创建”、“开发”和“领导”等主动动词，而不是被动动词。使用 ChatGPT 改进你的用词。
-
你越资深，你的领域就越重要。 具体描述你的技能至关重要。对于高级职位，比如自然语言处理，仅列出一般的机器学习经验是不够的。你的简历必须提到 NLP、NLU 和 LLMs 才能通过初筛阶段。市场对更有经验的人更有利，但这并不意味着新入行的人就没有机会。你也必须这样做才能脱颖而出。
-
让你的简历接受审阅。 让同行评审你的简历，询问相关领域的朋友意见。这有助于简化内容，使其更易于阅读。再一次，利用 ChatGPT。只需确保提供适当的背景信息，例如：“假装你是一名招聘人员在筛选[插入职位名称]职位的简历，这是职位描述[插入职位描述]，我该如何改进我的简历以适应这一职位，确保我的经验使用主动语言来回答我做了什么，我是如何做到的，同时用数字量化影响？”。
-
可以考虑获取 Tech Resume。
求职者来自各种背景，让你的背景闪耀出来。
-
如果你是自学成才，强调你通过实际项目获得的动手经验。
-
如果你是应届毕业生，突出你的学术成就和任何实习经历，以展示你的基础知识。
-
如果你从其他职业转行，重点介绍你过去的经验或学习如何使你适合数据科学角色。
-
如果你在这个领域已经有经验，详细描述你在数据科学中的具体成就很重要。
你通往数据科学的独特路径是你的优势，所以要确保它在你的求职策略中闪耀出来。
常见的简历错误，你需要避免
-
未针对职位进行定制： 未能根据特定职位及其要求定制简历。我之前提过，但这是第一位！
-
难以扫描的简历： 过于复杂的布局使招聘人员难以快速找到重要信息。
-
过于花哨的简历： 过多的图形或颜色可能会分散对内容的注意力。
-
格式不一致： 字体、大小或样式的不统一会使简历显得不专业。
-
过度加粗： 使用过多的粗体文字可能会让读者感到不适，并稀释关键点的影响。
-
忽略你的受众： 使用内部缩写或行话可能会让读者感到困惑。
-
混乱的 phrases： 结构不良的句子或语法错误可能会留下负面印象。
-
长句子： 冗长的文字块或句子可能使简历难以阅读。
-
不可点击的链接： 简历中的链接应可点击，确保没有损坏。只包括更新和相关的 GitHub、LinkedIn 或个人网站链接。
-
链接中的完整网址： 避免使用长网址，确保链接简洁明了。
优秀的简历能够立即吸引读者的注意。它们展示了该职位所需的具体、稀有且有价值的技能或经验。它们不是通用的，而是专门为你申请的职位量身定制的。它们展示了与该职位相关的清晰而坚实的职业发展。
为了帮助你入门，查看我的一步步制作优秀简历的指南 ⬇ (我展示了让我在 Spotify 找到工作的确切简历)
## 我如何通过简历进入 Spotify 的 0.1%被录取者
成功申请科技行业数据科学职位的秘诀
[towardsdatascience.com
3. 网络建立将快速引导你到招聘人员桌前
“单靠申请是没有用的。你必须建立人际网络。”
不仅仅是我在说，这也是数据科学社区的共识。你可以自己看看。
这是一个在 Blind 上匿名讨论职业相关话题的帖子。我偶然发现了这篇关于数据科学当前就业市场的帖子，里面有经验丰富的数据科学家分享了宝贵的建议。我将这些建议与我从认识的数据科学家那里了解到的相似经历混合在本文中分享——作者图像
“我那家小创业公司发布一个职位招聘——结果是什么？一天内收到 700 份简历。”
你如何期望在 700 份简历中脱颖而出（这只是一天的数量）？目前就业市场充斥着候选人。虽然需求依然存在并且在增长，但供应增长的速度更快，因此出现了“过度饱和”现象。
这种不平衡是现在找工作变得更加困难的原因。这是任何蓬勃发展的领域的常见轨迹，数据科学也不例外。它只是进入了生命周期的这一阶段。
“但是你知道谁立即获得了电话面试吗？是我的推荐人。”
网络不仅仅是关于你认识谁；它还涉及到建立能够带来推荐的联系，这大大增加了你获得面试的机会。
推荐是求职市场的圣杯
我的一个朋友最近获得了苹果一项高度竞争职位的推荐。她只是通过 LinkedIn 联系了一位苹果员工，询问关于职位的信息。结果怎么样？她获得了推荐。她甚至没有要求！那个人非常友善，主动提供了帮助。
这是我朋友告诉我她在苹果获得了推荐的消息——图片由作者提供
推荐很重要，因为它们几乎总是能保证你的简历会得到详细审阅。至少，现在会有人看到。那已经让你比 90%的申请者处于更好的位置。
如果你申请了符合你资格的职位，完善了你的简历，并获得了推荐，你被招聘者看到的机会会显著增加。
但获取推荐并不是直接的。你显然不能只是在 LinkedIn 上随便找人问推荐。他们很可能不会回复，因为他们不认识你。
这时，磨练你的网络技能变得至关重要。许多技术人员可能对网络交流感到不舒服，但这在现在是必需的。
我已经在一篇详细的文章中分享了网络策略如何帮助我获得育碧数据科学家的角色。务必查看一下 ⬇
## 网络帮我在技术行业找到了工作，即使在被拒绝之后，这里是我做到这一点的方式
为什么你的 LinkedIn 网络游戏可能会阻碍你的成功（以及如何解决这个问题）
[towardsdatascience.com
这是如何在 LinkedIn 上进行网络交流的总结
-
检查一下你是否在 LinkedIn 上有在你目标的 5-6 家公司中的第一或第二度连接。
-
主动联系他们，看看他们是否愿意聊天或喝咖啡。许多人通常愿意提供帮助。如果你的网络不强，现在是开始建立它的时候——永远不嫌晚。
-
从这些人那里了解公司的重点领域。这些知识在面试中也很有帮助。面试官会欣赏那些对公司进行了彻底研究的候选人。
图片由作者提供
由于他们接受了网络交流的不适，资格较低的人会比你更容易获得工作
在我之前的一篇文章中，我讨论了网络的重要性，这引发了一些激烈的反响。很多人对网络在求职市场中的作用超过纯粹能力的不公平性感到沮丧。
评论中的一个例子提到了一位科技影响者，他在比自己更有技能的候选人中获得了职位。这看似不公平，但这再次展示了如今网络的重要性。
我理解那种愤慨的感觉，但职业世界并不总是公平的。这是一个难以接受的事实，但与之和解是至关重要的。关键是要聪明地应对这一系统，而不是成为其受害者。
如果一个技能平平的人能通过人脉找到工作，那么拥有扎实技能和强大网络的人有多大潜力呢？为什么要让技能不如你的那些人因为更擅长利用社交资源而超越你呢？
科技行业不再是那些喜欢避免人际互动的人的避风港。就像其他行业一样，建立联系和网络已经成为求职过程中不可或缺的一部分。
适应这一市场环境至关重要。你无法改变游戏规则，所以只能相应地制定策略。当然，你也可以选择拒绝。如果你不愿意走出舒适区来拥抱网络建设，你可以选择忽略这个建议，并对 1%的成功率感到满意。
最后一言 + 总结
现在，我坐在 Spotify 的办公桌前，气氛相当紧张。我们正在经历一次大规模裁员——1500 人，约占我们 17%，正在失去工作。这是我在短短 10 个月内经历的第三次大裁员。
不确定这次我是否能度过难关。我在法国，这里的裁员方式不同。可能需要几周甚至几个月才能知道我是否会加入裁员行列。这是一个既奇怪又艰难的境地，但让我深刻反思现在找工作有多么困难。
我写这份路线图时也考虑了自己。如果我的工作被裁掉，我会按这种方式操作。我以前使用过这种策略，并且有效。如果需要，我还会再次使用。没有其他的秘密。
TL;DR
-
对齐期望与现实：了解求职过程中可控和不可控的因素。重点瞄准与你的技能水平匹配的职位，并相应优化你的简历。
-
简历制作大师：专门为 5 到 6 家挑选的公司量身定制简历。在简历中强调可量化的影响，并使用积极的语言。避免常见的简历错误。
-
拥抱网络：认识到网络在获取职位推荐和面试中的关键作用（即使在科技行业也是如此）。利用你的 LinkedIn 连接寻找潜在的推荐。通过网络了解公司的关注点，提升你的面试准备。
这个策略假设你符合你申请的工作的资格。如果不符合，那就回去学习，等准备好后再尝试。
否则，缩小你的搜索范围。我不建议这样做；这应该作为最后的选择，我总是鼓励读者将求职与职业目标对齐，而不是反过来。
请记住，虽然就业市场正在慢慢恢复，但确实在改善，未来会更好，祝你好运，朋友！
我有礼物送给你🎁！
注册我的newsletter K’s DataLadder，你将自动获得我每天在大科技公司工作中使用的终极 SQL 备忘单以及另一个神秘礼物！
我每周分享作为技术领域数据科学家的经验，包括实用技巧、技能和故事，旨在帮助你提升自我——因为没人真正了解，直到他们亲身经历！
如果你还没这样做的话
很快见！
扫一扫
1769
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
