Python-数学应用第二版-二-

Python 数学应用第二版（二）

原文：annas-archive.org/md5/f7fb298b3279144e9a56cc1cc6fa87cf

译者：飞龙

协议：CC BY-NC-SA 4.0

第四章：与随机性和概率打交道

在本章中，我们将讨论随机性和概率。我们将从简要探索概率的基本概念开始，通过从数据集中选择元素来学习概率。然后，我们将学习如何使用 Python 和 NumPy 生成（伪）随机数，并如何根据特定的概率分布生成样本。最后，我们将通过讨论一些高级主题，涵盖随机过程和贝叶斯技术，并使用马尔科夫链蒙特卡洛（MCMC）方法来估计简单模型的参数，来结束本章。

概率是对特定事件发生的可能性进行量化。我们经常直观地使用概率，尽管有时正式的理论可能相当违反直觉。概率论旨在描述随机变量的行为，随机变量的值未知，但该随机变量的值的概率会在已知的某些（范围的）值中取值。这些概率通常呈现为几种概率分布之一。可以说，最著名的概率分布之一是正态分布，它例如可以描述某一特征在一个大群体中的分布情况。

在第六章《数据与统计学应用》中，我们将再次讨论概率，并且将介绍统计学。在这里，我们将应用概率论来量化误差，并构建分析数据的系统化理论。

本章将涵盖以下实例：

• 随机选择项目

• 生成随机数据

• 更改随机数生成器

• 生成正态分布的随机数

• 与随机过程打交道

• 使用贝叶斯技术分析转化率

• 使用蒙特卡洛模拟估计参数

技术要求

本章需要使用标准的科学 Python 包：NumPy、Matplotlib 和 SciPy。我们还需要 PyMC 包来完成最终的实例。你可以通过你喜欢的包管理工具安装，比如pip：

python3.10 -m pip install pymc

此命令将安装最新版本的 PyMC，目前版本为 4.0.1\。此包提供了概率编程功能，涉及通过随机生成的数据进行多次计算，以了解解决方案的可能分布。

注意

在上一版中，当前版本的 PyMC 是 3.9.2，但自那时以来，PyMC 4.0 版本已经发布，并且在此更新中名称恢复为 PyMC，而不是 PyMC3\。

本章的代码可以在 GitHub 仓库的Chapter 04文件夹中找到，链接：github.com/PacktPublishing/Applying-Math-with-Python-2nd-Edition/tree/main/Chapter%2004。

随机选择项目

概率和随机性的核心思想是从某种集合中选择一个项目。正如我们所知道的，从集合中选择一个项目的概率量化了该项目被选择的可能性。随机性描述了根据概率从集合中选择项目，而没有任何额外的偏差。随机选择的对立面可以被描述为确定性选择。一般来说，使用计算机复制一个纯随机过程是非常困难的，因为计算机及其处理本质上是确定性的。然而，我们可以生成伪随机数序列，当这些序列正确构造时，能够展示出一个合理的随机性近似。

在本配方中，我们将从一个集合中选择项目，并学习一些与概率和随机性相关的关键术语，这些术语将在本章中贯穿始终。

准备工作

Python 标准库包含一个用于生成（伪）随机数的模块，叫做 random，但在本配方和整个章节中，我们将使用 NumPy 的 random 模块。NumPy random 模块中的例程可以用于生成随机数数组，并且比标准库中的相应模块更具灵活性。像往常一样，我们会使用 np 别名导入 NumPy。

在继续之前，我们需要固定一些术语。样本空间是一个集合（一个没有重复元素的集合），事件是样本空间的一个子集。事件发生的概率用 表示，记作 ，它是一个介于 0 和 1 之间的数字。概率为 0 表示事件永远不会发生，而概率为 1 表示事件肯定会发生。整个样本空间的概率必须为 1。

当样本空间是离散的时，概率只是介于 0 和 1 之间的数字，关联到每个元素，其中所有这些数字的总和为 1。这样就给出了从一个集合中选择单个项目（一个包含单一元素的事件）的概率的意义。我们将在这里考虑从离散集合中选择项目的方法，并将在生成正态分布的随机 数的配方中处理连续的情况。

如何做...

执行以下步骤从容器中随机选择项目：

1. 第一步是设置随机数生成器。此时，我们将使用 NumPy 的默认随机数生成器，这在大多数情况下都是推荐的。我们可以通过调用 NumPy random 模块中的 default_rng 例程来实现，这将返回一个随机数生成器实例。我们通常在不指定种子的情况下调用这个函数，但在本配方中，我们将添加一个 12345 种子，以便使我们的结果具有可重复性：

   rng = np.random.default_rng(12345) 
   

   # changing seed for repeatability
   

2. 接下来，我们需要创建要选择的数据和概率。如果你已经存储了数据，或者希望选择具有相等概率的元素，可以跳过此步骤：

   data = np.arange(15)
   

   probabilities = np.array(
   

        [0.3, 0.2, 0.1, 0.05, 0.05, 0.05, 0.05, 0.025,
   

        0.025, 0.025, 0.025, 0.025, 0.025, 0.025, 0.025]
   

   )
   

作为快速的基本检查，我们可以使用断言来检查这些概率的和确实为 1：

assert round(sum(probabilities), 10) == 1.0,
    "Probabilities must sum to 1"

现在，我们可以在随机数生成器 rng 上使用 choice 方法，根据刚才创建的概率从 data 中选择样本。对于这个选择，我们希望开启替换，因此多次调用该方法可以从整个 data 中进行选择：

selected = rng.choice(data,p=probabilities,replace=True)

要从 data 中选择多个项目，我们还可以提供 size 参数，指定要选择的数组形状。这个参数和其他 NumPy 数组创建函数中的 shape 关键字参数起到相同的作用。提供给 size 的参数可以是一个整数，也可以是一个整数元组：

selected_array = rng.choice(data, p=probabilities,  replace=True, size=(5, 5))
#array([[ 1, 6, 4, 1, 1],
#         [ 2, 0, 4, 12, 0],
#         [12, 4, 0, 1, 10],
#         [ 4, 1, 5, 0, 0],
#         [ 0, 1, 1, 0, 7]])

我们可以看到，在采样数据中，0 和 1 的出现次数似乎更多，我们分别为其分配了 0.3 和 0.2 的概率。有趣的是，尽管 12 出现的概率是 2 的一半，结果我们只看到了一个 2，而出现了两个 12。这并不是问题；较大的概率并不保证某个特定数字在样本中出现的次数，只有在大量样本中，我们才会预期大约看到两倍于 12 的 2。

它是如何工作的…

default_rng 函数创建了一个新的 random 模块。然而，通常建议显式地使用 default_rng 创建生成器，或者自己创建一个 Generator 实例。以这种方式更加明确，这也更符合 Python 的风格，并且应该能产生更可重复的结果（在某种程度上）。

种子是传递给随机数生成器的一个值，用于生成随机值。生成器基于种子以完全确定的方式生成数字序列。这意味着，如果提供相同的种子，两个相同的 PRNG 实例将生成相同的随机数序列。如果没有提供种子，生成器通常会生成一个依赖于用户系统的种子。

Generator 类是 NumPy 中一个对低级伪随机比特生成器的封装，实际的随机数就是在这里生成的。在 NumPy 的最新版本中，默认的 PRNG 算法是 128 位的置换同余生成器。相比之下，Python 内置的 random 模块使用的是梅森旋转算法 PRNG。有关不同 PRNG 算法选项的更多信息，请参阅更改随机数 生成器的食谱。

Generator实例上的choice方法根据底层BitGenerator生成的随机数进行选择。可选的p关键字参数指定了与提供的数据中每个项目相关的概率。如果没有提供此参数，则假定均匀概率，即每个项目被选中的概率相等。replace关键字参数指定选择是否应有放回或无放回。我们启用了放回选项，以便同一个元素可以被选中多次。choice方法使用生成器提供的随机数来进行选择，这意味着两个使用相同种子并且类型相同的伪随机数生成器（PRNG）在使用choice方法时会选中相同的项目。

从袋子中选择可能的点是理解离散概率的一个好方法。在这里，我们为每个有限数量的点分配一个特定的权重——例如，点数的倒数——这些权重的总和为 1。采样是根据概率分配的权重随机选择点的过程（我们也可以为无限集合分配离散概率，但这会更复杂，因为权重必须总和为 1，并且这在计算上也不实用）。

还有更多...

choice方法还可以通过传递replace=False参数来创建指定大小的随机样本。这可以确保从数据中选择不同的项目，这对于生成随机样本非常有用。例如，可能会用它来从整个用户组中选择用户来测试新版本的界面；大多数样本统计技术依赖于随机选择的样本。

生成随机数据

许多任务涉及生成大量的随机数，这些随机数在最基本的形式下，可能是整数或浮动精度（双精度）的浮点数，范围在 之间。理想情况下，这些数字应该均匀选择，这样如果我们绘制大量这些数字，它们将在范围 内大致均匀分布。

在这个食谱中，我们将展示如何使用 NumPy 生成大量随机整数和浮点数，并通过直方图展示这些数字的分布。

准备工作

在开始之前，我们需要从 NumPy 的random模块导入default_rng方法，并创建一个默认的随机数生成器实例来在食谱中使用：

from numpy.random import default_rng
rng = default_rng(12345) # changing seed for reproducibility

我们在随机选择项目食谱中已经讨论过这个过程。

我们还通过plt别名导入了 Matplotlib 的pyplot模块。

如何操作...

执行以下步骤以生成均匀的随机数据并绘制直方图以理解其分布：

1. 要生成介于 0 和 1 之间的随机浮点数（包括 0 但不包括 1），我们使用rng对象上的random方法：

   random_floats = rng.random(size=(5, 5))
   

   # array([[0.22733602, 0.31675834, 0.79736546, 0.67625467, 0.39110955],
   

   #           [0.33281393, 0.59830875, 0.18673419, 0.67275604, 0.94180287],
   

   #           [0.24824571, 0.94888115, 0.66723745, 0.09589794, 0.44183967],
   

   #           [0.88647992, 0.6974535 , 0.32647286, 0.73392816, 0.22013496],
   

   #           [0.08159457, 0.1598956 , 0.34010018, 0.46519315, 0.26642103]])
   

2. 要生成随机整数，我们使用 rng 对象上的 integers 方法。这样会返回指定范围内的整数：

   random_ints = rng.integers(1, 20, endpoint=True, size=10)
   

   # array([12, 17, 10, 4, 1, 3, 2, 2, 3, 12])
   

3. 为了检查随机浮点数的分布，我们首先需要生成一个大数组的随机数，正如我们在步骤 1中所做的那样。虽然这并不是严格必要的，但较大的样本能更清楚地显示分布。我们按以下方式生成这些数字：

   dist = rng.random(size=1000)
   

4. 为了显示我们生成的数字的分布，我们绘制了数据的直方图：

   fig, ax = plt.subplots()
   

   ax.hist(dist, color="k", alpha=0.6)
   

   ax.set_title("Histogram of random numbers")
   

   ax.set_xlabel("Value")
   

   ax.set_ylabel("Density")
   

结果的图形如图 4.1所示。正如我们所看到的，数据在整个范围内大致均匀分布：

](tos-cn-i-73owjymdk6/e0b34b2ff3cb44b39740c370195593a1)

图 4.1 – 随机生成的介于 0 和 1 之间的随机数的直方图

随着采样点数量的增加，我们预计这些条形图会“平整”并越来越像我们期望的均匀分布的平线。可以将其与图 4.2中的 10,000 个随机点的直方图进行比较：

](github.com/OpenDocCN/f…)

图 4.2 – 10,000 个均匀分布的随机数的直方图

我们可以看到，尽管分布并不是完全平坦的，但在整个范围内分布更加均匀。

它是如何工作的...

Generator 接口提供了三个简单的方法来生成基本的随机数，但不包括我们在随机选择项目一节中讨论的 choice 方法。除了用于生成随机浮点数的 random 方法和用于生成随机整数的 integers 方法外，还有一个用于生成原始随机字节的 bytes 方法。这些方法都会调用底层 BitGenerator 实例的相关方法。每个方法还可以改变生成的数字的数据类型，例如，从双精度浮点数变为单精度浮点数。

还有更多...

Generator 类上的 integers 方法通过添加 endpoint 可选参数，结合了旧版 RandomState 接口中的 randint 和 random_integers 方法的功能（在旧版接口中，randint 方法排除了上限，而 random_integers 方法包括了上限）。Generator 上的所有随机数据生成方法都允许定制生成数据的数据类型，这在旧接口中是无法做到的（该接口在 NumPy 1.17 中引入）。

在图 4.1中，我们可以看到我们生成的数据的直方图在范围内大致均匀分布 。也就是说，所有的条形图大致平齐（由于数据的随机性质，它们不会完全平齐）。这是我们期望的均匀分布随机数的特征，比如通过random方法生成的随机数。我们将在生成正态分布随机数的食谱中更详细地解释随机数的分布。

更换随机数生成器

NumPy 中的random模块提供了多个备用的伪随机数生成器，它们的默认 PRNG 使用 128 位置换同余生成器。虽然这是一个很好的通用随机数生成器，但它可能不足以满足你的特定需求。例如，这种算法与 Python 内部随机数生成器使用的算法差异很大。我们将遵循 NumPy 文档中设定的最佳实践指南，进行可重复但适当随机的仿真运行。

在这个食谱中，我们将向你展示如何切换到备用的伪随机数生成器（PRNG），以及如何在程序中有效地使用种子。

准备工作

像往常一样，我们通过np别名导入 NumPy。由于我们将使用random包中的多个项目，我们也通过以下代码从 NumPy 导入该模块：

from numpy import random

你需要选择 NumPy 提供的一个备用随机数生成器（或者定义你自己的；请参阅本食谱中的*还有更多...*部分）。在这个食谱中，我们将使用MT19937随机数生成器，它使用基于梅森旋转算法的生成器，类似于 Python 内部的随机数生成器所使用的算法。

如何操作...

以下步骤展示了如何以可重复的方式生成种子和不同的随机数生成器：

1. 我们将生成一个SeedSequence对象，该对象可以从给定的熵源可重复地生成新的种子。我们可以像提供default_rng的种子一样，提供我们自己的整数作为熵，或者我们也可以让 Python 从操作系统收集熵。这里我们选择后一种方法来演示它的使用。为此，我们在创建SeedSequence对象时不提供任何额外的参数：

   seed_seq = random.SeedSequence()
   

2. 现在我们有了为整个会话生成随机数生成器种子的手段，接下来记录熵，以便以后在必要时能够重现这个会话。以下是熵的示例；你的结果可能会有所不同：

   print(seed_seq.entropy)
   

   # 9219863422733683567749127389169034574
   

3. 现在，我们可以创建底层的BitGenerator实例，它将为包装的Generator对象提供随机数：

   bit_gen = random.MT19937(seed_seq)
   

4. 接下来，我们围绕这个BitGenerator实例创建一个包装的Generator对象，从而生成一个可用的随机数生成器：

   rng = random.Generator(bit_gen)
   

一旦创建，你就可以像我们在之前的食谱中看到的那样使用这个随机数生成器。

它是如何工作的...

如在 随机选择项 食谱中提到的，Generator 类是一个封装底层 BitGenerator 的类，后者实现了给定的伪随机数算法。NumPy 通过 BitGenerator 类的各种子类提供了几种伪随机数算法的实现：PCG64（默认）；MT19937（如本食谱所示）；Philox；和 SFC64。这些位生成器是在 Cython 中实现的。

PCG64 生成器应提供高性能的随机数生成，并且具有良好的统计质量（在 32 位系统上可能不是这样）。MT19937 生成器比更现代的伪随机数生成器要慢，并且生成的随机数不具有良好的统计属性。然而，这是 Python 标准库 random 模块使用的随机数生成算法。Philox 生成器相对较慢，但生成的随机数质量非常高，而 SFC64 生成器速度较快，质量也相对较好，但其统计属性不如其他生成器。

本食谱中创建的 SeedSequence 对象是一种以独立且可重复的方式为随机数生成器创建种子的方法。特别是，如果你需要为多个并行进程创建独立的随机数生成器，但仍需要能够在稍后重建每个会话以调试或检查结果时，这非常有用。该对象上存储的熵是一个 128 位整数，来自操作系统，并作为随机种子的来源。

SeedSequence 对象允许我们为每个独立的进程或线程创建一个独立的随机数生成器，从而消除可能导致结果不可预测的数据竞争问题。它还生成彼此非常不同的种子值，这有助于避免一些伪随机数生成器（例如 MT19937，它可能会使用两个相似的 32 位整数种子值生成非常相似的流）的问题。显然，当我们依赖这些值的独立性时，拥有两个生成相同或非常相似值的独立随机数生成器会导致问题。

还有更多...

BitGenerator 类作为生成原始随机整数的生成器的通用接口。前面提到的类是在 NumPy 中实现的，并采用 BitGenerator 接口。你也可以创建自己的 BitGenerator 子类，尽管这需要在 Cython 中实现。

注意

有关更多信息，请参考 NumPy 文档 numpy.org/devdocs/reference/random/extending.html#new-bit-generators。

生成正态分布的随机数

在生成随机数据的示例中，我们生成了遵循均匀分布（在 0 和 1 之间，但不包括 1）的随机浮点数。然而，在大多数需要随机数据的情况下，我们需要遵循几种不同的分布之一。粗略来说，分布函数是一个函数，，它描述了一个随机变量的值小于的概率。从实际应用上讲，分布描述了随机数据在一个范围内的分布情况。特别地，如果我们绘制一个遵循特定分布的数据的直方图，那么它应该大致与该分布函数的图形相似。通过实例来看，这一点最为明显。

最常见的分布之一是正态分布，它在统计学中经常出现，并且是我们在第六章《与数据和统计一起工作》中的许多统计方法的基础。在这个示例中，我们将展示如何生成遵循正态分布的数据，并绘制该数据的直方图，以查看分布的形态。

准备工作

与生成随机数据示例一样，我们从 NumPy 的random模块中导入default_rng函数，并创建一个带有种子生成器的Generator实例用于演示：

from numpy.random import default_rng
rng = default_rng(12345)

像往常一样，我们导入 Matplotlib 的pyplot模块为plt，并导入 NumPy 为np。

如何实现...

在接下来的步骤中，我们生成遵循正态分布的随机数据：

1. 我们使用Generator实例上的normal方法生成符合normal分布的随机数据。正态分布有两个参数：位置和尺度。还有一个可选的size参数，用于指定生成数据的形状（有关size参数的更多信息，请参阅生成随机数据示例）。我们生成一个包含 10,000 个值的数组，以获得一个合理大小的样本：

   mu = 5.0 # mean value
   

   sigma = 3.0 # standard deviation
   

   rands = rng.normal(loc=mu, scale=sigma, size=10000)
   

2. 接下来，我们绘制这个数据的直方图。我们增加了直方图中的bins数量。这并非严格必要，因为默认数量（10）已经足够，但它确实能稍微更好地展示数据分布：

   fig, ax = plt.subplots()
   

   ax.hist(rands, bins=20, color="k", alpha=0.6)
   

   ax.set_title("Histogram of normally distributed data")
   

   ax.set_xlabel("Value")
   

   ax.set_ylabel("Density")
   

3. 接下来，我们创建一个函数，生成一系列值的预期密度。这是通过将正态分布的概率密度函数与样本数量（10,000）相乘来实现的：

   def normal_dist_curve(x):
   

       return 10000*np.exp(
   

           -0.5*((x-mu)/sigma)**2)/(sigma*np.sqrt(2*np.pi))
   

4. 最后，我们将预期分布与数据的直方图一起绘制：

   x_range = np.linspace(-5, 15)
   

   y = normal_dist_curve(x_range)
   

   ax.plot(x_range, y, "k--")
   

结果如图 4.3所示。我们可以看到，样本数据的分布与正态分布曲线预期的分布非常接近：

图 4.3 – 从正态分布中抽取的数据的直方图，叠加了预期密度

再次强调，如果我们取更大更大的样本，我们会预期样本的粗糙度会开始平滑并接近预期的密度（如图 4.3中的虚线所示）。

它是如何工作的...

正态分布的概率密度函数由以下公式定义：

这与正态分布函数相关，，根据以下公式：

这个概率密度函数在均值处达到峰值，该均值与位置参数一致，而钟形曲线的宽度由尺度参数决定。我们可以在图 4.3中看到，Generator对象上使用normal方法生成的数据的直方图与预期分布非常吻合。

Generator类使用 256 步锯齿法生成正态分布的随机数据，这比 NumPy 中其他可用的 Box-Muller 或逆 CDF 实现要快。

还有更多...

正态分布是连续概率分布的一个例子，因为它是为实数定义的，且分布函数由积分（而非求和）定义。正态分布（以及其他连续概率分布）的一个有趣特性是，选择任何给定实数的概率为 0。这是合理的，因为只有在测量一个值是否位于给定范围内时，才有意义去计算这个分布中选中该值的概率。

正态分布在统计学中非常重要，主要得益于中心极限定理。简而言之，这个定理指出，具有共同均值和方差的独立同分布（IID）随机变量的和，最终将趋近于一个具有共同均值和方差的正态分布。无论这些随机变量的实际分布如何，这一结论始终成立。这使得我们在许多情况下，即使变量的实际分布不一定是正态分布，仍然可以使用基于正态分布的统计检验（但我们在引用中心极限定理时需要非常谨慎）。

除了正态分布之外，还有许多其他连续概率分布。我们已经遇到过 0 到 1 范围内的均匀分布。更一般地，均匀分布在范围内，其概率密度函数由以下方程给出：

其他常见的连续概率密度函数包括指数分布、贝塔分布和伽马分布。每种分布在Generator类中都有对应的方法，用来从该分布生成随机数据。这些方法通常根据分布的名称命名，全部使用小写字母，因此对于上述分布，对应的方法是exponential、beta和gamma。这些分布每个都有一个或多个参数，例如正态分布的地点和尺度参数，它们决定了分布的最终形状。你可能需要查阅 NumPy 文档（numpy.org/doc/1.18/reference/random/generator.html#numpy.random.Generator）或其他资料，以了解每种分布所需的参数。NumPy 文档还列出了可以生成随机数据的概率分布。

处理随机过程

在这个教程中，我们将探讨一个简单的随机过程示例，它模拟了公交车到达一个站点的数量随时间的变化。这个过程被称为泊松过程。泊松过程，，有一个单一的参数，，通常被称为强度或速率，且在给定时间，取值的概率由以下公式给出：

这个方程描述了到时间为止，辆公交车已经到达的概率。从数学上讲，这个方程意味着服从参数为的泊松分布。然而，有一种简单的方法可以通过对服从指数分布的到达间隔时间求和来构建泊松过程。例如，设为第()-次到达和第-次到达之间的时间，这些时间服从参数为的指数分布。现在，我们得到以下方程：

这里，数字是最大值，使得。这是我们将在这个教程中逐步构建的内容。我们还将通过取到达间隔时间的均值来估计这个过程的强度。

准备工作

在开始之前，我们从 NumPy 的random模块中导入default_rng函数，并创建一个新的随机数生成器，并为演示目的设置种子：

from numpy.random import default_rng
rng = default_rng(12345)

除了随机数生成器外，我们还导入了 NumPy 库作为np，并将 Matplotlib 的pyplot模块导入为plt。我们还需要确保安装了 SciPy 库。

如何操作…

以下步骤展示了如何使用泊松过程模拟公交车的到达：

1. 我们的第一项任务是通过从指数分布中采样数据来创建样本到达时间间隔。NumPy Generator类中的exponential方法需要一个scale参数，公式为 ，其中 是速率。我们选择速率为 4，并创建 50 个样本到达时间间隔：

   rate = 4.0
   

   inter_arrival_times = rng.exponential(
   

       scale=1./rate, size=50)
   

2. 接下来，我们使用 NumPy add通用函数的accumulate方法计算实际到达时间。我们还创建了一个包含整数 0 到 49 的数组，表示每个时刻的到达次数：

   arrivals = np.add.accumulate(inter_arrival_times)
   

   count = np.arange(50)
   

3. 接下来，我们使用step绘图方法绘制到达次数随时间变化的图：

   fig1, ax1 = plt.subplots()
   

   ax1.step(arrivals, count, where="post")
   

   ax1.set_xlabel("Time")
   

   ax1.set_ylabel("Number of arrivals")
   

   ax1.set_title("Arrivals over time")
   

结果如图 4.4所示，在图中每条水平线的长度表示到达时间间隔：

图 4.4 – 随时间变化的到达次数，且时间间隔服从指数分布

1. 接下来，我们定义一个函数，用于评估单位时间内的计数概率分布，这里我们取1作为单位时间。这个函数使用了我们在本配方介绍中给出的泊松分布公式：

   def probability(events, time=1, param=rate):
   

       return ((param*time)**events/factorial(
   

           events))*np.exp(- param*time)
   

2. 现在，我们绘制单位时间内的计数概率分布，因为在之前的步骤中我们选择了time=1。我们将在后续的步骤中继续完善这个图：

   fig2, ax2 = plt.subplots()
   

   ax2.plot(N, probability(N), "k", label="True distribution")
   

   ax2.set_xlabel("Number of arrivals in 1 time unit")
   

   ax2.set_ylabel("Probability")
   

   ax2.set_title("Probability distribution")
   

3. 现在，我们继续通过样本数据估计速率。我们通过计算到达时间间隔的均值来实现这一点，对于指数分布来说，均值是尺度参数 的估计值：

   estimated_scale = np.mean(inter_arrival_times)
   

   estimated_rate = 1.0/estimated_scale
   

4. 最后，我们绘制基于这个估计速率的单位时间内计数的概率分布。我们将其绘制在我们在步骤 5中得到的真实概率分布之上：

   ax2.plot(N, probability(
   

       N, param=estimated_rate),
   

       "k--",label="Estimated distribution")
   

   ax2.legend()
   

结果图如图 4.5所示，在图中我们可以看到，除了一个小的偏差，估计的分布与真实分布非常接近：

图 4.5 – 每单位时间到达次数的分布，估计值与真实值

图 4.5中展示的分布遵循了本配方介绍中描述的泊松分布。你可以看到，单位时间内的适中到达次数比大量到达次数更为常见。最可能的到达次数由速率参数 确定，在这个例子中速率为 4.0。

工作原理...

随机过程无处不在。大致来说，随机过程是相关随机变量的系统，通常根据时间对连续随机过程进行索引 ，或根据自然数对离散随机过程进行索引 。许多（离散）随机过程满足马尔可夫性质，这使它们成为马尔可夫链。马尔可夫性质的表述是，过程是无记忆的，即只有当前值对于下一个值的概率是重要的。

泊松过程是一种计数过程，用于计算在一段时间内发生的事件（例如公交车到达）的数量，前提是事件在时间上是随机分布的，并且具有固定参数的指数分布。我们通过从指数分布中抽样到达间隔时间来构建泊松过程，遵循我们在介绍中描述的构建方法。然而，事实证明，这一特性（即到达间隔时间是指数分布的）是所有泊松过程的属性，当它们按照概率的正式定义给出时。

在这个方法中，我们从具有给定rate参数的指数分布中抽取了 50 个点。我们必须做一个小的转换，因为 NumPy 的Generator方法用于从指数分布中抽样时，使用的是相关的scale参数，它是1除以rate参数。拿到这些点后，我们创建一个包含这些指数分布数值的累积和的数组。这就生成了我们的到达时间。实际的泊松过程是图 4**.4中显示的过程，它是到达时间与在该时间点发生的事件数量的组合。

指数分布的均值（期望值）与尺度参数相同，因此从指数分布中抽取样本的均值是一种估计尺度（rate）参数的方法。这个估计不会是完美的，因为我们的样本相对较小。这就是为什么在图 4**.5中的两幅图之间存在小的偏差。

还有更多内容...

随机过程有很多种类型，用于描述各种各样的现实场景。在这个方法中，我们使用泊松过程对到达时间进行了建模。泊松过程是一个连续随机过程，这意味着它是由一个连续变量 参数化的，而不是由一个离散变量 参数化的。泊松过程实际上是马尔可夫链，按照对马尔可夫链的适当推广定义，也是一种更新过程的例子。更新过程是描述在一段时间内发生的事件数量的过程。这里描述的泊松过程就是更新过程的一个例子。

许多马尔可夫链除了满足其定义的马尔可夫性质外，还满足一些其他性质。例如，如果下列等式对所有、和值成立，则该马尔可夫链是齐次的：

简单来说，这意味着随着我们增加步骤数，从一个状态到另一个状态的单步转移概率不会改变。这对于检查马尔可夫链的长期行为非常有用。

构建简单的均匀马尔可夫链非常容易。假设我们有两个状态，和。在任何给定的步骤中，我们可以处于状态或状态之一。我们根据一定的概率在这些状态之间移动。例如，假设从状态到状态的转移概率为 0.4，从状态到状态的转移概率为 0.6。同样，假设从状态到状态的转移概率为 0.2，从状态到状态的转移概率为 0.8。请注意，在这两种情况下，转移的概率和保持相同状态的概率之和都为 1。我们可以用矩阵形式表示从每个状态转移的概率，在这种情况下，使用以下方程：

这个矩阵叫做转移矩阵。这里的思路是，经过一步后处于特定状态的概率是通过将包含状态和（分别是位置 0 和 1）的概率的向量相乘得到的。例如，如果我们从状态开始，则概率向量在索引 0 处为 1，在索引 1 处为 0。那么，经过一步后处于状态的概率为 0.4，处于状态的概率为 0.6。这是我们之前概述的概率所预期的结果。然而，我们也可以使用矩阵公式来写出这个计算：

要获取经过两步后处于任一状态的概率，我们再次将右侧向量与转移矩阵相乘，得到如下结果：

我们可以将这个过程继续进行无限次，以获得一系列状态向量，这些向量构成了我们的马尔可夫链。如果需要，可以通过增加更多的状态来应用这一构建方法，用以建模许多简单的现实世界问题。

使用贝叶斯技术分析转化率

贝叶斯概率允许我们通过考虑数据系统地更新我们对情况的理解（从概率的角度）。用更技术化的语言来说，我们使用数据更新先验分布（我们的当前理解），以获得后验分布。例如，在研究用户浏览网站后购买产品的比例时，这非常有用。我们从先验信念分布开始。为此，我们将使用beta分布，它根据观察到的成功次数（完成购买）与失败次数（未购买）来模拟成功的概率。在这个示例中，我们假设我们的先验信念是我们期望在 100 次浏览中有 25 次成功（75 次失败）。这意味着我们的先验信念遵循 beta(25, 75)分布。假设我们想计算真实成功率至少为 33%的概率。

我们的方法大致分为三个步骤。首先，我们需要理解我们对转化率的先验信念，我们已经决定它遵循 beta(25, 75)分布。通过对先验分布的概率密度函数从 0.33 到 1 进行数值积分，我们计算出转化率至少为 33%的概率。接下来的步骤是应用贝叶斯推理，通过新的信息更新我们的先验信念。然后，我们可以对后验（更新后的）信念进行相同的积分，查看在此新信息下转化率至少为 33%的概率。

在本示例中，我们将展示如何使用贝叶斯技术基于新信息更新我们假设网站的先验信念。

准备工作

和往常一样，我们需要导入 NumPy 和 Matplotlib 包，分别命名为np和plt。我们还需要导入 SciPy 包，命名为sp。

如何操作...

以下步骤展示了如何使用贝叶斯推理估算和更新转化率估算：

1. 第一步是设置先验分布。为此，我们使用来自 SciPy stats 模块的beta分布对象，该对象提供了多种处理贝塔分布的方法。我们从stats模块导入beta分布对象，并将其命名为beta_dist，然后创建一个便捷函数，用于计算概率密度函数：

   from scipy.stats import beta as beta_dist
   

   beta_pdf = beta_dist.pdf
   

2. 接下来，我们需要计算在先验信念分布下，成功率至少为 33%的概率。为此，我们使用 SciPy integrate模块中的quad例程，该例程执行函数的数值积分。我们使用它来对先前步骤中导入的贝塔分布的概率密度函数进行积分，结合我们的先验参数。我们根据先验分布将概率打印到控制台：

   prior_alpha = 25
   

   prior_beta = 75
   

   args = (prior_alpha, prior_beta)
   

   prior_over_33, err = sp.integrate.quad(
   

       beta_pdf, 0.33, 1, args=args)
   

   print("Prior probability", prior_over_33)
   

   # 0.037830787030165056
   

3. 现在，假设我们收到了关于新时间段内成功和失败的信息。例如，在这段时间内，我们观察到 122 次成功和 257 次失败。我们创建新的变量来反映这些值：

   observed_successes = 122
   

   observed_failures = 257
   

4. 要获得使用贝塔分布的后验分布的参数值，我们只需将观察到的成功次数和失败次数分别加到prior_alpha和prior_beta参数中：

   posterior_alpha = prior_alpha + observed_successes
   

   posterior_beta = prior_beta + observed_failures
   

5. 现在，我们重复数值积分，以计算使用后验分布（使用我们先前计算的参数）成功率现在是否超过 33%的概率。我们再次将此概率输出到终端：

   args = (posterior_alpha, posterior_beta)
   

   posterior_over_33, err2 = sp.integrate.quad(
   

       beta_pdf, 0.33, 1, args=args)
   

   print("Posterior probability", posterior_over_33)
   

   # 0.13686193416281017
   

6. 我们可以看到，更新后的后验分布给出的新概率是 14%，而先验分布给出的概率是 4%。这是一个显著的差异，尽管根据这些值，我们仍然无法确信转化率超过 33%。现在，我们将先验和后验分布绘制出来，以可视化这种概率的增加。首先，我们创建一个值的数组，并基于这些值评估我们的概率密度函数：

   p = np.linspace(0, 1, 500)
   

   prior_dist = beta_pdf(p, prior_alpha, prior_beta)
   

   posterior_dist = beta_pdf(
   

       p, posterior_alpha, posterior_beta)
   

7. 最后，我们将第 6 步中计算的两个概率密度函数绘制到新的图表上：

   fig, ax = plt.subplots()
   

   ax.plot(p, prior_dist, "k--", label="Prior")
   

   ax.plot(p, posterior_dist, "k", label="Posterior")
   

   ax.legend()
   

   ax.set_xlabel("Success rate")
   

   ax.set_ylabel("Density")
   

   ax.set_title("Prior and posterior distributions for success rate")
   

结果图如图 4**.6所示，我们可以看到后验分布比先验分布更加狭窄，并且集中在先验分布的右侧：

图 4.6 – 成功率的先验和后验分布，遵循贝塔分布

我们可以看到，后验分布在约 0.3 处达到峰值，但分布的大部分质量都集中在这个峰值附近。

它是如何工作的……

贝叶斯技术通过采取一个先验信念（概率分布），然后利用贝叶斯定理将先验信念与在该先验信念下我们的数据的可能性结合起来，从而形成一个后验（更新后的）信念。这类似于我们在现实生活中的理解方式。例如，当你在某天早上醒来时，可能会有一个信念（来自天气预报或其他途径），认为外面有 40%的机会下雨。打开窗帘后，你看到外面非常多云，这可能表明雨水的可能性更大，因此我们根据这个新数据更新我们的信念，认为下雨的机会是 70%。

要理解这个过程，我们需要了解条件概率。条件概率处理的是给定另一个事件已经发生的情况下，一个事件发生的概率。用符号表示，事件 在事件 已经发生的条件下发生的概率可以写作：

贝叶斯定理是一个强大的工具，可以用以下方式（符号化）表示：

概率 代表我们的先验信念。事件 代表我们已收集的数据，因此 是在我们的先验信念下，数据产生的可能性。概率 代表我们的数据产生的概率， 代表在给定数据后，我们的后验信念。实际上，概率 可能很难计算或估计，因此通常会将上面的强等式替换为贝叶斯定理的比例版本：

在这个例子中，我们假设我们的先验信念服从贝塔分布。贝塔分布的概率密度函数由以下公式给出：

这里， 是伽马函数。该似然性呈二项分布，其概率密度函数由以下公式给出：

这里， 是观察次数， 是成功的次数之一。在这个例子中，我们观察到 次成功和 次失败，这给出了 和 。为了计算后验分布，我们可以利用贝塔分布是二项分布的共轭先验这一事实，来观察贝叶斯定理的比例形式右侧是一个贝塔分布，参数为 和 。这就是我们在例子中使用的内容。贝塔分布作为二项随机变量的共轭先验，使其在贝叶斯统计中非常有用。

我们在这个例子中展示的方法是使用贝叶斯方法的一个相对基础的示例，但它仍然对我们在系统地获取新数据时更新先验信念非常有用。

还有更多…

贝叶斯方法可以用于各种任务，成为一个强大的工具。在这个例子中，我们使用贝叶斯方法根据我们对网站表现的先验信念和用户收集的额外数据来建模网站的成功率。这个例子比较复杂，因为我们将先验信念建模为贝塔分布。这里是另一个使用贝叶斯定理的例子，用简单的概率（介于 0 和 1 之间的数字）来检验两个竞争的假设。

假设你每天回家时都会把钥匙放在同一个地方，但有一天早晨你醒来发现钥匙不在这个地方。经过短暂的寻找后，你找不到它们，于是得出它们一定已经从这个世界消失的结论。我们称这个假设为。现在，无疑解释了数据，，即你找不到钥匙——因此，的可能性（如果你的钥匙消失了，那么你根本不可能找到它们）。另一个假设是你回家时把钥匙放到了别的地方。我们称这个假设为。现在，这个假设同样能解释数据，所以，但实际上，比更为可信。假设你的钥匙完全消失的概率是百万分之一——这显然是个过高的估计，但我们需要保持数字的合理性——而你估计你把它们放到别处的概率是千分之一。计算后验概率，我们得到如下结果：

这突显了这样一个现实：你丢失钥匙的可能性比它们消失的可能性要大 10,000 倍。果不其然，你很快就发现钥匙已经在你的口袋里，因为你早晨曾把它放进过那里。

使用蒙特卡洛模拟估计参数

蒙特卡洛方法广泛描述了利用随机抽样来解决问题的技术。当问题涉及某种不确定性时，这些技术尤其强大。一般方法包括执行大量模拟，每次根据给定的概率分布抽取不同的输入，然后汇总结果，从而比任何单一的样本解提供更好的真实解近似。

MCMC（马尔可夫链蒙特卡洛）是一种特定类型的蒙特卡洛模拟，在这种方法中，我们构造一个马尔可夫链，逐步得到我们所寻求的真实分布的更好的近似。这个过程通过基于精心选择的接受概率，在每个阶段接受或拒绝一个随机采样的拟议状态，目的是构造一个其唯一平稳分布恰好是我们想要找到的未知分布的马尔可夫链。

在这个例子中，我们将使用 PyMC 包和 MCMC 方法来估计一个简单模型的参数。该包将处理运行模拟的大部分技术细节，因此我们无需进一步深入了解不同 MCMC 算法的实际工作原理。

准备工作

如常，我们导入 NumPy 包和 Matplotlib 的 pyplot 模块，分别命名为 np 和 plt。我们还导入并创建一个默认的随机数生成器，使用种子来进行演示，代码如下：

from numpy.random import default_rng
rng = default_rng(12345)

我们还需要 SciPy 包中的一个模块，以及 PyMC 包，它是一个用于概率编程的包。我们将 PyMC 包导入，并命名为 pm：

import pymc as pm

让我们来看一下如何使用 PyMC 包来估计给定观测到的噪声样本的模型参数。

如何做到...

执行以下步骤，使用 MCMC 模拟来估计简单模型的参数，基于样本数据：

1. 我们的第一个任务是创建一个表示我们希望识别的基础结构的函数。在本例中，我们将估计一个二次方程（一个二次多项式）的系数。该函数接受两个参数，一个是固定的范围内的点，另一个是我们希望估计的变量参数：

   def underlying(x, params):
   

        return params[0]*x**2 + params[1]*x + params[2]
   

2. 接下来，我们设置 true 参数和 size 参数，这将决定我们生成的样本中有多少个点：

   size = 100
   

   true_params = [2, -7, 6]
   

3. 我们生成将用于估计参数的样本数据。这个样本由基础数据组成，基础数据是由我们在步骤 1中定义的 underlying 函数生成的，再加上一些符合正态分布的随机噪声。我们首先生成一个范围的 值，这些值在整个过程里保持不变，然后使用 underlying 函数和随机数生成器的 normal 方法来生成样本数据：

   x_vals = np.linspace(-5, 5, size)
   

   raw_model = underlying(x_vals, true_params)
   

   noise = rng.normal(loc=0.0, scale=10.0, size=size)
   

   sample = raw_model + noise
   

4. 在开始分析之前，最好先绘制采样数据，并将基础数据叠加在其上。我们使用 scatter 绘图方法仅绘制数据点（不连接线条），然后用虚线绘制基础的二次结构：

   fig1, ax1 = plt.subplots()
   

   ax1.scatter(x_vals, sample,
   

       label="Sampled data", color="k", 
   

       alpha=0.6)
   

   ax1.plot(x_vals, raw_model,
   

       "k--", label="Underlying model")
   

   ax1.set_title("Sampled data")
   

   ax1.set_xlabel("x")
   

   ax1.set_ylabel("y")
   

结果是 图 4.7，在图中可以看到，即使存在噪声，基础模型的形状仍然可见，尽管该模型的确切参数已不再明显：

图 4.7 – 采样数据与基础模型叠加

1. PyMC 编程的基本对象是 Model 类，通常使用上下文管理器接口来创建。我们还为参数创建了先验分布。在本例中，我们假设我们的先验参数服从均值为 1、标准差为 1 的正态分布。我们需要三个参数，因此提供了 shape 参数。Normal 类创建将用于蒙特卡洛模拟的随机变量：

   with pm.Model() as model:
   

       params = pm.Normal(
   

           "params", mu=1, sigma=1, shape=3)
   

2. 我们为基础数据创建一个模型，可以通过将我们在步骤 6中创建的随机变量param传递给步骤 1中定义的underlying函数来完成。我们还创建了一个处理我们观测值的变量。为此，我们使用Normal类，因为我们知道我们的噪声在基础数据y周围呈正态分布。我们设置标准差为2，并将我们的观测sample数据传递给observed关键字参数（这也在Model上下文中）：

   y = underlying(x_vals, params)
   

   y_obs = pm.Normal("y_obs",
   

       mu=y, sigma=2, observed=sample)
   

3. 为了运行模拟，我们只需要在Model上下文中调用sample例程。我们传递cores参数来加速计算，但保留所有其他参数为默认值：

       trace = pm.sample(cores=4)
   

这些模拟应该执行起来很快。

1. 接下来，我们使用来自 PyMC 的plot_posterior例程绘制后验分布。该例程使用从采样步骤中获得的trace结果，我们已提前使用plt.subplots例程创建了自己的图形和坐标轴，但这不是严格必要的。我们在一个图形上使用了三个子图，并将axs2元组的Axes传递给绘图例程的ax关键字参数：

   fig2, axs2 = plt.subplots(1, 3, tight_layout=True)
   

   pm.plot_posterior(trace, ax=axs2, color="k")
   

结果的图示显示在图 4**.8中，您可以看到这些分布大致是正态分布，均值与真实参数值相似：

图 4.8 – 估计参数的后验分布

1. 现在，从trace结果中提取每个估计参数的均值。我们通过访问trace上的后验属性来获取估计的参数，然后对params项目使用mean方法（使用axes=(0,1)在所有链和所有样本上进行平均），并将其转换为 NumPy 数组。我们在终端中打印这些估计的参数：

   estimated_params = trace.posterior["params"].mean(
   

       axis=(0, 1)). to_numpy()
   

   print("Estimated parameters", estimated_params)
   

   # Estimated parameters [ 2.03220667 -7.09727509  5.27548983]
   

2. 最后，我们使用估计的参数生成我们的估计基础数据，通过将的值和估计的参数传递给步骤 1中定义的underlying函数。然后，我们将这个估计的基础数据与真实的基础数据一起绘制在同一坐标轴上：

   estimated = underlying(x_vals, estimated_params)
   

   fig3, ax3 = plt.subplots()
   

   ax3.plot(x_vals, raw_model, "k", label="True model")
   

   ax3.plot(x_vals, estimated, "k--", label="Estimated model")
   

   ax3.set_title("Plot of true and estimated models")
   

   ax3.set_xlabel("x")
   

   ax3.set_ylabel("y")
   

   ax3.legend()
   

结果的图示显示在图 4**.9中，在这个范围内，这两个模型之间只有一个小的差异：

图 4.9 – 真实模型和估计模型绘制在同一坐标轴上

在图 4**.9中，我们可以看到真实模型和估计模型之间存在一个小的差异。

它是如何工作的……

这个配方中有趣的部分可以在Model上下文管理器中找到。这个对象跟踪随机变量，协调模拟并保持状态。上下文管理器为我们提供了一种方便的方式来将概率变量与周围的代码分离。

我们从为表示参数的随机变量（共有三个）的分布提出先验分布开始。我们提出了正态分布，因为我们知道这些参数不能偏离值 1 太远（例如，通过查看我们在步骤 4中生成的图形，我们可以看出这一点）。使用正态分布将给那些接近当前值的数值更高的概率。接下来，我们加入与观察数据相关的细节，这些数据用于计算接受概率，决定是接受还是拒绝某个状态。最后，我们使用sample例程启动采样器。这将构建马尔可夫链并生成所有的步长数据。

sample例程根据将要模拟的变量类型设置采样器。由于正态分布是一个连续变量，sample例程选择了无转弯采样器（NUTS）。这是一种适用于连续变量的通用采样器。NUTS 的常见替代方法是 Metropolis 采样器，后者虽然在某些情况下比 NUTS 速度更快，但可靠性较差。PyMC 文档建议尽可能使用 NUTS。

一旦采样完成，我们绘制了马尔可夫链给出的状态（即后验分布）的图形，以查看我们生成的近似值的最终形态。我们可以看到，所有三个随机变量（参数）都呈现出大致正确值附近的正态分布。

在背后，PyMC 使用 Aesara——PyMC3 中 Theano 的继任者——来加速计算。这使得 PyMC 可以在图形处理单元（GPU）上进行计算，而不是在中央处理单元（CPU）上，从而大大提升计算速度。

还有更多内容...

蒙特卡罗方法非常灵活，我们这里给出的例子是它可以应用的一种特定情况。蒙特卡罗方法的一个更典型的基本应用例子是估算积分的值——通常被称为蒙特卡罗积分。一个非常有趣的蒙特卡罗积分应用是估算！的值。让我们简要看一下这是如何实现的。

首先，我们取一个单位圆，其半径为 1，因此面积为 。我们可以将这个圆形放置在一个边长为 2 的正方形内部，其顶点分别位于 、、 和 。该正方形的面积为 4，因为边长为 2。现在，我们可以在这个正方形上均匀生成随机点。当我们这样做时，任何一个随机点位于给定区域内的概率与该区域的面积成正比。因此，可以通过将位于该区域内的随机生成点的比例乘以正方形的总面积来估算区域的面积。特别地，我们可以通过将位于圆内的随机生成点的数量乘以 4，并除以我们生成的点的总数，来估算该圆的面积（当半径为 1 时，这个面积是 ）。

我们可以轻松地用 Python 编写一个执行此计算的函数，可能如下所示：

import numpy as np
from numpy.random import default_rng
def estimate_pi(n_points=10000):
    rng = default_rng()
    points = rng.uniform(-1, 1, size=(2, n_points))
    inside = np.less(points[0, :]**2 + points[1, :]**2, 1)
    return 4.0*inside.sum() / n_points

仅运行一次此函数即可给出π的合理近似值：

estimate_pi()  # 3.14224

我们可以通过使用更多的点来提高估计的准确性，但我们也可以多次运行这个过程并对结果进行平均。我们来运行这个模拟 100 次，并计算结果的平均值（我们将使用并发未来（concurrent futures）来并行化这个过程，这样如果需要，我们可以运行更大数量的样本）：

from statistics import mean
results = list(estimate_pi() for _ in  range(100))
print(mean(results))

运行一次此代码会打印出估计值 ，其值为 3.1415752，这是对真实值的一个更精确估计。

另见

PyMC 包含许多功能，这些功能通过大量示例进行了文档说明（docs.pymc.io/）。还有另一个基于 TensorFlow 的概率编程库（www.tensorflow.org/probability）。

进一步阅读

一本关于概率和随机过程的好书是：

• Grimmett, G. 和 Stirzaker, D. (2009). Probability and random processes. 第 3 版. 牛津：牛津 大学出版社。

贝叶斯定理和贝叶斯统计的简单介绍如下：

• Kurt, W. (2019). Bayesian statistics the fun way. San Francisco, CA: No Starch Press, Inc.

第五章：处理树与网络

网络是包含 节点 和节点对之间 边 的对象。它们可以用来表示各种现实世界的情况，如分配和调度。数学上，网络对于可视化组合问题非常有用，并且它有着丰富且迷人的理论。

当然，网络有多种类型。我们将主要处理简单网络，其中边连接两个不同的节点（因此没有自环），每两个节点之间最多只有一条边，并且所有边都是双向的。树 是一种特殊的网络，其中没有环；也就是说，没有节点列表，每个节点都通过边连接到下一个节点，最后一个节点又连接到第一个节点。树在理论上特别简单，因为它用最少的边连接多个节点。完全网络 是一种每个节点都通过边与其他每个节点相连的网络。

网络可以是有向的，其中每条边都有一个源节点和一个目标节点，或者可以承载额外的属性，例如权重。加权网络在某些应用中尤其有用。还有一些网络允许在两个给定节点之间有多条边。

在本章中，我们将学习如何创建、操作和分析网络，并应用网络算法解决各种问题。

注意

在文献中，特别是在数学文本中，网络更常被称为 图。节点有时被称为 顶点。我们更倾向于使用“网络”这一术语，以避免与更常见的将图表示为函数图像的用法混淆。

本章将涵盖以下内容：

• 在 Python 中创建网络

• 可视化网络

• 获取网络的基本特征

• 生成网络的邻接矩阵

• 创建有向和加权网络

• 查找网络中的最短路径

• 定量化网络中的聚类

• 给网络上色

• 查找最小生成树和支配集

让我们开始吧！

技术要求

在本章中，我们将主要使用 NetworkX 包来处理树和网络。您可以通过您喜欢的包管理器（例如 pip）来安装此包：

python3.10 -m pip install networkx

我们通常使用 nx 别名导入此模块，遵循官方 NetworkX 文档中建立的惯例（networkx.org/documentation/stable/），使用以下 import 语句：

import networkx as nx

本章的代码可以在本书的 GitHub 仓库中的 Chapter 05 文件夹找到，链接为 github.com/PacktPublishing/Applying-Math-with-Python-2nd-Edition/tree/main/Chapter%2005。

在 Python 中创建网络

为了解决可以表示为网络问题的各种问题，我们需要一种在 Python 中创建网络的方式。为此，我们将利用 NetworkX 包及其提供的例程和类来创建、操作和分析网络。

在本配方中，我们将创建一个表示网络的 Python 对象，并向该对象添加节点和边。

准备工作

正如我们在技术要求部分提到的，我们需要将 NetworkX 包以nx别名导入。我们可以使用以下import语句来完成：

import networkx as nx

如何实现...

按照以下步骤创建一个简单图的 Python 表示：

1. 我们需要创建一个新的Graph对象，用于存储构成图的节点和边：

   G = nx.Graph()
   

2. 接下来，我们需要使用add_node方法为网络添加节点：

   G.add_node(1)
   

   G.add_node(2)
   

3. 为了避免重复调用该方法，我们可以使用add_nodes_from方法从可迭代对象（如列表）中添加节点：

   G.add_nodes_from([3, 4, 5, 6])
   

4. 接下来，我们需要使用add_edge或add_edges_from方法为我们已添加的节点之间添加边，这些方法分别是添加单条边或一组边（作为元组）：

   G.add_edge(1, 2)  # edge from 1 to 2
   

   G.add_edges_from([(2, 3),(3, 4)(3, 5),(3, 6),
   

       (4,5),(5,6)])
   

5. 最后，我们必须通过访问nodes和edges属性，分别获取图中当前节点和边的视图：

   print(G.nodes)
   

   print(G.edges)
   

   # [1, 2, 3, 4, 5, 6]
   

   # [(1, 2), (2, 3), (3, 4), (3, 5), (3, 6), (4, 5), (5, 6)]
   

它是如何工作的...

NetworkX 包添加了多个类和例程，用于使用 Python 创建、操作和分析网络。Graph类是表示不包含多个边的网络的最基本类，其中边是无向的（双向的）。

一旦创建了一个空白的Graph对象，我们就可以使用本配方中描述的方法添加新的节点和边。在本配方中，我们创建了包含整数值的节点。然而，节点可以包含任何可哈希的 Python 对象，除了None。此外，还可以通过传递关键字参数给add_node方法来为节点添加关联数据。在使用add_nodes_from方法时，可以通过提供包含节点对象和属性字典的元组列表来添加属性。add_nodes_from方法对于批量添加节点非常有用，而add_node方法则适用于将单个节点附加到现有网络中。

网络中的一条边是一个包含两个（不同）节点的元组。在一个简单的网络中，例如由基本的Graph类表示的网络，任何两个给定节点之间最多只能有一条边。这些边是通过add_edge或add_edges_from方法添加的，分别是添加单个边或添加一组边。与节点类似，边也可以通过属性字典来保存任意的关联数据。特别是，当添加边时，可以通过提供weight属性来添加权重。我们将在创建有向和加权 网络的配方中提供有关加权图的更多细节。

nodes 和 edges 属性分别包含构成网络的节点和边。nodes 属性返回一个 NodesView 对象，这是一个类似字典的接口，提供节点及其关联数据的访问。同样，edges 属性返回一个 EdgeView 对象。我们可以使用这个对象检查个别边及其相关数据。

还有更多...

Graph 类表示 简单网络，它是节点最多由一条边连接且边不具有方向的网络。我们将在 创建有向和加权网络 食谱中讨论有向网络。还有一个单独的类用于表示可以在一对节点之间有多条边的网络，称为 MultiGraph。所有网络类型都允许自环，这在文献中的 简单网络 中有时是不允许的，简单网络通常指的是没有自环的无向网络。

所有网络类型都提供了多种方法来添加节点和边，并检查当前的节点和边。还有方法可以将网络复制到其他类型的网络中，或提取子网络。此外，NetworkX 包还提供了多个实用例程，用于生成标准网络并将子网络添加到现有网络中。

NetworkX 还提供了多种例程，用于将网络读写到不同的文件格式中，例如 GraphML、JSON 和 YAML。例如，我们可以使用 nx.write_graphml 例程将网络写入 GraphML 文件，并使用 nx.read_graphml 例程读取它。

可视化网络

分析网络的常见第一步是绘制网络，这有助于我们识别网络的一些显著特征。（当然，图示可能会产生误导，因此我们不应过于依赖它们进行分析。）

在这个食谱中，我们将描述如何使用 NetworkX 包中的网络绘制功能来可视化网络。

准备工作

对于这个食谱，我们需要按 技术要求 部分的描述，导入 NetworkX 包并使用 nx 别名。我们还需要 Matplotlib 包。为此，我们通常需要使用以下 import 语句导入 pyplot 模块，并将其命名为 plt：

import matplotlib.pyplot as plt

如何操作...

以下步骤概述了如何使用 NetworkX 的绘图例程绘制一个简单的网络对象：

1. 首先，我们将创建一个简单的示例网络进行绘制：

   G = nx.Graph()
   

   G.add_nodes_from(range(1, 7))
   

   G.add_edges_from([
   

       (1, 2), (2, 3), (3, 4), (3, 5),
   

       (3, 6), (4, 5), (5, 6)
   

   ])
   

2. 接下来，我们将为它创建新的 Matplotlib Figure 和 Axes 对象，准备好使用 plt 中的 subplots 例程来绘制网络：

   fig, ax = plt.subplots()
   

3. 现在，我们可以创建一个布局，用于在图形中定位节点。对于这个图形，我们将使用 shell_layout 例程来创建一个壳布局：

   layout = nx.shell_layout(G)
   

4. 我们可以使用draw函数将网络绘制到图形上。由于我们已经创建了一个 Matplotlib Figure和Axes，我们可以提供ax关键字参数。我们还将使用with_labels关键字参数为节点添加标签，并通过pos参数指定我们刚刚创建的布局：

   nx.draw(G, ax=ax, pos=layout, with_labels=True)
   

   ax.set_title("Simple network drawing")
   

结果图像可以在以下图中看到：

图 5.1 – 使用 shell 布局排列的简单网络图

由于本示例中的节点数量相对较少，它们被安排在一个圆圈内。边缘通过线条表示。

工作原理...

draw函数是一个专门的绘图函数，专门用于绘制网络。我们创建的布局指定了每个节点将被放置的坐标。我们使用了shell 布局，它将节点按同心圆的方式排列（此示例中仅使用了一个圆），该布局由网络的节点和边缘决定。默认情况下，draw函数会创建一个随机布局。

draw函数有多个关键字参数，用于定制绘制网络的外观。在本示例中，我们添加了with_labels关键字参数，根据节点持有的对象为节点添加标签。节点持有整数值，这就是为什么前面图中的节点按整数标注的原因。

我们还使用plt.subplots函数单独创建了一组坐标轴。严格来说，这不是必须的，因为draw函数会在没有提供坐标轴的情况下自动创建一个新的图形和坐标轴。

还有更多...

NetworkX 包提供了几个布局生成函数，类似于我们在本示例中使用的shell_layout函数。此布局实际上是一个字典，以节点为索引，其元素是节点应该绘制的位置的 x 和 y 坐标。NetworkX 的布局生成函数表示了适用于大多数情况的常见排列，但如果需要，您也可以创建自定义布局。NetworkX 文档中提供了不同布局创建函数的完整列表。此外，还有快捷的绘图函数，它们会使用特定的布局，而无需单独创建布局；例如，draw_shell函数会绘制使用 shell 布局的网络，这等同于本示例中draw函数的调用。

draw函数接受多个关键字参数来定制图形的外观。例如，您可以使用关键字参数控制节点的大小、颜色、形状和透明度。我们还可以添加箭头（用于有向边）和/或仅绘制网络中的特定节点和边。

获取网络的基本特征

网络除了节点和边的数量之外，还有许多有助于分析图的基本特性。例如，节点的度数是指从该节点出发（或到达）的边的数量。度数越高，表示该节点与网络的其他部分连接越紧密。

在本教程中，我们将学习如何访问网络的基本属性，并计算与网络相关的各种基本度量。

准备工作

和往常一样，我们需要导入 nx 别名下的 NetworkX 包。我们还需要导入 Matplotlib 的 pyplot 模块，别名为 plt。

如何操作...

按照以下步骤访问网络的各种基本特征：

1. 创建我们将在本教程中分析的示例网络，代码如下：

   G = nx.Graph()
   

   G.add_nodes_from(range(10))
   

   G.add_edges_from([
   

       (0, 1), (1, 2), (2, 3), (2, 4),
   

       (2, 5), (3, 4), (4, 5), (6, 7),
   

       (6, 8), (6, 9), (7, 8), (8, 9)
   

   ])
   

2. 接下来，最好绘制网络并将节点排列成圆形布局：

   fig, ax = plt.subplots()
   

   nx.draw_circular(G, ax=ax, with_labels=True)
   

   ax.set_title("Simple network")
   

结果图可以在下图中看到。正如我们所见，网络被分成了两个不同的部分：

图 5.2 – 一个简单的网络，采用圆形排列，包含两个不同的组件

1. 接下来，我们必须打印Graph对象以显示网络的基本信息：

   print(G)
   

   # Name:
   

   # Type: Graph
   

   # Number of nodes: 10
   

   # Number of edges: 12
   

   # Average degree: 2.4000
   

2. 现在，我们可以使用Graph对象的degree属性来获取特定节点的度数：

   for i in [0, 2, 7]:
   

       degree = G.degree[i]
   

       print(f"Degree of {i}: {degree}")
   

   # Degree of 0: 1
   

   # Degree of 2: 4
   

   # Degree of 7: 2
   

3. 我们可以使用 connected_components 例程获取网络的连通组件，该例程返回一个生成器，我们将其转换为列表：

   components = list(nx.connected_components(G))
   

   print(components)
   

   # [{0, 1, 2, 3, 4, 5}, {8, 9, 6, 7}]
   

4. 我们可以使用密度例程计算网络的密度，该例程返回一个介于 0 和 1 之间的浮动值。它表示连接到该节点的边与该节点所有可能的边之间的比例：

   density = nx.density(G)
   

   print("Density", density)
   

   # Density 0.26666666666666666
   

5. 最后，我们可以使用 check_planarity 例程来确定网络是否是平面的——即没有两条边需要交叉——：

   is_planar, _ = nx.check_planarity(G)
   

   print("Is planar", is_planar)
   

   # Is planar True
   

回头看看图 5.2，我们可以看到，确实可以在不交叉任何两条边的情况下绘制这个图。

它是如何工作的...

info 例程生成网络的小结，包括网络的类型（在本教程中是简单的 Graph 类型）、节点和边的数量以及网络中节点的平均度数。可以使用 degree 属性访问网络中节点的实际度数，它提供了一种类似字典的接口来查找每个节点的度数。

一组节点如果该组中的每个节点都通过一条边或一系列边与其他节点相连，则称这组节点为连通的。网络的连通分量是最大的一组相互连通的节点。任何两个不同的连通分量是相互不相交的。每个网络都可以分解为一个或多个连通分量。我们在这个示例中定义的网络有两个连通分量，分别是{0, 1, 2, 3, 4, 5}和{8, 9, 6, 7}。这些连通分量在前面的图中可见，其中第一个连通分量位于第二个连通分量的上方。在这张图中，我们可以沿着网络的边从组件中的任何节点到达任何其他节点；例如，从 0 到 5。

网络的密度衡量的是网络中边的数量与网络中节点数量所能形成的边的总数之间的比例。完全网络的密度为 1，但通常情况下，密度会小于 1。

如果一个网络可以在平面上绘制而没有交叉边，则称该网络为平面网络。五节点的完全网络是最简单的非平面网络。最多包含四个节点的完全网络是平面网络。通过一些尝试绘制这些网络，你将发现一个没有交叉边的绘图。此外，任何包含至少五个节点的完全图的网络都不是平面网络。平面网络在理论中很重要，因为它们相对简单，但在实际应用中不太常见。

还有更多内容...

除了网络类中的方法外，NetworkX 包中还有其他多个常用函数，可以用来访问网络中节点和边的属性。例如，nx.get_node_attributes可以获取网络中每个节点的指定属性。

为网络生成邻接矩阵

分析图的一个强大工具是邻接矩阵，如果存在从节点 到节点 的边，则矩阵中的值为 ，否则为 0。对于大多数网络，邻接矩阵通常是稀疏的（大多数值为 0）。对于无向网络，邻接矩阵还会是对称的（）。可以与网络关联的其他矩阵也有很多。我们将在这个示例的*还有更多内容...*部分简要讨论这些矩阵。

在这个示例中，我们将生成一个网络的邻接矩阵，并学习如何从这个矩阵中获取网络的一些基本属性。

准备工作

对于这个示例，我们将需要导入 NetworkX 包，并使用nx别名，同时导入 NumPy 模块，使用np别名。

如何实现...

以下步骤概述了如何为网络生成邻接矩阵，并从该矩阵推导出网络的一些简单属性：

1. 首先，我们将生成一个网络，在整个示例中使用。我们将生成一个具有五个节点和五条边的随机网络，并使用种子确保可重复性：

   G = nx.dense_gnm_random_graph(5, 5, seed=12345)
   

2. 为了生成邻接矩阵，我们可以使用 NetworkX 中的 adjacency_matrix 函数。默认情况下，这将返回一个稀疏矩阵，因此我们还将通过 todense 方法将其转换为一个完整的 NumPy 数组进行演示：

   matrix = nx.adjacency_matrix(G).todense()
   

   print(matrix)
   

   # [[0 0 1 0 0]
   

   #  [0 0 1 1 0]
   

   #  [1 1 0 0 1]
   

   #  [0 1 0 0 1]
   

   #  [0 0 1 1 0]]
   

3. 对邻接矩阵进行 次幂运算，可以得到从一个节点到另一个节点的路径数量，路径长度为 ：

   paths_len_4 = np.linalg.matrix_power(matrix, 4)
   

   print(paths_len_4)
   

   # [[ 3 5  0  0 5]
   

   #  [ 5 9  0  0 9]
   

   #  [ 0 0 13 10 0]
   

   #  [ 0 0 10  8 0]
   

   #  [ 5 9  0  0 9]]
   

步骤 2 中的邻接矩阵和步骤 3 中的四次幂都是对称矩阵。此外，注意到 paths_len_4 中非零条目的位置与邻接矩阵中的零位置相对应。这是因为存在两个不同的节点组，而奇数长度的路径在这两个组之间交换，偶数长度的路径则返回到起始组。

它是如何工作的……

dense_gnm_random_graph 函数生成一个（密集）随机网络，网络从所有具有 个节点和 条边的网络族中均匀选择。在这个示例中， 和 。dense 前缀表示该函数使用的算法应比替代的 gnm_random_graph 在边与节点比值较大的密集网络中速度更快。

网络的邻接矩阵很容易生成，特别是在图较小时，使用稀疏形式尤为方便。对于较大的网络，这可能是一个昂贵的操作，因此可能不太实际，特别是当你像本示例中那样将其转换为完整矩阵时。通常你不需要这样做，因为我们可以简单地使用 adjacency_matrix 函数生成的稀疏矩阵，并使用 SciPy sparse 模块中的稀疏线性代数工具。

矩阵的幂次提供了关于给定长度路径数量的信息。通过矩阵乘法的定义，可以很容易地看出这一点。记住，当两个节点之间存在边（路径长度为 1）时，邻接矩阵中的条目为 1。

还有更多……

网络的邻接矩阵的特征值提供了关于网络结构的额外信息，比如网络的着色数界限。（有关着色网络的更多信息，请参见 着色网络 这一部分。）有一个单独的函数用于计算邻接矩阵的特征值。例如，adjacency_spectrum 函数可以生成网络邻接矩阵的特征值。涉及到与网络相关的矩阵特征值的方法通常被称为 谱方法。

网络中还有其他矩阵，如关联矩阵和拉普拉斯矩阵。网络的关联矩阵是一个矩阵，其中是节点数，是边的数量。如果节点出现在边中，则该矩阵的项为 1，否则为 0。网络的拉普拉斯矩阵定义为矩阵，其中是包含网络中节点度数的对角矩阵，是网络的邻接矩阵。这两种矩阵对于网络分析非常有用。

创建有向加权网络

像之前食谱中描述的简单网络，对于描述边的方向不重要且边权相等的网络非常有用。实际上，大多数网络携带额外的信息，如权重或方向。

在本食谱中，我们将创建一个有向加权网络，并探索此类网络的一些基本属性。

准备工作

对于本食谱，我们需要使用 NetworkX 包，通过别名nx导入（如同往常一样），Matplotlib 的pyplot模块以plt导入，以及 NumPy 包以np导入。

如何操作…

以下步骤概述了如何创建一个带权重的有向网络，以及如何探索我们在之前的食谱中讨论的一些属性和技术：

1. 要创建一个有向网络，我们可以使用 NetworkX 中的DiGraph类，而不是简单的Graph类：

   G = nx.DiGraph()
   

2. 如同往常一样，我们必须使用add_node或add_nodes_from方法向网络中添加节点：

   G.add_nodes_from(range(5))
   

3. 要添加带权重的边，我们可以使用add_edge方法并提供weight关键字参数，或者使用add_weighted_edges_from方法：

   G.add_edge(0, 1, weight=1.0)
   

   G.add_weighted_edges_from([
   

       (1, 2, 0.5), (1, 3, 2.0), (2, 3, 0.3), (3, 2, 0.3),
   

       (2, 4, 1.2), (3, 4, 0.8)
   

   ])
   

4. 接下来，我们必须绘制网络，并用箭头表示每条边的方向。我们还必须为此图提供位置：

   fig, ax = plt.subplots()
   

   pos = {0: (-1, 0), 1: (0, 0), 2: (1, 1), 3: (1, -1),
   

       4:(2, 0)}
   

   nx.draw(G, ax=ax, pos=pos, with_labels=True)
   

   ax.set_title("Weighted, directed network")
   

结果图如下所示：

图 5.3 – 一个带权有向网络

1. 有向矩阵的邻接矩阵与简单网络创建方式相同，但结果矩阵将不再是对称的：

   adj_mat = nx.adjacency_matrix(G).todense()
   

   print(adj_mat)
   

   # [[0\. 1\. 0\. 0\. 0\. ]
   

   # [0\. 0\. 0.5 2\. 0\. ]
   

   # [0\. 0\. 0\. 0.3 1.2]
   

   # [0\. 0\. 0.3 0\. 0.8]
   

   # [0\. 0\. 0\. 0\. 0\. ]]
   

与其说是两个给定节点之间的边数，邻接矩阵包含的是这些节点之间边的权重之和。

它是如何工作的…

DiGraph类表示一个有向网络，其中添加边时节点的顺序很重要。在本食谱中，我们添加了两条边来连接节点 2 和节点 3，每个方向一条。在简单网络（Graph类）中，添加第二条边不会增加另一条边。然而，在有向网络（DiGraph类）中，添加边时节点的顺序决定了边的方向。

加权边没有特别之处，除了附加在边上的weight属性。（可以通过关键字参数将任意数据附加到网络中的边或节点。）add_weighted_edges_from方法只是将对应的权重值（元组中的第三个值）添加到相关的边上。权重可以添加到任何网络中的任何边，不仅仅是本食谱中展示的有向网络。

draw例程在绘制有向网络时会自动为边添加箭头。可以通过传递arrows=False关键字参数来关闭此行为。有向或加权网络的邻接矩阵也不同于简单网络的邻接矩阵。在有向网络中，矩阵通常不是对称的，因为边可能存在于一个方向上但不在另一个方向上。对于加权网络，矩阵的条目可能不是 1 或 0，而是对应边的权重。

还有更多内容...

加权网络出现在许多应用中，例如描述交通网络中的距离或速度。你还可以使用网络来检查流量通过网络的情况，方法是为网络中的边提供容量（作为权重或其他属性）。NetworkX 提供了多个分析网络流量的工具，例如通过nx.maximum_flow例程找到网络中的最大流量。

有向网络向网络中添加了方向信息。许多实际应用产生了具有单向边的网络，例如工业过程或供应链网络中的那些边。这种额外的方向信息对许多处理网络的算法有影响，正如我们在本章中将看到的那样。

在网络中寻找最短路径

网络出现的一个常见问题是找到两个节点之间的最短路径——或者更准确地说，找到最高奖励的路径。例如，这可能是两个城市之间的最短距离，其中节点代表城市，边代表连接城市对的道路。在这种情况下，边的权重就是它们的长度。

在本食谱中，我们将找到一个加权网络中两个节点之间的最短路径。

准备工作

对于本食谱，我们将按常规导入 NetworkX 包，并使用nx别名，导入 Matplotlib 的pyplot模块作为plt，以及从 NumPy 导入随机数生成器对象：

from numpy.random import default_rng
rng = default_rng(12345) # seed for reproducibility

如何操作...

按照以下步骤在网络中找到两个节点之间的最短路径：

1. 首先，我们将使用gnm_random_graph和一个seed来创建一个随机网络，作为本演示的基础：

   G = nx.gnm_random_graph(10, 17, seed=12345)
   

2. 接下来，我们将使用圆形布局绘制网络，以便查看节点之间的连接方式：

   fig, ax = plt.subplots()
   

   nx.draw_circular(G, ax=ax, with_labels=True)
   

   ax.set_title("Random network for shortest path finding")
   

结果图形可以在下图中看到。这里，我们可以看到节点 7 和节点 9 之间没有直接的边：

图 5.4 – 一个随机生成的网络，包含 10 个节点和 17 条边。

1. 现在，我们需要给每条边添加一个权重，以便在最短路径方面有些路线比其他路线更具优势：

   for u, v in G.edges:
   

       G.edges[u, v]["weight"] = rng.integers(5, 15)
   

2. 接下来，我们将使用nx.shortest_path例程计算从节点 7 到节点 9 的最短路径：

   path = nx.shortest_path(G, 7, 9, weight="weight")
   

   print(path)
   

   # [7, 5, 2, 9]
   

3. 我们可以使用nx.shortest_path_来找到这条最短路径的长度。

4. length例程：

   length = nx.shortest_path_length(G, 7, 9,
   

       weight="weight")
   

   print("Length", length)
   

   # Length 32
   

这里路径的长度是沿最短路径遍历的边的权重总和。如果网络没有权重，那么它将等于沿路径遍历的边的数量。

它是如何工作的……

shortest_path例程计算每一对节点之间的最短路径。或者，当提供源节点和目标节点时（这就是我们在本示例中所做的），它计算两个指定节点之间的最短路径。我们提供了可选的weight关键字参数，这使得算法根据边的权重属性来寻找最短路径。这个参数改变了最短的定义，默认情况下是最少边数。

查找两个节点之间最短路径的默认算法是 Dijkstra 算法，它是计算机科学和数学课程中的基础算法。它是一个通用的算法，但效率并不是特别高。其他的路线寻找算法包括 A算法。通过使用 A算法并加入额外的启发式信息来指导节点选择，可以提高效率。

还有更多……

有许多算法可以用来寻找网络中两个节点之间的最短路径。也有一些变体可以寻找最大权重路径。

关于网络中路径查找，有几个相关的问题，比如旅行推销员问题和路径检查问题。在旅行推销员问题中，我们需要找到一个环路（一个从同一节点出发并返回的路径），它遍历网络中的每个节点，且总权重最小（或最大）。在路径检查问题中，我们寻找遍历网络中每条边并返回起点的最短环路（按权重计算）。旅行推销员问题是已知的 NP 难题，但路径检查问题可以在多项式时间内解决。

图论中的一个著名问题是哥尼斯堡的桥梁问题，它要求在网络中找到一条路径，使得每条边恰好经过一次。正如欧拉所证明的那样，哥尼斯堡桥问题中找到这样一条路径是不可能的。这样一条恰好经过每条边一次的路径被称为 欧拉回路。一个网络如果包含欧拉回路，则称为 欧拉网络。只有当每个节点的度数都是偶数时，网络才是欧拉网络。哥尼斯堡桥问题的网络表示可以在下图中看到。图中的边代表河流上的不同桥梁，而节点代表不同的陆地块。我们可以看到，所有四个节点的度数都是奇数，这意味着不存在一条路径能恰好经过每一条边一次：

图 5.5 – 代表哥尼斯堡桥问题的网络

这些边代表连接不同陆地块的桥梁，这些陆地块由节点表示。

网络中的聚类度量

网络中有各种量度用于衡量网络的特征。例如，节点的聚类系数衡量了附近节点之间的互联性（此处，附近指的是通过边连接的节点）。实际上，它衡量了相邻节点形成一个完整网络或 团体 的接近程度。

节点的聚类系数衡量了邻近节点通过边连接的比例；也就是说，两个相邻节点与给定节点形成一个三角形。我们计算三角形的数量，并将其除以根据节点的度数可以形成的三角形的总数。从数值上讲，在一个简单的无权网络中，节点的聚类系数 由以下公式给出：

这里， 是在 处的三角形数量，分母是 处可能形成的三角形总数。如果 的度数（即从 出发的边数）为 0 或 1，则我们将 设为 0。

在这个示例中，我们将学习如何计算网络中节点的聚类系数。

准备工作

对于这个示例，我们将需要导入 nx 别名下的 NetworkX 包和作为 plt 导入的 Matplotlib pyplot 模块。

如何进行...

以下步骤展示了如何计算网络中节点的聚类系数：

1. 首先，我们需要创建一个样本网络来进行操作：

   G = nx.Graph()
   

   complete_part = nx.complete_graph(4)
   

   cycle_part = nx.cycle_graph(range(4, 9))
   

   G.update(complete_part)
   

   G.update(cycle_part)
   

   G.add_edges_from([(0, 8), (3, 4)])
   

2. 接下来，我们必须绘制网络，以便我们能够比较我们将要计算的聚类系数。这样，我们可以看到这些节点在网络中的表现：

   fig, ax = plt.subplots()
   

   nx.draw_circular(G, ax=ax, with_labels=True)
   

   ax.set_title("Network with different clustering behavior")
   

结果图表可以在下图中看到：

图 5.6 – 测试聚类的示例网络

1. 现在，我们可以使用 nx.clustering 例程计算网络中节点的聚类系数：

   cluster_coeffs = nx.clustering(G)
   

2. nx.clustering 例程的输出是网络中各个节点的字典。因此，我们可以按如下方式打印一些选定的节点：

   for i in [0, 2, 6]:
   

       print(f"Node {i}, clustering {cluster_coeffs[i]}")
   

   # Node 0, clustering 0.5
   

   # Node 2, clustering 1.0
   

   # Node 6, clustering 0
   

3. 所有网络节点的平均聚类系数可以通过 nx.average_clustering 例程计算：

   av_clustering = nx.average_clustering(G)
   

   print(av_clustering)
   

   # 0.3333333333333333
   

这个平均聚类系数表明，平均而言，节点大约拥有 1/3 的所有可能连接。

它是如何工作的...

节点的聚类系数衡量的是该节点的邻域有多接近于一个完整的网络（所有节点都相互连接）。在这个例子中，我们计算了三个不同的值：节点 0 的聚类系数为 0.5，节点 2 的聚类系数为 1.0，节点 6 的聚类系数为 0。这意味着与节点 2 相连的节点构成了一个完整的网络，这是因为我们设计网络时就是这样设计的。（节点 0 到 4 按设计形成了一个完整的网络。）节点 6 的邻域离完整网络非常远，因为它的两个邻居之间没有任何相互连接的边。

平均聚类值是对网络中所有节点的聚类系数的简单平均。它与全局聚类系数（通过 NetworkX 中的 nx.transitivity 例程计算）略有不同，但它确实能给我们一个网络整体接近完整网络的程度的概念。全局聚类系数衡量的是三角形的数量与三元组的数量之间的比率——三元组是由至少两条边连接的三个节点组成的集合——覆盖整个网络。

全局聚类和平均聚类之间的差异相当微妙。全局聚类系数衡量的是整个网络的聚类程度，而平均聚类系数衡量的是网络在局部的聚类程度。这个差异在风车型网络中最为明显，该网络由一个中心节点和一个偶数个节点围成的圆圈组成。所有节点都连接到中心，但圆圈上的节点只按交替的模式连接。外部节点的局部聚类系数为 1，而中心节点的局部聚类系数为 ，其中 表示连接中心节点的三角形数量。然而，全局聚类系数是 。

还有更多...

聚类系数与网络中的 团体（cliques）有关。团体是一个子网络，其中所有节点都通过边连接。网络理论中的一个重要问题是寻找网络中的最大团体，这通常是一个非常困难的问题（这里，最大意味着 不能更大）。

给网络着色

网络在调度问题中也非常有用，在这些问题中，你需要将活动安排到不同的时间段，以避免冲突。例如，我们可以使用网络来安排课程，以确保选择不同课程的学生不会同时上两门课。在这种情况下，节点将表示不同的课程，边将表示学生同时选修两门课程。我们用来解决这类问题的过程叫做网络着色。这个过程涉及为网络中的节点分配最少的颜色，以确保没有两个相邻节点有相同的颜色。

在这个配方中，我们将学习如何着色一个网络以解决简单的调度问题。

准备工作

对于这个配方，我们需要导入 NetworkX 包并使用nx别名，同时导入 Matplotlib 的pyplot模块并命名为plt。

如何操作...

按照以下步骤来解决网络着色问题：

1. 首先，我们将创建一个示例网络，以便在这个配方中使用：

   G = nx.complete_graph(3)
   

   G.add_nodes_from(range(3, 7))
   

   G.add_edges_from([
   

       (2, 3), (2, 4), (2, 6), (0, 3), (0, 6), (1, 6),
   

       (1, 5), (2, 5), (4, 5) ])
   

2. 接下来，我们将绘制网络，以便在生成着色时理解它。为此，我们将使用draw_circular例程：

   fig, ax = plt.subplots()
   

   nx.draw_circular(G, ax=ax, with_labels=True)
   

   ax.set_title("Scheduling network")
   

结果图可以在下图中看到：

](tos-cn-i-73owjymdk6/68abac282a4f4e7898339ad15dccd94a)

图 5.7 – 简单调度问题的示例网络

1. 我们将使用nx.greedy_color例程生成着色：

   coloring = nx.greedy_color(G)
   

   print("Coloring", coloring)
   

   # Coloring {2: 0, 0: 1, 1: 2, 5: 1, 6: 3, 3: 2, 4: 2}
   

2. 为了查看实际使用的颜色，我们将从coloring字典中生成一组值：

   different_colors = set(coloring.values())
   

   print("Different colors", different_colors)
   

   # Different colors {0, 1, 2, 3}
   

请注意，着色中的颜色数量不能更少，因为节点 0、1、2 和 6 构成了一个完全网络——这些节点彼此相连，因此每个节点都需要一个单独的颜色。

它是如何工作的...

nx.greedy_color例程使用几种可能的策略之一来着色网络。默认情况下，它按度数从大到小的顺序工作。在我们的例子中，它首先将颜色 0 分配给度数为 6 的节点 2，然后将颜色 1 分配给度数为 4 的节点 0，以此类推。在这个序列中，为每个节点选择第一个可用的颜色。这不一定是最有效的着色算法。

任何网络都可以通过为每个节点分配不同的颜色来着色，但在大多数情况下，只需要更少的颜色。在这个配方中，网络有七个节点，但只需要四种颜色。所需的最小颜色数量称为网络的色度数。

我们在这里描述的问题是节点着色问题。还有一个相关的问题叫做边着色。我们可以通过考虑一个网络，其中的节点是原始网络的边，并且当两个节点之间存在共同的原始网络节点时，在这两个节点之间添加一条边，从而将边着色问题转化为节点着色问题。

还有更多...

网络着色问题有多种变体。一个这样的变体是列表着色问题，其中我们为网络中的每个节点从预定义的颜色列表中选择一种颜色进行着色。这个问题比一般的着色问题更为复杂。

一般的着色问题有一些令人惊讶的结果。例如，每个平面网络最多可以用四种不同的颜色进行着色。这是图论中的一个著名定理，称为四色定理，由 Appel 和 Haken 在 1977 年证明。这个定理表明，每个平面图的色数不超过 4。

查找最小生成树和支配集

网络在许多问题中都有应用。通信和分配是两个显而易见的应用领域。例如，我们可能希望找到一种将商品分配到多个城市（节点）的方法，前提是它能覆盖从某个特定点出发的最小距离。对于此类问题，我们需要研究最小生成树和支配集。

在本食谱中，我们将找到网络中的最小生成树和支配集。

准备工作

对于本食谱，我们需要导入nx别名下的 NetworkX 包和 Matplotlib 的pyplot模块，别名为plt。

如何操作...

按照以下步骤为网络查找最小生成树和支配集：

1. 首先，我们将创建一个示例网络来进行分析：

   G = nx.gnm_random_graph(15, 22, seed=12345)
   

2. 接下来，和往常一样，我们将在进行任何分析之前先绘制网络：

   fig, ax = plt.subplots()
   

   pos = nx.circular_layout(G)
   

   nx.draw(G, pos=pos, ax=ax, with_labels=True, style="--")
   

   ax.set_title("Network with minimum spanning tree overlaid")
   

3. 最小生成树可以使用nx.minimum_``spanning_tree例程来计算：

   min_span_tree = nx.minimum_spanning_tree(G)
   

   print(list(min_span_tree.edges))
   

   # [(0, 13), (0, 7), (0, 5), (1, 13), (1, 11),
   

   #    (2, 5), (2, 9), (2, 8), (2, 3), (2, 12),
   

   #    (3, 4), (4, 6), (5, 14), (8, 10)]
   

4. 接下来，我们将把最小生成树的边叠加到图上：

   nx.draw_networkx_edges(min_span_tree, pos=pos,
   

                          ax=ax,width=2.)
   

5. 最后，我们将使用nx.dominating_set例程为网络找到一个支配集——一个网络中每个节点都与集合中的至少一个节点相邻的集合：

   dominating_set = nx.dominating_set(G)
   

   print("Dominating set", dominating_set)
   

   # Dominating set {0, 1, 2, 4, 10, 14}
   

可以在以下图中看到叠加了最小生成树的网络图：

图 5.8 – 显示叠加最小生成树的网络

最小生成树中使用的边是粗体的未断开的线，而原始网络中的边是虚线。尽管最小生成树实际上是一个树，但由于布局的原因这一点稍显模糊，但我们可以轻松地追踪并看到没有任何两个节点连接到同一个父节点的情况。

它是如何工作的...

网络的生成树是包含网络中所有节点的树。最小生成树是包含最少边的生成树——或具有最低总权重的生成树。最小生成树在网络分发问题中非常有用。找到最小生成树的一个简单算法是选择边（如果网络是加权的，则优先选择最小权重的边），以避免形成环路，直到无法再避免为止。

网络的支配集是一个顶点集合，其中网络中的每个节点都与支配集中的至少一个节点相邻。支配集在通信网络中有应用。我们通常希望找到最小的支配集，但这在计算上是困难的。测试是否存在一个比给定大小更小的支配集是 NP 完全问题。然而，对于某些类型的图，存在一些高效的算法来找到最小的支配集。非正式地说，问题在于，一旦你找到了一个最小大小支配集的候选集，你还必须验证是否没有比它更小的支配集。如果你事先不知道所有可能的支配集，这会变得非常困难。

进一步阅读

关于图论，有几本经典的著作，包括 Bollobás 和 Diestel 的书：

• Diestel, R., 2010. 图论。第 3 版。柏林：Springer。

• Bollobás, B., 2010. 现代图论。纽约，NY：Springer。

第六章：数据和统计学的应用

对于需要分析数据的人来说，Python 最吸引人的特点之一是其庞大的数据处理和分析包生态系统，以及活跃的 Python 数据科学社区。Python 易于使用，同时提供了非常强大、快速的库，使得即使是相对初学的程序员也能够快速、轻松地处理大量数据。许多数据科学包和工具的核心是 pandas 库。pandas 提供了两种基于 NumPy 数组的数据容器类型，并且对标签（不仅仅是简单的整数）有很好的支持。这些数据容器使得处理大型数据集变得非常简单。

数据和统计学是现代世界中无处不在的一部分。Python 在试图理解每天产生的大量数据中处于领先地位，通常，这一切都从 pandas 开始——Python 用于处理数据的基础库。首先，我们将了解一些使用 pandas 处理数据的基本技巧。然后，我们将讨论统计学的基础知识，这将为我们提供通过观察一个小样本来理解整个群体的系统方法。

本章中，我们将学习如何利用 Python 和 pandas 处理大规模数据集并进行统计检验。

本章包括以下内容：

• 创建 Series 和 DataFrame 对象

• 从 DataFrame 加载和存储数据

• 在 DataFrame 中操作数据

• 从 DataFrame 中绘制数据

• 从 DataFrame 中获取描述性统计数据

• 使用抽样理解总体

• 在 DataFrame 中对分组数据进行操作

• 使用 t 检验测试假设

• 使用 ANOVA 检验假设

• 对非参数数据进行假设检验

• 使用 Bokeh 创建交互式图表

什么是统计学？

统计学是通过数学——特别是概率——理论来系统地研究数据。统计学有两个主要方面。第一个方面是数据总结。在这里，我们找到描述一组数据的数值，包括数据的中心（均值或中位数）和分布（标准差或方差）等特征。这些值被称为描述性统计。我们在做的事情是拟合一个概率分布，用来描述某一特征在总体中出现的可能性。在这里，总体指的是某一特征的所有测量值的完整集合——例如，地球上所有活着的人的身高。

统计学的第二个——也是可以说是更重要的——方面是推断。在这里，我们通过计算该总体相对较小样本的数值来估计描述总体数据的分布。我们不仅尝试估计总体的分布，还试图量化我们的近似值有多准确。通常这以置信区间的形式呈现。置信区间是一个值的范围，我们可以有信心地认为真值位于该范围内，基于我们观察到的数据。我们通常会为估计值提供 95%或 99%的置信区间。

推断还包括测试两个或多个样本数据集是否来自同一总体。这就是假设检验的领域。在这里，我们比较两个数据集的可能分布，以确定它们是否可能是相同的。许多假设检验要求数据符合正态分布，或者更可能的是我们可以应用中心极限定理。这些检验有时被描述为参数检验，包括 t 检验和方差分析（ANOVA）。然而，如果你的数据不符合足够的正态性，以至于无法应用中心极限定理，那么一些检验不需要假设正态性。这些检验被称为非参数检验。

技术要求

对于本章，我们将主要使用 pandas 库进行数据处理，它提供了类似 R 的数据结构，例如用于存储、组织和操作数据的Series和DataFrame对象。在本章的最后一节食谱中，我们还将使用 Bokeh 数据可视化库。这些库可以通过你喜欢的包管理器安装，例如pip：

python3.10 -m pip install pandas bokeh

我们还将使用 NumPy 和 SciPy 包。

本章的代码可以在本书的 GitHub 仓库中的Chapter 06文件夹找到：github.com/PacktPublishing/Applying-Math-with-Python-2nd-Edition/tree/main/Chapter%2006。

创建 Series 和 DataFrame 对象

在 Python 中，大多数数据处理是通过 pandas 库完成的，pandas 基于 NumPy 提供类似 R 的结构来存储数据。这些结构允许使用字符串或其他 Python 对象（除了整数）进行轻松的行列索引。一旦数据被加载到 pandas DataFrame或Series中，它就可以像在电子表格中一样轻松操作。这使得 Python 与 pandas 结合时，成为处理和分析数据的强大工具。

在本食谱中，我们将看到如何创建新的 pandas Series和DataFrame对象，并访问其中的项目。

准备工作

对于本食谱，我们将使用以下命令将 pandas 库导入为pd：

import pandas as pd

NumPy 包是np。我们还必须从 NumPy 中创建一个（带种子的）随机数生成器，如下所示：

from numpy.random import default_rng
rng = default_rng(12345)

如何操作...

以下步骤概述了如何创建保存数据的Series和DataFrame对象：

1. 首先，创建我们将存储在Series和DataFrame对象中的随机数据：

   diff_data = rng.normal(0, 1, size=100)
   

   cumulative = diff_data.cumsum()
   

2. 接下来，创建一个包含diff_data的Series对象。我们将打印Series以查看数据：

   data_series = pd.Series(diff_data)
   

   print(data_series)
   

3. 现在，创建一个包含两列的DataFrame对象：

   data_frame = pd.DataFrame({
   

       "diffs": data_series,
   

       "cumulative": cumulative
   

   })
   

4. 打印DataFrame对象以查看它所包含的数据：

   print(data_frame)
   

打印出的对象如下；左边是Series对象，右边是DataFrame对象：

                                     diffs  cumulative
0    -1.423825                0  -1.423825   -1.423825
1     1.263728                1   1.263728   -0.160097
2    -0.870662                2  -0.870662   -1.030758
3    -0.259173                3  -0.259173   -1.289932
4    -0.075343                4  -0.075343   -1.365275
        ...                  ..       ...         ...
95   -0.061905               95 -0.061905   -1.107210
96   -0.359480               96 -0.359480   -1.466690
97   -0.748644               97 -0.748644   -2.215334
98   -0.965479               98 -0.965479   -3.180813
99    0.360035               99  0.360035   -2.820778
Length: 100, dtype: float64  [100 rows x 2 columns]

如预期的那样，Series和DataFrame都包含 100 行。由于系列中的数据是单一类型的——这通过它只是一个 NumPy 数组来保证——数据类型显示为float64。DataFrame有两列，通常这两列可能有不同的数据类型（尽管在这里，它们都有float64）。

它是如何工作的……

pandas 包提供了Series和DataFrame类，它们的功能和特性与 R 语言中的类似对象相对应。Series用于存储一维数据，如时间序列数据，而DataFrame用于存储多维数据；你可以将DataFrame对象视为一个“电子表格”。

Series与简单的 NumPy ndarray的区别在于Series对其项的索引方式。NumPy 数组是通过整数索引的，这也是Series对象的默认索引方式。然而，Series可以通过任何可哈希的 Python 对象进行索引，包括字符串和datetime对象。这使得Series在存储时间序列数据时非常有用。Series可以通过多种方式创建。在这个示例中，我们使用了 NumPy 数组，但也可以使用任何 Python 可迭代对象，如列表。

DataFrame对象中的每一列都是包含行的系列，就像传统的数据库或电子表格一样。在这个示例中，DataFrame对象通过字典的键来构建时，列会被赋予标签。

DataFrame和Series对象在打印时会创建一个数据摘要。这包括列名、行数和列数，以及框架（序列）的前五行和后五行。这有助于快速获取对象及其数据分布的概况。

还有更多内容……

Series对象的单独行（记录）可以通过常规的索引符号访问，只需提供相应的索引。我们也可以通过特殊的iloc属性对象按其数字位置访问行。这允许我们像在 Python 列表或 NumPy 数组中一样，通过数字（整数）索引访问行。

DataFrame对象中的列可以通过常规的索引符号访问，提供列名即可。这样得到的结果是一个包含所选列数据的Series对象。DataFrame还提供了两个可以用于访问数据的属性。loc属性通过索引访问单独的行，无论该对象是什么。而iloc属性通过数字索引访问行，就像在Series对象中一样。

您可以向 loc 提供选择标准（或仅使用对象的索引表示法）来选择数据。这包括单个标签、标签列表、标签切片或布尔数组（大小适当）。iloc 选择方法接受类似的标准。

从 Series 或 DataFrame 对象中选择数据的方式有很多，超出了我们在此所描述的简单方法。例如，我们可以使用 at 属性访问对象中指定行（和列）的单个值。

有时候，pandas 的 Series 或 DataFrame 不能充分描述数据，因为它们本身是低维的。xarray 包基于 pandas 接口并提供对带标签的多维数组（即 NumPy 数组）的支持。在 第十章 的 从 NetCDF 文件加载和存储数据 示例中，我们将学习关于 xarray 的内容。更多关于 xarray 的信息可以在文档中找到：docs.xarray.dev/en/stable/index.html。

另见

pandas 文档包含了创建和索引 DataFrame 或 Series 对象的不同方式的详细描述：pandas.pydata.org/docs/user_guide/indexing.html。

从 DataFrame 加载和存储数据

在 Python 会话中从原始数据创建 DataFrame 对象是相当不寻常的。实际上，数据通常来自外部来源，例如现有的电子表格或 CSV 文件、数据库或 API 接口。因此，pandas 提供了许多用于加载和存储数据到文件的工具。pandas 默认支持从 CSV、Excel（xls 或 xlsx）、JSON、SQL、Parquet 和 Google BigQuery 加载和存储数据。这使得将数据导入 pandas 并使用 Python 操作和分析这些数据变得非常简单。

在本示例中，我们将学习如何加载和存储 CSV 文件中的数据。加载和存储其他文件格式的数据的操作方法将类似。

准备工作

对于这个示例，我们需要导入 pandas 包并使用 pd 别名，同时导入 NumPy 库并命名为 np。我们还必须使用以下命令从 NumPy 创建一个默认的随机数生成器：

from numpy.random import default_rng
rng = default_rng(12345) # seed for example

让我们学习如何存储数据并从 DataFrame 中加载数据。

如何操作...

按照以下步骤将数据存储到文件中，然后将数据加载回 Python：

1. 首先，我们将使用随机数据创建一个示例的 DataFrame 对象。然后，我们将打印这个 DataFrame 对象，以便我们可以将其与稍后读取的数据进行比较：

   diffs = rng.normal(0, 1, size=100)
   

   cumulative = diffs.cumsum()
   

   data_frame = pd.DataFrame({
   

       "diffs": diffs, 
   

       "cumulative": cumulative
   

   })
   

   print(data_frame)
   

2. 我们将通过在 DataFrame 对象上使用 to_csv 方法，将数据存储到 sample.csv 文件中。我们将使用 index=False 关键字参数，以确保索引不被存储在 CSV 文件中：

   data_frame.to_csv("sample.csv", index=False)
   

3. 现在，我们可以使用 pandas 中的read_csv方法将sample.csv文件读入一个新的DataFrame对象。我们将打印这个对象，以展示结果：

   df = pd.read_csv("sample.csv", index_col=False)
   

   print(df)
   

两个打印出来的 DataFrame 并排显示。第一步中的DataFrame对象在左侧，第三步中的DataFrame对象在右侧：

    diffs      cumulative          diffs       cumulative
0  -1.423825   -1.423825        0  -1.423825   -1.423825
1   1.263728   -0.160097        1   1.263728   -0.160097
2  -0.870662   -1.030758        2  -0.870662   -1.030758
3  -0.259173   -1.289932        3  -0.259173   -1.289932
4  -0.075343   -1.365275        4  -0.075343   -1.365275
..         ...            ...        ..         ...            ...
95 -0.061905   -1.107210        95 -0.061905   -1.107210
96 -0.359480   -1.466690        96 -0.359480   -1.466690
97 -0.748644   -2.215334        97 -0.748644   -2.215334
98 -0.965479   -3.180813        98 -0.965479   -3.180813
99  0.360035   -2.820778        99  0.360035   -2.820778
[100 rows x 2 columns]           [100 rows x 2 columns]

如我们从行中所看到的，这两个 DataFrame 是完全相同的。

它是如何工作的...

本教程的核心是 pandas 中的read_csv方法。该方法接受路径或类似文件的对象作为参数，并将文件内容读取为 CSV 数据。我们可以通过sep关键字参数定制分隔符，默认是逗号（,）。还可以定制列头以及每列的类型。

DataFrame或Series中的to_csv方法将内容存储到 CSV 文件中。我们在这里使用了index关键字参数，以确保索引不会被打印到文件中。这意味着 pandas 会根据 CSV 文件中的行号推断索引。如果数据是按整数索引的，这种行为是可取的，但如果数据是按时间或日期索引的，情况可能不同。我们还可以使用这个关键字参数指定 CSV 文件中作为索引的列。

另请参阅

请参阅 pandas 文档，了解支持的文件格式：pandas.pydata.org/docs/reference/io.html。

操作 DataFrame 中的数据

一旦我们有了DataFrame中的数据，我们通常需要对数据进行一些简单的转换或筛选，然后才能进行分析。这可能包括，例如，筛选出缺失数据的行，或者对单独的列应用函数。

在这个教程中，我们将学习如何对DataFrame对象进行一些基本操作，以准备数据进行分析。

准备工作

对于这个教程，我们需要导入 pandas 库，别名为pd，导入 NumPy 库，别名为np，并通过以下命令创建一个默认的 NumPy 随机数生成器对象：

from numpy.random import default_rng
rng = default_rng(12345)

让我们学习如何对DataFrame中的数据进行一些简单的操作。

如何做到...

以下步骤展示了如何对 pandas DataFrame进行一些基本的筛选和操作：

1. 首先，我们将使用随机数据创建一个示例DataFrame：

   three = rng.uniform(-0.2, 1.0, size=100)
   

   three[three < 0] = np.nan
   

   data_frame = pd.DataFrame({
   

       "one": rng.random(size=100),
   

       "two": rng.normal(0, 1, size=100).cumsum(),
   

       "three": three
   

   })
   

2. 接下来，我们将从现有的列生成一个新列。这个新列将存储True，如果对应的"one"列的值大于0.5，否则存储False：

   data_frame["four"] = data_frame["one"] > 0.5
   

3. 现在，让我们创建一个新函数，我们将应用到我们的DataFrame上。这个函数将把行 "two" 的值乘以行 "one" 和0.5中的最大值（还有更简洁的方式来写这个函数）：

   def transform_function(row):
   

       if row["four"]:
   

           return 0.5*row["two"]
   

           return row["one"]*row["two"]
   

4. 现在，我们将把之前定义的函数应用到 DataFrame 的每一行，生成一个新列。我们还将打印更新后的 DataFrame，方便后续对比：

   data_frame["five"] = data_frame.apply(
   

       transform_function, axis=1)
   

   print(data_frame)
   

5. 最后，我们需要过滤掉包含非数字（NaN）值的 DataFrame 行。我们将打印出处理后的 DataFrame：

   df = data_frame.dropna()
   

   print(df)
   

步骤 4 中 print 命令的输出如下：

         one       two     three   four      five
0   0.168629  1.215005  0.072803  False  0.204885
1   0.240144  1.431064  0.180110  False  0.343662
2   0.780008  0.466240  0.756839   True  0.233120
3   0.203768 -0.090367  0.611506  False -0.018414
4   0.552051 -2.388755  0.269331   True -1.194377
..         ...         ...         ...     ...         ...
95  0.437305  2.262327  0.254499  False  0.989326
96  0.143115  1.624691  0.131136  False  0.232517
97  0.422742  2.769187  0.959325  False  1.170652
98  0.764412  1.128285         NaN   True  0.564142
99  0.413983 -0.185719  0.290481  False -0.076885
[100 rows x 5 columns]

在第 98 行可见一个 NaN 值。正如预期的那样，我们总共有 100 行数据和 5 列数据。现在，我们可以将其与 步骤 6 中 print 命令的输出进行对比：

         one       two     three   four      five
0   0.168629  1.215005  0.072803  False  0.204885
1   0.240144  1.431064  0.180110  False  0.343662
2   0.780008  0.466240  0.756839   True  0.233120
3   0.203768 -0.090367  0.611506  False -0.018414
4   0.552051 -2.388755  0.269331   True -1.194377
..         ...         ...         ...     ...         ...
94  0.475131  3.065343  0.330151  False  1.456440
95  0.437305  2.262327  0.254499  False  0.989326
96  0.143115  1.624691  0.131136  False  0.232517
97  0.422742  2.769187  0.959325  False  1.170652
99  0.413983 -0.185719  0.290481  False -0.076885
[88 rows x 5 columns]

正如预期的那样，行数减少了 12 行，因为我们已删除了所有包含 NaN 值的行。（请注意，第 98 行在第 3 列中不再包含 NaN。）

工作原理...

可以通过简单地将新列分配给新的列索引来向现有的 DataFrame 中添加新列。然而，这里需要特别注意。在某些情况下，pandas 会创建 DataFrame 对象的“视图”而不是复制对象，在这种情况下，将新列分配给 DataFrame 可能不会产生预期的效果。有关这一点，可以参阅 pandas 文档（pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy）。

pandas 的 Series 对象（即 DataFrame 中的列）支持丰富的比较运算符，比如等于、小于和大于（在本示例中，我们使用了大于运算符）。这些比较运算符返回一个 Series，其中包含与比较结果为真或假的位置相对应的布尔值。接着，可以利用这个布尔值 Series 来索引原始的 Series，从而获取比较结果为真的行。在这个示例中，我们只是将这个布尔值 Series 添加到原始的 DataFrame 中。

apply 方法接受一个函数（或其他可调用的函数），并将其应用于 DataFrame 对象中的每一列。在这个示例中，我们希望将函数应用于每一行，因此使用了 axis=1 关键字参数，将函数应用于 DataFrame 对象中的每一行。无论哪种方式，函数都会获得一个按行（列）索引的 Series 对象。我们还对每一行应用了一个函数，该函数返回一个基于每一行数据计算出的值。实际上，如果 DataFrame 对象包含大量行，这样的操作会非常慢。如果可能的话，你应该操作整个列，使用专为 NumPy 数组设计的函数，这样效率更高。尤其是在对 DataFrame 中列的值进行简单算术运算时，这一点尤为重要。就像 NumPy 数组一样，Series 对象也实现了标准的算术运算，这可以大大提高处理大规模 DataFrame 的速度。

在本教程的最后一步，我们使用了dropna方法来快速选择仅包含非 NaN 值的 DataFrame 行。pandas 使用 NaN 来表示DataFrame中的缺失数据，因此此方法选择那些不包含缺失值的行。该方法返回原始DataFrame对象的视图，但也可以通过传递inplace=True关键字参数来修改原始DataFrame。正如本教程所示，这大致相当于使用索引数组选择行，数组中包含布尔值。

注意

在直接修改原始数据时，您应始终小心，因为可能无法返回该数据以便以后重复分析。如果确实需要直接修改数据，您应确保数据已备份，或者确保修改不会删除以后可能需要的数据。

还有更多...

大多数 pandas 例程都能合理处理缺失数据（NaN）。然而，如果确实需要在DataFrame中删除或替换缺失数据，则有几种方法可以做到。在本教程中，我们使用dropna方法简单地删除了缺失数据的行。作为替代方案，我们可以使用fillna方法将所有缺失值填充为特定值，或者使用interpolate方法根据周围的值插值缺失值。

更一般来说，我们可以使用replace方法将特定的（非 NaN）值替换为其他值。此方法可以同时处理数值和字符串类型的值，包括使用正则表达式进行模式匹配。

DataFrame类有许多有用的方法。我们这里只介绍了最基本的方法，但还有另外两个方法也值得提及。这些方法是agg方法和merge方法。

agg方法对给定轴上的一个或多个操作的结果进行聚合。这允许我们通过应用聚合函数快速生成每列（或每行）的汇总信息。输出是一个DataFrame，其中包含应用函数的名称作为行，所选轴（例如列标签）的标签作为列。

merge方法在两个 DataFrame 之间执行类似 SQL 的连接操作。这将生成一个新的DataFrame，包含连接操作的结果。可以通过how关键字参数指定要执行的连接类型，默认是inner。连接操作的列名或索引应通过on关键字参数传递——如果两个DataFrame对象包含相同的键——或者通过left_on和right_on传递。以下是一个简单的 DataFrame 连接示例：

rng = default_rng(12345)
df1 = pd.DataFrame({
    "label": rng.choice(["A", "B", "C"], size=5),
    "data1": rng.standard_normal(size=5)
})
df2 = pd.DataFrame({
    "label": rng.choice(["A", "B", "C", "D"], size=4),
    "data2": rng.standard_normal(size=4)
})
df3 = df1.merge(df2, how="inner", on="label")

这将生成一个包含label、data1和data2的DataFrame，这些数据对应于df1和df2中具有相同标签的行。我们打印这三个 DataFrame 来查看结果：

>>> print(df1)                          >>> print(df2)
  label      data1                        label      data2
0      C -0.259173                     0      D  2.347410
1      A -0.075343                     1      A  0.968497
2      C -0.740885                     2      C -0.759387
3      A -1.367793                     3      C  0.902198
4      A  0.648893
>>> df3
  label      data1      data2
0      C -0.259173 -0.759387
1      C -0.259173  0.902198
2      C -0.740885 -0.759387
3      C -0.740885  0.902198
4      A -0.075343  0.968497
5      A -1.367793  0.968497
6      A  0.648893  0.968497

在这里，你可以看到每一组来自df1和df2的data1和data2值组合，配有相应的标签，都会有一行出现在df3中。此外，df2中标签为D的行未被使用，因为df1中没有标签为D的行。

绘制来自 DataFrame 的数据

正如许多数学问题一样，找到一种可视化问题和所有信息的方法的第一步通常是制定策略。对于数据相关问题，这通常意味着生成数据的图表，并通过目视检查来寻找趋势、模式和潜在结构。由于这是一个常见的操作，pandas 提供了一个快速简单的接口，用于直接从Series或DataFrame绘制数据，默认情况下使用 Matplotlib 作为后台。

在本教程中，我们将学习如何直接从DataFrame或Series绘制数据，以理解其背后的趋势和结构。

准备工作

对于本教程，我们将需要导入pandas库（作为pd），NumPy库（作为np），Matplotlib的pyplot模块（作为plt），并使用以下命令创建默认的随机数生成器实例：

from numpy.random import default_rng
rng = default_rng(12345)

如何操作...

按照以下步骤，使用随机数据创建一个简单的DataFrame并绘制其包含的数据图表：

1. 使用随机数据创建一个示例DataFrame：

   diffs = rng.standard_normal(size=100)
   

   walk = diffs.cumsum()
   

   df = pd.DataFrame({
   

       "diffs": diffs,
   

       "walk": walk
   

   })
   

2. 接下来，我们需要创建一个空白图形，包含两个子图以备绘制：

   fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True)
   

3. 我们需要将walk列绘制为标准的折线图。这可以通过在Series（列）对象上使用plot方法，且不需要额外的参数来完成。我们将通过传递ax=ax1关键字参数强制在ax1上绘制：

   df["walk"].plot(ax=ax1, title="Random walk", color="k")
   

   ax1.set_xlabel("Index")
   

   ax1.set_ylabel("Value")
   

4. 现在，我们需要通过将kind="hist"关键字参数传递给plot方法，绘制diffs列的直方图：

   df["diffs"].plot(kind="hist", ax=ax2, 
   

       title="Histogram of diffs", color="k", alpha=0.6)
   

   ax2.set_xlabel("Difference")
   

结果图形如下所示：

图 6.1 – 来自 DataFrame 的行走值和差异的直方图

在这里，我们可以看到，差异的直方图接近标准正态分布（均值 0，方差 1）。随机游走图展示了差异的累计和，并且围绕 0 上下摆动（相当对称）。

它是如何工作的...

Series（或DataFrame）上的plot方法是一种快速的方式，用于将其包含的数据与行索引进行绘制。kind关键字参数用于控制绘制的图表类型，默认情况下是折线图。还有许多绘图类型的选项，包括bar用于垂直条形图，barh用于水平条形图，hist用于直方图（在本教程中也有提到），box用于箱线图，scatter用于散点图。还有其他几个关键字参数可用于自定义生成的图表。在本教程中，我们还提供了title关键字参数，为每个子图添加标题。

由于我们希望将两个图表并排放置在同一图形上，使用之前已创建的子图，因此我们使用 ax 关键字参数将相应的坐标轴句柄传递给绘图例程。即使你让 plot 方法构建图形，你可能仍然需要使用 plt.show 例程来显示图形并应用某些设置。

还有更多...

我们可以通过 pandas 接口生成几种常见的图表类型。除了本食谱中提到的图表外，还包括散点图、条形图（水平条形图和垂直条形图）、区域图、饼图和箱型图。plot 方法还接受各种关键字参数，用于自定义图表的外观。

从 DataFrame 获取描述性统计数据

描述性统计数据，或称汇总统计数据，是与一组数据相关的简单值，例如均值、中位数、标准差、最小值、最大值和四分位数。这些值从不同角度描述了数据集的位置和分布。均值和中位数是数据中心（位置）的度量，而其他值衡量数据从均值和中位数的分散程度。这些统计数据对于理解数据集至关重要，并且构成了许多分析技术的基础。

在本食谱中，我们将学习如何为 DataFrame 中的每一列生成描述性统计数据。

准备工作

对于本食谱，我们需要导入 pandas 包并命名为 pd，导入 NumPy 包并命名为 np，导入 Matplotlib 的 pyplot 模块并命名为 plt，并通过以下命令创建默认的随机数生成器：

from numpy.random import default_rng
rng = default_rng(12345)

如何实现...

以下步骤展示了如何为 DataFrame 中的每一列生成描述性统计数据：

1. 首先，我们将创建一些示例数据供我们分析：

   uniform = rng.uniform(1, 5, size=100)
   

   normal = rng.normal(1, 2.5, size=100)
   

   bimodal = np.concatenate([rng.normal(0, 1, size=50), 
   

       rng.normal(6, 1, size=50)])
   

   df = pd.DataFrame({
   

       "uniform": uniform, 
   

       "normal": normal, 
   

       "bimodal": bimodal
   

   })
   

2. 接下来，我们将绘制数据的直方图，以便了解 DataFrame 对象中数据的分布：

   fig, (ax1, ax2, ax3) = plt.subplots(1, 3,
   

                                       tight_layout=True)
   

   df["uniform"].plot(kind="hist",
   

       title="Uniform", ax=ax1, color="k", alpha=0.6)
   

   df["normal"].plot(kind="hist",
   

       title="Normal", ax=ax2, color="k", alpha=0.6)
   

3. 为了更好地观察 bimodal 数据的分布，我们将直方图的箱数更改为 20：

   df["bimodal"].plot(kind="hist", title="Bimodal",
   

       ax=ax3, bins=20, color="k", alpha=0.6)
   

4. pandas DataFrame 对象有一个方法，可以获取每一列的几种常见描述性统计数据。describe 方法会创建一个新的 DataFrame，该 DataFrame 的列标题与原始对象相同，每一行包含不同的描述性统计数据：

   descriptive = df.describe()
   

5. 我们还必须计算 峰度 并将其添加到我们刚刚获得的新 DataFrame 对象中。我们还必须将描述性统计数据打印到控制台，以查看这些值是什么：

   descriptive.loc["kurtosis"] = df.kurtosis()
   

   print(descriptive)
   

   #             uniform      normal     bimodal
   

   # count     100.000000 100.000000 100.000000
   

   # mean         2.813878   1.087146   2.977682
   

   # std           1.093795   2.435806   3.102760
   

   # min           1.020089  -5.806040  -2.298388
   

   # 25%           1.966120  -0.498995   0.069838
   

   # 50%           2.599687   1.162897   3.100215
   

   # 75%           3.674468   2.904759   5.877905
   

   # max           4.891319   6.375775   8.471313
   

   # kurtosis  -1.055983   0.061679  -1.604305
   

6. 最后，我们必须在直方图中添加垂直线，以显示每种情况下均值的值：

   uniform_mean = descriptive.loc["mean", "uniform"]
   

   normal_mean = descriptive.loc["mean", "normal"]
   

   bimodal_mean = descriptive.loc["mean", "bimodal"]
   

   ax1.vlines(uniform_mean, 0, 20, "k")
   

   ax2.vlines(normal_mean, 0, 25, "k")
   

   ax3.vlines(bimodal_mean, 0, 20, "k")
   

结果的直方图如下所示：

图 6.2 – 三组数据的直方图及其均值标示

在这里，我们可以看到平均值对于正态分布的数据（中间）是中心的，但对于均匀分布的数据（左侧），分布的“质量”稍微偏向于左侧的较低值。对于双峰日（右侧），平均线恰好位于质量的两个组成部分之间。

工作原理...

describe 方法返回一个包含以下数据描述统计的 DataFrame：计数、平均值、标准差、最小值、25% 四分位数、中位数（50% 四分位数）、75% 四分位数和最大值。计数、最小值和最大值都很容易理解。平均值和中位数是数据的两种不同的平均值，大致代表数据的中心值。平均值的定义应该很熟悉，即所有值的总和除以值的数量。我们可以用以下公式表示这个数量：

这里， 值代表数据值， 是数值（计数）的数量。在这里，我们还采用了用横线表示平均值的常见符号。中位数是当所有数据排序时的“中间值”（如果值的数量为奇数，则取两个中间值的平均值）。25% 和 75% 处的四分位数同样定义，但是取有序值的 25%或 75%处的值。你也可以将最小值视为 0%四分位数，最大值视为 100%四分位数。

标准差是数据与平均值之间的分布范围的度量，与统计学中经常提到的另一个量方差有关。方差是标准差的平方，定义如下：

在这里你可能也会看到分数中出现的，这是在从样本估计总体参数时对偏差的校正。我们将在下一个示例中讨论总体参数及其估计。标准差、方差、四分位数、最大值和最小值描述了数据的分布。例如，如果最大值为 5，最小值为 0，25%四分位数为 2，75%四分位数为 4，则表示大部分（实际上至少 50%的值）的数据集中在 2 和 4 之间。

峰度是衡量数据在分布的“尾部”（远离均值的地方）集中程度的指标。这种指标不像我们在本教程中讨论的其他量那样常见，但它确实在某些分析中会出现。我们在这里主要是为了演示如何计算在describe方法返回的DataFrame对象中未出现的汇总统计量——这里是kurtosis方法。对于计算均值（mean）、标准差（std）以及describe方法中的其他量，当然也有单独的方法。

注意

当 pandas 计算本教程中描述的各项量时，它会自动忽略任何由 NaN 表示的“缺失值”。这也会反映在描述性统计中报告的计数中。

还有更多...

我们在统计中包含的第三个数据集说明了观察数据的重要性，以确保我们计算的数值合理。事实上，我们计算出的均值约为2.9，但从直方图来看，显然大多数数据远离这个值。我们应该始终检查我们计算的汇总统计量是否能准确地反映样本数据的特征。仅仅引用均值可能会对样本产生不准确的描述。

使用抽样理解总体

统计学中的一个核心问题是进行估计——并量化这些估计的准确性——即在只有一个小（随机）样本的情况下，估计整个总体的分布。一个经典的例子是估计一个国家所有人群的平均身高，方法是测量一个随机选取的人的身高。这类问题在真实的总体分布无法实测时特别有趣，通常我们指的总体均值就是指整体人口的均值。在这种情况下，我们必须依赖统计学知识和一个（通常要小得多的）随机选取的样本来估计总体的真实均值和标准差，同时也量化我们估计的准确性。正是后者导致了公众对统计学的困惑、误解和错误解读。

在本教程中，我们将学习如何估算总体均值，并为这些估算提供置信区间。

准备工作

对于本教程，我们需要导入pandas包并使用别名pd，导入 Python 标准库中的math模块，并导入 SciPy 的stats模块，使用以下命令：

from scipy import stats

让我们学习如何使用 SciPy 的统计函数构建置信区间。

如何操作...

在接下来的步骤中，我们将基于随机选取的 20 人的样本，估算英国男性的平均身高：

1. 首先，我们必须将样本数据加载到 pandas Series中：

   sample_data = pd.Series(
   

       [172.3, 171.3, 164.7, 162.9, 172.5, 176.3, 174.8,
   

       171.9,176.8, 167.8, 164.5, 179.7, 157.8, 170.6,
   

       189.9, 185., 172.7, 165.5, 174.5, 171.5]
   

   )
   

2. 接下来，我们必须计算样本的均值和标准差：

   sample_mean = sample_data.mean()
   

   sample_std = sample_data.std()
   

   print(f"Mean {sample_mean}, st. dev {sample_std}")
   

   # Mean 172.15, st. dev 7.473778724383846
   

3. 然后，我们必须计算标准误差，如下所示：

   N = sample_data.count()
   

   std_err = sample_std/math.sqrt(N)
   

4. 我们必须从学生 t 分布中计算出我们希望的置信度所需的临界值：

   cv_95, cv_99 = stats.t.ppf([0.975, 0.995], df=N-1)
   

5. 现在，我们可以使用以下代码计算真实总体均值的 95%和 99%置信区间：

   pm_95 = cv_95*std_err
   

   conf_interval_95 = [sample_mean - pm_95,
   

       sample_mean + pm_95]
   

   pm_99 = cv_99*std_err
   

   conf_interval_99 = [sample_mean - pm_99,
   

       sample_mean + pm_99]
   

   print("95% confidence", conf_interval_95)
   

   # 95% confidence [168.65216388659374, 175.64783611340627]
   

   print("99% confidence", conf_interval_99)
   

   # 99% confidence [167.36884119608774, 176.93115880391227]
   

它是如何工作的...

参数估计的关键是正态分布，我们在第四章《与随机性和概率打交道》中讨论过它。如果我们找到标准正态分布的临界值！，使得一个标准正态分布随机数小于该值的概率为 97.5%，那么该随机数位于！和！之间的概率为 95%（每尾 2.5%）。这个临界值！的结果是 1.96，保留到小数点后两位。也就是说，我们可以 95%确定一个标准正态分布随机数的值介于！和！之间。类似地，99%置信度的临界值是 2.58（保留到小数点后两位）。

如果我们的样本是“大的”，我们可以调用 SciPy stats模块中的stats.t.ppf例程。

学生 t 分布与正态分布相关，但它有一个参数——自由度——改变分布的形态。随着自由度的增加，学生 t 分布将越来越像正态分布。你认为分布足够相似的临界点取决于你的应用和数据。一个常见的经验法则是样本量为 30 时，可以调用中心极限定理，直接使用正态分布，但这并不一定是一个好规则。你在基于样本做出推论时要非常小心，尤其是当样本相对于总体来说非常小的时候。

一旦我们得到了临界值，就可以通过将临界值乘以样本的标准误差，然后将结果加减到样本均值上来计算真实总体均值的置信区间。标准误差是给定样本量的样本均值分布的扩展的近似值，相对于真实总体均值。这就是为什么我们使用标准误差来提供我们对总体均值估计的置信区间。当我们将标准误差乘以从学生 t 分布（在此情况下）获得的临界值时，我们得到了观察到的样本均值与真实总体均值之间的最大差异估计值，该差异值对应于给定的置信水平。

在本节中，这意味着我们有 95% 的信心认为英国男性的平均身高位于 168.7 cm 和 175.6 cm 之间，且我们有 99% 的信心认为英国男性的平均身高位于 167.4 cm 和 176.9 cm 之间。我们的样本来自一个均值为 175.3 cm，标准差为 7.2 cm 的总体。这个真实的均值（175.3 cm）确实位于我们的两个置信区间内，但只是勉强。

另请参见

有一个非常有用的包叫做 uncertainties，用于处理涉及带有不确定性的值的计算。更多信息请参见 第十章 中的 计算中的不确定性处理 配方，位于 提高生产力 一章。

在 DataFrame 中对分组数据执行操作

pandas DataFrame 的一个重要特性是可以通过特定列中的值对数据进行分组。例如，我们可以按生产线 ID 和班次 ID 对装配线数据进行分组。能够以符合人体工程学的方式对这些分组数据进行操作非常重要，因为数据通常是为了分析而聚合的，但在预处理时需要进行分组。

在本节中，我们将学习如何对 DataFrame 中的分组数据执行操作。我们还将借此机会展示如何对（分组的）数据的滚动窗口进行操作。

准备工作

对于本节，我们需要导入 NumPy 库并命名为 np，导入 Matplotlib 的 pyplot 接口并命名为 plt，还需要导入 pandas 库并命名为 pd。我们还需要如下创建一个默认随机数生成器的实例：

rng = np.random.default_rng(12345)

在开始之前，我们还需要设置 Matplotlib 绘图设置，以便在本节中更改绘图样式。我们将改变在同一坐标轴上生成多个图形时，循环使用的绘图样式机制，这通常会导致不同的颜色。为了实现这一点，我们将其更改为生成黑色线条并使用不同的线条样式：

from matplotlib.rcsetup import cycler
plt.rc("axes", prop_cycle=cycler(
    c=["k"]*3, ls=["-", "--", "-."]))

现在，让我们学习如何使用 pandas DataFrame 的分组功能。

如何操作...

按照以下步骤，学习如何在 pandas DataFrame 中对分组数据执行操作：

1. 首先，我们需要在 DataFrame 中生成一些示例数据。对于这个示例，我们将生成两个标签列和一个数值数据列：

   labels1 = rng.choice(["A", "B", "C"], size=50)
   

   labels2 = rng.choice([1, 2], size=50)
   

   data = rng.normal(0.0, 2.0, size=50)
   

   df = pd.DataFrame({"label1": labels1, "label2": labels2, "data": data})
   

2. 接下来，让我们添加一个新列，该列由按第一个标签 "label1" 分组的 "data" 列的累积和组成：

   df[“first_group”] = df.groupby(“label1”)[“data”].cumsum()
   

   print(df.head())
   

现在，df 的前五行如下：

  label1  label2      data  first_group
0      C       2  0.867309     0.867309
1      A       2  0.554967     0.554967
2      C       1  1.060505     1.927814
3      A       1  1.073442     1.628409
4      A       1  1.236700     2.865109

在这里，我们可以看到 "first_group" 列包含了 "label1" 列中每个标签的累积和。例如，第 0 行和第 1 行的和只是 "data" 列中的值。第 2 行的新条目是第 0 行和第 2 行数据的总和，因为这两行是第一个带有标签 "C" 的行。

1. 现在，让我们对 "label1" 和 "label2" 列同时进行分组：

   grouped = df.groupby(["label1", "label2"])
   

2. 现在，我们可以使用 transform 和 rolling 方法在分组数据上计算每组内连续条目的滚动均值：

   df["second_group"] = grouped["data"].transform(lambda d:
   

       d.rolling(2, min_periods=1).mean())
   

   print(df.head())
   

   print(df[df["label1"]=="C"].head())
   

前五行打印结果如下：

  label1  label2      data  first_group  second_group
0      C       2  0.867309     0.867309      0.867309
1      A       2  0.554967     0.554967      0.554967
2      C       1  1.060505     1.927814      1.060505
3      A       1  1.073442     1.628409      1.073442
4      A       1  1.236700     2.865109      1.155071

与之前一样，前几行都表示不同的组，因此 "second_group" 列中的值与 "data" 列中的相应值相同。第 4 行的值是第 3 行和第 4 行数据值的平均值。接下来的五行是标签为 C 的行：

  label1  label2      data  first_group  second_group
0      C       2  0.867309     0.867309      0.867309
2      C       1  1.060505     1.927814      1.060505
5      C       1 -1.590035     0.337779     -0.264765
7      C       1 -3.205403    -2.867624     -2.397719
8      C       1  0.533598    -2.334027     -1.335903

在这里，我们可以更清楚地看到滚动平均值和累计和。除了第一行外，其他行的标签都相同。

1. 最后，让我们绘制由 "label1" 列分组的 "first_group" 列的值：

   fig, ax = plt.subplots()
   

   df.groupby("label1")["first_group"].plot(ax=ax)
   

   ax.set(title="Grouped data cumulative sums",     xlabel="Index", ylabel="value")
   

   ax.legend()
   

结果图如 图 6.3 所示：

图 6.3 – 按 label1 组绘制的累计和图

在这里，我们可以看到每个组在图上都产生了一条不同的线。这是一种快速简便的方法，可以从 DataFrame 中绘制分组数据的图。（请记住，在 准备工作 部分中我们更改了默认的样式周期，使得图表风格在页面上更加鲜明。）

它是如何工作的...

groupby 方法为 DataFrame 创建了一个代理，其索引是从请求的列中生成的。我们可以在这个代理对象上执行操作。在这种情况下，我们使用 cumsum 方法来生成每个组内 "data" 列的累计和。我们可以使用这种方法以相同的方式生成分组数据的汇总统计信息。这对于数据探索非常有用。

在这个食谱的第二部分，我们根据两个不同的标签列进行了分组，并对每个组计算了滚动平均值（窗口长度为 2）。请注意，我们使用 transform 方法“包装”了这个计算，而不是直接在分组的 DataFrame 上调用 rolling。这样做是为了确保结果具有正确的索引，可以放回到 df 中。否则，mean 的输出将继承分组索引，我们将无法将结果放回到 df 中。我们在 rolling 上使用了 min_periods 可选参数，以确保所有行都有值。否则，窗口大小之前的行将被分配为 NaN。

这个食谱的最后一部分使用了 plot 函数，绘制了按 "label1" 分组的数据。这是一种快速简便的方法，可以从同一个 DataFrame 对象中绘制多个数据流。不幸的是，在这种情况下，定制绘图稍显困难，尽管可以通过使用 Matplotlib 中的 rcparams 设置来完成。

使用 t 检验进行假设检验

统计学中最常见的任务之一是测试一个关于正态分布人口均值的假设有效性，前提是你已经从该人口中收集了样本数据。例如，在质量控制中，我们可能希望检验一个生产厂的板材厚度是否为 2 毫米。为了进行此测试，我们可以随机选择样本板材并测量其厚度，从而获得样本数据。然后，我们可以使用 stats 模块计算 t 统计量和 值。如果 值低于 0.05，则我们以 5% 的显著性水平（95% 的置信度）接受零假设。如果 值大于 0.05，则我们必须拒绝零假设，支持替代假设。

在这个食谱中，我们将学习如何使用 t 检验来检验给定样本的假设人口均值是否有效。

准备工作

对于这个食谱，我们需要导入 pandas 包并命名为 pd，同时导入 SciPy stats 模块，使用以下命令：

from scipy import stats

让我们学习如何使用 SciPy stats 模块进行 t 检验。

如何操作...

按照以下步骤使用 t 检验来检验给定一些样本数据的假设人口均值的有效性：

1. 首先，我们必须将数据加载到 pandas Series 中：

   sample = pd.Series([
   

       2.4, 2.4, 2.9, 2.6, 1.8, 2.7, 2.6, 2.4, 2.8, 
   

       2.4, 2.4, 2.4, 2.7, 2.7, 2.3, 2.4, 2.4, 3.2, 
   

       2.9, 2.5, 2.5, 3.2, 2\. , 2.3, 3\. , 1.5, 3.1,
   

       2.5, 2.2, 2.5, 2.1,1.8, 3.1, 2.4, 3\. , 2.5,
   

       2.7, 2.1, 2.3, 2.2, 2.5, 2.6, 2.5, 2.8, 2.5,
   

       2.9, 2.1, 2.8, 2.1, 2.3
   

   ])
   

2. 现在，让我们设置假设的人口均值和我们将进行检验的显著性水平：

   mu0 = 2.0
   

   significance = 0.05
   

3. 接下来，我们将使用 SciPy stats 模块中的 ttest_1samp 例程来生成 t 统计量和 值：

   t_statistic, p_value = stats.ttest_1samp(sample, mu0)
   

   print(f"t stat: {t_statistic}, p value: {p_value}")
   

   # t stat: 9.752368720068665, p value: 4.596949515944238e-13
   

4. 最后，让我们测试 值是否小于我们选择的显著性水平：

   if p_value <= significance:
   

       print("Reject H0 in favour of H1: mu != 2.0")
   

   else:
   

       print("Accept H0: mu = 2.0")
   

   # Reject H0 in favour of H1: mu != 2.0
   

我们可以有 95% 的信心得出结论：样本数据来源的人口均值不等于 2（由于样本中大多数数值大于 2，这并不令人惊讶）。鉴于这里的 值如此之小，我们可以非常有信心地认为这一点是正确的。

它是如何工作的...

t 统计量通过以下公式计算：

这里， 是假设的均值（来自零假设）， 是样本均值， 是样本标准差， 是样本的大小。t 统计量是观察到的样本均值与假设人口均值 之间差异的估算值，经过标准误差标准化。如果假设人口呈正态分布，t 统计量将遵循 t 分布，具有 自由度。通过查看 t 统计量在相应的学生 t 分布中的位置，我们可以大致了解观察到的样本均值来自假设均值人口的可能性。这个值用 来表示。

值是观察到的样本均值比假设总体均值更极端的概率，前提是假设总体均值等于。如果值小于我们选择的显著性值，则不能期望真实的总体均值是我们假设的值。在这种情况下，我们接受备择假设，认为真实的总体均值不等于。

还有更多……

本配方中演示的检验是 t 检验的最基本应用。在这里，我们将样本均值与假设的总体均值进行比较，以决定整个总体均值是否合理地是这个假设值。更一般来说，我们可以使用 t 检验来比较两个独立的总体，前提是从每个总体中抽取了样本，使用双样本 t 检验，或者比较数据以某种方式配对的总体，使用配对 t 检验。这使得 t 检验成为统计学家重要的工具。

显著性和置信度是统计学中经常出现的两个概念。统计学上显著的结果有很高的正确概率。在许多情况下，我们认为任何错误概率低于某个阈值（通常是 5%或 1%）的结果都是统计显著的。置信度是我们对结果有多确定的量化表示。结果的置信度是 1 减去显著性。

不幸的是，结果的显著性往往被误用或误解。说一个结果在 5%的显著性水平上是显著的，意味着我们有 5%的概率错误地接受了零假设。换句话说，如果我们对来自总体的 20 个其他样本重复相同的测试，我们期望至少其中一个样本会给出相反的结果。然而，这并不意味着其中一个样本必然会如此。

高显著性表明我们更确定我们得出的结论是正确的，但这并不能保证情况确实如此。这个配方中的结果就是证据；我们使用的样本来自一个均值为2.5，标准差为0.35的总体。（样本在创建后进行了某些四舍五入处理，这稍微改变了分布。）这并不是说我们的分析是错误的，或者我们从样本中得出的结论不是正确的。

重要的是要记住，t 检验仅在基础总体遵循正态分布，或至少大致符合时才有效。如果不是这种情况，您可能需要使用非参数检验。我们将在非参数 数据假设检验的章节中讨论这一点。

使用 ANOVA 进行假设检验

假设我们设计了一个实验，测试两个新流程与当前流程的差异，并且我们想要测试这些新流程的结果是否与当前流程不同。在这种情况下，我们可以使用方差分析（ANOVA）来帮助我们判断三组结果的均值是否存在差异（为此，我们需要假设每个样本来自具有共同方差的正态分布）。

在本食谱中，我们将学习如何使用 ANOVA 来比较多个样本之间的差异。

准备就绪

对于这个食谱，我们需要使用 SciPy 的 stats 模块。我们还需要通过以下命令创建一个默认的随机数生成器实例：

from numpy.random import default_rng
rng = default_rng(12345)

如何操作…

按照以下步骤执行（单因素）ANOVA 检验，以测试三个不同流程之间的差异：

1. 首先，我们将创建一些样本数据，然后进行分析：

   current = rng.normal(4.0, 2.0, size=40)
   

   process_a = rng.normal(6.2, 2.0, size=25)
   

   process_b = rng.normal(4.5, 2.0, size=64)
   

2. 接下来，我们将设置我们检验的显著性水平：

   significance = 0.05
   

3. 然后，我们将使用 SciPy stats 模块中的 f_oneway 例程生成 F 统计量和 值：

   F_stat, p_value = stats.f_oneway(
   

       current, process_a, process_b)
   

   print(f"F stat: {F_stat}, p value: {p_value}")
   

   # F stat: 9.949052026027028, p value: 9.732322721019206e-05
   

4. 现在，我们必须测试 值是否足够小，以判断我们是否应接受或拒绝原假设，即所有均值相等：

   if p_value <= significance:
   

       print("Reject H0: there is a difference between means")
   

   else:
   

       print("Accept H0: all means equal")
   

   # Reject H0: there is a difference between means
   

这里， 值非常小（在 级别），因此差异在 95% 置信度下显著（即，），并且在 99% 置信度下也显著（）。

工作原理…

ANOVA 是一种强大的技术，可以同时比较多个样本之间的差异。它通过比较样本之间的变异性与总体变异性来工作。ANOVA 在比较三个或更多样本时尤其强大，因为它不会因进行多次检验而产生累积误差。不幸的是，如果 ANOVA 检测到并非所有均值相等，则无法通过检验信息确定哪些样本之间存在显著差异。为此，您需要使用额外的检验来找出差异。

f_oneway SciPy stats 模块中的例程执行单因素方差分析（ANOVA）检验——ANOVA 中生成的检验统计量遵循 F 分布。再次强调， 值是检验中最关键的信息。如果 值小于我们预先设定的显著性水平（在本食谱中为 5%），我们接受原假设；否则拒绝原假设。

还有更多…

ANOVA 方法非常灵活。我们在这里介绍的一元 ANOVA 检验是最简单的情况，因为这里只有一个因子需要检验。二元 ANOVA 检验可以用来检验两个不同因子之间的差异。例如，在药物临床试验中，我们可以进行对照组的测试，同时还衡量性别（例如）对结果的影响。不幸的是，SciPy 在 stats 模块中没有执行二元 ANOVA 的例程。你需要使用其他包，如 statsmodels 包。

如前所述，ANOVA 只能检测是否存在差异。如果存在显著差异，它无法检测这些差异发生在哪里。例如，我们可以使用 Durnett 检验来测试其他样本的均值是否与对照样本的均值不同，或者使用 Tukey 范围检验来测试每个组的均值与其他所有组的均值。

非参数数据的假设检验

t 检验和 ANOVA 的一个主要缺点是：所采样的总体必须遵循正态分布。在许多应用中，这并不太具限制性，因为许多实际的总体值遵循正态分布，或者某些定理（例如中心极限定理）允许我们分析相关数据。然而，实际上并非所有的总体值都以合理的方式遵循正态分布。对于这些（幸运的是，较为少见的）情况，我们需要使用一些替代的检验统计量，作为 t 检验和 ANOVA 的替代方法。

在本方案中，我们将使用 Wilcoxon 秩和检验和 Kruskal-Wallis 检验来检验两个（或更多，在后一种情况下）总体之间的差异。

准备工作

对于此方案，我们将需要导入 pandas 包作为 pd，SciPy 的 stats 模块，以及使用以下命令创建一个默认的随机数生成器实例：

from numpy.random import default_rng
rng = default_rng(12345)

让我们学习如何使用 SciPy stats 中的非参数假设检验工具。

如何操作…

按照以下步骤比较两个或更多非正态分布的样本：

1. 首先，我们将生成一些样本数据来进行分析：

   sample_A = rng.uniform(2.5, 3.5, size=25)
   

   sample_B = rng.uniform(3.0, 4.4, size=25)
   

   sample_C = rng.uniform(3.1, 4.5, size=25)
   

2. 接下来，我们将设定分析中使用的显著性水平：

   significance = 0.05
   

3. 现在，我们将使用 stats.kruskal 例程生成检验统计量和 值，用于检验原假设，即总体具有相同的中位数值：

   statistic, p_value = stats.kruskal(sample_A, sample_B,
   

       sample_C)
   

   print(f"Statistic: {statistic}, p value: {p_value}")
   

   # Statistic: 40.22214736842102, p value: 1.8444703308682906e-09
   

4. 我们将使用条件语句打印测试结果的陈述：

   if p_value <= significance:
   

       print("There are differences between population  medians")
   

   else:
   

       print("Accept H0: all medians equal")
   

   # There are differences between population medians
   

5. 现在，我们将使用 Wilcoxon 秩和检验获取每对样本比较的 值。这些检验的原假设是它们来自同一分布：

   _, p_A_B = stats.ranksums(sample_A, sample_B)
   

   _, p_A_C = stats.ranksums(sample_A, sample_C)
   

   _, p_B_C = stats.ranksums(sample_B, sample_C)
   

6. 接下来，我们将使用条件语句打印出那些表示显著差异的比较结果：

   if p_A_B <= significance:
   

       print("Significant differences between A and B,
   

           p value", p_A_B)
   

   # Significant differences between A and B, p value
   

   1.0035366080480683e-07
   

   if p_A_C <= significance:
   

       print("Significant differences between A and C,
   

           p value", p_A_C)
   

   # Significant differences between A and C, p value
   

   2.428534673701913e-08
   

   if p_B_C <= significance:
   

       print("Significant differences between B and C,
   

           p value", p_B_C)
   

   else:
   

       print("No significant differences between B and C,
   

           p value", p_B_C)
   

   # No significant differences between B and C, p value
   

   0.3271631660572756
   

这些打印出来的线条显示我们的测试发现，A 组和 B 组、A 组和 C 组之间存在显著差异，但 B 组和 C 组之间没有差异。

它是如何工作的...

我们说数据是非参数的，如果数据来源的总体不符合可以通过少量参数描述的分布。这通常意味着总体不是正态分布的，但范围比这更广。在本示例中，我们从均匀分布中抽样，但这仍然是一个比我们通常需要使用非参数检验时更为结构化的例子。非参数检验可以且应该用于任何我们不确定基础分布的情况。这样做的代价是，检验的效力稍微较弱。

任何（真实的）分析的第一步应该是绘制数据的直方图并直观地检查分布。如果你从一个正态分布的总体中抽取随机样本，你也可以预期样本呈正态分布（在本书中我们已经多次看到这一点）。如果你的样本显示出正态分布的特征性钟形曲线，那么该总体很可能是正态分布的。你还可以使用kind="kde"。如果你仍然不确定总体是否符合正态分布，可以使用统计检验，例如 D'Agostino 的 K-squared 检验或皮尔逊卡方检验来检验正态性。这两种检验被结合成一个名为normaltest的常规方法，用于在 SciPy 的stats模块中进行正态性检验，此外还包括其他几种正态性检验。

Wilcoxon 秩和检验是两样本 t 检验的非参数替代方法。与 t 检验不同，秩和检验不是通过比较样本均值来量化总体是否具有不同的分布。相反，它将样本数据合并并按大小排序。检验统计量由样本中元素最少的那一组的秩和生成。从这里开始，像往常一样，我们生成值，用于检验两个总体是否具有相同分布的原假设。

Kruskal-Wallis 检验是单因素方差分析（ANOVA）检验的非参数替代方法。与秩和检验类似，它利用样本数据的排序生成检验统计量和值，用以检验所有总体的中位数是否相同的原假设。与单因素方差分析一样，我们只能检测所有总体的中位数是否相同，而不能指出差异所在。要做这个，我们还需要使用额外的检验方法。

在本食谱中，我们使用了 Kruskal-Wallis 测试来判断我们的三个样本对应的群体之间是否存在显著差异。我们发现了一个差异，具有非常小的 值。接着我们使用秩和检验来确定群体之间差异的具体位置。在这里，我们发现样本 A 与样本 B 和 C 有显著差异，但 B 与 C 之间没有显著差异。考虑到这些样本的生成方式，这并不令人惊讶。

注意

不幸的是，由于我们在本食谱中使用了多个测试，我们对结论的整体信心并不像我们预期的那样高。我们进行了四次 95%置信度的测试，这意味着我们对结论的总体信心只有大约 81%。这是因为多次测试的误差会累积，从而降低总体信心。为了解决这个问题，我们需要调整每个测试的显著性阈值，使用 Bonferroni 校正（或类似方法）。

使用 Bokeh 创建交互式图表

测试统计和数值推理对于系统地分析数据集非常有用。然而，它们并不能像图表那样给我们一个整体数据的清晰图像。数值值是确定性的，但在统计学中往往难以理解，而图表则能立即直观地展示数据集之间的差异和趋势。因此，存在大量的库，可以用更富创意的方式绘制数据。一个特别有趣的数据绘图库是 Bokeh，它通过利用 JavaScript 库，允许我们在浏览器中创建交互式图表。

在本食谱中，我们将学习如何使用 Bokeh 创建一个可以在浏览器中显示的交互式图表。

准备工作

对于本食谱，我们需要导入 pandas 包作为 pd，导入 NumPy 包作为 np，并使用以下代码构造默认的随机数生成器实例，还需要导入 Bokeh 的 plotting 模块，并用 bk 别名导入：

from bokeh import plotting as bk
from numpy.random import default_rng
rng = default_rng(12345)

如何操作...

这些步骤展示了如何使用 Bokeh 在浏览器中创建交互式图表：

1. 首先，我们需要创建一些样本数据来绘制图表：

   date_range = pd.date_range("2020-01-01", periods=50)
   

   data = rng.normal(0, 3, size=50).cumsum()
   

   series = pd.Series(data, index=date_range)
   

2. 接下来，我们必须使用 output_file 函数指定一个输出文件，用于存储图表的 HTML 代码：

   bk.output_file("sample.html")
   

3. 现在，我们将创建一个新的图形，并设置标题和坐标轴标签，同时将 轴类型设置为 datetime，以便我们的日期索引能够正确显示：

   fig = bk.figure(title="Time series data", 
   

                          x_axis_label="date",
   

                          x_axis_type="datetime",
   

                          y_axis_label="value")
   

4. 我们将数据作为一条线添加到图形中：

   fig.line(date_range, series)
   

5. 最后，我们可以使用 show 函数或 save 函数来保存或更新指定输出文件中的 HTML。在这里我们使用 show 来使图表在浏览器中打开：

   bk.show(fig)
   

Bokeh 图表不是静态对象，应该通过浏览器进行交互。为了进行比较，这里使用matplotlib重新创建了 Bokeh 图表中将显示的数据：

图 6.4 – 使用 Matplotlib 创建的时间序列数据图

图 6.4 – 使用 Matplotlib 创建的时间序列数据图

Bokeh 的真正优势在于其能够将动态、交互式的图表插入到网页和文档中（例如，Jupyter 笔记本），这样读者就可以深入查看绘制数据的细节。

它是如何工作的…

Bokeh 使用 JavaScript 库在浏览器中渲染图表，数据由 Python 后端提供。这样做的优点是，用户可以自行检查图表。例如，我们可以放大以查看可能隐藏的图表细节，或自然地平移数据。此食谱中给出的示例只是使用 Bokeh 可以实现的功能的一个尝试。

figure例程创建一个表示图表的对象，我们向其中添加元素——例如，通过数据点的线——就像我们在 Matplotlib Axes对象中添加图表一样。在这个食谱中，我们创建了一个简单的 HTML 文件，其中包含 JavaScript 代码来渲染数据。每次保存或调用show例程时，这段 HTML 代码会被导出到指定文件。实际上，值越小，我们就越能确信假设的总体均值是正确的。

还有更多…

Bokeh 的功能远远超出了这里的描述。Bokeh 图表可以嵌入到诸如 Jupyter 笔记本之类的文件中，而这些文件也会在浏览器中渲染，或者嵌入到现有的网站中。如果你使用的是 Jupyter 笔记本，应该使用output_notebook例程，而不是output_file例程，将图表直接打印到笔记本中。它具有多种不同的绘图样式，支持在图表之间共享数据（例如，数据可以在一个图表中选择并在其他图表中突出显示），并支持流数据。

进一步阅读

统计学和统计理论方面有大量的教材。以下书籍是本章所涉及的统计学内容的良好参考：

• Mendenhall, W., Beaver, R., 和 Beaver, B.（2006 年），概率与统计导论。第 12 版，（加利福尼亚州贝尔蒙特：Thomson Brooks/Cole）。

• Freedman, D., Pisani, R., 和 Purves, R.（2007 年），统计学。纽约：W.W. Norton。

pandas 文档（pandas.pydata.org/docs/index.html）和以下的 pandas 书籍是学习 pandas 的好参考：

• McKinney, W.（2017 年），Python 数据分析。第二版，（Sebastopol：O'Reilly Media, Inc，美国）。

SciPy 文档 (docs.scipy.org/doc/scipy/tutorial/stats.html) 还包含了关于本章多次使用的统计模块的详细信息。