tnkstat3e-merge-0

统计思维（程序员的概率统计）第三版（一）

原文：allendowney.github.io/ThinkStats/index.html

译者：飞龙

协议：CC BY-NC-SA 4.0

统计思维第三版

原文：allendowney.github.io/ThinkStats/index.html

Think Stats 是面向 Python 程序员的概率论与数理统计入门。如果你有基本的 Python 技能，你可以使用它们来学习概率论与数理统计的概念以及处理数据的相关技能。

第三版现在可以从 Bookshop.org 和 Amazon（这些是联盟链接）购买。如果你正在享受免费的在线版本，考虑买我一杯咖啡。

这是 Green Tea Press 第三版的着陆页面.

• 本书强调了一些简单的技术，你可以使用这些技术来探索真实数据集并回答重要的问题。

• 它包括使用来自美国国立卫生研究院和其他来源的数据集的案例研究。

• 许多练习都使用简短的程序来运行实验，帮助读者发展理解。

本书在 Creative Commons 许可下可用，这意味着你可以自由地复制、分发和修改它，只要您注明出处并且不将其用于商业目的。

笔记本

Jupyter 简介:

• 如果你不太熟悉在 Colab 上运行 Jupyter 笔记本，你可能想从这篇介绍开始。

第一章:

• 在此处运行 Colab 上的第一章

第二章:

• 在此处运行 Colab 上的第二章

第三章:

• 在此处运行 Colab 上的第三章

第四章:

• 在此处运行 Colab 上的第四章

第五章:

• 在此处运行 Colab 上的第五章

第六章:

• 在此处运行 Colab 上的第六章

第七章:

• 在此处运行 Colab 上的第七章

第八章:

• 在此处运行 Colab 上的第八章

第九章:

• 在此处运行 Colab 上的第九章

第十章:

• 在此处运行第十章的 Colab

第十一章:

• 在此处运行第十一章的 Colab

第十二章:

• 在此处运行第十二章的 Colab

第十三章:

• 在此处运行第十三章的 Colab

第十四章:

• 在此处运行第十四章的 Colab

前言

原文：allendowney.github.io/ThinkStats/chap00.html

从统计学最早的起源开始，就有两种关于统计学是什么的观点。一种观点认为，统计学是数学的一个分支，其目标是建立概率和统计推断的理论基础。另一种观点认为，统计学是一套处理数据、回答问题和做出更好决策的工具和实践。许多统计学入门课程基于第一种观点。本书基于第二种观点。

Think Stats 是一本介绍探索和可视化数据、发现关系和趋势以及传达结果的实用方法的入门书。本书的组织结构遵循我开始处理数据集时使用的流程：

• 导入和清理：无论数据格式如何，通常需要一些时间和精力来读取数据、清理和转换它，并检查是否一切在翻译过程中完好无损。

• 单变量探索：我通常一次检查一个变量，找出变量的含义，查看值的分布，并选择适当的汇总统计量。

• 对偶探索：为了识别变量之间可能的关系，我会查看表格和散点图，并计算相关性和线性拟合。

• 多变量分析：如果有变量之间存在明显的关系，我会使用多元回归来添加控制变量并研究更复杂的关系。

• 估计和假设检验：在报告统计结果时，重要的是要回答三个问题：效果有多大？如果我们再次进行相同的测量，我们应该期望多少可变性？明显的效应是否可能是偶然的？

• 可视化：在探索过程中，可视化是寻找可能的关系和效应的重要工具。然后，如果明显的效应经得起推敲，可视化是传达结果的有效方式。

这本书采用计算方法，与更数学的处理方法相比有几个优点：

• 我通常使用 Python 代码而不是数学符号来展示大多数观点。一般来说，Python 代码更易读——此外，因为它可执行，读者可以运行它并修改它来发展洞察力。

• 每章都包括读者可以做的练习，以检查和巩固他们的学习。当你编写程序时，你在代码中表达你的理解——当你调试程序时，你也在检查你的理解。

• 一些练习涉及实验来测试统计行为。例如，你可以通过生成随机样本并计算它们的总和来探索中心极限定理（CLT）。结果的可视化显示了 CLT 为什么有效以及何时无效。

• 一些在数学上难以理解的概念，通过模拟可以更容易地理解。例如，我们通过运行随机模拟来近似 p 值，这加强了假设检验的意义。

• 由于本书基于通用编程语言（Python），读者可以从几乎任何来源导入数据。他们不受特定统计工具已清洗和格式化数据集的限制。

为了展示我的统计分析方法，示例和练习使用了来自多个来源的数据，包括：

• 美国疾病控制与预防中心（CDC）进行的全国家庭成长调查（NSFG），旨在收集有关“家庭生活、婚姻和离婚、怀孕、不孕、避孕药具的使用以及男性和女性健康”的信息。（见www.cdc.gov/nchs/nsfg/index.htm）

• 行为风险因素监测系统（BRFSS），由国家慢性疾病预防和健康促进中心进行，旨在“追踪美国的健康状况和风险行为。”（见cdc.gov/BRFSS/）

• 帕默企鹅数据集，该数据集包括来自南极洲帕默站附近企鹅样本的测量数据（见allisonhorst.github.io/palmerpenguins/）。

• 美国能源信息署（EIA）关于美国可再生能源发电的数据。

我感谢收集这些数据并使其可用的个人和机构——并且我希望与来自各个领域的真实数据一起工作可以使本书对读者更具吸引力。

新增内容？

对于这一版，我首先将本书移入 Jupyter 笔记本。这种变化有一个直接的好处——你可以在一个地方阅读文本、运行代码和完成练习。而且笔记本设计用于在 Google Colab 上工作，因此你可以不安装任何东西就开始使用。

转向笔记本还有另一个好处——代码更易于查看。在前两版中，一些代码在书中，一些代码在在线提供的支持文件中。回顾起来，很明显，以这种方式分割材料并不理想，并且使代码比必要的更复杂。在第三版中，我能够简化代码并使其更易于阅读。

自从上一版出版以来，我开发了一个名为empiricaldist的库，它提供了表示统计分布的对象。这个库现在更加成熟，因此更新的代码更好地利用了它。

当我开始这个项目时，NumPy 和 SciPy 并没有像现在这样广泛使用，Pandas 更是如此，因此原始代码使用了 Python 的数据结构，如列表和字典。这一版广泛使用了数组和 Pandas 结构，并更多地使用了这些库提供的函数。

第三版涵盖了与原版几乎相同的内容，顺序也几乎相同，但文本进行了大幅修订。一些示例是新的；其他示例则更新了新的数据。我开发了新的练习，修改了一些旧的练习，并删除了一些。我认为更新的练习与示例的联系更好，而且更有趣。

自第一版以来，这本书基于这样一个论点：许多难以用数学解释的想法用代码更容易解释。在本版中，我进一步强化了这个观点，以至于几乎不再有数学符号。

总体来说，我认为这些变化使《Think Stats》成为一本更好的书。希望你们喜欢它！

使用代码

本书使用的代码和数据可以从AllenDowney/ThinkStats获取，这是一个 GitHub 上的 Git 仓库。Git 是一个版本控制系统，有助于跟踪构成项目的文件。Git 控制下的文件集合被称为仓库。GitHub 是一个托管服务，为 Git 仓库提供存储空间并提供方便的 Web 界面。

对于本书的每一章，仓库提供了一个 Jupyter 笔记本，这是一个包含文本、代码和代码运行结果的文档。你可以使用这些笔记本来运行代码并完成练习。

你有两种方法可以运行笔记本。到目前为止，最简单的方法是使用 Colab，这是由 Google 提供的一项服务，你可以在 Web 浏览器中运行笔记本而无需在电脑上安装任何东西。如果你从allendowney.github.io/ThinkStats/的Think Stats主页开始，你会找到笔记本的链接，包括一个介绍 Colab 和 Jupyter 笔记本的链接。

如果你不想使用 Colab，你可以下载笔记本并在你的电脑上运行它们，但那样的话，你必须安装 Python、Jupyter 以及书中使用的库，包括 NumPy、SciPy 和 StatsModels。如果你有安装软件的经验，设置一个可以运行笔记本的环境并不困难。但如果你没有这样的经验，你的第一次尝试可能会很具挑战性，有时甚至令人沮丧。在这种情况下，这可能会成为充分利用本书的障碍。如果你想学习 Python 中的探索性数据分析，你不想把时间和认知能力花在安装软件上！

因此，我强烈建议你在 Colab 上至少运行前几章。然后，如果你想设置自己的环境，你可以在不影响书中进度的情况下完成。最后，还有一个建议：如果你在安装软件时遇到任何问题，可以利用像 ChatGPT 这样的工具——它们通常在这些话题上提供很好的指导。

我写这本书假设读者熟悉核心 Python，包括面向对象特性。如果你熟悉 NumPy 和 Pandas，那会有帮助，但不是必需的——我会解释你需要知道的内容。我假设读者了解基本的数学，包括对数，例如，以及求和。你不需要知道线性代数或微积分。有一个地方我提到了导数和积分，但如果你不熟悉这些概念，它们完全是可选的。最后，我不假设你了解任何关于统计学的内容。

致谢

感谢为本书前几版提供纠错和建议的读者，以及那些忍受了粗糙草稿的奥利恩学院的同学们。

感谢本版的审稿人：Zachary del Rosario、Jerzy Wieczorek、Thomas Nield、Walter Paczkowski 和 Peter Bruce。

还要感谢奥莱利媒体的所有人，特别是编辑 Sara Hunter 和 Aaron Black，以及…待定

《Think Stats：Python 中的探索性数据分析，第 3 版

版权所有 2024 艾伦·B·唐尼

代码许可：MIT 许可证

文本许可：Creative Commons Attribution-NonCommercial-ShareAlike 4.0 国际

探索性数据分析

原文：allendowney.github.io/ThinkStats/chap01.html

这本书的论点是我们可以使用数据来回答问题、解决辩论和做出更好的决策。

本章介绍了我们将要使用的步骤：加载数据和验证数据、探索数据，并选择衡量我们感兴趣内容的统计方法。作为一个例子，我们将使用国家家庭成长调查（NSFG）的数据来回答当我和妻子期待我们第一个孩子时听到的问题：第一个孩子是否倾向于晚出生？

点击此处运行此笔记本在 Colab 上。

from  os.path  import basename, existsdef  download(url):filename = basename(url)if not exists(filename):from  urllib.request  import urlretrievelocal, _ = urlretrieve(url, filename)print("Downloaded " + local)download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/thinkstats.py") 

try:import  empiricaldist
except ImportError:%pip install empiricaldist 

import  numpy  as  np
import  pandas  as  pd
import  matplotlib.pyplot  as  pltfrom  IPython.display  import HTML
from  thinkstats  import decorate 

证据

你可能听说过第一个孩子更有可能晚出生。如果你用这个问题在网络上搜索，你会找到很多讨论。有些人声称这是真的，有些人说这是一个神话，还有一些人说情况正好相反：第一个孩子会提前出生。

在许多这样的讨论中，人们提供数据来支持他们的主张。我发现了很多这样的例子：

“我最近有两个朋友生了第一个孩子，他们两个在分娩或被诱导分娩前都几乎迟到了两周。”

“我的第一个孩子迟到了两周，现在我觉得第二个孩子可能会提前两周出生！！！”

“我不认为这是真的，因为我的姐姐是我母亲第一个孩子，她很早就出生了，就像我的许多表亲一样。”

这样的报告被称为轶事证据，因为它们基于未发表的数据，通常是个人数据。在闲聊中，轶事并没有什么不妥，所以我并不是要批评我引用的人。

但我们可能需要更有说服力的证据和更可靠的答案。按照这些标准，轶事证据通常是不够的，因为：

• 观察数量少：如果第一个孩子的怀孕时间更长，那么与自然变异相比，差异可能很小。在这种情况下，我们可能需要比较大量的怀孕案例来了解是否存在差异。

• 选择偏误：加入关于这个问题讨论的人可能是因为他们的第一个孩子晚出生而感兴趣的。在这种情况下，选择数据的过程可能会使结果产生偏差。

• 确认偏误：相信这个说法的人可能更倾向于提供证实它的例子。怀疑这个说法的人更有可能引用反例。

• 不准确性：轶事通常是个人故事，可能会被错误地记住、错误地表述、不准确地重复等。

为了解决轶事的局限性，我们将使用统计工具，这包括：

• 数据收集：我们将使用来自一个大型国家调查的数据，该调查明确旨在生成关于美国人口的统计有效推断。

• 描述性统计：我们将生成总结数据的统计量，并评估不同的数据可视化方法。

• 探索性数据分析：我们将寻找模式、差异和其他特征，以解决我们感兴趣的问题。同时，我们将检查不一致性并确定局限性。

• 估计：我们将使用样本数据来估计总体特征。

• 假设检验：当我们看到明显的效应，如两组之间的差异时，我们将评估这种效应是否可能偶然发生。

通过小心执行这些步骤以避免陷阱，我们可以得出更加合理且更有可能正确的结论。

全国家庭成长调查

自 1973 年以来，美国疾病控制与预防中心（CDC）一直在进行全国家庭成长调查（NSFG），旨在收集有关“家庭生活、婚姻与离婚、怀孕、不孕、避孕药具的使用以及男女健康”的信息。调查结果被用于……规划健康服务和健康教育计划，以及进行家庭、生育和健康方面的统计分析。

你可以在cdc.gov/nchs/nsfg.htm了解更多关于 NSFG 的信息。

我们将使用此调查收集的数据来调查是否第一个孩子倾向于晚出生，以及其他问题。为了有效地使用这些数据，我们必须了解研究的设计。

通常，统计研究的目标是关于总体的结论。在 NSFG 中，目标总体是美国 15-44 岁的人群。

理想情况下，调查将收集总体中每个成员的数据，但这很少可能实现。相反，我们收集来自称为样本的总体子集的数据。参与调查的人被称为受访者。

NSFG 是一项横断面研究，这意味着它捕捉了在某一时间点的总体快照。NSFG 已经进行了几次；每次部署被称为周期。我们将使用第 6 周期的数据，该周期从 2002 年 1 月持续到 2003 年 3 月。

通常，横断面研究旨在是代表性的，这意味着样本在所有对研究目的重要方面都与目标总体相似。在实践中实现这一理想很困难，但进行调查的人会尽可能接近。

美国国家家庭生育调查（NSFG）不具有代表性；相反，它是分层的，这意味着它故意过度抽样某些群体。该研究的制定者招募了三个群体——西班牙裔、非裔美国人和青少年——其在美国人口中的代表性高于其比例，以确保每个群体中受访者的数量足够多，以便得出有效的结论。过度抽样的缺点是，基于样本的统计数据，基于人口得出结论并不那么容易。我们稍后会回到这个点。

当处理这类数据时，熟悉代码簿非常重要，它记录了研究的设计、调查问题和响应的编码。

NSFG 数据的代码簿和用户指南可在www.cdc.gov/nchs/nsfg/nsfg_cycle6.htm获取。

数据读取

在下载 NSFG 数据之前，您必须同意使用条款：

任何对个人或机构的故意识别或披露都违反了向信息提供者提供的保密保证。因此，用户将：

• 仅使用此数据集中的数据进行统计分析。
• 不要试图了解这些数据中包含的任何个人或机构的身份。
• 不要将此数据集与来自其他 NCHS 或非 NCHS 数据集的个人可识别数据链接。
• 不要参与评估用于保护个人和机构的披露方法或任何关于个人和机构重新识别方法的研究。

如果您同意遵守这些条款，本章笔记本中提供了下载数据的说明。

数据文件可以直接从 NSFG 网站获取，网址为www.cdc.gov/nchs/data_access/ftp_dua.htm?url_redirect=ftp://ftp.cdc.gov/pub/Health_Statistics/NCHS/Datasets/NSFG，但我们将从本书的存储库中下载它们，该存储库提供了数据文件的压缩版本。

以下单元格下载数据文件并安装statadict，这是我们读取数据所需的。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dct")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dat.gz") 

try:import  statadict
except ImportError:%pip install statadict 

数据存储在两个文件中，一个“字典”描述了数据的格式，一个数据文件。

dct_file = "2002FemPreg.dct"
dat_file = "2002FemPreg.dat.gz" 

本章笔记本定义了一个函数，用于读取这些文件。它被称为read_stata，因为这种数据格式与名为 Stata 的统计软件包兼容。

以下函数将这些文件名作为参数，读取字典，并使用结果读取数据文件。

from  statadict  import parse_stata_dictdef  read_stata(dct_file, dat_file):stata_dict = parse_stata_dict(dct_file)resp = pd.read_fwf(dat_file,names=stata_dict.names,colspecs=stata_dict.colspecs,compression="gzip",)return resp 

这是我们的使用方法。

preg = read_stata(dct_file, dat_file) 

结果是一个DataFrame，这是 Pandas 数据结构，用于表示行和列的表格数据。这个DataFrame为每个受访者报告的怀孕包含一行，为每个变量包含一列。一个变量可以包含对调查问题的回答或基于一个或多个问题的回答计算出的值。

除了数据之外，DataFrame还包含变量名及其类型，并提供访问和修改数据的方法。DataFrame有一个名为shape的属性，其中包含行数和列数。

preg.shape 

(13593, 243) 

该数据集包含 243 个变量，涉及 13,593 次怀孕的信息。DataFrame提供了一个名为head的方法，用于显示前几行。

preg.head() 

| | caseid | pregordr | howpreg_n | howpreg_p | moscurrp | nowprgdk | pregend1 | pregend2 | nbrnaliv | multbrth | ... | poverty_i | laborfor_i | religion_i | metro_i | basewgt | adj_mod_basewgt | finalwgt | secu_p | sest | cmintvw |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 0 | 1 | 1 | NaN | NaN | NaN | NaN | 6.0 | NaN | 1.0 | NaN | ... | 0 | 0 | 0 | 0 | 3410.389399 | 3869.349602 | 6448.271112 | 2 | 9 | 1231 |
| 1 | 1 | 2 | NaN | NaN | NaN | NaN | 6.0 | NaN | 1.0 | NaN | ... | 0 | 0 | 0 | 0 | 3410.389399 | 3869.349602 | 6448.271112 | 2 | 9 | 1231 |
| 2 | 2 | 1 | NaN | NaN | NaN | NaN | 5.0 | NaN | 3.0 | 5.0 | ... | 0 | 0 | 0 | 0 | 7226.301740 | 8567.549110 | 12999.542264 | 2 | 12 | 1231 |
| 3 | 2 | 2 | NaN | NaN | NaN | NaN | 6.0 | NaN | 1.0 | NaN | ... | 0 | 0 | 0 | 0 | 7226.301740 | 8567.549110 | 12999.542264 | 2 | 12 | 1231 |
| 4 | 2 | 3 | NaN | NaN | NaN | NaN | 6.0 | NaN | 1.0 | NaN | ... | 0 | 0 | 0 | 0 | 7226.301740 | 8567.549110 | 12999.542264 | 2 | 12 | 1231 |

5 行 × 243 列

左侧列是DataFrame的索引，其中包含每行的标签。在这种情况下，标签是从 0 开始的整数，但它们也可以是字符串和其他类型。

DataFrame有一个名为columns的属性，其中包含变量的名称。

preg.columns 

Index(['caseid', 'pregordr', 'howpreg_n', 'howpreg_p', 'moscurrp', 'nowprgdk','pregend1', 'pregend2', 'nbrnaliv', 'multbrth',...'poverty_i', 'laborfor_i', 'religion_i', 'metro_i', 'basewgt','adj_mod_basewgt', 'finalwgt', 'secu_p', 'sest', 'cmintvw'],dtype='object', length=243) 

列名包含在一个Index对象中，这是另一个 Pandas 数据结构。要从DataFrame中访问列，可以使用列名作为键。

pregordr = preg["pregordr"]
type(pregordr) 

pandas.core.series.Series 

结果是 Pandas 的Series，它表示一系列值。Series也提供了head方法，用于显示前几个值及其标签。

pregordr.head() 

0    1
1    2
2    1
3    2
4    3
Name: pregordr, dtype: int64 

最后一行包含Series的名称和dtype，即值的类型。在这个例子中，int64表示值是 64 位整数。

NSFG 数据集总共包含 243 个变量。以下是本书探索中将使用的一些变量。

• caseid是受访者的整数 ID。

• pregordr是怀孕序列号：受访者第一次怀孕的代码是 1，第二次怀孕是 2，依此类推。

• prglngth是怀孕的整数持续时间（周）。

• outcome 是妊娠结果的整数代码。代码 1 表示活产。

• birthord 是活产的序列号：受访者的第一个孩子的代码是 1，依此类推。对于除活产之外的结果，此字段为空。

• birthwgt_lb 和 birthwgt_oz 包含婴儿出生体重的磅和盎司部分。

• agepreg 是妊娠结束时的母亲年龄。

• finalwgt 是与受访者相关的统计权重。它是一个浮点值，表示该受访者代表美国人口中的多少人。

如果仔细阅读代码簿，你会看到许多变量是重新编码，这意味着它们不是调查收集的原始数据的一部分——它们是使用原始数据计算得出的。

例如，对于活产，如果可用，prglngth 等于原始变量 wksgest（妊娠周数）；否则，它使用 mosgest * 4.33（妊娠月数乘以每月平均周数）来估算。

重新编码通常基于检查数据的完整性和准确性的逻辑。一般来说，当可用时使用重新编码是一个好主意，除非有充分的理由自行处理原始数据。

验证

当数据从一个软件环境导出并导入到另一个软件环境时，可能会引入错误。当你熟悉新的数据集时，你可能会错误地解码数据或误解其含义。如果你投入时间验证数据，你可以在以后节省时间并避免错误。

验证数据的一种方法是通过计算基本统计量并将它们与已发布的结果进行比较。例如，NSFG 代码簿包括总结每个变量的表格。以下是 outcome 的表格，它编码了每个妊娠的结果。

值	标签	总计
1	活产	9148
2	人工流产	1862
3	死产	120
4	流产	1921
5	宫外孕	190
6	当前妊娠	352
总计		13593

“总计”列表示具有每种结果的妊娠数量。为了检查这些总计，我们将使用 value_counts 方法，该方法计算每个值出现的次数，以及 sort_index，该方法根据 Index（左侧列）中的值对结果进行排序。

preg["outcome"].value_counts().sort_index() 

outcome
1    9148
2    1862
3     120
4    1921
5     190
6     352
Name: count, dtype: int64 

将结果与已发布的表格进行比较，我们可以确认 outcome 中的值是正确的。同样，以下是 birthwgt_lb 的已发布表格。

值	标签	总计
.	不适用	4449
0-5	低于 6 磅	1125
6	6 磅	2223
7	7 磅	3049
8	8 磅	1889
9-95	9 磅或以上	799
97	未确定	1
98	拒绝	1
99	未知	57
总计		13593

出生重量仅记录在以活产结束的怀孕中。表格表明，有 4449 个案例中这个变量不适用。此外，还有一个案例中未提问，一个案例中受访者未回答，以及 57 个案例中他们不知道。

再次，我们可以使用value_counts来比较数据集中的计数与代码簿中的计数。

counts = preg["birthwgt_lb"].value_counts(dropna=False).sort_index()
counts 

birthwgt_lb
0.0        8
1.0       40
2.0       53
3.0       98
4.0      229
5.0      697
6.0     2223
7.0     3049
8.0     1889
9.0      623
10.0     132
11.0      26
12.0      10
13.0       3
14.0       3
15.0       1
51.0       1
97.0       1
98.0       1
99.0      57
NaN     4449
Name: count, dtype: int64 

参数dropna=False意味着value_counts不会忽略“NA”或“不适用”的值。这些值在结果中显示为NaN，代表“不是一个数字”——并且这些值的计数与代码簿中不适用的案例计数一致。

6、7 和 8 磅的计数与代码簿一致。为了检查 0 至 5 磅范围内的计数，我们可以使用一个名为loc的属性——它是“位置”的缩写——以及一个切片索引来选择计数的一个子集。

counts.loc[0:5] 

birthwgt_lb
0.0      8
1.0     40
2.0     53
3.0     98
4.0    229
5.0    697
Name: count, dtype: int64 

我们还可以使用sum方法将它们加起来。

counts.loc[0:5].sum() 

np.int64(1125) 

总数与代码簿一致。

值 97、98 和 99 代表出生重量未知的情况。我们可能有几种处理缺失数据的方法。一个简单的选项是将这些值替换为NaN。同时，我们还将替换一个明显错误的值，51 磅。

我们可以使用这种方法使用replace方法：

preg["birthwgt_lb"] = preg["birthwgt_lb"].replace([51, 97, 98, 99], np.nan) 

第一个参数是要替换的值的列表。第二个参数np.nan从 NumPy 获取NaN值。

当你以这种方式读取数据时，你通常需要检查错误并处理特殊值。这种操作被称为数据清洗。

转换

作为另一种数据清洗方式，有时我们必须将数据转换为不同的格式，并执行其他计算。

例如，agepreg包含怀孕结束时的母亲年龄。根据代码簿，它是一个以百分之一年（百分之一年）为单位的整数，正如我们可以通过使用mean方法计算其平均值来得知。

preg["agepreg"].mean() 

np.float64(2468.8151197039497) 

为了将其转换为年，我们可以将其除以 100。

preg["agepreg"] /= 100.0
preg["agepreg"].mean() 

np.float64(24.6881511970395) 

现在平均值更可信了。

作为另一个例子，birthwgt_lb和birthwgt_oz包含出生重量，磅和盎司分别在不同的列中。将它们合并成一个包含磅和磅分数的单列将更方便。

首先，我们将像处理birthwgt_lb一样清洗birthwgt_oz。

preg["birthwgt_oz"].value_counts(dropna=False).sort_index() 

birthwgt_oz
0.0     1037
1.0      408
2.0      603
3.0      533
4.0      525
5.0      535
6.0      709
7.0      501
8.0      756
9.0      505
10.0     475
11.0     557
12.0     555
13.0     487
14.0     475
15.0     378
97.0       1
98.0       1
99.0      46
NaN     4506
Name: count, dtype: int64 

preg["birthwgt_oz"] = preg["birthwgt_oz"].replace([97, 98, 99], np.nan) 

现在我们可以使用清洗后的值来创建一个新的列，该列将磅和盎司合并成一个单一的数量。

preg["totalwgt_lb"] = preg["birthwgt_lb"] + preg["birthwgt_oz"] / 16.0
preg["totalwgt_lb"].mean() 

np.float64(7.265628457623368) 

结果的平均值看起来是合理的。

摘要统计

统计量是从数据集中派生出的一个数字，通常旨在量化数据的某个方面。例如包括计数、平均值、方差和标准差。

Series对象有一个count方法，它返回非nan值的数量。

weights = preg["totalwgt_lb"]
n = weights.count()
n 

np.int64(9038) 

它还提供了一个返回值总和的sum方法——我们可以用它来计算平均值，如下所示。

mean = weights.sum() / n
mean 

np.float64(7.265628457623368) 

但正如我们已经看到的，还有一个 mean 方法可以做到同样的事情。

weights.mean() 

np.float64(7.265628457623368) 

在这个数据集中，平均出生体重约为 7.3 磅。

方差是一种统计量，用于量化一组值的分布。它是平方偏差的平均值，即每个点与平均值之间的距离。

squared_deviations = (weights - mean) ** 2 

我们可以这样计算平方偏差的平均值。

var = squared_deviations.sum() / n
var 

np.float64(1.983070989750022) 

如您所预期，Series 提供了一个 var 方法，它几乎做的是同样的事情。

weights.var() 

np.float64(1.9832904288326545) 

结果略有不同，因为当 var 方法计算平方偏差的平均值时，它除以 n-1 而不是 n。这是因为根据你试图做什么，有两种计算样本方差的方法。我将在 第八章 中解释这种差异——但在实践中通常无关紧要。如果你更喜欢分母中有 n 的版本，你可以通过将 ddof=0 作为关键字参数传递给 var 方法来获得它。

weights.var(ddof=0) 

np.float64(1.983070989750022) 

在这个数据集中，出生体重的方差大约是 1.98，但这个值很难解释——一方面，它是以磅平方为单位的。方差在某些计算中很有用，但不是描述数据集的好方法。更好的选择是 标准差，它是方差的平方根。我们可以这样计算它。

std = np.sqrt(var)
std 

np.float64(1.40821553384062) 

或者，我们可以使用 std 方法。

weights.std(ddof=0) 

np.float64(1.40821553384062) 

在这个数据集中，出生体重的标准差约为 1.4 磅。非正式地说，距离平均值一个或两个标准差的值很常见——距离平均值更远的值很少见。

解释

为了有效地处理数据，你必须同时从两个层面思考：统计层面和背景层面。例如，让我们选择怀孕文件中 caseid 为 10229 的行。query 方法接受一个字符串，可以包含列名、比较运算符和数字等。

subset = preg.query("caseid == 10229")
subset.shape 

(7, 244) 

结果是一个只包含查询为 True 的行的 DataFrame。这位受访者报告了七次怀孕——以下是他们的结果，这些结果按时间顺序记录。

subset["outcome"].values 

array([4, 4, 4, 4, 4, 4, 1]) 

结果代码 1 表示活产。代码 4 表示流产——即怀孕丢失，通常没有已知的医学原因。

从统计学的角度来看，这位受访者并不异常。怀孕丢失很常见，还有其他受访者报告了同样多的实例。但考虑到背景，这些数据讲述了一个女人怀孕六次，每次都以流产告终的故事。她的第七次也是最最近的一次怀孕以活产结束。如果我们带着同情心考虑这些数据，那么这个故事自然会触动人心。

NSFG 数据集中的每一行代表一个对许多个人和困难问题给出诚实回答的人。我们可以使用这些数据来回答关于家庭生活、生育和健康的统计问题。同时，我们有责任考虑数据所代表的人，并给予他们尊重和感激。

术语表

每章的结尾提供了一个定义在章节中的词汇表。

• 轶事证据：从少量个体案例中非正式收集的数据，通常没有系统抽样。

• 横断面研究：在某一时间点或时间间隔内从代表性样本中收集数据的调查。

• 周期：在一个在多个时间间隔收集数据的调查中，一个数据收集间隔。

• 总体：研究主题的整个个体或项目群体。

• 样本：总体的一部分，通常随机选择。

• 受访者：参与调查并回答问题的人。

• 代表性：如果一个样本在研究目的上重要的方面与总体相似，则该样本是具有代表性的。

• 分层抽样：如果一个样本故意对某些群体进行过抽样，通常是为了确保包含足够的成员以支持有效的结论，则该样本是分层的。

• 过抽样：如果一个群体的成员在样本中出现的概率更高，则该群体是过抽样的。

• 变量：在调查数据中，变量是对问题的回答或从回答中计算出的值的集合。

• 代码簿：描述数据集中变量的文档，并提供有关数据的其他信息。

• 重新编码：基于数据集中其他变量的变量。

• 原始数据：收集后未经处理的原始数据。

• 数据清洗：识别和纠正数据集中错误的过程，处理缺失值，并计算重新编码。

• 统计量：描述或总结样本属性的值。

• 标准差：一个统计量，用于量化数据围绕平均值的分布。

练习

本章的练习基于 NSFG 怀孕文件。

练习 1.1

从preg中选择birthord列，打印值计数，并与ftp.cdc.gov/pub/Health_Statistics/NCHS/Dataset_Documentation/NSFG/Cycle6Codebook-Pregnancy.pdf中发布的代码簿结果进行比较。

练习 1.2

创建一个名为totalwgt_kg的新列，其中包含以千克为单位的出生体重（每千克大约有 2.2 磅）。计算新列的平均值和标准差。

练习 1.3

对于caseid为 2298 的受访者，其怀孕时长是多少？

caseid为 5013 的受访者第一个孩子的出生体重是多少？提示：您可以使用and在查询中检查多个条件。

Think Stats: Python 中的探索性数据分析，第 3 版

版权所有 2024 艾伦·B·唐尼

代码许可：MIT 许可证

文本许可：Creative Commons 知识共享署名-非商业性使用-相同方式共享 4.0 国际

分布

原文链接

本章介绍了统计学中最基本的思想之一，即分布。我们将从频率表开始——它表示数据集中的值及其出现的次数——并使用它们来探索来自国家家庭增长调查（NSFG）的数据。我们还将寻找极端或错误值，称为异常值，并考虑处理它们的方法。

在此处运行此笔记本。

from  os.path  import basename, existsdef  download(url):filename = basename(url)if not exists(filename):from  urllib.request  import urlretrievelocal, _ = urlretrieve(url, filename)print("Downloaded " + local)download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/thinkstats.py") 

try:import  empiricaldist
except ImportError:%pip install empiricaldist 

import  numpy  as  np
import  matplotlib.pyplot  as  plt
from  thinkstats  import decorate 

频率表

描述变量的一个方法是一个频率表，它包含变量的值及其频率——也就是说，每个值出现的次数。这种描述被称为变量的分布。

为了表示分布，我们将使用一个名为 empiricaldist 的库。在这个上下文中，“empirical”意味着这些分布是基于数据而不是数学模型。empiricaldist 提供了一个名为 FreqTab 的类，我们可以使用它来计算和绘制频率表。我们可以这样导入它。

from  empiricaldist  import FreqTab 

为了展示它是如何工作的，我们将从一个小的值列表开始。

t = [1.0, 2.0, 2.0, 3.0, 5.0] 

FreqTab 提供了一个名为 from_seq 的方法，它接受一个序列并将其转换为 FreqTab 对象。

ftab = FreqTab.from_seq(t)
ftab 

	freqs
1.0	1
2.0	2
3.0	1
5.0	1

FreqTab 对象是一种包含值及其频率的 Pandas Series。在这个例子中，值 1.0 对应频率 1，值 2.0 对应频率 2，等等。

FreqTab 提供了一个名为 bar 的方法，它将频率表绘制为条形图。

ftab.bar()
decorate(xlabel="Value", ylabel="Frequency") 

因为 FreqTab 是 Pandas Series，我们可以使用方括号运算符来查找值并获取其频率。

ftab[2.0] 

np.int64(2) 

但与 Pandas Series 不同，我们还可以像函数一样调用 FreqTab 对象来查找值。

ftab(2.0) 

np.int64(2) 

如果我们查找 FreqTab 中不存在的值，函数语法将返回 0。

ftab(4.0) 

0 

FreqTab 对象有一个名为 qs 的属性，它返回一个值数组——qs 代表数量，尽管技术上并非所有值都是数量。

ftab.qs 

array([1., 2., 3., 5.]) 

FreqTab 还有一个名为 fs 的属性，它返回一个频率数组。

ftab.fs 

array([1, 2, 1, 1]) 

FreqTab 提供了一个 items 方法，我们可以使用它来遍历数量-频率对：

for x, freq in ftab.items():print(x, freq) 

1.0 1
2.0 2
3.0 1
5.0 1 

我们将在继续的过程中看到更多的 FreqTab 方法。

NSFG 分布

当你开始使用一个新的数据集时，我建议你逐个探索你计划使用的变量，一个好的开始方式是查看频率表。

作为例子，让我们看看国家家庭增长调查（NSFG）的数据。在前一章中，我们下载了这个数据集，将其读入 Pandas DataFrame，并对几个变量进行了清理。我们用来加载数据和清理数据的代码在一个名为nsfg.py的模块中——关于安装此模块的说明在本章的笔记本中。

以下单元格下载数据文件并安装statadict，这是我们读取数据所需的。

try:import  statadict
except ImportError:%pip install statadict 

download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/nsfg.py")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dct")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dat.gz") 

我们可以这样导入它并读取妊娠文件。

from  nsfg  import read_fem_pregpreg = read_fem_preg() 

在本章的例子中，我们将关注以活产结束的妊娠。我们可以使用query方法来选择outcome为 1 的行。

live = preg.query("outcome == 1") 

在传递给query的字符串中，像outcome这样的变量名指的是DataFrame中的列名。这个字符串还可以包含像==这样的运算符和像1这样的操作数。

现在我们可以使用FreqTab.from_seq来计算每个数量在birthwgt_lb（出生体重的磅部分）中出现的次数。name参数给FreqTab对象一个名称，当我们在绘图时用作标签。

ftab_lb = FreqTab.from_seq(live["birthwgt_lb"], name="birthwgt_lb") 

这里是分布的形状。

ftab_lb.bar()
decorate(xlabel="Pounds", ylabel="Frequency") 

看到这样的分布，我们首先注意到的是形状，它类似于著名的钟形曲线，更正式地称为正态分布或高斯分布。分布的另一个显著特征是峰值，即最常见的值。要找到峰值，我们可以使用idxmax方法，它找到与最高频率相关的数量。

ftab_lb.idxmax() 

np.float64(7.0) 

FreqTab提供了一个名为mode的方法，它做的是同样的事情。

ftab_lb.mode() 

np.float64(7.0) 

在这个分布中，峰值在 7 磅。

作为另一个例子，这是birthwgt_oz的频率表，它是出生体重的盎司部分。

ftab_oz = FreqTab.from_seq(live["birthwgt_oz"], name="birthwgt_oz")
ftab_oz.bar()
decorate(xlabel="Ounces", ylabel="Frequency") 

因为自然界不知道磅和盎司，我们可能会预期birthwgt_oz的所有值出现的可能性都是相同的——也就是说，这个分布应该是均匀分布。但看起来0比其他数值更常见，而1和15则不那么常见，这表明受访者会将接近整数磅的出生体重四舍五入。

作为另一个例子，让我们看看agepreg的频率表，这是妊娠结束时母亲的年龄。

ftab_age = FreqTab.from_seq(live["agepreg"], name="agepreg") 

在 NSFG 中，年龄以年和月为单位记录，因此比我们之前看到的分布有更多的唯一值。因此，我们将width=0.1作为关键字参数传递给bar方法，该方法调整条形的宽度，以便它们不会重叠太多。

ftab_age.bar(width=0.1)
decorate(xlabel="Age", ylabel="Frequency") 

分布非常粗略地呈钟形，但它是右偏的——也就是说，尾巴比左边延伸得更远。

最后，让我们看看prglngth的频率表，这是怀孕的周数。xlim参数将 x 轴的极限设置为从 20 周到 50 周的范围——这个范围之外没有很多值，它们可能是错误的。

ftab_length = FreqTab.from_seq(live["prglngth"], name="prglngth")
ftab_length.bar()
decorate(xlabel="Weeks", ylabel="Frequency", xlim=[20, 50]) 

最常见的数量是 39 周。左尾比右尾长——早产的婴儿很常见，但怀孕很少超过 43 周，如果确实如此，医生通常会干预。

异常值

查看频率表，很容易识别分布的形状和最常见的数量，但稀有数量并不总是显而易见。在继续之前，检查异常值是个好主意，这些异常值可能是测量或记录错误，也可能是对稀有事件的准确报告。

为了识别异常值，以下函数接受一个FreqTab对象和一个整数n，并使用切片索引来选择n个最小的数量及其频率。

def  smallest(ftab, n=10):return ftab[:n] 

在prglngth的频率表中，这里有 10 个最小的值。

smallest(ftab_length) 

prglngth
0     1
4     1
9     1
13    1
17    2
18    1
19    1
20    1
21    2
22    7
Name: prglngth, dtype: int64 

由于我们选择了活产行，怀孕长度小于 10 周肯定是错误。最可能的解释是结果没有被正确编码。长度超过 30 周可能是合法的。在 10 到 30 周之间，很难确定——一些数量可能是错误，但一些是正确记录的早产儿。

以下函数从FreqTab对象中选择最大值。

def  largest(ftab, n=10):return ftab[-n:] 

这是数据集中最长的怀孕长度。

largest(ftab_length) 

prglngth
40    1116
41     587
42     328
43     148
44      46
45      10
46       1
47       1
48       7
50       2
Name: prglngth, dtype: int64 

再次强调，其中一些值可能是错误。大多数医生建议如果怀孕超过 41 周就进行引产，所以 50 周似乎不太可能是正确的。但是，在确定错误值和可能为稀有事件正确报告的值之间没有明确的界限。

处理异常值最好的方法取决于“领域知识”——也就是说，关于数据来源及其含义的信息。并且它取决于你计划执行的分析。

在这个例子中，激励问题是初生儿是否比其他婴儿更早或更晚。所以我们将使用不会因少量错误值而受到太大影响的统计数据。

初生儿

现在让我们比较初生儿和其他人的怀孕长度分布。我们可以使用query方法来选择代表初生儿和其他人的行。

firsts = live.query("birthord == 1")
others = live.query("birthord != 1") 

为每个组制作一个FreqTab的怀孕长度表。

ftab_first = FreqTab.from_seq(firsts["prglngth"], name="firsts")
ftab_other = FreqTab.from_seq(others["prglngth"], name="others") 

以下函数可以并排绘制两个频率表。

def  two_bar_plots(ftab1, ftab2, width=0.45):ftab1.bar(align="edge", width=-width)ftab2.bar(align="edge", width=width, alpha=0.5) 

这是它们的分布情况。

two_bar_plots(ftab_first, ftab_other)
decorate(xlabel="Weeks", ylabel="Frequency", xlim=[20, 50]) 

分布形状或异常值之间没有明显的差异。看起来更多的非初产妇在 39 周出生，但数据集中有更多的非初产妇，所以我们不应该直接比较计数。

firsts["prglngth"].count(), others["prglngth"].count() 

(np.int64(4413), np.int64(4735)) 

比较分布的均值，初产妇的平均出生时间似乎稍晚一些。

first_mean = firsts["prglngth"].mean()
other_mean = others["prglngth"].mean()
first_mean, other_mean 

(np.float64(38.60095173351461), np.float64(38.52291446673706)) 

但差异仅为 0.078 周，大约是 13 小时。

diff = first_mean - other_mean
diff, diff * 7 * 24 

(np.float64(0.07803726677754952), np.float64(13.11026081862832)) 

这种明显的差异可能有几个可能的原因：

• 首胎和其他人之间在平均怀孕长度上可能存在实际差异。

• 我们在这个数据集中看到的这种明显的差异可能是抽样过程中的偏差的结果——也就是说，调查受访者的选择。

• 这种明显的差异可能是测量误差的结果——例如，初产妇或其他人的自我报告的怀孕长度可能更准确。

• 这种明显的差异可能是抽样过程中的随机变化的结果。

在后面的章节中，我们将更仔细地考虑这些可能的原因，但现在我们将这个结果视为事实：在这个数据集中，这些组之间在怀孕长度上存在微小差异。

效应量

这种差异有时被称为“效应”。有几种方法可以量化效应的大小。最简单的是报告绝对差异——在这个例子中，差异是 0.078 周。

另一种方法是报告相对差异。例如，我们可能会说，初产妇的平均怀孕时间比其他人长 0.2%。

diff / live["prglngth"].mean() * 100 

np.float64(0.20237586646738304) 

另一种选择是报告一个标准化的效应量，这是一个旨在以不同数量和不同组之间可比较的方式量化效应大小的统计量。

标准化意味着我们以标准差的倍数来表示差异。因此，我们可能会倾向于写一些像这样的事情。

diff / live["prglngth"].std() 

np.float64(0.028877623375210403) 

但请注意，我们使用了两组来计算标准差。如果两组有显著差异，当我们把它们放在一起时，标准差会比任何一个组都大，这可能会使效应量看起来很小。

另一种选择是只使用一个组的标准差，但不确定是哪一个。因此，我们可以取两个标准差的平均值，但如果组的大小不同，这会给一个组过多的权重，而给另一个组不足的权重。

一种常见的解决方案是使用合并标准差，它是合并方差的平方根，而合并方差是各组方差的加权总和。为了计算它，我们将从方差开始。

group1, group2 = firsts["prglngth"], others["prglngth"] 

v1, v2 = group1.var(), group2.var() 

这里是加权总和，组大小作为权重。

n1, n2 = group1.count(), group2.count()
pooled_var = (n1 * v1 + n2 * v2) / (n1 + n2) 

最后，这是合并标准差。

np.sqrt(pooled_var) 

np.float64(2.7022108144953862) 

合并标准差位于各组标准差之间。

firsts["prglngth"].std(), others["prglngth"].std() 

(np.float64(2.7919014146687204), np.float64(2.6158523504392375)) 

使用合并标准差的标准效应量被称为Cohen 效应量。这里有一个计算它的函数。

def  cohen_effect_size(group1, group2):diff = group1.mean() - group2.mean()v1, v2 = group1.var(), group2.var()n1, n2 = group1.count(), group2.count()pooled_var = (n1 * v1 + n2 * v2) / (n1 + n2)return diff / np.sqrt(pooled_var) 

下面是平均怀孕长度差异的效应量。

cohen_effect_size(firsts["prglngth"], others["prglngth"]) 

np.float64(0.028879044654449834) 

在本例中，差异为 0.029 个标准差，这是一个小的差异。为了更好地理解这一点，男性和女性之间身高的差异大约是 1.7 个标准差。

报告结果

我们已经看到了几种描述第一胎和其他婴儿之间怀孕长度差异（如果有的话）的方法。我们应该如何报告这些结果？

答案取决于提问者是谁。科学家可能对任何（真实）效应都感兴趣，无论大小如何。医生可能只关心那些实际有意义的效应——也就是说，在实践中有意义的差异。一个孕妇可能对与她相关的研究结果感兴趣，比如早产的几率或晚产的几率。

你如何报告结果也取决于你的目标。如果你试图证明一个效应的重要性，你可能会选择强调差异的汇总统计量。如果你试图安慰一个病人，你可能会选择将差异置于背景中的统计量。

当然，你的决定也应该受到专业伦理的指导。说服人是可以的——你应该设计清晰讲述故事的统计报告和可视化。但你也应该尽力使你的报告诚实，并承认不确定性和局限性。

术语表

• 分布：数据集中值及其出现频率的集合。

• 频率表：从值到频率的映射。

• 频率：一个值在样本中出现的次数。

• 偏态：如果一个分布是不对称的，且极端数值在一个方向上比另一个方向延伸得更远，则该分布是偏态的。

• 众数：样本中最频繁的数值，或最频繁的数值之一。

• 均匀分布：一个分布，其中所有数值具有相同的频率。

• 异常值：分布中的一个极端数值。

• 标准化：如果一个统计量以不同数据集和领域可比的术语表示，则该统计量是标准化的。

• 合并标准差：一个统计量，它结合了两个或更多组的数据来计算一个共同的标准差。

• Cohen 效应量：一个标准化的统计量，用于量化两组均值之间的差异。

• 实际有意义的：如果一个效应足够大，在实践中有意义，则该效应是实际有意义的。

练习

对于本章的练习，我们将加载 NSFG 女性受访者文件，该文件包含每个女性受访者的一个记录。下载数据和代码簿的说明在本书的笔记本中。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemResp.dct")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemResp.dat.gz") 

女性受访者文件代码簿位于ftp.cdc.gov/pub/Health_Statistics/NCHS/Dataset_Documentation/NSFG/Cycle6Codebook-Female.pdf。

nsfg.py模块提供了一个函数，用于读取女性受访者文件，清理一些变量，并返回一个DataFrame。

from  nsfg  import read_fem_respresp = read_fem_resp()
resp.shape 

(7643, 3092) 

这个DataFrame包含 3092 列，但我们只会使用其中的一小部分。

练习 2.1

我们将从totincr开始，它记录受访者的家庭总收入，用 1 到 14 的值进行编码。你可以阅读受访者文件的代码簿，以查看每个值代表什么收入水平。

创建一个表示这个变量分布的FreqTab对象，并将其绘制为条形图。

练习 2.2

制作一个记录每个受访者所生育子女数量的parity列频率表。你将如何描述这个分布的形状？

使用largest函数找到parity的最大值。你认为是否有任何值是错误的？

练习 2.3

让我们调查收入较高或较低的女性生育的孩子数量。使用查询方法选择收入最高的受访者（第 14 级）。仅绘制高收入受访者的parity频率表。

比较高收入受访者和其他人的平均parity。

计算这个差异的 Cohen 效应量。它与第一胎和其他人的怀孕长度差异相比如何？

这些结果是否表明收入较高的人生育的孩子较少，或者你能想到其他解释这个明显差异的原因吗？

Think Stats: Python 中的探索性数据分析，第 3 版

版权所有 2024 艾伦·B·唐尼

代码许可：MIT 许可

文本许可：Creative Commons Attribution-NonCommercial-ShareAlike 4.0 国际

概率质量函数

原文：allendowney.github.io/ThinkStats/chap03.html

在上一章中，我们使用FreqTab对象表示分布，它包含一组值及其频率——即每个值出现的次数。在这一章中，我们将介绍另一种描述分布的方法，即概率质量函数（PMF）。

为了表示概率质量函数（PMF），我们将使用一个名为Pmf的对象，它包含一组值及其概率。我们将使用Pmf对象来计算分布的均值和方差，以及偏度，这表明它是左偏还是右偏。最后，我们将探讨一个称为“检查悖论”的现象如何导致样本给出对分布的偏见观点。

在此处运行此笔记本。

from  os.path  import basename, existsdef  download(url):filename = basename(url)if not exists(filename):from  urllib.request  import urlretrievelocal, _ = urlretrieve(url, filename)print("Downloaded " + local)download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/thinkstats.py") 

try:import  empiricaldist
except ImportError:%pip install empiricaldist 

import  numpy  as  np
import  pandas  as  pd
import  matplotlib.pyplot  as  pltfrom  thinkstats  import decorate 

PMFs

Pmf对象就像一个包含概率而不是频率的FreqTab。因此，创建Pmf的一种方法是从一个FreqTab开始。例如，这里有一个FreqTab，它表示一个短序列中值的分布。

from  empiricaldist  import FreqTabftab = FreqTab.from_seq([1, 2, 2, 3, 5])
ftab 

	freqs
1	1
2	2
3	1
5	1

频率的总和是原始序列的大小。

n = ftab.sum()
n 

np.int64(5) 

如果我们将频率除以n，它们代表比例，而不是计数。

pmf = ftab / n
pmf 

	probs
1	0.2
2	0.4
3	0.2
5	0.2

这个结果表示序列中有 20%的值是 1，40%是 2，依此类推。

我们也可以从以下意义上将这些比例视为概率：如果我们从原始序列中随机选择一个值，选择值为 1 的概率是 0.2，选择值为 2 的概率是 0.4，依此类推。

因为我们将除以n，所以概率的总和是 1，这意味着这个分布是归一化的。

pmf.sum() 

np.float64(1.0) 

归一化的FreqTab对象代表一个概率质量函数（PMF），之所以称为“概率质量函数”，是因为与离散值相关的概率也称为“概率质量”。

empiricaldist库提供了一个表示概率质量函数的Pmf对象，因此，我们不需要创建FreqTab对象然后对其进行归一化，我们可以直接创建Pmf对象。

from  empiricaldist  import Pmfpmf = Pmf.from_seq([1, 2, 2, 3, 5])
pmf 

	probs
1	0.2
2	0.4
3	0.2
5	0.2

Pmf是归一化的，所以总概率是 1。

pmf.sum() 

np.float64(1.0) 

Pmf和FreqTab对象在许多方面都很相似。为了查找与一个值相关的概率，我们可以使用方括号操作符。

pmf[2] 

np.float64(0.4) 

或者使用括号来像函数一样调用Pmf。

pmf(2) 

np.float64(0.4) 

要给一个值分配一个概率，你必须使用方括号操作符。

pmf[2] = 0.2
pmf(2) 

np.float64(0.2) 

你可以通过增加与一个值相关的概率来修改现有的Pmf：

pmf[2] += 0.3
pmf[2] 

np.float64(0.5) 

或者你可以将一个概率乘以一个因子：

pmf[2] *= 0.5
pmf[2] 

np.float64(0.25) 

如果你修改一个Pmf，结果可能不会归一化——也就是说，概率可能不再加起来等于 1。

pmf.sum() 

np.float64(0.8500000000000001) 

normalize 方法通过除以总和来重新标准化 Pmf，并返回总和。

pmf.normalize() 

np.float64(0.8500000000000001) 

Pmf 对象提供了一个 copy 方法，这样你就可以创建并修改一个副本，而不会影响原始对象。

pmf.copy() 

	probs
1	0.235294
2	0.294118
3	0.235294
5	0.235294

与 FreqTab 对象一样，Pmf 对象有一个 qs 属性，可以访问数量，还有一个 ps 属性，可以访问概率。

它还有一个 bar 方法，可以将 Pmf 绘制为条形图，还有一个 plot 方法，可以将它绘制为折线图。

总结 PMF

在 第一章 中，我们通过将元素相加并除以元素数量来计算样本的平均值。这里有一个简单的例子。

seq = [1, 2, 2, 3, 5]n = len(seq)
mean = np.sum(seq) / n
mean 

np.float64(2.6) 

现在假设我们计算序列值的 PMF。

pmf = Pmf.from_seq(seq) 

给定 Pmf，我们仍然可以计算平均值，但过程不同——我们必须将概率和数量相乘，然后将乘积相加。

mean = np.sum(pmf.ps * pmf.qs)
mean 

np.float64(2.6) 

注意，我们不需要除以 n，因为我们已经在标准化 Pmf 时做了这件事。Pmf 对象有一个 mean 方法，做同样的事情。

pmf.mean() 

np.float64(2.6) 

给定一个 Pmf，我们可以通过计算每个数量与平均值的偏差来计算方差。

deviations = pmf.qs - mean 

然后我们将平方偏差乘以概率，并将乘积相加。

var = np.sum(pmf.ps * deviations**2)
var 

np.float64(1.84) 

var 方法做同样的事情。

pmf.var() 

np.float64(1.84) 

从方差中，我们可以按照通常的方式计算标准差。

np.sqrt(var) 

np.float64(1.3564659966250536) 

或者 std 方法做同样的事情。

pmf.std() 

np.float64(1.3564659966250536) 

Pmf 还提供了一个 mode 方法，可以找到概率最高的值。

pmf.mode() 

np.int64(2) 

随着我们继续前进，我们将看到更多的方法，但这些都足以开始。

班级规模悖论

作为使用 Pmf 对象可以做到的事情的例子，让我们考虑一个我称之为“班级规模悖论”的现象。

在许多美国大学和学院中，学生与教师的比例大约是 10:1。但学生常常惊讶地发现，他们中的许多课程的学生人数超过 10 人，有时甚至更多。这种差异有两个原因：

• 学生通常每学期选修 4 或 5 门课程，但教授通常只教 1 或 2 门。

• 小班的学生人数少，大班的学生人数多。

第一个影响是明显的，至少一旦指出；第二个影响更为微妙。让我们看一个例子。假设某大学在某个学期提供了 65 门课程，并且我们给出了以下大小范围内的课程数量。

ranges = pd.interval_range(start=5, end=50, freq=5, closed="left")
ranges.name = "class size"data = pd.DataFrame(index=ranges)
data["count"] = [8, 8, 14, 4, 6, 12, 8, 3, 2]
data 

	count
class size
---	---
[5, 10)	8
[10, 15)	8
[15, 20)	14
[20, 25)	4
[25, 30)	6
[30, 35)	12
[35, 40)	8
[40, 45)	3
[45, 50)	2

Pandas 函数 interval_range 创建一个 Index，其中每个标签代表一个值范围。符号 [5, 10) 表示 5 包含在区间内，而 10 不包含。由于我们不知道每个区间内班级的大小，让我们假设所有大小都是范围的中点。

sizes = ranges.left + 2
sizes 

Index([7, 12, 17, 22, 27, 32, 37, 42, 47], dtype='int64') 

现在让我们创建一个代表班级规模分布的Pmf。因为我们知道规模及其频率，我们可以直接创建一个Pmf，将计数、规模和名称作为参数传递。当我们对新的Pmf进行归一化时，结果是计数的总和。

counts = data["count"]
actual_pmf = Pmf(counts, sizes, name="actual")
actual_pmf.normalize() 

np.int64(65) 

如果你向大学询问平均班级规模，他们会报告这个分布的平均值，即 23.7。

actual_pmf.mean() 

np.float64(23.692307692307693) 

但如果你对一组学生进行调查，询问他们班级中有多少学生，并计算平均值，平均数会更大。让我们看看会大多少。

以下函数接受实际的班级规模Pmf，并创建一个新的Pmf，该Pmf代表学生观察到的班级规模。两个分布中的数量是相同的，但分布中的概率是乘以数量的，因为在规模为x的班级中，有x名学生观察到了这个班级。因此，观察到一个班级的概率与其规模成正比。

def  bias(pmf, name):# multiply each probability by class sizeps = pmf.ps * pmf.qs# make a new Pmf and normalize itnew_pmf = Pmf(ps, pmf.qs, name=name)new_pmf.normalize()return new_pmf 

现在我们可以计算学生观察到的有偏差的Pmf。

observed_pmf = bias(actual_pmf, name="observed") 

这里展示了这两个分布的形状。

from  thinkstats  import two_bar_plotstwo_bar_plots(actual_pmf, observed_pmf, width=2)
decorate(xlabel="Class size", ylabel="PMF") 

在观察到的分布中，小班级较少，大班级较多。有偏差的平均值是 29.1，比实际平均值高出近 25%。

observed_pmf.mean() 

np.float64(29.123376623376622) 

此操作也可以逆转。假设你想找到大学中班级规模分布，但你无法获得可靠的数据。一个选择是随机抽取学生样本，并询问他们班级中有多少学生。

由于我们刚刚看到的原因，结果会有偏差，但你可以用它来估计实际的分布。下面是一个通过除以大小来对Pmf进行去偏的函数。

def  unbias(pmf, name):# divide each probability by class sizeps = pmf.ps / pmf.qsnew_pmf = Pmf(ps, pmf.qs, name=name)new_pmf.normalize()return new_pmf 

下面是结果。

debiased_pmf = unbias(observed_pmf, "debiased")
debiased_pmf.mean() 

np.float64(23.692307692307693) 

去偏Pmf的均值与最初实际分布的均值相同。

如果你认为这个例子很有趣，你可能喜欢《很可能是在过度思考》的第二章，其中包含了这个以及其他几个被称为“检查悖论”的例子。

NSFG 数据

在上一章中，我们绘制了第一胎和其他婴儿的怀孕长度频率表。但由于组的大小不同，我们无法直接比较频率表。因为 PMFs 是归一化的，我们可以比较它们。所以让我们再次加载 NSFG 数据，并创建Pmf对象来表示怀孕长度的分布。

以下单元格下载数据文件并安装statadict，这是我们读取数据所需的。

try:import  statadict
except ImportError:%pip install statadict 

download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/nsfg.py")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dct")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dat.gz") 

nsfg模块提供了一个read_nsfg_groups函数，该函数读取数据，选择代表活产行的行，并将活产行分为第一胎和其他。它返回三个DataFrame对象。

from  nsfg  import get_nsfg_groupslive, firsts, others = get_nsfg_groups() 

我们可以使用firsts和others来为每个组中的怀孕长度创建一个Pmf。

first_pmf = Pmf.from_seq(firsts["prglngth"], name="firsts")
other_pmf = Pmf.from_seq(others["prglngth"], name="others") 

这里是第一胎和其他婴儿的 PMFs，以条形图的形式展示。

two_bar_plots(first_pmf, other_pmf)
decorate(xlabel="Weeks", ylabel="Probability", xlim=[20, 50]) 

通过绘制 PMF 而不是频率表，我们可以比较两个分布，而不会因为样本大小的差异而误导。根据此图，第一胎婴儿似乎比其他人更不可能准时到达（第 39 周）并且更有可能迟到（第 41 周和第 42 周）。

其他可视化

FreqTabograms 和 PMFs 在探索数据并尝试识别模式和关系时非常有用。一旦你有了想法，下一步就是设计一个可视化，尽可能清晰地展示你已识别的模式。

在 NSFG 数据中，分布的最大差异接近峰值。因此，在图表的这一部分进行放大，并选择 35 周到 46 周的数据是有意义的。

当我们像调用函数一样调用Pmf对象时，我们可以查找一系列数量并获取一系列概率。

weeks = range(35, 46)
first_pmf(weeks) 

array([0.03602991, 0.03897575, 0.04713347, 0.06163608, 0.4790392 ,0.12145932, 0.08157716, 0.04645366, 0.01971448, 0.00521187,0.00135962]) 

other_pmf(weeks) 

array([0.03210137, 0.03146779, 0.05216473, 0.07074974, 0.54466737,0.12249208, 0.04794087, 0.02597677, 0.01288279, 0.00485744,0.00084477]) 

因此，我们可以这样计算概率的差异。

diffs = first_pmf(weeks) - other_pmf(weeks)
diffs 

array([ 0.00392854,  0.00750796, -0.00503126, -0.00911366, -0.06562817,-0.00103276,  0.03363629,  0.02047689,  0.00683169,  0.00035443,0.00051485]) 

这是它们的模样，乘以 100 以表示百分比点的差异。

plt.bar(weeks, diffs * 100)
decorate(xlabel="Weeks", ylabel="Difference (percentage points)") 

此图使模式更清晰：第一胎婴儿在 39 周出生的可能性较小，而在 41 周和 42 周出生的可能性略大。

当我们在样本中看到这种模式时，我们不能确定它也在总体中成立——我们也不知道我们是否会在来自同一总体的另一个样本中看到它。我们将在第九章（chap09.html#chapter-hypothesis-testing）中重新审视这个问题。

术语表

本章的新术语不如前几章多。

• 归一化：如果一组概率的总和为 1，则称其为归一化。

• 概率质量函数（PMF）：一个将每个数量映射到其概率的函数。

练习

对于本章的练习，我们将使用 NSFG 受访者文件，其中每行代表一个受访者。下载数据的说明在本书的笔记本中。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemResp.dct")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemResp.dat.gz") 

此数据集的代码簿位于ftp.cdc.gov/pub/Health_Statistics/NCHS/Dataset_Documentation/NSFG/Cycle6Codebook-Female.pdf。

nsfg.py模块提供了一个函数，用于读取受访者文件并返回一个DataFrame。

from  nsfg  import read_fem_respresp = read_fem_resp()
resp.shape 

(7643, 3092) 

此DataFrame包含 7643 行和 3092 列。

练习 3.1

选择numbabes列，该列记录每个受访者的“出生存活婴儿数量”。创建一个FreqTab对象并显示该列值的频率。检查它们是否与代码簿中的频率一致。是否有任何特殊值应该替换为NaN？

然后创建一个Pmf对象并将其作为条形图绘制。分布是对称的，左偏斜还是右偏斜？

pmf.bar()
decorate(xlabel="Number of babies", ylabel="PMF") 

练习 3.2

与均值标识分布中的中心点，方差量化其分散程度一样，还有一个称为偏度的统计量，它表明分布是向左还是向右偏斜。

给定一个样本，我们可以通过计算立方偏差之和并除以标准差立方来计算偏度。例如，以下是计算numbabes偏度的方法。

numbabes = resp["numbabes"].replace(97, np.nan) 

deviations = numbabes - numbabes.mean()
skewness = np.mean(deviations**3) / numbabes.std(ddof=0) ** 3
skewness 

np.float64(1.7018914266755958) 

正值表示分布向右偏斜，负值表示它向左偏斜。

如果你给定的是一个Pmf而不是一系列值，你可以这样计算偏度：

1. 计算每个数量在Pmf中与均值的偏差。

2. 将偏差值立方，乘以Pmf中的概率，然后将乘积相加。

3. 将总和除以标准差立方。

编写一个名为pmf_skewness的函数，该函数接受一个Pmf对象并返回其偏度。

使用你的函数和numbabes的Pmf来计算偏度，并确认你得到我们上面计算出的相同结果。

练习 3.3

如果你调查儿童并询问他们家庭中有多少孩子，就会出现类似班级规模悖论的情况。有多个孩子的家庭更有可能出现在你的样本中，而没有孩子的家庭根本不可能出现在样本中。

从resp中选择numkdhh，它记录了每个受访者家庭中 18 岁以下孩子的数量。为这个列中的值创建一个Pmf。

使用bias函数计算如果我们调查儿童并询问他们家庭中有多少 18 岁以下的儿童（包括他们自己）时我们会看到的分布。

绘制实际和有偏分布，并计算它们的均值。

Think Stats: Exploratory Data Analysis in Python, 3rd Edition

版权所有 2024 艾伦·B·唐尼

代码许可：MIT License

文本许可：Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

累积分布函数

原文：allendowney.github.io/ThinkStats/chap04.html

频率表和概率质量函数是表示分布最熟悉的方式，但正如我们将在本章中看到的，它们有局限性。一个替代方案是累积分布函数（CDF），它对于计算百分位数很有用，尤其是在比较分布时特别有用。

在本章中，我们还将计算基于百分位数的统计量来量化分布的位置、范围和偏度。

在此处运行此笔记本。

from  os.path  import basename, existsdef  download(url):filename = basename(url)if not exists(filename):from  urllib.request  import urlretrievelocal, _ = urlretrieve(url, filename)print("Downloaded " + local)download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/thinkstats.py") 

try:import  empiricaldist
except ImportError:%pip install empiricaldist 

import  numpy  as  np
import  pandas  as  pd
import  matplotlib.pyplot  as  pltfrom  thinkstats  import decorate 

百分位数和百分位数排名

如果你参加过标准化考试，你可能会以原始分数和百分位数排名的形式收到你的成绩。在这个上下文中，百分位数排名是指得分与你相同或更低的人的百分比。所以如果你是“90 百分位数”，那么你的表现与参加考试的人中的 90%一样好或更好。

要理解百分位数和百分位数排名，让我们考虑一个基于跑步速度的例子。几年前我参加了马萨诸塞州的詹姆斯·乔伊斯漫步 10 公里公路赛跑。赛后，我下载了成绩来查看我的时间与其他跑者的比较。

下载数据的说明在本章的笔记本中。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/relay.py")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/Apr25_27thAn_set1.shtml"
) 

relay.py 模块提供了一个函数，用于读取结果并返回 Pandas DataFrame。

from  relay  import read_resultsresults = read_results()
results.head() 

	排名	分组/总数	分组	总时间	净时间	每英里分钟	英里每小时
0	1	1/362	M2039	30:43	30:42	4:57	12.121212
1	2	2/362	M2039	31:36	31:36	5:06	11.764706
2	3	3/362	M2039	31:42	31:42	5:07	11.726384
3	4	4/362	M2039	32:28	32:27	5:14	11.464968
4	5	5/362	M2039	32:52	32:52	5:18	11.320755

results 包含每个完成比赛的 1633 名跑者的一个行。我们将使用 MPH 这一列来量化性能，它包含每位跑者的平均速度（每小时英里）。我们将选择这一列并使用 values 提取速度作为一个 NumPy 数组。

speeds = results["MPH"].values 

我以 42:44 的成绩完成比赛，因此我们可以这样找到我的行。

my_result = results.query("Nettime == '42:44'")
my_result 

	排名	分组/总数	分组	总时间	净时间	每英里分钟	英里每小时
96	97	26/256	M4049	42:48	42:44	6:53	8.716707

我的行索引是 96，因此我们可以这样提取我的速度。

my_speed = speeds[96] 

我们可以使用 sum 来计算与我速度或更慢的跑者的数量。

(speeds <= my_speed).sum() 

np.int64(1537) 

我们可以使用 mean 来计算与我速度或更慢的跑者的百分比。

(speeds <= my_speed).mean() * 100 

np.float64(94.12124923453766) 

结果是我的场内百分位数排名，大约是 94%。

更普遍地，以下函数计算序列中特定值的百分位数排名。

def  percentile_rank(x, seq):"""Percentile rank of x.x: valueseq: sequence of valuesreturns: percentile rank 0-100"""return (seq <= x).mean() * 100 

在 results 中，Division 列表示每个跑者所在的分组，通过性别和年龄范围来识别——例如，我所在的分组是 M4049，包括 40 至 49 岁的男性跑者。我们可以使用 query 方法来选择我所在分组的行并提取他们的速度。

my_division = results.query("Division == 'M4049'")
my_division_speeds = my_division["MPH"].values 

现在我们可以使用 percentile_rank 来计算我在我所在分组的百分位数排名。

percentile_rank(my_speed, my_division_speeds) 

np.float64(90.234375) 

反过来，如果我们给定一个百分位数排名，以下函数会找到序列中相应的值。

def  percentile(p, seq):n = len(seq)i = (1 - p / 100) * (n + 1)return seq[round(i)] 

n 是序列中的元素数量；i 是具有给定百分位数排名的元素的索引。当我们查找百分位数排名时，相应的值被称为 百分位数。

percentile(90, my_division_speeds) 

np.float64(8.591885441527447) 

在我的分组中，第 90 个百分位数大约是 8.6 英里每小时。

现在，在那场比赛多年之后，我进入了 M5059 分组。那么让我们看看我需要跑多快才能在新分组中获得相同的百分位数排名。我们可以通过将我在 M4049 分组中的百分位数排名，大约 90.2%，转换为 M5059 分组中的速度来回答这个问题。

next_division = results.query("Division == 'M5059'")
next_division_speeds = next_division["MPH"].valuespercentile(90.2, next_division_speeds) 

np.float64(8.017817371937639) 

与我具有相同百分位数排名的 M5059 分组的人跑了略超过 8 英里每小时。我们可以使用 query 来找到他。

next_division.query("MPH > 8.01").tail(1) 

	Place	Div/Tot	Division	Guntime	Nettime	Min/Mile	MPH
222	223	18/171	M5059	46:30	46:25	7:29	8.017817

他以 46:25 的成绩完成比赛，在他所在分组中排名第 18，共有 171 人。

通过对百分位数和百分位数的介绍，我们为累积分布函数做好了准备。

CDFs

累积分布函数，或称为 CDF，是描述一组值分布的另一种方式，与频率表或 PMF 一起。给定一个值 x，CDF 计算小于或等于 x 的值的比例。作为一个例子，我们将从一个简短的序列开始。

t = [1, 2, 2, 3, 5] 

计算 CDF 的一种方法是从 PMF 开始。这里有一个 Pmf 对象，它代表了 t 中值的分布。

from  empiricaldist  import Pmfpmf = Pmf.from_seq(t)
pmf 

	probs
1	0.2
2	0.4
3	0.2
5	0.2

如前一章所述，我们可以使用方括号操作符在 Pmf 中查找值。

pmf[2] 

np.float64(0.4) 

结果是序列中等于给定值的值的比例。在这个例子中，五个值中有两个等于 2，所以结果是 0.4。我们也可以将这个比例视为从序列中随机选择一个值等于 2 的概率。

Pmf 有一个 make_cdf 方法，它计算 Pmf 中概率的累积和。

cdf = pmf.make_cdf()
cdf 

	probs
1	0.2
2	0.6
3	0.8
5	1.0

结果是一个 Cdf 对象，它是一种 Pandas Series 类型。我们可以使用方括号操作符来查找值。

cdf[2] 

np.float64(0.6000000000000001) 

结果是序列中小于或等于给定值的值的比例。在这个例子中，序列中有三个值小于或等于2，所以结果是 0.6。我们也可以将这个比例视为从序列中随机选择一个值小于或等于2的概率。

我们可以使用括号来像函数一样调用Cdf对象。

cdf(3) 

array(0.8) 

累积分布函数对所有数字都定义，而不仅仅是序列中出现的数字。

cdf(4) 

array(0.8) 

为了可视化累积分布函数（Cdf），我们可以使用step方法，该方法将累积分布函数绘制为阶梯函数。

cdf.step()
decorate(xlabel="x", ylabel="CDF") 

作为第二个例子，让我们创建一个表示上一节中跑步速度分布的累积分布函数（Cdf）。Cdf类提供了一个from_seq函数，我们可以用它从序列中创建一个Cdf对象。

from  empiricaldist  import Cdfcdf_speeds = Cdf.from_seq(speeds) 

下面是这个样子——垂直线是我的速度。

cdf_speeds.step()
plt.axvline(my_speed, ls=":", color="gray")
decorate(xlabel="Speed (mph)", ylabel="CDF") 

如果我们查看我的速度，结果是与我速度或更慢的跑者的比例。如果我们乘以 100，我们得到我的百分位数排名。

cdf_speeds(my_speed) * 100 

np.float64(94.12124923453766) 

因此，这是思考累积分布函数的一种方式——给定一个值，它计算类似于百分位数排名的东西，只不过它是在 0 到 1 之间的比例，而不是 0 到 100 之间的百分比。

Cdf提供了一个inverse方法，它计算累积分布函数的逆——给定 0 到 1 之间的比例，它找到相应的值。

例如，如果有人说他们跑得和 50%的选手一样快或更快，我们可以这样找到他们的速度。

cdf_speeds.inverse(0.5) 

array(6.70391061) 

如果你有一个比例，并且使用逆 CDF 找到相应的值，结果被称为分位数——因此逆 CDF 有时被称为分位数函数。

如果你有一个分位数，并且使用 CDF 找到相应的比例，结果实际上并没有一个名字，奇怪的是。为了与百分位数和百分位数排名保持一致，它可以被称为“分位数排名”，但据我所知，没有人这样称呼它。最常见的是，它只是被称为“累积概率”。

比较累积分布函数

累积分布函数在比较分布时特别有用。作为一个例子，让我们比较第一个孩子和其他人的出生体重的分布。我们将再次加载 NSFG 数据集，并将其分为三个DataFrames：所有活产婴儿、第一个孩子和其他人。

下面的单元格下载数据文件并安装statadict，这是我们读取数据所需的。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/nsfg.py")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dct")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dat.gz") 

try:import  statadict
except ImportError:%pip install statadict 

from  nsfg  import get_nsfg_groupslive, firsts, others = get_nsfg_groups() 

从firsts和others中，我们将选择总出生体重（以磅为单位），使用dropna删除nan值。

first_weights = firsts["totalwgt_lb"].dropna()
first_weights.mean() 

np.float64(7.201094430437772) 

other_weights = others["totalwgt_lb"].dropna()
other_weights.mean() 

np.float64(7.325855614973262) 

看起来，第一个孩子平均体重略轻。但是，有几种方式会导致这种差异——例如，可能有一小部分第一个孩子特别轻，或者一小部分其他孩子特别重。在这些情况下，分布会有不同的形状。另一种可能性是，分布可能有相同的形状，但位置不同。

要比较分布，我们可以尝试绘制概率质量函数（PMF）。

from  empiricaldist  import Pmffirst_pmf = Pmf.from_seq(first_weights, name="first")
other_pmf = Pmf.from_seq(other_weights, name="other") 

但正如我们可以在以下图中看到的，它并不奏效。

from  thinkstats  import two_bar_plotstwo_bar_plots(first_pmf, other_pmf, width=0.06)
decorate(xlabel="Weight (pounds)", ylabel="PMF") 

我调整了柱子的宽度和透明度，以尽可能清楚地显示分布，但很难进行比较。有许多峰和谷，以及一些明显的变化，但很难判断哪些特征是有意义的。此外，很难看到整体模式；例如，从视觉上不明显哪个分布的均值更高。

这些问题可以通过数据分箱来解决——也就是说，将数量的范围分成不重叠的区间，并计算每个箱中数量的数量。分箱可能很有用，但很难正确确定箱的大小。如果它们足够大以平滑噪声，它们也可能平滑掉有用的信息。

一个好的替代方案是绘制累积分布函数（CDF）。

first_cdf = first_pmf.make_cdf()
other_cdf = other_pmf.make_cdf() 

它们看起来是这样的。

first_cdf.plot(ls="--")
other_cdf.plot(alpha=0.5)
decorate(xlabel="Weight (pounds)", ylabel="CDF") 

这张图使得分布的形状和它们之间的差异更加清晰。第一个孩子的曲线始终位于其他人的曲线左侧，这表明第一个孩子在整个分布中略轻——在中间点以上的差异更大。

百分位统计量

在第三章中，我们计算了算术平均值，它确定了分布中的中心点，以及标准差，它量化了分布的分散程度。在之前的练习中，我们计算了偏度，它表明分布是左偏还是右偏。所有这些统计量的一个缺点是它们对异常值敏感。数据集中单个极端值可以对均值、标准差和偏度产生重大影响。

一个替代方案是使用基于分布百分位的统计量，这些统计量通常更稳健，这意味着它们对异常值不太敏感。为了演示，让我们再次加载 NSFG 数据，但不进行任何数据清理。

from  nsfg  import read_statadct_file = "2002FemPreg.dct"
dat_file = "2002FemPreg.dat.gz"preg = read_stata(dct_file, dat_file) 

记住，出生体重记录在两列中，一列是磅，另一列是盎司。

birthwgt_lb = preg["birthwgt_lb"]
birthwgt_oz = preg["birthwgt_oz"] 

如果我们用birthwgt_oz中的值创建一个Hist对象，我们可以看到它包括特殊的值 97、98 和 99，这些值表示缺失数据。

from  empiricaldist  import HistHist.from_seq(birthwgt_oz).tail(5) 

	频率
birthwgt_oz
---	---
14.0	475
15.0	378
97.0	1
98.0	1
99.0	46

birthwgt_lb列包括相同的特殊值；它还包括值 51，这肯定是一个错误。

Hist.from_seq(birthwgt_lb).tail(5) 

	freqs
birthwgt_lb
---	---
15.0	1
51.0	1
97.0	1
98.0	1
99.0	57

现在让我们想象两种情况。在一种情况下，我们通过用nan替换缺失和无效值来清理这些变量，然后计算以磅为单位的总重量。将birthwgt_oz_clean除以 16 将其转换为十进制磅。

birthwgt_lb_clean = birthwgt_lb.replace([51, 97, 98, 99], np.nan)
birthwgt_oz_clean = birthwgt_oz.replace([97, 98, 99], np.nan)total_weight_clean = birthwgt_lb_clean + birthwgt_oz_clean / 16 

在另一种情况下，我们忽略了数据清理，并意外地用这些假数据计算了总重量。

total_weight_bogus = birthwgt_lb + birthwgt_oz / 16 

假数据集只包含 49 个假数据，大约是数据的 0.5%。

count1, count2 = total_weight_bogus.count(), total_weight_clean.count()
diff = count1 - count2diff, diff / count2 * 100 

(np.int64(49), np.float64(0.5421553441026776)) 

现在让我们计算两种情况下的数据平均值。

mean1, mean2 = total_weight_bogus.mean(), total_weight_clean.mean()
mean1, mean2 

(np.float64(7.319680587652691), np.float64(7.265628457623368)) 

假数据对平均值有适度的影响。如果我们认为清洁数据的平均值是正确的，那么假数据的平均值偏差不到 1%。

(mean1 - mean2) / mean2 * 100 

np.float64(0.74394294099376) 

这样的错误可能不会被察觉——但现在让我们看看标准差会发生什么变化。

std1, std2 = total_weight_bogus.std(), total_weight_clean.std()
std1, std2 

(np.float64(2.096001779161835), np.float64(1.4082934455690173)) 

(std1 - std2) / std2 * 100 

np.float64(48.83274403900607) 

假数据的标准差偏差接近 50%，因此这更加明显。最后，这是两个数据集的偏度。

def  skewness(seq):"""Compute the skewness of a sequenceseq: sequence of numbersreturns: float skewness"""deviations = seq - seq.mean()return np.mean(deviations**3) / seq.std(ddof=0) ** 3 

skew1, skew2 = skewness(total_weight_bogus), skewness(total_weight_clean)
skew1, skew2 

(np.float64(22.251846195422484), np.float64(-0.5895062687577697)) 

# how much is skew1 off by?
(skew1 - skew2) / skew2 

np.float64(-38.74658112171128) 

假数据集的偏度偏差接近 40 倍，并且符号错误！添加了异常值后，分布强烈向右偏斜，正如大正偏度所指示的。但有效数据的分布略微向左偏斜，正如小负偏度所指示的。

这些结果表明，少数异常值对平均值有适度的影响，对标准差有强烈的影响，对偏度有灾难性的影响。

另一种方法是使用基于百分比的统计。具体来说：

• 中位数，即第 50 百分位数，像平均值一样，在分布中确定一个中心点。

• 四分位距，即第 25 百分位数和第 75 百分位数之间的差值，量化了分布的分散程度，就像标准差一样。

• 四分位距偏度使用分布的四分位数（第 25、第 50 和第 75 百分位数）来量化偏度。

Cdf对象提供了一种高效的方法来计算这些基于百分位的统计量。为了演示，让我们从假数据和清洁数据集中创建一个Cdf对象。

cdf_total_weight_bogus = Cdf.from_seq(total_weight_bogus)
cdf_total_weight_clean = Cdf.from_seq(total_weight_clean) 

以下函数接受一个Cdf并使用其inverse方法来计算第 50 百分位数，即中位数（至少，这是定义数据集中位数的一种方法）。

def  median(cdf):m = cdf.inverse(0.5)return m 

现在我们可以计算两个数据集的中位数。

median(cdf_total_weight_bogus), median(cdf_total_weight_clean) 

(array(7.375), array(7.375)) 

结果是相同的，因此在这种情况下，异常值对中位数没有任何影响。一般来说，异常值对中位数的影响小于对平均值的影响。

四分位距（IQR）是第 75 百分位数和第 25 百分位数之间的差值。以下函数接受一个Cdf并返回 IQR。

def  iqr(cdf):low, high = cdf.inverse([0.25, 0.75])return high - low 

下面是两个数据集的四分位距。

iqr(cdf_total_weight_bogus), iqr(cdf_total_weight_clean) 

(np.float64(1.625), np.float64(1.625)) 

通常，异常值对四分位距的影响小于对标准差的影响——在这种情况下，它们完全没有影响。

最后，这是一个计算四分位数偏斜的函数，它依赖于三个统计数据：

• 中位数，

• 第 25 百分位数和第 75 百分位数的中间值，

• 半四分位距，即四分位距的一半。

def  quartile_skewness(cdf):low, median, high = cdf.inverse([0.25, 0.5, 0.75])midpoint = (high + low) / 2semi_iqr = (high - low) / 2return (midpoint - median) / semi_iqr 

这里是两个数据集的四分位数偏斜。

qskew1 = quartile_skewness(cdf_total_weight_bogus)
qskew2 = quartile_skewness(cdf_total_weight_clean)
qskew1, qskew2 

(np.float64(-0.07692307692307693), np.float64(-0.07692307692307693)) 

这些例子中的少量异常值对这些四分位数偏斜没有影响。这些例子表明，基于百分位数的统计对异常值和数据中的错误不太敏感。

随机数

Cdf 对象提供了一种高效的方法来从分布中生成随机数。首先，我们从 0 到 1 之间的均匀分布中生成随机数。然后，在这些点上评估累积分布函数（CDF）的逆函数。以下函数实现了此算法。

def  sample_from_cdf(cdf, n):ps = np.random.random(size=n)return cdf.inverse(ps) 

为了演示，让我们生成一个跑步速度的随机样本。

sample = sample_from_cdf(cdf_speeds, 1001) 

为了确认它是否有效，我们可以比较样本和原始数据集的累积分布函数（CDF）。

cdf_sample = Cdf.from_seq(sample)cdf_speeds.plot(label="original", ls="--")
cdf_sample.plot(label="sample", alpha=0.5)decorate(xlabel="Speed (mph)", ylabel="CDF") 

样本遵循原始数据的分布。为了理解此算法的工作原理，考虑以下问题：假设我们从跑步速度的总体中随机选择一个样本，并查找样本中速度的百分位数排名。现在假设我们计算百分位数排名的累积分布函数（CDF）。你认为它会是什么样子？

让我们来探究一下。以下是我们在生成的样本中的百分位数排名。

percentile_ranks = cdf_speeds(sample) * 100 

这里是百分位数排名的累积分布函数（CDF）。

cdf_percentile_rank = Cdf.from_seq(percentile_ranks)
cdf_percentile_rank.plot()decorate(xlabel="Percentile rank", ylabel="CDF") 

百分位数排名的累积分布函数（CDF）在 0 到 1 之间接近一条直线。这是有道理的，因为在任何分布中，百分位数排名小于 50%的比例是 0.5；百分位数排名小于 90%的比例是 0.9，依此类推。

Cdf 提供了一个使用此算法的 sample 方法，因此我们也可以生成这样的样本。

sample = cdf_speeds.sample(1001) 

术语表

• 百分位数排名：在分布中小于或等于给定数量的值的百分比。

• 百分位数：与给定百分位数排名相关联的分布中的值。

• 累积分布函数（CDF）：一个将值映射到小于或等于该值的分布比例的函数。

• 分位数：在分布中大于或等于给定比例的值的值。

• 稳健：如果一个统计量受极端值或异常值的影响较小，则称其为稳健。

• 四分位距（IQR）：第 75 百分位数和第 25 百分位数之间的差值，用于衡量分布的离散程度。

练习

练习 4.1

你出生时体重是多少？如果你不知道，打电话给你的母亲或其他人。如果没有人知道，你可以使用我的出生体重，8.5 磅，来进行这个练习。

使用 NSFG 数据（所有活产婴儿），计算出生体重的分布，并使用它来找到你的百分位数排名。如果你是第一个孩子，找到第一个孩子分布中的百分位数排名。否则使用其他人的分布。如果你在 90 个百分位数或更高，给你的母亲打个电话道歉。

from  nsfg  import get_nsfg_groupslive, firsts, others = get_nsfg_groups() 

练习 4.2

在 NSFG 数据集中，对于活产婴儿，babysex 列表示婴儿是男性还是女性。我们可以使用 query 来选择男性和女性婴儿的行。

male = live.query("babysex == 1")
female = live.query("babysex == 2")
len(male), len(female) 

(4641, 4500) 

创建代表男性和女性婴儿出生体重分布的 Cdf 对象。绘制这两个 CDF。这两个分布的形状和位置有什么不同？

如果一个男婴体重为 8.5 磅，他的百分位数排名是多少？具有相同百分位数排名的女婴体重是多少？

练习 4.3

从 NSFG 数据集的怀孕数据中，选择 agepreg 列，创建一个 Cdf 来表示每次怀孕的怀孕年龄分布。使用 CDF 来计算年龄小于或等于 20 岁的百分比，以及年龄小于或等于 30 岁的百分比。使用这些结果来计算 20 到 30 岁之间的百分比。

from  nsfg  import read_fem_pregpreg = read_fem_preg() 

练习 4.4

这里是之前在本章中提到的詹姆斯·乔伊斯漫步完成者的跑步速度。

speeds = results["MPH"].values 

创建一个代表这些速度分布的 Cdf，并使用它来计算中位数、四分位数范围和四分位数偏斜。分布是向左还是向右偏斜？

练习 4.5

np.random.random 生成的数字应该在 0 和 1 之间是均匀的，这意味着样本的 CDF 应该是一条直线。让我们看看这是否正确。这里有一个 1001 个数字的样本。绘制这个样本的 CDF。它看起来像一条直线吗？

t = np.random.random(1001) 

Think Stats: 使用 Python 进行探索性数据分析，第 3 版

版权 2024 艾伦·B·唐尼

代码许可：MIT 许可证

文本许可：Creative Commons 知识共享署名-非商业性使用-相同方式共享 4.0 国际