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

布客飞龙
91 阅读42分钟

Python 数学应用第二版（四）

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

译者：飞龙

协议：CC BY-NC-SA 4.0

第九章：寻找最优解

在这一章中，我们将探讨寻找给定情境下最佳结果的各种方法。这被称为优化，通常涉及最小化或最大化目标函数。目标函数是一个包含一个或多个自变量的函数，返回一个标量值，表示给定参数选择的成本或收益。最小化和最大化函数的问题实际上是等价的，因此我们将在本章中只讨论最小化目标函数。最小化一个函数，，等同于最大化函数。关于这一点的更多细节将在我们讨论第一个食谱时提供。

用于最小化给定函数的算法取决于函数的性质。例如，包含一个或多个变量的简单线性函数与具有多个变量的非线性函数所适用的算法是不同的。线性函数的最小化属于线性规划的范畴，这是一个成熟的理论。线性函数可以通过标准的线性代数技术解决。对于非线性函数，我们通常利用函数的梯度来寻找最小点。我们将讨论几种用于最小化不同类型函数的方法。

寻找单变量函数的最小值和最大值尤其简单，如果已知函数的导数，则可以轻松完成。如果不知道导数，那么可以使用适当食谱中描述的方法。最小化非线性函数食谱中的备注提供了一些额外的细节。

我们还将简要介绍一下博弈论。广义来说，博弈论是围绕决策制定的理论，具有广泛的应用，尤其在经济学等学科中尤为重要。具体来说，我们将讨论如何将简单的双人博弈表示为 Python 中的对象，计算与某些选择相关的收益，并计算这些博弈的纳什均衡。

我们将首先学习如何最小化包含一个或多个变量的线性和非线性函数。然后，我们将继续探讨使用梯度下降法和最小二乘法的曲线拟合方法。最后，我们将通过分析双人博弈和纳什均衡来总结这一章。

本章将涵盖以下食谱：

  • 最小化简单线性函数

  • 最小化非线性函数

  • 在优化中使用梯度下降法

  • 使用最小二乘法对数据进行曲线拟合

  • 分析简单的双人博弈

  • 计算纳什均衡

让我们开始吧！

技术要求

在本章中，我们将需要 NumPy、SciPy 和 Matplotlib 包，和往常一样。对于最后两个食谱，我们还需要 Nashpy 包。可以使用你喜欢的包管理器（如pip）安装这些包：

python3.10 -m pip install numpy scipy matplotlib nashpy

本章节的代码可以在 GitHub 仓库中的Chapter 09文件夹找到，地址为 github.com/PacktPublishing/Applying-Math-with-Python-2nd-Edition/tree/main/Chapter%2009。

最小化一个简单的线性函数

我们在优化中面临的最基本问题是找到一个函数的最小值所在的参数。通常，这个问题会受到一些关于参数可能值的约束，这增加了问题的复杂性。显然，如果我们要最小化的函数本身也很复杂，那么问题的复杂性会进一步增加。因此，我们必须首先考虑线性函数，其形式如下：

要解决这种问题，我们需要将约束转换成计算机能够处理的形式。在这种情况下，我们通常将其转换为线性代数问题（矩阵和向量）。完成这一步后，我们可以利用 NumPy 和 SciPy 中线性代数包的工具来找到我们需要的参数。幸运的是，由于这种问题非常常见，SciPy 已经有处理这种转换和后续求解的例程。

在本例中，我们将使用 SciPy optimize 模块中的例程解决以下约束线性最小化问题：

这将受到以下条件的限制：

让我们看看如何使用 SciPy 的optimize例程来解决这个线性规划问题。

准备就绪

对于这个示例，我们需要导入 NumPy 包（别名为np），Matplotlib 的 pyplot 模块（别名为plt），以及 SciPy 的 optimize 模块。我们还需要从 mpl_toolkits.mplot3d 导入 Axes3D 类，以便进行 3D 绘图：

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

让我们看看如何使用 optimize 模块中的例程来最小化一个有约束的线性系统。

如何实现...

按照以下步骤使用 SciPy 解决一个有约束的线性最小化问题：

  1. 将系统设置为 SciPy 能够识别的形式：

    A = np.array([

        [2, 1], # 2*x0 + x1 <= 6

        [-1, -1] # -x0 - x1 <= -4

    ])

    b = np.array([6, -4])

    x0_bounds = (-3, 14) # -3 <= x0 <= 14

    x1_bounds = (2, 12)      # 2 <= x1 <= 12

    c = np.array([1, 5])

  2. 接下来，我们需要定义一个例程来评估在线性函数值为 时的函数值，该值是一个向量（NumPy 数组）：

    def func(x):

        return np.tensordot(c, x, axes=1)

  3. 然后，我们创建一个新图形，并添加一组 3d 坐标轴，以便在其上绘制函数：

    fig = plt.figure()

    ax = fig.add_subplot(projection="3d")

    ax.set(xlabel="x0", ylabel="x1", zlabel="func")

    ax.set_title("Values in Feasible region")

  4. 接下来，我们创建一个覆盖问题区域的值网格，并在该区域上绘制函数值：

    X0 = np.linspace(*x0_bounds)

    X1 = np.linspace(*x1_bounds)

    x0, x1 = np.meshgrid(X0, X1)

    z = func([x0, x1])

    ax.plot_surface(x0, x1, z, cmap="gray",

        vmax=100.0, alpha=0.3)

  5. 现在，我们绘制函数值平面中对应于临界线2*x0 + x1 == 6的直线，并在我们的图上绘制落在该范围内的值：

    Y = (b[0] - A[0, 0]*X0) / A[0, 1]

    I = np.logical_and(Y >= x1_bounds[0], Y <= x1_bounds[1])

    ax.plot(X0[I], Y[I], func([X0[I], Y[I]]), 

        "k", lw=1.5, alpha=0.6)

  6. 我们对第二条临界线x0 + x1 == -4重复进行绘图步骤：

    Y = (b[1] - A[1, 0]*X0) / A[1, 1]

    I = np.logical_and(Y >= x1_bounds[0], Y <= x1_bounds[1])

    ax.plot(X0[I], Y[I], func([X0[I], Y[I]]), 

        "k", lw=1.5, alpha=0.6)

  7. 接下来，我们对位于两条临界线之间的区域进行着色，这对应于最小化问题的可行区域：

    B = np.tensordot(A, np.array([x0, x1]), axes=1)

    II = np.logical_and(B[0, ...] <= b[0], B[1, ...] <= b[1])

    ax.plot_trisurf(x0[II], x1[II], z[II], 

        color="k", alpha=0.5)

在下图中可以看到函数值在可行区域上的绘图：

图 9.1 – 线性函数的值，并突出显示了可行区域

图 9.1 – 线性函数的值，并突出显示了可行区域

如我们所见，位于这个阴影区域内的最小值出现在两条关键线的交点处。

  1. 接下来，我们使用linprog来解决在步骤 1中创建的约束最小化问题。我们在终端中打印出结果对象：

    res = optimize.linprog(c, A_ub=A, b_ub=b,

        bounds= (x0_bounds, x1_bounds))

    print(res)

  2. 最后，我们将在可行区域上方绘制最小函数值：

    ax.plot([res.x[0]], [res.x[1]], [res.fun], "kx")

更新后的图示可以在下图中看到：

图 9.2 – 最小值绘制在可行区域上

图 9.2 – 最小值绘制在可行区域上

在这里，我们可以看到linprog程序确实找到了最小值，它位于两条关键线的交点处。

它是如何工作的…

约束线性最小化问题在经济场景中非常常见，在这些问题中，你试图在保持其他参数的同时最小化成本。事实上，优化理论中的很多术语都反映了这一点。解决这些问题的一种非常简单的算法叫做单纯形法，它通过一系列阵列操作来找到最小解。从几何角度看，这些操作代表了在单纯形的不同顶点之间变化（我们在这里不会定义单纯形），正是这种变化赋予了算法这个名字。

在我们继续之前，我们将简要概述单纯形法用于求解约束线性优化问题的过程。该问题，呈现给我们的并不是一个矩阵方程问题，而是一个矩阵不等式问题。我们可以通过引入松弛变量来解决这个问题，将不等式转化为等式。例如，通过引入松弛变量，第一条约束不等式可以重写为如下形式：。

只要不为负数，这样就满足了所需的不等式。第二个约束不等式是“大于或等于”型的不等式，我们必须首先将其转化为“小于或等于”型的不等式。我们通过将所有项乘以-1 来实现这一点。这就得到了我们在步骤中定义的A矩阵的第二行。在引入第二个松弛变量后，，我们得到了第二个方程：

从这里，我们可以构建一个矩阵，其列包含两个参数变量的系数， 和 ，以及两个松弛变量， 和 。该矩阵的行表示两个边界方程和目标函数。现在，这个方程组可以通过对该矩阵进行基本的行操作来求解，以得到最小化目标函数的 和 的值。由于求解矩阵方程既简单又快速，这意味着我们可以迅速高效地最小化线性函数。

幸运的是，我们不需要记住如何将我们的不等式系统化简为线性方程组，因为像linprog这样的例程已经为我们做了这件事。我们只需提供边界不等式作为一个矩阵和向量对，包括每个不等式的系数，并提供一个定义目标函数的单独向量。linprog例程会负责公式化并解决最小化问题。

实际上，linprog例程并未使用单纯形法来最小化函数。相反，linprog使用了一个内部点算法，这种算法效率更高。（该方法实际上可以通过提供适当的方法名称和method关键字参数设置为simplex或revised-simplex。在打印的结果输出中，我们可以看到，达到解决方案仅需五次迭代。）该例程返回的结果对象包含在最小值处的参数值，存储在x属性中，最小值处函数的值存储在fun属性中，此外还包含关于求解过程的各种信息。如果方法失败，则status属性将包含一个数值代码，描述方法失败的原因。

在本食谱的第 2 步中，我们创建了一个表示此问题目标函数的函数。该函数接受一个包含函数应评估的参数空间值的单一数组作为输入。在这里，我们使用了 NumPy 的tensordot例程（axes=1）来评估系数向量的点积，，与每个输入，。在这里我们必须非常小心，因为我们传入函数的值在后续步骤中将会是一个 2 × 50 × 50 的数组。普通的矩阵乘法（np.dot）在这种情况下无法给出我们想要的 50 × 50 的数组输出。

在第 5 步和第 6 步中，我们根据以下方程计算了临界线上的点：

然后我们计算了相应的值，以便绘制在目标函数定义的平面上的线条。我们还需要修剪这些值，确保只包含那些在问题中指定的范围内的值。通过在代码中构造标记为 I 的索引数组来完成这一步，该数组由位于边界值内的点组成。

还有更多...

本节介绍了受限最小化问题以及如何使用 SciPy 解决它。然而，相同的方法也可以用来解决受限的最大化问题。这是因为最大化和最小化是对偶的，从某种意义上来说，最大化一个函数与最小化函数并取其负值是一样的。事实上，我们在本节中就利用了这一点，将第二个约束不等式从 ≥ 改为 ≤。

在本节中，我们解决了一个只有两个参数变量的问题，但相同的方法也适用于涉及更多变量的问题（除了绘图步骤）。我们只需要为每个数组增加更多的行和列来处理增加的变量数量——这包括传递给函数的边界元组。当处理大量变量时，该函数还可以在适当的情况下与稀疏矩阵一起使用，以提高效率。

linprog 函数得名于线性规划，用于描述此类问题——找到满足某些矩阵不等式的的值，同时满足其他条件。由于矩阵理论与线性代数之间有着紧密的联系，因此线性规划问题有许多非常快速和高效的技术，而这些技术在非线性上下文中是不可用的。

最小化非线性函数

在上一节中，我们展示了如何最小化一个非常简单的线性函数。不幸的是，大多数函数并非线性，通常也不具备我们希望的良好特性。对于这些非线性函数，我们无法使用为线性问题开发的快速算法，因此我们需要设计出可以在这些更一般化情况下使用的新方法。我们将在这里使用的算法叫做 Nelder-Mead 算法，它是一种健壮的通用方法，用于寻找函数的最小值，并且不依赖于函数的梯度。

在本节中，我们将学习如何使用 Nelder-Mead 单纯形法来最小化一个包含两个变量的非线性函数。

准备工作

在本节中，我们将使用导入的 NumPy 包（命名为 np）、Matplotlib 的 pyplot 模块（命名为 plt）、从 mpl_toolkits.mplot3d 导入的 Axes3D 类来启用 3D 绘图，以及 SciPy 的 optimize 模块：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy import optimize

让我们看看如何使用这些工具来解决非线性优化问题。

如何进行...

以下步骤展示了如何使用 Nelder-Mead 单纯形法找到一般非线性目标函数的最小值：

  1. 定义我们要最小化的目标函数：

    def func(x):

        return ((x[0] - 0.5)**2 + (

            x[1] + 0.5)**2)*np.cos(0.5*x[0]*x[1])

  2. 接下来，创建一个值的网格，以便我们可以在上面绘制目标函数：

    x_r = np.linspace(-1, 1)

    y_r = np.linspace(-2, 2)

    x, y = np.meshgrid(x_r, y_r)

  3. 现在，我们在这个点阵上评估该函数：

    z = func([x, y])

  4. 接下来，我们创建一个包含3d坐标轴对象的新图，并设置轴标签和标题：

    fig = plt.figure(tight_layout=True)

    ax = fig.add_subplot(projection="3d")

    ax.tick_params(axis="both", which="major", labelsize=9)

    ax.set(xlabel="x", ylabel="y", zlabel="z")

    ax.set_title("Objective function")

  5. 现在，我们可以将目标函数作为表面绘制在刚才创建的坐标轴上：

    ax.plot_surface(x, y, z, cmap="gray",

        vmax=8.0, alpha=0.5)

  6. 我们选择一个初始点，作为最小化例程开始迭代的起点，并将其绘制在表面上：

    x0 = np.array([-0.5, 1.0])

    ax.plot([x0[0]], [x0[1]], func(x0), "k*")

目标函数表面的绘图，以及初始点，可以在以下图示中看到。这里，我们可以看到最小值似乎出现在 x 轴大约 0.5，y 轴大约-0.5 的位置：

图 9.3 – 一个包含起始点的非线性目标函数

图 9.3 – 一个包含起始点的非线性目标函数

  1. 现在，我们使用optimize包中的minimize例程来寻找最小值，并打印出它生成的result对象：

    result = optimize.minimize(

        func, x0, tol=1e-6, method= "Nelder-Mead")

    print(result)

  2. 最后，我们将minimize例程找到的最小值绘制在目标函数的表面上：

    ax.plot([result.x[0]], [result.x[1]], [result.fun], "kx")

更新后的目标函数图，包含minimize例程找到的最小点，可以在以下图示中看到：

图 9.4 – 一个包含起始点和最小点的目标函数

图 9.4 – 一个包含起始点和最小点的目标函数

这表明该方法确实在从初始点（左上角）开始的区域内找到了最小点（右下角）。

它是如何工作的...

Nelder-Mead 单纯形法——与线性优化问题的单纯形法不同——是一种简单的算法，用于找到非线性函数的最小值，即使该目标函数没有已知的导数也能有效工作。（对于这个示例中的函数并非如此；使用基于梯度的方法唯一的收益是收敛速度的提升。）该方法通过比较单纯形顶点处目标函数的值来工作，单纯形在二维空间中是一个三角形。具有最大函数值的顶点会被反射至对边，并执行适当的扩展或收缩，从而使单纯形向下移动。

来自 SciPy optimize模块的minimize例程是许多非线性函数最小化算法的入口点。在这个例子中，我们使用了 Nelder-Mead 单纯形算法，但也有许多其他算法可供选择。许多这些算法需要了解函数的梯度，梯度可能由算法自动计算。通过提供适当的名称给method关键字参数，就可以使用该算法。

minimize例程返回的result对象包含了关于找到的解（如果没有找到解则为错误）的大量信息。特别地，计算出的最小值对应的目标参数存储在结果的x属性中，而函数值存储在fun属性中。

minimize例程需要目标函数和x0的起始值。在这个示例中，我们还提供了一个容差值，要求最小值使用tol关键字参数进行计算。更改该值将会修改计算解的精度。

还有更多内容……

Nelder-Mead 算法是一个无梯度最小化算法的例子，因为它不需要任何关于目标函数梯度（导数）的信息。有几种类似的算法，它们通常涉及在多个点上评估目标函数，然后利用这些信息向最小值移动。一般来说，无梯度方法的收敛速度通常比梯度下降模型慢。然而，它们可以用于几乎任何目标函数，即使在无法准确计算梯度或通过近似计算梯度的情况下也能使用。

优化单变量函数通常比多维情况更容易，并且在 SciPy 的optimize库中有一个专门的函数。minimize_scalar例程用于执行单变量函数的最小化，在这种情况下应替代minimize使用。

在优化中使用梯度下降方法

在之前的示例中，我们使用了 Nelder-Mead 单纯形算法来最小化一个包含两个变量的非线性函数。这是一种相当稳健的方法，即使对目标函数了解甚少，它也能有效工作。然而，在许多情况下，我们对目标函数有更多的了解，这使我们能够设计出更快速、更高效的算法来最小化该函数。我们可以通过利用诸如函数梯度等属性来实现这一点。

一个多变量函数的梯度描述了该函数在各个分量方向上的变化速率。这是该函数关于每个变量的偏导数向量。通过这个梯度向量，我们可以推断出函数增长最快的方向，反之，也可以推断出函数下降最快的方向。从任意给定位置出发，这为我们提供了梯度下降方法的基础，用于最小化一个函数。算法非常简单：给定一个起始位置，，我们计算在的梯度以及梯度下降最快的对应方向，然后朝着该方向迈出一个小步伐。经过几次迭代后，这将使我们从起始位置移动到函数的最小值。

在这个食谱中，我们将学习如何实现一个基于最速下降法的算法，在一个有界区域内最小化目标函数。

准备开始

对于这个食谱，我们需要导入 NumPy 包作为np，Matplotlib 的pyplot模块作为plt，以及从mpl_toolkits.mplot3d导入Axes3D对象：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

让我们实现一个简单的梯度下降算法，并用它来解决前面食谱中描述的最小化问题，看看它是如何工作的。

如何做…

在接下来的步骤中，我们将实现一个简单的梯度下降法，用来最小化一个已知梯度函数的目标函数（实际上我们将使用一个生成器函数，以便在方法运行时看到它的工作过程）：

  1. 我们将从定义一个descend例程开始，它将执行我们的算法。函数声明如下：

    def descend(func,x0,grad,bounds,tol=1e-8,max_iter=100):

  2. 接下来，我们需要实现这个例程。我们从定义在方法运行过程中存储迭代值的变量开始：

        xn = x0

        previous = np.inf

        grad_xn = grad(x0)

  3. 然后我们开始我们的循环，它将执行迭代。在继续之前，我们立即检查是否在取得有意义的进展：

        for i in range(max_iter):

            if np.linalg.norm(xn - previous) < tol:

                break

  4. 方向是梯度向量的负值。我们计算一次并将其存储在direction变量中：

            direction = -grad_xn

  5. 现在，我们更新前一个和当前的值，分别是xnm1和xn，为下一次迭代做好准备。这就结束了descend例程的代码：

            previous = xn

            xn = xn + 0.2*direction

  6. 现在，我们可以计算当前值的梯度并生成所有适当的值：

            grad_xn = grad(xn)

            yield i, xn, func(xn), grad_xn

这就结束了descend例程的定义。

  1. 我们现在可以定义一个示例目标函数进行最小化：

    def func(x):

        return ((x[0] - 0.5)**2 + (

            x[1] + 0.5)**2)*np.cos(0.5*x[0]*x[1])

  2. 接下来，我们创建一个网格，用于评估并绘制目标函数：

    x_r = np.linspace(-1, 1)

    y_r = np.linspace(-2, 2)

    x, y = np.meshgrid(x_r, y_r)

  3. 一旦创建了网格，我们可以评估我们的函数并将结果存储在z变量中：

    z = func([x, y])

  4. 接下来，我们创建目标函数的三维曲面图：

    surf_fig = plt.figure(tight_layout=True)

    surf_ax = surf_fig.add_subplot(projection="3d")

    surf_ax.tick_params(axis="both", which="major",

        labelsize=9)

    surf_ax.set(xlabel="x", ylabel="y", zlabel="z")

    surf_ax.set_title("Objective function")

    surf_ax.plot_surface(x, y, z, cmap="gray", 

        vmax=8.0, alpha=0.5)

  5. 在开始最小化过程之前，我们需要定义一个初始点x0。我们将这个点绘制在我们之前创建的目标函数图上：

    x0 = np.array([-0.8, 1.3])

    surf_ax.plot([x0[0]], [x0[1]], func(x0), "k*")

目标函数的表面图以及初始值可以在下图中看到：

图 9.5 – 目标函数的表面图和初始位置

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

图 9.5 – 目标函数的表面图和初始位置

  1. 我们的descend例程需要一个评估目标函数梯度的函数，因此我们将定义一个：

    def grad(x):

        c1 = x[0]**2 - x[0] + x[1]**2 + x[1] + 0.5

        cos_t = np.cos(0.5*x[0]*x[1])

        sin_t = np.sin(0.5*x[0]*x[1])

        return np.array([

            (2*x[0]-1)*cos_t - 0.5*x[1]*c1*sin_t,

            (2*x[1]+1)*cos_t - 0.5*x[0]*c1*sin_t

        ])

  2. 我们将把迭代次数绘制在等高线图上，因此我们按以下方式进行设置：

    cont_fig, cont_ax = plt.subplots()

    cont_ax.set(xlabel="x", ylabel="y")

    cont_ax.set_title("Contour plot with iterates")

    cont_ax.contour(x, y, z, levels=25, cmap="gray",

        vmax=8.0, opacity=0.6)

  3. 现在，我们创建一个变量，用来保存在和方向上的边界，这些边界是一个元组的元组。这些边界与第 10 步中的linspace调用的边界相同：

    bounds = ((-1, 1), (-2, 2))

  4. 现在，我们可以使用for循环驱动descend生成器来生成每次迭代，并将步骤添加到等高线图中：

    xnm1 = x0

    for i, xn, fxn, grad_xn in descend(func, x0, grad, bounds):

        cont_ax.plot([xnm1[0], xn[0]], [xnm1[1], xn[1]],           	        "k*--")

        xnm1, grad_xnm1 = xn, grad_xn

  5. 循环完成后，我们将最终的值打印到终端：

    print(f"iterations={i}")

    print(f"min val at {xn}")

    print(f"min func value = {fxn}")

之前print语句的输出如下：

iterations=37
min val at [ 0.49999999 -0.49999999]
min func value = 2.1287163880894953e-16

在这里，我们可以看到我们的例程使用了 37 次迭代，找到了大约在(0.5, -0.5)的最小值，这是正确的。

带有迭代次数的等高线图可以在下图中看到：

图 9.6 – 目标函数的等高线图，梯度下降迭代到最小值

图 9.6 – 目标函数的等高线图，显示梯度下降迭代到最小值

在这里，我们可以看到每次迭代的方向——由虚线表示——是目标函数下降最迅速的方向。最终的迭代位于目标函数的碗的中心，即最小值所在的位置。

它是如何工作的...

这个食谱的核心是descend例程。该例程中定义的过程是梯度下降法的一个非常简单的实现。在给定点计算梯度由grad参数处理，然后通过direction = -grad来推导迭代的行进方向。我们将这个方向乘以一个固定的比例因子（有时称为0.2*direction），然后加到当前位置。

这个食谱中的解法在收敛时花费了 37 次迭代，相较于最小化非线性函数食谱中的 Nelder-Mead 单纯形算法（该算法花费了 58 次迭代），这是一个轻微的改进。（这不是一个完美的比较，因为我们更改了这个食谱的起始位置。）该性能在很大程度上依赖于我们选择的步长。在这种情况下，我们将最大步长固定为方向向量大小的 0.2 倍。这使得算法保持简单，但效率并不是特别高。

在这个例子中，我们选择将算法实现为生成器函数，这样可以在每一步迭代中看到输出，并将其绘制在等高线图上。当我们逐步执行迭代时，我们可以看到每一步的结果。实际上，我们可能不想这样做，而是希望在迭代完成后直接返回计算得到的最小值。为了实现这一点，我们只需移除yield语句，并在函数的最后（即主函数的缩进位置，不在循环内部）将其替换为return xn。如果你想防止不收敛的情况，可以使用for循环的else特性来捕捉循环因到达迭代器末尾而没有触发break关键字的情况。这个else块可以引发异常，表示算法未能收敛到一个解。在这个例子中，我们用来结束迭代的条件并不能保证方法已经达到最小值，但通常情况下它是有效的。

还有更多…

实际上，你通常不会自己实现梯度下降算法，而是会使用像 SciPy optimize 模块这样的库中的通用函数。我们可以使用在前一个例子中使用的相同minimize函数，通过多种不同的算法执行最小化，包括几种梯度下降算法。这些实现通常会比像这种自定义实现具有更高的性能和更强的鲁棒性。

我们在这个例子中使用的梯度下降方法是一个非常简单的实现，可以通过允许算法在每一步选择步长来大大改进。（允许自己选择步长的方法有时被称为自适应方法。）这种改进的难点在于如何选择在此方向上采取的步长。为此，我们需要考虑一个单变量的函数，该函数由以下方程给出：

在这里，表示当前点，表示当前方向，是一个参数。为了简便起见，我们可以使用一个名为minimize_scalar的最小化例程来处理标量值函数，来自 SciPy optimize 模块。不幸的是，这并不是像传入这个辅助函数并找到最小值那样简单。我们需要给![](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/b7e32bdd82774789972f4e1cde1ef50a~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1770865377&x-signature=DEGzIuAHdOJIakOUTx5FIiyWzt8%3D)的可能值设置边界，以便计算得到的最小化点，![](https://p9-xtjj-sign.byteimg.com/tos-cn-i-73owjymdk6/dafdd1e607344c8fa1c8f2a261cc3c8c~tplv-73owjymdk6-jj-mark-v1:0:0:0:0:5o6Y6YeR5oqA5pyv56S-5Yy6IEAg5biD5a6i6aOe6b6Z:q75.awebp?rk3s=f64ab15b&x-expires=1770865377&x-signature=Ea6bqK%2BthvfR649xiH%2F%2FQ7FIeg4%3D)，位于我们感兴趣的区域内。

要理解我们如何约束 的值，我们必须先从几何角度看构建过程。我们引入的辅助函数评估了沿给定方向的单条直线上的目标函数。我们可以将其想象为对表面进行单一的截面切割，切割面通过当前的 点，并沿 方向延伸。算法的下一步是寻找步长 ，它最小化沿该直线的目标函数值——这是一个标量函数，比较容易最小化。然后，边界值应该是 的值范围，在此范围内，直线将位于由 和 边界值定义的矩形内。我们确定这条直线与 和 边界线交叉的四个值，其中两个是负数，两个是正数。（这是因为当前点必须位于矩形内。）我们取两个正值中的最小值和两个负值中的最大值，并将这些边界传递给标量最小化程序。可以使用以下代码实现此过程：

alphas = np.array([
    (bounds[0][0] - xn[0]) / direction[0],
    # x lower
   (bounds[1][0] - xn[1]) / direction[1],
    # y lower
    (bounds[0][1] - xn[0]) / direction[0],
    # x upper
    (bounds[1][1] - xn[1]) / direction[1] 
    # y upper
])
alpha_max = alphas[alphas >= 0].min()
alpha_min = alphas[alphas < 0].max()
result = minimize_scalar(lambda t: 
    func(xn + t*direction),
    method="bounded",
    bounds=(alpha_min, alpha_max))
amount = result.x

一旦步长被选定，剩下的唯一步骤是更新当前的 xn 值，如下所示：

xn = xn + amount * direction

使用这种自适应步长增加了程序的复杂性，但性能得到了显著提高。使用这种改进的程序，方法仅在三次迭代内收敛，这比本食谱中朴素代码使用的迭代次数（37 次）或上一食谱中 Nelder-Mead 单纯形算法的迭代次数（58 次）要少得多。迭代次数的减少正是我们通过提供梯度函数这一信息而预期的结果。

我们创建了一个函数，返回给定点处的函数梯度。我们在开始之前手动计算了这个梯度，但这并不总是容易甚至可能的。相反，通常会将这里使用的解析梯度替换为一个数值计算梯度，该梯度通过有限差分或类似算法估算得到。这会对性能和精度产生影响，就像所有的近似方法一样，但鉴于梯度下降方法在收敛速度上的改进，这些问题通常是微不足道的。

梯度下降类算法在机器学习应用中尤其流行。大多数流行的 Python 机器学习库——包括 PyTorch、TensorFlow 和 Theano——都提供了用于自动计算数据数组梯度的工具。这使得梯度下降方法可以在后台使用，从而提高性能。

梯度下降法的一个流行变体是随机梯度下降法，其中梯度是通过随机抽样估算的，而不是使用整个数据集。这可以显著减轻该方法的计算负担——虽然代价是收敛速度较慢——特别是在高维问题中，如机器学习应用中常见的问题。随机梯度下降法通常与反向传播结合，成为训练人工神经网络在机器学习中的基础。

基本的随机梯度下降算法有几种扩展。例如，动量算法将前一个增量融入到下一个增量的计算中。另一个例子是自适应梯度算法，它结合了每个参数的学习率，以提高涉及大量稀疏参数问题的收敛速度。

使用最小二乘法拟合曲线到数据

最小二乘法是一种强大的技术，用于从相对较小的潜在函数家族中找到最能描述特定数据集的函数。该技术在统计学中尤为常见。例如，最小二乘法用于线性回归问题——在这里，潜在函数的家族是所有线性函数的集合。通常，我们尝试拟合的函数家族具有相对较少的可调参数，用以解决问题。

最小二乘法的理念相对简单。对于每个数据点，我们计算残差的平方——即点的值与给定函数的期望值之间的差异——并尝试将这些残差平方的和尽可能小（因此称为最小二乘法）。

在本教程中，我们将学习如何使用最小二乘法拟合一条曲线到样本数据集。

准备工作

对于本教程，我们将像往常一样导入 NumPy 包，命名为 np，并导入 Matplotlib pyplot 模块，命名为 plt：

import numpy as np
import matplotlib.pyplot as plt

我们还需要从 NumPy random 模块导入默认的随机数生成器实例，如下所示：

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

最后，我们需要 SciPy optimize 模块中的 curve_fit 函数：

from scipy.optimize import curve_fit

让我们看看如何使用这个函数将一条非线性曲线拟合到一些数据上。

如何实现...

以下步骤展示了如何使用 curve_fit 函数将一条曲线拟合到一组数据：

  1. 第一步是创建样本数据：

    SIZE = 100

    x_data = rng.uniform(-3.0, 3.0, size=SIZE)

    noise = rng.normal(0.0, 0.8, size=SIZE)

    y_data = 2.0*x_data**2 - 4*x_data + noise

  2. 接下来，我们生成数据的散点图，以查看是否能识别数据中的潜在趋势：

    fig, ax = plt.subplots()

    ax.scatter(x_data, y_data)

    ax.set(xlabel="x", ylabel="y",

        title="Scatter plot of sample data")

我们所生成的散点图如下所示。这里，我们可以看到数据显然不遵循线性趋势（直线）。既然我们知道趋势是多项式，那么我们接下来的猜测是二次趋势。这就是我们在这里使用的：

图 9.7 – 样本数据的散点图 – 我们可以看到数据并没有遵循线性趋势

图 9.7 – 样本数据的散点图——我们可以看到它没有遵循线性趋势

  1. 接下来，我们创建一个表示我们希望拟合的模型的函数：

    def func(x, a, b, c):

        return a*x**2 + b*x + c

  2. 现在，我们可以使用 curve_fit 例程将模型函数拟合到样本数据中：

    coeffs, _ = curve_fit(func, x_data, y_data)

    print(coeffs)

    # [ 1.99611157 -3.97522213 0.04546998]

  3. 最后，我们将最佳拟合曲线绘制到散点图上，以评估拟合曲线对数据的描述效果：

    x = np.linspace(-3.0, 3.0, SIZE)

    y = func(x, coeffs[0], coeffs[1], coeffs[2])

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

更新后的散点图可以在下图中看到：

图 9.8 – 一个使用叠加最小二乘法找到的最佳拟合曲线的散点图

图 9.8 – 一个使用叠加最小二乘法找到的最佳拟合曲线的散点图

在这里，我们可以看到我们找到的曲线拟合数据的效果相当好。系数并不完全等于真实模型——这是由于添加的噪声效应。

它是如何工作的...

curve_fit 例程执行最小二乘法拟合，将模型的曲线拟合到样本数据中。在实践中，这相当于最小化以下目标函数：

这里，配对的 是样本数据中的点。在这种情况下，我们在一个三维参数空间中进行优化，每个维度对应一个参数。例程返回估计的系数——在参数空间中最小化目标函数的点——以及一个包含拟合协方差矩阵估计值的第二个变量。在这个例子中，我们忽略了这一部分。

从 curve_fit 例程返回的估计协方差矩阵可以用来为估计的参数提供置信区间。这是通过取对角元素的平方根并除以样本大小（本例中为 100）来实现的。这会给出估计的标准误差，当乘以对应于置信度的适当值时，就能得到置信区间的大小。（我们在第六章 数据与统计学应用 中讨论了置信区间。）

你可能已经注意到，curve_fit 例程估计的参数与我们在步骤 1 中用来定义样本数据的参数接近，但并不完全相等。之所以没有完全相等，是因为我们向数据中添加了正态分布噪声。在这个例子中，我们知道数据的潜在结构是二次型的——也就是二次多项式——而不是其他更复杂的函数。实际上，我们不太可能如此清楚数据的潜在结构，这也是我们向样本中添加噪声的原因。

还有更多内容...

SciPy optimize模块中有另一个执行最小二乘拟合的例程，称为least_squares。这个例程的签名稍微不那么直观，但确实返回一个包含更多优化过程信息的对象。然而，这个例程的设置方式可能更类似于我们在*如何运作...*部分中构造底层数学问题的方式。要使用这个例程，我们定义目标函数如下：

def func(params, x, y):
    return y -(
        params[0]*x**2 + params[1]*x + params[2])

我们将这个函数和在参数空间中的初始估计x0一起传递，例如(1, 0, 0)。目标函数func的附加参数可以通过args关键字参数传递——例如，我们可以使用args=(x_data, y_data)。这些参数会被传递到目标函数的x和y参数中。总结来说，我们可以通过以下方式调用least_squares来估计参数：

results = least_squares(func, [1, 0, 0], args=(x_data, y_data))

从least_squares例程返回的results对象实际上与本章描述的其他优化例程返回的对象相同。它包含诸如使用的迭代次数、过程是否成功、详细的错误信息、参数值以及目标函数在最小值处的值等细节。

分析简单的二人游戏

博弈论是研究决策和策略分析的数学分支。它在经济学、生物学和行为科学中有广泛应用。许多看似复杂的情境可以简化为一个相对简单的数学博弈，经过系统分析后能够找到最优解决方案。

博弈论中的经典问题是囚徒困境，其原始形式如下：两个同谋被抓获，必须决定是否保持沉默或指证对方。如果两人都保持沉默，他们都服刑 1 年；如果一个人作证而另一个人不作证，作证的人被释放，另一个人服刑 3 年；如果两人都相互作证，他们都服刑 2 年。那么每个同谋应该怎么做呢？事实证明，考虑到对对方的合理不信任，每个同谋做出的最佳选择是作证。采取这种策略，他们要么不服刑，要么最多服刑 2 年。

由于这本书是关于 Python 的，我们将用这个经典问题的变种来说明这个问题的普遍性。考虑以下问题：你和你的同事必须为客户编写一些代码。你认为你可以用 Python 写得更快，但你的同事认为他们可以用 C 语言写得更快。问题是，你们应该为项目选择哪种语言？

你认为自己可以写出比 C 快四倍的 Python 代码，所以你以速度 1 写 C，速度 4 写 Python。你的同事说他们写 C 稍微比写 Python 快一些，所以他们写 C 以速度 3，写 Python 以速度 2。如果你们都同意一个语言，那么你按照自己预测的速度写代码；但是如果你们意见不合，那么较快的程序员的生产力会减少 1。我们可以总结如下：

同事/你CPython
C3/13/2
Python2/12/4

图 9.9 – 预测各种配置下的工作速度表

在这个教程中，我们将学习如何在 Python 中构建一个对象来表示这个简单的两人游戏，并对游戏的结果进行一些基础分析。

准备工作

对于这个教程，我们需要导入 NumPy 包为 np，以及导入 Nashpy 包为 nash：

import numpy as np
import nashpy as nash

让我们看看如何使用 nashpy 包来分析一个简单的两人游戏。

如何实现...

以下步骤将向你展示如何使用 Nashpy 创建并进行一些简单的两人博弈分析：

  1. 首先，我们需要创建存储每个玩家回报信息的矩阵（在这个例子中是你和你的同事）：

    you = np.array([[1, 3], [1, 4]])

    colleague = np.array([[3, 2], [2, 2]])

  2. 接下来，我们创建一个 Game 对象，该对象包含了由这些回报矩阵表示的游戏：

    dilemma = nash.Game(you, colleague)

  3. 我们通过索引符号计算给定选择的效用：

    print(dilemma[[1, 0], [1, 0]])      # [1 3]

    print(dilemma[[1, 0], [0, 1]])      # [3 2]

    print(dilemma[[0, 1], [1, 0]])      # [1 2]

    print(dilemma[[0, 1], [0, 1]])      # [4 2]

  4. 我们还可以根据选择某个特定决策的概率来计算期望效用：

    print(dilemma[[0.1, 0.9], [0.5, 0.5]]) # [2.45 2.05]

这些期望效用表示我们期望（平均来说）在多次重复游戏时看到的结果，按照指定的概率。

它是如何工作的...

在这个教程中，我们构建了一个 Python 对象，表示一个非常简单的两人战略游戏。这里的思路是有两个玩家需要做出决策，每种玩家选择的组合会给出一个特定的回报值。我们在这里的目标是找到每个玩家可以做出的最佳选择。假设两个玩家同时做出一个决定，在这个过程中，任何一个玩家都不知道对方的选择。每个玩家都有一个策略来决定他们的选择。

在步骤 1中，我们为每个玩家创建两个矩阵，每个矩阵对应每种选择组合的回报值。这两个矩阵被 Nashpy 的 Game 类包装，提供了一个便捷且直观（从博弈论角度看）的界面来处理游戏。我们可以通过传入选择并使用索引符号，快速计算给定选择组合的效用。

我们还可以基于一种策略计算预期效用，其中根据某种概率分布随机选择。语法与前面描述的确定性情况相同，只是我们为每个选择提供一个概率向量。我们根据您 90%的时间选择 Python，而您的同事选择 Python 的概率为 50%来计算预期速度，分别为 2.45 和 2.05。

还有更多内容...

在 Python 中，计算游戏理论的另一种方法是使用 Gambit 项目。Gambit 项目是一个工具集，用于计算博弈论，具有 Python 接口（www.gambit-project.org/）。这是一个成熟的项目，构建在 C 库的基础上，并提供比 Nashpy 更高的性能。

计算纳什均衡

纳什均衡是一个两人策略游戏——类似于我们在分析简单的两人游戏配方中看到的——它代表了一个稳定状态，在这种状态下每个玩家都看到了最佳可能的结果。然而，这并不意味着与纳什均衡相关联的结果在总体上是最好的。纳什均衡比这更微妙。一个非正式的纳什均衡定义如下：在这个动作配置中，假设所有其他玩家遵循这个配置，没有一个个体玩家可以改善他们的结果。

我们将通过经典的剪刀石头布游戏探讨纳什均衡的概念。规则如下。每个玩家可以选择以下选项之一：剪刀，石头或布。石头打败剪刀，但输给布；纸打败石头，但输给剪刀；剪刀打败纸，但输给石头。如果两个玩家选择相同的选项，则为平局。数值上，我们用+1 表示赢，-1 表示输，0 表示平局。从中，我们可以构建一个两人游戏并计算该游戏的纳什均衡。

在这个配方中，我们将计算经典剪刀石头布游戏的纳什均衡。

准备工作

对于这个配方，我们需要导入 NumPy 包作为np，并导入 Nashpy 包作为nash：

import numpy as np
import nashpy as nash

让我们看看如何使用nashpy包来计算两人策略游戏的纳什均衡。

怎么做...

下面的步骤展示了如何计算一个简单两人游戏的纳什均衡：

  1. 首先，我们需要为每个玩家创建一个收益矩阵。我们从第一个玩家开始：

    rps_p1 = np.array([

        [ 0, -1, 1], # rock payoff

        [ 1, 0, -1], # paper payoff

        [-1, 1, 0] # scissors payoff

    ])

  2. 第二个玩家的收益矩阵是rps_p1的转置：

    rps_p2 = rps_p1.transpose()

  3. 接下来，我们创建代表游戏的Game对象：

    rock_paper_scissors = nash.Game(rps_p1, rps_p2)

  4. 我们使用支持枚举算法计算游戏的纳什均衡：

    equilibria = rock_paper_scissors.support_enumeration()

  5. 我们遍历均衡，并打印每个玩家的策略：

    for p1, p2 in equilibria:

        print("Player 1", p1)

        print("Player 2", p2)

这些打印语句的输出如下：

Player 1 [0.33333333 0.33333333 0.33333333]
Player 2 [0.33333333 0.33333333 0.33333333]

工作原理...

纳什均衡在博弈论中极其重要，因为它们使我们能够分析战略博弈的结果，并识别有利的位置。它们最早由约翰·纳什在 1950 年描述，并在现代博弈论中发挥了关键作用。一个二人博弈可能有多个纳什均衡，但任何有限的二人博弈至少必须有一个纳什均衡。问题在于找到一个给定博弈的所有可能的纳什均衡。

在这个案例中，我们使用了支持集枚举法，它有效地枚举了所有可能的策略，并筛选出纳什均衡策略。在这个案例中，支持集枚举算法仅找到了一个纳什均衡，它是一种混合策略。这意味着唯一没有改进的策略是随机选择其中一个选项，每个选项的概率为 1/3。对于玩过石头剪子布的人来说，这并不令人惊讶，因为无论我们做出什么选择，我们的对手都有 1/3 的概率（随机）选择一个能够击败我们选择的动作。同样，我们也有 1/3 的机会平局或获胜，因此我们在所有这些可能性下的期望值如下：

在不知道我们的对手会选择哪个选项的情况下，没有办法改进这个预期结果。

还有更多...

Nashpy 包还提供了其他计算纳什均衡的算法。具体来说，vertex_enumeration方法在Game对象上使用顶点枚举算法，而lemke_howson_enumeration方法使用Lemke-Howson算法。这些替代算法具有不同的特性，可能在某些问题中更高效。

另见

Nashpy 包的文档包含关于所涉及算法和博弈论的更详细信息。这其中包括许多关于博弈论的参考文献。该文档可以在nashpy.readthedocs.io/en/latest/找到。

进一步阅读

和往常一样，Numerical Recipes书是数值算法的一个很好的来源。第十章，函数的最小化与最大化，讨论了函数的最大化和最小化：

  • Press, W.H., Teukolsky, S.A., Vetterling, W.T., 和 Flannery, B.P., 2017. 数值计算的艺术：数值算法。第三版。剑桥：剑桥大学出版社。

关于优化的更多具体信息可以在以下书籍中找到：

  • Boyd, S.P. 和 Vandenberghe, L., 2018. 凸优化。剑桥：剑桥大学出版社。

  • Griva, I., Nash, S., 和 Sofer, A., 2009. 线性和非线性优化。第二版。费城：工业与应用数学学会。

最后，以下书籍是博弈论的一个很好的入门书籍：

  • Osborne, M.J., 2017. 博弈论导论。牛津：牛津大学出版社。

第十章：提高生产力

本章将讨论一些不属于本书前几章分类的主题。这些主题大多数与简化计算过程和优化代码执行有关。还有一些涉及处理特定类型的数据或文件格式。

本章的目的是为您提供一些工具，虽然它们本质上并非严格的数学工具，但在数学问题中经常出现。这些工具包括分布式计算和优化等主题——它们帮助您更快速地解决问题，验证数据和计算结果，从常见的科学计算文件格式加载和存储数据，并且融入其他有助于提高代码生产力的主题。

在前两种配方中，我们将介绍帮助跟踪计算中的单位和不确定性的包。这些对处理具有直接物理应用的数据的计算非常重要。在接下来的配方中，我们将讨论如何从 网络通用数据格式（NetCDF）文件中加载和存储数据。NetCDF 是一种通常用于存储天气和气候数据的文件格式。在第四个配方中，我们将讨论如何处理地理数据，如可能与天气或气候数据相关的数据。接下来，我们将讨论如何从终端运行 Jupyter 笔记本，而无需启动交互式会话。然后，我们将转向数据验证，并在接下来的配方中关注使用 Cython 和 Dask 等工具的性能。最后，我们将简要概述一些编写数据科学可重现代码的技术。

本章将涵盖以下配方：

  • 使用 Pint 跟踪单位

  • 计算中的不确定性

  • 从 NetCDF 文件加载和存储数据

  • 处理地理数据

  • 将 Jupyter 笔记本作为脚本执行

  • 数据验证

  • 使用 Cython 加速代码

  • 使用 Dask 分布式计算

  • 为数据科学编写可重现的代码

让我们开始吧！

技术要求

由于本章包含的内容类型，本章需要许多不同的包。我们需要的包清单如下：

  • Pint

  • uncertainties

  • netCDF4

  • xarray

  • Pandas

  • Scikit-learn

  • GeoPandas

  • Geoplot

  • Jupyter

  • Papermill

  • Cerberus

  • Cython

  • Dask

所有这些包都可以通过您喜欢的包管理工具（如 pip）进行安装：

python3.10 -m pip install pint uncertainties netCDF4 xarray pandas scikit-learn geopandas geoplot jupyter papermill cerberus cython

要安装 Dask 包，我们需要安装与该包相关的各种附加功能。可以通过在终端中使用以下 pip 命令来完成：

python3.10 -m pip install dask[complete]

除了这些 Python 包，我们还需要安装一些支持软件。对于处理地理数据配方，GeoPandas 和 Geoplot 库有许多较低级别的依赖项，可能需要单独安装。详细的安装说明请参考 GeoPandas 包文档，网址：geopandas.org/install.html。

对于使用 Cython 加速代码配方，我们需要安装 C 编译器。有关如何获得GNU C 编译器（GCC）的说明，请参见 Cython 文档：cython.readthedocs.io/en/latest/src/quickstart/install.html。

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

使用 Pint 追踪单位

正确追踪计算中的单位可能非常困难，特别是当不同单位可以互换使用时。例如，非常容易忘记转换不同的单位——英尺/英寸转化为米，或者公制前缀——例如将 1 千米转换为 1,000 米。

在这个配方中，我们将学习如何使用 Pint 包在计算中追踪单位。

准备工作

对于这个配方，我们需要安装 Pint 包，可以通过以下方式导入：

import pint

如何操作...

以下步骤展示了如何使用 Pint 包在计算中追踪单位：

  1. 首先，我们需要创建一个UnitRegistry对象：

    ureg = pint.UnitRegistry(system="mks")

  2. 要创建一个带有单位的量，我们将数字乘以注册对象的相应属性：

    distance = 5280 * ureg.feet

  3. 我们可以使用其中一个可用的转换方法来更改量的单位：

    print(distance.to("miles"))

    print(distance.to_base_units())

    print(distance.to_base_units().to_compact())

这些print语句的输出如下：

0.9999999999999999 mile
1609.3439999999998 meter
1.6093439999999999 kilometer

  1. 我们封装一个例程，使其期望接受秒为单位的参数，并输出以米为单位的结果：

    @ureg.wraps(ureg.meter, ureg.second)

    def calc_depth(dropping_time):

        # s = u*t + 0.5*a*t*t

        # u = 0, a = 9.81

        return 0.5*9.81*dropping_time*dropping_time

  2. 现在，当我们用minute单位调用calc_depth例程时，它会自动转换为秒进行计算：

    depth = calc_depth(0.05 * ureg.minute)

    print("Depth", depth)

    # Depth 44.144999999999996 meter

它是如何工作的...

Pint 包提供了一个包装类，用于数值类型，给类型添加了单位元数据。这个包装类型实现了所有标准的算术操作，并在整个计算过程中追踪单位。例如，当我们将长度单位除以时间单位时，我们会得到一个速度单位。这意味着，您可以使用 Pint 来确保在复杂计算后单位是正确的。

UnitRegistry对象追踪会话中所有存在的单位，并处理不同单位类型之间的转换等操作。它还维护着一个度量参考系统，在本配方中，是以米、千克和秒为基本单位的国际标准系统，称为mks。

wraps 功能允许我们声明一个例程的输入和输出单位，这使得 Pint 能够对输入函数进行自动单位转换——在这个示例中，我们将分钟转换为秒。试图使用没有关联单位的数量，或者使用不兼容单位调用包装函数，会抛出异常。这使得在运行时可以验证参数，并自动转换成例程所需的正确单位。

还有更多内容...

Pint 包提供了一个包含大量预编程单位的测量单位列表，覆盖了大多数全球使用的系统。单位可以在运行时定义，也可以从文件中加载。这意味着你可以定义特定于应用程序的自定义单位或单位系统。

单位也可以在不同的上下文中使用，这使得在通常不相关的不同单位类型之间轻松转换成为可能。这可以在计算过程中需要频繁转换单位时节省大量时间。

在计算中考虑不确定性

大多数测量设备并非 100% 准确，而是具有一定的准确度，通常在 0 到 10% 之间。例如，一个温度计的准确度可能为 1%，而一把数显卡尺的准确度可能为 0.1%。在这两种情况下，真实值不太可能与报告值完全一致，尽管它们会非常接近。追踪一个值的不确定性是非常困难的，特别是当多个不确定性以不同的方式结合时。与其手动跟踪这些，不如使用一个一致的库来代替。这正是 uncertainties 包的功能。

在这个示例中，我们将学习如何量化变量的不确定性，并观察这些不确定性是如何在计算中传播的。

准备工作

对于这个示例，我们将需要 uncertainties 包，其中我们将导入 ufloat 类和 umath 模块：

from uncertainties import ufloat, umath

如何操作...

以下步骤展示了如何在计算中量化数值的不确定性：

  1. 首先，我们创建一个不确定浮动值 3.0 加减 0.4：

    seconds = ufloat(3.0, 0.4)

    print(seconds)      # 3.0+/-0.4

  2. 接下来，我们进行涉及该不确定值的计算，以获得一个新的不确定值：

    depth = 0.5*9.81*seconds*seconds

    print(depth)      # 44+/-12

  3. 接下来，我们创建一个新的不确定浮动值，并应用 umath 模块中的 sqrt 例程，进行前一个计算的反向操作：

    other_depth = ufloat(44, 12)

    time = umath.sqrt(2.0*other_depth/9.81)

    print("Estimated time", time)

    # Estimated time 3.0+/-0.4

如我们所见，第一次计算（步骤 2）的结果是一个不确定浮动值，其数值为 44，并且存在 系统误差。这意味着真实值可能在 32 到 56 之间。根据我们现有的测量数据，我们无法得出更精确的结果。

它是如何工作的...

ufloat类封装了float对象，并在整个计算过程中追踪不确定性。该库利用线性误差传播理论，使用非线性函数的导数来估算计算中的误差传播。该库还正确处理相关性，因此从自身减去一个值时，结果为零且没有误差。

为了追踪标准数学函数中的不确定性，你需要使用umath模块中提供的版本，而不是使用 Python 标准库或第三方包如 NumPy 中定义的版本。

还有更多...

uncertainties包为 NumPy 提供支持，而之前例子中提到的 Pint 包可以与uncertainties结合使用，以确保单位和误差范围正确地归属于计算的最终值。例如，我们可以计算本例第2 步中计算的单位，如下所示：

import pint
from uncertainties import ufloat
ureg = pint.UnitRegistry(system="mks")
g = 9.81*ureg.meters / ureg.seconds ** 2
seconds = ufloat(3.0, 0.4) * ureg.seconds
depth = 0.5*g*seconds**2
print(depth)

正如预期的那样，最后一行的print语句输出了44+/-12米。

从 NetCDF 文件加载和存储数据

许多科学应用要求我们从大量多维数据开始，并使用稳健的格式。NetCDF 是气象和气候行业开发的一种数据格式示例。不幸的是，数据的复杂性意味着我们不能简单地使用例如 Pandas 包中的工具来加载这些数据进行分析。我们需要netcdf4包来读取并将数据导入到 Python 中，但我们还需要使用xarray。与 Pandas 库不同，xarray可以处理更高维度的数据，同时仍提供类似 Pandas 的接口。

在本例中，我们将学习如何从 NetCDF 文件加载数据并将数据存储到 NetCDF 文件中。

准备工作

对于本例，我们将需要导入 NumPy 包为np，Pandas 包为pd，Matplotlib 的pyplot模块为plt，以及 NumPy 中默认随机数生成器的实例：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from numpy.random import default_rng
rng = default_rng(12345)

我们还需要导入xarray包，并以xr为别名。你还需要安装 Dask 包，如技术要求部分所述，以及netCDF4包：

import xarray as xr

我们不需要直接导入这两个包。

如何操作...

按照以下步骤加载并存储示例数据到 NetCDF 文件：

  1. 首先，我们需要创建一些随机数据。这些数据包含日期范围、位置代码列表以及随机生成的数字：

    dates = pd.date_range("2020-01-01", periods=365, name="date")

    locations = list(range(25))

    steps = rng.normal(0, 1, size=(365,25))

    accumulated = np.add.accumulate(steps)

  2. 接下来，我们创建一个包含数据的 xarray Dataset对象。日期和位置是索引，而steps和accumulated变量是数据：

    data_array = xr.Dataset({

        "steps": (("date", "location"), steps),

        "accumulated": (("date", "location"), accumulated)     },

        {"location": locations, "date": dates}

    )

这里显示了print语句的输出：

<xarray.Dataset>
Dimensions: (date: 365, location: 25)
Coordinates:
* location (location) int64 0 1 2 3 4 5 6 7 8 ... 17 18 19 20 21 22 23 24
* date (date) datetime64[ns] 2020-01-01 2020-01-02 ... 2020-12-30
Data variables:
steps (date, location) float64 geoplot.pointplot(cities, ax=ax, fc="r", marker="2")
ax.axis((-180, 180, -90, 90))-1.424 1.264 ... -0.4547 -0.4873
accumulated (date, location) float64 -1.424 1.264 -0.8707 ... 8.935 -3.525

  1. 接下来，我们计算每个时间索引下所有位置的平均值：

    means = data_array.mean(dim="location")

  2. 现在，我们在新的坐标轴上绘制平均累计值：

    fig, ax = plt.subplots()

    means["accumulated"].to_dataframe().plot(ax=ax)

    ax.set(title="Mean accumulated values", 

        xlabel="date", ylabel="value")

生成的图像如下所示：

![图 10.1 - 累计均值随时间变化的图表]

](github.com/OpenDocCN/f…)

图 10.1 - 随时间变化的累积均值图

  1. 使用to_netcdf方法将此数据集保存为一个新的 NetCDF 文件：

    data_array.to_netcdf("data.nc")

  2. 现在，我们可以使用xarray中的load_dataset例程加载新创建的 NetCDF 文件：

    new_data = xr.load_dataset("data.nc")

    print(new_data)

上述代码的输出如下：

<xarray.Dataset>
Dimensions: (date: 365, location: 25)
Coordinates:
            * location (location) int64 0 1 2 3 4 5 6 7 8 ... 17 18 19 20 21 22 23 24
            * date (date) datetime64[ns] 2020-01-01 2020-01-02 ... 2020-12-30
Data variables:
            steps (date, location) float64 -1.424 1.264 ... -0.4547 -0.4873
            accumulated (date, location) float64 -1.424 1.264 -0.8707 ... 8.935 -3.525

输出显示加载的数组包含了我们在之前步骤中添加的所有数据。关键步骤是5和6，在这两个步骤中，我们存储并加载了这个"data.nc"数据。

它是如何工作的...

xarray包提供了DataArray和DataSet类，它们大致上是 Pandas Series和DataFrame对象的多维等效物。在这个示例中，我们使用数据集，因为每个索引——一个包含日期和位置的元组——都关联着两条数据。这两个对象都暴露了与它们的 Pandas 等效物相似的接口。例如，我们可以使用mean方法计算某个轴上的均值。DataArray和DataSet对象还有一个方便的方法，可以将其转换为 Pandas DataFrame，这个方法叫做to_dataframe。我们在这个食谱中使用它，将means Dataset中的累积列转换为DataFrame以便绘图，虽然这并不是必须的，因为xarray本身就内置了绘图功能。

本食谱的真正重点是to_netcdf方法和load_dataset例程。前者将一个DataSet对象存储为 NetCDF 格式的文件。这需要安装netCDF4包，因为它允许我们访问解码 NetCDF 格式文件所需的相关 C 库。load_dataset例程是一个通用的例程，用于从各种文件格式加载数据到DataSet对象中，包括 NetCDF 格式（同样，这需要安装netCDF4包）。

还有更多...

xarray包除了支持 NetCDF 格式外，还支持多种数据格式，如 OPeNDAP、Pickle、GRIB 等，这些格式也被 Pandas 支持。

处理地理数据

许多应用涉及处理地理数据。例如，在跟踪全球天气时，我们可能想要绘制世界各地传感器测量的温度，并标注它们在地图上的位置。为此，我们可以使用 GeoPandas 包和 Geoplot 包，它们都允许我们操作、分析和可视化地理数据。

在本食谱中，我们将使用 GeoPandas 和 Geoplot 包加载并可视化一些示例地理数据。

准备工作

对于这个食谱，我们将需要导入 GeoPandas 包、Geoplot 包以及 Matplotlib 的pyplot包，命名为plt：

import geopandas
import geoplot
import matplotlib.pyplot as plt

如何实现...

按照以下步骤，使用示例数据在世界地图上绘制首都城市的简单图：

  1. 首先，我们需要加载 GeoPandas 包中的示例数据，这些数据包含了全球的几何信息：

    world = geopandas.read_file(

        geopandas.datasets.get_path("naturalearth_lowres")

    )

  2. 接下来，我们需要加载包含全球每个首都城市名称和位置的数据：

    cities = geopandas.read_file(

        geopandas.datasets.get_path("naturalearth_cities")

    )

  3. 现在，我们可以创建一个新的图形，并使用 polyplot 例程绘制世界地理轮廓：

    fig, ax = plt.subplots()

    geoplot.polyplot(world, ax=ax, alpha=0.7)

  4. 最后，我们使用 pointplot 例程将首都城市的位置叠加到世界地图上。我们还设置了坐标轴的限制，以确保整个世界都能显示出来：

    geoplot.pointplot(cities, ax=ax, fc="k", marker="2")

    ax.axis((-180, 180, -90, 90))

最终得到的图像展示了世界首都城市的位置：

图 10.2 - 世界各国首都城市在地图上的分布

图 10.2 - 世界各国首都城市在地图上的分布

该图展示了世界不同国家的粗略轮廓。每个首都城市的位置通过标记表示。从这个视角看，很难区分中欧地区的各个城市。

工作原理...

GeoPandas 包是 Pandas 的一个扩展，专门用于处理地理数据，而 Geoplot 包是 Matplotlib 的一个扩展，用于绘制地理数据。GeoPandas 包提供了一些示例数据集，我们在本节中使用了这些数据集。naturalearth_lowres 包含描述世界各国边界的几何图形。这些数据的分辨率不高，正如其名称所示，因此地图上可能缺少一些地理特征的细节（一些小岛根本没有显示）。naturalearth_cities 包含世界各国首都城市的名称和位置。我们使用 datasets.get_path 例程来检索这些数据集在包数据目录中的路径。read_file 例程将数据导入到 Python 会话中。

Geoplot 包提供了一些额外的绘图例程，专门用于绘制地理数据。polyplot 例程绘制来自 GeoPandas DataFrame 的多边形数据，这些数据描述了国家的地理边界。pointplot 例程则在一组坐标轴上绘制来自 GeoPandas DataFrame 的离散点，这些点在本例中表示首都城市的位置。

将 Jupyter notebook 作为脚本执行

Jupyter notebook 是一种流行的工具，广泛用于编写科学计算和数据分析应用中的 Python 代码。Jupyter notebook 实际上是一系列存储在 .ipynb 格式文件中的代码块。每个代码块可以是不同类型的，比如代码块或 Markdown 块。这些 notebook 通常通过 Web 应用访问，Web 应用解释这些代码块并在后台的内核中执行代码，结果再返回给 Web 应用。如果你在个人电脑上工作，这种方式非常方便，但如果你想在远程服务器上执行 notebook 中的代码呢？在这种情况下，可能连 Jupyter Notebook 提供的 Web 界面都无法访问。papermill 包允许我们在命令行中对 notebook 进行参数化和执行。

在本节中，我们将学习如何使用 papermill 从命令行执行 Jupyter notebook。

准备工作

对于此方法，我们需要安装 papermill 包，并且在当前目录中有一个示例 Jupyter notebook。我们将使用本章代码库中存储的 sample.ipynb notebook 文件。

如何操作...

按照以下步骤使用 papermill 命令行界面远程执行 Jupyter notebook：

  1. 首先，我们从本章的代码库中打开示例 notebook 文件 sample.ipynb。该 notebook 包含三个代码单元，其中包含以下代码：

    import matplotlib.pyplot as plt

    from numpy.random import default_rng

    rng = default_rng(12345)

    uniform_data = rng.uniform(-5, 5, size=(2, 100))

    fig, ax = plt.subplots(tight_layout=True)

    ax.scatter(uniform_data[0, :], uniform_data[1, :])

    ax.set(title="Scatter plot", xlabel="x", ylabel="y")

  2. 接下来，我们在终端中打开包含 Jupyter notebook 的文件夹，并使用以下命令：

    papermill --kernel python3 sample.ipynb output.ipynb

  3. 现在，我们打开输出文件 output.ipynb，该文件应该包含已经更新的 notebook，其中包含执行代码的结果。最终生成的散点图如下所示：

图 10.3 - 在 Jupyter notebook 内生成的随机数据的散点图

图 10.3 - 在 Jupyter notebook 内生成的随机数据的散点图

请注意，papermill 命令的输出是一个全新的 notebook，它复制了原始代码和文本内容，并填充了运行命令后的输出。这对于“冻结”用于生成结果的确切代码非常有用。

工作原理...

papermill 包提供了一个简单的命令行界面，它解释并执行 Jupyter notebook，并将结果存储在新的 notebook 文件中。在此方法中，我们将第一个参数——输入的 notebook 文件——设为 sample.ipynb，第二个参数——输出的 notebook 文件——设为 output.ipynb。然后，该工具执行 notebook 中的代码并生成输出。notebook 的文件格式会跟踪最后一次运行的结果，因此这些结果会被添加到输出 notebook 中并存储在指定位置。在这个示例中，我们将其存储为一个简单的本地文件，但 papermill 也可以将其存储在云端位置，如 Amazon Web Services (AWS) S3 存储或 Azure 数据存储。

在 步骤 2 中，我们在使用 papermill 命令行界面时添加了 --kernel python3 选项。该选项允许我们指定用于执行 Jupyter notebook 的内核。如果 papermill 尝试使用与编写 notebook 时不同的内核执行 notebook，可能会导致错误，因此此选项可能是必要的。可以通过在终端中使用以下命令来查看可用的内核列表：

jupyter kernelspec list

如果执行 notebook 时遇到错误，您可以尝试更换为其他内核。

还有更多内容...

Papermill 还提供了一个 Python 接口，使得你可以在 Python 应用程序中执行 notebook。这对于构建需要在外部硬件上执行长时间计算并且结果需要存储在云中的 Web 应用程序可能很有用。它还可以向 notebook 提供参数。为此，我们需要在 notebook 中创建一个标记为参数的块，并设置默认值。然后，可以通过命令行接口使用-p标志提供更新的参数，后跟参数名称和相应的值。

验证数据

数据通常以原始形式呈现，可能包含异常、错误或格式不正确的数据，这显然会给后续的处理和分析带来问题。通常，在处理管道中构建一个验证步骤是个好主意。幸运的是，Cerberus 包为 Python 提供了一个轻量级且易于使用的验证工具。

对于验证，我们必须定义一个模式，它是数据应该是什么样子以及应该执行哪些检查的技术描述。例如，我们可以检查类型并设定最大和最小值的边界。Cerberus 验证器还可以在验证步骤中执行类型转换，这使我们能够将直接从 CSV 文件加载的数据插入到验证器中。

在这个示例中，我们将学习如何使用 Cerberus 验证从 CSV 文件加载的数据。

准备就绪

对于这个示例，我们需要从 Python 标准库中导入csv模块（docs.python.org/3/library/csv.html），以及 Cerberus 包：

import csv
import cerberus

我们还需要此章节中的sample.csv文件，来自代码仓库（github.com/PacktPublishing/Applying-Math-with-Python/tree/master/Chapter%2010）。

如何操作...

在接下来的步骤中，我们将验证一组已从 CSV 加载的数据，使用 Cerberus 包进行验证：

  1. 首先，我们需要构建一个描述我们期望数据的模式。为此，我们必须定义一个简单的浮点数模式：

    float_schema = {"type": "float", "coerce": float, 

        "min": -1.0, "max": 1.0}

  2. 接下来，我们为单个项构建模式。这些将是我们数据的行：

    item_schema = {

        "type": "dict",

        "schema": {

            "id": {"type": "string"},

            "number": {"type": "integer",

            "coerce": int},

        "lower": float_schema,

        "upper": float_schema,

        }

    }

  3. 现在，我们可以为整个文档定义一个模式，该模式将包含一个项目列表：

    schema = {

        "rows": {

            "type": "list",

            "schema": item_schema

        }

    }

  4. 接下来，我们创建一个带有我们刚定义的模式的Validator对象：

    validator = cerberus.Validator(schema)

  5. 然后，我们使用csv模块中的DictReader加载数据：

    with open("sample.csv") as f:

        dr = csv.DictReader(f)

        document = {"rows": list(dr)}

  6. 接下来，我们使用validate方法在validator上验证文档：

    validator.validate(document)

  7. 然后，我们从validator对象中获取验证过程中的错误：

    errors = validator.errors["rows"][0]

  8. 最后，我们可以打印出现的任何错误信息：

    for row_n, errs in errors.items():

                print(f"row {row_n}: {errs}")

错误消息的输出如下：

row 11: [{'lower': ['min value is -1.0']}]
row 18: [{'number': ['must be of integer type',      "field 'number' cannot be coerced: invalid literal for int() with base 10: 'None'"]}]
row 32: [{'upper': ['min value is -1.0']}]
row 63: [{'lower': ['max value is 1.0']}]

这已识别出四行不符合我们设定的模式，这限制了“lower”和“upper”中的浮动值仅限于-1.0到1.0之间。

它是如何工作的...

我们创建的架构是所有需要检查数据标准的技术描述。通常，它会被定义为一个字典，字典的键是项目名称，值是一个包含属性（如类型或值范围）的字典。例如，在第 1 步中，我们为浮点数定义了一个架构，限制这些数字的值在-1 和 1 之间。请注意，我们包含了coerce键，它指定了在验证过程中值应该转换为的类型。这使得我们可以传入从 CSV 文档中加载的数据，尽管它仅包含字符串，而无需担心其类型。

validator对象负责解析文档，以便验证并检查它们包含的数据是否符合架构中描述的所有标准。在这个食谱中，我们在创建validator对象时提供了架构。然而，我们也可以将架构作为第二个参数传递给validate方法。错误信息被存储在一个嵌套的字典中，字典的结构与文档的结构相对应。

使用 Cython 加速代码

Python 经常被批评为一种慢编程语言——这一点常常被辩论。许多批评可以通过使用具有 Python 接口的高性能编译库来解决，比如科学 Python 栈，从而大大提高性能。然而，有些情况下，我们无法避免 Python 不是一种编译语言的事实。在这些（相对罕见）情况下，改善性能的一种方式是编写 C 扩展（甚至将代码完全用 C 重写），以加速关键部分。这肯定会让代码运行得更快，但可能会让维护该包变得更加困难。相反，我们可以使用 Cython，它是 Python 语言的扩展，通过将 Python 代码转换为 C 并编译，从而实现极大的性能提升。

例如，我们可以考虑一些用来生成曼德尔布罗集图像的代码。为了对比，我们假设纯 Python 代码作为起点，如下所示：

# mandelbrot/python_mandel.py
import numpy as np
def in_mandel(cx, cy, max_iter):
    x = cx
    y = cy
    for i in range(max_iter):
        x2 = x**2
        y2 = y**2
        if (x2 + y2) >= 4:
            return i
        y = 2.0*x*y + cy
        x = x2 - y2 + cx
    return max_iter
def compute_mandel(N_x, N_y, N_iter):
    xlim_l = -2.5
    xlim_u = 0.5
    ylim_l = -1.2
    ylim_u = 1.2
    x_vals = np.linspace(xlim_l, xlim_u,
        N_x, dtype=np.float64)
y_vals = np.linspace(ylim_l, ylim_u,
        N_y, dtype=np.float64)
    height = np.empty((N_x, N_y), dtype=np.int64)
    for i in range(N_x):
        for j in range(N_y):
        height[i, j] = in_mandel(
		    x_vals[i], y_vals[j], N_iter)
    return height

这段代码在纯 Python 中相对较慢的原因相当明显：嵌套的循环。为了演示目的，假设我们无法使用 NumPy 向量化这段代码。初步测试表明，使用这些函数生成曼德尔布罗集图像，使用 320 × 240 的点和 255 步大约需要 6.3 秒。你的测试时间可能会有所不同，取决于你的系统。

在这个食谱中，我们将使用 Cython 大幅提升前述代码的性能，以生成曼德尔布罗集图像。

准备工作

对于这个食谱，我们需要安装 NumPy 包和 Cython 包。你还需要在系统上安装一个 C 编译器，如 GCC。例如，在 Windows 上，你可以通过安装 MinGW 来获得 GCC 版本。

如何做到...

按照以下步骤使用 Cython 大幅提高生成曼德尔布罗特集图像代码的性能：

  1. 在mandelbrot文件夹中开始一个名为cython_mandel.pyx的新文件。在此文件中，我们将添加一些简单的导入和类型定义：

    # mandelbrot/cython_mandel.pyx

    import numpy as np

    cimport numpy as np

    cimport cython

    ctypedef Py_ssize_t Int

    ctypedef np.float64_t Double

  2. 接下来，我们使用 Cython 语法定义in_mandel例程的新版本。在该例程的前几行添加一些声明：

    cdef int in_mandel(Double cx, Double cy, int max_iter):

        cdef Double x = cx

        cdef Double y = cy

        cdef Double x2, y2

        cdef Int i

  3. 其余部分与 Python 版本的函数完全相同：

        for i in range(max_iter):

            x2 = x**2

            y2 = y**2

            if (x2 + y2) >= 4:

                return i

            y = 2.0*x*y + cy

            x = x2 - y2 + cx

        return max_iter

  4. 接下来，我们定义compute_mandel函数的新版本。我们向该函数添加了两个来自 Cython 包的装饰器：

    @cython.boundscheck(False)

    @cython.wraparound(False)

    def compute_mandel(int N_x, int N_y, int N_iter):

  5. 然后，我们定义常量，就像在原始例程中一样：

        cdef double xlim_l = -2.5

        cdef double xlim_u = 0.5

        cdef double ylim_l = -1.2

        cdef double ylim_u = 1.2

  6. 我们使用与 Python 版本完全相同的方式调用 NumPy 包中的linspace和empty例程。唯一的不同是，我们声明了i和j变量，它们是Int类型：

        cdef np.ndarray x_vals = np.linspace(xlim_l,

            xlim_u, N_x, dtype=np.float64)

        cdef np.ndarray y_vals = np.linspace(ylim_l,

            ylim_u, N_y, dtype=np.float64)

        cdef np.ndarray height = np.empty(

            (N_x, N_y),dtype=np.int64)

        cdef Int i, j

  7. 定义的其余部分与 Python 版本完全相同：

        for i in range(N_x):

            for j in range(N_y):

                height[i, j] = in_mandel(

                    xx_vals[i], y_vals[j], N_iter)

            return height

  8. 接下来，我们在mandelbrot文件夹中创建一个新的文件，命名为setup.py，并在文件顶部添加以下导入：

    # mandelbrot/setup.py

    import numpy as np

    from setuptools import setup, Extension

    from Cython.Build import cythonize

  9. 之后，我们定义一个扩展模块，源文件指向原始的python_mandel.py文件。将此模块的名称设置为hybrid_mandel：

    hybrid = Extension(

        "hybrid_mandel",

        sources=["python_mandel.py"],

        include_dirs=[np.get_include()],

        define_macros=[("NPY_NO_DEPRECATED_API",

            "NPY_1_7_API_VERSION")]

    )

  10. 现在，我们定义第二个扩展模块，源文件设置为刚才创建的cython_mandel.pyx文件：

    cython = Extension(

        "cython_mandel",

        sources=["cython_mandel.pyx"],

        include_dirs=[np.get_include()],

        define_macros=[("NPY_NO_DEPRECATED_API",

            "NPY_1_7_API_VERSION")]

    )

  11. 接下来，我们将这两个扩展模块添加到一个列表中，并调用setup例程来注册这些模块：

    extensions = [hybrid, cython]

    setup(

        ext_modules = cythonize(

            extensions, compiler_directives={

    		    "language_level": "3"}),

    )

  12. 在mandelbrot文件夹中创建一个名为__init__.py的新空文件，使其成为一个可以导入到 Python 中的包。

  13. 打开mandelbrot文件夹中的终端，使用以下命令构建 Cython 扩展模块：

    python3.8 setup.py build_ext --inplace

  14. 现在，开始一个名为run.py的新文件，并添加以下import语句：

    # run.py

    from time import time

    from functools import wraps

    import matplotlib.pyplot as plt

  15. 从我们定义的每个模块中导入不同的compute_mandel例程：python_mandel为原始版本；hybrid_mandel为 Cython 化的 Python 代码；cython_mandel为编译后的纯 Cython 代码：

    from mandelbrot.python_mandel import compute_mandel

                as compute_mandel_py

    from mandelbrot.hybrid_mandel import compute_mandel

                as compute_mandel_hy

    from mandelbrot.cython_mandel import compute_mandel

                as compute_mandel_cy

  16. 定义一个简单的计时器装饰器，我们将用它来测试例程的性能：

    def timer(func, name):

    	@wraps(func)

    	def wrapper(*args, **kwargs):

    		t_start = time()

    		val = func(*args, **kwargs)

    		t_end = time()

    		print(f"Time taken for {name}:

    			{t_end - t_start}")

    		return val

    	return wrapper

  17. 将timer装饰器应用到每个导入的例程，并为测试定义一些常量：

    mandel_py = timer(compute_mandel_py, "Python")

    mandel_hy = timer(compute_mandel_hy, "Hybrid")

    mandel_cy = timer(compute_mandel_cy, "Cython")

    Nx = 320

    Ny = 240

    steps = 255

  18. 使用我们之前设置的常量运行每个已装饰的例程。将最后一次调用（Cython 版本）的输出记录在vals变量中：

    mandel_py(Nx, Ny, steps)

    mandel_hy(Nx, Ny, steps)

    vals = mandel_cy(Nx, Ny, steps)

  19. 最后，绘制 Cython 版本的输出，以检查该例程是否正确计算曼德尔布罗特集：

    fig, ax = plt.subplots()

    ax.imshow(vals.T, extent=(-2.5, 0.5, -1.2, 1.2))

    plt.show()

运行run.py文件将会打印每个例程的执行时间到终端，如下所示：

Time taken for Python: 11.399756908416748
Time taken for Hybrid: 10.955225229263306
Time taken for Cython: 0.24534869194030762

注意

这些时间不如第一版那样好，这可能是由于作者的 PC 上 Python 安装的方式。你的时间可能会有所不同。

曼德尔布罗特集的图像可以在下图中看到：

图 10.4 - 使用 Cython 代码计算的曼德尔布罗特集图像

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

图 10.4 - 使用 Cython 代码计算的曼德尔布罗特集图像

这就是我们对曼德尔布罗集合的预期效果。细节部分在边界处有所显示。

它是如何工作的...

这个食谱中有很多内容需要说明，因此我们先从解释整体过程开始。Cython 将用 Python 语言扩展编写的代码编译成 C 代码，然后生成一个可以导入到 Python 会话中的 C 扩展库。实际上，你甚至可以使用 Cython 将普通的 Python 代码直接编译成扩展，尽管结果不如使用修改后的语言时好。食谱中的前几步定义了用修改后的语言编写的 Python 代码的新版本（保存为 .pyx 文件），其中除了常规的 Python 代码外，还包含了类型信息。为了使用 Cython 构建 C 扩展，我们需要定义一个设置文件，然后创建一个文件运行它来生成结果。

Cython 编译后的最终版本比其 Python 等效版本运行得要快得多。Cython 编译后的 Python 代码（我们在本食谱中称之为混合代码）比纯 Python 代码稍微快一些。这是因为生成的 Cython 代码仍然需要与 Python 对象进行交互，且必须考虑所有相关问题。通过在 .pyx 文件中的 Python 代码中添加类型信息，我们开始看到性能有了显著的提升。这是因为 in_mandel 函数现在有效地被定义为一个 C 级别的函数，不再与 Python 对象交互，而是直接操作原始数据类型。

Cython 代码与 Python 等效代码之间有一些小但非常重要的差异。在 步骤 1 中，你可以看到我们照常导入了 NumPy 包，但我们还使用了 cimport 关键字将一些 C 级别的定义引入作用域。在 步骤 2 中，我们在定义 in_mandel 函数时使用了 cdef 关键字，而不是 def 关键字。这意味着 in_mandel 函数被定义为一个 C 级别的函数，不能从 Python 级别调用，这样在调用该函数时（这会发生很多次）就节省了大量开销。

关于该函数定义的唯一其他实际差异是签名中和函数的前几行中包含了一些类型声明。我们在这里应用的两个装饰器禁用了访问列表（数组）元素时的边界检查。boundscheck 装饰器禁用了检查索引是否有效（即是否在 0 和数组大小之间），而 wraparound 装饰器禁用了负索引。这两个装饰器都能在执行过程中带来适度的速度提升，尽管它们会禁用 Python 内置的一些安全功能。在这个食谱中，禁用这些检查是可以接受的，因为我们正在遍历数组的有效索引。

设置文件是我们告诉 Python（因此也告诉 Cython）如何构建 C 扩展的地方。Cython 中的 cythonize 例程是关键，它触发了 Cython 构建过程。在 步骤 9 和 步骤 10 中，我们使用 setuptools 中的 Extension 类定义了扩展模块，以便为构建定义一些额外的细节；具体来说，我们为 NumPy 编译设置了一个环境变量，并添加了 NumPy C 头文件的 include 文件。这是通过 Extension 类的 define_macros 关键字参数完成的。我们在 步骤 13 中使用的终端命令使用 setuptools 构建 Cython 扩展，并且添加 --inplace 标志意味着编译后的库将被添加到当前目录，而不是放置在集中位置。这对于开发来说非常方便。

运行脚本非常简单：从每个已定义的模块中导入例程——其中有两个实际上是 C 扩展模块——并测量它们的执行时间。我们需要在导入别名和例程名称上稍作创意，以避免冲突。

还有更多……

Cython 是一个强大的工具，可以提高代码某些方面的性能。然而，在优化代码时，你必须始终谨慎，明智地分配时间。可以使用 Python 标准库中提供的配置文件（如 cProfile）来找到代码中性能瓶颈发生的地方。在这个案例中，性能瓶颈的位置相当明显。Cython 是解决这个问题的一个好办法，因为它涉及到在（双重）for 循环中重复调用一个函数。然而，它并不是解决所有性能问题的万能方法，通常情况下，重构代码使其能够利用高性能库，往往能显著提高代码的性能。

Cython 与 Jupyter Notebook 紧密集成，可以无缝地在笔记本的代码块中使用。当 Cython 使用 Anaconda 发行版安装时，它也包含在 Anaconda 中，因此不需要额外的设置就可以在 Jupyter 笔记本中使用 Cython。

还有其他替代方案可以将 Python 代码编译为机器码。例如，Numba 包 (numba.pydata.org/) 提供了一个 即时编译（JIT）编译器，通过简单地在特定函数上添加装饰器，在运行时优化 Python 代码。Numba 旨在与 NumPy 和其他科学 Python 库一起使用，也可以用来利用 GPU 加速代码。

通过 pyjion 包（www.trypyjion.com/）还可以使用一个通用的 JIT 编译器。这可以在各种场景中使用，不像 Numba 库主要用于数值代码。jax 库中也有一个内置的 JIT 编译器，如在第三章所讨论的，但它也仅限于数值代码。

使用 Dask 进行分布式计算

Dask 是一个用于在多个线程、进程甚至计算机之间分布式计算的库，旨在有效地进行大规模计算。即使在单台笔记本电脑上工作，它也能大幅提升性能和吞吐量。Dask 提供了 Python 科学计算栈中大多数数据结构的替代品，比如 NumPy 数组和 Pandas DataFrame。这些替代品具有非常相似的接口，但底层是为分布式计算而构建的，可以在多个线程、进程或计算机之间共享。在许多情况下，切换到 Dask 就像改变 import 语句那么简单，可能还需要添加几个额外的方法调用来启动并发计算。

在这个实例中，我们将学习如何使用 Dask 在 DataFrame 上进行一些简单的计算。

准备工作

对于这个实例，我们需要从 Dask 包中导入 dataframe 模块。按照 Dask 文档中的约定，我们将以 dd 别名导入此模块：

import dask.dataframe as dd

我们还需要本章代码库中的 sample.csv 文件。

如何操作...

按照以下步骤使用 Dask 对 DataFrame 对象执行一些计算：

  1. 首先，我们需要将数据从 sample.csv 加载到 Dask DataFrame 中。number 列的类型设置为 "object"，因为否则 Dask 的类型推断会失败（因为此列包含 None，但其余部分为整数）：

    data = dd.read_csv("sample.csv", dtype={

        "number":"object"})

  2. 接下来，我们对 DataFrame 的列执行标准计算：

    sum_data = data.lower + data.upper

    print(sum_data)

与 Pandas DataFrame 不同，结果不是一个新的 DataFrame。print 语句为我们提供了以下信息：

Dask Series Structure:
npartitions=1
             float64
                               ...
dtype: float64
Dask Name: add, 4 graph layers

  1. 要实际获取结果，我们需要使用 compute 方法：

    result = sum_data.compute()

    print(result.head())

结果现在如预期所示：

0      -0.911811
1       0.947240
2      -0.552153
3      -0.429914
4       1.229118
dtype:  float64

  1. 我们计算最后两列的均值，方法和使用 Pandas DataFrame 完全相同，只是需要调用 compute 方法来执行计算：

    means = data[["lower", "upper"]].mean().compute()

    print(means)

结果如我们所预期的那样打印出来：

lower -0.060393
upper -0.035192
dtype: float64

它是如何工作的...

Dask 为计算构建了一个任务图，该图描述了需要在数据集合上执行的各种操作和计算之间的关系。这将计算步骤拆解开来，以便可以在不同的工作节点上按照正确的顺序进行计算。这个任务图会被传递给调度器，调度器将实际任务发送到工作节点进行执行。Dask 提供了几种不同的调度器：同步、线程化、进程化和分布式。调度器的类型可以在调用 compute 方法时选择，或者全局设置。如果没有指定，Dask 会选择一个合理的默认值。

同步、线程化和进程化调度器适用于单台机器，而分布式调度器则用于集群。Dask 允许你以相对透明的方式在调度器之间进行切换，尽管对于小任务来说，由于设置更复杂的调度器所带来的开销，你可能不会获得任何性能提升。

compute 方法是本食谱的关键。通常会对 Pandas DataFrame 执行计算的方法，现在仅仅是设置一个通过 Dask 调度器执行的计算。直到调用 compute 方法，计算才会开始。这类似于 Future（如来自 asyncio 标准库的 Future）的方式，它作为异步函数调用结果的代理，直到计算完成时，才会返回实际结果。

还有更多...

Dask 提供了用于 NumPy 数组的接口，以及在本食谱中展示的 DataFrame 接口。还有一个名为 dask_ml 的机器学习接口，提供与 scikit-learn 包类似的功能。一些外部包，如 xarray，也有 Dask 接口。Dask 还可以与 GPU 协同工作，以进一步加速计算并从远程源加载数据，这在计算分布在集群中的时候非常有用。

编写可重复的数据科学代码

科学方法的基本原则之一是，结果应该是可重复的并且可以独立验证的。遗憾的是，这一原则常常被低估，反而更重视“新颖”的想法和结果。作为数据科学的从业者，我们有责任尽力使我们的分析和结果尽可能具有可重复性。

由于数据科学通常完全在计算机上进行——也就是说，它通常不涉及测量中的仪器误差——有些人可能会认为所有的数据科学工作本质上都是可重复的。但事实并非如此。在使用随机化的超参数搜索或基于随机梯度下降的优化时，容易忽略一些简单的细节，比如种子设置（参见第三章）。此外，一些更微妙的非确定性因素（例如使用线程或多进程）如果不加以注意，可能会显著改变结果。

在这个例程中，我们将展示一个基本数据分析管道的示例，并实施一些基本步骤，以确保您可以重复结果。

准备工作

对于这个例程，我们将需要 NumPy 包，通常以np导入，Pandas 包，以pd导入，Matplotlib 的pyplot接口，以plt导入，以及从scikit-learn包中导入以下内容：

from sklearn.metrics import ConfusionMatrixDisplay, accuracy_score
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeClassifier

我们将模拟我们的数据（而不是从其他地方获取数据），因此需要使用一个具有种子值的默认随机数生成器实例（以确保可重复性）：

rng = np.random.default_rng(12345)

为了生成数据，我们定义了以下例程：

def get_data():
	permute = rng.permutation(200)
	data = np.vstack([
		rng.normal((1.0, 2.0, -3.0), 1.0,
		size=(50, 3)),
		rng.normal((-1.0, 1.0, 1.0), 1.0,
		size=(50, 3)),
		rng.normal((0.0, -1.0, -1.0), 1.0,
		size=(50, 3)),
		rng.normal((-1.0, -1.0, -2.0), 1.0,
		size=(50, 3))
		])
	labels = np.hstack(
		[[1]*50, [2]*50, [3]*50,[4]*50])
	X = pd.DataFrame(
		np.take(data, permute, axis=0),
		columns=["A", "B", "C"])
	y = pd.Series(np.take(labels, permute, axis=0))
	return X, y

我们使用这个函数来代替其他将数据加载到 Python 中的方法，例如从文件中读取或从互联网上下载。

如何实现……

请按照以下步骤创建一个非常简单且可重复的数据科学管道：

  1. 首先，我们需要使用之前定义的get_data例程“加载”数据：

    data, labels = get_data()

  2. 由于我们的数据是动态获取的，最好将数据与我们生成的任何结果一起存储。

    data.to_csv("data.csv")

    labels.to_csv("labels.csv")

  3. 现在，我们需要使用scikit-learn中的train_test_split例程将数据分为训练集和测试集。我们将数据按 80/20 的比例进行划分，并确保设置了随机状态，以便可以重复这个过程（尽管我们将在下一步保存索引以供参考）：

    X_train, X_test, y_train, y_test = train_test_split(

        data,labels, test_size=0.2, random_state=23456)

  4. 现在，我们确保保存训练集和测试集的索引，以便我们确切知道每个样本中的观测值。我们可以将这些索引与步骤 2中存储的数据一起使用，以便稍后完全重建这些数据集：

    X_train.index.to_series().to_csv("train_index.csv",

        index=False, header=False)

    X_test.index.to_series().to_csv("test_index.csv",

        index=False, header=False)

  5. 现在，我们可以设置并训练分类器。在这个示例中，我们使用一个简单的DecisionTreeClassifier，但这个选择并不重要。由于训练过程涉及一些随机性，请确保将random_state关键字参数设置为种子值，以便控制这种随机性：

    classifier = DecisionTreeClassifier(random_state=34567)

    classifer.fit(X_train, y_train)

  6. 在继续之前，最好先收集一些关于训练模型的信息，并将其与结果一起存储。不同模型的有趣信息会有所不同。对于这个模型，特征重要性信息可能很有用，因此我们将其记录在一个 CSV 文件中：

    feat_importance = pd.DataFrame(

    	classifier.feature_importances_,

    	index=classifier.feature_names_in_,

    	columns=["Importance"])

    feat_importance.to_csv("feature_importance.csv")

  7. 现在，我们可以继续检查模型的表现。我们将在训练数据和测试数据上评估模型，稍后我们会将其与真实标签进行比较：

    train_predictions = classifier.predict(X_train)

    test_predictions = classifier.predict(X_test)

  8. 始终保存此类预测任务（或回归任务，或其他将以某种方式成为报告一部分的最终结果）的结果。我们首先将它们转换为Series对象，以确保索引设置正确：

    pd.Series(train_predictions,index=X_train.index,

        name="Predicted label").to_csv(

    		"train_predictions.csv")

    pd.Series(test_predictions,index=X_test.index,

        name="Predicted label").to_csv(

    		"test_predictions.csv")

  9. 最后，我们可以生成任何图形或度量，以帮助我们决定如何继续进行分析。在这里，我们将为训练和测试队列分别生成一个混淆矩阵图，并打印出一些准确度总结评分：

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

    ax1.set_title("Confusion matrix for training data")

    ax2.set_title("Confusion matrix for test data")

    ConfusionMatrixDisplay.from_predictions(

    	y_train, train_predictions,

    	ax=ax1 cmap="Greys", colorbar=False)

    ConfusionMatrixDisplay.from_predictions(

    	y_test, test_predictions,

    	ax=ax2 cmap="Greys", colorbar=False)

    print(f"Train accuracy {accuracy_score(y_train, train_predictions)}",

    	f"Test accuracy {accuracy_score(y_test, test_predictions)}",

    	sep="\n")

    # Train accuracy 1.0

    # Test accuracy 0.65

结果的混淆矩阵见图 10.5：

图 10.5 - 简单分类任务的混淆矩阵

图 10.5 - 简单分类任务的混淆矩阵

这个例子的测试结果并不出色，这不应该令人惊讶，因为我们没有花时间选择最合适的模型或进行调优，并且我们的样本量相当小。为这些数据生成一个准确的模型并不是目标。在当前目录中（脚本运行所在的目录），应该会有一些新的 CSV 文件，包含我们写入磁盘的所有中间数据：data.csv、labels.csv、train_index.csv、test_index.csv、feature_importance.csv、train_predictions.csv 和 test_predictions.csv。

它是如何工作的…

在可重复性方面没有明确的正确答案，但肯定有错误的答案。我们这里只触及了如何使代码更具可重复性的一些想法，但还有很多可以做的事情。（参见更多内容…）。在这个过程中，我们实际上更专注于存储中间值和结果，而不是其他任何东西。这一点常常被忽视，大家更倾向于生成图表和图形——因为这些通常是展示结果的方式。然而，我们不应该为了更改图表的样式而重新运行整个流程。存储中间值可以让你审计流程中的各个部分，检查你做的事情是否合理和适当，并确保你能从这些中间值中重现结果。

一般来说，数据科学流程包含五个步骤：

  1. 数据获取

  2. 数据预处理和特征选择

  3. 模型和超参数调优

  4. 模型训练

  5. 评估和结果生成

在本教程中，我们将数据获取步骤替换为一个随机生成数据的函数。如引言中所述，这一步通常会涉及从磁盘加载数据（来自 CSV 文件或数据库）、从互联网下载数据，或直接从测量设备采集数据。我们缓存了数据获取的结果，因为我们假设这是一个开销较大的操作。当然，这并非总是如此；如果你直接从磁盘加载所有数据（例如通过 CSV 文件），那么显然不需要存储这份数据的第二份副本。然而，如果你通过查询一个大型数据库生成数据，那么存储数据的平面副本将大大提高你在数据管道上迭代的速度。

我们的预处理仅包括将数据分割成训练集和测试集。再次说明，我们在这一步之后存储了足够的数据，以便稍后可以独立地重建这些数据集——我们只存储了与每个数据集对应的 ID。由于我们已经存储了这些数据集，因此在 train_test_split 函数中种子随机数并非绝对必要，但通常来说这是一个好主意。如果你的预处理涉及更为密集的操作，你可能会考虑缓存处理过的数据或在数据管道中使用的生成特征（我们稍后将更详细地讨论缓存）。如果你的预处理步骤涉及从数据的列中选择特征，那么你应该绝对保存这些选择的特征到磁盘，并与结果一起存储。

我们的模型非常简单，且没有任何（非默认）超参数。如果你做过超参数调整，你应该将这些参数以及任何可能需要用来重建模型的元数据进行存储。存储模型本身（通过序列化或其他方式）是有用的，但请记住，序列化后的模型可能无法被其他方读取（例如，如果他们使用的是不同版本的 Python）。

你应该始终存储来自模型的数值结果。当你检查后续运行时结果是否相同，比较图表和其他摘要图形几乎是不可能的。此外，这样做可以让你在之后快速重新生成图形或数值，如果有需要的话。例如，如果你的分析涉及二元分类问题，那么存储用于生成接收者操作特征（ROC）曲线的值是个好主意，即使你已经绘制了 ROC 曲线的图形并报告了曲线下的面积。

还有更多……

我们在这里还没有讨论很多内容。首先，让我们处理一个显而易见的问题。Jupyter notebooks 是一种常用的数据科学工作流工具。这没问题，但用户应该了解，这种格式有几个缺点。首先，最重要的一点是，Jupyter notebooks 可以无序运行，后面的单元格可能会依赖于前面的单元格，这些依赖关系可能不是很直观。为了解决这个问题，确保你每次都在一个干净的内核上完整地运行整个 notebook，而不是仅仅在当前内核中重新运行每个单元格（例如，使用 Executing a Jupyter notebook as a script 这一配方中的 Papermill 工具）。其次，notebook 中存储的结果可能与代码单元格中编写的代码不一致。这种情况发生在 notebook 被运行过并且代码在事后被修改，但没有重新运行的情况下。一个好的做法是，保留一份没有任何存储结果的主副本，然后创建多个副本，其中包含结果，并且这些副本不再进行修改。最后，Jupyter notebooks 通常在一些难以正确缓存中间步骤结果的环境中执行。这部分问题由 notebook 内部的缓存机制部分解决，但它并不总是完全透明的。

现在，让我们讨论两个与可重复性相关的一般性问题：配置和缓存。配置是指用于控制工作流设置和执行的值的集合。在这个配方中，我们没有明显的配置值，除了在 train_test_split 例程中使用的随机种子和模型（以及数据生成，暂时不考虑这些），以及训练/测试集拆分时所取的百分比。这些值是硬编码在配方中的，但这可能不是最好的做法。至少，我们希望能够记录每次运行分析时使用的配置。理想情况下，配置应该从文件中加载（仅加载一次），然后在工作流运行之前被最终确定并缓存。这意味着，完整的配置应该从一个或多个来源（配置文件、命令行参数或环境变量）加载，汇总为一个“真实来源”，然后将其序列化为机器可读和人类可读的格式，如 JSON，并与结果一起存储。这样你就可以确切知道是使用了什么配置生成了这些结果。

缓存是存储中间结果的过程，以便在随后的运行中可以重用，从而减少运行时间。在本食谱中，我们确实存储了中间结果，但我们没有建立机制来重用已存储的数据（如果它存在且有效）。这是因为实际的检查和加载缓存值的机制较为复杂，并且有些依赖于具体的配置。由于我们的项目非常小，所以缓存值未必有意义。然而，对于有多个组件的大型项目来说，缓存确实能带来差异。在实现缓存机制时，您应该建立一个系统来检查缓存是否有效，例如，可以使用代码文件及其依赖的任何数据源的 SHA-2 哈希值。

在存储结果时，通常的好做法是将所有结果一起存储在一个带有时间戳的文件夹或类似的地方。我们在本食谱中没有这样做，但其实很容易实现。例如，使用标准库中的datetime和pathlib模块，我们可以轻松创建一个用于存储结果的基础路径：

from pathlib import Path
from datetime import datetime
RESULTS_OUT = Path(datetime.now().isoformat())
...
results.to_csv(RESULTS_OUT / "name.csv")

如果您使用多进程并行运行多个分析任务，必须小心，因为每个新进程都会生成一个新的RESULTS_OUT全局变量。更好的选择是将其纳入配置过程，这样用户还可以自定义输出路径。

除了我们迄今为止讨论的脚本中的实际代码外，在项目层面上还有很多可以做的事情，以提高代码的可重现性。第一步，也是最重要的一步，是尽可能使代码可用，这包括指定代码可以共享的许可证（如果可以的话）。此外，好的代码应该足够健壮，以便可以用于分析多个数据集（显然，这些数据应该与最初使用的数据类型相同）。同样重要的是使用版本控制（如 Git、Subversion 等）来跟踪更改。这也有助于将代码分发给其他用户。最后，代码需要有良好的文档说明，并且理想情况下应有自动化测试，以检查管道在示例数据集上的预期表现。

参见...

以下是一些关于可重现编码实践的额外信息来源：

这结束了本书的第十章，也是最后一章。请记住，我们才刚刚触及了使用 Python 做数学时可能实现的表面，您应该阅读本书中提到的文档和资源，以了解这些包和技术能够实现的更多信息。

avatar
文章
阅读
粉丝