MIT 6.0001：Python 计算和编程入门第三版（三）

原文：zh.z-lib.gs/md5/b81f9400901fb07c6e4e456605c4cd1f

译者：飞龙

协议：CC BY-NC-SA 4.0

第十五章：动态规划

动态规划是理查德·贝尔曼在 1950 年代发明的。不要试图从其名称中推测任何关于该技术的内容。正如贝尔曼所描述的，名称“动态规划”被选择是为了掩盖政府赞助者“我实际上是在做数学… [动态规划这个短语]是连国会议员都无法反对的东西。” ⁹⁴

动态规划是一种高效解决具有重叠子问题和最优子结构特征的问题的方法。幸运的是，许多优化问题具有这些特征。

如果一个全局最优解可以通过结合局部子问题的最优解来找到，则该问题具有最优子结构。我们已经查看了许多此类问题。例如，归并排序利用了列表可以通过先排序子列表然后合并解决方案来进行排序这一事实。

如果一个问题的最优解涉及多次解决同一问题，则该问题具有重叠子问题。归并排序并不具备这个特性。尽管我们多次进行归并，但每次归并的列表都是不同的。

这并不是显而易见的，但 0/1 背包问题同时具备这两个特性。首先，我们稍作偏离，来看一个最优子结构和重叠子问题更加明显的问题。

15.1 重新审视斐波那契数列

在第四章中，我们查看了斐波那契函数的一个简单递归实现：

def fib(n):
    """Assumes n is an int >= 0
       Returns Fibonacci of n"""
    if n == 0 or n == 1:
        return 1
    else:
        return fib(n-1) + fib(n-2)

尽管这个递归实现显然是正确的，但它非常低效。例如，尝试运行fib(120)，但不要等它完成。该实现的复杂度有点难以推导，但大致为O(fib(n))。也就是说，它的增长与结果值的增长成比例，而斐波那契序列的增长率是相当大的。例如，fib(120)是8,670,007,398,507,948,658,051,921。如果每个递归调用花费一纳秒，fib(120)将需要大约250,000年才能完成。

让我们试着找出为什么这个实现耗时如此之长。考虑到fib主体中的代码量很小，很明显问题在于fib调用自身的次数。例如，查看与调用fib(6)相关的调用树。

图 15-1 递归斐波那契的调用树

请注意，我们在反复计算相同的值。例如，fib被调用了三次，且每次调用都会引发四次额外的fib调用。想出一个好主意，即记录第一次调用返回的值，然后查找而不是每次都计算，并不需要天才。这是动态规划背后的关键思想。

动态规划有两种方法

• 记忆化从自顶向下解决原始问题。它从原始问题开始，将其分解为子问题，再将子问题分解为子问题，依此类推。每次解决一个子问题时，它都会将答案存储在表中。每次需要解决一个子问题时，它会首先尝试在表中查找答案。

• 表格法是一种自底向上的方法。它从最小的问题开始，并将这些问题的答案存储在表中。然后，它将这些问题的解决方案结合起来，以解决下一个最小的问题，并将这些答案存储在表中。

图 15-2 包含使用每种动态规划方法实现的斐波那契数列。函数 fib_memo 有一个参数 memo，用于跟踪它已经评估的数字。该参数有一个默认值，即空字典，因此 fib_memo 的客户不必担心提供 memo 的初始值。当 fib_memo 被调用时，若 n > 1，它会尝试在 memo 中查找 n。如果不存在（因为这是第一次以该值调用 fib_memo），则会引发异常。当这种情况发生时，fib_memo 使用正常的斐波那契递归，然后将结果存储在 memo 中。

函数 fib_tab 非常简单。它利用了所有斐波那契的子问题都是提前已知且容易以有用的顺序列举的事实。

图 15-2 使用备忘录实现斐波那契数列

如果你尝试运行 fib_memo 和 fib_tab，你会看到它们确实非常快：fib(120) 几乎瞬间返回。这个函数的复杂度是什么？fib_memo 针对从 0 到 n 的每个值调用 fib 恰好一次。因此，在假设字典查找可以在常数时间内完成的前提下，fib_memo(n) 的时间复杂度为 O(n)。而 fib_tab 的复杂度更明显也是 O(n)。

如果解决原始问题需要解决所有子问题，通常使用表格方法更好。它更简单，速度更快，因为它没有递归调用的开销，并且可以预先分配适当大小的表，而不是动态增长备忘录。如果仅需要解决某些子问题（这通常是情况），记忆化通常更高效。

手指练习：使用表格方法实现一个符合规范的动态规划解决方案

 def make_change(coin_vals, change):
    """coin_vals is a list of positive ints and coin_vals[0] = 1
       change is a positive int,
       return the minimum number of coins needed to have a set of
          coins the values of which sum to change. Coins may be used
          more than once. For example, make_change([1, 5, 8], 11)
          should return 3."""

15.2 动态规划与 0/1 背包问题

我们在第十四章中研究的优化问题之一是 0/1 背包问题。回想一下，我们研究了一种在 n log n 时间内运行的贪心算法，但不能保证找到最优解。我们还研究了一种保证找到最优解的暴力算法，但其运行时间为指数级。最后，我们讨论了这个问题在输入规模上本质上是指数级的。在最坏的情况下，如果不查看所有可能的答案，就无法找到最优解。

幸运的是，情况并不像看起来那么糟。动态规划提供了一种在合理时间内解决大多数 0/1 背包问题的实用方法。作为推导此类解决方案的第一步，我们从基于穷举枚举的指数解决方案开始。关键思想是通过构建一个根节点二叉树来探索满足重量约束的所有可能解的空间。

根节点二叉树是一个无环有向图，其中。

• 恰好有一个没有父节点的节点。这被称为根。

• 每个非根节点恰好有一个父节点。

• 每个节点最多有两个子节点。没有子节点的节点称为叶子。

0/1 背包问题的搜索树中的每个节点都带有一个四元组，表示背包问题的部分解。四元组的元素包括：

• 一组要取的物品。

• 尚未做出决策的物品列表。

• 要取的物品集合的总价值（这仅仅是一个优化，因为价值可以从集合中计算得出）。

• 背包中剩余的空间。（同样，这是一个优化，因为它仅仅是允许的重量与迄今为止所有物品重量的差值。）

树是自上而下构建的，从根节点开始。⁹⁵ 从待考虑的物品中选择一个元素。如果该物品可以放入背包中，就构建一个节点，反映选择这个物品的后果。根据约定，我们将该节点绘制为左子节点。右子节点显示选择不带这个物品的后果。然后递归应用此过程，直到背包满或没有更多物品可考虑。因为每条边代表一个决策（选择带入或不带入一个物品），这样的树被称为决策树。⁹⁶

图 15-3 是描述一组物品的表。

图 15-3 物品的价值和重量表。

图 15-4 是一个决策树，用于决定在假设背包最大重量为5的情况下应选择哪些物品。树的根节点（节点 0）标记为<{}, [a,b,c,d], 0, 5>，表示没有物品被选中，所有物品仍需考虑，已选物品的价值为0，并且仍有5的重量可用。节点1表示物品a已被选中，剩余需考虑的物品为[b,c,d]，已选物品的价值为6，而背包还可以容纳2磅。节点1没有左子节点，因为重量为3磅的物品b无法放入背包。

在图 15-4 中，每个节点前的冒号之前的数字表示节点生成的一种顺序。这种特定的顺序称为左优先深度优先。在每个节点，我们尝试生成左节点。如果不可能，则尝试生成右节点。如果也不可能，则向上回退一个节点（到父节点）并重复该过程。最终，我们生成了根节点的所有子孙，过程停止。当过程结束时，所有可以放入背包的物品组合都已生成，任何具有最大价值的叶节点都代表了一个最佳解。注意，对于每个叶节点，第二个元素要么是空列表（表示没有更多物品需要考虑），要么第四个元素为 0（表示背包没有剩余空间）。

图 15-4 背包问题的决策树

毫不奇怪（特别是如果你阅读了第十四章），深度优先树搜索的自然实现是递归的。图 15-5 包含这样的实现。它使用Item类及在图 14-2 中定义的函数。

函数max_val返回两个值，所选物品的集合及这些物品的总价值。它被调用时有两个参数，分别对应树中节点标签的第二个和第四个元素：

• to_consider。那些在树中较高节点（对应于递归调用栈中的早期调用）尚未考虑的物品。

• avail。仍然可用的空间量。

注意max_val的实现并不是构建决策树再寻找最佳节点。相反，它使用局部变量result记录到目前为止找到的最佳解。可以使用图 15-6 中的代码来测试max_val。

当运行small_test（使用图 15-3 中的值）时，它打印出一个结果，表明图 15-4 中的节点8是一个最佳解：

<c, 8, 2>
<b, 7, 3>
Total value of items taken = 15

图 15-5 使用决策树解决背包问题

函数build_many_items和big_test可用于在随机生成的项目集上测试max_val。尝试big_test(10, 40)。这并没有花费太长时间。现在尝试big_test(40, 100)。当你等得不耐烦时，停止计算，问问自己发生了什么。

让我们思考一下我们正在探索的树的大小。由于在树的每一层我们都在决定保留或不保留一个项目，因此树的最大深度为len(items)。在层级0时，我们只有一个节点；在层级1时最多有两个节点；在层级2时最多有四个节点；在层级3时最多有八个节点。在层级39时，最多有 2³⁹个节点。难怪运行起来需要很长时间！

让我们看看动态规划是否能提供帮助。

在图 15-4 和图 15-5 中都可以看到最优子结构。每个父节点结合其子节点达到的解决方案，以推导出该父节点根树的最优解决方案。这在图 15-5 中通过注释#Choose better branch后的代码得以体现。

图 15-6 测试基于决策树的实现

是否也存在重叠子问题？乍一看，答案似乎是“没有”。在树的每一层，我们都有不同的可用项目集可供考虑。这意味着，如果确实存在公共子问题，它们必须在树的同一层级。实际上，在树的每一层中，每个节点都有相同的项目集可供考虑。然而，通过查看图 15-4 中的标签，我们可以看到，每个层级的节点表示关于更高层树中所考虑项目的不同选择集。

思考每个节点正在解决的问题：在给定剩余可用重量的情况下，从剩下的可考虑项目中找到最佳项目。可用重量依赖于到目前为止所取项目的总重量，但与所取项目或所取项目的总价值无关。因此，例如，在图 15-4 中，节点 2 和 7 实际上在解决同一个问题：在可用重量为2的情况下，决定应该取[c,d]中的哪些元素。

图 15-7 中的代码利用了最优子结构和重叠子问题，提供了一种基于记忆化的动态规划解决方案，来解决 0/1 背包问题。

图 15-7 背包问题的动态规划解决方案

添加了一个额外的参数memo来跟踪已经解决的子问题的解。它是使用一个字典实现的，键由to_consider的长度和可用重量构成。表达式len(to_consider)是一种简洁的方式来表示仍需考虑的物品。这个方法有效，因为物品总是从列表to_consider的同一端（前端）移除。

现在，将对max_val的调用替换为对fast_max_val的调用，并尝试运行big_test(40, 100)。它几乎瞬间返回问题的最优解。

图 15-8 显示了当我们在具有不同物品数量和最大重量为 100 的情况下运行代码时的调用次数。增长量很难量化，但显然远低于指数级。⁹⁷但这怎么可能呢？因为我们知道 0/1 背包问题在物品数量上本质上是指数级的？我们是否找到了一种推翻宇宙基本法则的方法？不，但我们发现计算复杂性可能是一个微妙的概念。⁹⁸

图 15-8 动态编程解决方案的性能

fast_max_val的运行时间由生成的不同对<to_consider, avail>的数量决定。这是因为关于下一步该做什么的决定仅依赖于仍然可用的物品和已选物品的总重量。

to_consider的可能值数量受len(items)的限制。avail的可能值更难以界定。它的上限是背包可以容纳的项目重量总和的最大不同值。如果背包最多可以容纳n个物品（根据背包的容量和可用物品的重量），则avail最多可以取2^n 个不同值。从原则上讲，这可能是一个相当大的数字。然而，实际上通常并不是如此。即使背包的容量很大，如果物品的重量是从合理小的重量集合中选择的，许多物品的组合将具有相同的总重量，从而大大减少运行时间。

这个算法属于一个被称为伪多项式的复杂度类。对此概念的仔细解释超出了本书的范围。粗略来说，fast_max_val在表示avail可能值所需的位数上是指数级的。

要查看当avail的值从一个相当大空间中选择时会发生什么，请在图 15-6 中的big_test函数里将对max_val的调用更改为

val, taken = fast_max_val(items, 1000)

现在，当物品数量为1024时，找到解决方案需要1,028,403次对fast_max_val的调用。

为了观察在选取来自一个巨大空间的权重时会发生什么，我们可以从正实数中选择可能的权重，而不是从正整数中选择。为此，替换该行，

items.append(Item(str(i),
                  random.randint(1, max_val),
                  random.randint(1, max_weight)))

在build_many_items中，由于这一行

items.append(Item(str(i),
                  random.randint(1, max_val),
                  random.randint(1, max_weight)*random.random()))

每次调用时，random.random()返回一个介于0.0和1.0之间的随机浮点数，因此，在所有实用意义上，可能的权重数量是无限的。不要指望等待这个最后的测试完成。动态规划在一般意义上可能是一种奇迹技术，⁹⁹，但在礼仪意义上并不能创造奇迹。

15.3 动态规划与分治法

与分治算法类似，动态规划基于解决独立的子问题，然后将这些解决方案组合起来。然而，仍然存在一些重要的区别。

分治算法基于寻找比原始问题小得多的子问题。例如，归并排序通过在每一步将问题规模减半来工作。相对而言，动态规划涉及解决仅比原始问题稍小的问题。例如，计算第十九个斐波那契数并不是一个比计算第二十个斐波那契数小得多的问题。

另一个重要的区别是，分治算法的效率并不依赖于将算法结构化以便重复解决相同的问题。相比之下，动态规划仅在不同子问题的数量显著少于总子问题数量时才高效。

15.4 章节中介绍的术语

• 动态规划

• 最优子结构

• 重叠子问题

• 备忘录法

• 表格法

• 有根二叉树

• 根

• 叶子

• 决策树

• 伪多项式复杂度

第十六章：随机游走与数据可视化

本书讲述的是如何使用计算来解决问题。到目前为止，我们关注的问题是可以通过确定性程序解决的问题。程序是确定性的，如果它在相同输入下运行时，产生相同输出。这种计算非常有用，但显然不足以应对某些类型的问题。我们生活的世界的许多方面只能准确地建模为随机过程。¹⁰⁰ 如果一个过程的下一个状态可能依赖于某个随机因素，那么它就是随机的。随机过程的结果通常是不确定的。因此，我们很少能对随机过程会做什么作出明确的陈述。相反，我们会对它们可能做什么作出概率性陈述。本书的其余部分主要讨论构建帮助理解不确定情况的程序。这些程序中的许多将是模拟模型。

模拟模仿真实系统的活动。例如，图 10-11 中的代码模拟了一个人进行一系列按揭付款。可以将这段代码视为一种实验设备，称为模拟 模型，它提供有关所建模系统可能行为的有用信息。除此之外，模拟通常用于预测物理系统的未来状态（例如，50年后地球的温度），并代替那些进行成本过高、耗时或危险的真实世界实验（例如，税法变化的影响）。

重要的是要始终记住，模拟模型和所有模型一样，仅仅是对现实的近似。我们永远无法确定实际系统是否会按照模型预测的方式表现。事实上，我们通常相当自信地认为实际系统不会完全按模型预测的方式表现。例如，并不是每个借款人都会按时支付所有的按揭款。常言道：“所有模型都是错误的，但有些是有用的。” ¹⁰¹

16.1 随机游走

1827 年，苏格兰植物学家罗伯特·布朗观察到悬浮在水中的花粉颗粒似乎在随机漂浮。他对后来的布朗运动没有合理的解释，也没有试图从数学上进行建模。¹⁰² 1900 年，路易斯·巴歇利耶在其博士论文投机理论中首次提出了这一现象的清晰数学模型。然而，由于该论文涉及当时不光彩的金融市场理解问题，因此被主流学术界大多忽视。五年后，年轻的阿尔伯特·爱因斯坦将这种随机思维带入物理学界，提出了几乎与巴歇利耶的模型相同的数学模型，并描述了如何利用该模型确认原子的存在。¹⁰³ 不知为何，人们似乎认为理解物理学比赚钱更重要，因此世界开始关注这方面。那时的时代确实不同。

布朗运动是随机漫步的一个例子。随机漫步广泛用于模拟物理过程（例如扩散）、生物过程（例如 DNA 对异源双链 RNA 位移的动力学）和社会过程（例如股市的运动）。

在本章中，我们出于三个原因研究随机漫步：

• 随机漫步本身具有内在趣味性，且应用广泛。

• 它们为我们提供了一个良好的示例，说明如何使用抽象数据类型和继承来构建一般程序，尤其是模拟模型。

• 这些提供了引入 Python 更多特性和展示一些额外绘图技巧的机会。

16.2 醉汉的漫步

让我们来看一个实际上涉及走动的随机漫步。一个醉汉农民站在田地中央，每秒钟随机朝一个方向迈出一步。她（或他）在1000秒后与原点的预期距离是多少？如果她走了很多步，是否更可能离原点越来越远，还是更可能一次又一次地徘徊回到原点，最终不远离起点？让我们写一个模拟来找出答案。

在开始设计程序之前，尝试对程序所要模拟的情况发展一些直觉总是个好主意。让我们先使用笛卡尔坐标系勾勒出一个简单的模型。假设农民正站在一片田地中，草神秘地被修剪成类似于方格纸的形状。进一步假设，农民每走一步的长度为一，并且与 x 轴或 y 轴平行。

图 16-1 左侧的图片描绘了一个站在田野中间的农夫¹⁰⁴。微笑脸表示农夫在一步之后可能到达的所有地方。请注意，在一步之后，她始终离起始点恰好一单位。假设她在第一步时朝东从她的初始位置游荡。她在第二步之后可能离她的初始位置有多远？

图 16-1 一个不寻常的农夫

看着右侧图片中的微笑脸，我们看到她以0.25的概率会离原点0单位，以0.25的概率会离原点2单位，以0.5的概率会离原点单位。¹⁰⁵因此，平均而言，经过两步后她会比经过一步后离原点更远。第三步呢？如果第二步是到顶部或底部的微笑脸，第三步将使农夫在一半的情况下更靠近原点，而在另一半的情况下更远离。如果第二步是到左侧的微笑脸（原点），第三步将远离原点。如果第二步是到右侧的微笑脸，第三步在四分之一的时间内将更靠近原点，而在四分之三的时间内则会更远离原点。

看起来醉汉走的步数越多，预期离原点的距离就越大。我们可以继续这种可能性的详尽枚举，也许可以很好地直观理解这种距离是如何随着步数的增加而增长的。然而，这变得有些繁琐，因此写一个程序来为我们完成这项工作似乎是个更好的主意。

让我们通过思考一些可能在构建此模拟以及其他类型随机行走的模拟中有用的数据抽象来开始设计过程。像往常一样，我们应该尝试发明与我们试图建模的情况中出现的事物类型相对应的类型。三个显而易见的类型是Location、Field和Drunk。在查看提供这些类型的类时，考虑它们可能对我们能够构建的模拟模型暗示的内容是很有价值的。

让我们从Location开始，图 16-2。这是一个简单的类，但它体现了两个重要的决策。它告诉我们模拟将最多涉及两个维度。这与上面的图片是一致的。此外，由于提供给delta_x和delta_y的值可以是浮点数而不是整数，因此这个类没有关于醉汉可能移动的方向集合的内置假设。这是对非正式模型的一个概括，其中每一步的长度为一，并且与 x 轴或 y 轴平行。

类Field，图 16-2，也相当简单，但它同样体现了一些显著的决策。它只是维护了醉汉与位置的映射。它对位置没有任何限制，因此可以推测Field是无限大小的。它允许多个醉汉随机添加到Field中的任意位置。它没有对醉汉的移动模式做出任何说明，也不禁止多个醉汉占据同一位置或穿过其他醉汉占据的空间。

图 16-2 Location和Field类

图 16-3 中的Drunk和Usual_drunk类定义了醉汉在场地上漫游的方式。特别是，Usual_drunk中step_choices的值引入了每一步的长度为一且平行于 x 轴或 y 轴的限制。由于函数random.choice返回传入序列中随机选择的成员，因此每种步伐的可能性是相等的，并且不受之前步伐的影响。稍后我们将查看具有不同行为的Drunk子类。

图 16-3 定义醉汉的类

下一步是使用这些类构建一个模拟，以回答最初的问题。图 16-4 包含用于此模拟的三个函数。

图 16-4 醉汉的行走（带有错误）

函数walk模拟一次num_steps步的行走。函数sim_walks调用walk模拟num_trials次每次num_steps步的行走。函数drunk_test调用sim_walks以模拟不同长度的行走。

sim_walks的参数d_class是class类型，并在代码的第一行用于创建适当子类的Drunk。稍后，当drunk.take_step从Field.move_drunk调用时，将自动选择适当子类的方法。

函数drunk_test还有一个参数d_class，类型为class。它被使用了两次，一次在调用sim_walks时，另一次在第一个print语句中。在print语句中，内置的class属性__name__用于获取类名的字符串。

当我们执行drunk_test((10, 100, 1000, 10000), 100, Usual_drunk)时，它打印了

Usual_drunk walk of 10 steps: Mean = 8.634, Max = 21.6, Min = 1.4
Usual_drunk walk of 100 steps: Mean = 8.57, Max = 22.0, Min = 0.0
Usual_drunk walk of 1000 steps: Mean = 9.206, Max = 21.6, Min = 1.4
Usual_drunk walk of 10000 steps: Mean = 8.727, Max = 23.5, Min = 1.4

这令人惊讶，因为我们之前形成的直觉是，平均距离应该随着步数的增加而增加。这可能意味着我们的直觉是错误的，或者可能意味着我们的模拟有问题，或者两者都有。

此时的第一步是对我们已经认为知道答案的值运行模拟，确保模拟产生的结果与预期结果匹配。我们来试试零步数的步行（对于它，离原点的平均、最小和最大距离都应该是0）和一步的步行（对于它，离原点的平均、最小和最大距离都应该是1）。

当我们运行drunk_test((0,1), 100, Usual_drunk)时，得到了高度可疑的结果

Usual_drunk walk of 0 steps: Mean = 8.634, Max = 21.6, Min = 1.4
Usual_drunk walk of 1 steps: Mean = 8.57, Max = 22.0, Min = 0.0

零步数的步行的平均距离怎么可能超过8？我们的模拟中肯定有至少一个错误。经过一些调查，问题变得清楚。在sim_walks中，函数调用walk(f, Homer, num_trials)应该是walk(f, Homer, num_steps)。

这里的道德非常重要：在查看模拟结果时，始终保持一些怀疑态度。首先要问结果是否合理（即，看起来可信）。并始终对你对结果有强烈直觉的参数进行烟雾测试¹⁰⁶。

当纠正后的模拟在我们的两个简单案例中运行时，得到了完全预期的结果：

Usual_drunk walk of 0 steps: Mean = 0.0, Max = 0.0, Min = 0.0
Usual_drunk walk of 1 steps: Mean = 1.0, Max = 1.0, Min = 1.0

当在更长的步行中运行时，它打印出

Usual_drunk walk of 10 steps: Mean = 2.863, Max = 7.2, Min = 0.0
Usual_drunk walk of 100 steps: Mean = 8.296, Max = 21.6, Min = 1.4
Usual_drunk walk of 1000 steps: Mean = 27.297, Max = 66.3, Min = 4.2
Usual_drunk walk of 10000 steps: Mean = 89.241, Max = 226.5, Min = 10.0

正如预期的那样，离原点的平均距离随着步数的增加而增长。

现在让我们看看从原点的平均距离的图表，图 16-5。为了展示距离增长的速度，我们在图上绘制了一条显示步数平方根的线（并将步数增加到100,000）。

图 16-5 从起点到已走步数的距离

指尖练习： 编写代码生成图 16-5 中的图表。

这个图表是否提供了关于醉汉最终可能位置的任何信息？它确实告诉我们，平均而言，醉汉会在一个以原点为中心，半径等于预期距离的圆上。然而，它对于我们在任何特定步行结束时可能找到醉汉的具体位置几乎没有提供信息。我们将在下一节回到这个话题。

16.3 偏置随机步行

现在我们有了一个有效的模拟，可以开始修改它来研究其他类型的随机步行。例如，假设我们想考虑一个厌恶寒冷的北半球醉汉农夫的行为，即使在醉酒状态下，他在朝南方向移动时速度也会快一倍。或者可能是一个光向性醉汉，总是朝着太阳移动（早上向东，下午向西）。这些都是偏置随机步行的例子。步行仍然是随机的，但结果中存在偏差。

图 16-6 定义了两个额外的Drunk子类。在每种情况下，专业化涉及为step_choices选择合适的值。函数sim_all迭代一系列Drunk子类，以生成有关每种类型行为的信息。

图 16-6 Drunk基类的子类

当我们运行时

`sim_all((Usual_drunk, Cold_drunk, EW_drunk), (100, 1000), 10)`

它打印了

Usual_drunk walk of 100 steps: Mean = 9.64, Max = 17.2, Min = 4.2
Usual_drunk walk of 1000 steps: Mean = 22.37, Max = 45.5, Min = 4.5
Cold_drunk walk of 100 steps: Mean = 27.96, Max = 51.2, Min = 4.1
Cold_drunk walk of 1000 steps: Mean = 259.49, Max = 320.7, Min = 215.1
EW_drunk walk of 100 steps: Mean = 7.8, Max = 16.0, Min = 0.0
EW_drunk walk of 1000 steps: Mean = 20.2, Max = 48.0, Min = 4.0

看起来我们的热寻求醉汉比其他两种醉汉更快地远离原点。然而，消化这段输出中的所有信息并不容易。是时候远离文本输出，开始使用图表了。

由于我们在同一个图表上展示不同类型的醉汉，我们将为每种醉汉关联一种独特的样式，以便于区分。样式将具有三个方面：

• 线条和标记的颜色

• 标记的形状

• 线条的类型，例如，实线或虚线。

类style_iterator，图 16-7 通过传递给style_iterator.__init__的参数旋转一系列样式。

图 16-7 迭代样式

图 16-8 中的代码在结构上与图 16-4 中的代码相似。

图 16-8 绘制不同醉汉的步态

sim_drunk和sim_all_plot中的print语句对仿真的结果没有贡献。它们的存在是因为这个仿真可能需要较长时间完成，偶尔打印一条指示进展的消息可以让用户感到安心，避免他们怀疑程序是否真的在运行。

图 16-9 中的图是通过执行生成的。

sim_all_plot((Usual_drunk, Cold_drunk, EW_drunk),
             (10, 100, 1000, 10000, 100000), 100)

图 16-9 不同类型醉汉的平均距离

普通醉汉和光向性醉汉（EW_drunk）似乎以大约相同的速度远离原点，但热寻求醉汉（Cold_drunk）似乎远离的速度要快几个数量级。这很有趣，因为平均而言，他只比其他人快25%（他每四步平均走五步）。

让我们构建一个不同的图表，以帮助我们更深入地了解这三个类的行为。代码在图 16-10 中绘制了单步数的最终位置分布，而不是随着步数增加绘制距离随时间的变化。

图 16-10 绘制最终位置

plot_locs的第一个操作是创建一个style_iterator实例，该实例有三种样式的标记。然后它使用plt.plot在每次试验结束时对应的位置放置标记。对plt.plot的调用设置了标记的颜色和形状，使用了style_iterator返回的值。

调用plot_locs((Usual_drunk, Cold_drunk, EW_drunk), 100, 200)生成了图 16-11 中的图。

图 16-11 醉汉停下的地方

首先要说的是，我们的醉汉似乎表现如预期。EW_drunk最终位于 x 轴上，Cold_drunk似乎向南推进，而Usual_drunk则显得漫无目的。

但为什么圆形标记似乎比三角形或+标记少得多呢？因为许多EW_drunk的行走最终到达了同一个地方。这并不令人惊讶，考虑到EW_drunk可能的终点数量只有200。而且，圆形标记在 x 轴上似乎分布相对均匀。

至少对我们来说，Cold_drunk为什么在平均情况下能比其他类型的醉汉走得更远，这一点仍然不是显而易见的。也许是时候关注一条单一行走的路径，而不是许多行走的终点。 图 16-12 中的代码生成了图 16-13 中的图。

图 16-12 行走轨迹追踪

图 16-13 行走轨迹

由于行走长度为200步，而EW_drunk的行走访问了不到30个不同的位置，因此显然他花费了很多时间在回溯自己的步骤上。对于Usual_drunk也是如此。相比之下，虽然Cold_drunk并不是直奔佛罗里达，但他相对花费更少的时间去访问他已经去过的地方。

这些模拟本身并没有什么有趣之处。（在第十八章中，我们将关注更内在有趣的模拟。）但有几点值得注意：

• 起初，我们将模拟代码分为四个独立的部分。其中三个是类（Location、Field和Drunk），对应于问题非正式描述中出现的抽象数据类型。第四部分是一组函数，这些函数使用这些类执行简单的模拟。

• 我们随后将Drunk详细划分为一个类的层次结构，以便观察不同类型的偏向随机行走。Location和Field的代码保持不变，但模拟代码被更改为迭代Drunk的不同子类。在此过程中，我们利用了类本身是一个对象，因此可以作为参数传递的事实。

• 最后，我们对模拟进行了一系列渐进的更改，而没有涉及表示抽象类型的类的任何变化。这些更改主要涉及引入旨在提供对不同游走的洞察的图表。这是模拟开发的典型方式。首先让基本模拟工作，然后再开始添加功能。

16.4 危险场

你是否玩过在美国被称为滑梯与梯子、在英国被称为蛇与梯子的桌面游戏？这个儿童游戏起源于印度（可能是在公元前二世纪），当时称为Moksha-patamu。落在代表美德（例如，慷慨）的方块上，会让玩家爬上梯子，进入更高的生活层次。落在代表邪恶（例如，欲望）的方块上，则会将玩家送回更低的生活层次。

我们可以通过创建一个带有虫洞的Field轻松地为我们的随机游走添加这种特性，¹⁰⁷，如图 16-14 所示，并将函数trace_walk中的第二行代码替换为

`f = Odd_field(1000, 100, 200)`

在Odd_field中，走入虫洞位置的醉汉会被传送到虫洞另一端的位置。

图 16-14 具有奇怪属性的场

当我们运行trace_walk((Usual_drunk, Cold_drunk, EW_drunk), 500)时，得到了图 16-15 中相当奇怪的图。

图 16-15 奇怪的游走

显然，更改场的属性产生了显著效果。然而，这并不是本示例的重点。主要观点是：

• 由于我们代码的结构方式，适应被建模情况的重大变化变得容易。正如我们可以添加不同类型的醉汉而不触动Field一样，我们也可以添加一种新的Field而不触动Drunk或其任何子类。（如果我们足够有远见，将场作为trace_walk的参数，我们就不必更改trace_walk了。）

• 尽管可以从分析上推导出简单随机游走甚至有偏随机游走的预期行为的不同信息，但一旦引入虫洞，这就变得很具挑战性。然而，改变模拟以适应新情况却极为简单。模拟模型通常相对于分析模型享有这一优势。

16.5 在章节中引入的术语确定性程序

• 随机过程

• 模拟模型

• 随机游走

• 烟雾测试

• 有偏随机游走

• 对数尺度

第十七章：随机程序、概率和分布

纽顿力学给人一种安慰。你在杠杆的一端施加压力，另一端就会上升。你把球抛向空中；它沿着抛物线轨迹移动，最终落下。简而言之，一切都是有原因的。物理世界是一个完全可预测的地方——所有物理系统的未来状态都可以从对其当前状态的知识中推导出来。

几个世纪以来，这一直是主流科学智慧；然后量子力学和哥本哈根公理出现了。以波尔和海森堡为首的公理支持者们争辩说，从最根本的层面来看，物理世界的行为是无法预测的。可以做出“x 很可能发生”的概率性陈述，但不能做出“x 必然发生”的陈述。其他著名物理学家，尤其是爱因斯坦和薛定谔， vehemently disagreed.

这场辩论搅动了物理学、哲学甚至宗教的世界。辩论的核心是因果非确定性的有效性，即并非每个事件都是由之前的事件引起的信念。爱因斯坦和薛定谔发现这种观点在哲学上不可接受，正如爱因斯坦经常重复的那句话，“上帝不会掷骰子。”他们所能接受的是预测非确定性，即我们的测量能力有限使得无法对未来状态做出准确预测。这一区别被爱因斯坦很好地总结，他说：“当代理论的本质统计特征完全归因于这一理论以不完整的物理系统描述为基础。”

因果非确定性的问题仍然没有解决。然而，我们无法预测事件的原因是因为它们真正不可预测，还是因为我们仅仅没有足够的信息去预测它们，这在实践中并无重要意义。

当波尔/爱因斯坦的辩论涉及如何理解物理世界的最低层次时，相同的问题也出现在宏观层面。或许马赛的结果、轮盘赌的旋转以及股票市场的投资是因果确定的。然而，有充足的证据表明，将它们视为可预测的确定性是危险的。

17.1 随机程序

如果一个程序是确定性的，那么每当它在相同输入上运行时，它产生相同的输出。注意，这与说输出完全由问题的规格定义并不相同。考虑，例如，square_root的规格：

`def square_root(x, epsilon):     """Assumes x and epsilon are of type float; x >= 0 and epsilon > 0        Returns float y such that x-epsilon <= y*y <= x+epsilon"""`

这个规格为函数调用square_root(2, 0.001)允许了许多可能的返回值。然而，我们在第三章看到的连续逼近算法将始终返回相同的值。该规格并不要求实现是确定性的，但它确实允许确定性实现。

不是所有有趣的规格都可以通过确定性实现来满足。例如，考虑实现一个骰子游戏的程序，比如西洋双陆棋或掷骰子。程序中会有一个函数来模拟一个公平的六面骰子的掷骰子。¹⁰⁹假设它的规格是这样的。

`def roll_die():     """Returns an int between 1 and 6"""` 

这将是个问题，因为它允许实现每次调用时返回相同的数字，这会使游戏变得相当无聊。更好的方法是指定roll_die“返回一个在 1 和 6 之间随机选择的整数”，这样就要求实现是随机的。

大多数编程语言，包括 Python，都提供了简单的方法来编写随机程序，即利用随机性的程序。图 17-1 中的小程序是一个模拟模型。我们不是让某个人多次掷骰子，而是写了一个程序来模拟这个活动。代码使用了从导入的 Python 标准库模块random中找到的几个有用函数之一。

import random
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate

所有这些都已执行。

函数random.choice以一个非空序列作为参数，并返回该序列中随机选择的成员。几乎所有random中的函数都是使用函数random.random构建的，该函数生成一个在0.0和1.0之间的随机浮点数。¹¹⁰

图 17-1 掷骰子

现在，想象一下运行roll_n(10)。你会对打印出1111111111还是5442462412感到更惊讶？换句话说，这两个序列哪个更随机？这是一个技巧性问题。这两个序列发生的可能性是一样的，因为每次掷骰子的值与之前的掷骰子值是独立的。在随机过程中的两个事件是独立的，如果一个事件的结果对另一个事件的结果没有影响。

如果我们通过考虑一个两面的骰子（也称为硬币），其值为 0 和 1 来简化情况，这样更容易理解。这样我们可以将roll_n的输出视为一个二进制数字。当我们使用二进制骰子时，roll_n可能返回的序列有2^n 种可能性。这些序列的可能性是一样的，因此每个序列发生的概率为(1/2)^n。

让我们回到六面骰子。长度为10的不同序列有多少种？6¹⁰。所以，连续掷到 10 个1的概率是1/6¹⁰。少于六千万分之一。相当低，但不低于任何其他序列的概率，例如5442462412。

17.2 计算简单概率

一般来说，当我们谈论某个结果具有某种属性的概率（例如，都是1）时，我们是在询问所有可能结果中有多少比例具有该属性。这就是为什么概率的范围是从 0 到 1。假设我们想知道掷骰子时得到的序列不是全部1的概率。它只是1 – (1/6¹⁰)，因为某件事发生的概率和同一事件不发生的概率之和必须为1。

假设我们想知道在 10 次掷骰子中没有掷到一个1的概率。回答这个问题的一种方法是将其转化为6¹⁰个可能序列中有多少个不包含1。这可以通过以下方式计算：

1. 1. 在任何一次掷骰子中，不掷到1的概率是5/6。

2. 2. 在第一次或第二次掷骰子时，不掷到1的概率是(5/6)*(5/6)，或(5/6)²。

3. 3. 所以，不连续掷到1 10 次的概率是(5/6)¹⁰，略高于0.16。

第 2 步是对独立概率的乘法法则的应用。例如，考虑两个独立事件A和B。如果A发生的概率为1/3，而B发生的概率为1/4，那么A和B同时发生的概率为1/4 of 1/3，即(1/3)/4或(1/3)*(1/4)。

至于掷到至少一个1的概率呢？它就是 1 减去不掷到至少一个1的概率，即1 - (5/6)¹⁰。注意，不能简单地说任何一次掷骰子的1的概率是1/6，因此至少掷到一个1的概率是10*(1/6)，即10/6。这显然是不正确的，因为概率不能大于1。

那么在 10 次掷骰子中掷到恰好两个1的概率是多少呢？这相当于问在6¹⁰个整数中，有多少个整数在其基数6表示中恰好有两个1。我们可以轻松编写程序生成所有这些序列，并计算包含恰好一个1的数量。用解析方法推导这个概率有点棘手，我们将其留到第 17.4.4 节讨论。

17.3 推论统计

正如你刚才看到的，我们可以使用系统化的过程基于已知的一个或多个简单事件的概率，推导出一些复杂事件的精确概率。例如，我们可以轻松计算投掷硬币得到 10 次连续正面的概率，前提是我们假设每次投掷是独立的，并且知道每次投掷出现正面的概率。然而，假设我们实际上并不知道相关的简单事件的概率。举例来说，假设我们不知道这枚硬币是否公平（即，正反面出现的概率相等）。

一切尚未失去。如果我们对硬币的行为有一些数据，我们可以将这些数据与我们的概率知识结合，以推导出真实概率的估计。我们可以使用推断统计来估计一次投掷出现正面的概率，然后用传统概率来计算这种行为的硬币连续出现 10 次正面的概率。

简而言之（因为这不是一本关于统计的书），推断统计的指导原则是，随机样本往往表现出与其来源总体相同的特性。

假设哈维·登特（也称为“双面人”）投掷了一枚硬币，结果是正面。你不会由此推断下一次投掷也会是正面。假设他投掷了两次，结果都是正面。你可能会推测，对于一枚公平的硬币，这种情况发生的概率是0.25，因此仍然没有理由假设下一次投掷会是正面。然而，假设100次投掷中都有100次是正面。事件的概率(1/2)¹⁰⁰（假设硬币是公平的）非常小，因此你可能会觉得推断这枚硬币的两面都是正面是安全的。

你对硬币是否公平的信念基于这样的直觉：单次100次投掷的行为与所有100次投掷样本的总体行为是相似的。当所有100次投掷都是正面时，这种信念似乎是合适的。假设有52次投掷是正面，48次是反面。你是否会对预测下一次100次投掷的正反面比例保持信心？更进一步，你会对预测下一次100次投掷中正面会多于反面感到有把握吗？花几分钟思考一下这个问题，然后尝试进行实验。或者，如果你没有现成的硬币，可以使用图 17-2 中的代码来模拟投掷。

图 17-2 中的函数flip模拟了一枚公平硬币投掷num_flips次，并返回其中出现正面的比例。对于每次投掷，调用random.choice(('H', ‘T'))随机返回'H'或'T'。

图 17-2 投掷硬币

尝试执行函数 flip_sim(10, 1) 几次。我们在前两次尝试 print('Mean =', flip_sim(10, 1)) 时看到的结果如下：

Mean = 0.2
Mean = 0.6

从单次 10 次抛掷的试验中假设太多（除了硬币有正反两面）似乎不合适。这就是为什么我们通常将模拟结构化为包含多个试验并比较结果。让我们尝试 flip_sim(10, 100) 几次：

Mean = 0.5029999999999999
Mean = 0.496

你对这些结果感觉更好了吗？当我们尝试 flip_sim(100, 100000) 时，我们得到了

Mean = 0.5005000000000038
Mean = 0.5003139999999954

这看起来真的很好（尤其是我们知道答案应该是 0.5——但这算作弊）。现在似乎我们可以安全地得出关于下一次抛掷的结论，即正反面大致同样可能。但我们为什么认为可以得出这样的结论呢？

我们所依赖的是 大数法则（也称为 伯努利定理¹¹¹）。该法则声明，在相同实际概率 p 的重复独立测试中（在这种情况下为抛掷），该结果的出现频率与 p 的差异概率随着试验次数趋向于无穷大而收敛为零。

值得注意的是，大数法则并不意味着，正如太多人所想，若出现与预期行为的偏差，这些偏差可能会被未来的相反偏差“抵消”。这种对大数法则的误用被称为 赌徒谬误。¹¹²

人们常常将赌徒谬误与回归均值混淆。回归均值¹¹³表示，在经历极端随机事件后，下一个随机事件可能会更不极端。如果你抛一枚公正的硬币六次并得到六个正面，回归均值暗示下一个六次抛掷的结果更接近三个正面的期望值。它并不意味着，正如赌徒谬误所暗示的，下一次抛掷的正面数量可能少于反面数量。

大多数事务的成功需要技能与运气的结合。技能成分决定了均值，而运气成分则解释了变异性。运气的随机性导致回归均值。

图 17-3 中的代码生成了一个图，图 17-4，说明了回归均值。函数 regress_to_mean 首先生成 num_trials 次试验，每次 num_flips 次抛硬币。然后，它识别出所有结果中正面比例低于 1/3 或高于 2/3 的试验，并将这些极端值绘制为圆圈。接着，对于每个这些点，它在与圆圈相同的列中绘制后续试验的值作为三角形。

0.5处的水平线，即预期均值，是使用axhline函数创建的。plt.xlim函数控制 x 轴的范围。函数调用plt.xlim(xmin, xmax)设置当前图形 x 轴的最小值和最大值。函数调用plt.xlim()返回一个元组，包含当前图形 x 轴的最小值和最大值。plt.ylim函数的工作方式相同。

图 17-3 回归到均值

注意到，在极端结果后，通常会跟随一个更接近均值的试验，但这并不总是发生——如框中的一对所示。手指练习： 安德烈在打高尔夫时每洞平均5杆。一天，她完成前九洞用了40杆。她的搭档推测她可能会回归均值，接下来用50杆完成后九洞。你同意她的搭档吗？ 图 17-4 均值回归的插图 图 17-5 包含一个函数flip_plot，生成两个图， 图 17-6，旨在展示大数法则的应用。第一个图显示头和尾的数量差的绝对值如何随翻转次数变化。第二个图比较头与尾的比率和翻转次数。由于 x 轴上的数字较大，我们使用plt.xticks(rotation = ‘vertical')来旋转它们。 图 17-5 掷硬币结果的绘图 random.seed(0)在图底部确保random.random使用的伪随机数生成器每次执行时生成相同的伪随机数序列。这在调试时很方便。函数random.seed可以用任何数字调用。如果没有参数调用，则随机选择种子。 图 17-6 大数法则的应用 左边的图似乎表明，头和尾的绝对差在开始时波动，随后迅速下跌，然后迅速上升。然而，我们需要记住，x = 300,000右侧只有两个数据点。plt.plot用线连接这些点可能会误导我们看到趋势，而实际上只有孤立点。这不是不常见的现象，所以在对图的含义做出任何结论之前，总是要询问图实际包含多少个点。 在右侧的图中，很难看出任何内容，主要是平坦的线。这也是误导。尽管有 16 个数据点，但大多数都挤在图左侧的小区域中，因此细节难以看清。这是因为绘制的点的 x 值是2⁴, 2⁵, 2⁶, …, 2²⁰，因此 x 轴的值范围从16到超过一百万，除非另行指示，plot会根据这些点与原点的相对距离来放置它们。这被称为**线性缩放**。由于大多数点的 x 值相对于2²⁰较小，它们看起来相对接近原点。 幸运的是，这些可视化问题在 Python 中很容易解决。正如我们在第十三章和本章早些时候看到的，我们可以轻松指示程序绘制不连接的点，例如，通过编写plt.plot(xAxis, diffs, 'ko')。 图 17-7 中的两个图都在 x 轴上使用对数尺度。由于flip_plot生成的 x 值为2^(minExp), 2^(minExp+1), …, 2^(maxExp)，使用对数 x 轴使得点沿 x 轴均匀分布——提供点之间的最大分离。 图 17-7 左侧的图在 y 轴上也使用对数尺度。该图的 y 值范围从接近0到约550。如果 y 轴是线性缩放，则很难看到图左端 y 值的相对小差异。另一方面，右侧的图 y 值相对较紧凑，因此我们使用线性 y 轴。![c17-fig-0007.jpg](https://gitee.com/OpenDocCN/geekdoc-ds-zh/tree/master/docs/intr-comp-prog-py/img/c17-fig-0007.jpg) 图 17-7 大数法则的应用 **手指练习：** 修改图 17-5 中的代码，使其生成如图 17-7 所示的图。 图 17-7 中的图比早期图更容易解释。右侧的图强烈暗示，随着翻转次数增大，头与尾的比率收敛于1.0。左侧图的意义不太清晰。似乎绝对差随着翻转次数增长，但这并不完全令人信服。 通过抽样实现完美的准确性永远是不可能的，除非抽样整个种群。无论我们检查多少样本，只有在检查种群的每个元素后，才能确定样本集是否典型（而且由于我们通常处理的是无限种群，例如所有可能的掷硬币序列，这往往是不可能的）。当然，这并不是说估计不能是精确正确的。我们可能掷硬币两次，得到一头一尾，得出每个结果的真实概率是0.5。我们会得出正确的结论，但我们的推理是错误的。 我们需要查看多少个样本，才能对我们计算的答案有合理的信心？这取决于**方差**在潜在分布中的情况。大致而言，方差是可能不同结果的扩散程度的度量。更正式地，值集合*X*的方差定义为：![c17-fig-5002.jpg](https://gitee.com/OpenDocCN/geekdoc-ds-zh/tree/master/docs/intr-comp-prog-py/img/c17-fig-5002.jpg) 其中|*X*|是集合的大小，*μ*（mu）是其均值。非正式地，方差描述有多少比例的值接近均值。如果许多值相对接近均值，方差相对较小。如果许多值相对远离均值，方差相对较大。如果所有值都相同，方差为零。 一组值的**标准差**是方差的平方根。虽然它包含与方差完全相同的信息，但标准差更易于解释，因为它与原始数据的单位相同。例如，“一个种群的平均身高是70英寸，标准差是4英寸”的说法比“一个种群的平均身高是70英寸，方差是16平方英寸”的说法更易于理解。图 17-8 包含方差和标准差的实现。¹¹⁵![c17-fig-0008.jpg](https://gitee.com/OpenDocCN/geekdoc-ds-zh/tree/master/docs/intr-comp-prog-py/img/c17-fig-0008.jpg) 图 17-8 方差和标准差 我们可以利用标准差的概念来思考我们查看的样本数量与我们在计算的答案中应有的信心程度之间的关系。图 17-10 包含一个修改版的flip_plot。它使用辅助函数run_trial运行每个掷硬币次数的多个试验，然后绘制头/尾比率的均值和标准差。辅助函数make_plot，图 17-9，包含生成图的代码。![c17-fig-0009.jpg](https://gitee.com/OpenDocCN/geekdoc-ds-zh/tree/master/docs/intr-comp-prog-py/img/c17-fig-0009.jpg) 图 17-9 掷硬币模拟的辅助函数 ![c17-fig-0010.jpg](https://gitee.com/OpenDocCN/geekdoc-ds-zh/tree/master/docs/intr-comp-prog-py/img/c17-fig-0010.jpg) 图 17-10 掷硬币模拟 让我们尝试flip_plot1(4, 20, 20)。它生成图 17-11 中的图。![c17-fig-0011.jpg](https://gitee.com/OpenDocCN/geekdoc-ds-zh/tree/master/docs/intr-comp-prog-py/img/c17-fig-0011.jpg) 图 17-11 头尾比率的收敛 这令人鼓舞。平均头/尾比率正在收敛于1，标准差的对数与每次试验的翻转次数的对数线性下降。当我们进行大约10⁶次掷硬币时，标准差（约10^(-3)）大约比均值（约1）小三个数量级，因此我们可以相当有信心预期头/尾比率非常接近1.0`。随着我们掷更多硬币，我们不仅得到更精确的答案，更重要的是，我们还有更多理由相信它

第十八章：蒙特卡罗模拟

在第十六章和第十七章中，我们探讨了在计算中使用随机性的不同方式。我们呈现的许多示例属于被称为蒙特卡罗模拟的计算类别。蒙特卡罗模拟是一种通过多次运行相同模拟并平均结果来近似事件概率的技术。

斯坦尼斯瓦夫·乌拉姆和尼古拉斯·梅特罗波利斯于 1949 年创造了“蒙特卡罗模拟”这一术语，以向摩纳哥公国赌场中的机会游戏致敬。乌拉姆以与爱德华·泰勒共同设计氢弹而闻名，他这样描述该方法的发明：

我对[蒙特卡罗方法]的首次思考和尝试，是在 1946 年我在康复期间玩接龙时想到的一个问题。这个问题是，一个用 52 张牌摆成的坎菲尔德接龙成功的机会是多少？在花了大量时间试图通过纯组合计算来估算它们之后，我想知道是否有比“抽象思维”更实际的方法，比如说摆一百次并简单观察和计算成功的次数。随着新一代快速计算机的出现，这已经是可以想象的了，**¹²⁴ 我立即想到了中子扩散问题和其他数学物理问题，以及更一般地如何将由某些微分方程描述的过程转变为可解释为随机操作序列的等效形式。后来…… [在 1946 年，我] 向约翰·冯·诺依曼描述了这个想法，我们开始规划实际计算。**¹²⁵

该技术在曼哈顿计划期间被用来预测核裂变反应会发生什么，但直到 1950 年代计算机变得更加普及和强大时，它才真正起飞。

乌拉姆并不是第一个想到使用概率工具来理解机会游戏的数学家。概率的历史与赌博的历史密切相关。正是这种不确定性使得赌博成为可能。而赌博的存在促使了许多必要的数学的发展，以便更好地推理不确定性。卡尔达诺、帕斯卡尔、费马、伯努利、德摩根和拉普拉斯对概率理论基础的贡献，都是出于希望更好地理解（并可能从中获利）机会游戏的愿望。

18.1 帕斯卡尔问题

早期概率理论的大部分研究围绕着使用骰子的游戏进行。传闻，帕斯卡对后来被称为概率理论的领域的兴趣开始于一个朋友问他，在 24 次掷骰中掷出双6是否会有利可图。这在十七世纪中期被视为一个棘手的问题。帕斯卡和费尔马这两位聪明人就如何解决这个问题交换了多封信件，但现在看来这是一个容易回答的问题：

• 在第一次掷骰中，每个骰子掷出6的概率为1/6，因此两个骰子同时掷出6的概率为1/36。

• 因此，第一次掷骰不掷出双6的概率为1 ‑ 1/36 = 35/36。

• 因此，连续 24 次掷骰不掷出6的概率为(35/36)²⁴，约为0.51，因此掷出双6的概率为1 - (35/36)²⁴，大约为0.49。从长远来看，押注在 24 次掷骰中掷出双6并不划算。

为了安全起见，我们写一个小程序，图 18-1，来模拟帕斯卡朋友的游戏，并确认我们得到的答案与帕斯卡相同。（本章中的所有代码都假设

import random
import numpy as np

在代码出现的文件开始时出现。当第一次运行时，调用check_pascal(1000000)打印了结果。

`Probability of winning = 0.490761`

这确实非常接近1 - (35/36)²⁴；在 Python shell 中输入1-(35.0/36.0)**24产生0.49140387613090342。

图 18-1 检查帕斯卡的分析

18.2 过线或不过线？

关于机会游戏的问题并非都那么容易回答。在掷骰子游戏中，掷骰者（投掷骰子的人）在“过线”或“不过线”投注之间选择。

• 过线：如果第一次掷骰结果是“自然数”（7或11），则掷骰者获胜；如果结果是“掷骰”（2、3或12），则掷骰者失败。如果掷出其他数字，该数字成为“点数”，掷骰者继续掷骰。如果掷骰者在掷出7之前掷出点数，掷骰者获胜；否则掷骰者失败。

• 不过线：如果第一次掷骰结果是7或11，则掷骰者失败；如果结果是2或3，则掷骰者获胜；如果结果是12，则平局（在赌博术语中称为“平推”）。如果掷出其他数字，该数字成为点数，掷骰者继续掷骰。如果掷骰者在掷出点数之前掷出7，掷骰者获胜；否则掷骰者失败。

这两者中哪个更划算？哪一个算是一个好赌注吗？可以通过分析得出这些问题的答案，但对我们来说，编写一个模拟掷骰游戏的程序并看看发生了什么似乎更简单。图 18-2 包含了这种模拟的核心。

图 18-2 Craps_game类

Craps_game类的实例变量记录自游戏开始以来通过线和不通过线的表现。观察方法pass_results和dp_results返回这些值。方法play_hand模拟一局游戏。一局的开始是当投掷者“出场”，这是在确定点数之前的投掷术语。一局在投掷者赢得或失去初始赌注时结束。play_hand中的大部分代码仅是上述规则的算法描述。请注意，else子句中有一个循环，对应于确定点数后发生的情况。当掷出七点或点数时，通过break语句退出循环。

图 18-3 包含一个使用Craps_game类模拟一系列掷骰子游戏的函数。

图 18-3 模拟掷骰子游戏

craps_sim的结构是许多模拟程序的典型：

1. 1. 它运行多个游戏（可以把每个游戏看作是我们之前模拟中的一次试验）并累积结果。每个游戏包括多局，因此有一个嵌套循环。

2. 2. 然后它生成并存储每个游戏的统计数据。

3. 3. 最后，它生成并输出汇总统计数据。在这种情况下，它打印每种投注线的预计投资回报率（ROI）和该 ROI 的标准差。

投资回报率由以下方程定义¹²⁷

由于通过线和不通过线支付的是等额奖金（如果你投注$1并赢得，你获得$1），因此投资回报率为

例如，如果你进行了100次通过线的投注并赢得了一半，那么你的 ROI 将是

如果你在不通过线上投注 100 次，并赢得 25 次和 5 次平局，那么 ROI 将是

让我们运行掷骰子游戏模拟，看看尝试craps_sim(20, 10)时会发生什么：¹²⁸

Pass: Mean ROI = -7.0% Std. Dev. = 23.6854%
Don't pass: Mean ROI = 4.0% Std Dev = 23.5372%

看起来避免通过线是个好主意——预计投资回报为7%的损失。但不通过线似乎是一个不错的赌注。或者说并非如此？

从标准差来看，不通过线似乎并不是如此安全的赌注。回想一下，在假设分布是正态的情况下，95%的置信区间由均值两侧的1.96个标准差包围。对于不通过线，95%的置信区间是[4.0–1.9623.5372, 4.0+1.9623.5372]——大约为[-43%, +51%]。这当然并不表明投注不通过线是万无一失的。

是时候将大数法则付诸实践；craps_sim(1000000, 10) 输出

Pass: Mean ROI = -1.4204% Std. Dev. = 0.0614%
Don't pass: Mean ROI = -1.3571% Std Dev = 0.0593%

我们现在可以相对安全地假设这两个选项都不是好的投注。¹²⁹ 看起来不通过线可能稍微好一些，但我们可能不应该依赖于此。如果通过线和不通过线的95%置信区间没有重叠，假设两个均值之间的差异是统计显著的将是安全的。¹³⁰ 然而，它们确实重叠，所以不能安全得出结论。

假设我们不是增加每局的手数，而是增加游戏的局数，例如，通过调用craps_sim(20, 1000000)：

Pass: Mean ROI = -1.4133% Std. Dev. = 22.3571%
Don't pass: Mean ROI = -1.3649% Std Dev = 22.0446%

标准差很高——这表明一局20手的结果高度不确定。

模拟的一个好处是它们使得进行“如果”实验变得简单。例如，如果一名玩家能够悄悄地引入一对有利于5而非2的作弊骰子（5和2在骰子的对面）呢？要测试这个，只需将roll_die的实现替换为类似的内容。

def roll_die():
    return random.choice([1,1,2,3,3,4,4,5,5,5,6,6])

骰子上的这个相对较小的变化在赔率上产生了显著差异。运行craps_sim(1000000, 10)得出的结果是

Pass: Mean ROI = 6.6867% Std. Dev. = 0.0993%
Don't pass: Mean ROI = -9.469% Std Dev = 0.1024%

难怪赌场费尽心思确保玩家不在游戏中引入自己的骰子！

手指练习：如果在掷出 6 之前掷出了 7，"大 6"的投注支付是平赔。假设每小时有 30 个$5 的投注，编写一个蒙特卡罗模拟来估算玩“大 6”投注的每小时成本和该成本的标准差。

18.3 使用表查找来提高性能

你可能不想在家尝试运行craps_sim(100000000, 10)。在大多数计算机上完成这个任务需要很长时间。这引出了一个问题：是否有简单的方法来加快模拟速度。

函数craps_sim的实现复杂度大约是θ(``play_hand``)hands_per_gamenum_games。play_hand的运行时间取决于循环执行的次数。从理论上讲，循环可以被执行无限次，因为掷出7或点数所需的时间没有上限。但实际上，我们有充分理由相信它总会终止。

不过请注意，play_hand的调用结果并不依赖于循环执行的次数，而只取决于达到的退出条件。对于每个可能的点，我们可以轻松计算在掷出7之前掷出该点的概率。例如，使用一对骰子我们可以用三种方式掷出4：<1, 3>, <3, 1>, 和 <2, 2>。我们可以用六种方式掷出7：<1, 6>, <6, 1>, <2, 5>, <5, 2>, <3, 4> 和 <4, 3>。因此，通过掷出7退出循环的可能性是通过掷出4的两倍。

图 18-4 包含了一个利用这种思维的play_hand实现。我们首先计算在掷出7之前，针对每个可能的点值形成该点的概率，并将这些值存储在字典中。例如，假设点值为8。投掷者会继续掷骰，直到掷出该点或掷出“掷坏”。掷出8有五种方式(<6,2>, <2,6>, <5,3>, <3,5>, <4,4>)，而掷出7有六种方式。因此，字典键8的值为表达式5/11的值。拥有这个表允许我们将包含不受限制掷骰次数的内部循环替换为对一次random.random调用的测试。这个版本的play_hand的渐进复杂度为O(1)。

用查表替代计算的思想具有广泛的适用性，通常在速度成为问题时使用。查表是以时间换空间这一一般思想的一个例子。正如我们在第十五章中看到的，这也是动态规划背后的关键思想。在我们对哈希分析中看到过这一技术的另一个例子：表越大，碰撞越少，平均查找速度越快。在这种情况下，表很小，所以空间成本微不足道。

图 18-4 使用查表来提高性能

18.4 寻找π

很容易看出，蒙特卡罗模拟在解决非确定性问题时是有用的。有趣的是，蒙特卡罗模拟（以及随机算法一般）也可以用来解决那些本质上并不是随机的问题，即没有结果不确定性的问题。

考虑π。数千年来，人们一直知道有一个常数（自 18 世纪以来称为π），使得圆的周长等于π * 直径，而圆的面积等于π * 半径²。他们不知道的是这个常数的值。

最早的估算之一，4(8/9)² = 3.16，可以在公元前 1650 年的埃及林德纸草书中找到。千年之后，旧约（列王纪上 7.23）在给出所罗门王某个建筑项目的规格时暗示了一个不同的π*值，

他造了一个铸造的海，从一边到另一边十肘：它周围是圆的，高五肘：周围有三十肘的线环绕着它。

解决π，10π = 30，因此π = 3。也许圣经只是错了，或者铸造的海并不完美圆形，或者周长是从墙外测量的而直径是从内部测量的，或者这只是诗意的许可。我们留给读者自己决定。

阿基米德（公元前 287-212 年）通过使用高阶多边形来近似圆形，从而得出了π的上下界。使用 96 边的多边形，他得出结论223/71 < π < 22/7。给出上下界在当时是一种相当复杂的方法。如果将他的两个界的平均值作为最佳估计，我们得到了3.1418，误差约为0.0002。不错！但大约 700 年后，中国数学家祖冲之使用一个有 24,576 个边的多边形得出结论3.1415962 < π < 3.1415927\。大约 800 年后，荷兰制图师阿德里安·安东尼兹（1527-1607）估算为355/113，大约为 3.1415929203539825。这一估算对大多数实际用途来说已经足够好，但并未阻止数学家继续研究这个问题。

在计算机发明之前，法国数学家布丰（1707-1788）和拉普拉斯（1749-1827）提出使用随机模拟来估算π的值。¹³¹ 想象一下在一条边长为2的正方形内画一个圆，使得圆的半径r为1。

图 18-5 单位圆被画在正方形内

根据π的定义，面积 = πr²。由于r为1，π = 面积。但圆的面积是多少呢？布丰建议通过在正方形附近抛投大量针来估算圆的面积（他认为针在下落时会随机移动）。落在正方形内的针尖数量与落在圆内的针尖数量的比率可以用来估算圆的面积。

如果针的位置是完全随机的，我们知道

并求解圆的面积，

记住，2乘以2的正方形面积是4，所以，

一般来说，要估算某个区域R的面积：

1. 1. 选择一个封闭区域E，使得E的面积容易计算，且R完全位于E内。

2. 2. 选择一组随机点，这些点位于E内。

3. 3. 设F为落在R内的点的比例。

4. 4. 将区域E的面积乘以F。

如果你尝试布丰的实验，你会很快意识到针落下的位置并不是完全随机的。此外，即使你能随机投放它们，也需要大量的针才能得到与Bible相当的π近似值。幸运的是，计算机可以以惊人的速度随机投放模拟针。¹³²

图 18-6 包含一个使用布丰-拉普拉斯方法估算π的程序。为了简化，程序仅考虑落在正方形右上象限的针。

图 18-6 估计 π

函数 throw_needles 通过首先使用 random.random 获取一对正的笛卡尔坐标（x 和 y 值），模拟扔掉一根针，这代表了针相对于正方形中心的位置。然后，它使用毕达哥拉斯定理计算以 x 为底、y 为高的直角三角形的斜边。这是针尖离原点（正方形的中心）的距离。由于圆的半径为 1，我们知道，只有当从原点的距离不大于 1 时，针才会落在圆内。我们利用这一事实来计算落在圆内的针的数量。

函数 get_est 使用 throw_needles 通过首先扔掉 num_needles 根针，然后对 num_trials 次试验的结果取平均，来找出 π 的估计值。它返回试验的均值和标准差。

函数 est_pi 以不断增加的针的数量调用 get_est，直到 get_est 返回的标准差不大于 precision/1.96。在假设误差服从正态分布的前提下，这意味着 95% 的值位于均值的 precision 范围内。

当我们运行 est_pi(0.01, 100) 时，它打印了

Est. = 3.14844, Std. dev. = 0.04789, Needles = 1000
Est. = 3.13918, Std. dev. = 0.0355, Needles = 2000
Est. = 3.14108, Std. dev. = 0.02713, Needles = 4000
Est. = 3.14143, Std. dev. = 0.0168, Needles = 8000
Est. = 3.14135, Std. dev. = 0.0137, Needles = 16000
Est. = 3.14131, Std. dev. = 0.00848, Needles = 32000
Est. = 3.14117, Std. dev. = 0.00703, Needles = 64000
Est. = 3.14159, Std. dev. = 0.00403, Needles = 128000

正如我们所预期的，随着样本数量的增加，标准差单调减少。最开始时，π 的估计值也稳步提高。有些高于真实值，有些低于真实值，但每次增加 num_needles 都会导致估计值的改善。在每次试验中使用 1000 个样本时，模拟的估计值已经优于圣经和林德纸草书。

有趣的是，当针的数量从 8,000 增加到 16,000 时，估计值反而变差，因为 3.14135 离真实的 π 值比 3.14143 更远。然而，如果我们查看每个均值周围一个标准差所定义的范围，这两个范围都包含真实的 π 值，且与较大样本量相关的范围更小。尽管使用 16,000 个样本生成的估计恰好离真实的 π 值更远，但我们应该对其准确性更有信心。这是一个极其重要的概念。仅仅给出一个好答案是不够的。我们必须有合理的理由相信它实际上是一个好答案。当我们扔掉足够多的针时，小标准差给我们提供了信心，表明我们得到了正确的答案，对吧？

不完全是。拥有小的标准差是对结果有效性有信心的必要条件，但不是充分条件。统计有效结论的概念永远不应与正确结论的概念混淆。

每个统计分析都始于一组假设。这里的关键假设是我们的模拟是现实的准确模型。回想一下，我们的布丰-拉普拉斯模拟的设计是从一些代数开始的，这些代数展示了我们如何利用两个区域的比率来找到π的值。然后我们将这个想法转化为依赖于一些几何知识和random.random随机性的代码。

让我们看看如果我们在这些方面出错会发生什么。例如，假设我们将函数throw_needles最后一行中的4替换为2，然后再次运行est_pi(0.01, 100)。这次它打印

Est. = 1.57422, Std. dev. = 0.02394, Needles = 1000
Est. = 1.56959, Std. dev. = 0.01775, Needles = 2000
Est. = 1.57054, Std. dev. = 0.01356, Needles = 4000
Est. = 1.57072, Std. dev. = 0.0084, Needles = 8000
Est. = 1.57068, Std. dev. = 0.00685, Needles = 16000
Est. = 1.57066, Std. dev. = 0.00424, Needles = 32000

仅仅32,000根针的标准差表明我们对这个估计应该有相当的信心。但这到底意味着什么呢？这意味着我们可以合理地相信，如果我们从同一分布中抽取更多样本，我们会得到一个类似的值。它并未说明这个值是否接近实际的π值。如果你要记住关于统计学的一件事，请记住这一点：统计上有效的结论不应与正确的结论混淆！

在相信模拟结果之前，我们需要对我们的概念模型的正确性以及我们是否正确实现了该模型有信心。尽可能地，你应该尝试将结果与现实进行验证。在这种情况下，你可以使用其他方法计算圆的面积的近似值（例如，物理测量），并检查计算得到的π值至少是否在正确的范围内。

18.5 关于模拟模型的一些结束语

在科学历史的大部分时间里，理论家们使用数学技术构建纯粹的解析模型，这些模型可以根据一组参数和初始条件预测系统的行为。这导致了重要数学工具的发展，从微积分到概率论。这些工具帮助科学家们对宏观物理世界形成了相对准确的理解。

随着二十世纪的进展，这种方法的局限性变得越来越明显。这些原因包括：

• 对社会科学（例如经济学）的兴趣增加，促使人们希望构建那些在数学上无法处理的系统的良好模型。

• 随着被建模的系统变得越来越复杂，似乎逐步完善一系列模拟模型比构建准确的解析模型要容易得多。

• 从模拟中提取有用的中间结果往往比从解析模型中提取更容易，例如进行“如果……会怎样”的游戏。

• 计算机的可用性使得运行大规模模拟成为可能。在二十世纪中叶现代计算机出现之前，模拟的实用性受限于手动计算所需的时间。

仿真模型是描述性的，而不是处方性的。它们描述系统在给定条件下如何工作，而不是如何安排条件以使系统表现最佳。仿真并不优化，它只是描述。这并不是说仿真不能作为优化过程的一部分。例如，仿真通常作为寻找最佳参数设置的一部分搜索过程。

仿真模型可以沿三个维度分类：

• 确定性与随机性

• 静态与动态

• 离散与连续

确定性仿真的行为完全由模型定义。重新运行仿真不会改变结果。确定性仿真通常用于被建模系统本身也是确定性的情况，但分析过于复杂，例如，处理器芯片的性能。随机仿真在模型中引入了随机性。对同一模型的多次运行可能会生成不同的值。这种随机因素迫使我们生成多个结果，以查看可能性的范围。生成10、1000或100,000个结果的问题是一个统计问题，如前所述。

在静态模型中，时间没有本质作用。本章中用于估计π的针落仿真就是一个静态仿真的例子。在动态模型中，时间或某种类似物起着重要作用。在第十六章中模拟的一系列随机行走中，所采取的步数被用作时间的替代。

在离散模型中，相关变量的值是可枚举的，例如，它们是整数。在连续模型中，相关变量的值范围在不可枚举的集合上，例如，实数。想象分析高速公路上的交通流。我们可能选择对每辆汽车进行建模，在这种情况下，我们有一个离散模型。或者，我们可能选择将交通视为一种流，其中流的变化可以用微分方程描述。这就导致了一个连续模型。在这个例子中，离散模型更接近物理情况（没有人开半辆车，尽管有些车的尺寸是其他车的一半），但计算复杂性大于连续模型。在实践中，模型往往同时具有离散和连续组件。例如，我们可能选择使用离散模型对血液流动进行建模（即，对单个血球建模），并使用连续模型对血压进行建模。

18.6 本章引入的术语

• 蒙特卡洛仿真

• 投资回报率（ROI）

• 表查找

• 时间/空间权衡

• 描述模型

• 处方模型

• 确定性仿真

• 随机仿真

• 静态模型

• 动态模型

• 离散模型

• 连续模型

第十九章：抽样与置信度

回想一下，推断统计涉及通过分析随机选择的样本来对总体进行推断。这个样本被称为样本。

抽样很重要，因为通常无法观察到整个感兴趣的总体。医生无法计算患者血液中某种细菌的数量，但可以测量患者血液的小样本中的细菌数量，从而推断整个总体的特征。如果你想知道十八岁美国人的平均体重，可以尝试把他们全部召集起来，放在一个非常大的秤上，然后除以人数。或者，你可以随机召集 50 个十八岁的人，计算他们的平均体重，并假设这个平均体重是整个十八岁人群平均体重的合理估计。

样本与感兴趣总体之间的对应关系至关重要。如果样本不能代表总体，再复杂的数学也无法得出有效的推断。50 个女性或 50 个亚裔美国人或 50 个足球运动员的样本不能用于对美国所有十八岁人群的平均体重进行有效推断。

在本章中，我们关注概率抽样。使用概率抽样，感兴趣总体的每个成员都有非零的被纳入样本的概率。在简单随机样本中，总体的每个成员被选中的机会是均等的。在分层抽样中，总体首先被划分为子群体，然后从每个子群体随机抽样以构建样本。分层抽样可以提高样本代表整个总体的概率。例如，确保样本中男性和女性的比例与总体中的比例相符，会增加样本均值（样本均值）是整个总体均值（总体均值）的良好估计的概率。

本章中的代码假设以下导入语句

import random
import numpy as np
import matplotlib.pyplot as plt
import scipy

19.1 抽样波士顿马拉松

自 1897 年以来，每年都有运动员（主要是跑步者，但自 1975 年以来有残疾人组别）聚集在马萨诸塞州参加波士顿马拉松。¹³³近年来，每年约有 20,000 名勇敢的人成功挑战42.195公里（26英里，385码）的赛道。

包含 2012 年比赛数据的文件可在与本书相关的网站上获得。文件bm_results2012.csv为逗号分隔格式，包含每位参与者的姓名、性别、¹³⁴年龄、组别、国家和时间。图 19-1 包含该文件内容的前几行。

图 19-1 bm_results2012.csv中的前几行

由于每场比赛的完整结果数据很容易获得，因此没有实际需要使用抽样来推导比赛的统计数据。然而，从教育角度来看，将从样本中得出的统计估计与实际估计值进行比较是有益的。

图 19-2 中的代码生成了图 19-3 中显示的图表。函数get_BM_data从包含比赛中每位竞争者信息的文件中读取数据。它返回一个包含六个元素的字典。每个键描述与该键关联的列表中元素的数据类型（例如，'name'或'gender'）。例如，data['time']是一个浮点数列表，包含每个竞争者的完成时间，data['name'][i]是第i位竞争者的名字，data['time'][i]是第i位竞争者的完成时间。函数make_hist生成完成时间的可视化表示。（在第二十三章中，我们将研究一个可以简化本章中很多代码的 Python 模块 Pandas，包括get_BM_data和make_hist。）

图 19-2 读取数据并生成波士顿马拉松的图表

代码

times = get_BM_data('bm_results2012.csv')['time']
make_hist(times, 20, '2012 Boston Marathon',
          'Minutes to Complete Race', 'Number of Runners')

在图 19-3 中生成图表。

图 19-3 波士顿马拉松完成时间

完成时间的分布类似于正态分布，但由于右侧的粗尾明显不正常。

现在，让我们假装没有关于所有竞争者的数据，而是想通过抽样一小部分随机选择的竞争者来估计整个参赛者完成时间的一些统计数据。

图 19-4 中的代码创建了times元素的简单随机样本，然后使用该样本来估计times的均值和标准差。函数sample_times使用random.sample(times, num_examples)来提取样本。调用random.sample返回一个大小为num_examples的列表，包含从列表times中随机选择的不同元素。在提取样本后，sample_times生成一个直方图，显示样本中值的分布。

图 19-4 抽样完成时间

正如图 19-5 所示，样本的分布与其抽取的分布相比，远离正态分布。这并不令人惊讶，因为样本大小较小。更令人惊讶的是，尽管样本大小为(40，约为21,000)小，但估计的均值与总体均值相差约3%。我们是运气好，还是有理由期待均值的估计会相当准确？换句话说，我们能否以量化的方式表达我们对估计的信心有多大？

图 19-5 分析小样本

正如我们在第十七章和第十八章讨论的那样，提供置信区间和置信水平以指示估计的可靠性通常是有用的。给定一个从更大总体中抽取的单个样本（任意大小），总体均值的最佳估计是样本的均值。估计所需达到期望置信水平的置信区间宽度则更复杂。这在一定程度上取决于样本的大小。

很容易理解样本大小为什么重要。大数法则告诉我们，随着样本大小的增加，样本值的分布更可能类似于其抽取的总体的分布。因此，随着样本大小的增加，样本均值和样本标准差更可能接近总体均值和总体标准差。

所以，越大越好，但多大才算足够？这取决于总体的方差。方差越高，需要的样本越多。考虑两个正态分布，一个均值为0，标准差为1，另一个均值为0，标准差为100。如果我们从这些分布中随机选择一个元素并用它来估计该分布的均值，那么该估计在任何期望精度∈内的真实均值(0)的概率，将等于在概率密度函数下方，范围在−∈和∈之间的面积（见第 17.4.1 节）。图 19-6 中的代码计算并打印了对于∈ = 3 分钟的这些概率。

图 19-6 方差对均值估计的影响

当运行图 19-6 中的代码时，它会打印出

Probability of being within 3 of true mean of tight dist. = 0.9973
Probability of being within 3 of true mean of wide dist. = 0.0239

图 19-7 中的代码绘制了来自两个正态分布的1000个样本，每个样本大小为40的均值。再次强调，每个分布的均值为0，但一个的标准差为1，另一个的标准差为100。

图 19-7 计算并绘制样本均值

图 19-8 的左侧显示了每个样本的均值。正如预期的那样，当总体标准差为 1 时，样本均值都接近总体均值0，这就是为什么看不到明显的圆圈——它们密集到合并成看似一条柱形的状态。相反，当总体标准差为100时，样本均值则以难以辨认的模式分散。

图 19-8 样本均值

然而，当我们查看标准差为100时的均值直方图时，在图 19-8 的右侧，出现了一个重要的现象：均值形成了一个类似于以0为中心的正态分布。右侧的图 19-8 看起来如此并非偶然。这是中心极限定理的结果，这是所有概率和统计中最著名的定理。

19.2 中心极限定理

中心极限定理解释了为什么可以使用从一个总体中抽取的单一样本来估计从同一总体中抽取的一组假设样本均值的变异性。

中心极限定理（简称CLT）的一个版本最早由拉普拉斯于 1810 年发表，并在 1820 年代由泊松进一步完善。但我们今天所知的 CLT 是 20 世纪上半叶一系列杰出数学家工作的成果。

尽管（或许正因为）有许多杰出的数学家参与其中，CLT 实际上相当简单。它的意思是

• 给定从同一总体中抽取的足够大的样本集，样本的均值（样本均值）将近似呈正态分布。

• 这种正态分布的均值将接近总体的均值。

• 样本均值的方差（使用numpy.var计算）将接近于总体方差除以样本大小。

让我们来看一个中心极限定理（CLT）实际应用的例子。想象一下，你有一个骰子，每次掷出的结果会产生一个在 0 到 5 之间的随机实数。图 19-9 中的代码模拟了多次掷骰子的过程，打印均值和方差（variance函数在图 17-8 中定义），并绘制了显示各个数字范围出现概率的直方图。它还模拟了多次掷100个骰子，并在同一图中绘制这些100个骰子的均值的直方图。hatch关键字参数用于在视觉上区分两个直方图。

图 19-9 估计连续骰子的均值

weights关键字绑定到与hist的第一个参数相同长度的数组，用于给第一个参数中的每个元素分配权重。在生成的直方图中，箱子中的每个值都贡献其相关的权重到箱子计数中（而不是通常的1）。在这个例子中，我们使用weights来缩放 y 值，以反映每个箱子的相对（而非绝对）大小。因此，对于每个箱子，y 轴上的值是均值落在该箱子内的概率。

运行代码后，生成了图 19-10 中的绘图，并打印了，

图 19-10 中心极限定理的示意图

Mean of rolling 1 die = 2.5003 Variance = 2.0814
Mean of rolling 100 dice = 2.4999 Variance = 0.0211

在每种情况下，均值都非常接近预期的均值2.5。由于我们的骰子是公平的，一个骰子的概率分布几乎是完全均匀的，¹³⁵即远非正态。然而，当我们查看100个骰子的平均值时，分布几乎是完全正态的，峰值包括预期的均值。此外，100次掷骰子的均值方差接近单次掷骰子的值除以100的方差。所有结果都如中心极限定理所预测。

中心极限定理似乎有效，但它有什么用呢？或许它可以帮助那些在特别书呆子的酒吧喝酒的人赢得酒吧赌注。然而，中心极限定理的主要价值在于它允许我们计算置信水平和区间，即使基础的总体分布不是正态分布。当我们在第 17.4.2 节中讨论置信区间时，指出经验法则是基于对所采样空间性质的假设。我们假设

• 均值估计误差为 0。

• 估计值的误差分布是正态的。

当这些假设成立时，针对正态分布的经验法则提供了一种便捷的方法，根据均值和标准差估计置信区间和水平。

让我们回到波士顿马拉松的例子。代码在图 19-11 中生成的绘图在图 19-12 中显示，针对各种样本大小绘制了 200 个简单随机样本。对于每个样本大小，它计算了这 200 个样本的均值；然后计算这些均值的均值和标准差。由于中心极限定理告诉我们样本均值将服从正态分布，我们可以使用标准差和经验法则为每个样本大小计算95%置信区间。

图 19-11 带误差条的绘图

图 19-12 带误差条的完成时间估计

如图 19-12 所示，所有估计值都与实际总体均值相当接近。然而，请注意，估计均值的误差并不是随着样本大小单调减少——使用700个例子的估计恰好比使用50个例子的估计更差。随着样本大小的增加，我们对均值估计的信心是单调增加的。当样本大小从100增加到1500时，置信区间从大约±15减少到大约±2.5。这非常重要。仅仅运气好并得到一个好的估计是不够的。我们需要知道对我们的估计要有多少信心。

19.3 均值的标准误差

我们刚刚看到，如果选择 200 个随机样本，每个样本1,500名竞争者，我们可以以95%的置信度，在大约五分钟的范围内估计均值完成时间。我们使用了样本均值的标准差。不幸的是，由于这涉及到使用比竞争者更多的总例子(200*1500 = 300,000)，这似乎并不是一个有用的结果。我们直接使用整个总体计算实际均值会更好。我们需要一种方法，通过单个例子估计置信区间。引入均值的标准误差（SE或SEM）的概念。

样本大小为n的 SEM 是从同一总体中抽取的无穷多个样本均值的标准差，样本大小为n。不出所料，它依赖于n和σ，总体的标准差：

图 19-13 将图 19-12 中使用的样本大小的 SEM 与我们为每个样本大小生成的 200 个样本的均值标准差进行了比较。

图 19-13 标准误差

我们的 200 个样本的均值的实际标准差与标准误差(SE)紧密相关。注意，SEM 和 SD 在开始时迅速下降，然后随着样本大小的增大而减缓。这是因为该值依赖于样本大小的平方根。换句话说，要将标准差减半，我们需要将样本大小增加四倍。

可惜的是，如果我们只有一个样本，就不知道总体的标准差。通常，我们假设样本的标准差是总体标准差的合理替代。这在总体分布不严重偏斜时是成立的。

图 19-14 中的代码创建了来自波士顿马拉松数据的100个不同大小的样本，并将每个大小样本的均值标准差与总体的标准差进行了比较。它生成了图 19-15 中的图。

图 19-14 样本标准偏差与总体标准偏差

图 19-15 样本标准偏差

当样本大小达到100时，样本标准偏差与总体标准偏差之间的差异相对较小（约为实际平均完成时间的 1.2%）。

实际上，人们通常使用样本标准偏差来代替（通常未知的）总体标准偏差以估计标准误差。如果样本大小足够大，¹³⁶并且总体分布与正态分布相差不大，那么使用这个估计值来计算基于经验法则的置信区间是安全的。

这意味着什么？如果我们抽取一个包含 200 名跑步者的单一样本，我们可以

• 计算该样本的均值和标准偏差。

• 使用该样本的标准偏差来估计标准误差（SE）。

• 使用估计的标准误差（SE）来生成围绕样本均值的置信区间。

图 19-16 中的代码执行此操作10,000次，然后打印样本均值与总体均值之间超过1.96个估计标准误差的次数比例。（请记住，对于正态分布，95%的数据落在均值的1.96个标准偏差范围内。）

图 19-16 估计总体均值 10,000 次

当代码运行时，它会打印，

Fraction outside 95% confidence interval = 0.0533

这基本上是理论的预测。中心极限定理胜出一分！

19.4 本章引入的术语

• 总体

• 样本

• 样本大小

• 概率抽样

• 简单随机样本

• 分层抽样

• 样本均值

• 总体均值

• 中心极限定理

• 标准误差（SE，SEM）

第二十章：理解实验数据

本章是关于理解实验数据的。我们将广泛使用绘图来可视化数据，并展示如何使用线性回归构建实验数据模型。我们还将讨论物理实验和计算实验之间的相互作用。我们将讨论如何得出有效的统计结论的内容推迟到第二十一章。

20.1 弹簧的行为

弹簧是非常奇妙的东西。当它们被某种力压缩或拉伸时，会储存能量。当这个力不再施加时，它们释放储存的能量。这一特性使它们能够在汽车中平滑行驶，帮助床垫适应我们的身体，收回安全带，发射弹丸。

1676 年，英国物理学家罗伯特·胡克制定了弹性胡克定律：Ut tensio, sic vis，用英语表达为F = -kx。换句话说，存储在弹簧中的力F与弹簧被压缩（或拉伸）的距离呈线性关系。（负号表示弹簧施加的力与位移方向相反。）胡克定律适用于各种材料和系统，包括许多生物系统。当然，它不适用于任意大的力。所有弹簧都有一个弹性极限，超过这个极限，定律就失效了。那些拉伸过度的 Slinky 玩具的人对此了解得太清楚了。

比例常数k称为弹簧常数。如果弹簧很硬（比如汽车悬挂中的弹簧或弓箭手的弓臂），k就大。如果弹簧很弱，比如圆珠笔中的弹簧，k就小。

知道特定弹簧的弹簧常数可能非常重要。简单秤和原子力显微镜的标定都依赖于知道组件的弹簧常数。DNA 链的机械行为与压缩它所需的力相关。弓发射箭矢的力与其弓臂的弹簧常数相关，等等。

代代物理学生通过使用类似于图 20-1 所示的实验装置来估计弹簧常数。

图 20-1 经典实验

我们从一个没有附加重量的弹簧开始，测量弹簧底部到支架顶部的距离。然后我们在弹簧上挂上已知质量的物体，并等待它停止移动。在此时，存储在弹簧中的力就是悬挂物体施加在弹簧上的力。这就是胡克定律中的F值。我们再次测量弹簧底部到支架顶部的距离。这个距离与挂上重量之前的距离之间的差值就是胡克定律中的x值。

我们知道施加在弹簧上的力F等于质量m乘以重力加速度g（9.81 m/s²是这个星球表面g的一个相当好的近似值），因此我们将mg代入F。通过简单的代数，我们知道k = -(mg)/x.

假设，例如，m = 1kg和x = 0.1m，那么

根据这个计算，拉伸弹簧一米需要*98.1*牛顿¹³⁷的力。

如果一切都很好

• 我们完全相信我们能够完美地进行这个实验。在这种情况下，我们可以进行一次测量，执行计算，并知道我们找到了k。不幸的是，实验科学几乎从不这样运作。

• 我们可以确保我们在弹簧的弹性极限以下进行操作。

一个更稳健的实验是将一系列越来越重的重物悬挂在弹簧上，每次测量弹簧的伸长并绘制结果。我们进行了这样的实验，并将结果输入到一个名为springData.csv的文件中：

Distance (m), Mass (kg)
0.0865,0.1
0.1015,0.15
…
0.4416,0.9
0.4304,0.95
0.437,1.0

图 20-2 中的函数从一个文件读取数据，例如我们保存的文件，并返回包含距离和质量的列表。

图 20-2 从文件中提取数据

图 20-3 中的函数使用get_data从文件中提取实验数据，然后生成图 20-4 中的图表。

图 20-3 绘制数据

图 20-4 弹簧的位移

这不是胡克定律所预测的。胡克定律告诉我们，距离应与质量线性增加，即点应位于一条直线上，其斜率由弹簧常数决定。当然，我们知道当我们进行真实测量时，实验数据很少与理论完全吻合。测量误差是可以预期的，因此我们应当期望点位于一条线附近，而不是在线上。

仍然，看到一条代表我们最佳猜测的线会很好，假如没有测量误差，点将会在何处。通常的做法是对数据进行线性拟合。

20.1.1 使用线性回归找到拟合

每当我们将任何曲线（包括直线）拟合到数据时，我们需要某种方法来判断哪条曲线是数据的最佳拟合。这意味着我们需要定义一个客观函数，以定量评估曲线与数据的拟合程度。一旦我们有了这样的函数，找到最佳拟合可以被表述为寻找一个最小化（或最大化）该函数值的曲线，即作为一个优化问题（见第 14 和 15 章）。

最常用的目标函数称为 最小二乘法。令 *observed* 和 *predicted* 为相同长度的向量，其中 *observed* 包含测量点，predicted 包含建议拟合的相应数据点。

然后定义目标函数为：

对观察值和预测值之间的差异进行平方处理，使观察值和预测值之间的大差异相对比小差异更为重要。平方差异还丢失了关于差异是正还是负的信息。

我们如何找到最佳的最小二乘拟合？一种方法是使用类似于第三章中牛顿–拉夫森算法的逐次逼近算法。或者，通常适用解析解。但我们不必实现牛顿–拉夫森或解析解，因为 numpy 提供了一个内置函数 polyfit，该函数可以找到最佳最小二乘拟合的近似值。该调用

`np.polyfit(observed_x_vals, observed_y_vals, n)`

查找提供最佳最小二乘拟合的多项式的系数，次数为 n，该多项式由两个数组 observed_x_vals 和 observed_y_vals 定义的点集给出。例如，该调用

`np.polyfit(observed_x_vals, observed_y_vals, 1)`

将找到由多项式 y = ax + b 描述的直线，其中 a 是直线的斜率，b 是 y 轴截距。在这种情况下，调用返回一个包含两个浮点值的数组。类似地，抛物线由二次方程 y = ax² + bx + c 描述。因此，该调用

`np.polyfit(observed_x_vals, observed_y_vals, 2)`

返回一个包含三个浮点值的数组。

polyfit 使用的算法称为 线性回归。这可能有些令人困惑，因为我们可以用它来拟合除直线以外的曲线。一些作者确实区分线性回归（当模型是直线时）和 多项式回归（当模型是次数大于 1 的多项式时），但大多数作者没有。¹³⁸

图 20-5 中的函数 fit_data 通过添加一条表示数据最佳拟合的直线来扩展 图 20-3 中的 plot_data 函数。它使用 polyfit 找到系数 a 和 b，然后利用这些系数生成每个力的预测弹簧位移。注意，forces 和 distance 的处理方式存在不对称性。forces 中的值（来自悬挂在弹簧上的质量）被视为独立变量，用于生成因悬挂该质量而产生的因变量 predicted_distances 的值。

图 20-5 拟合数据曲线

该函数还计算弹簧常数 k。直线的斜率 a 是 Δdistance/Δforce。另一方面，弹簧常数 k 是 Δforce/Δdistance。因此，k 是 a 的倒数。

调用fit_data('springData.csv')生成图 20-6 中的图表。

图 20-6 测量点和线性模型

有趣的是，实际上很少有点落在最小二乘拟合上。这是合理的，因为我们试图最小化平方误差的总和，而不是最大化落在直线上的点的数量。不过，这看起来似乎并不是一个很好的拟合。让我们通过向fit_data添加代码来尝试三次拟合

#find cubic fit
fit = np.polyfit(forces, distances, 3)
predicted_distances = np.polyval(fit, forces)
plt.plot(forces, predicted_distances, 'k:', label = 'cubic fit')

在这段代码中，我们使用了函数polyval来生成与三次拟合相关的点。这个函数接受两个参数：一组多项式系数和一组用于计算多项式值的自变量。代码片段

fit = np.polyfit(forces, distances, 3)
predicted_distances = np.polyval(fit, forces)

和

a,b,c,d = np.polyfit(forces, distances, 3)
predicted_distances = a*(forces**3) + b*forces**2 + c*forces + d

是等价的。

这生成了图 20-7 中的图表。三次拟合似乎比线性拟合更好地描述了数据，但真的如此吗？可能不是。

图 20-7 线性与三次拟合

技术性和大众化文章经常包含这样的图表，展示原始数据和与数据拟合的曲线。然而，作者常常假设拟合曲线就是对实际情况的描述，而原始数据仅仅是实验误差的指示。这是很危险的。

请记住，我们开始时的理论是x和y值之间应该存在线性关系，而不是三次关系。让我们看看如果使用我们的线性和三次拟合来预测对应于悬挂1.5kg重量的点会落在哪里，图 20-8。

图 20-8 使用模型进行预测

现在三次拟合看起来不太好。特别是，通过在弹簧上挂一个大重量，弹簧会升高到（y 值为负）其悬挂的杆上方，这似乎极不可能。我们遇到的是过拟合的例子。过拟合通常发生在模型过于复杂时，例如，相对于数据量，它有太多参数。当这种情况发生时，拟合可以捕捉到数据中的噪声，而不是有意义的关系。过拟合的模型通常具有较差的预测能力，正如这个例子所示。

手指练习： 修改图 20-5 中的代码，以便生成图 20-8 中的图表。

让我们回到线性拟合。此刻，忘掉直线，研究原始数据。它有什么奇怪的地方吗？如果我们对最右侧的六个点拟合一条直线，它将几乎与 x 轴平行。这似乎与胡克定律相矛盾——直到我们记起胡克定律仅在某个弹性极限内有效。也许这个极限在7N（大约0.7kg）附近。

让我们看看如果通过替换fit_data的第二和第三行来消除最后六个点会发生什么。

`distances = np.array(distances[:-6]) masses = np.array(masses[:-6])`

正如图 20-9 所示，去掉那些点确实会产生影响：k值显著下降，线性和立方拟合几乎无法区分。但我们怎么知道哪条线性拟合更好地代表了我们的弹簧在其弹性极限内的表现呢？我们可以使用某种统计检验来确定哪条线更适合数据，但那并不是重点。这个问题无法通过统计来回答。毕竟，我们可以抛弃所有数据，只保留任意两个点，并知道polyfit会找到一条完全适合这两个点的线。单纯为了获得更好的拟合而抛弃实验结果是完全不合适的。¹³⁹ 在这里，我们通过引用胡克定律的理论，即弹簧具有弹性极限，来合理化抛弃最右侧的点。这个理由不应该被用于其他地方的数据点。

图 20-9 弹性极限的模型

20.2 发射体的行为

对于仅仅拉伸弹簧感到厌倦，我们决定用我们的一个弹簧制作一个能够发射发射体的装置。¹⁴⁰ 我们使用该装置四次，将发射体发射到距离发射点30码（1080英寸）的目标上。每次，我们都测量发射体在距离发射点不同距离时的高度。发射点和目标在同一高度，我们在测量中将其视为0.0。

数据存储在一个文件中，部分数据如图 20-10 所示。第一列包含发射体距离目标的距离。其他列包含在该距离下四次实验中发射体的高度。所有测量均以英寸为单位。

图 20-10 来自发射体实验的数据

图 20-11 中的代码用于绘制四次实验中发射体的平均高度与发射点之间的距离关系。它还绘制了这些点的最佳线性和二次拟合。（如果你忘记了将列表与整数相乘的含义，表达式[0]*len(distances)将生成一个包含len(distances)个0的列表。）

图 20-11 绘制发射体的轨迹

通过快速查看图 20-12 中的图表，很明显，二次拟合远优于线性拟合。¹⁴¹ 但从绝对意义上讲，这条线的拟合有多糟糕，二次拟合又有多好呢？

图 20-12 轨迹图

20.2.1 决定系数

当我们为一组数据拟合曲线时，我们是在寻找一个将自变量（本例中从发射点水平距离的英寸数）与因变量的预测值（本例中发射点以上的英寸数）关联的函数。询问拟合优度相当于询问这些预测的准确性。请记住，拟合是通过最小化均方误差来找到的。这表明我们可以通过查看均方误差来评估拟合优度。该方法的问题在于，虽然均方误差有下界（0），但没有上界。这意味着尽管均方误差对于比较同一数据的两个拟合的相对优度是有用的，但它对于获取拟合的绝对优度的感觉并不特别有用。

我们可以使用决定系数计算拟合的绝对优度，通常写作R²。¹⁴² 设y[i]为第i^(th)个观察值，p[i]为模型预测的对应值，μ为观察值的均值。

通过将估计误差（分子）与原始值的变异性（分母）进行比较，R²旨在捕捉统计模型所解释的数据集中相对于均值的变异比例。当评估的模型由线性回归产生时，R²的值始终介于0和1之间。如果R²= 1，则模型与数据完全拟合。如果R²= 0，则模型预测的值与数据围绕均值的分布之间没有关系。

图 20-13 中的代码提供了这个统计度量的直接实现。它的紧凑性源于对 numpy 数组操作的表达能力。表达式(predicted - measured)**2从一个数组的元素中减去另一个数组的元素，然后对结果中的每个元素进行平方。表达式(measured - mean_of_measured)**2从数组measured的每个元素中减去标量值mean_of_measured，然后对结果中的每个元素进行平方。

图 20-13 计算 R²

当代码行

print('r**2 of linear fit =', r_squared(mean_heights, altitudes))

和

print('r**2 of quadratic fit =', r_squared(mean_heights, altitudes))

在process_trajectories中适当调用plt.plot之后插入时，它们会打印

r**2 of linear fit = 0.0177433205440769
r**2 of quadratic fit = 0.9857653692869693

粗略地说，这告诉我们，线性模型只能解释测量数据中不到2%的变异，但二次模型可以解释超过98%的变异。

20.2.2 使用计算模型

现在我们有了一个似乎是我们数据的良好模型，我们可以利用这个模型帮助回答关于我们原始数据的问题。一个有趣的问题是弹道在击中目标时的水平速度。我们可以使用以下思路来设计一个计算，以回答这个问题：

1. 1. 我们知道，弹道的轨迹由形式为y = ax² + bx + c的公式给出，即它是一个抛物线。由于每个抛物线在其顶点处是对称的，我们知道其最高点发生在发射点与目标之间的中间位置；将此距离称为xMid。因此，峰值高度 yPeak 由 yPeak = axMid² + bxMid + c 给出。

2. 2. 如果我们忽略空气阻力（记住没有模型是完美的），我们可以计算弹道从yPeak降到目标高度所需的时间，因为这纯粹是重力的函数。它由方程给出。¹⁴³这也是弹道从xMid移动到目标所需的时间，因为一旦它到达目标就停止移动。

3. 3. 给定从xMid到目标的时间，我们可以轻松计算该时间段内弹道的平均水平速度。如果我们假设弹道在该时间段内在水平方向上既没有加速也没有减速，我们可以将平均水平速度作为弹道击中目标时水平速度的估计。

图 20-14 实现了这种估算弹道水平速度的技术。¹⁴⁴

图 20-14 计算弹道的水平速度

当在process_trajectories (图 20-11)末尾插入行get_horizontal_speed(fit, distances[-1], distances[0])时，它会打印。

Horizontal speed = 136 feet/sec

我们刚刚经历的步骤序列遵循一个常见模式。

1. 1. 我们首先进行了一项实验，以获取有关物理系统行为的数据。

2. 2. 然后我们使用计算来寻找并评估系统行为模型的质量。

3. 3. 最后，我们使用一些理论和分析设计了一个简单的计算，以推导出模型的一个有趣结果。

手指练习： 在真空中，物体下落的速度由方程v = v0 + gt定义，其中v0是物体的初始速度，t是物体下落的秒数，g是重力常数，地球表面约为9.8 m/sec²，火星上为3.711 m/sec²。一位科学家在一个未知星球上测量下落物体的速度。她通过在不同时间点测量物体的下落速度来实现。在时间0时，物体的速度为未知的v0。实现一个函数，将模型拟合到时间和速度数据上，并估计该星球的g和实验的v0。它应该返回g和v0的估计值，以及模型的 r 平方值。

20.3 拟合指数分布数据

Polyfit使用线性回归来找到某个给定次数的多项式，该多项式是某些数据的最佳最小二乘拟合。如果数据可以直接用多项式近似，它工作得很好。但这并不总是可能。例如，考虑简单的指数增长函数y = 3^x。图 20-15 中的代码拟合了一个五次多项式到前十个点，并按图 20-16 中的结果绘制。它使用函数调用np.arange(10)，返回一个包含整数0-9的array。参数设置markeredgewidth = 2设置了标记中使用的线条宽度。

图 20-15 拟合指数分布的多项式曲线

图 20-16 拟合指数分布

拟合对于这些数据点显然是好的。然而，让我们看看模型对3²⁰的预测。当我们添加代码时

print('Model predicts that 3**20 is roughly',
      np.polyval(fit, [3**20])[0])
print('Actual value of 3**20 is', 3**20)

到图 20-15 的末尾，它打印出，

Model predicts that 3**20 is roughly 2.4547827637212492e+48
Actual value of 3**20 is 3486784401

哎呀！尽管拟合了数据，但polyfit产生的模型显然不好。这是因为5不是正确的次数吗？不。这是因为没有多项式能很好地拟合指数分布。这是否意味着我们无法使用polyfit来建立指数分布的模型？幸运的是，不是这样，因为我们可以使用polyfit找到一个适合原始独立值和依赖值对数的曲线。

考虑指数序列[1, 2, 4, 8, 16, 32, 64, 128, 256, 512]。如果我们对每个值取以2为底的对数，我们得到序列[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]，即线性增长的序列。实际上，如果一个函数y = f(x)表现出指数增长，则f(x)的对数（以任何底数）都是线性增长的。这可以通过绘制具有对数 y 轴的指数函数来可视化。代码

x_vals, y_vals = [], []
for i in range(10):
    x_vals.append(i)
    y_vals.append(3**i)
plt.plot(x_vals, y_vals, 'k')
plt.semilogy()

产生的图在图 20-17 中。

图 20-17 半对数图上的指数

对指数函数取对数会产生线性函数这一事实可以用来构建指数分布的数据点模型，如图 20-18 中的代码所示。我们使用 polyfit 找到适合 x 值和 y 值对数的曲线。请注意，我们还使用了另一个 Python 标准库模块 math，它提供了一个 log 函数。（我们本可以使用 np.log2，但想指出 math 有一个更通用的对数函数。）

图 20-18 使用 polyfit 拟合指数函数

运行时，代码

x_vals = range(10)
f = lambda x: 3**x
y_vals = create_data(f, x_vals)
plt.plot(x_vals, y_vals, 'ko', label = 'Actual values')
fit, base = fit_exp_data(x_vals, y_vals)
predictedy_vals = []
for x in x_vals:
    predictedy_vals.append(base**np.polyval(fit, x))
plt.plot(x_vals, predictedy_vals, label = 'Predicted values')
plt.title('Fitting an Exponential Function')
plt.legend(loc = 'upper left')
#Look at a value for x not in original data
print('f(20) =', f(20))
print('Predicted value =', int(base**(np.polyval(fit, [20]))))

生成的图在图 20-19 中，实际值与预测值重合。此外，当模型在一个未用于生成拟合值的值（20）上进行测试时，它输出为

f(20) = 3486784401
Predicted value = 3486784401 

图 20-19 指数函数的拟合

使用 polyfit 找到数据模型的方法适用于关系可以用形式为 y = base^(ax+b) 的方程描述的情况。如果用于不适合此模型的数据，将会产生糟糕的结果。

为了看到这一点，我们使用

f = lambda x: 3**x + x

该模型现在的预测效果很差，输出为

f(20) = 3486784421
Predicted value = 2734037145

20.4 当理论缺失时

在本章中，我们强调了理论、实验和计算科学之间的相互作用。然而，有时我们发现自己拥有大量有趣的数据，却几乎没有理论。在这种情况下，我们通常 resort to 使用计算技术通过构建一个似乎符合数据的模型来发展理论。

在理想情况下，我们会进行一次受控实验（例如，从弹簧上挂重物），研究结果，然后回顾性地制定与结果一致的模型。接着我们会进行新的实验（例如，从同一个弹簧上挂不同的重物），并将实验结果与模型预测进行比较。

不幸的是，在许多情况下，甚至无法进行一次受控实验。举个例子，假设建立一个模型来揭示利率如何影响股价。我们中的很少有人有能力设定利率并观察结果。另一方面，相关的历史数据却是源源不断的。

在这种情况下，我们可以通过将现有数据分为训练集和保留集来模拟一组实验，用作测试集。在不查看保留集的情况下，我们建立一个似乎能解释训练集的模型。例如，我们找到一个对训练集具有合理 R² 的曲线。然后我们在保留集上测试该模型。大多数情况下，该模型会比保留集拟合训练集得更好。但是如果模型是好的，它应该能合理地拟合保留集。如果没有，那么该模型可能应该被舍弃。

我们如何选择训练集？我们希望它能代表整个数据集。一种方法是随机选择训练集的样本。如果数据集足够大，这通常效果不错。

检查模型的另一种相关但稍微不同的方法是对原始数据的多个随机选择子集进行训练，并查看模型之间的相似程度。如果它们非常相似，那么我们可以比较有信心。这种方法称为交叉验证。

交叉验证在第二十一章和第二十四章中有更详细的讨论。

20.5 在章节中引入的术语

• 胡克定律

• 弹性极限

• 弹簧常数

• 曲线拟合

• 最小二乘法

• 线性回归

• 多项式回归

• 过拟合

• 拟合优度

• 决定系数 (R²)

• 训练集

• 测试集

• 保留集

• 交叉验证

第二十一章：随机试验与假设检验

X 博士发明了一种药物 PED-X，旨在帮助职业自行车手骑得更快。当他试图推销时，自行车手们坚持要求 X 博士证明他的药物优于 PED-Y，这种他们使用多年的禁药。X 博士从一些投资者那里筹集资金并启动了随机试验。

他说服200名职业自行车手参与试验。然后他将他们随机分为两个组：治疗组和对照组。治疗组的每位成员都接受了 PED-X 的剂量。对照组的成员被告知他们接受了 PED-X 的剂量，但实际上他们接受的是 PED-Y 的剂量。

每位自行车手被要求尽可能快地骑行50英里。每组的完成时间呈正态分布。治疗组的平均完成时间为118.61分钟，而对照组为120.62分钟。图 21-1 显示了每位自行车手的时间。

图 21-1 自行车手的完成时间

X 博士非常兴奋，直到遇到一位统计学家，她指出两个组中几乎总会有一个组的均值低于另一个组，或许均值之间的差异仅仅是随机现象。当她看到科学家失落的表情时，统计学家提出要教他如何检查研究的统计显著性。

21.1 检查显著性

在任何涉及从人群中随机抽样的实验中，观察到的效应总有可能纯粹是偶然的。图 21-2 可视化了 2020 年 1 月的温度与 1981 年至 2010 年 1 月的平均温度的变化。现在，想象你通过选择地球上的 20 个随机地点来构建一个样本，然后发现样本的平均温度变化为+1摄氏度。观察到的平均温度变化是因为你恰好抽取的地点的伪影，而不是表明整个星球正在变暖的概率是多少？回答这种问题就是统计显著性的核心所在。

图 21-2 2020 年 1 月与 1981-2010 年平均气温差异¹⁴⁵

在 20 世纪初，罗纳德·费舍尔开发了一种统计假设检验的方法，已成为评估观察到的效果纯粹是偶然发生的概率的最常见方法。费舍尔说他发明了这种方法，是回应布里斯托-罗奇博士的一个声明，她声称当她喝加奶茶时，可以判断是茶先倒入茶杯还是牛奶。费舍尔向她挑战进行“茶测试”，她被提供了八杯茶（每种加茶和牛奶的顺序各四杯），要求识别那些茶先倒入的杯子。她完美地完成了。费舍尔随后计算了她纯粹偶然做到这一点的概率。正如我们在第 17.4.4 节中看到的，，即有70种方法从8杯中选择4杯。由于只有这70种组合中有一种包含所有“茶先倒入”的4杯，费舍尔计算出布里斯托-罗奇博士纯靠运气选择正确的概率是。由此，他得出结论，她的成功极不可能归因于运气。

费舍尔的显著性检验方法可以总结为

1. 1. 陈述原假设和备择假设。原假设是“处理”没有有趣的效果。对于“茶测试”，原假设是布里斯托-罗奇博士无法品尝出差异。备择假设只有在原假设为假时才能成立，例如，布里斯托-罗奇博士能够品尝出差异。¹⁴⁶

2. 2. 理解被评估样本的统计假设。对于“茶测试”，费舍尔假设布里斯托-罗奇博士为每杯茶做出了独立的决定。

3. 3. 计算相关的检验统计量。在这种情况下，检验统计量是布里斯托-罗奇博士给出的正确答案的比例。

4. 4. 在原假设下推导该检验统计量的概率。在这种情况下，这是偶然得到所有杯子的概率，即0.014。

5. 5. 决定该概率是否足够小，以至于你愿意假设原假设为假，即拒绝原假设。拒绝水平的常见值，应该提前选择，是0.05和0.01。

回到我们的骑行者，假设治疗组和对照组的时间是从 PED-X 用户和 PED-Y 用户的无限完成时间总体中抽取的样本。这个实验的原假设是这两个更大总体的均值相同，即治疗组的总体均值与对照组的总体均值之间的差为 0。备择假设是它们不相同，即均值之差不等于 0。

接下来，我们开始尝试拒绝原假设。我们为统计显著性选择一个阈值α，并尝试证明数据来自与原假设一致的分布的概率小于α。然后我们说可以以置信度α拒绝原假设，并以概率1 – α接受原假设的否定。

α的选择影响我们犯错的类型。α越大，我们越可能拒绝实际上真实的原假设。这被称为第一类错误。当α较小的时候，我们更可能接受实际上是错误的原假设。这被称为第二类错误。

通常，人们选择α = 0.05。然而，根据错误的后果，选择一个较小或较大的α可能更为合适。例如，假设原假设是，服用 PED-X 和服用 PED-Y 之间的早逝率没有差异。我们可能希望选择一个较小的α，比如0.001，作为拒绝该假设的基础，然后再决定哪种药物更安全。另一方面，如果原假设是 PED-X 和 PED-Y 在增强表现的效果上没有差异，我们可能会舒适地选择一个相对较大的α。¹⁴⁷

下一步是计算检验统计量。最常见的检验统计量是t 统计量。t 统计量告诉我们，从数据中得出的估计值与原假设之间的差异有多大，以标准误的单位来衡量。t 统计量越大，越有可能拒绝原假设。在我们的例子中，t 统计量告诉我们两个均值的差异（118.44 – 119.82 = -1.38）距离0有多少个标准误。我们的 PED-X 示例的 t 统计量是-2.11（稍后你会看到如何计算这个值）。这意味着什么？我们如何使用它？

我们使用 t 统计量的方式与使用均值的标准差数量来计算置信区间的方法非常相似（参见第 17.4.2 节）。请记住，对于所有正态分布，样本落在均值的固定标准差数量内的概率是固定的。在这里，我们做了一些稍微复杂的事情，考虑了用于计算标准误的样本数量。我们假设的是t 分布，而不是正态分布。

t 分布最早由威廉·戈塞特（William Gosset）于 1908 年描述，他是一位为亚瑟·吉尼斯与儿子酿酒厂工作的统计学家。¹⁴⁸ t 分布实际上是一系列分布，因为分布的形状取决于样本的自由度。

自由度描述了用于推导 t 统计量的独立信息量。一般来说，我们可以将自由度视为样本中可用于估计某个关于样本所抽取的总体的统计量的独立观察数量。

t 分布类似于正态分布，自由度越大，越接近正态分布。对于小自由度，t 分布的尾部比正态分布明显更胖。对于自由度在 30 或以上的情况，t 分布非常接近正态分布。

现在，让我们用样本方差来估计总体方差。考虑一个包含三个例子 100、200 和 300 的样本。回想一下

所以我们的样本方差是

看起来我们似乎在使用三个独立的信息片段，但实际上并不是。分子中的三个项并不是彼此独立的，因为这三个观察值都用于计算200名骑行者样本的均值。自由度为2，因为如果我们知道均值和三个观察值中的任何两个，则第三个观察值的值是固定的。

自由度越大，样本统计量代表总体的概率越高。从单个样本计算的 t 统计量的自由度比样本大小少一，因为在计算 t 统计量时使用了样本的均值。如果使用两个样本，则自由度比样本大小之和少二，因为在计算 t 统计量时使用了每个样本的均值。例如，对于 PED-X/PED-Y 实验，自由度为198。

在给定自由度的情况下，我们可以绘制一个显示适当 t 分布的图，并查看我们为 PED-X 示例计算的 t 统计量在分布中的位置。图 21-3 中的代码执行了这一点，并生成了图 21-4 中的图。该代码首先使用函数scipy.random.standard_t生成许多从自由度为198的 t 分布中抽取的样本。然后，它在 PED-X 样本的 t 统计量及其负值处绘制白线。

图 21-3 绘制 t 分布

图 21-4 可视化 t 统计量

白线左右的直方图区域的分数之和等于在原假设为真时，观察值至少极端的概率。

• 样本代表总体，并且

• 原假设为真。

我们需要查看两个尾部，因为我们的原假设是总体均值相等。因此，如果治疗组的均值在统计上显著大于或小于对照组的均值，检验应该失败。

在原假设成立的假设下，得到至少与观察值一样极端的值的概率被称为p 值。对于我们的 PED-X 示例，p 值是指在假设治疗组和对照组的实际总体均值相同的情况下，观察到的均值差异至少与实际差异一样大的概率。

如果原假设成立，p 值似乎告诉我们事件发生的概率，这看起来有些奇怪，因为我们通常希望原假设不成立。然而，这与经典的科学方法并没有太大区别，科学方法是基于设计能够反驳假设的实验。图 21-5 中的代码计算并打印我们的两个样本的 t 统计量和 p 值，一个样本包含对照组的时间，另一个样本包含治疗组的时间。库函数 scipy.stats.ttest_ind 执行双尾两样本t 检验并返回 t 统计量和 p 值。将参数 equal_var 设置为 False 表示我们不知道这两个总体是否具有相同的方差。

图 21-5 计算并打印 t 统计量和 p 值

当我们运行代码时，它报告

Treatment mean - control mean = -1.38 minutes
The t-statistic from two-sample test is -2.11
The p-value from two-sample test is 0.04

“是的，”X 博士兴奋地说，“看来 PED-X 不比 PED-Y 更好的概率只有 4%，因此 PED-X 有效的概率是 96%。让现金注册机开始响起来。”可惜，他的兴奋仅持续到他阅读本章的下一部分。

21.2 警惕 P 值

过于简单地将某种含义读入 p 值是非常容易的。将 p 值视为原假设为真的概率是很有诱惑力的，但这并不是它的实际含义。

零假设类似于英美刑事司法系统中的被告。该系统基于一种称为“无罪推定”的原则，即在未被证明有罪之前被认为是无罪。类似地，我们假设零假设是真实的，除非我们看到足够的证据反对它。在审判中，陪审团可以裁定被告是“有罪”或“无罪”。“无罪”裁决意味着证据不足以使陪审团相信被告“超出合理怀疑”是有罪的。¹⁴⁹可以将其视为“没有证明有罪”。“无罪”裁决并不意味着证据足以使陪审团相信被告是无罪的。它也没有说明如果陪审团看到不同的证据会得出什么结论。可以将 p 值视为陪审团裁决，其中标准“超出合理怀疑”由 α 定义，而证据则是构建 t 统计量的数据。

小 p 值表明，如果零假设为真，则特定样本不太可能出现。它类似于陪审团得出结论，认为如果被告是无罪的，就不太可能会呈现出这一组证据，因此做出了有罪的裁决。当然，这并不意味着被告实际上有罪。也许陪审团看到了误导性的证据。类似地，低 p 值可能是因为零假设实际上是错误的，或者只是因为样本并不能代表其来源的人群，即证据是误导性的。

正如你所预料的，Dr. X 坚决声称他的实验表明零假设可能是错误的。Dr. Y 坚持认为低 p 值可能是由于不具代表性的样本，并资助了与 Dr. X 的实验规模相同的另一项实验。当用她的实验样本计算统计数据时，代码打印了

Treatment mean - control mean = 0.18 minutes
The t-statistic from two-sample test is -0.27
The p-value from two-sample test is 0.78

这个 p 值比 Dr. X 的实验得到的值大了 17 倍以上，确实没有理由怀疑零假设。混乱盛行。但我们可以澄清它！

你可能不会惊讶地发现这并不是一个真实的故事——毕竟，骑自行车的人服用增强表现的药物的想法确实令人怀疑。实际上，实验的样本是通过图 21-6 中的代码生成的。

图 21-6 生成竞赛示例的代码

由于实验纯粹是计算性的，我们可以多次运行它以获得许多不同的样本。当我们生成了 10,000 对样本（每个分布各一个）并绘制 p 值的概率时，得到了图 21-7。

图 21-7 p 值的概率

由于约10%的 p 值低于0.04，因此一个实验恰好在4%水平上显示显著性并不令人感到惊讶。另一方面，第二个实验得出完全不同的结果也并不令人意外。看起来令人惊讶的是，考虑到我们知道两个分布的均值实际上是不同的，我们得到的结果在5%水平上只有约12%的时间是显著的。大约88%的时间，我们未能在5%水平上拒绝错误的原假设。

p 值可能是不可靠的指标，这也是许多出现在科学文献中的结果无法被其他科学家重复的原因之一。一个问题是，研究能力（样本的大小）与统计发现的可信度之间存在强关系。¹⁵⁰ 如果我们将样本大小增加到3000，我们仅约1%的时间未能拒绝错误的原假设。

为什么这么多研究的统计能力不足？如果我们真的在进行一项对人进行实验（而不是模拟），那么提取大小为2000的样本将比提取大小为100的样本贵 20 倍。

样本大小的问题是所谓频率主义统计方法的一个内在特征。在 21.7 节中，我们讨论一种试图缓解这一问题的替代方法。

21.3 单尾和单样本检验

到目前为止，我们在本章中只讨论了双尾双样本检验。有时，使用单尾和/或单样本t 检验更为合适。

首先，让我们考虑一个单尾双样本检验。在我们对 PED-X 和 PED-Y 相对有效性的双尾检验中，我们考虑了三种情况：1）它们同样有效，2）PED-X 比 PED-Y 更有效，以及 3）PED-Y 比 PED-X 更有效。目标是通过论证如果原假设（情况 1）为真，那么看到 PED-X 和 PED-Y 样本均值之间如此大的差异是不太可能的，从而拒绝原假设。

然而，假设 PED-X 的成本显著低于 PED-Y。为了为他的化合物找到市场，X 博士只需证明 PED-X 至少与 PED-Y 一样有效。可以这样理解，我们希望拒绝均值相等或 PED-X 均值更大的假设。请注意，这比均值相等的假设要严格弱。(假设A严格弱于假设B，如果B为真时A也为真，但反之不然。)

为此，我们从一个双样本检验开始，原始原假设由图 21-5 中的代码计算得出。它打印了

Treatment mean - control mean = -1.38 minutes
The t-statistic from two-sample test is -2.11
The p-value from two-sample test is 0.04

允许我们在约4%的水平上拒绝原假设。

我们的较弱假设怎么样？回想一下图 21-4。我们观察到在零假设成立的假设下，白线左侧和右侧直方图区域的分数之和等于获得至少与观察值一样极端的值的概率。然而，为了拒绝我们的较弱假设，我们不需要考虑左尾下的区域，因为这对应于 PED-X 比 PED-Y 更有效（一个负时间差），而我们只对拒绝 PED-X 不如 PED-Y 的假设感兴趣。也就是说，我们可以进行单尾检验。

由于 t 分布是对称的，因此为了获得单尾检验的值，我们将双尾检验的 p 值减半。因此，单尾检验的 p 值为0.02。这使我们能够在2%的水平上拒绝我们较弱的假设，而这在使用双尾检验时无法做到。

因为单尾检验提供了更强的效能来检测效应，所以每当有关于效应方向的假设时，使用单尾检验是很诱人的。但这通常不是一个好主意。仅当未检验方向上漏掉效应的后果微不足道时，单尾检验才合适。

现在让我们看看单样本检验。假设经过多年使用 PED-Y 的经验，赛车手在 PED-Y 下完成 50 英里赛道的平均时间是120分钟。为了发现 PED-X 是否与 PED-Y 有不同的效应，我们将检验零假设，即单个 PED-X 样本的平均时间等于120。我们可以使用函数scipy.stats.ttest_1samp来实现，该函数接受一个样本和与之比较的总体均值作为参数。它返回一个包含 t 统计量和 p 值的元组。例如，如果我们在图 21-5 的代码末尾附加代码。

one_sample_test = scipy.stats.ttest_1samp(treatment_times, 120)
print('The t-statistic from one-sample test is', one_sample_test[0])
print('The p-value from one-sample test is', one_sample_test[1])

它打印。

The t-statistic from one-sample test is -2.9646117910591645
The p-value from one-sample test is 0.0037972083811954023 

p 值小于我们使用双样本双尾检验得到的值并不令人惊讶。通过假设我们知道两个均值之一，我们消除了一个不确定性来源。

那么，经过这一切，我们从 PED-X 和 PED-Y 的统计分析中学到了什么？尽管 PED-X 和 PED-Y 用户的预期表现存在差异，但没有任何有限的 PED-X 和 PED-Y 用户样本能保证揭示这一差异。此外，由于预期均值的差异很小（不到千分之五），因此不太可能像 Dr. X 进行的实验（每组100名骑手）会产生足够的证据，使我们能在95%的置信水平下得出均值存在差异的结论。我们可以通过使用单尾检验来提高在95%水平上获得统计显著结果的可能性，但那将是误导性的，因为我们没有理由假设 PED-X 的效果不低于 PED-Y。

21.4 显著性与否？

林德赛和约翰在过去几年中浪费了大量时间玩一个叫做“与朋友单词”的游戏。他们之间共进行了 1273 场比赛，林德赛赢了 666 场，因此她得意洋洋地说：“我在这个游戏中比你强多了。”约翰坚称林德赛的说法毫无意义，并认为胜利的差异完全应该归因于运气。

最近读过一本关于统计学的书的约翰，提出了一种方法来判断是否合理将林德赛的相对成功归因于技能：

• 将每场1,273场比赛视为一次实验，如果林德赛获胜则返回1，否则返回0。

• 选择零假设，即这些实验的平均值为0.5。

• 对该零假设执行双尾单样本检验。

当他运行代码时

num_games = 1273
lyndsay_wins = 666
outcomes = [1.0]*lyndsay_wins + [0.0]*(num_games - lyndsay_wins)
print('The p-value from a one-sample test is',
      scipy.stats.ttest_1samp(outcomes, 0.5)[1])

它打印了

The p-value from a one-sample test is 0.0982205871243577 

促使约翰声称差异甚至没有接近5%水平的显著性。

林德赛没有学习过统计学，但读过本书第十八章，对此并不满意。“让我们进行一次蒙特卡罗模拟，”她建议，并提供了图 21-8 中的代码。

图 21-8 林德赛的游戏模拟

当林德赛的代码运行时，它打印了，

Probability of result at least this extreme by accident = 0.0491 

促使她声称约翰的统计检验完全无效，并且胜利的差异在5%水平上是统计显著的。

“不，”约翰耐心解释，“是你的模拟有问题。它假设你是更好的玩家，并进行了相当于单尾检验的操作。你模拟的内部循环是错误的。你应该进行相当于双尾检验，测试在模拟中是否有任一玩家赢得了你在实际比赛中赢得的666场比赛。”约翰随后运行了图 21-9 中的模拟。

图 21-9 游戏的正确模拟

约翰的模拟打印了

Probability of result at least this extreme by accident = 0.0986

“这与我的双尾检验预测的非常接近，”约翰得意地说。林德赛不雅的反应不适合出现在家庭导向的书中。

手指练习：一位调查记者发现，林德赛不仅使用了可疑的统计方法，还将其应用于她仅仅虚构的数据上。¹⁵¹ 实际上，约翰击败了林德赛 479 次，输掉了 443 次。这个差异在统计上有多显著？

21.5 哪个 N？

一位教授想知道上课是否与他所在系的成绩相关。他招募了40名新生，并给他们所有人都佩戴了脚踝手环，以便追踪他们的去向。一半的学生不被允许参加他们所有课程的任何讲座，¹⁵²而另一半则被要求参加所有讲座。¹⁵³ 在接下来的四年中，每位学生参加了40门不同的课程，为每组学生提供了800个成绩。

当教授对这两个样本（每个样本大小为800）的均值进行双尾 t 检验时，p 值约为0.01。这让教授感到失望，他原本希望没有统计显著性效应——这样他就能少一些因取消讲座而感到的愧疚，可以去海滩。出于绝望，他查看了两组的平均 GPA，发现差异非常小。他想知道，均值如此微小的差异为何会在那个水平上显著？

当样本大小足够大时，即使是微小的效应也可能具有高度统计显著性。N很重要，影响很大。 图 21-10 绘制了1000次试验的平均 p 值与试验中使用的样本大小的关系。对于每个样本大小和每次试验，我们生成了两个样本。每个样本均来自标准差为5的高斯分布。其中一个样本的均值为100，另一个的均值为100.5。平均 p 值随着样本大小线性下降。当样本大小达到约2000时，均值之间0.5%的差异在5%水平上变得显著，而当样本大小接近3000时，在1%水平上显著。

图 21-10 样本大小对 p 值的影响

回到我们的例子，教授在他的研究中使用N为800的每个组别是否合理？换句话说，是否真的有800个独立的样本对应每组20名学生？可能不是。每个样本有800个成绩，但只有20名学生，而与每名学生相关的40个成绩可能不应视为独立样本。毕竟，有些学生始终能获得好成绩，而有些学生的成绩则令人失望。

教授决定以不同的方式查看数据。他计算了每位学生的 GPA。当他对这两个样本（每个样本大小为20）进行双尾 t 检验时，p 值约为0.3。于是他去了海滩。

21.6 多重假设

在第十九章中，我们使用波士顿马拉松的数据进行了抽样。 图 21-11 中的代码读取了 2012 年比赛的数据，并查找来自少数国家的女性完赛时间的统计显著性差异。它使用了图 19-2 中定义的get_BM_data函数。

图 21-11 比较所选国家的平均完赛时间

当代码运行时，它打印

ITA and JPN have significantly different means, p-value = 0.025

看起来意大利或日本都可以声称其女性跑者比对方更快。¹⁵⁴然而，这样的结论将是相当脆弱的。虽然一组跑者的平均时间确实比另一组快，但样本大小（20和32）较小，可能并不能代表各国女性马拉松选手的能力。

更重要的是，我们的实验构建方式存在缺陷。我们检查了10个零假设（每对国家各一个），并发现其中一个在5%水平上可以被拒绝。可以认为我们实际上在检查零假设：“所有国家对之间的女性马拉松选手的平均完赛时间相同。”拒绝这个零假设可能没问题，但这并不等同于拒绝意大利和日本的女性马拉松选手速度相同的零假设。

这一点在图 21-12 的例子中表现得非常明显。在那个例子中，我们从同一总体中抽取 50 对样本，每个样本大小为200，并测试每对样本的均值是否存在统计差异。

图 21-12 检查多个假设

由于所有样本都来自同一总体，我们知道零假设是正确的。然而，当我们运行代码时，它打印

# of statistically significantly different (p < 0.05) pairs  = 2 

表示可以拒绝两个配对的零假设。

这并不特别令人惊讶。回想一下，p 值为0.05表示如果零假设成立，观察到两个样本之间的均值差异至少与此差异一样大的概率为0.05。因此，如果我们检查 50 对样本，其中两对样本的均值在统计上显著不同，这并不奇怪。运行大量相关实验，然后挑选你喜欢的结果，可以被温和地描述为马虎。不友好的人可能会称之为其他东西。

回到我们的波士顿马拉松实验，我们检查是否可以拒绝零假设（均值无差异）对于10对样本。当进行涉及多个假设的实验时，最简单且最保守的方法是使用邦费罗尼校正。其直觉很简单：在检查m个假设的家庭时，维持适当的家庭错误率的一种方法是以的水平测试每个单独的假设。使用邦费罗尼校正检查意大利和日本之间的差异在α = 0.05水平上是否显著，我们应该检查 p 值是否小于0.05/10，即0.005——而它并不是。

如果有许多测试或者测试统计量之间正相关，博恩费罗尼校正是保守的（即，它比必要时更频繁地未能拒绝原假设）。另一个问题是缺乏一个普遍接受的“假设家族”的定义。显然，图 21-12 中代码生成的假设是相关的，因此需要进行校正。但情况并不总是如此明确。

21.7 条件概率与贝叶斯统计

到目前为止，我们采取了所谓的频率主义统计方法。我们完全基于数据的频率或比例从样本中得出结论。这是最常用的推理框架，并导致本书前面讨论的统计假设检验和置信区间等成熟方法。从原则上讲，它的优点在于无偏。结论完全基于观察到的数据。

然而在某些情况下，另一种统计方法贝叶斯统计更为合适。考虑图 21-13 中的漫画。¹⁵⁵

图 21-13 太阳爆炸了吗？

这里发生了什么？频率主义者知道只有两种可能性：机器掷出一对六且在说谎，或者没有掷出一对六且在讲真话。由于不掷出一对六的概率为35/36（97.22%），频率主义者得出结论，机器可能在讲真话，因此太阳可能爆炸了。¹⁵⁶

贝叶斯使用额外信息来建立她的概率模型。她同意机器不太可能掷出一对六；然而，她认为这种情况发生的概率需要与太阳未爆炸的先验概率进行比较。她得出结论，太阳未爆炸的可能性甚至高于97.22%，并决定押注“太阳明天会出来”。

21.7.1 条件概率

贝叶斯推理的关键思想是条件概率。

在我们之前讨论概率时，依赖于事件独立的假设。例如，我们假设抛硬币结果为正面或反面与之前的结果无关。这在数学上是方便的，但生活并不总是这样。在许多实际情况下，独立性是一个糟糕的假设。

考虑随机选择的美国成年人是男性且体重超过197磅的概率。男性的概率约为0.5，而体重超过197磅（美国的平均体重¹⁵⁷）的概率也约为0.5。¹⁵⁸ 如果这些事件是独立的，则选定的人同时是男性且体重超过197磅的概率将是0.25。然而，这些事件并不是独立的，因为平均美国男性的体重比平均女性重约30磅。因此，更好的问题是 1）选定的人是男性的概率是多少，以及 2）在选定的人是男性的情况下，该人体重超过197磅的概率是多少？条件概率的符号使得这一点变得容易表达。

符号P(A|B)表示在假设B为真的情况下A为真的概率。通常读作“在B的条件下A的概率”。因此，公式

正好表达了我们要寻找的概率。如果P(A)和P(B)是独立的，则P(A|B) = P(A)。在上述例子中，B是男性，而*A*是体重> 197。

一般而言，如果P(B) ≠ 0，

像常规概率一样，条件概率总是在0和1之间。此外，如果Ā表示不 A，则P(A|B) + P(Ā|B) = 1。人们常常错误地假设P(A|B)等于P(B|A)。没有理由期望这种情况成立。例如，P(Male|Maltese)的值大约为0.5，但P(Maltese|Male)约为0.000064。¹⁵⁹

指尖练习： 估计随机选择的美国人同时是男性且体重超过197磅的概率。假设50%的人口是男性，男性的体重呈正态分布，平均为210磅，标准差为30磅。（提示：考虑使用经验法则。）

公式P(A|B, C)表示在B和C均为真的条件下* A为真的概率。假设B和C*相互独立，条件概率的定义和独立概率的乘法法则表明

公式P(A, B, C)代表* A、B和C*都为真的概率。

同样，P(A, B|C)表示在C条件下* A 和 B的概率。假设 A和B*相互独立

21.7.2 贝叶斯定理

假设一位四十多岁的无症状女性去做乳腺 X 光检查，并收到了坏消息：乳腺 X 光检查结果“阳性”。¹⁶⁰

患有乳腺癌的女性在乳腺 X 光检查中得到真实阳性结果的概率是0.9。没有乳腺癌的女性在乳腺 X 光检查中得到误报的概率是0.07。

我们可以使用条件概率来表达这些事实。设

canc = has breast cancer
TP = true positive
FP = false positive

使用这些变量，我们写出条件概率

P(TP | canc) = 0.9
P(FP | not Canc) = 0.07

考虑到这些条件概率，四十多岁女性在阳性乳腺 X 光检查后该多担心？她实际患乳腺癌的概率是多少？是0.93吗，因为误报率是7%？更多？更少？

这是一个技巧性问题：我们没有提供足够的信息让你以合理的方式回答问题。要做到这一点，你需要知道四十多岁女性患乳腺癌的先验概率。四十多岁女性患乳腺癌的比例是0.008 (1000 中有 8)。因此，她们没有乳腺癌的比例为1 – 0.008 = 0.992：

P(canc | woman in her 40s) = 0.008
P(not canc | woman in her 40s) = 0.992

现在我们拥有了解决四十多岁女性该多担心的所有信息。要计算她患乳腺癌的概率，我们使用称为贝叶斯定理的东西¹⁶¹（通常称为贝叶斯法则或贝叶斯规则）：

在贝叶斯世界中，概率测量的是信念的程度。贝叶斯定理将考虑证据前后的信念程度联系起来。等号左侧的公式P(A|B)是后验概率，即在考虑B后对A的信念程度。后验概率以先验P(A)和证据B对A提供的支持来定义。支持是B在A成立时的概率与B独立于A的概率之比，即。

如果我们使用贝叶斯定理来估计这位女性实际患乳腺癌的概率，我们会得到（其中canc在我们的贝叶斯定理中扮演A的角色，pos则扮演B的角色）

阳性检测的概率是

所以

也就是说，约90%的阳性乳腺 X 光检查结果是误报。贝叶斯定理在这里帮助了我们，因为我们准确估计了四十多岁女性患乳腺癌的先验概率。

请记住，如果我们开始时使用了错误的先验，将该先验纳入我们的概率估计中，会使估计结果更糟而非更好。例如，如果我们开始时的先验是

P(canc | women in her 40's) = 0.6

我们会得出误报率大约为5%，即，四十多岁女性在阳性乳腺 X 光检查中患乳腺癌的概率大约为0.95。

指尖练习： 你正在森林中漫步，看到一片看起来美味的蘑菇田。你用篮子装满了它们，准备回家烹饪并给丈夫端上。不过，在你烹饪之前，他要求你查阅一本关于当地蘑菇种类的书，检查它们是否有毒。书中说当地森林中 80%的蘑菇是有毒的。然而，你将你的蘑菇与书中的图片进行比较，决定你有 95%的把握认为这些蘑菇是安全的。你在给丈夫端上这些蘑菇时应该有多安心（假设你宁愿不成为寡妇）？

21.8 章节中引入的术语

• 随机试验

• 治疗组

• 对照组

• 统计显著性

• 假设检验

• 零假设

• 备择假设

• 检验统计量

• 假设拒绝

• 第一类错误

• 第二类错误

• t 统计量

• t 分布

• 自由度

• p 值

• 科学方法

• t 检验

• 研究的效能

• 双尾 p 检验

• 单尾 p 检验

• 研究的组别

• 樱桃采摘

• 邦费罗尼校正

• 家族型错误率

• 频率统计

• 贝叶斯统计

• 条件概率

• 真阳性

• 假阳性

• 先验概率

• 贝叶斯定理

• 信度

• 后验概率

• 先验

• 支持

第二十二章：谎言、可怕的谎言和统计数据

“如果你无法证明你想证明的东西，就展示其他东西并假装它们是同一回事。在统计与人类思维碰撞后的迷茫中，几乎没有人会注意到差异。”**¹⁶³

任何人都可以通过简单地编造虚假的统计数据来撒谎。用准确的统计数据讲谎话更具挑战性，但仍然不难。

统计思维是相对较新的发明。在大多数记录历史中，事物是以定性而非定量的方式评估的。人们必然对一些统计事实有直观的认识（例如，女性通常比男性矮），但他们没有数学工具可以从轶事证据得出统计结论。这种情况在十七世纪中叶开始改变，尤其是约翰·格朗特的《自然与政治观察——关于死亡率的记录》的出版。这部开创性的著作利用统计分析从死亡记录中估算伦敦人口，并试图提供一个可用于预测瘟疫传播的模型。

可惜，自那时以来，人们在使用统计数据时，既用于误导也用于告知。有些人故意使用统计数据进行误导；其他人则仅仅是无能。在本章中，我们讨论了一些人可能被引导到从统计数据中得出不恰当推论的方式。我们相信你会仅将这些信息用于善良的目的——成为更好的消费者和更诚实的统计信息传播者。

22.1 垃圾进垃圾出（GIGO）

“我曾两次被[国会议员]问到，‘请问，巴贝奇先生，如果你把错误的数据输入机器，能否得到正确的答案？’我无法正确理解引发这样问题的思维混乱。”——查尔斯·巴贝奇**¹⁶⁴

这里传达的信息很简单。如果输入数据严重缺陷，任何数量的统计调整都无法产生有意义的结果。

1840 年的美国人口普查显示，自由黑人和混血者的精神病发生率大约是被奴役的黑人和混血者的十倍。结论显而易见。正如美国参议员（前副总统和未来国务卿）约翰·C·卡尔霍恩所说：“这个普查所揭示的精神病数据是无可争议的。我们的国家必须得出结论，废除奴隶制对非洲人而言将是一个诅咒。”更别提很快就清楚普查数据充满了错误。正如卡尔霍恩据说向约翰·昆西·亚当斯解释的那样：“错误太多，以至于互相抵消，导致的结论就像所有数据都正确一样。”

卡尔霍恩对亚当斯的（或许是故意的）错误反应基于一个经典错误，即独立性假设。如果他的数学素养更高，他可能会说：“我相信测量误差是无偏的，且相互独立，因此在均值两侧均匀分布。”事实上，后来的分析显示，误差严重偏向，以至于无法得出统计有效的结论。¹⁶⁵

垃圾进垃圾出（GIGO）在科学文献中是一个特别有害的问题——因为它难以被检测到。2020 年 5 月，世界上最负盛名的医学期刊之一（柳叶刀）发表了一篇关于当时肆虐的 Covid-19 大流行的论文。该论文依赖于来自六大洲近 700 家医院的 96,000 名患者的数据。在审核过程中，审稿人检查了论文中报告的分析的合理性，但没有检查分析所依据的数据的合理性。发表不到一个月后，因发现数据存在缺陷，该论文被撤回。

22.2 测试是不完美的

每个实验都应被视为一个潜在存在缺陷的测试。我们可以对化学物质、现象、疾病等进行测试。然而，我们测试的事件不一定与测试结果相同。教授们设计考试的目的是为了了解学生对某些学科内容的掌握情况，但考试结果不应与学生实际理解的程度混淆。每个测试都有固有的错误率。想象一下，一个学习第二语言的学生被要求学习100个单词的含义，但他只学习了 80 个单词的含义。他的理解率是80%，但他在一个包含20个单词的测试中得分 80%的概率当然不是1。

测试可能同时存在假阴性和假阳性。如我们在第 21.7 节中看到的，阴性的乳腺 X 光检查并不能保证没有乳腺癌，而阳性的乳腺 X 光检查也不能保证其存在。此外，测试概率和事件概率并不是同一回事。当测试稀有事件（例如，稀有疾病的存在）时，这一点尤其相关。如果假阴性的代价很高（例如，错过了一个严重但可治愈的疾病），则测试应设计得具有高敏感性，即使因此导致许多假阳性。

22.3 图片可能具有误导性

毫无疑问，图形在快速传达信息方面是非常有用的。然而，当使用不当（或恶意使用）时，图表可能会产生极具误导性的效果。例如，考虑图 22-1 中描绘的美国中西部州的房价图表。

图 22-1 美国中西部的房价

从图 22-1 左侧的图表来看，似乎在 2006-2009 年间房价相对稳定。但等一下！在 2008 年底，美国住宅房地产崩溃并引发全球金融危机不是吗？确实如此，正如右侧图表所示。

这两个图表展示了完全相同的数据，但传达了截然不同的印象。左侧的图表旨在给人房价稳定的印象。在 y 轴上，设计者使用了一个范围极低的房屋平均价格为$1,000到一个极高的平均价格$500,000。这最小化了房价变动区域所占空间，给人变化相对较小的印象。右侧的图表则旨在展示房价的剧烈波动，随后崩溃。设计者使用了一个狭窄的价格范围，从而夸大了变化的幅度。

图 22-2 中的代码生成了我们之前查看的两个图表，以及一个旨在准确反映房价变动的图表。它使用了我们尚未见过的两种绘图工具。

图 22-2 绘制房价

调用plt.bar(quarters, prices, width)会生成一个条形图，条形宽度为给定值。条形的左边缘对应列表quarters的元素值，而条形的高度对应列表prices的相应元素值。函数调用plt.xticks(quarters+width/2, labels)描述了与条形相关的标签。第一个参数指定每个标签的放置位置，第二个参数则是标签的文本。函数yticks的行为类似。调用plot_housing('fair')会生成图 22-3 中的图表。

图 22-3 房价的不同视角

手指练习：相对基线绘制图形有时会启发思考，如图 22-4 所示。修改plot_housing以生成这样的图表。基线以下的条形应为红色。提示：使用plt.bar的bottom关键字参数。

图 22-4 相对于$200,000 的房价

对数 y 轴为制作具有误导性的图表提供了一个绝佳的工具。考虑图 22-5 中的条形图。左侧的图表更准确地展示了跟随 khemric 和 katyperry 的人数差异。右侧图表中跟随人数稀少的 jguttag 使 y 轴不得不将更大比例的长度分配给较小的值，从而留下更少的距离来区分 khemric 和 katyperry 的粉丝数量。¹⁶⁶

图 22-5 比较 Instagram 粉丝数量

22.4 因果关系与相关关系¹⁶⁷

研究表明，定期上课的大学生平均成绩高于偶尔上课的学生。我们这些教授这些课程的人希望相信，这是因为学生从我们教授的课程中学到了东西。当然，学生们能取得更好成绩也可能同样因为那些更有可能上课的学生也更有可能努力学习。

相关性是衡量两个变量朝同一方向移动程度的指标。如果x与y朝同一方向移动，变量是正相关的；如果朝相反方向移动，则是负相关的。如果没有关系，则相关性为0。人们的身高与父母的身高正相关。吸烟与寿命之间的相关性是负的。

当两件事相关时，人们容易假设其中一件事导致了另一件事。考虑北美的流感发生率。病例数呈可预测的模式上下波动。夏季几乎没有病例；病例数在初秋开始上升，接着在夏季来临时开始下降。现在考虑上学的儿童人数。夏季上学的儿童非常少；入学人数在初秋开始上升，然后在夏季来临时下降。

学校开学与流感发生率之间的关联是无可争辩的。这使得一些人得出结论，上学是流感传播的重要因果因素。这可能是对的，但我们不能仅基于相关性得出这一结论。相关性并不意味着因果关系！毕竟，相关性同样可以用来证明流感爆发导致学校开学的信念。或者，也许在任何一个方向上都没有因果关系，而是我们尚未考虑的某个潜在变量导致了两者的发生。实际上，流感病毒在凉爽干燥的空气中存活的时间明显长于温暖潮湿的空气，在北美，流感季节与学校开学时间都与较凉爽和干燥的天气相关。

只要有足够的回溯数据，总是可以找到两个相关的变量，正如图 22-6 中的图表所示。¹⁶⁸

图 22-6 墨西哥柠檬能救命吗？

当发现这样的相关性时，首先要问的是是否有一个合理的理论来解释这种相关性。

陷入cum hoc ergo propter hoc谬论可能是相当危险的。2002 年初，大约有六百万美国女性在被开处方激素替代疗法（HRT），她们相信这将大大降低她们的心血管疾病风险。这一信念得到了几项高度可信的已发表研究的支持，这些研究表明，使用 HRT 的女性心血管死亡的发生率降低。

当《美国医学协会杂志》发表一篇文章声称激素替代疗法（HRT）实际上增加了心血管疾病的风险时，许多女性及其医生感到十分惊讶。¹⁶⁹ 这怎么可能发生？

对一些早期研究的重新分析表明，进行 HRT 的女性往往来自饮食和锻炼习惯优于平均水平的群体。或许进行 HRT 的女性在健康意识上普遍高于研究中的其他女性，因此 HRT 和心脏健康改善是共同原因的巧合效应。

指尖练习：在过去的 100 年里，加拿大每年的死亡人数与每年肉类消费量呈正相关。有什么潜在变量可以解释这一点？

22.5 统计指标并不能全面反映事实

从数据集中可以提取出大量不同的统计数据。通过仔细选择这些数据，可以传达关于同一数据的不同印象。一个好的解毒剂是查看数据集本身。

1973 年，统计学家 F.J.安斯科姆发表了一篇论文，表格见图 22-7，通常称为安斯科姆四重奏。它包含了来自四个数据集的点的<x, y>坐标。这四个数据集的x（9.0）和y（7.5）的均值相同，x（10.0）和y（3.75）的方差相同，以及x与y之间的相关性（0.816）相同。如果我们使用线性回归为每个数据集拟合一条线，结果都是y = 0.5x + 3。

图 22-7 安斯科姆四重奏的统计数据

这是否意味着没有明显的方法可以区分这些数据集？不。我们只需绘制数据，就会发现这些数据集并不相同（图 22-8）。

图 22-8 安斯科姆四重奏的数据

道理很简单：如果可能的话，始终查看一些原始数据的代表。

22.6 抽样偏差

在第二次世界大战期间，每当一架盟军飞机从欧洲的任务中返回时，都会检查飞机，看高射炮弹是在哪些地方击中的。根据这些数据，机械师会加强那些看起来最有可能被击中的飞机部位。

这有什么问题呢？他们没有检查那些未能从任务中返回的飞机，因为它们是被高射炮击落的。也许这些未检查的飞机正是因为在高射炮造成最大伤害的地方受到了攻击而未能返回。这种特定的错误被称为非响应偏差。在调查中相当常见。例如，在许多大学，学生会在学期末的一次讲座中被要求填写一份表格，评价教授讲座的质量。尽管这些调查的结果常常不尽如人意，但情况可能更糟。那些认为讲座糟糕到不值得参加的学生并没有被纳入调查中。¹⁷⁰

如第十九章所讨论，所有统计技术都是基于这样一个假设：通过抽样某一人群的子集，我们可以推断出整个群体的情况。如果使用随机抽样，我们可以对样本与整个群体之间的期望关系做出精确的数学陈述。不幸的是，许多研究，特别是在社会科学领域，是基于所谓的便利（或偶然）抽样。这涉及到根据样本获取的便利性选择样本。为什么那么多心理学研究使用本科生作为研究对象？因为他们在大学校园中容易找到。便利样本可能是具有代表性的，但我们无法确定它实际上是否具有代表性。

手指练习：一种疾病的感染致死率是感染该疾病的人数与因该疾病死亡的人数之比。疾病的病例致死率是被诊断为该疾病的人数与因该疾病死亡的人数之比。哪一种更容易准确估计，为什么？

22.7 上下文很重要

在阅读数据时，尤其是将数据放在上下文之外时，很容易对数据的含义解读得过于深刻。2009 年 4 月 29 日，CNN 报道说：“墨西哥卫生官员怀疑，猪流感疫情已导致超过159人死亡和大约2,500人感染。”听起来相当可怕——直到我们将其与每年在美国季节性流感中约36,000的死亡人数进行比较。

一个常被引用且准确的统计数据是大多数汽车事故发生在离家10英里以内。那么这又意味着什么呢？大多数驾驶都是在离家10英里以内完成的！此外，在这个上下文中，“家”意味着什么？该统计数据是以注册汽车的地址作为“家”来计算的。你是否可以通过将汽车注册在某个遥远的地方来降低发生事故的概率？

反对美国政府减少枪支普及率的倡议者喜欢引用这样的统计数据：大约99.8%的美国枪支在任何给定年份都不会用于暴力犯罪。但没有一些背景，很难知道这意味着什么。这是否意味着美国的枪支暴力并不严重？全国步枪协会报告称，美国大约有3亿支私人拥有的枪支——0.2%的3亿是600,000！

22.8 比较苹果与橙子

快速查看图 22-9 中的图像。

图 22-9 福利与全职工作

这给你留下了什么印象？领取福利的美国人是否比工作的人要多得多？

左边的柱子比右边的柱子高约 500%。然而，柱子上的数字告诉我们 y 轴已经被截断。如果没有被截断，左边的柱子只会高出 6.8%。不过，想到领取福利的人比工作的人多 6.8%还是让人震惊的，震惊且误导。

“领取福利的人数”来源于美国人口普查局对参与条件性援助项目人数的统计。这项统计包括任何居住在至少一人领取福利的家庭中。例如，考虑一个由两个父母和三个孩子组成的家庭，其中一个父母有全职工作，另一个有兼职工作。如果该家庭领取了食品券，则该家庭将为“领取福利的人数”统计增加五人，为全职工作统计增加一人。

这两个数字都是“正确”的，但它们不可比较。这就像是得出结论，奥尔加比马克更会种田，因为她每英亩种植 20 吨土豆，而马克每英亩只种 3 吨蓝莓。

22.9 摘樱桃

既然我们谈到了水果，摘樱桃和比较苹果与橙子一样糟糕。选择性取材涉及选择特定的数据，而忽略其他数据，以支持某种立场。

请考虑图 22-10 中的图表。趋势很明显，但如果我们希望用这些数据来争辩地球没有变暖，我们可以引用 2013 年 4 月冰量比 1988 年 4 月更多这一事实，而忽略其他数据。

图 22-10 北极海冰

22.10 小心推断

从数据中推断是非常容易的。我们在第 20.1.1 节中就这样做了，当时我们将从线性回归中得出的拟合延伸到了回归中使用的数据之外。只有在有合理的理论依据时，才应进行外推。尤其要警惕直线外推。

考虑一下图 22-11 左侧的图表。它显示了 1994 年至 2000 年间美国互联网使用的增长。正如你所看到的，直线提供了相当不错的拟合。

图 22-11 美国互联网使用的增长。

图 22-11 右侧的图表利用这一拟合预测了未来几年使用互联网的美国人口百分比。这个预测让人难以相信。似乎到 2009 年，美国的每个人都在使用互联网是非常不太可能的，而到 2015 年，美国使用互联网的人口超过140%更是不可能。

22.11 德克萨斯神枪手谬论。

想象一下，你在德克萨斯的乡村公路上行驶。你看到一座有六个靶子画在上面的谷仓，每个靶子的正中心都有一个子弹孔。“是的，先生，”谷仓的主人说，“我从不失手。” “没错，”他的配偶说，“没有一个德克萨斯州的人比他用刷子更准确。”明白了吗？他开了六枪，然后在他们周围画了靶子。

图 22-12 教授困惑于学生们扔粉笔的准确性。

这一类型的经典案例出现在 2001 年。¹⁷¹它报告了一项研究团队在阿伯丁的皇家康希尔医院发现，“厌食症女性最有可能是在春季或初夏出生……在三月到六月之间，出生的厌食症患者比平均多了13%，而六月份出生的则多了30%。”

让我们看看那些在六月份出生的女性的令人担忧的统计数据。该团队研究了446名被诊断为厌食症的女性，因此每月的平均出生人数略多于37。这表明六月出生的人数为48 (37*1.3)。让我们写一个短程序(图 22-13)来估计这纯属偶然发生的概率。

图 22-13 在六月出生 48 名厌食症患者的概率。

当我们运行june_prob(10000)时，它打印了。

`Probability of at least 48 births in June = 0.0427`

看起来在六月份出生至少有48个婴儿的概率大约是4.25%。所以也许阿伯丁的那些研究者发现了一些东西。嗯，如果他们从假设更多将成为厌食症患者的婴儿在六月份出生开始，然后进行一项旨在检验该假设的研究，那他们可能确实发现了一些东西。

但他们并没有这样做。相反，他们查看了数据，然后模仿德克萨斯神枪手，在六月周围画了一个圈。正确的统计问题是，在12个月中，至少有48名婴儿出生的概率是什么。图 22-14 中的程序回答了这个问题。

图 22-14 48 名厌食症患者出生于某个月的概率

调用any_prob(10000)打印

`Probability of at least 48 births in some month = 0.4357` 

看起来，研究中报告的结果并非如此不太可能，这些结果反映了一个偶然事件，而不是出生月份与厌食症之间的真实关联。人们不必来自德克萨斯州就能成为德克萨斯神枪手谬论的受害者。

结果的统计显著性取决于实验的进行方式。如果阿伯丁团队最初假设在六月出生的厌食症患者更多，那么他们的结果是值得考虑的。但如果他们假设存在一个异常比例的厌食症患者出生的月份，那么他们的结果就不是很有说服力。实际上，他们是在测试多个假设并选择性地挑选结果。他们可能应该应用 Bonferroni 校正（见第 21.6 节）。

阿伯丁团队可能采取了什么后续步骤来验证他们新发现的假设？一种可能性是进行前瞻性研究。在前瞻性研究中，研究者从一组假设开始，在研究对象发展出感兴趣的结果（在本例中为厌食症）之前招募对象，然后在一段时间内跟踪这些对象。如果该小组进行了具有特定假设的前瞻性研究，并得到了类似的结果，我们可能会信服。

前瞻性研究可能既昂贵又耗时。在回顾性研究中，必须以减少获取误导性结果的方式分析现有数据。一种常见的技术，如第 20.4 节所述，是将数据分成训练集和保留测试集。例如，他们可以随机选择446/2名女性作为训练集，并统计每个月的出生人数。然后，他们可以将这与剩余女性（保留集）每个月的出生人数进行比较。

22.12 百分比可能会引起混淆

一位投资顾问打电话给客户，报告他的股票投资组合在上个月上涨了16%。顾问承认，过去一年中有一些起伏，但很高兴地报告平均每月变动为+0.5%。想象一下，当客户收到一年的报表并注意到他的投资组合价值在一年内下降时的惊讶。

他打电话给他的顾问，并指责他是个骗子。“在我看来，”他说，“我的投资组合下降了大约8%，而你告诉我它上涨了0.5%每月。” “我没有，”财务顾问回答。“我告诉你平均每月变化是+0.5%。”当他检查每月的账单时，投资者意识到他并没有被欺骗，只是被误导了。他的投资组合在上半年每个月下降了15%，然后在下半年每个月上涨了16%。

在考虑百分比时，我们始终需要关注计算百分比的基数。在这种情况下，15%的下降是基于一个比16%的上升更高的平均基础。

当应用于小基数时，百分比可能特别具有误导性。你可能会看到关于某种药物的报道，称其副作用使某种疾病的发生率增加了200%。但如果该疾病的基线发生率非常低，比如1,000,000中有一个，你可能会认为服用该药物的风险被其积极效果所抵消。

手指练习：2020 年 5 月 19 日，纽约时报报道称美国航空旅行在一个月内增加了 123%（从 95,161 名乘客增至 212,508 名乘客）。它还报道了这一增长是在最近一次航空旅行下降 96%之后发生的。那么总的净百分比变化是多少？

22.13 回归谬误

回归谬误发生在当人们未能考虑事件的自然波动时。

所有运动员都有好日子和坏日子。当他们有好日子时，他们会尽量不做任何改变。然而，当他们经历一系列异常糟糕的日子时，他们通常会尝试进行改变。即使这些改变实际上并没有帮助，均值回归（第 17.3 节）使得在接下来的几天里，运动员的表现可能会比之前异常糟糕的表现更好。这可能会误导运动员，以为有治疗效果，即将改善的表现归因于他或她所做的改变。

诺贝尔奖得主心理学家丹尼尔·卡尼曼讲述了一个故事，关于一位以色列空军飞行教官拒绝了卡尼曼的主张，即“对改善表现的奖励比对错误的惩罚更有效。”教官的论点是：“我曾多次表扬飞行学员清晰执行某些特技动作。下次他们尝试相同动作时，通常会表现得更糟。另一方面，我常常对学员的错误执行大声训斥，通常他下次会表现得更好。” ¹⁷² 人类自然会想象治疗效果，因为我们喜欢因果推理。但有时这仅仅是运气问题。

想象一个不存在的治疗效果可能是危险的。这可能导致人们相信疫苗有害，蛇油能治愈所有疼痛，或仅投资于“去年超过市场”的基金是一个好策略。

22.14 统计显著差异可能是微不足道的

在毛伊岛科技学院（MIT），一位招生官希望向世界证明 MIT 的招生过程是“性别盲”的，宣称：“在 MIT，男性和女性的平均成绩没有显著差异。”同一天，一位热情的女性沙文主义者宣称“在 MIT，女性的平均成绩显著高于男性。”一位困惑的学生报社记者决定调查数据以揭露这个说谎者。但当她最终设法从大学获得数据时，她得出结论：两者都在说真话。

句子“在麻省理工学院，女性的平均成绩显著高于男性”实际上意味着什么？没有学习过统计学的人（大多数人群）可能会得出女性和男性的 GPA 之间存在“有意义”的差异。相比之下，最近学习过统计的人可能只会得出 1）女性的平均 GPA 高于男性的平均 GPA，2）可以在5%水平上拒绝因随机性导致的 GPA 差异的零假设。

假设，例如，有2,500名女性和2,500名男性在麻省理工学院学习。进一步假设男性的平均 GPA 为3.5，女性的平均 GPA 为3.51，男性和女性的 GPA 标准差均为0.25。大多数理智的人会认为 GPA 的差异“微不足道”。然而，从统计学的角度看，这一差异在接近2%水平上是“显著的”。这种奇怪二分法的根源是什么？正如我们在第 21.5 节中所示，当研究具有足够的效能——即足够的样本时，即使微不足道的差异也可以统计显著。

当研究规模非常小的时候，相关问题就会出现。假设你投了两次硬币，结果都是正面。现在，让我们使用在第 21.3 节中看到的双尾单样本 t 检验来检验硬币是公正的零假设。如果我们假设正面的值为1，反面的值为0，我们可以使用代码得到 p 值。

`scipy.stats.ttest_1samp([1, 1], 0.5)[1]`

它返回一个 p 值为0，这表明如果硬币是公正的，那么得到两个连续正面的概率为零。如果我们采用贝叶斯方法，从硬币是公正的先验出发，我们会得到不同的答案。

22.15 只需小心

填充几百页统计滥用的历史会很简单，也很有趣。但现在你可能已经明白：用数字撒谎和用文字撒谎一样简单。在你得出结论之前，确保你了解实际测量的内容以及那些“统计显著”结果是如何计算的。正如诺贝尔经济学奖获得者罗纳德·科斯所说：“如果你足够折磨数据，它会承认任何事情。”

22.16 本章引入的术语

• 垃圾进垃圾出（GIGO）

• 独立性假设

• 条形图

• 相关性

• 因果关系

• 潜在变量

• 非应答偏差

• 便利性（偶然）抽样

• 感染致死率

• 病例致死率

• 精挑细选

• 前瞻性研究

• 回顾性研究

• 回归谬误

• 治疗效果

第二十三章：使用 PANDAS 探索数据

本书后半部分的大部分内容都集中在构建各种计算模型上，这些模型可以用来从数据中提取有用的信息。在接下来的章节中，我们将快速了解使用机器学习从数据构建模型的简单方法。

然而，在这样做之前，我们将查看一个流行的库，它可以快速帮助我们熟悉数据集，然后再深入到更详细的分析中。Pandas¹⁷³ 是建立在 numpy 之上的。Pandas 提供了促进

• 组织数据

• 计算数据的简单统计

• 以便于未来分析的格式存储数据

23.1 DataFrames 和 CSV 文件

Pandas 中的一切都是围绕 DataFrame 类型构建的。DataFrame 是一个可变的二维表格数据结构，具有标记的轴（行和列）。可以将其视为一种增强版的电子表格。

虽然可以使用 Python 代码从头构建 DataFrame，但创建 DataFrame 更常见的方法是通过读取 CSV 文件。正如我们在第十九章中看到的，CSV 文件的每一行由一个或多个值组成，这些值由逗号分隔。¹⁷⁴ CSV 文件通常用于以纯文本格式存储表格数据。在这种情况下，行通常具有相同数量的字段。由于它们是纯文本，因此经常用于将数据从一个应用程序移动到另一个应用程序。例如，大多数电子表格程序允许用户将电子表格的内容写入 CSV 文件。

图 23-1 显示了一个 DataFrame，其中包含关于 2019 年 FIFA 女子世界杯后期轮次的信息。每一列代表一个称为 series 的东西。每一行都有一个 index 相关联。默认情况下，索引是连续的数字，但它们并不一定需要如此。每一列都有一个 name 相关联。正如我们将看到的，这些名称的作用类似于字典中的键。

图 23-1 一个绑定到变量 wwc 的示例 Pandas DataFrame

在 图 23-1 中显示的 DataFrame 是使用下面的代码和 图 23-2 中所示的 CSV 文件生成的。

import pandas as pd
`wwc = pd.read_csv('wwc2019_q-f.csv') print(wwc)`

图 23-2 一个示例 CSV 文件

在导入 Pandas 后，代码使用 Pandas 的函数 read_csv 来读取 CSV 文件，然后以 图 23-1 中所示的表格形式打印出来。如果 DataFrame 具有大量的行或列，print 将在 DataFrame 的中心用省略号替换列和/或行。这可以通过首先使用 DataFrame 方法 to_string 将 DataFrame 转换为字符串来避免。

一行索引和一个列标签一起指示一个数据单元（如在电子表格中）。我们将在第 23.3 节讨论如何访问单个单元和单元组。通常，但并不总是，列中的单元格都是同一种类型。在 图 23-1 的 DataFrame 中，Round、Winner 和 Loser 列中的每个单元格都是 str 类型。W Goals 和 L Goals 列中的单元格是 numpy.int64 类型。如果你将它们视为 Python 的整型，你就不会有问题。

我们可以直接使用属性 index、columns 和 values 访问 DataFrame 的三个组成部分。

index 属性的类型为 RangeIndex。例如，wwc.index 的值为 RangeIndex(start=0, stop=8, step=1)。因此，代码

for i in wwc.index:
    print(i)

将按升序打印整数 0-7。

columns 属性的类型为 Index。例如，值 wwc.columns 是 Index(['Round', 'Winner', 'W Goals', 'Loser', 'L Goals'], dtype='object')，而代码

for c in wwc.columns:
    print(c)

打印

Round
Winner
W Goals
Loser
L Goals

values 属性的类型为 numpy.ndarray。在第十三章中，我们介绍了类型 numpy.array。事实证明，array 是 ndarray 的特殊情况。虽然数组是一维的（像其他序列类型一样），但 ndarrays 可以是多维的。ndarray 的维数和项目数称为其 shape，并用一个表示每个维度大小的非负整数元组表示。值 wwc.values 是二维 ndarrays

[['Quarters' 'England' 3 'Norway' 0]
 ['Quarters' 'USA' 2 'France' 1]
 ['Quarters' 'Netherlands' 2 'Italy' 0]
 ['Quarters' ‘Sweden' 2 'Germany' 1]
 ['Semis' 'USA' 2 'England' 1]
 ['Semis' 'Netherlands' 1 ‘Sweden' 0]
 ['3rd Place' ‘Sweden' 2 'England' 1]
 ['Championship' 'USA' 2 'Netherlands' 0]]

由于它有八行五列，因此它的形状为 (8, 5)。

23.2 创建系列和 DataFrame

实际上，Pandas 的 DataFrame 通常是通过加载存储为 SQL 数据库、CSV 文件或与电子表格应用程序相关的格式的数据集来创建的。然而，有时使用 Python 代码构建系列和 DataFrame 是有用的。

表达式 pd.DataFrame() 生成一个空的 DataFrame，而语句 print(pd.DataFrame()) 产生的输出

Empty DataFrame
Columns: []
Index: []

创建一个非空 DataFrame 的简单方法是传入一个列表。例如，代码

rounds = ['Semis', ‘Semis', '3rd Place', 'Championship']
print(pd.DataFrame(rounds))

打印

              0
0         Semis
1         Semis
2     3rd Place
3  Championship

请注意，Pandas 为 DataFrame 的唯一一列自动生成了一个标签，尽管这个标签并不是特别描述性。为了获得更具描述性的标签，我们可以传入一个字典而不是一个列表。例如，代码 print(pd.DataFrame({'Round': rounds})) 会打印

          Round
0         Semis
1         Semis
2     3rd Place
3  Championship

要直接创建一个具有多列的 DataFrame，我们只需传入一个包含多个条目的字典，每个条目由一个列标签作为键和一个与每个键相关联的列表作为值。这些列表必须具有相同的长度。例如，代码

rounds = ['Semis', ‘Semis', '3rd Place', 'Championship']
teams = ['USA', 'Netherlands', ‘Sweden', 'USA']
df = pd.DataFrame({'Round': rounds, 'Winner': teams})
print(df)

打印

          Round       Winner
0         Semis          USA
1         Semis  Netherlands
2     3rd Place       Sweden
3  Championship          USA

一旦创建了 DataFrame，就很容易添加列。例如，语句 df['W Goals'] = [2, 1, 0, 0] 会修改 df，使其值变为

          Round       Winner  W Goals
0         Semis          USA        2
1         Semis  Netherlands        1
2     3rd Place       Sweden        0
3  Championship          USA        0

就像字典中与键关联的值可以被替换一样，与列关联的值也可以被替换。例如，在执行语句df['W Goals'] = [2, 1, 2, 2]后，df的值变为

          Round       Winner  W Goals
0         Semis          USA        2
1         Semis  Netherlands        1
2     3rd Place       Sweden        2
3  Championship          USA        2

从 DataFrame 中删除列也是很简单的。函数调用print(df.drop('Winner', axis = 'columns'))会打印出

          Round  W Goals
0         Semis        2
1         Semis        1
2     3rd Place        2
3  Championship        2

并且df保持不变。如果在调用drop时没有包含axis = 'columns'（或等效地axis = 1），则轴将默认为'rows'（等效于axis = 0），这将导致生成异常KeyError: "['Winner'] not found in axis."

如果一个 DataFrame 很大，以这种方式使用drop是低效的，因为这需要复制 DataFrame。通过将drop的inplace关键字参数设置为True，可以避免复制。调用df.drop('Winner', axis = 'columns', inplace = True)会修改df并返回None。

可以使用DataFrame构造函数将行添加到 DataFrame 的开头或末尾，然后使用concat函数将新 DataFrame 与现有 DataFrame 组合。例如，代码

quarters_dict = {'Round': ['Quarters']*4,
                 'Winner': ['England', 'USA', 'Netherlands', ‘Sweden'],
                 'W Goals': [3, 2, 2, 2]}
df = pd.concat([pd.DataFrame(quarters_dict), df], sort = False)

将df设置为

          Round       Winner  W Goals
0      Quarters      England        3
1      Quarters          USA        2
2      Quarters  Netherlands        2
3      Quarters       Sweden        2
0         Semis          USA        2
1         Semis  Netherlands        1
2     3rd Place       Sweden        2
3  Championship          USA        2

如果将关键字参数sort设置为True，concat也会根据列标签的字典序改变列的顺序。也就是说，

pd.concat([pd.DataFrame(quarters_dict), df], sort = True)

交换最后两列的位置并返回 DataFrame。

          Round  W Goals       Winner
0      Quarters        3      England
1      Quarters        2          USA
2      Quarters        2  Netherlands
3      Quarters        2       Sweden
0         Semis        2          USA
1         Semis        1  Netherlands
2     3rd Place        2       Sweden
3  Championship        2          USA

如果没有提供sort的值，则默认为False。

注意每个连接的 DataFrame 的索引保持不变。因此，会有多个具有相同索引的行。可以使用reset_index方法重置索引。例如，表达式df.reset_index(drop = True)的值为

             Round       Winner  W Goals
0      Quarters      England        3
1      Quarters          USA        2
2      Quarters  Netherlands        2
3      Quarters       Sweden        2
4         Semis          USA        2
5         Semis  Netherlands        1
6     3rd Place       Sweden        2
7  Championship          USA        2

如果reset_index被调用时drop = False，则会向 DataFrame 添加一个包含旧索引的新列。该列被标记为index。

你可能会想知道为什么 Pandas 甚至允许重复索引。原因是使用语义上有意义的索引来标记行通常是有帮助的。例如，df.set_index('Round')的值为

                   Winner  W Goals
Round                             
Quarters          England        3
Quarters              USA        2
Quarters      Netherlands        2
Quarters           Sweden        2
Semis                 USA        2
Semis         Netherlands        1
3rd Place          Sweden        2
Championship          USA        2

23.3 选择列和行

对于 Python 中的其他复合类型，方括号是选择 DataFrame 部分的主要机制。要选择 DataFrame 的单列，只需将列的标签放在方括号之间。例如，wwc['Winner']的值为

0        England
1            USA
2    Netherlands
3         Sweden
4            USA
5    Netherlands
6         Sweden
7            USA

该对象的类型是**Series**，即它不是 DataFrame。Series 是一维值序列，每个值都有一个索引标签。要从 Series 中选择单个项目，我们在系列后面的方括号中放入一个索引。因此，wwc['Winner'][3]的值为字符串Sweden。

我们可以使用for循环遍历一个序列。例如，

winners = ''
for w in wwc['Winner']:
    winners += w + ','
print(winners[:-1])

打印出England,USA,Netherlands,Sweden,USA,Netherlands,Sweden,USA。

指尖练习： 编写一个函数，返回获胜者所进球的总和。

方括号也可以用来从 DataFrame 中选择多个列。这是通过在方括号内放置列标签的列表来完成的。这将产生一个 DataFrame 而不是 Series。例如，wwc[['Winner', 'Loser']]产生的 DataFrame 为：

        Winner        Loser
0      England       Norway
1          USA       France
2  Netherlands        Italy
3       Sweden      Germany
4          USA      England
5  Netherlands       Sweden
6       Sweden      England
7          USA  Netherlands

选择方括号中的标签列表的顺序不必与原始 DataFrame 中的标签顺序相同。这使得通过选择来重新组织 DataFrame 变得方便。例如，wwc[['Round','Winner','Loser','W Goals','L Goals']]返回的 DataFrame 为：

          Round       Winner        Loser  W Goals  L Goals
0      Quarters      England       Norway        3        0
1      Quarters          USA       France        2        1
2      Quarters  Netherlands        Italy        2        0
3      Quarters       Sweden      Germany        2        1
4         Semis          USA      England        2        1
5         Semis  Netherlands       Sweden        1        0
6     3rd Place       Sweden      England        2        1
7  Championship          USA  Netherlands        2        0

请注意，尝试通过将行的索引放入方括号中来选择一行将不起作用。这将产生一个KeyError异常。然而，奇怪的是，我们可以通过切片选择行。因此，虽然wwc[1]会导致异常，wwc[1:2]会产生一个包含单行的 DataFrame。

 Round Winner  W Goals   Loser  L Goals
 1  Quarters    USA        2  France        1

我们将在下一小节讨论其他选择行的方法。

23.3.1 使用 loc 和 iloc 进行选择

**loc**方法可以用于从 DataFrame 中选择行、列或行列的组合。重要的是，所有选择都是通过标签进行的。这一点值得强调，因为某些标签（例如索引）看起来可能像数字。

如果df是一个 DataFrame，则表达式df.loc[label]返回与df中label关联的行对应的 Series。例如，wwc.loc[3]返回的 Series 为：

Round      Quarters
Winner       Sweden
W Goals           2
Loser       Germany
L Goals           1

请注意，wwc的列标签是 Series 的索引标签，与这些标签相关联的值是wwc中标签为3的行对应列的值。

要选择多行，我们只需在.loc后面的方括号内放置标签列表（而不是单个标签）。这样做时，表达式的值是一个 DataFrame 而不是一个 Series。例如，表达式wwc.loc[[1,3,5]]产生：

      Round       Winner  W Goals    Loser  L Goals
1  Quarters          USA        2   France        1
3  Quarters       Sweden        2  Germany        1
5     Semis  Netherlands        1   Sweden        0

请注意，新 DataFrame 中每一行的索引是旧 DataFrame 中该行的索引。

切片提供了另一种选择多行的方法。一般形式为df.loc[first:last:step]。如果未提供first，则默认为 DataFrame 中的第一个索引。如果未提供last，则默认为 DataFrame 中的最后一个索引。如果未提供step，则默认为1。表达式wwc.loc[3:7:2]产生的 DataFrame 为：

          Round       Winner  W Goals        Loser  L Goals
3      Quarters       Sweden        2      Germany        1
5         Semis  Netherlands        1       Sweden        0
7  Championship          USA        2  Netherlands        0

作为一名 Python 程序员，你可能会惊讶于标签为7的行被包含在内。对于其他 Python 数据容器（如列表），切片时最后一个值会被排除，但对于 DataFrame 则不是这样。¹⁷⁵ 表达式wwc.loc[6:]产生的 DataFrame 为：

          Round  Winner  W Goals        Loser  L Goals
6     3rd Place  Sweden        2      England        1
7  Championship     USA        2  Netherlands        0

表达式wwc.loc[:2]产生：

      Round       Winner  W Goals   Loser  L Goals
0  Quarters      England        3  Norway        0
1  Quarters          USA        2  France        1
2  Quarters  Netherlands        2   Italy        0

指尖练习： 写一个表达式，选择wwc中所有偶数编号的行。

如前所述，loc可以用来同时选择行和列的组合。这是通过类似以下形式的表达式完成的：

df.loc[*row_selector*, *column_selector*]

行和列选择器可以使用之前讨论过的任何机制编写，即单个标签、标签列表或切片表达式。例如，wwc.loc[0:2, 'Round':'L Goals':2]生成

 Round  W Goals  L Goals 
0  Quarters        3        0 
1  Quarters        2        1 
2  Quarters        2        0

手指练习： 写一个生成数据框的表达式

 Round       Winner  W Goals   Loser  L Goals
1  Quarters          USA        2  France        1
2  Quarters  Netherlands        2   Italy        0

到目前为止，如果你把索引标签视为整数，你是不会错的。让我们看看当 1)标签不是数字型，2)有多行具有相同标签时，选择是如何工作的。让wwc_by_round成为数据框

                   Winner  W Goals        Loser  L Goals
Round                                                   
Quarters          England        3       Norway        0
Quarters              USA        2       France        1
Quarters      Netherlands        2        Italy        0
Quarters           Sweden        2      Germany        1
Semis                 USA        2      England        1
Semis         Netherlands        1       Sweden        0
3rd Place          Sweden        2      England        1
Championship          USA        2  Netherlands        0

你认为表达式wwc_by_round.loc['Semis']的结果是什么？它选择所有标签为Semis的行，以返回

            Winner  W Goals    Loser  L Goals
Round                                        
Semis          USA        2  England        1
Semis  Netherlands        1   Sweden        0

同样，wwc_by_round.loc[['Semis', 'Championship']]选择所有标签为Semis或Championship的行：

                   Winner  W Goals        Loser  L Goals
Round                                                   
Semis                 USA        2      England        1
Semis         Netherlands        1       Sweden        0
Championship          USA        2  Netherlands        0

切片也适用于非数字索引。表达式

`wwc_by_round.loc['Quarters':'Semis':2]`

通过选择标记为Quarters的第一行，然后选择每隔一行，直到经过标记为Semis的行，生成

               Winner  W Goals    Loser  L Goals
Round                                           
Quarters      England        3   Norway        0
Quarters  Netherlands        2    Italy        0
Semis             USA        2  England        1

现在，假设我们想选择标记为Quarters的第二行和第三行。我们不能简单地写wwc_by_round.loc['Quarters']，因为那样会选择所有四行标记为Quarters。使用iloc方法。

**iloc**方法类似于loc，但它是基于整数而不是标签（因此有i在iloc中）。数据框的第一行是iloc 0，第二行为iloc 1，依此类推。因此，要选择标记为Quarters的第二行和第三行，我们写wwc_by_round.iloc[[1,2]]。

23.3.2 按组选择

将数据框分成子集，并对每个子集分别应用一些聚合或转换，通常很方便。groupby方法使得这种操作变得简单。

假设，例如，我们想知道每轮中获胜和失利球队总进球数。代码

grouped_by_round = wwc.groupby('Round')

将group_by_round绑定到类型为DataFrameGroupBy的对象。我们可以对该对象应用聚合器sum以生成一个数据框。代码

grouped_by_round = wwc.groupby('Round')
print(grouped_by_round.sum())

打印

               W Goals  L Goals
Round                         
3rd Place           2        1
Championship        2        0
Quarters            9        2
Semis               3        1

代码print(wwc.groupby('Winner').mean())打印出

             W Goals   L Goals
Winner                        
England          3.0  0.000000
Netherlands      1.5  0.000000
Sweden           2.0  1.000000
USA              2.0  0.666667

从中我们可以很容易地看出，英格兰在赢得的比赛中平均进了三球，同时零封了对手。

代码print(wwc.groupby(['Loser', 'Round']).mean())打印出

`                          W Goals  L Goals Loser       Round                          England     3rd Place           2        1             Semis               2        1 France      Quarters            2        1 Germany     Quarters            2        1 Italy       Quarters            2        0 Netherlands Championship        2        0 Norway      Quarters            3        0 Sweden      Semis               1        0`

从中我们可以很容易地看出，英格兰在输掉的比赛中平均进了一球，而丢了两个。

23.3.3 按内容选择

假设我们想从数据框中选择瑞典赢得的所有比赛的行，如图 23-1。由于这个数据框比较小，我们可以查看每一行并找到对应比赛的行索引。当然，这种方法不适用于大型数据框。幸运的是，使用称为布尔索引的东西可以轻松根据内容选择行。

基本思想是编写一个逻辑表达式，引用 DataFrame 中的值。该表达式随后在 DataFrame 的每一行上进行评估，评估结果为 True 的行将被选中。表达式 wwc.loc[wwc['Winner'] == 'Sweden'] 评估为 DataFrame

       Round  Winner  W Goals    Loser  L Goals
3   Quarters  Sweden        2  Germany        1
6  3rd Place  Sweden        2  England        1

提取所有涉及瑞典的比赛稍微复杂一些。逻辑运算符 **&**（对应于与），**|**（对应于或）和 **–**（对应于非）可用于形成表达式。表达式 wwc.loc[(wwc['Winner'] == 'Sweden') | (wwc['Loser'] == 'Sweden')] 返回

 Round       Winner  W Goals    Loser  L Goals
3   Quarters       Sweden        2  Germany        1
5      Semis  Netherlands        1   Sweden        0
6  3rd Place       Sweden        2  England        1

请注意，逻辑表达式两个子项周围的括号是必要的，因为在 Pandas 中，| 的优先级高于 ==。

练习: 编写一个表达式，返回一个包含美国参赛但法国没有参赛的比赛的 DataFrame。

如果我们预计会进行许多查询以选择某个国家参与的比赛，定义一个函数可能会很方便。

   def get_country(df, country):
    """df a DataFrame with series labeled Winner and Loser
       country a str
       returns a DataFrame with all rows in which country appears
       in either the Winner or Loser column"""
    return df.loc[(df['Winner'] == country) | (df['Loser'] == country)]

由于 get_country 返回一个 DataFrame，因此通过组合两次调用 get_country 很容易提取成对球队之间的比赛。例如，评估 get_country(get_country(wwc, 'Sweden'),'Germany') 提取了这两支球队之间的一场比赛（球队在淘汰赛阶段最多相互对阵一次）。

假设我们想要对 get_country 进行概括，使其接受国家列表作为参数并返回列表中任何国家参加的所有比赛。我们可以使用 isin 方法来实现：

  def get_games(df, countries):
    return df[(df['Winner'].isin(countries)) |
            (df['Loser'].isin(countries))]

isin 方法通过在指定列中选择仅包含指定值（或指定值集合中的元素）的行来过滤 DataFrame。在 get_games 的实现中，表达式 df['Winner'].isin(countries) 选择了 df 中 Winner 列包含 countries 列表中的元素的行。

练习: 打印一个只包含瑞典与德国或荷兰比赛的 DataFrame。

23.4 在 DataFrame 中操作数据

我们现在已经查看了一些创建和选择 DataFrame 部分的简单方法。创建 DataFrame 的一个原因是能够轻松提取聚合信息。让我们首先看看如何从 DataFrame wwc 中提取聚合信息，如 图 23-1 所示。

DataFrame 的列可以以类似于我们对 numpy 数组操作的方式进行操作。例如，类似于表达式 2*np.array([1,2,3]) 评估为数组 [2 4 6]，表达式 2*wwc['W Goals'] 评估为系列

0    6
1    4
2    4
3    4
4    4
5    2
6    4
7    4

表达式 wwc['W Goals'].sum() 对 W Goals 列中的值进行求和，得到值 16。类似地，表达式

(wwc[wwc['Winner'] == ‘Sweden']['W Goals'].sum() +
wwc[wwc['Winner'] == ‘Sweden']['L Goals'].sum())

计算瑞典队进球总数为 6，并且表达式

`(wwc['W Goals'].sum() - wwc['L Goals'].sum())/len(wwc['W Goals'])`

计算 DataFrame 中比赛的平均进球差为 1.5。

手指练习： 写一个表达式，计算所有轮次中进球的总数。

手指练习： 写一个表达式，计算四分之一决赛中输球队进球的总数。

假设我们想添加一列，包含所有比赛的进球差异，并添加一行，总结所有包含数字的列的总计。添加列很简单。我们只需执行wwc['G Diff'] = wwc['W Goals'] - wwc['L Goals']。添加行则更复杂。我们首先创建一个包含所需行内容的字典，然后使用该字典创建一个只包含新行的新数据框。接着，我们使用concat函数将wwc和新数据框连接起来。

#Add new column to wwc
wwc['G Diff'] = wwc['W Goals'] - wwc['L Goals']
#create a dict with values for new row
new_row_dict = {'Round': ['Total'],
           'W Goals': [wwc['W Goals'].sum()],
           'L Goals': [wwc['L Goals'].sum()],
           'G Diff': [wwc['G Diff'].sum()]}
#Create DataFrame from dict, then pass it to concat
new_row = pd.DataFrame(new_row_dict)
wwc = pd.concat([wwc, new_row], sort = False).reset_index(drop = True)

这段代码生成了数据框。

          Round       Winner  W Goals        Loser  L Goals  G Diff
0      Quarters      England        3       Norway        0       3
1      Quarters          USA        2       France        1       1
2      Quarters  Netherlands        2        Italy        0       2
3      Quarters       Sweden        2      Germany        1       1
4         Semis          USA        2      England        1       1
5         Semis  Netherlands        1       Sweden        0       1
6     3rd Place       Sweden        2      England        1       1
7  Championship          USA        2  Netherlands        0       2
8         Total          NaN       16          NaN        4      12

注意，当我们尝试对不包含数字的列中的值求和时，Pandas 没有产生异常。相反，它提供了特殊值NaN（非数字）。

除了提供简单的算术操作，如求和和均值，Pandas 还提供计算各种有用统计函数的方法。其中最有用的是corr，用于计算两个系列之间的相关性。

相关性是一个介于-1 和 1 之间的数字，提供有关两个数值之间关系的信息。正相关表明一个变量的值增加时，另一个变量的值也增加；负相关表明一个变量的值增加时，另一个变量的值减少。相关性为零表示变量之间没有关系。

最常用的相关性度量是皮尔逊相关性。皮尔逊相关性衡量两个变量之间线性关系的强度和方向。除了皮尔逊相关性外，Pandas 还支持其他两种相关性度量，斯皮尔曼和肯德尔。这三种度量之间存在重要差异（例如，斯皮尔曼对异常值的敏感性低于皮尔逊，但仅对单调关系有用），但讨论何时使用哪一种超出了本书的范围。

要打印W Goals、L Goals和G Diff的皮尔逊成对相关性（并排除包含总计的行），我们只需执行。

print(wwc.loc[wwc['Round'] != ‘Total'].corr(method = 'pearson'))

这产生了。

          W Goals   L Goals    G Diff
W Goals  1.000000  0.000000  0.707107
L Goals  0.000000  1.000000 -0.707107
G Diff   0.707107 -0.707107  1.000000

对角线上的值都是 1，因为每个系列与自身完全正相关。不出所料，进球差异与获胜队的进球数强正相关，而与输球队的进球数强负相关。获胜者和输者进球之间较弱的负相关在职业足球中也有其道理。¹⁷⁶

23.5 扩展示例

在本节中，我们将查看两个数据集，一个包含 21 个美国城市的历史温度数据，另一个包含全球化石燃料使用的历史数据。

23.5.1 温度数据

代码

pd.set_option('display.max_rows', 6)
pd.set_option('display.max_columns', 5)
temperatures = pd.read_csv('US_temperatures.csv')
print(temperatures)

打印

           Date  Albuquerque  ...  St Louis  Tampa
0      19610101        -0.55  ...     -0.55  15.00
1      19610102        -2.50  ...     -0.55  13.60
2      19610103        -2.50  ...      0.30  11.95
        ...          ...  ...       ...    ...
20085  20151229        -2.15  ...      1.40  26.10
20086  20151230        -2.75  ...      0.60  25.55
20087  20151231        -0.75  ...     -0.25  25.55
[20088 rows x 22 columns]

前两行代码设置了默认选项，以限制打印 DataFrame 时显示的行数和列数。这些选项的作用类似于我们用于设置各种绘图默认值的rcParams。函数reset_option可以用来将选项恢复为系统默认值。

这个 DataFrame 以一种便于查看特定日期不同城市天气的方式组织。例如，查询

temperatures.loc[temperatures['Date']==19790812][['New York','Tampa']]

告诉我们在 1979 年 8 月 12 日，纽约的温度为 15°C，坦帕为 25.55°C。

手指练习： 写一个表达式，如果 2000 年 10 月 31 日凤凰城比坦帕温暖则评估为True，否则评估为False。

手指练习： 编写代码提取凤凰城温度为 41.4°C 的日期。¹⁷⁷

不幸的是，查看 21 个城市在 20,088 个日期的数据并不能直接洞察与温度趋势相关的大问题。让我们开始添加提供每日温度汇总信息的列。代码

temperatures['Max T'] = temperatures.max(axis = 'columns')
temperatures['Min T'] = temperatures.min(axis = 'columns')
temperatures['Mean T'] = round(temperatures.mean(axis = 'columns'), 2)
print(temperatures.loc[20000704:20000704])

打印

 Date  Albuquerque  ...  Min T      Mean T

14429  20000704        26.65  ...  15.25  1666747.37

2000 年 7 月 4 日，那 21 个城市的平均温度真的比太阳表面的温度高很多吗？可能不是。更可能的是我们的代码存在 bug。问题在于我们的 DataFrame 将日期编码为数字，而这些数字用于计算每行的平均值。从概念上讲，考虑日期作为温度系列的索引可能更有意义。因此，让我们将 DataFrame 中的日期改为索引。代码

temperatures.set_index('Date', drop = True, inplace = True)
temperatures['Max'] = temperatures.max(axis = 'columns')
temperatures['Min'] = temperatures.min(axis = 'columns')
temperatures['Mean T'] = round(temperatures.mean(axis = 'columns'), 2)
print(temperatures.loc[20000704:20000704])

打印更可信的

          Albuquerque  Baltimore  ...  Min T  Mean T
Date                              ...               
20000704        26.65      25.55  ...  15.25   24.42

顺便提一下，由于Date不再是列标签，我们不得不使用不同的打印语句。我们为什么要用切片选择单行？因为我们想创建一个 DataFrame 而不是一个系列。

我们现在可以开始绘制一些显示各种趋势的图表。例如，

plt.figure(figsize = (14, 3)) #set aspect ratio for figure
plt.plot(list(temperatures['Mean T']))
plt.title('Mean Temp Across 21 US Cities')
plt.xlabel('Days Since 1/1/1961')
plt.ylabel('Degrees C')

生成的图表显示了美国温度的季节性。请注意，在绘制平均温度之前，我们将系列转换为列表。如果直接绘制系列，它将使用系列的索引（代表日期的整数）作为 x 轴。这将产生一种相当奇怪的图，因为 x 轴上的点会奇怪地间隔。例如，1961 年 12 月 30 日和 1961 年 12 月 31 日之间的距离为 1，但 1961 年 12 月 31 日和 1962 年 1 月 1 日之间的距离为 8870（19620,101 – 19611231）。

通过放大几年的数据并使用调用plt.plot(list(temperatures['Mean T'])[0:3*365])生成图表，我们可以更清楚地看到季节模式。

在过去的几十年里，关于地球变暖的共识已经形成。让我们看看这些数据是否与这一共识一致。由于我们正在研究一个关于长期趋势的假设，可能不应关注每日或季节性的温度变化。相反，让我们看一下年度数据。

作为第一步，让我们使用temperatures中的数据构建一个新的数据框，其中行代表年份而不是天。执行此操作的代码包含在图 23-3 和图 23-4 中。大部分工作在函数get_dict中完成， 图 23-3 返回一个字典，将年份映射到一个字典，该字典给出与不同标签相关的该年份的值。get_dict的实现使用iterrows遍历temperatures中的行。该方法返回一个迭代器，每一行返回一个包含索引标签和该行内容的系列对。可以使用列标签选择生成系列的元素。¹⁷⁸

图 23-3 构建一个将年份映射到温度数据的字典

如果test是数据框。

          Max T  Min T  Mean T
Date                          
19611230  24.70 -13.35    3.35
19611231  24.75 -10.25    5.10
19620101  25.55 -10.00    5.70
19620102  25.85  -4.45    6.05

调用get_dict(test, ['Max', 'Min'])将返回字典。

{'1961': {'Max T': [24.7, 24.75], 'Min T': [-13.35, -10.25], 'Mean T': [3.35, 5.1]}, '1962': {'Max T': [25.55, 25.85], 'Min T': [-10.0, -4.45], 'Mean T': [5.7, 6.05]}}

图 23-4 围绕年份构建数据框

在图 23-4 中调用get_dict后的代码构建了一个包含出现在temperatures中的每一年的列表，以及包含这些年的最低、最高和平均温度的附加列表。最后，它使用这些列表构建数据框yearly_temps：

    Year  Min T  Max T  Mean T
0   1961 -17.25  38.05   15.64
1   1962 -21.65  36.95   15.39
2   1963 -24.70  36.10   15.50
..   ...    ...    ...     ...
52  2013 -15.00  40.55   16.66
53  2014 -22.70  40.30   16.85
54  2015 -18.80  40.55   17.54

现在我们已经以便于处理的格式获得了数据，让我们生成一些图表来可视化温度随时间的变化。 图 23-5 中的代码生成了图 23-6 中的图表。

图 23-5 生成与年份相关的温度测量图

图 23-6 年度平均和最低温度

图 23-6 左侧的图表显示了一个不可否认的趋势；¹⁷⁹这 21 个城市的平均温度随着时间的推移而上升。右侧的图表则不那么清晰。极端的年度波动使得很难看出趋势。通过绘制温度的移动平均，可以生成更具启示性的图表。

Pandas 方法rolling用于对系列的多个连续值执行操作。评估表达式yearly_temps['Min T'].rolling(7).mean()会生成一个系列，其中前 6 个值为NaN，对于每个大于 6 的 i，系列中的第 i 个值为yearly_temps['Min'][i-6:i+1]的平均值。将该系列与年份绘制在一起会生成图 图 23-7，这确实暗示了一个趋势。

图 23-7 滚动平均最低温度

虽然可视化两个系列之间的关系可以提供信息，但通常更有用的是以更定量的方式观察这些关系。让我们先看一下年份与七年滚动平均最低、最高和平均温度之间的相关性。在计算相关性之前，我们首先更新yearly_temps中的系列以包含滚动平均值，然后将年份值从字符串转换为整数。代码

num_years = 7
for label in ['Min T', 'Max T', 'Mean T']:
    yearly_temps[label] = yearly_temps[label].rolling(num_years).mean()
yearly_temps['Year'] = yearly_temps['Year'].apply(int)
print(yearly_temps.corr())

打印

            Year     Min T     Max T    Mean T
Year    1.000000  0.713382  0.918975  0.969475
Min T   0.713382  1.000000  0.629268  0.680766
Max T   0.918975  0.629268  1.000000  0.942378
Mean T  0.969475  0.680766  0.942378  1.000000

所有汇总的温度值与年份呈正相关，其中平均温度的相关性最强。这引发了一个问题：年份解释了平均温度滚动平均值方差的多少。以下代码打印决定系数（第 20.2.1 节）。

indices = np.isfinite(yearly_temps['Mean T'])
model = np.polyfit(list(yearly_temps['Year'][indices]),
                  list(yearly_temps['Mean T'][indices]), 1)
print(r_squared(yearly_temps['Mean T'][indices],
              np.polyval(model, yearly_temps['Year'][indices])))

由于Mean系列中的一些值为NaN，我们首先使用函数np.isfinite获取yearly_temps['Mean']中非NaN值的索引。然后我们构建一个线性模型，最后使用r_squared函数（见图 20-13）将模型预测的结果与实际温度进行比较。与年份相关的七年滚动平均温度几乎解释了 94%的方差。

指尖练习： 找到平均年温度而非滚动平均和十年滚动平均的决定系数（r²）。

如果你恰好住在美国或计划前往美国，你可能更有兴趣按城市而不是按年份查看数据。让我们首先生成一个新的 DataFrame，以提供每个城市的汇总数据。为了考虑到我们的美国读者，我们通过对city_temps中的所有值应用转换函数，将所有温度转换为华氏度。倒数第二行添加了一列，显示温度变化的极端程度。执行此代码会生成 图 23-8 中的 DataFrame。¹⁸⁰

temperatures = pd.read_csv('US_temperatures.csv')
temperatures.drop('Date', axis = 'columns', inplace = True)
means = round(temperatures.mean(), 2)
maxes = temperatures.max()
mins = temperatures.min()
city_temps = pd.DataFrame({'Min T':mins, 'Max T':maxes,
                           'Mean T':means})
city_temps = city_temps.apply(lambda x: 1.8*x + 32)
city_temps['Max-Min'] = city_temps['Max T'] - city_temps['Min T']
print(city_temps.sort_values('Mean T', ascending = False).to_string())

图 23-8 部分城市的平均温度

为了可视化城市之间的差异，我们使用代码生成了图 图 23-9

plt.plot(city_temps.sort_values('Max-Min', ascending=False)
        ['Max-Min'], 'o')
plt.figure()
plt.plot(city_temps.sort_values('Max-Min', ascending=False)['Min T'],
        'b∧', label = 'Min T')
plt.plot(city_temps.sort_values('Max-Min', ascending=False)['Max T'],
        'kx', label = 'Max T')
plt.plot(city_temps.sort_values('Max-Min', ascending=False)['Mean T'],
        'ro', label = 'Mean T')
plt.xticks(rotation = ‘vertical')
plt.legend()
plt.title('Variation in Extremal Daily\nTemperature 1961-2015')
plt.ylabel('Degrees F')

请注意，我们对所有三个系列使用了排序顺序Max - Min。使用ascending = False会逆转默认的排序顺序。

图 23-9 温度极端变化

看这个图我们可以看到，除了其他内容之外，

• 在城市之间，最低温度的差异远大于最高温度。因此，最大值减去最小值（排序顺序）与最低温度之间存在强正相关。

• 旧金山或西雅图从未变得非常炎热。

• 圣胡安的温度几乎保持不变。

• 芝加哥的温度并不接近恒定。这个风城的温度既会变得相当炎热，也会变得令人畏惧的寒冷。

• 凤凰城和拉斯维加斯都变得异常炎热。

• 旧金山和阿尔伯克基的平均气温大致相同，但最低和最高温度差异显著。

23.5.2 化石燃料消费

文件 global-fossil-fuel-consumption.csv 包含 1965 年至 2015 年地球上化石燃料年消费的数据。代码

emissions = pd.read_csv('global-fossil-fuel-consumption.csv')
print(emissions)

打印

    Year         Coal    Crude Oil   Natural Gas
0   1965  16151.96017  18054.69004   6306.370076
1   1966  16332.01679  19442.23715   6871.686791
2   1967  16071.18119  20830.13575   7377.525476
..   ...          ...          ...           ...
50  2015  43786.84580  52053.27008  34741.883490
51  2016  43101.23216  53001.86598  35741.829870
52  2017  43397.13549  53752.27638  36703.965870

现在，让我们将显示每种燃料消费的列替换为两列，一列显示三者的总和，另一列显示五年滚动平均的总和。

emissions['Fuels'] = emissions.sum(axis = 'columns')
emissions.drop(['Coal', 'Crude Oil', 'Natural Gas'], axis = 'columns',
              inplace = True)
num_years = 5
emissions['Roll F'] =\
    emissions['Fuels'].rolling(num_years).mean()
emissions = emissions.round()

我们可以使用

plt.plot(emissions['Year'], emissions['Fuels'],
         label = 'Consumption')
plt.plot(emissions['Year'], emissions['Roll F'],
         label = str(num_years) + ' Year Rolling Ave.')
plt.legend()
plt.title('Consumption of Fossil Fuels')
plt.xlabel('Year')
plt.ylabel('Consumption')

要获取图 23-10 中的情节。

图 23-10 全球化石燃料消费

尽管消费量在少数小幅下滑（例如，2008 年金融危机期间），但上升趋势显而易见。

科学界已达成共识，燃料消费的上升与地球平均温度的上升之间存在关联。让我们看看它与我们在 23.5.1 节中查看的 21 个美国城市的温度之间的关系。

请记住，yearly_temps被绑定到数据框中。

    Year  Min T  Max T  Mean T
0   1961 -17.25  38.05   15.64
1   1962 -21.65  36.95   15.39
2   1963 -24.70  36.10   15.50
..   ...    ...    ...     ...
52  2013 -15.00  40.55   16.66
53  2014 -22.70  40.30   16.85
54  2015 -18.80  40.55   17.54

如果有一种简单的方法来组合yearly_temps和emissions，那该多好啊？Pandas 的merge函数正是如此。代码

yearly_temps['Year'] = yearly_temps['Year'].astype(int)
merged_df = pd.merge(yearly_temps, emissions,
                      left_on = 'Year', right_on = 'Year')
print(merged_df)

打印数据框

    Year  Min T  ...     Fuels    Roll F
0   1965  -21.7  ...   42478.0       NaN
1   1966  -25.0  ...   44612.0       NaN
2   1967  -17.8  ...   46246.0       NaN
..   ...    ...  ...       ...       ...
48  2013  -15.0  ...  131379.0  126466.0
49  2014  -22.7  ...  132028.0  129072.0
50  2015  -18.8  ...  132597.0  130662.0

数据框包含出现在yearly_temps和emissions中的列的并集，但仅包括来自yearly_temps和emissions中具有相同Year值的行构建的行。

现在我们在同一个数据框中拥有排放和温度信息，轻松查看它们之间的相关性。代码print(merged_df.corr().round(2).to_string())打印

 Year  Min T  Max T  Mean T  Fuels  Roll F
Year    1.00   0.37   0.72    0.85   0.99    0.98
Min T   0.37   1.00   0.22    0.49   0.37    0.33
Max T   0.72   0.22   1.00    0.70   0.75    0.66
Mean T  0.85   0.49   0.70    1.00   0.85    0.81
Fuels   0.99   0.37   0.75    0.85   1.00    1.00
Roll F  0.98   0.33   0.66    0.81   1.00    1.00 

我们看到，前几年的全球燃料消费确实与这些美国城市的平均温度和最高温度高度相关。这是否意味着增加的燃料消费导致温度上升？并不是。请注意，两者都与年份高度相关。也许某个潜在变量也与年份相关并且是因果因素。从统计学的角度来看，我们可以说，数据并不反驳广泛接受的科学假设，即化石燃料的增加使用产生的温室气体导致温度上升。

这就是我们对 Pandas 的简要介绍。我们只是在它所提供的内容上轻轻触及了一下。我们将在书中后续部分使用它，并介绍更多功能。如果你想了解更多，有许多在线资源以及一些优秀的廉价书籍。网站 [www.dataschool.io/best-python-pandas-resources/](https://www.dataschool.io/best-python-pandas-resources/) 列出了其中的一些。

23.6 在章节中引入的术语

• 数据框

• 行

• 序列

• 索引

• 名称

• CSV 文件

• ndarray 的形状

• 布尔索引

• 序列的相关性

• 移动（滚动）平均