数据科学学习指南（一）

原文：zh.annas-archive.org/md5/9950f82165bee426119f99afc3ab612d

译者：飞龙

协议：CC BY-NC-SA 4.0

序言

数据科学是令人兴奋的工作。从混乱的数据中获得洞察力对于业务、医学、政策等各种决策具有重要价值。本书《学习数据科学》旨在为读者做好数据科学的准备。为此，我们设计了以下特别功能：

焦点放在基础上

技术变迁无常。尽管本书中使用特定技术，但我们的目标是为读者提供数据科学的基础构建块。我们通过揭示如何思考数据科学问题和挑战，以及涵盖各个技术背后的基础知识来实现这一点。我们的目标是即使技术变化，也能为读者提供帮助。

覆盖整个数据科学生命周期

我们不仅仅关注单一主题，比如如何处理数据表或如何应用机器学习技术，而是覆盖整个数据科学生命周期——提出问题、获取数据、理解数据和理解世界的过程。完成整个生命周期通常是数据科学家工作中最难的部分。

使用真实数据

为了准备处理实际问题，我们认为从使用真实数据的例子中学习至关重要，这些数据有其缺陷和优点。我们在选择本书中呈现的数据集时，精心挑选了对实际数据分析产生影响的数据，而不是使用过于精炼或合成的数据。

通过案例研究应用概念

本书中我们包含了一些扩展案例研究，这些案例研究追随或扩展了其他数据科学家的分析。这些案例研究向读者展示了如何在实际环境中应用数据科学生命周期。

结合计算和推断思维

在工作中，数据科学家需要预见他们编写代码时的决策如何以及数据集的大小可能会如何影响统计分析。为了为读者未来的工作做好准备，《学习数据科学》整合了计算和统计思维。我们还通过模拟研究而非数学证明来激励统计概念。

本书的文本和代码都是开源的，在 GitHub 上可获取。

预期的背景知识

我们期望读者精通 Python，并了解如何使用内置的数据结构如列表、字典和集合；导入和使用其他包中的函数和类；并且能够从头开始编写函数。我们还使用numpy Python 包而不进行介绍，但不要求读者具有太多使用经验。

如果读者对概率、微积分和线性代数有一些了解，他们将从本书中受益更多，但我们的目标是直观地解释数学思想。

本书的组织结构

本书共有 21 章，分为六个部分：

第一部分（章节 1–5）

第一部分 描述了生命周期是什么，以基本水平通过了整个生命周期，并介绍了本书中使用的术语。该部分以一个关于公交到站时间的短期案例研究结束。

第二部分（第 6 – 7 章）

第二部分 介绍数据框架和关系以及如何编写代码使用 pandas 和 SQL 操作数据。

第三部分（第 8 – 12 章）

第三部分 关注数据获取，发现其特征和发现问题。在理解这些概念后，读者可以获取一个数据文件，并向他人描述数据集的有趣特征。本部分以一个关于空气质量的案例研究结束。

第四部分（第 13 – 14 章）

第四部分 探讨了广泛使用的替代数据源，如文本，二进制和来自网络的数据。

第五部分（第 15 – 18 章）

第五部分 专注于使用数据理解世界。除了模型拟合，特征工程和模型选择之外，它还涵盖了推断主题，如置信区间和假设检验。本部分以一个关于为肯尼亚兽医预测驴重的案例研究结束。

第六部分（第 19 – 21 章）

第六部分 完成了我们对使用逻辑回归和优化进行监督学习的研究。它以一个案例研究结束，预测新闻文章是否做出真实或虚假声明。

本书末尾，我们提供了有关本书引入的许多主题的更多学习资源，并提供了使用的完整数据集列表。

本书中使用的约定

本书使用以下排版约定：

斜体

表示新术语，URL，电子邮件地址，文件名和文件扩展名。

等宽字体

用于程序清单，以及段落内引用程序元素，如变量或函数名，数据库，数据类型，环境变量，语句和关键字。

等宽粗体

显示用户应按字面输入的命令或其他文本。

等宽斜体

显示应用程序中应替换为用户提供值或根据上下文确定值的文本。

注意

此元素表示一般注释。

提示

此元素表示提示。

警告

此元素表示警告或注意。

使用代码示例

补充材料（代码示例，练习等）可在 learningds.org 下载。

如果您有技术问题或在使用代码示例时遇到问题，请发送电子邮件至 bookquestions@oreilly.com。

本书旨在帮助您完成工作。一般情况下，如果本书提供示例代码，您可以在自己的程序和文档中使用它。除非您复制了大量代码片段，否则无需征得我们的许可。例如，编写使用本书多个代码片段的程序不需要许可。售卖或分发来自 O’Reilly 书籍的示例需要许可。引用本书并引用示例代码回答问题不需要许可。将本书大量示例代码整合到您产品的文档中需要许可。

我们感谢您的署名。署名通常包括标题、作者、出版商和 ISBN。例如：“学习数据科学 by Sam Lau, Joseph Gonzalez, and Deborah Nolan (O’Reilly). Copyright 2023 Sam Lau, Joseph Gonzalez, and Deborah Nolan, 978-1-098-11300-1.”

如果您认为使用代码示例超出了公平使用范围或上述权限，请随时通过 bookquestions@oreilly.com 与我们联系。

O’Reilly 在线学习

注意

40 多年来，O’Reilly Media 提供技术和商业培训，为企业的成功提供知识和见解。

我们独特的专家和创新者网络通过书籍、文章和我们的在线学习平台分享他们的知识和专业知识。O’Reilly 的在线学习平台为您提供按需访问的现场培训课程、深入学习路径、交互式编码环境以及来自 O’Reilly 和 200 多个其他出版商的大量文本和视频。欲了解更多信息，请访问https://oreilly.com.

如何联系我们

请将有关本书的评论和问题发送给出版商：

• O’Reilly Media, Inc.

• 1005 Gravenstein Highway North

• Sebastopol, CA 95472

• 800-889-8969（美国或加拿大）

• 707-829-7019（国际或本地）

• 707-829-0104（传真）

• support@oreilly.com

• https://www.oreilly.com/about/contact.html

我们为本书设立了网页，列出勘误、示例和任何额外信息。您可以访问https://oreil.ly/learning-data-science获取这些信息。

获取有关我们的书籍和课程的新闻和信息，请访问https://oreilly.com.

在 LinkedIn 上找到我们：https://linkedin.com/company/oreilly-media.

在 Twitter 上关注我们：https://twitter.com/oreillymedia.

在 YouTube 上观看我们：https://youtube.com/oreillymedia.

致谢

这本书源于我们共同设计和教授的“数据科学原理与技术”课程，这是加州大学伯克利分校的本科课程。我们首次在 2017 年春季教授了“数据 100”，以响应学生对第二门数据科学课程的需求；他们希望这门课程能为他们在数据科学的深造和职场生涯做好准备。

自那次首次开设以来，我们教授的成千上万名学生对我们是一种激励。我们还受益于与其他讲师的合作教学，包括 Ani Adhikari、Andrew Bray、John DeNero、Sandrine Dudoit、Will Fithian、Joe Hellerstein、Josh Hug、Anthony Joseph、Scott Lee、Fernando Perez、Alvin Wan、Lisa Yan 和 Bin Yu。我们特别感谢 Joe Hellerstein 在数据整理方面的见解，Fernando Perez 鼓励我们包含像 NetCDF 这样的更复杂的数据结构，Josh Hug 提出了 PurpleAir 案例研究的想法，Duncan Temple Lang 在课程早期版本上的合作。我们还感谢伯克利的学生担任助教，并特别提到那些为本书的早期版本做出贡献的人员：Ananth Agarwal、Ashley Chien、Andrew Do、Tiffany Jann、Sona Jeswani、Andrew Kim、Jun Seo Park、Allen Shen、Katherine Yen 和 Daniel Zhu。

本书的核心部分是我们整理和分析的许多数据集，我们非常感谢那些使他们的数据向我们开放和可用的个人和组织。在本书的结尾，我们列出这些贡献者以及原始数据来源，以及相关的研究论文、博客文章和报告。

最后，我们感谢 O’Reilly 团队的工作，将这本书从课堂笔记转化为出版物，特别是 Melissa Potter、Jess Haberman、Aaron Black、Danny Elfanbaum 和 Mike Loukides。我们还要感谢技术审阅人员的评论，这些评论改进了本书：Sona Jeswani、Thomas Nield、Siddharth Yadav 和 Abhijit Dasgupta。

第一部分：数据科学生命周期

第一章：数据科学生命周期

数据科学是一个快速发展的领域。在撰写本文时，人们仍在努力确定数据科学究竟是什么，数据科学家做什么，以及数据科学家应该具备哪些技能。然而，我们知道的是，数据科学利用统计学和计算机科学的方法和原则结合在一起，从数据中获取洞见。学习计算机科学和统计学的结合使我们成为更好的数据科学家。我们还知道，我们获取的任何洞见都需要在我们正在解决的问题的背景下进行解释。

本书涵盖了数据科学家需要帮助做出各种重要决策的基本原理和技能。凭借技术技能和概念理解，我们可以处理以数据为中心的问题，例如评估疫苗的有效性、自动过滤假新闻、校准空气质量传感器，并为分析师在政策变更方面提供建议。

为了帮助您掌握更大的局面，我们围绕一种称为数据科学生命周期的工作流程组织了主题。在本章中，我们介绍这个生命周期。与其他数据科学书籍不同，这些书籍倾向于专注于生命周期的某一部分，或者只涉及计算或统计主题，我们涵盖了从头到尾的整个周期，并考虑了统计和计算两个方面。

生命周期的阶段

图 1-1 显示了数据科学生命周期，分为四个阶段：提出问题、获取数据、理解数据和理解世界。我们特意使这些阶段变得广泛。根据我们的经验，生命周期的机制经常变化。计算机科学家和统计学家继续开发用于处理数据的新软件包和编程语言，并开发更专业的新方法。

图 1-1. 数据科学生命周期的四个高级阶段，箭头表示这些阶段如何相互关联

尽管发生了这些变化，我们发现几乎每个数据项目都包括以下四个阶段：

提出问题

提出好问题是数据科学的核心，识别不同类型的问题指导我们的分析。我们涵盖四类问题：描述性、探索性、推论性和预测性。例如，“房价随时间的变化是如何的？”是描述性的，而“房屋的哪些方面与售价相关？”是探索性的。将一个广泛的问题缩小到可以用数据回答的问题是生命周期中这个第一阶段的关键元素。这可能涉及咨询参与研究的人员、弄清如何测量某事以及设计数据收集协议。清晰而专注的研究问题帮助我们确定我们需要的数据、要查找的模式以及如何解释结果。它还可以帮助我们完善问题、识别所提出的问题类型，并计划生命周期的数据收集阶段。

获取数据

当数据昂贵且难以收集时，当我们的目标是从数据中推广到世界时，我们的目标是为收集数据定义精确的协议。其他时候，数据是廉价且易于获取的。这在在线数据源中尤其如此。例如，Twitter让人们可以快速下载数百万个数据点。当数据丰富时，我们可以通过获取和探索数据，然后进一步明确研究问题来开始分析。在这两种情况下，大多数数据都有缺失或异常值以及其他我们需要考虑的异常情况。无论数据来自何处，我们都需要检查数据质量。同样重要的是考虑数据的范围；例如，我们确定数据的代表性以及在收集过程中可能存在的潜在偏差来源。这些考虑有助于我们确定我们对发现的信任程度。通常，我们必须在更正式地分析数据之前操作数据。我们可能需要修改结构、清理数据值并转换测量值以准备分析。

理解数据

在获取和准备数据后，我们希望仔细检查它们，而探索性数据分析通常是关键。在我们的探索中，我们制作图表以揭示有趣的模式并以视觉方式总结数据。我们还继续寻找数据中的问题。当我们搜索模式和趋势时，我们使用汇总统计和构建统计模型，如线性和逻辑回归。根据我们的经验，这个生命周期的阶段是高度迭代的。理解数据还可能使我们回到数据科学生命周期的早期阶段。我们可能发现需要修改或重新进行数据清洗和操作，获取更多数据以补充我们的分析，或者鉴于数据的限制，重新定义我们的研究问题。我们在这个阶段进行的描述性和探索性分析可能已经足够回答我们的问题，或者我们可能需要进入下一个阶段以便对我们的数据进行推广。

理解世界

当我们的目标纯粹是描述性或探索性时，分析就止步于生命周期的“理解数据”阶段。在其他时候，我们的目标是量化我们发现的趋势在我们数据之外的泛化程度。我们可能希望使用我们拟合到数据的模型来从样本推断到总体，以便对世界进行推断或预测未来的观察。为了从样本到总体进行推断，我们使用统计技术如 A/B 测试和置信区间。为了对未来的观察进行预测，我们创建预测区间并使用数据的训练-测试分割。

对于生命周期的每个阶段，我们解释理论概念，介绍数据技术和统计方法，并展示它们在实际示例中的运用。在整个过程中，我们依赖其他数据科学家的真实数据和分析，而不是虚构的数据，因此您可以学习如何进行自己的数据获取、清洗、探索和正式分析，并得出合理的结论。本书的每一章通常都集中于数据科学生命周期的一个阶段，但我们还包括案例研究章节，展示整个生命周期的实际应用。

注意

理解探索、推断、预测和因果之间的差异可能是一个挑战。我们很容易在数据中找到的相关性与因果关系混淆。例如，探索性或推断性分析可能会针对“暴露于空气污染较多的人是否具有更高的肺病发病率？”这样的问题寻找相关性，而因果性问题可能会问：“给维基百科贡献者颁发奖是否会增加生产力？”除非我们进行了随机实验（或近似实验），通常我们无法回答因果性问题。我们在整本书中都强调这些重要的区别。

生命周期的示例

本书中放置了几个案例研究，涵盖整个数据科学生命周期。这些案例研究有双重作用。它们专注于生命周期中的一个阶段，以提供本书部分中的具体示例，并演示整个循环。

第五章的重点是问题的相互作用以及数据如何用来回答问题。简单的问题“为什么我的公交总是迟到？”提供了一个丰富的案例研究，对于初学者来说足够基础，可以追踪生命周期的各个阶段，但又足够微妙，可以展示我们如何应用统计和计算思维来回答问题。在这个案例研究中，我们进行了模拟研究，以了解乘客等待时间的分布。我们还拟合了一个简单的模型来总结等待时间的统计数据。这个案例研究还展示了作为数据科学家，您如何收集自己的数据来回答您感兴趣的问题。

第十二章研究了在美国各地使用的大众市场空气传感器的准确性。我们设计了一种利用环保局维护的高精度传感器数据来改进廉价传感器读数的方法。这个案例研究展示了如何通过来自严格维护和精确的政府监测设备的数据来改进众包开放数据。在这个过程中，我们专注于清理和合并来自多个来源的数据，但我们也拟合模型来调整和改进空气质量测量。

在第十八章中，我们关注模型构建和预测。但我们涵盖了整个生命周期，并看到问题的兴趣如何影响我们构建的模型。我们的目标是帮助肯尼亚农村的兽医，他们无法访问用于称重驴的称重器，以便为生病的动物开药。当我们学习研究设计，清理数据，并在简单性与准确性之间取得平衡时，我们评估了我们模型的预测能力，并展示了科学家如何与面临实际问题的人合作，并协助他们找到解决方案。

最后，在第二十一章中，我们试图通过算法区分手工分类的新闻故事中的假新闻和真实新闻。在这个案例研究中，我们再次看到，可读性强的信息为数据科学家开发新技术并研究当今重要问题创造了惊人的机会。这些数据来自于网上新闻报道，并由阅读这些报道的人员将其分类为假新闻或真实新闻。我们还看到，数据科学家如何以创造性的方式将一般信息（如新闻文章的内容）转化为可分析的数据，以解决当前的热点问题。

概要

数据科学生命周期为本书提供了一个组织结构。在处理来自科学、医学、政治、社交媒体和政府等多个来源的数据集时，我们始终牢记这一生命周期。第一次使用数据集时，我们提供了数据收集背景的上下文，对研究数据感兴趣的问题以及理解数据所需的描述。通过这种方式，我们旨在贯彻整本书的良好数据科学实践。

生命周期的第一阶段——提出问题——通常被视为需要应用技术以获得数字的问题，比如“这个 A/B 测试的p-值是多少？”或者在实践中经常见到的模糊问题，比如“我们能恢复美国梦吗？”回答第一类问题几乎没有练习开发研究问题的实践经验。而回答第二类问题则很难在没有指导的情况下进行，指导需要将兴趣的一般领域转化为可以用数据回答的问题。提出问题与理解数据局限性以回答问题之间的相互作用是下一章的主题。

第二章：问题与数据范围

作为数据科学家，我们使用数据来回答问题，数据收集过程的质量可以显著影响数据的有效性和准确性，我们从分析中得出的结论的强度以及我们做出的决策。在本章中，我们描述了理解数据收集和评估数据在解决感兴趣问题的有用性方面的一般方法。理想情况下，我们希望数据能够代表我们研究的现象，无论该现象是人群特征、物理模型还是某种社会行为。通常情况下，我们的数据并不包含完整信息（范围在某种程度上受限），但我们希望使用数据准确描述人口、估算科学数量、推断特征之间的关系形式或预测未来结果。在所有这些情况下，如果我们的数据不能代表我们研究对象，那么我们的结论可能会受到限制，可能具有误导性，甚至是错误的。

为了激励对这些问题的思考，我们首先以大数据的力量和可能出现的问题为例子。然后，我们提供一个框架，可以帮助您将研究目标（您的问题）与数据收集过程联系起来。我们称之为数据范围^(1)，并提供术语来帮助描述数据范围，以及来自调查、政府数据、科学仪器和在线资源的示例。在本章后期，我们考虑数据准确性的含义。在那里，我们介绍不同形式的偏差和变异，并描述可能发生这些情况的情形。整个过程中，例子涵盖了您作为数据科学家可能使用的数据种类的光谱；这些例子来自科学、政治、公共卫生和在线社区。

大数据和新机遇

公开可用数据的巨大增加为数据科学创造了新的角色和机会。例如，数据记者在数据中寻找有趣的故事，就像传统报道记者寻找新闻故事一样。数据记者的生命周期始于寻找可能有趣故事的现有数据，而不是从研究问题开始，然后找出如何收集新数据或使用现有数据来解决问题。

公民科学项目是另一个例子。它们吸引许多人（和仪器）参与数据收集。这些数据集体上供研究人员组织项目使用，并且通常在公共存储库中供大众进一步调查使用。

行政和组织数据的可用性创造了其他机会。研究人员可以将从科学研究中收集的数据与例如医疗数据进行链接；这些行政数据之所以收集，并不是因为直接源于感兴趣的问题，但它们在其他情境中可能很有用。这种链接可以帮助数据科学家扩展其分析的可能性，并交叉检查其数据的质量。此外，找到的数据可以包括数字痕迹，例如您的网络浏览活动、社交媒体上的发布以及您的在线朋友和熟人网络，它们可能非常复杂。

当我们拥有大量的行政数据或广泛的数字痕迹时，很容易将它们视为比传统的小型研究更具决定性的数据。我们甚至可能认为这些大型数据集可以取代科学研究，本质上是人口普查。这种过度自信被称为“大数据傲慢”。具有广泛范围的数据并不意味着我们可以忽略数据代表性的基础问题，也不能忽略测量、依赖性和可靠性问题。（甚至很容易仅仅因为巧合而发现毫无意义或荒谬的关系。）一个众所周知的例子就是谷歌流感趋势追踪系统。

示例：谷歌流感趋势

数字流行病学是流行病学的一个新分支，利用公共卫生系统之外产生的数据来研究人群中的疾病模式和健康动态。谷歌流感趋势（Google Flu Trends，GFT）追踪系统是数字流行病学的早期示例之一。2007 年，研究人员发现，统计人们搜索与流感相关术语的次数可以准确估计流感病例数。这一明显的成功引起了轰动，许多研究人员对大数据的可能性感到兴奋。然而，GFT 未能达到预期，并在 2015 年被放弃。

发生了什么问题？毕竟，GFT 利用了来自在线查询的数百万数字痕迹来预测流感活动。尽管最初取得了成功，在 2011-2012 年流感季节中，谷歌的数据科学家们发现，GFT 并不能替代由美国疾病控制与预防中心（CDC）从全国各地实验室收集的三周前统计数据。相比之下，GFT 高估了 CDC 的数据，有 108 周中的 100 周。周复一周，尽管基于大数据，GFT 对流感病例的估计仍然过高：

在这幅图中，从第 412 周到第 519 周，GFT（实线）高估了实际的 CDC 报告（虚线）100 倍。还在此处绘制了基于三周前 CDC 数据和季节性趋势的模型预测（点线），这比 GFT 更接近实际情况。

数据科学家发现，一个简单的模型，基于过去的 CDC 报告，使用三周前的 CDC 数据和季节性趋势，比 GFT 在预测流感流行方面做得更好。GFT 忽视了可以通过基本统计方法提取的大量信息。这并不意味着从在线活动中获取的大数据是无用的。事实上，研究人员已经表明，GFT 数据与 CDC 数据的结合可以显著改善 GFT 预测和基于 CDC 的模型。通常情况下，结合不同的方法会比单一方法带来改进。

GFT 的例子告诉我们，即使我们拥有大量信息，数据与提出问题之间的关联至关重要。理解这一框架可以帮助我们避免回答错误的问题，将不适当的方法应用于数据，并夸大我们的发现。

注意

在大数据时代，我们往往会被诱惑收集更多数据来精确回答问题。毕竟，人口普查给了我们完美的信息，那么大数据不应该几乎是完美的吗？不幸的是，这通常并非如此，特别是在管理数据和数字痕迹方面。你想研究的人群的一小部分无法接触（参见第三章中的 2016 年选举反转）或者测量过程本身（就像这个 GFT 的例子）可能会导致预测不准确。重要的是考虑数据的范围，因为它与正在调查的问题相关联。

需要牢记的一个关键因素是数据的范围。范围包括考虑我们想要研究的人群，如何获取有关该人群的信息，以及我们实际测量的内容。深思熟虑这些要点可以帮助我们发现方法上的潜在漏洞。我们将在下一节中详细探讨这一点。

目标人群、访问框架和样本

在数据生命周期的重要初始步骤是将感兴趣的问题表达在主题领域的背景下，并考虑问题与收集的数据之间的联系。在甚至考虑分析或建模步骤之前这样做是一个很好的做法，因为它可能会揭示问题和数据之间的不一致。作为将数据收集过程与调查主题联系起来的一部分，我们识别人群、访问人群的手段、测量工具和收集过程中使用的附加协议。这些概念——目标人群、访问框架和样本——帮助我们理解数据的范围，无论我们的目标是获取关于人群、科学数量、物理模型、社会行为或其他内容的知识。

目标人群

目标人群包括构成你最终意图描述和得出结论的人群的元素集合。元素可以是人群中的一个人，选民中的一名选民，一系列推文中的一条推文，或者一个州中的一个县。我们有时将一个元素称为单位或原子。

访问框架

访问框架是可供测量和观察的元素集合。这些单位是你可以研究目标人群的对象。理想情况下，访问框架和人群完全对齐，意味着它们由完全相同的元素组成。然而，访问框架中的单位可能只是目标人群的一个子集；此外，框架可能包括不属于人群的单位。例如，为了了解选民在选举中的投票意向，你可能会通过电话联系人们。你打电话的人可能不是选民，因此他们在你的框架中但不在人群中。另一方面，从未接听未知号码的选民无法联系到，因此他们在人群中但不在你的框架中。

样本

样本是从访问框架中取出的单位子集，用于观察和测量。样本为你提供了分析数据，从而进行关于感兴趣人群的预测或概括的基础。当资源被用于跟进非响应者和追踪难以找到的单位时，一个小样本比一个大样本或试图对被忽视人群子集进行人口普查更为有效。

访问框架的内容与目标人群相比，以及从框架中选择单位到样本中的方法，是确定数据是否可以被视为目标人群代表性的重要因素。如果访问框架不代表目标人群，那么样本数据很可能也不具有代表性。如果单位的抽样方式存在偏差，也会引发代表性问题。

你还需要考虑数据范围中的时间和地点。例如，在某个疾病流行的地区测试药物试验的有效性可能与在感染率较低的世界其他地区进行的试验不相称（见 第三章）。此外，为了研究随时间变化的数据收集，如大气中二氧化碳（CO[2]）的月度测量（见 第九章）和预测流感趋势的谷歌搜索的每周报告，具有时间结构，我们在检查数据时需要注意。在其他时候，数据可能存在空间模式。例如，本节稍后描述的环境健康数据是针对加利福尼亚州每个人口普查区报告的，我们可能会制作地图以寻找空间相关性。

如果你没有收集数据，你需要考虑是谁收集了数据以及出于什么目的。现在更多的数据是被动收集而不是出于特定目的而收集的，这一点尤为重要。仔细审视发现的数据，并问自己这些数据是否可以用来解决你的问题，可以避免进行无效的分析或得出不恰当的结论。

对于以下小节中的示例，我们从一个一般性问题开始，将其缩小为可以用数据回答的问题，并在此过程中确定目标人群、访问框架和样本。这些概念在图表中用圆圈和矩形表示，这些形状的重叠配置有助于揭示范围的关键方面。在每个示例中，我们还描述了数据范围的相关时间和空间特征。

例如：什么使在线社区的成员活跃？

Wikipedia 上的内容是由属于 Wikipedia 社区的志愿者编写和编辑的。这个在线社区对 Wikipedia 的成功和活力至关重要。为了理解如何激励在线社区的成员，研究人员进行了一项关于 Wikipedia 贡献者的实验。一个更具体的问题是：奖励是否会增加 Wikipedia 贡献者的活跃度？在这个实验中，目标人群是前一个月内对 Wikipedia 贡献最活跃的顶尖贡献者——即前 1% 的最活跃贡献者。接入框架排除了那些在该月已经接受过奖励（奖励）的人群。接入框架故意排除了人群中的一些贡献者，因为研究人员想要衡量奖励的影响，而那些已经接受了一种奖励的人可能会有不同的行为（见 图 2-1）。

图 2-1. Wikipedia 实验中范围的表现

样本是从框架中随机选择的 200 名贡献者。观察了这些贡献者 90 天，并收集了他们在 Wikipedia 上的数字活动迹象。请注意，贡献者群体并非静态；有定期的更替。在研究开始前的一个月内，超过 144,000 名志愿者为 Wikipedia 创作内容。从这群顶尖贡献者中选择限制了研究结果的普遍性，但鉴于顶尖贡献者群体的规模，如果能通过非正式奖励影响他们以维持或增加其贡献，这仍然是一个有价值的发现。

在许多实验和研究中，我们无法包括所有的人口单位在框架内。通常情况下，接入框架包括愿意参与研究/实验的志愿者。

例子：谁会赢得选举？

2016 年美国总统选举的结果让许多人和许多民意调查员感到意外。大多数选举前的民意调查预测希拉里·克林顿会击败唐纳德·特朗普。政治民意调查是一种试图在选举之前衡量人们会投票给谁的公众意见调查。由于意见随时间变化，焦点缩小到了一个“赛马”问题，即受访者被问及如果明天举行选举，他们会投票给哪位候选人：候选人 A 还是候选人 B。

民意调查定期在总统竞选期间进行，随着选举日的临近，我们预计民意调查会更好地预测结果，因为偏好会稳定下来。民意调查通常也会在州内进行，然后再结合起来预测总体赢家。因此，民意调查的时间和地点都很重要。民意调查机构也很重要，有些机构一直比其他机构更准确。

在这些选举前调查中，目标人群包括将在选举中投票的人，例如此例中的 2016 年美国总统选举。然而，调查员只能猜测某人是否会在选举中投票，因此访问框架包括那些被认为可能是选民的人（通常基于过去的投票记录，但也可能使用其他因素）。由于人们通过电话联系，因此访问框架仅限于拥有座机或移动电话的人群。样本由框架中按随机拨号方案选择的人组成（见 图 2-2）。

图 2-2. 2016 年总统选举调查中范围的表示

在 第三章 中，我们讨论了人们不愿接听电话或参与民意调查对选举预测的影响。

示例：环境危害如何影响个体健康？

为了解决这个问题，加利福尼亚环境保护局（CalEPA）、加利福尼亚州环境健康危害评估办公室（OEHHA）和公众共同开发了 CalEnviroScreen 项目。该项目利用从美国人口普查的人口统计数据、加利福尼亚州卫生保健访问和信息部门的健康统计数据，以及加利福尼亚州空气资源局维护的空气监测站收集的污染测量数据，研究加利福尼亚社区中人口健康与环境污染的关系。

理想情况下，我们希望研究加利福尼亚州的人群，并评估这些环境危害对个体健康的影响。然而，在这种情况下，数据只能在人口普查区一级别获取。访问框架由居住在同一人口普查区的群体组成。因此，框架中的单位是人口普查区，样本是一个普查，即州内所有普查区，因为提供了所有普查区的数据（见 图 2-3）。

图 2-3. CalEnviroScreen 项目的范围；访问框架中的网格代表人口普查区

不幸的是，我们无法分解一个普查区内的信息以研究个人。这种聚合影响了我们可以探讨的问题以及我们可以得出的结论。例如，我们可以询问加利福尼亚社区因哮喘住院率与空气质量之间的关系。但是我们无法回答最初提出的有关个体健康的问题。

这些示例展示了目标、访问框架和样本的可能配置。当一个框架无法覆盖所有人时，我们应考虑这些缺失信息可能如何影响我们的发现。同样地，我们要思考当一个框架包括非人口中的成员时可能会发生什么。此外，抽样技术可以影响样本对总体的代表性。当你考虑推广你的数据发现时，你还需要考虑收集数据所用仪器和程序的质量。如果你的样本是一次完全符合目标的普查，但信息收集不当，那么你的发现将毫无价值。这是下一节的主题。

仪器和协议

当我们考虑数据的范围时，也要考虑用于进行测量和测量过程的仪器和协议，我们称之为协议。对于调查而言，仪器通常是样本中的个体填写的问卷。调查的协议包括如何选择样本，如何跟进非回应者，访谈者培训，保护机密性等。

良好的仪器和协议对所有类型的数据收集都很重要。如果我们想测量自然现象，如大气中的二氧化碳，我们需要量化仪器的准确性。校准仪器和进行测量的协议对于获取准确的测量结果至关重要。仪器可能会失校，测量结果可能随时间漂移，导致测量极不准确。

在实验中，协议也至关重要。理想情况下，任何可能影响实验结果的因素都应受到控制。例如，温度、时间、医疗记录的机密性，甚至测量顺序都需要保持一致，以排除这些因素可能带来的影响。

随着数字痕迹的增加，支持在线活动的算法是动态的，并不断进行重新设计。例如，谷歌的搜索算法不断调整以提高用户服务和广告收入。搜索算法的变化可以影响从搜索中生成的数据，进而影响基于这些数据构建的系统，例如谷歌流感趋势跟踪系统。这种变化环境可能使得维护数据收集协议变得不可行，并且难以复制研究结果。

许多数据科学项目涉及将来自多个来源的数据进行关联。每个来源应通过这种数据范围构建来检查，考虑到各个来源的差异。此外，用于合并多个来源数据的匹配算法需要被清楚理解，以便比较来源的人口和框架。

用于研究自然现象的仪器测得的测量可以在目标、访问框架和样本的范围图中进行描述。这种方法有助于理解仪器的准确性。

测量自然现象

引入的用于观察目标群体的范围图可以扩展到我们想要测量数量的情况，比如空气中的粒子计数、化石的年龄或光速。在这些情况下，我们将要测量的数量视为一个未知的确切值。（这个未知值通常被称为参数。）我们可以根据这种情况调整我们的范围图：我们将目标缩小到代表未知的一个点；仪器的准确性作为访问框架；样本由仪器测得的测量组成。你可以把框架想象成一个飞镖靶，仪器是扔飞镖的人，而飞镖落在靶心周围。飞镖的散布对应仪器测量的结果。目标点不被飞镖投手看到，但理想情况下它与靶心重合。

为了说明测量误差及其与抽样误差的关系，我们考虑测量空气中 CO[2] 浓度的问题。

例子：空气中 CO[2] 的水平是多少？

CO[2] 是全球变暖的重要信号，因为它在地球大气中捕获热量。没有 CO[2]，地球将会异常寒冷，但这是一个微妙的平衡。CO[2] 浓度的增加推动全球变暖，并威胁到我们地球的气候。为了解决这个问题，自 1958 年以来在毛纳罗亚观测站对 CO[2] 浓度进行了监测。这些数据为理解全球变暖的威胁提供了关键的基准。

在考虑数据的范围时，我们要考虑数据收集的地点和时间。科学家选择在毛纳罗亚火山测量 CO[2]，因为他们希望在空气中测量 CO[2] 的背景水平。毛纳罗亚位于太平洋中，远离污染源，观测站位于一个被裸露的火山岩石包围的山顶上，远离可以从空气中去除 CO[2] 的植物。

测量 CO[2] 的仪器尽可能准确非常重要。严格的协议被制定用来保持仪器处于最佳状态。例如，毛纳罗亚定期使用不同类型的设备对空气样本进行测量，并将其他样本送往实验室进行更精确的测量。这些测量有助于确定仪器的准确性。此外，每小时对一个参考气体进行 5 分钟的测量，每天对另外两个参考气体进行 15 分钟的测量。这些参考气体具有已知的 CO[2] 浓度。将测量浓度与已知值进行比较有助于识别仪器中的偏差。

尽管毛纳洛亚岛上背景空气中的二氧化碳相对稳定，但在任何小时内测量的五分钟平均浓度与小时平均值会有偏差。这些偏差反映了仪器的精确度和气流的变化。

数据收集的范围可以总结如下：在毛纳洛亚岛上空的特定位置，在特定的一个小时内，存在真实的二氧化碳浓度背景；这是我们的目标（参见图 2-4）。仪器进行测量并报告五分钟平均值。这些读数构成一个样本，包含在访问帧中，即飞镖板。如果仪器工作正常，飞镖靶心与目标（一小时平均二氧化碳浓度）重合，测量值集中在靶心周围，偏差约为 0.30 百万分之一（ppm）。二氧化碳的测量是每百万份干燥空气中二氧化碳分子的数量，因此测量单位是 ppm。

图 2-4. 访问帧代表仪器的精确度；星星代表感兴趣的真实值

我们在下一节继续使用飞镖板的类比，介绍偏差和变异的概念，描述样本可能不代表总体的常见方式，并将精度与协议联系起来。

精确度

在人口普查中，访问帧与人口匹配，样本捕捉整个人口。在这种情况下，如果我们进行了设计良好的问卷调查，那么我们对人口有完整准确的信息，范围是完美的。类似地，在测量大气中的二氧化碳浓度时，如果我们的仪器具有完美的精确度并正确使用，那么我们可以测量二氧化碳浓度的确切值（忽略空气波动）。这些情况是罕见的，如果不是不可能的话。在大多数情况下，我们需要量化测量的准确度，以便将我们的发现推广到未被观察到的领域。例如，我们经常使用样本来估计人群的平均值，从测量中推断科学未知值的值，或预测新个体的行为。在每种情况下，我们还希望有一个可量化的精度度量。我们想知道我们的估计、推断和预测与真实情况有多接近。

早些时候引入的飞镖投掷到飞镖板的类比在理解精度方面很有帮助。我们将精度分为两个基本部分：偏差和精确度（也称为变异性）。我们的目标是让飞镖击中飞镖板上的靶心，并使靶心与看不见的目标对齐。飞镖在板上的分布代表我们测量中的精度，从靶心到我们正在瞄准的未知值的差距则代表偏差。

图 2-5 显示了低偏差和高偏差以及精度的组合。在这些图中，点代表采取的测量，星代表真实的、未知的参数值。点在访问框架内形成一个散射，在这个散射中，点代表采取的测量，星代表真实的、未知的参数值。当访问框架的靶心大致位于星星上方时（顶行），测量点在感兴趣的值周围散布，偏差较小。较大的飞镖板（右侧）表示测量中的更广泛的分布（较低的精度）。

图 2-5. 低和高测量偏差和精度的组合

代表性数据将我们置于图表的顶行，其中偏差小，意味着未知的目标与靶心对齐。理想情况下，我们的仪器和协议将我们置于图表的左上部，那里的变化也很小。底行中点的模式系统地错过了目标值。增加样本量不会纠正这种偏差。

偏差类型

偏差有多种形式。我们在这里描述了一些经典类型，并将它们与我们的目标-访问-样本框架联系起来：

覆盖偏差

当访问框架不包括目标人群中的所有人时发生。例如，基于电话调查的调查无法接触到没有电话的人。在这种情况下，不能接触到的人可能在重要方面与访问框架中的人不同。

选择偏差

当用于选择样本单位的机制倾向于选择某些单位的频率超过它们应该被选择的频率时，就会出现覆盖偏差。例如，方便采样选择最容易获得的单位。当那些容易接触到的人与那些难以接触到的人在重要方面有所不同时，问题就会出现。另一个例子是，观察研究和实验通常依赖于志愿者（选择参与的人），而这种自我选择可能会导致偏差，如果志愿者与目标人群在重要方面不同。

非响应偏差

有两种形式：单位和项目。当被选样本中的某人不愿意参与时发生单位无响应（他们可能永远不会接听来自陌生人的电话）。当某人接听电话但拒绝回答特定的调查问题时发生项目无响应。如果选择不回答的人与选择回答的人有系统性差异，那么非响应可能会导致偏差。

测量偏差

当一个仪器系统地错过了一个方向的目标时发生。例如，低湿度可能会导致我们对空气污染的测量错误地偏高。此外，测量设备随时间变化可能变得不稳定并漂移，从而产生系统误差。在调查中，当问题措辞不清或引导性，或者受访者可能不愿诚实回答时，可能会出现测量偏差。

这些偏见中的每一种都可能导致数据未集中在未知的目标值上。通常情况下，我们无法评估偏见的潜在幅度，因为关于那些不在访问范围之外、不太可能被选中作为样本或不愿意回应的人几乎没有信息可用。协议是减少这些偏见来源的关键。从框架中选择样本或将单位分配给实验条件的机会性机制可以消除选择偏见。非响应跟进协议可以鼓励参与，从而减少非响应偏见。试点调查可以改善问题措辞，从而减少测量偏见。校准仪器的程序和随机顺序进行测量的协议可以减少测量偏见。

在 2016 年美国总统选举中，非响应偏见和测量偏见是预测胜者不准确的关键因素。几乎所有选民在选举前预测希拉里将击败特朗普。特朗普意外获胜令人惊讶。选举后，许多民意调查专家试图诊断民意调查中的问题所在。美国公共舆论研究协会发现预测有两个关键原因出现了偏差：

• 大学受教育选民被过度代表。大学受教育选民比受教育较少的人更有可能参与调查，并且在 2016 年，他们更有可能支持希拉里。更高的受教育选民的响应率使样本存在偏见，并且高估了对希拉里的支持。

• 选民在选举前几天临时决定或改变他们的偏好。由于民意调查是静态的，只能直接测量当前的信仰，它无法反映态度的变化。

很难弄清楚人们是忍住了他们的偏好还是改变了他们的偏好，以及这种偏见有多大。然而，选后民调帮助民意调查专家了解事后发生了什么。它们表明，在密集竞争的州（例如密歇根州），许多选民在竞选的最后一周作出了选择，而且这群人大多支持特朗普。

并非在所有情况下都需要避免偏见。如果仪器精度高（低方差）且偏见小，则该仪器可能优于另一种具有更高方差但无偏见的仪器。例如，偏倚研究可能对试点调查仪器或捕获大型研究设计中的有用信息有用。许多时候，我们最多只能招募研究志愿者。鉴于这种限制，仍然有可能将这些志愿者纳入研究，并使用随机分配将其分为治疗组。这就是随机对照实验背后的理念。

无论是否存在偏差，数据通常也表现出变异。可以通过使用偶然机制选择样本来有意引入变异，也可以通过仪器的精确性自然发生。在下一节中，我们将识别三种常见的变异来源。

变异类型

以下类型的变异由偶然机制产生，并具有可量化的优势：

抽样变异

使用偶然选择样本的结果。在这种情况下，原则上我们可以计算特定元素集被选为样本的机会。

分配变异

在控制实验中发生，当我们将单位随机分配到处理组时。在这种情况下，如果我们以不同方式分割单位，那么我们可以从实验中获得不同的结果。这种分配过程使我们能够计算特定组分配的机会。

测量误差

测量过程的结果。如果用于测量的仪器没有漂移或偏差，并且具有可靠的误差分布，那么当我们对同一对象进行多次测量时，我们会得到围绕真实值中心的随机测量变化。

瓮模型 是一个简单的抽象，对理解变异很有帮助。这个模型建立了一个容器（一个瓮，类似于花瓶或桶），里面装满了标记过的相同弹珠，我们使用从瓮中抽取弹珠的简单动作来推理取样方案、随机对照实验和测量误差。对于这些变异类型，瓮模型帮助我们使用概率或模拟估计变异的大小（参见第三章）。选择维基百科贡献者接受非正式奖励的例子提供了瓮模型的两个示例。

回顾维基百科实验，从 1440 位顶级贡献者中随机选择了 200 位贡献者。然后，这 200 位贡献者再次随机分成两组，每组 100 人。一组获得非正式奖励，另一组没有。我们使用瓮模型来描述这个选择和分组过程：

1. 想象一个装满 1440 颗相同形状和大小的弹珠的瓮，每颗弹珠上写着 1440 个 Wikipedia 用户名之一。（这是访问框架。）

2. 仔细将弹珠在瓮中混合，选择一个弹珠并将其放在一边。

3. 重复混合并选择弹珠以获得 200 个弹珠。

从样本中抽取的弹珠。接下来，为了确定哪 200 位贡献者将获得奖励，我们使用另一个瓮：

1. 在第二个瓮中，放入前述样本中的 200 个弹珠。

2. 仔细混合这些弹珠，选择一个弹珠并将其放在一边。

3. 重复，选择 100 个弹珠。也就是说，一次选择一个弹珠，在中间混合，并将选择的弹珠放在一边。

从 100 颗抽取的弹珠被分配到治疗组，并对应于获奖的贡献者。留在瓮中的 100 颗弹珠形成对照组，他们不会获得奖励。

无论是样本选择还是奖励接受者的选择都使用了机会机制。如果我们再次重复第一次取样活动，将所有 1,440 颗弹珠放回原来的瓮中，那么我们很可能会得到一个不同的样本。这种变化是抽样变异的来源。同样地，如果我们再次重复随机分配过程（保持 200 个样本不变），那么我们会得到一个不同的治疗组。分配变异是由这第二个机会过程引起的。

维基百科实验提供了抽样和分配变异的例子。在这两种情况下，研究人员对数据收集过程施加了机会机制。测量误差有时也可以被视为遵循瓮模型的机会过程。例如，我们可以用这种方式表征马瓦努阿洛阿的 CO[2]监测仪的测量误差。

如果我们可以在数据的变化和瓮模型之间进行准确的类比，瓮模型为我们提供了估算变化规模的工具（见第三章）。这是非常理想的，因为我们可以为数据的变化提供具体的值。然而，确认瓮模型是否合理地描述了变化来源至关重要。否则，我们的准确性声明可能存在严重缺陷。我们需要尽可能了解数据范围的所有内容，包括数据收集中使用的仪器、协议和机会机制，以应用这些瓮模型。

总结

无论您处理的数据类型是什么，在进行清理、探索和分析之前，请花点时间了解数据的来源。如果您没有收集数据，请问自己：

• 谁收集了数据？

• 为什么要收集数据？

这些问题的答案可以帮助确定这些找到的数据是否可以用来解决您感兴趣的问题。

考虑数据的范围。关于数据收集的时间和空间方面的问题可以提供宝贵的见解：

• 数据是何时收集的？

• 数据是在哪里收集的？

这些问题的答案帮助您确定您的发现是否与您感兴趣的情况相关，或者您的情况可能与其他地方和时间不可比较。

对于范围概念的核心是回答以下问题：

• 目标人群是什么（或未知参数值是什么）？

• 目标是如何访问的？

• 选择样本/进行测量的方法是什么？

• 使用了哪些仪器，它们是如何校准的？

尽可能回答这些问题可以为您提供宝贵的见解，了解您对发现的信任程度以及是否可以从中进行概括。

本章为您提供了术语和思考并回答这些问题的框架。该章节还概述了如何识别可能影响您发现准确性的偏差和方差的方法。为了帮助您理解偏差和方差，我们介绍了以下图表和概念：

• 范围图示，显示目标人群、访问框架和样本之间的重叠。

• 飞镖板，描述仪器的偏差和方差。

• 在使用概率机制从访问框架中选择样本、将群体分为实验处理组或从良好校准的仪器中进行测量时，使用乌尔恩模型。

这些图表和模型试图概括理解如何识别数据限制并评估数据在回答问题中的有用性所需的关键概念。第三章继续发展乌尔恩模型，更正式地量化准确性并设计模拟研究。

^(1) “范围”概念改编自约瑟夫·赫勒斯坦（Joseph Hellerstein）的课程笔记，涉及范围、时间性和忠实度。

第三章：模拟与数据设计

在本章中，我们开发了理解数据抽样及其对偏差和方差影响的基础理论。我们建立这个基础不是基于经典统计的干涩方程，而是基于一个装满大理石的乌尔恩的故事。我们使用模拟的计算工具来推理从乌尔恩中选择大理石的属性以及它们对现实世界数据收集的启示。我们将模拟过程与常见的统计分布联系起来（干涩的方程），但是模拟的基本工具使我们能够超越仅仅使用方程直接建模的范围。

例如，我们研究了民意调查员未能预测 2016 年美国总统选举结果的失败。我们的模拟研究使用了宾夕法尼亚州实际投票情况。我们模拟了这六百万选民民意调查的抽样变异，揭示了回应偏差如何扭曲民意调查，并看到仅仅收集更多数据是无济于事的。

在第二个模拟研究中，我们研究了一项控制实验，证明了一种 COVID-19 疫苗的有效性，但也引发了关于疫苗相对有效性的激烈争论。将实验抽象为一个乌尔恩模型为我们提供了一个工具，用于研究随机对照实验中的分配变化。通过模拟，我们找到了临床试验的预期结果。我们的模拟，连同对数据范围的仔细检查，驳斥了疫苗无效的说法。

第三个例子使用模拟来模拟测量过程。当我们将我们人工测量的空气质量的波动与真实测量进行比较时，我们可以评估乌尔恩模型模拟测量空气质量波动的适当性。这种比较为我们校准紫外线空气质量监测器提供了背景，使它们可以更准确地在低湿度时期（如火灾季节）测量空气质量。

然而，在我们解决我们时代一些最重要的数据争论之前，我们首先从一个非常小的故事开始，这个故事是关于几颗大理石坐在一个乌尔恩中的故事。

乌尔恩模型

乌尔恩模型是由雅各布·伯努利在 18 世纪初开发的，作为模拟从人群中选择项目的一种方式。图 3-1 显示的乌尔恩模型提供了一个随机从乌尔恩中取样大理石过程的视觉描述。乌尔恩最初有五颗大理石：三颗黑色和两颗白色。图表显示进行了两次抽取：首先抽出一颗白色大理石，然后抽出一颗黑色大理石。

图 3-1. 两颗弃不重复从乌尔恩中抽出的大理石的图表

要建立一个乌尔恩模型，我们首先需要做出几个决定：

• 乌尔恩中的大理石数量

• 每颗大理石的颜色（或标签）

• 从乌尔恩中抽出的大理石的数量

最后，我们还需要决定抽样过程。对于我们的过程，我们将弹珠混合在容器中，当我们选择一个弹珠进行样本时，我们可以选择记录其颜色并将其放回容器（有放回抽样），或者将其设置为不能再次被选中（无放回抽样）。

这些决策构成了我们模型的参数。我们可以通过选择这些参数来调整弹珠模型，以描述许多现实世界的情况。举例来说，考虑 图 3-1 中的示例。我们可以使用numpy的random.choice方法在两次取样之间模拟从我们的容器中取出两个弹珠（无放回）。numpy库支持数组的函数，对于数据科学特别有用：

`import` `numpy` `as` `np`

`urn` `=` `[``"``b``"``,` `"``b``"``,` `"``b``"``,` `"``w``"``,` `"``w``"``]`
`print``(``"``Sample 1:``"``,` `np``.``random``.``choice``(``urn``,` `size``=``2``,` `replace``=``False``)``)`
`print``(``"``Sample 2:``"``,` `np``.``random``.``choice``(``urn``,` `size``=``2``,` `replace``=``False``)``)`

Sample 1: ['b' 'w']
Sample 2: ['w' 'b']

注意我们将replace参数设置为False，以表明一旦我们取出一个弹珠，就不会放回到容器中。

有了这个基本设置，我们可以大致回答关于我们期望看到的样本种类的问题。我们的样本中仅包含一种颜色的弹珠的机会是多少？如果我们在选择后将每个弹珠放回，机会会改变吗？如果我们改变了容器中的弹珠数量会怎样？如果我们从容器中抽取更多的弹珠会怎样？如果我们重复这个过程很多次会发生什么？

这些问题的答案对我们理解数据收集至关重要。我们可以基于这些基本技能来模拟弹珠，并将模拟技术应用于经典概率方程难以解决的现实问题。

例如，我们可以使用模拟来轻松估计我们抽取的两个弹珠中匹配颜色的样本比例。在下面的代码中，我们运行了 10,000 轮从我们的容器中抽取两个弹珠的过程。利用这些样本，我们可以直接计算具有匹配弹珠的样本比例：

`n` `=` `10_000`
`samples` `=` `[``np``.``random``.``choice``(``urn``,` `size``=``2``,` `replace``=``False``)` `for` `_` `in` `range``(``n``)``]`
`is_matching` `=` `[``marble1` `==` `marble2` `for` `marble1``,` `marble2` `in` `samples``]`
`print``(``f``"``Proportion of samples with matching marbles:` `{``np``.``mean``(``is_matching``)``}``"``)`

Proportion of samples with matching marbles: 0.4032

我们刚刚进行了一项模拟研究。我们对np.random.choice的调用模拟了从容器中无放回地抽取两个弹珠的概率过程。每次对np.random.choice的调用给我们一个可能的样本。在模拟研究中，我们重复这个概率过程多次（在这种情况下是10_000次），以获得大量样本。然后，我们利用这些样本的典型行为来推理出我们可能从概率过程中获得的结果。虽然这可能看起来像是一个假设性的例子（确实如此），但请考虑如果我们将弹珠替换为在线约会服务中的人，用更复杂的属性替换颜色，并可能使用神经网络来评分匹配，你就能开始看到更复杂分析的基础了。

到目前为止，我们关注的是样本，但我们通常对我们可能观察到的样本与它可以告诉我们关于最初在容器中的“种群”弹珠之间的关系感兴趣。

我们可以将数据范围的类比与第二章联系起来：从坛子中抽取的一组大理石是一个样本，而放置在坛子中的所有大理石是接触框架，在这种情况下，我们认为与种群相同。这种模糊了访问框架和种群之间的差异指向了模拟与现实之间的差距。模拟往往简化模型。尽管如此，它们可以为现实世界的现象提供有用的见解。

在我们不在抽样过程中替换大理石的坛子模型中，有一个常见的选择方法称为简单随机样本。我们接下来描述这种方法及其他基于它的抽样技术。

抽样设计

从坛子中不替换地抽取大理石的过程等同于一个简单随机样本。在简单随机样本中，每个样本被选中的机会相同。虽然方法名称中含有“简单”一词，但构建一个简单随机样本通常并不简单，在许多情况下也是最佳的抽样过程。此外，如果我们诚实地说，它有时也会有些令人困惑。

为了更好地理解这种抽样方法，我们回到坛子模型。考虑一个有七个大理石的坛子。我们不给大理石着色，而是用字母A到G对每个大理石进行标记。由于每个大理石都有不同的标签，我们可以更清楚地识别我们可能得到的所有可能样本。让我们从坛子中不替换地选择三个大理石，并使用itertools库生成所有组合的列表：

`from` `itertools` `import` `combinations`

`all_samples` `=` `[``"``"``.``join``(``sample``)` `for` `sample` `in` `combinations``(``"``ABCDEFG``"``,` `3``)``]`
`print``(``all_samples``)`
`print``(``"``Number of Samples:``"``,` `len``(``all_samples``)``)`

['ABC', 'ABD', 'ABE', 'ABF', 'ABG', 'ACD', 'ACE', 'ACF', 'ACG', 'ADE', 'ADF', 'ADG', 'AEF', 'AEG', 'AFG', 'BCD', 'BCE', 'BCF', 'BCG', 'BDE', 'BDF', 'BDG', 'BEF', 'BEG', 'BFG', 'CDE', 'CDF', 'CDG', 'CEF', 'CEG', 'CFG', 'DEF', 'DEG', 'DFG', 'EFG']
Number of Samples: 35

我们的列表显示，有 35 个唯一的三大理石集。我们可以以六种不同的方式抽取每组集合。例如，集合{ A , B , C }可以被抽样：

`from` `itertools` `import` `permutations`

`print``(``[``"``"``.``join``(``sample``)` `for` `sample` `in` `permutations``(``"``ABC``"``)``]``)`

['ABC', 'ACB', 'BAC', 'BCA', 'CAB', 'CBA']

在这个小例子中，我们可以全面了解我们可以从坛子中抽取任意三个大理石的所有方法。

由于从七个种群中选择的每组三个大理石同等可能地发生，任何一个特定样本的机会必须是1 / 35：

P ( A B C ) = P ( ABD ) = ⋯ = P ( EFG ) = 1 35

我们使用特殊符号P表示“概率”或“机会”，并且我们将语句P ( A B C )读作“样本包含标记为 A、B 和 C 的大理石的机会，无论顺序如何”。

我们可以使用从坛子中所有可能样本的枚举来回答关于这一随机过程的其他问题。例如，要找出大理石A在样本中的机会，我们可以加总所有包含A的样本的机会。共有 15 个，因此机会是：

P ( A i s i n t h e s a m p l e ) = 15 35 = 3 7

当难以列举和计算所有可能的样本时，我们可以使用模拟来帮助理解这一概率过程。

注意

许多人错误地认为简单随机样本的定义特性是每个单位有相同的抽样机会。然而，情况并非如此。从人口中抽取n个单位的简单随机样本意味着每个N个单位的所有可能的n个集合被选择的机会相同。稍有变体的是带放回的简单随机样本，在这种情况下，单位/彩球在每次抽取后都被放回罐子。这种方法也具有每个N个单位人口的n个样本被选中的属性。不同之处在于，因为同一彩球可以在样本中出现多次，所以可能的n个单位集合更多。

简单随机样本（及其对应的罐子）是更复杂的调查设计的主要构建块。我们简要描述了两种更广泛使用的设计：

分层抽样

将人口分成互不重叠的群体，称为层（一个群体称为层，多个称为层），然后从每个群体中简单随机抽取样本。这就像每个层都有一个单独的罐子，并且独立地从每个罐子中抽取彩球。这些层的大小可以不同，并且我们不需要从每个层中取相同数量的彩球。

簇抽样

将人口分成互不重叠的子群体，称为簇，从簇中简单随机抽取样本，并将簇中的所有单位包括在样本中。我们可以将这看作是从一个包含大彩球的罐子中简单随机抽样，而这些大彩球本身是装有小彩球的容器。（大彩球的数量不一定相同。）当打开时，大彩球样本变成小彩球样本。（簇通常比层小。）

例如，我们可以将标记为A-G的七个彩球，组织成三个簇：(A,B)，(C,D)和(E,F,G)。然后，大小为一的簇样本有同等机会从这三个簇中抽取任何一个。在这种情况下，每个彩球被抽取为样本的机会相同：

P ( A in sample ) = P ( cluster ( A , B ) chosen ) = 1 3 P ( B in sample ) = P ( cluster ( A , B ) chosen ) = 1 3 ⋮ P ( G in sample ) = P ( cluster ( E , F , G ) chosen ) = 1 3

但并非每种元素组合的发生概率都相等：样本不可能同时包含A和C，因为它们位于不同的簇中。

经常，我们对样本的总结感兴趣；换句话说，我们对统计量感兴趣。对于任何样本，我们可以计算统计量，而罐模型帮助我们找到该统计量可能具有的值的分布。接下来，我们检查我们简单示例的统计量分布。

统计量的抽样分布

假设我们有兴趣测试新燃料箱设计的失败压力，这样做成本高昂，因为我们需要摧毁燃料箱，并且可能需要测试多个燃料箱以解决制造上的变化。

我们可以使用罐模型选择要测试的原型，并且可以通过失败测试的原型比例总结我们的测试结果。罐模型为我们提供了每个样本被选择的相同机会，因此压力测试结果代表了整体群体。

为了简单起见，假设我们有七个与之前弹珠类似标记的燃料箱。让我们看看如果选择时，如何处理坦克A、B、D和F未通过压力测试，而选择通过测试的C、E和G。

对于每三颗弹珠的样本，我们可以根据这些四个次品原型中有多少个来计算失败的比例。我们举几个计算示例：

样本	ABC	BCE	BDF	CEG
失败比例	2/3	1/3	1	0

由于我们从罐中抽取三颗弹珠，唯一可能的样本比例是0、1 / 3、2 / 3和1，对于每个三元组，我们可以计算其相应的比例。有四个样本使我们得到全部测试失败（样本比例为 1）。它们分别是A B D、A B F、A D F和B D F，因此观察到样本比例为 1 的机会为4 / 35。我们可以将样本比例的分布总结为一张表格，称为样本比例的抽样分布：

失败比例	样本数	样本比例
0	1	1/35 ≈ 0.03
1/3	12	12/35 ≈ 0.34
2/3	18	18/35 ≈ 0.51
1	4	4/35 ≈ 0.11
总数	35	1

尽管这些计算相对直观，我们可以通过仿真研究来近似它们。为此，我们反复从我们的总体中取三个样本，比如说一万次。对于每个样本，我们计算失败比例。这给我们 10,000 个模拟样本比例。仿真比例表应接近抽样分布。我们通过仿真研究确认这一点。

模拟抽样分布

仿真可以是理解复杂随机过程的强大工具。在我们七个燃料箱的例子中，我们能够考虑来自相应盒模型的所有可能样本。然而，在具有大量人口和样本以及更复杂抽样过程的情况下，直接计算某些结果的概率可能并不可行。在这些情况下，我们经常转向仿真，以提供我们无法直接计算的数量的准确估计。

让我们设定一个问题，即寻找三个燃料箱的简单随机样本中失败比例的抽样分布，作为一个盒子模型。由于我们关心罐是否故障，我们使用 1 表示失败和 0 表示通过，这给我们一个标记如下的罐子：

`urn` `=` `[``1``,` `1``,` `0``,` `1``,` `0``,` `1``,` `0``]`

我们已经使用 1 表示失败和 0 表示通过对罐头 A 至 G 进行编码，因此我们可以取样本的平均值来获得样本中的失败比例：

`sample` `=` `np``.``random``.``choice``(``urn``,` `size``=``3``,` `replace``=``False``)`
`print``(``f``"``Sample:` `{``sample``}``"``)`
`print``(``f``"``Prop Failures:` `{``sample``.``mean``(``)``}``"``)`

Sample: [1 0 0]
Prop Failures: 0.3333333333333333

在仿真研究中，我们重复抽样过程数千次以获得数千个比例，然后从我们的仿真中估计比例的抽样分布。在这里，我们构建了 10,000 个样本（因此有 10,000 个比例）：

`samples` `=` `[``np``.``random``.``choice``(``urn``,` `size``=``3``,` `replace``=``False``)` `for` `_` `in` `range``(``10_000``)``]`
`prop_failures` `=` `[``s``.``mean``(``)` `for` `s` `in` `samples``]`

我们可以研究这 10,000 个样本比例，并将我们的发现与我们已经使用所有 35 个可能样本的完全枚举计算的结果进行匹配。我们预计仿真结果与我们之前的计算非常接近，因为我们重复了许多次抽样过程。也就是说，我们想比较 10,000 个样本比例的分数是否为 0、1 / 3 、2 / 3 和 1，与我们确切计算的那些分数相匹配；这些分数是 1 / 35 、12 / 35 、18 / 35 和 4 / 35 ，约为 0.03 、0.34 、0.51 和 0.11 ：

`unique_els``,` `counts_els` `=` `np``.``unique``(``prop_failures``,` `return_counts``=``True``)`
`pd``.``DataFrame``(``{`
    `"``Proportion of failures``"``:` `unique_els``,`
    `"``Fraction of samples``"``:` `counts_els` `/` `10_000``,`
`}``)`

	失败比例	样本分数比
0	0.00	0.03
1	0.33	0.35
2	0.67	0.51
3	1.00	0.11

仿真结果非常接近我们之前计算的准确概率。

注意

模拟研究利用随机数生成器从随机过程中采样许多结果。从某种意义上讲，模拟研究将复杂的随机过程转化为我们可以使用本书中涵盖的广泛计算工具进行分析的数据。虽然模拟研究通常不提供特定假设的确凿证据，但它们可以提供重要的证据。在许多情况下，模拟是我们拥有的最准确的估计过程。

从一个瓮中抽取 0 和 1 的球是理解随机性的一个流行框架，这种机会过程已被正式命名为超几何分布，大多数软件都提供快速进行此过程模拟的功能。在下一节中，我们将模拟燃料箱例子中的超几何分布。

使用超几何分布进行模拟

而不是使用 random.choice，我们可以使用 numpy 的 random.hypergeometric 来模拟从瓮中取球并计算失败次数。random.hypergeometric 方法针对 0-1 瓮进行了优化，并允许我们在一次调用中请求 10,000 次模拟。为了完整起见，我们重复了我们的模拟研究并计算了经验比例：

`simulations_fast` `=` `np``.``random``.``hypergeometric``(`
    `ngood``=``4``,` `nbad``=``3``,` `nsample``=``3``,` `size``=``10_000`
`)`
`print``(``simulations_fast``)`

[1 1 2 ... 1 2 2]

（我们并不认为通过的是“坏”的；这只是一种命名惯例，将你想要计数的类型称为“好”的，其他的称为“坏”的。）

我们统计了 10,000 个样本中有 0、1、2 或 3 次失败的比例：

`unique_els``,` `counts_els` `=` `np``.``unique``(``simulations_fast``,` `return_counts``=``True``)`
`pd``.``DataFrame``(``{`
    `"``Number of failures``"``:` `unique_els``,`
    `"``Fraction of samples``"``:` `counts_els` `/` `10_000``,`
`}``)`

	失败次数	样本分数比
0	0	0.03
1	1	0.34
2	2	0.52
3	3	0.11

也许你已经在想：既然超几何分布如此受欢迎，为什么不提供可能值的确切分布呢？事实上，我们可以精确计算这些：

`from` `scipy``.``stats` `import` `hypergeom`

`num_failures` `=` `[``0``,` `1``,` `2``,` `3``]`
`pd``.``DataFrame``(``{`
    `"``Number of failures``"``:` `num_failures``,`
    `"``Fraction of samples``"``:` `hypergeom``.``pmf``(``num_failures``,` `7``,` `4``,` `3``)``,`
`}``)`

	失败次数	样本分数比
0	0	0.03
1	1	0.34
2	2	0.51
3	3	0.11

注意

在可能的情况下，最好使用第三方包中提供的功能来模拟命名分布，比如 numpy 中提供的随机数生成器，而不是编写自己的函数。最好利用他人开发的高效和准确的代码。话虽如此，偶尔从头开始构建可以帮助你理解算法，所以我们建议试一试。

或许最常见的两种机会过程是那些由从 0-1 瓮中抽取 1 的数量产生的过程：不替换抽取是超几何分布，替换抽取是二项式分布。

虽然这个模拟过程非常简单，我们本可以直接使用hypergeom.pmf来计算我们的分布，但我们希望展示模拟研究能够揭示的直觉。我们在本书中采用的方法是基于模拟研究开发对机会过程的理解。然而，在第十七章中，我们确实正式地定义了统计数据的概率分布的概念。

现在我们把模拟作为了解准确性的工具之一，可以重新审视第二章中的选举例子，并进行一次选举后研究，看看选民民意调查可能出了什么问题。这个模拟研究模仿了从 600 万个选民中抽取超过一千个弹珠（参与民意调查的选民）。我们可以检查偏见的潜在来源和民意调查结果的变化，并进行“如果”分析，看看如果从选民中抽取更多的弹珠会对预测产生怎样的影响。

示例：模拟选举民意调查的偏见和方差

2016 年，几乎所有关于美国总统选举结果的预测都是错误的。这是一个历史性的预测误差水平，震惊了统计学和数据科学界。在这里，我们探讨为什么几乎每一个政治民意调查都是如此自信，却又如此错误。这个故事既展示了模拟的力量，也揭示了数据的傲慢和偏见挑战的难度。

美国总统是由选举人团选出的，而不是由普通选民的投票决定。根据各州的人口大小，每个州被分配一定数量的选举人票数。通常情况下，谁在某州赢得了普选，谁就会获得该州所有的选举人票数。在选举前进行的民意调查帮助下，评论家确定了“争夺”州，在这些州中选举预计会非常接近，选举人票数可能会左右选举结果。

2016 年，民意调查机构正确预测了 50 个州中的 46 个的选举结果。这并不差！毕竟，对于那 46 个州来说，唐纳德·特朗普获得了 231 张选举人票，希拉里·克林顿获得了 232 张选举人票——几乎是平局，克林顿领先微弱。不幸的是，剩下的四个州，佛罗里达、密歇根、宾夕法尼亚和威斯康辛，被认定为争夺州，并且合计 75 张选举人票。这四个州的普选投票比例非常接近。例如，在宾夕法尼亚州，特朗普获得了 6,165,478 票中的 48.18%，而克林顿获得了 47.46%。在这些州，由于民意调查使用的样本量较小，很难预测选举结果。但是，在调查过程本身也存在更大的挑战。

许多专家研究了 2016 年选举结果，以分析并确定出了什么问题。根据美国公共舆论研究协会的说法，一项在线自愿参与的民意调查对受访者的教育程度进行了调整，但只使用了三个广泛的类别（高中或以下、部分大学和大学毕业）。民意调查人员发现，如果他们将具有高级学位的受访者与具有大学学位的受访者分开，那么他们将会将克林顿的估计百分比降低 0.5 个百分点。换句话说，在事后，他们能够确定受过教育程度较高的选民更愿意参与投票。这种偏见很重要，因为这些选民也倾向于喜欢克林顿而不是特朗普。

现在我们知道人们实际投票的方式，我们可以进行类似Manfred te Grotenhuis 等人的模拟研究，模拟不同情况下的选举民意调查，以帮助形成对准确性、偏见和方差的直觉。我们可以模拟并比较宾夕法尼亚州的两种情况下的民意调查：

• 受访者没有改变他们的想法，也没有隐藏他们投票给谁，并且代表了在选举日投票的人群。

• 受过高等教育的人更有可能回答，这导致了对克林顿的偏见。

我们的最终目标是了解在样本收集过程中完全没有偏见和存在少量非响应偏见的情况下，民意调查错误地将选举归因于希拉里·克林顿的频率。我们首先为第一种情况建立瓮模型。

宾夕法尼亚州瓮模型

我们对宾夕法尼亚州选民进行民意调查的瓮模型是一种事后情况，我们使用选举结果。这个瓮中有 6,165,478 个弹珠，每个选民一个。就像我们的小样本一样，我们在每个弹珠上写上他们投票给的候选人，从瓮中抽取 1,500 个弹珠（1,500 个是这些调查的典型大小），并统计特朗普、克林顿和任何其他候选人的选票。通过统计，我们可以计算特朗普相对于克林顿的领先优势。

由于我们只关心特朗普相对于克林顿的领先优势，我们可以将其他候选人的所有选票合并在一起。这样，每个弹珠都有三种可能的选票：特朗普、克林顿或其他。我们不能忽略“其他”类别，因为它会影响领先优势的大小。让我们将选民数分配给这三个群体：

`proportions` `=` `np``.``array``(``[``0.4818``,` `0.4746``,` `1` `-` `(``0.4818` `+` `0.4746``)``]``)`               
`n` `=` `1_500` 
`N` `=` `6_165_478`
`votes` `=` `np``.``trunc``(``N` `*` `proportions``)``.``astype``(``int``)`
`votes`

array([2970527, 2926135,  268814])

这个版本的瓮模型有三种类型的弹珠。它比超几何分布复杂一些，但仍然足够常见以具有命名分布：多元超几何分布。在 Python 中，具有两种以上弹珠类型的瓮模型是通过scipy.stats.multivariate_hypergeom.rvs方法实现的。该函数返回从瓮中抽取的每种类型的弹珠数量。我们调用函数如下：

`from` `scipy``.``stats` `import` `multivariate_hypergeom`

`multivariate_hypergeom``.``rvs``(``votes``,` `n``)`

array([727, 703,  70])

每次调用multivariate_hypergeom.rvs时，我们都会得到一个不同的样本和计数：

`multivariate_hypergeom``.``rvs``(``votes``,` `n``)`

array([711, 721,  68])

我们需要计算每个样本的特朗普领先： ( n T − n C ) / n ，其中 n T 是特朗普的选票数，n C 是克林顿的选票数。如果领先是正数，则样本显示特朗普胜出。

我们知道实际领先是 0.4818 – 0.4746 = 0.0072. 为了了解民意调查的变化情况，我们可以模拟反复从瓮中取出的机会过程，并检查我们得到的值。现在，我们可以模拟在宾夕法尼亚投票中的 1,500 名选民的 100,000 次调查：

`def` `trump_advantage``(``votes``,` `n``)``:`
    `sample_votes` `=` `multivariate_hypergeom``.``rvs``(``votes``,` `n``)`
    `return` `(``sample_votes``[``0``]` `-` `sample_votes``[``1``]``)` `/` `n`

`simulations` `=` `[``trump_advantage``(``votes``,` `n``)` `for` `_` `in` `range``(``100_000``)``]` 

平均而言，民调结果显示特朗普领先接近 0.7%，这与投票结果的构成相符：

`np``.``mean``(``simulations``)`

0.007177066666666666

然而，很多时候样本的领先是负数，这意味着在该选民样本中克林顿是赢家。下图显示了在宾夕法尼亚 1,500 名选民样本中特朗普优势的抽样分布。在 0 处的垂直虚线显示，特朗普更常被提及，但在 1,500 人的调查中，有很多次克林顿处于领先地位：

在 100,000 次模拟的调查中，我们发现特朗普大约 60%的时间是胜利者：

`np``.``mean``(``np``.``array``(``simulations``)` `>` `0``)`

0.60613

换句话说，一个样本即使没有任何偏见地收集，也将大约 60%的时间正确预测特朗普的胜利。而这种无偏样本将在 40%的时间错误。

我们使用了瓮模型来研究简单民意调查的变化，并找出了如果选择过程没有偏差时民意调查的预测可能是什么样子（弹珠是无法区分的，六百多万个弹珠中的每一种可能的 1,500 个弹珠集合是同等可能的）。接下来，我们看看当一点偏差进入混合时会发生什么。

具有偏差的瓮模型

根据 Grotenhuis 的说法，“在完美的世界里，民意调查从选民群体中取样，他们会明确表达自己的政治偏好，然后相应地投票。”^(2) 这就是我们刚刚进行的模拟研究。然而，在现实中，往往很难控制每一种偏见来源。

我们在这里研究了教育偏见对民调结果的影响。具体来说，我们检查了对克林顿有利的 0.5%偏见的影响。这种偏见实质上意味着我们在民意调查中看到了选民偏好的扭曲图片。克林顿的选票不是 47.46%，而是 47.96%，而特朗普是 48.18 – 0.5 = 47.68%。我们调整了瓮中弹珠的比例以反映这一变化：

`bias` `=` `0.005`
`proportions_bias` `=` `np``.``array``(``[``0.4818` `-` `bias``,` `0.4747` `+` `bias``,` 
                             `1` `-` `(``0.4818` `+` `0.4746``)``]``)`
`proportions_bias`

array([0.48, 0.48, 0.04])

`votes_bias` `=` `np``.``trunc``(``N` `*` `proportions_bias``)``.``astype``(``int``)`
`votes_bias`

array([2939699, 2957579,  268814])

当我们再次进行模拟研究时，这次使用有偏的瓮，我们发现结果大不相同：

`simulations_bias` `=` `[``trump_advantage``(``votes_bias``,` `n``)` `for` `_` `in` `range``(``100_000``)``]` 

`np``.``mean``(``np``.``array``(``simulations_bias``)` `>` `0``)` 

0.44967

现在，特朗普在大约 45%的民意调查中领先。请注意，两次模拟的直方图形状相似。它们对称，并且尾部长度合理。也就是说，它们似乎大致遵循正态曲线。第二个直方图稍微向左移动，反映了我们引入的非响应偏倚。增加样本量会有所帮助吗？我们接下来研究这个话题。

进行更大规模的民调

通过我们的模拟研究，我们可以了解更大样本对样本领先的影响。例如，我们可以尝试使用 12,000 个样本，是实际民调规模的 8 倍，并针对无偏和有偏情况运行 100,000 次模拟：

`simulations_big` `=` `[``trump_advantage``(``votes``,` `12_000``)` `for` `_` `in` `range``(``100_000``)``]`
`simulations_bias_big` `=` `[``trump_advantage``(``votes_bias``,` `12_000``)`
                        `for` `_` `in` `range``(``100_000``)``]`

`scenario_no_bias` `=` `np``.``mean``(``np``.``array``(``simulations_big``)` `>` `0``)`
`scenario_bias` `=` `np``.``mean``(``np``.``array``(``simulations_bias_big``)` `>` `0``)`
`print``(``scenario_no_bias``,` `scenario_bias``)`

0.78968 0.36935

模拟显示，在有偏情况下，只有大约三分之一的模拟中检测到特朗普的领先地位。这些结果的直方图分布比仅有 1,500 位选民的情况更窄。不幸的是，它已经偏离了正确的数值。我们并没有克服偏倚；我们只是对偏倚情况有了更准确的了解。大数据并没有拯救我们。此外，更大规模的民调还有其他问题。它们通常更难进行，因为民调员在有限的资源下工作，本来可以用来改善数据范围的努力被重新定向到扩展民调上：

事后，通过多个相同选举的民调，我们可以检测到偏倚。在对 600 个州级、州长、参议院和总统选举的 4,000 多次民意调查进行的选举后分析中，研究人员发现，平均来说，选举民调显示出大约 1.5 个百分点的偏倚，这有助于解释为什么许多民调预测都错了。

当胜利的边际相对较小时，就像 2016 年那样，更大的样本量可以减少抽样误差，但不幸的是，如果存在偏倚，那么预测结果接近偏倚估计。如果偏倚使得预测从一个候选人（特朗普）转向另一个（克林顿），那么我们就会看到一场“意外”的颠覆。民调员开发选民选择方案，试图减少偏倚，比如按教育水平分离选民的偏好。但是，就像这个案例一样，很难，甚至不可能考虑到新的、意外的偏倚来源。民调仍然有用，但我们需要承认偏倚问题，并在减少偏倚方面做得更好。

在这个例子中，我们使用了乌尔恩模型来研究民调中的简单随机样本。乌尔恩模型在随机对照实验中也是常见的应用之一。

例如：模拟疫苗随机试验

在药物试验中，试验志愿者接受的是新治疗方法或者安慰剂（一种假的治疗方法），研究人员控制志愿者分配到治疗组和安慰剂组的过程。在随机对照实验中，他们使用随机过程来进行这种分配。科学家们基本上使用乌尔恩模型来选择接受治疗和安慰剂（即接受安慰剂的那些人）的对象。我们可以模拟乌尔恩的随机机制，以更好地理解实验结果的变化和临床试验中效力的含义。

2021 年 3 月，底特律市长迈克·达根（Mike Duggan）拒绝接受超过 6000 剂强生公司的疫苗 国家新闻，并表示他的市民应该“得到最好的”。市长指的是疫苗的有效率，据报道约为 66%。相比之下，Moderna 和 Pfizer 的疫苗报告的有效率约为 95%。

乍看之下，达根市长的理由似乎合理，但是三个临床试验的范围并不可比，这意味着直接比较试验结果是有问题的。此外，CDC 认为 66%的有效率相当不错，这也是为什么它获得了紧急批准。

让我们依次考虑范围和效力的要点。

范围

记住，当我们评估数据的范围时，我们考虑研究的对象、时间和地点。对于强生公司的临床试验，参与者：

• 包括 18 岁及以上成年人，其中大约 40%具有与患严重 COVID-19 风险增加相关的既往病史。

• 从 2020 年 10 月到 11 月，参与者被招募入研究。

• 来自八个国家，涵盖三大洲，包括美国和南非。

Moderna 和 Pfizer 试验的参与者主要来自美国，大约 40%有既往病史，试验在 2020 年夏季早些时候进行。试验的时间和地点使得它们难以比较。COVID-19 病例在美国夏季时期处于低点，但在晚秋期间迅速上升。同时，当时南非的一种更具传染性的病毒变种正在迅速传播，这也是 J&J 试验的时候。

每个临床试验旨在通过将受试者随机分配到治疗组和对照组来测试疫苗在没有疫苗的情况下的效果。尽管每个试验的范围有所不同，但试验内的随机化使得治疗组和对照组的范围基本相似。这使得在同一试验中组之间可以进行有意义的比较。三个疫苗试验的范围差异足够大，使得直接比较这三个试验的结果变得复杂。

在进行 J&J 疫苗 的试验中，有 43,738 人参与了。这些参与者被随机分成两组。一半人接受新疫苗，另一半接受安慰剂，比如生理盐水。然后每个人都被随访 28 天，看是否感染了 COVID-19。

关于每位患者都记录了大量信息，比如他们的年龄、种族和性别，以及他们是否感染了 COVID-19，包括疾病的严重程度。28 天后，研究人员发现了 468 例 COVID-19 病例，其中治疗组有 117 例，对照组有 351 例。

将患者随机分配到治疗组和对照组为科学家提供了评估疫苗有效性的框架。典型的推理如下：

1. 从疫苗无效的假设开始。

2. 因此，468 名感染 COVID-19 的人无论是否接种疫苗都会感染。

3. 试验中剩下的 43,270 人，无论是否接种疫苗，都不会生病。

4. 在治疗组有 117 人，对照组有 351 人的拆分，完全是由于将参与者分配到治疗或对照组的偶然过程。

我们可以建立一个反映这种情况的瓮模型，然后通过模拟研究实验结果的行为。

用于随机分配的瓮模型

我们的瓮有 43,738 个弹珠，代表临床试验中的每个人。由于其中有 468 例 COVID-19 病例，我们用 1 标记了 468 个弹珠，剩下的 43,270 个用 0 标记。我们从瓮中抽出一半的弹珠（21,869 个）接受治疗，剩下的一半接受安慰剂。实验的关键结果仅仅是从瓮中随机抽取的标记为 1 的弹珠数量。

我们可以模拟这个过程，以了解在这些假设下从瓮中最多抽取 117 个标记为 1 的弹珠的可能性有多大。由于我们从瓮中抽出一半的弹珠，我们预计大约有一半的 468 个弹珠，即 234 个，会被抽出。模拟研究给出了随机分配过程可能产生的变异的概念。也就是说，模拟可以给出试验中治疗组病毒病例如此之少的近似机会。

注意

这个瓮模型涉及到几个关键假设，比如疫苗无效的假设。跟踪这些假设的依赖很重要，因为我们的模拟研究给出了仅在这些关键假设下观察到的罕见结果的近似。

与以往一样，我们可以使用超几何概率分布模拟瓮模型，而不必从头编写偶然过程的程序：

`simulations_fast` `=` `np``.``random``.``hypergeometric``(``ngood``=``468``,` `nbad``=``43270``,` 
                                            `nsample``=``21869``,` `size``=``500000``)`

在我们的模拟中，我们重复了 500,000 次随机分配到治疗组的过程。事实上，我们发现 500,000 次模拟中没有一次发生 117 例或更少的情况。如果疫苗真的没有效果，看到如此少的 COVID-19 病例将是一个极其罕见的事件。

在解释了比较具有不同范围和预防 COVID-19 严重病例有效性的药物试验之后，达根市长撤回了他最初的声明，说：“我完全相信强生公司的疫苗既安全又有效。”^(3)

本例表明

• 使用随机过程将受试者分配到临床试验治疗组中，可以帮助我们回答各种假设情景。

• 考虑数据范围可以帮助我们确定是否合理比较不同数据集的数据。

从瓮中抽取弹珠的模拟是研究调查样本和控制实验可能结果的有用抽象。这种模拟有效，因为它模拟了用于选择样本或将人员分配到治疗中的机会机制。在我们测量自然现象的设置中，我们的测量倾向于遵循类似的机会过程。正如第二章中所述，仪器通常具有与其相关的误差，我们可以使用瓮来表示测量对象的变异性。

例子：测量空气质量

在美国各地，用于测量空气污染的传感器被广泛使用，包括个人、社区组织以及州和地方空气监测机构。例如，2020 年 9 月的两天，约 60 万加利福尼亚人和 50 万俄勒冈人在紫外空气地图上查看了火灾蔓延和疏散计划。(紫外空气从传感器收集的众包数据制作空气质量地图。)

传感器测量空气中直径小于 2.5 微米的颗粒物质的数量（测量单位为每立方米的微克数：μg/m³）。记录的测量值是两分钟内的平均浓度。例如，尽管颗粒物质的水平在一天中会发生变化，比如人们通勤上下班，但是在一天的某些时间，比如午夜，我们预计两分钟的平均值在半小时内变化不大。如果我们检查这些时间段内的测量值，我们可以了解仪器记录的综合变异性和空气中颗粒物的混合情况。

任何人都可以访问 PurpleAir 网站上的传感器测量值。该网站提供下载工具，并且对 PurpleAir 地图上出现的任何传感器都可用数据。我们下载了一个传感器在 24 小时内的数据，并选择了一天中分布在整个时间段内读数大致稳定的三个半小时时间间隔。这为我们提供了三组 15 个两分钟平均值，总共 45 个测量值：

	aq2.5	time	hour	meds	diff30
0	6.14	2022-04-01 00:01:10 UTC	0	5.38	0.59
1	5.00	2022-04-01 00:03:10 UTC	0	5.38	-0.55
2	5.29	2022-04-01 00:05:10 UTC	0	5.38	-0.26
...	...	...	...	...	...
42	7.55	2022-04-01 19:27:20 UTC	19	8.55	-1.29
43	9.47	2022-04-01 19:29:20 UTC	19	8.55	0.63
44	8.55	2022-04-01 19:31:20 UTC	19	8.55	-0.29

45 rows × 5 columns

线图可以让我们感受到测量值的变化。在一个 30 分钟的时间段内，我们期望测量值大致相同，除了空气中颗粒物的轻微波动和仪器的测量误差：

图表显示了一天中空气质量如何恶化，但在每个半小时间隔中，午夜、上午 11 点和下午 7 点时，空气质量分别大约为 5.4、6.6 和 8.6 μg/m³。我们可以将数据范围想象为：在特定位置的特定半小时时间间隔内，环绕传感器的空气中有一个平均颗粒浓度。这个浓度是我们的目标，我们的仪器——传感器，进行了许多形成样本的测量。 （请参见第二章中有关这一过程的飞镖板类比。）如果仪器工作正常，测量值将集中在目标上：30 分钟平均值。

为了更好地了解半小时间隔内的变化，我们可以检查测量值与相应半小时的中位数之间的差异。这些“误差”的分布如下：

直方图显示，测量值的典型波动通常小于 0.5 μg/m³，很少大于 1 μg/m³。使用仪器时，我们经常考虑其相对标准误差，即标准偏差占平均值的百分比。这 45 个偏差的标准偏差为：

`np``.``std``(``pm``[``'``diff30``'``]``)`

0.6870817156282193

鉴于每小时的测量值范围在 5 至 9 μg/m³之间，相对误差为 8%至 12%，这是相当精确的。

我们可以使用乌尔恩模型来模拟这个测量过程的变异性。我们将所有 45 次测量与它们的 30 分钟中位数的偏差放入乌尔恩，并通过从中抽取 15 次（有放回）并将抽取的偏差加到一个假设的 30 分钟平均值来模拟一个 30 分钟的空气质量测量序列：

`urn` `=` `pm``[``"``diff30``"``]`

`np``.``random``.``seed``(``221212``)`
`sample_err` `=` `np``.``random``.``choice``(``urn``,` `size``=``15``,` `replace``=``True``)`
`aq_imitate` `=` `11` `+` `sample_err`

我们可以为这组人工测量数据添加一条线图，并与之前的三个真实数据进行比较：

从模拟数据的线图形状看，它与其他数据相似，这表明我们对测量过程的模型是合理的。不幸的是，我们不知道这些测量是否接近真实的空气质量。为了检测仪器的偏差，我们需要与更精确的仪器进行比较，或者在空气含有已知颗粒物数量的受保护环境中进行测量。事实上，研究人员 发现低湿度会使读数偏高。在第十二章中，我们对 PurpleAir 传感器数据进行了更全面的分析，并校准仪器以提高其准确性。

总结

在本章中，我们使用了从乌尔恩抽取彩球的类比来模拟从人群中随机抽样和在实验中随机分配受试者到治疗组的过程。这个框架使我们能够进行针对假设调查、实验或其他随机过程的模拟研究，以研究它们的行为。我们发现了在假设治疗无效的情况下观察到特定临床试验结果的概率，并且基于实际选举投票结果的样本研究了对克林顿和特朗普的支持情况。这些模拟研究使我们能够量化随机过程中的典型偏差，并近似总结统计的分布，例如特朗普领先克林顿的情况。这些模拟研究揭示了统计量的抽样分布，并帮助我们回答关于在乌尔恩模型下观察到类似结果的可能性问题。

球罐模型简化为几个基本要素：罐中弹珠的数量，每个弹珠上的内容，从罐中抽取的弹珠数量，以及抽取过程中是否替换。从这些基础出发，我们可以模拟越来越复杂的数据设计。然而，罐模型的实用性关键在于将数据设计映射到罐子中。如果样本不是随机抽取的，受试者没有随机分配到治疗组，或者测量不是在校准良好的设备上进行的，那么这个框架在帮助我们理解数据并做出决策方面就会显得力不从心。另一方面，我们也需要记住，罐模型是对实际数据收集过程的简化。如果现实中数据收集存在偏差，那么我们在模拟中观察到的随机性并不能完整地捕捉到全貌。太多时候，数据科学家们忽略这些烦恼，只关注罐模型描述的变异性。这是预测 2016 年美国总统选举结果的调查中的主要问题之一。

在这些例子中，我们学习的总结统计数据是作为例子的一部分给出的。在下一章中，我们将讨论如何选择一个总结统计数据来代表这些数据。

^(1) Manfred te Grotenhuis 等人，《更好的民意抽样将更加怀疑希拉里·克林顿在 2016 年选举中可能获胜的潜力》，伦敦政治经济学院，2018 年 2 月 1 日。

^(2) Grotenhuis 等人，《更好的民意抽样将更加怀疑希拉里·克林顿在 2016 年选举中可能获胜的潜力》。

^(3) 不幸的是，尽管疫苗的有效性，美国食品药品监督管理局由于增加罕见且潜在致命的血栓风险，于 2022 年 5 月限制了 J&J 疫苗的使用。

第四章：模型与总结统计

我们在第二章中看到了数据范围的重要性，在第三章中看到了数据生成机制的重要性，例如可以用一个瓮模型来表示的机制。瓮模型解决了建模的一个方面：它描述了偶然变化，并确保数据代表了目标。良好的范围和代表性数据为从数据中提取有用信息奠定了基础，这在建模的另一部分中经常被称为数据中的信号。我们使用模型来近似这个信号，其中最简单的模型之一是常数模型，其中信号由一个单一数字（如均值或中位数）来近似。其他更复杂的模型总结了数据中特征之间的关系，例如空气质量中的湿度和颗粒物质(第十二章)，社区中的上升流动性和通勤时间(第十五章)，以及动物的身高和体重(第十八章)。这些更复杂的模型也是从数据中构建的近似值。当模型很好地适合数据时，它可以提供对世界的有用近似描述或仅仅是数据的有用描述。

在本章中，我们通过一个损失的形式介绍了模型拟合的基础知识。我们演示了如何通过考虑由简单总结描述数据引起的损失来建模数据中的模式，即常数模型。我们在第十六章深入探讨了瓮模型与拟合模型之间的联系，其中我们检查了拟合模型时信号和噪声之间的平衡，并在第十七章中讨论了推断、预测和假设检验的主题。

常数模型让我们可以从损失最小化的角度在简单情境中介绍模型拟合，它帮助我们将总结统计（如均值和中位数）与后续章节中的更复杂建模场景联系起来。我们从一个例子开始，该例子使用关于公交车晚点的数据来介绍常数模型。

常数模型

一个乘客，Jake，经常在西雅图市中心第三大道和派克街交界处的北行 C 路公交车站乘坐公交车。^(1) 这辆公交车应该每 10 分钟到达一次，但是 Jake 注意到有时候他等车的时间很长。他想知道公交车通常晚到多久。Jake 成功获取了从华盛顿州交通中心获得的公交车的预定到达时间和实际到达时间。根据这些数据，他可以计算每辆公交车晚点到达他所在站点的分钟数：

`times` `=` `pd``.``read_csv``(``'``data/seattle_bus_times_NC.csv``'``)`
`times`

	路线	方向	预定时间	实际时间	晚到分钟数
0	C	北行	2016-03-26 06:30:28	2016-03-26 06:26:04	-4.40
1	C	往北	2016-03-26 01:05:25	2016-03-26 01:10:15	4.83
2	C	往北	2016-03-26 21:00:25	2016-03-26 21:05:00	4.58
...	...	...	...	...	...
1431	C	往北	2016-04-10 06:15:28	2016-04-10 06:11:37	-3.85
1432	C	往北	2016-04-10 17:00:28	2016-04-10 16:56:54	-3.57
1433	C	往北	2016-04-10 20:15:25	2016-04-10 20:18:21	2.93

1434 rows × 5 columns

数据表中的 minutes_late 列记录了每辆公交车的迟到时间。请注意，有些时间是负数，这意味着公交车提前到达。让我们来看一下每辆公交车迟到时间的直方图：

`fig` `=` `px``.``histogram``(``times``,` `x``=``'``minutes_late``'``,` `width``=``450``,` `height``=``250``)`
`fig``.``update_xaxes``(``range``=``[``-``12``,` `60``]``,` `title_text``=``'``Minutes late``'``)`
`fig`

我们已经可以在数据中看到一些有趣的模式。例如，许多公交车提前到达，但有些车晚到超过 20 分钟。我们还可以看到明显的众数（高点），在 0 附近，意味着许多公交车大致按时到达。

要了解这条路线上公交车通常晚到多久，我们希望通过一个常数来总结迟到情况 —— 这是一个统计量，一个单一的数字，比如均值、中位数或众数。让我们找到数据表中 minutes_late 列的每个这些摘要统计量。

从直方图中，我们估计数据的众数是 0，并使用 Python 计算均值和中位数：

mean:    1.92 mins late
median:  0.74 mins late
mode:    0.00 mins late

自然地，我们想知道这些数字中哪一个最能代表迟到的摘要情况。我们不想依赖经验法则，而是采取更正式的方法。我们为公交车迟到建立一个常数模型。让我们称这个常数为 θ （在建模中，θ 通常被称为参数）。例如，如果我们考虑 θ = 5 ，那么我们的模型大致认为公交车通常晚到五分钟。

现在，θ = 5 并不是一个特别好的猜测。从迟到时间的直方图中，我们看到有更多的点接近 0 而不是 5。但是目前还不清楚 θ = 0（众数）是否比 θ = 0.74（中位数）、θ = 1.92（均值）或者完全不同的其他值更好。为了在不同的 θ 值之间做出选择，我们希望给每个 θ 值分配一个评分，以衡量这个常数如何与数据匹配。换句话说，我们希望评估用常数近似数据所涉及的损失，比如 θ = 5 。而理想情况下，我们希望选择最能匹配我们数据的常数，也就是具有最小损失的常数。在下一节中，我们将更正式地描述损失，并展示如何使用它来拟合模型。

最小化损失

我们想要通过一个常数来模拟北向 C 路线的延迟，这个常数我们称之为θ，并且我们想要利用每辆公交车实际延迟的分钟数的数据来找出一个合适的θ值。为此，我们使用一个损失函数，这个函数衡量我们的常数θ与实际数据之间的差距。

损失函数是一个数学函数，接受θ和数据值y作为输入。它输出一个单一的数字，损失，用来衡量θ和y之间的距离。我们将损失函数写成l ( θ , y )。

按照惯例，损失函数对于较好的θ值输出较低的值，对于较差的θ值输出较大的值。为了使常数适应我们的数据，我们选择产生所有θ选择下平均损失最低的特定θ。换句话说，我们找到最小化数据平均损失的θ，其中y 1 , … , y n 。更正式地，我们将平均损失写为L ( θ , y 1 , y 2 , … , y n )，其中：

L ( θ , y 1 , y 2 , … , y n ) = mean { l ( θ , y 1 ) , l ( θ , y 2 ) , … , l ( θ , y n ) } = 1 n ∑ i = 1 n l ( θ , y i )

作为简写，我们经常使用向量y = [ y 1 , y 2 , … , y n ]。然后我们可以将平均损失写为：

L ( θ , y ) = 1 n ∑ i = 1 n l ( θ , y i )

注意

注意，l ( θ , y ) 告诉我们模型对单个数据点的损失，而L ( θ , y ) 给出模型对所有数据点的平均损失。大写的L帮助我们记住平均损失结合了多个较小的l值。

一旦我们定义了一个损失函数，我们可以找到产生最小平均损失的 θ 值。我们称这个最小化的值为 θ ^ 。换句话说，对于所有可能的 θ 值，θ ^ 是为我们的数据产生最小平均损失的那个值。我们称这个优化过程为模型拟合；它找到了适合我们数据的最佳常数模型。

接下来，我们看一下两个特定的损失函数：绝对误差和平方误差。我们的目标是拟合模型并找到每个损失函数的 θ ^ 。

平均绝对误差

我们从绝对误差损失函数开始。这里是绝对损失背后的理念。对于某个 θ 值和数据值 y ：

1. 找到误差，y − θ 。

2. 取误差的绝对值，| y − θ | 。

因此，损失函数为 l ( θ , y ) = | y − θ | 。

取误差的绝对值是将负误差转换为正误差的简单方法。例如，点 y = 4 距离 θ = 2 和 θ = 6 同样远，因此误差是同样“糟糕”的。

平均绝对误差被称为平均绝对误差（MAE）。MAE 是每个单独绝对误差的平均值：

L ( θ , y ) = 1 n ∑ i = 1 n | y i − θ |

注意，MAE 的名称告诉你如何计算它：取误差的绝对值的平均值，{ y i − θ } 。

我们可以编写一个简单的 Python 函数来计算这个损失：

`def` `mae_loss``(``theta``,` `y_vals``)``:`
    `return` `np``.``mean``(``np``.``abs``(``y_vals` `-` `theta``)``)`

让我们看看当我们只有五个数据点 [ – 1 , 0 , 2 , 5 , 10 ] 时，这个损失函数的表现。我们可以尝试不同的 θ 值，并查看每个值对应的 MAE 输出：

我们建议通过手工验证一些这些损失值，以确保您理解如何计算 MAE。

在我们尝试的θ值中，我们发现θ=2具有最低的平均绝对误差。对于这个简单的例子，2 是数据值的中位数。这不是巧合。现在让我们来检查公交晚点时间原始数据的平均损失是多少，当我们将θ设置为分钟数的众数、中位数和平均数时，分别得到的 MAE 为：

我们再次看到中位数（中间图）比众数和平均数（左图和右图）有更小的损失。事实上，对于绝对损失，最小化的θ^是中位数{y1,y2,…,yn}。

到目前为止，我们通过简单尝试几个值并选择最小损失的值来找到了θ的最佳值。为了更好地理解θ的 MAE 作为函数的情况，我们可以尝试更多的θ值，并绘制一条曲线，显示L(θ,y)随θ变化的情况。我们为前述的五个数据值[–1,0,2,5,10]绘制了这条曲线：

前面的图表显示，实际上θ=2是这五个值的最佳选择。请注意曲线的形状。它是分段线性的，线段在数据值（–1, 0, 2 和 5）的位置连接。这是绝对值函数的特性。对于大量数据，平坦部分不那么明显。我们的公交数据有超过 1400 个数据点，MAE 曲线看起来更加平滑：

我们可以利用这个图来确认数据的中位数是最小化值；换句话说，θ^=0.74。这个图不是真正的证明，但希望它足够令你信服。

接下来，让我们看看另一个平方误差的损失函数。

均方误差

我们已经将常数模型拟合到我们的数据中，并发现使用均方误差时，最小化器是中位数。现在我们将保持模型不变，但切换到不同的损失函数：平方误差。我们不再取每个数据值y和常数θ之间的绝对差值，而是将误差平方。也就是说，对于某个θ值和数据值y：

1. 找到误差，y − θ。

2. 计算误差的平方，( y − θ ) 2。

这给出了损失函数l ( θ , y ) = ( y − θ ) 2。

与以往一样，我们希望利用所有数据来找到最佳的θ，因此我们计算均方误差，简称为 MSE：

L ( θ , y ) = L ( θ , y 1 , y 2 , … , y n ) = 1 n ∑ i = 1 n ( y i − θ ) 2

我们可以编写一个简单的 Python 函数来计算 MSE：

`def` `mse_loss``(``theta``,` `y_vals``)``:`
    `return` `np``.``mean``(``(``y_vals` `-` `theta``)` `*``*` `2``)`

让我们再次尝试均值、中位数和众数作为 MSE 的潜在最小化器：

现在，当我们使用 MSE 损失来拟合常数模型时，我们发现均值（右图）的损失小于中位数和众数（左图和中图）。

让我们绘制给定数据的不同θ值的 MSE 曲线。曲线显示最小化值θ ^接近 2：

这条曲线的一个特点是，与 MAE 相比，MSE 增长得非常迅速（注意纵轴上的范围）。这种增长与平方误差的性质有关；它对远离θ的数据值施加了更高的损失。如果θ = 10且y = 110，则平方损失为( 10 − 110 ) 2 = 10 , 000，而绝对损失为| 10 − 110 | = 100。因此，MSE 对异常大的数据值更为敏感。

从均方误差曲线来看，最小化的 θ ^ 看起来是 y 的均值。同样，这不是巧合；数据的均值总是与平方误差的 θ ^ 相符。我们展示了这是如何从均方误差的二次特性推导出来的。在此过程中，我们展示了平方损失作为方差和偏差项之和的常见表示，这是模型拟合中的核心。首先，我们在损失函数中添加和减去 y ¯ ，并展开平方如下：

L ( θ , y ) = 1 n ∑ i = 1 n ( y i − θ ) 2 = 1 n ∑ i = 1 n [ ( y i − y ¯ ) + ( y ¯ − θ ) ] 2 = 1 n ∑ i = 1 n [ ( y i − y ¯ ) 2 + 2 ( y i − y ¯ ) ( y ¯ − θ ) + ( y ¯ − θ ) 2 ]

接下来，我们将均方误差分解为这三个项的和，并注意中间项为 0，这是由于平均数的简单性质：∑ ( y i − y ¯ ) = 0 ：

1 n ∑ i = 1 n ( y i − y ¯ ) 2 + 1 n ∑ i = 1 n 2 ( y i − y ¯ ) ( y ¯ − θ ) + 1 n ∑ i = 1 n ( y ¯ − θ ) 2 = 1 n ∑ i = 1 n ( y i − y ¯ ) 2 + 2 ( y ¯ − θ ) 1 n ∑ i = 1 n ( y i − y ¯ ) + 1 n ∑ i = 1 n ( y ¯ − θ ) 2 = 1 n ∑ i = 1 n ( y i − y ¯ ) 2 + ( y ¯ − θ ) 2

在剩余的两个术语中，第一个不涉及 θ 。你可能认识到它是数据的方差。第二个术语始终为非负。它称为偏差平方。第二个术语，偏差平方，在 θ = y ¯ 时为 0，因此 θ ^ 给出任何数据集的最小均方误差。

我们已经看到，对于绝对损失，最好的常数模型是中位数，但对于平方误差，是均值。选择损失函数是模型拟合的一个重要方面。

选择损失函数

现在我们已经处理了两个损失函数，我们可以回到最初的问题：我们如何选择使用中位数、均值或模式？由于这些统计量最小化不同的损失函数，^(2) 我们可以等价地问：对于我们的问题，什么是最合适的损失函数？为了回答这个问题，我们看一下问题的背景。

与平均绝对误差（MAE）相比，均方误差（MSE）在公交车迟到（或提前）很多时会导致特别大的损失。希望了解典型迟到时间的公交车乘客会使用 MAE 和中位数（晚 0.74 分钟），但是讨厌意外大迟到时间的乘客可能会用均方误差和均值（晚 1.92 分钟）来总结数据。

如果我们想进一步优化模型，我们可以使用更专业的损失函数。例如，假设公交车提前到达时会在站点等待直到预定出发时间；那么我们可能希望将早到视为 0 损失。如果一个非常迟到的公交车比一个中度迟到的公交车更加令人恼火，我们可能会选择一个非对称损失函数，对超级迟到的惩罚更大。

实质上，在选择损失函数时上下文很重要。通过仔细考虑我们计划如何使用模型，我们可以选择一个有助于我们做出良好数据驱动决策的损失函数。

总结

我们介绍了恒定模型：一个通过单一值汇总数据的模型。为了拟合恒定模型，我们选择了一个度量给定常数与数据值匹配程度的损失函数，并计算所有数据值的平均损失。我们发现，根据损失函数的选择，我们得到不同的最小化值：我们发现平均值最小化平均平方误差（MSE），中位数最小化平均绝对误差（MAE）。我们还讨论了如何结合问题的上下文和知识来选择损失函数。

将模型拟合到损失最小化的概念将简单的汇总统计数据（如平均数、中位数和众数）与更复杂的建模情况联系起来。我们在建模数据时采取的步骤适用于许多建模场景：

1. 选择模型的形式（例如恒定模型）。

2. 选择一个损失函数（如绝对误差）。

3. 通过最小化所有数据的损失来拟合模型（如平均损失）。

在本书的其余部分，我们的建模技术扩展到这些步骤的一个或多个。我们引入新模型、新损失函数和新的最小化损失技术。第五章 重新审视了公交车晚点到达站点的研究。这一次，我们将问题呈现为案例研究，并访问数据科学生命周期的所有阶段。通过经历这些阶段，我们做出了一些不同寻常的发现；当我们通过考虑数据范围并使用瓮来模拟乘客到达公交车站时，我们发现建模公交车迟到不同于建模乘客等待公交车的经验。

^(1) 我们（作者）最初从名为杰克·范德普拉斯的数据科学家的分析中了解到公交到达时间数据。我们以他的名义命名本节的主角。

^(2) 众数最小化了一个称为 0-1 损失的损失函数。尽管我们尚未涵盖这种特定损失，但该过程是相同的：选择损失函数，然后找到最小化损失的内容。

第五章：案例研究：为什么我的公交车总是迟到？

Jake VanderPlas 的博客Pythonic Perambulations，提供了一个现代数据科学家的生活示例。作为数据科学家，我们在工作中、日常生活中以及个人生活中都看到数据，我们对这些数据可能带来的见解往往很好奇。在这第一个案例研究中，我们借鉴了 Pythonic Perambulations 上的一篇文章"等待时间悖论，或者，为什么我的公交车总是迟到？"，来模拟在西雅图街角等公交车的情景。我们涉及数据生命周期的每个阶段，但在这第一个案例研究中，我们的重点是思考问题、数据和模型的过程，而不是数据结构和建模技术。通过持续的模型和模拟研究，我们可以更好地理解这些问题。

VanderPlas 的文章受到他等公交车的经历的启发。等待时间总是比预期的长。这种经历与以下推理不符：如果每 10 分钟就来一辆公交车，并且你在随机时间到达车站，那么平均等待时间应该约为 5 分钟。作者借助华盛顿州交通中心提供的数据，能够研究这一现象。我们也做同样的研究。

我们应用在前几章介绍的概念，从总体问题开始，为什么我的公交车总是迟到？并将这个问题细化为更接近我们目标且可以通过数据调查的问题。然后我们考虑数据的范围，例如这些数据是如何收集的以及潜在的偏倚来源，并准备数据进行分析。我们对数据范围的理解帮助我们设计一个等车模型，我们可以模拟来研究这一现象。

问题和范围

我们的研究起源于一个经常乘坐公交车的乘客的经历，他想知道为什么他的公交车总是迟到。我们不是在寻找实际导致公交车迟到的原因，比如交通拥堵或维护延误。相反，我们想研究实际到达时间与计划到达时间之间的差异模式。这些信息将帮助我们更好地理解等待公交车的感受。

公交线路在世界各地甚至城市内部都有所不同，因此我们将调查范围缩小到西雅图市的一个公交车站。我们的数据是关于西雅图快速通行线路 C、D 和 E 在第三大道和派克街的停靠时间，华盛顿州交通中心提供了 2016 年 3 月 26 日至 5 月 27 日这三条公交线路所有实际和计划停靠时间的数据。

考虑到我们狭窄的范围仅限于两个月内一个特定停靠点的公交车，并且我们可以访问在这段时间窗口内收集的所有行政数据，人口、访问框架和样本是相同的。然而，我们可以想象我们的分析可能对西雅图及其以外的其他地点和其他时段也有用。如果我们幸运的话，我们发现的想法或采取的方法可能对他人有用。目前，我们保持狭窄的焦点。

让我们来看看这些数据，以更好地理解它们的结构。

数据整理

在开始分析之前，我们检查数据的质量，在可能的情况下简化结构，并导出可能有助于分析的新测量数据。我们在第九章中涵盖了这些操作类型，因此现在不必担心代码的细节。而是专注于清理数据时数据表之间的差异。我们首先将数据加载到 Python 中。

数据表的前几行显示如下：

`bus``.``head``(``3``)`

	OPD_DATE	VEHICLE_ID	RTE	DIR	...	STOP_ID	STOP_NAME	SCH_STOP_TM	ACT_STOP_TM
0	2016-03-26	6201	673	S	...	431	第三大道与派克街 (431)	01:11:57	01:13:19
1	2016-03-26	6201	673	S	...	431	第三大道与派克街 (431)	23:19:57	23:16:13
2	2016-03-26	6201	673	S	...	431	第三大道与派克街 (431)	21:19:57	21:18:46

3 rows × 9 columns

（原始数据以逗号分隔值的形式存储在文件中，我们已将其加载到此表中；详情请参见第八章。）

看起来表格中一些列可能是多余的，比如标记为STOP_ID和STOP_NAME的列。我们可以查找唯一值的数量及其计数来确认这一点：

`bus``[``[``'``STOP_ID``'``,``'``STOP_NAME``'``]``]``.``value_counts``(``)`

STOP_ID  STOP_NAME              
578      3RD AVE & PIKE ST (578)    19599
431      3RD AVE & PIKE ST (431)    19318
dtype: int64

有两个名为 3RD AVE & PIKE ST 的停靠点名称。我们想知道它们是否与公交车的行驶方向有关，可以通过方向、停靠点 ID 和停靠点名称的可能组合来检查：

`bus``[``[``'``DIR``'``,``'``STOP_ID``'``,``'``STOP_NAME``'``]``]``.``value_counts``(``)`

DIR  STOP_ID  STOP_NAME              
N    578      3RD AVE & PIKE ST (578)    19599
S    431      3RD AVE & PIKE ST (431)    19318
dtype: int64

实际上，北方向对应的是停靠点 ID 578，南方向对应的是停靠点 ID 431。由于我们只研究一个停靠点，我们实际上不需要比方向更多的信息。

我们还可以检查唯一路线名称的数量：

673    13228
674    13179
675    12510
Name: RTE, dtype: int64

这些路线是按编号排列，并不符合问题描述原文中的 C、D 和 E 的名称。这个问题涉及到数据整理的另一个方面：我们需要找到连接路线字母和数字的信息。我们可以从西雅图公共交通网站获取这些信息。数据整理的另一个部分是将数据值翻译成更容易理解的形式，因此我们用字母替换路线号码：

`def` `clean_stops``(``bus``)``:`
    `return` `bus``.``assign``(`
        `route``=``bus``[``"``RTE``"``]``.``replace``(``{``673``:` `"``C``"``,` `674``:` `"``D``"``,` `675``:` `"``E``"``}``)``,`
        `direction``=``bus``[``"``DIR``"``]``.``replace``(``{``"``N``"``:` `"``northbound``"``,` `"``S``"``:` `"``southbound``"``}``)``,`
    `)`

我们还可以在表格中创建新的列，帮助我们进行调查。例如，我们可以使用计划和实际到达时间计算公交车的晚点时间。这需要一些日期和时间格式的处理，这在第九章有介绍。

让我们检查这个新量的值，确保我们的计算是正确的：

 smallest amount late: -12.87 minutes
 greatest amount late:  150.28 minutes
 median amount late:    0.52 minutes

公交车晚点有负值有点令人惊讶，但这只是意味着公交车比计划时间提前到达。虽然中位数的晚点只有大约半分钟，但有些公交车晚到 2.5 小时！让我们看一下公交车晚点分钟数的直方图：

`px``.``histogram``(``bus``,` `x``=``"``minutes_late``"``,` `nbins``=``120``,` `width``=``450``,` `height``=``300``,` 
             `labels``=``{``'``minutes_late``'``:``'``Minutes late``'``}``)`

我们在第四章也看到了类似形状的直方图。公交车晚点情况的分布高度右偏，但许多车辆准时到达。

最后，我们通过创建数据表的简化版本来结束数据整理工作。因为我们只需要追踪路线、方向、计划和实际到达时间以及公交车的晚点时间，所以我们创建了一个更小的表格，并给列取了一些更易读的名称：

`bus` `=` `bus``[``[``"``route``"``,` `"``direction``"``,` `"``scheduled``"``,` `"``actual``"``,` `"``minutes_late``"``]``]`
`bus``.``head``(``)`

	路线	方向	计划时间	实际时间	迟到分钟
0	C	南行	2016-03-26 01:11:57	2016-03-26 01:13:19	1.37
1	C	南行	2016-03-26 23:19:57	2016-03-26 23:16:13	-3.73
2	C	南行	2016-03-26 21:19:57	2016-03-26 21:18:46	-1.18
3	C	南行	2016-03-26 19:04:57	2016-03-26 19:01:49	-3.13
4	C	南行	2016-03-26 16:42:57	2016-03-26 16:42:39	-0.30

这些表格操作在第六章有详细介绍。

在我们开始建模公交车晚点之前，我们想要深入探索和学习这些数据。我们接下来会做这件事。

探索公交时间

在我们清理和简化数据的过程中，我们对数据有了很多了解，但在开始建模等待时间之前，我们希望深入挖掘，更好地理解公交车晚点现象。我们将焦点缩小到了一个站点（第三大道和派克街口）在两个月内的公交活动上。我们看到公交车晚点的分布呈右偏态，确实有些公交车非常晚。在这个探索阶段，我们可能会问：

• 不同的三条公交线路的晚点分布看起来一样吗？

• 公交车是往北行驶还是往南行驶是否重要？

• 白天时间如何影响公交车的晚点情况？

• 公交车是否按照全天的规律间隔到达？

回答这些问题有助于我们更好地确定建模方法。

回顾一下 第四章 中我们发现的公交车晚点的中位数时间是 3/4 分钟。但这与我们为所有公交线路和方向计算的中位数（1/2 分钟）不符。我们来检查一下是否是由于该章节关注的是北行 C 线路的原因。我们创建每个六种公交线路和方向组合的延误直方图来解决这个问题和我们列表上的前两个问题：

y 轴上的比例尺（或密度）使得比较直方图变得更加容易，因为我们不会被不同组中的计数误导。x 轴上的范围在六个图中是相同的，这样更容易检测到分布的不同中心和扩展。（这些概念在 第十一章 中有描述。）

每条线路的南行和北行分布都不同。当我们深入了解背景时，我们了解到 C 线路起源于北部，而其他两条线路起源于南部。直方图暗示在公交路线的后半段到达时间的变异性较大，这是合理的，因为随着一天的进展，延误会逐渐累积。

接下来，为了探索不同时段的延误情况，我们需要推导一个新的量：公交车计划到达的小时。鉴于我们刚刚看到的路线和方向的变化导致的公交车晚点情况，我们再次为每条线路和方向创建单独的图表：

的确，似乎存在交通高峰时间的影响，而且晚高峰似乎比早高峰更严重。北行 C 线路看起来受到的影响最大。

最后，为了检查公交车的计划频率，我们需要计算计划到达时间之间的间隔。我们在表格中创建一个新列，其中包含北行 C 线路公交车的计划到达时间之间的时间：

`minute` `=` `pd``.``Timedelta``(``'``1 minute``'``)`
`bus_c_n` `=` `(`
    `bus``[``(``bus``[``'``route``'``]` `==` `'``C``'``)` `&` `(``bus``[``'``direction``'``]` `==` `'``northbound``'``)``]`
    `.``sort_values``(``'``scheduled``'``)`
    `.``assign``(``sched_inter``=``lambda` `x``:` `x``[``'``scheduled``'``]``.``diff``(``)` `/` `minute``)`
`)`
`bus_c_n``.``head``(``3``)`

	线路	方向	计划时间	实际时间	延迟分钟数	计划间隔
19512	C	北行	2016-03-26 00:00:25	2016-03-26 00:05:01	4.60	NaN
19471	C	北行	2016-03-26 00:30:25	2016-03-26 00:30:19	-0.10	30.0
19487	C	北行	2016-03-26 01:05:25	2016-03-26 01:10:15	4.83	35.0

让我们来看一下这些公交车的到站间隔时间分布直方图：

`fig` `=` `px``.``histogram``(``bus_c_n``,` `x``=``'``sched_inter``'``,` 
                   `title``=``"``Bus line C, northbound``"``,`
                   `width``=``450``,` `height``=``300``)`

`fig``.``update_xaxes``(``range``=``[``0``,` `40``]``,` `title``=``"``Time between consecutive buses``"``)`
`fig``.``update_layout``(``margin``=``dict``(``t``=``40``)``)`

我们可以看到，公交车在一天中的不同时间段被安排以不同的间隔到达。在这两个月的时间内，大约有 1500 辆公交车被安排在前一辆车后 12 分钟到达，大约有 1400 辆公交车被安排在前一辆车后 15 分钟到达。

在探索数据的过程中，我们学到了很多，并且现在更有能力拟合模型。尤其是，如果我们想要清楚地了解等待公交车的经历，我们需要考虑公交车预定到达之间的间隔，以及公交线路和方向。

建模等待时间

我们对建模等待公交车的人的经历很感兴趣。我们可以开发一个涉及预定到达间隔、公交线路和方向的复杂模型。相反，我们采取了更简单的方法，将焦点缩小到一个线路、一个方向和一个预定间隔。我们检查了预定间隔为 12 分钟的 C 线向北站点：

`bus_c_n_12` `=` `bus_c_n``[``bus_c_n``[``'``sched_inter``'``]` `==` `12``]`

复杂模型和狭窄方法都是合法的，但我们尚未掌握处理复杂模型的工具（请参阅第十五章了解更多建模细节）。

到目前为止，我们已经检查了公交车晚点的分钟数分布。我们为我们正在分析的数据子集（预定在前一辆公交车到达后 12 分钟的 C 线向北站点）创建了另一个此延迟的直方图：

`fig` `=` `px``.``histogram``(``bus_c_n_12``,` `x``=``'``minutes_late``'``,`
                   `labels``=``{``'``minutes_late``'``:``'``Minutes late``'``}``,`
                   `nbins``=``120``,` `width``=``450``,` `height``=``300``)`

`fig``.``add_annotation``(``x``=``20``,` `y``=``150``,`  `showarrow``=``False``,`
  `text``=``"``Line C, northbound<br>Scheduled arrivals: 12 minutes apart``"` `)`
`fig``.``update_xaxes``(``range``=``[``-``13``,` `40``]``)`
`fig``.``show``(``)`

现在让我们计算晚点的最小值、最大值和中位数：

smallest amount late:  -10.20 minutes
greatest amount late:  57.00 minutes
median amount late:    -0.50 minutes

有趣的是，C 线向北行驶的公交车，间隔 12 分钟的情况下，更多时候会提前抵达而不是晚点！

现在让我们重新审视我们的问题，确认我们正在答案的正确方向上。仅总结公交车晚点的情况并不能完全解释等待公交车的人的经历。当有人到达公交车站时，他们需要等待下一辆公交车到达。图 5-1 展示了乘客和公交车到达公交车站时时间流逝的理想化图像。如果人们在随机时间到达公交车站，注意到他们更有可能在公交车晚点的时间段到达，因为公交车之间的间隔更长。这种到达模式是大小偏倚抽样的一个例子。因此，要回答等待公交车的人们经历了什么这个问题，我们需要做的不仅仅是总结公交车晚点的情况。

图 5-1. 理想化时间轴，显示公交车到达（矩形）、乘客到达（圆圈）和乘客等待下一辆公交车到达的时间（大括号）

我们可以设计一个模拟，模拟一天中等待公交车的过程，使用来自第三章的思想。为此，我们设置了一系列从早上 6 点到午夜预定间隔为 12 分钟的公交车到达时间：

`scheduled` `=` `12` `*` `np``.``arange``(``91``)`
`scheduled`

array([   0,   12,   24, ..., 1056, 1068, 1080])

针对每个预定到达时间，我们通过添加每辆公交车晚点的随机分钟数来模拟其实际到达时间。为此，我们从实际公交车晚点的分布中选择晚点的分钟数。注意，在我们的模拟研究中，我们已经通过使用实际公交车延迟的分布来整合真实数据，这些公交车的间隔为 12 分钟：

`minutes_late` `=` `bus_c_n_12``[``'``minutes_late``'``]`
`actual` `=` `scheduled` `+` `np``.``random``.``choice``(``minutes_late``,` `size``=``91``,` `replace``=``True``)`

我们需要对这些到达时间进行排序，因为当一辆公交车迟到很久时，可能会有另一辆公交车在它之前到达：

`actual``.``sort``(``)`
`actual`

array([  -1.2 ,   25.37,   32.2 , ..., 1051.02, 1077\.  , 1089.43])

我们还需要随机模拟人们在一天中的不同时间到达公交车站的情况。我们可以使用另一个不同的罐子模型来描述乘客的到达情况。对于乘客，我们在罐子里放入带有时间的彩球。这些时间从时间 0 开始，代表上午 6 点，到午夜最后一辆公交车的到达，即从上午 6 点起计算到夜间 1,068 分钟。为了与我们数据中的公交车时间测量方式匹配，我们将这些时间分成每分钟 1/100 的间隔：

`pass_arrival_times` `=` `np``.``arange``(``100``*``1068``)`
`pass_arrival_times` `/` `100`

array([   0\.  ,    0.01,    0.02, ..., 1067.97, 1067.98, 1067.99])

现在我们可以模拟一天中，例如，五百人在公交车站的到达情况。我们从这个罐子中抽取五百次，每次抽取后替换掉彩球：

`sim_arrival_times` `=` `(`
    `np``.``random``.``choice``(``pass_arrival_times``,` `size``=``500``,` `replace``=``True``)` `/` `100`
`)`
`sim_arrival_times``.``sort``(``)`
`sim_arrival_times`

array([   2.06,    3.01,    8.54, ..., 1064\.  , 1064.77, 1066.42])

要了解每个人等待多长时间，我们寻找他们采样时间后最快到达的公交车。这两个时间的差值（个人采样时间和其后最快公交车到达时间）就是这个人的等待时间：

`i` `=` `np``.``searchsorted``(``actual``,` `sim_arrival_times``,` `side``=``'``right``'``)`
`sim_wait_times` `=` `actual``[``i``]` `-` `sim_arrival_times`
`sim_wait_times`

array([23.31, 22.36, 16.83, ..., 13\.  , 12.23, 10.58])

我们可以建立一个完整的模拟，例如模拟两百天的公交车到达，而每天我们模拟五百人在一天中的随机时间到达公交车站。总计是 100,000 次模拟等待时间：

`sim_wait_times` `=` `[``]`

`for` `day` `in` `np``.``arange``(``0``,` `200``,` `1``)``:`
    `bus_late` `=` `np``.``random``.``choice``(``minutes_late``,` `size``=``91``,` `replace``=``True``)`
    `actual` `=` `scheduled` `+` `bus_late`
    `actual``.``sort``(``)`
    `sim_arrival_times` `=` `(`
        `np``.``random``.``choice``(``pass_arrival_times``,` `size``=``500``,` `replace``=``True``)` `/` `100`
    `)`
    `sim_arrival_times``.``sort``(``)`
    `i` `=` `np``.``searchsorted``(``actual``,` `sim_arrival_times``,` `side``=``"``right``"``)`
    `sim_wait_times` `=` `np``.``append``(``sim_wait_times``,` `actual``[``i``]` `-` `sim_arrival_times``)`

让我们制作这些模拟等待时间的直方图，以检查其分布情况：

`fig` `=` `px``.``histogram``(``x``=``sim_wait_times``,` `nbins``=``40``,`
                   `histnorm``=``'``probability density``'``,`
                   `width``=``450``,` `height``=``300``)`

`fig``.``update_xaxes``(``title``=``"``Simulated wait times for 100,000 passengers``"``)`
`fig``.``update_yaxes``(``title``=``"``proportion``"``)`
`fig``.``show``(``)`

正如我们预期的那样，我们发现了一个偏斜的分布。我们可以用一个常数模型来描述这一点，其中我们使用绝对损失来选择最佳常数。我们在第四章中看到，绝对损失给出了中位数等待时间：

`print``(``f``"``Median wait time:` `{``np``.``median``(``sim_wait_times``)``:``.2f``}` `minutes``"``)`

Median wait time: 6.49 minutes

中位数约为六分半钟，看起来并不太长。虽然我们的模型捕捉了典型等待时间，但我们也想提供一个过程变异性的估计。这个话题在第十七章中有所涉及。我们可以计算等待时间的上四分位数，以帮助我们了解过程的变异性：

`print``(``f``"``Upper quartile:` `{``np``.``quantile``(``sim_wait_times``,` `0.75``)``:``.2f``}` `minutes``"``)`

Upper quartile: 10.62 minutes

上四分位数相当大。当你等待超过 10 分钟的公交车时，这无疑是令人难忘的，因为公交车本应每 12 分钟到达一次，但每四次乘车中就有一次会发生这种情况！

摘要

在我们的第一个案例研究中，我们已经遍历了数据建模的整个生命周期。也许你会觉得，这样一个简单的问题并不能立即用收集的数据回答。我们需要将公交车的预定到达时间和实际到达时间的数据与乘客在随机时间到达公交车站的模拟研究结合起来，以揭示乘客的等待体验。

这个模拟简化了公交乘车中的许多真实模式。我们关注的是单向行驶的一条公交线路，公交车每隔 12 分钟到达一次。此外，数据的探索显示，迟到的模式与一天中的时间相关，这在我们的分析中尚未考虑。尽管如此，我们的发现仍然是有用的。例如，它们证实了典型的等待时间长于计划间隔的一半。等待时间的分布具有右长尾，意味着乘客的体验可能受到过程中变化的影响。

我们还看到了如何衍生新的量，例如公交车的迟到时间和公交车之间的时间，并探索数据在建模中的实用性。我们的直方图显示，特定的公交线路和方向是重要的，需要加以考虑。我们还发现，随着一天中时间的变化，许多公交车在另一辆车到达后 10、12 和 15 分钟到达，有些则到达频率更高或更分散。这一观察进一步为建模阶段提供了信息。

最后，我们使用了pandas和plotly等数据工具库，这些将在后面的章节中介绍。我们的重点不在于如何操作表格或创建图表，而是专注于生命周期，将问题与数据连接到建模到结论。在下一章中，我们将转向处理数据表格的实际问题。

第二部分：矩形数据

第六章：使用 pandas 处理数据框

数据科学家处理存储在表格中的数据。本章介绍了“数据框架”，这是表示数据表的最常用方式之一。我们还介绍了pandas，这是处理数据框架的标准 Python 包。以下是一个包含有关流行狗品种信息的数据框架示例：

	grooming	food_cost	kids	size
breed
---	---	---	---	---
拉布拉多寻回犬	每周	466.0	高	中型
德国牧羊犬	每周	466.0	中等	大型
比格犬	每天	324.0	高	小型
金毛寻回犬	每周	466.0	高	中型
约克夏梗	每天	324.0	低	小型
英国斗牛犬	每周	466.0	中等	中型
拳师犬	每周	466.0	高	中型

在数据框中，每行代表一个单独的记录——在本例中是一个狗品种。每列代表记录的一个特征——例如，grooming 列表示每个狗品种需要多频繁地梳理。

数据框架具有列和行的标签。例如，此数据框架具有一个标记为 grooming 的列和一个标记为德国牧羊犬的行。数据框架的列和行是有序的——我们可以将拉布拉多寻回犬行称为数据框架的第一行。

在列内，数据具有相同的类型。例如，食物成本包含数字，狗的大小由类别组成。但是在行内，数据类型可以不同。

由于这些属性，数据框架使得各种有用的操作成为可能。

注意

数据科学家经常发现自己与使用不同背景的人合作。例如，计算机科学家说数据框中的列表示数据的“特征”，而统计学家则称之为“变量”。

有时，人们使用同一个术语来指代略有不同的概念。在编程中，“数据类型”指的是计算机如何在内部存储数据。例如，Python 中的size列具有字符串数据类型。但从统计学的角度来看，size列的类型是有序分类数据（序数数据）。我们在第十章中详细讨论了这种具体区别。

在本章中，我们介绍了常见的数据框架操作。数据科学家在 Python 中处理数据框架时使用pandas库。首先，我们解释了pandas提供的主要对象：DataFrame和Series类。然后，我们展示如何使用pandas执行常见的数据操作任务，如切片、过滤、排序、分组和连接。

子集

本节介绍了对数据框进行子集操作的操作。当数据科学家首次读取数据框时，他们通常希望获取计划使用的特定数据子集。例如，数据科学家可以从包含数百列的数据框中切片出 10 个相关特征。或者，他们可以过滤数据框以删除包含不完整数据的行。在本章的其余部分，我们使用一个婴儿姓名的数据框来演示数据框操作。

数据范围和问题

有一篇2021 年《纽约时报》文章讨论了哈里王子和梅根·马克尔为他们的新生女儿选择的独特姓名“莉莉贝特”。文章中采访了婴儿姓名专家帕梅拉·雷德蒙德，她谈到了人们如何给孩子取名的有趣趋势。例如，她说以字母“L”开头的名字近年来变得非常流行，而以字母“J”开头的名字在上世纪 70 年代和 80 年代最受欢迎。这些说法在数据中反映出来吗？我们可以使用pandas来找出答案。

首先，我们将包导入为pd，这是它的常用缩写：

`import` `pandas` `as` `pd`

我们有一个婴儿姓名数据集，存储在名为babynames.csv的逗号分隔值（CSV）文件中。我们使用pd.read_csv函数将文件读取为pandas.DataFrame对象：

`baby` `=` `pd``.``read_csv``(``'``babynames.csv``'``)`
`baby`

	名称	性别	计数	年份
0	利亚姆	M	19659	2020
1	诺亚	M	18252	2020
2	奥利弗	M	14147	2020
...	...	...	...	...
2020719	维罗纳	F	5	1880
2020720	维尔蒂	F	5	1880
2020721	威尔玛	F	5	1880

2020722 rows × 4 columns

baby表中的数据来自美国社会保障局（SSA），记录了出生证明目的的婴儿姓名和出生性别。SSA 将婴儿姓名数据提供在其网站上。我们已将此数据加载到baby表中。

SSA 网站有一个页面，详细描述了数据。在本章中，我们不会深入讨论数据的限制，但我们将指出网站上有关此相关信息：

所有姓名均来自于 1879 年后在美国出生的社会保障卡申请。请注意，许多 1937 年之前出生的人从未申请过社会保障卡，因此他们的姓名不包含在我们的数据中。对于那些申请过的人，我们的记录可能不显示出生地，同样，他们的姓名也不包含在我们的数据中。

所有数据都来自我们 2021 年 3 月社会保障卡申请记录的 100％样本。

在撰写本文时，重要的一点是指出，SSA 数据集仅提供了男性和女性的二进制选项。我们希望将来像这样的国家数据集能够提供更多包容性选项。

数据框和索引

让我们更详细地查看 baby 数据帧。数据帧有行和列。每行和列都有一个标签，如 图 6-1 中所示。

图 6-1. baby 数据帧为行和列都设置了标签（框起来的）

默认情况下，pandas 分配的行标签从 0 开始递增。在这种情况下，标记为 0 的行和标记为 Name 的列的数据为 'Liam'。

数据帧的行标签也可以是字符串。图 6-2 显示了一个狗数据的数据帧，其中行标签是字符串。

图 6-2. 数据帧中的行标签也可以是字符串，例如本例中，每行都使用狗品种名称进行标记

行标签有一个特殊的名称。我们称之为数据帧的 索引，pandas 将行标签存储在特殊的 pd.Index 对象中。我们暂时不讨论 pd.Index 对象，因为不常操作索引本身。但现在重要的是要记住，即使索引看起来像数据的一列，索引实际上代表行标签，而不是数据。例如，狗品种的数据帧有四列数据，而不是五列，因为索引不算作列。

切片

切片 是通过从另一个数据帧中取出部分行或列来创建一个新数据帧的操作。想象一下切西红柿——切片可以在垂直和水平方向上进行。在 pandas 中进行数据帧切片时，我们使用 .loc 和 .iloc 属性。让我们从 .loc 开始。

这是完整的 baby 数据帧：

`baby`

	Name	Sex	Count	Year
0	Liam	M	19659	2020
1	Noah	M	18252	2020
2	Oliver	M	14147	2020
...	...	...	...	...
2020719	Verona	F	5	1880
2020720	Vertie	F	5	1880
2020721	Wilma	F	5	1880

2020722 rows × 4 columns

.loc 允许我们使用它们的标签选择行和列。例如，要获取标记为 1 的行和标记为 Name 的列中的数据：

`#        The first argument is the row label`
`#        ↓`
`baby``.``loc``[``1``,` `'``Name``'``]`
`#            ↑`
`#            The second argument is the column label`

'Noah'

警告

注意，.loc 需要使用方括号；运行 baby.loc(1, 'Name') 将导致错误。

为了切片多行或列，我们可以使用 Python 切片语法，而不是单独的值：

`baby``.``loc``[``0``:``3``,` `'``Name``'``:``'``Count``'``]`

	Name	Sex	Count
0	Liam	M	19659
1	Noah	M	18252
2	Oliver	M	14147
3	Elijah	M	13034

要获取整列数据，我们可以将空切片作为第一个参数传递：

`baby``.``loc``[``:``,` `'``Count``'``]`

0          19659
1          18252
2          14147
           ...  
2020719        5
2020720        5
2020721        5
Name: Count, Length: 2020722, dtype: int64

请注意，这个输出看起来不像一个数据帧，因为它不是。选择数据帧的单行或单列会产生一个 pd.Series 对象：

`counts` `=` `baby``.``loc``[``:``,` `'``Count``'``]`
`counts``.``__class__``.``__name__`

'Series'

什么是pd.Series对象和pd.DataFrame对象的区别？本质上，pd.DataFrame是二维的——它有行和列，代表着数据表。pd.Series是一维的——它代表着数据列表。pd.Series和pd.DataFrame对象有许多共同的方法，但它们实际上代表着两种不同的东西。混淆两者可能会导致错误和混乱。

要选择数据帧的特定列，请将列表传递给.loc。以下是原始数据帧：

`baby`

	名字	性别	数量	年份
0	Liam	M	19659	2020
1	Noah	M	18252	2020
2	Oliver	M	14147	2020
...	...	...	...	...
2020719	Verona	F	5	1880
2020720	Vertie	F	5	1880
2020721	Wilma	F	5	1880

2020722 rows × 4 columns

`# And here's the dataframe with only Name and Year columns`
`baby``.``loc``[``:``,` `[``'``Name``'``,` `'``Year``'``]``]`
`#           └──────┬───────┘`
`#         list of column labels`

	名字	年份
0	Liam	2020
1	Noah	2020
2	Oliver	2020
...	...	...
2020719	Verona	1880
2020720	Vertie	1880
2020721	Wilma	1880

2020722 rows × 2 columns

选择列非常常见，所以有一种简便写法：

`# Shorthand for baby.loc[:, 'Name']`
`baby``[``'``Name``'``]`

0            Liam
1            Noah
2          Oliver
            ...  
2020719    Verona
2020720    Vertie
2020721     Wilma
Name: Name, Length: 2020722, dtype: object

`# Shorthand for baby.loc[:, ['Name', 'Count']]`
`baby``[``[``'``Name``'``,` `'``Count``'``]``]`

	名字	数量
0	Liam	19659
1	Noah	18252
2	Oliver	14147
...	...	...
2020719	Verona	5
2020720	Vertie	5
2020721	Wilma	5

2020722 rows × 2 columns

使用.iloc切片与.loc类似，不同之处在于.iloc使用行和列的位置而不是标签。当数据帧索引具有字符串时，演示时最容易显示.iloc和.loc之间的差异，所以让我们看一个有关狗品种信息的数据帧：

`dogs` `=` `pd``.``read_csv``(``'``dogs.csv``'``,` `index_col``=``'``breed``'``)`
`dogs`

	美容	食品成本	孩子	尺寸
品种
---	---	---	---	---
拉布拉多猎犬	每周	466.0	高	中等
德国牧羊犬	每周	466.0	中等	大型
猎犬	每日	324.0	高	小型
金毛寻回犬	每周	466.0	高	中等
约克夏梗	每日	324.0	低	小型
斗牛犬	每周	466.0	中等	中等
拳师犬	每周	466.0	高	中等

要通过位置获取前三行和前两列，请使用.iloc：

`dogs``.``iloc``[``0``:``3``,` `0``:``2``]`

	美容	食品成本
品种
---	---	---
拉布拉多猎犬	每周	466.0
德国牧羊犬	每周	466.0
猎犬	每日	324.0

使用.loc进行相同操作需要使用数据帧标签：

`dogs``.``loc``[``'``Labrador Retriever``'``:``'``Beagle``'``,` `'``grooming``'``:``'``food_cost``'``]`

	美容	食品成本
品种
---	---	---
拉布拉多猎犬	每周	466.0
德国牧羊犬	每周	466.0
猎犬	每日	324.0

接下来，我们将看看如何过滤行。

过滤行

到目前为止，我们已经展示了如何使用.loc和.iloc使用标签和位置来切片数据帧。

然而，数据科学家通常希望筛选行——他们希望使用某些条件获取行的子集。假设我们要找到 2020 年最流行的婴儿名字。为此，我们可以筛选行，仅保留Year为 2020 的行。

要筛选，我们想要检查Year列中的每个值是否等于 1970，然后仅保留这些行。

要比较Year中的每个值，我们切出列并进行布尔比较（这与我们在numpy数组中所做的类似）。以下是参考数据框：

`baby`

	名称	性别	计数	年份
0	Liam	M	19659	2020
1	Noah	M	18252	2020
2	Oliver	M	14147	2020
...	...	...	...	...
2020719	Verona	F	5	1880
2020720	Vertie	F	5	1880
2020721	Wilma	F	5	1880

2020722 rows × 4 columns

`# Get a Series with the Year data`
`baby``[``'``Year``'``]`

0          2020
1          2020
2          2020
           ... 
2020719    1880
2020720    1880
2020721    1880
Name: Year, Length: 2020722, dtype: int64

`# Compare with 2020`
`baby``[``'``Year``'``]` `==` `2020`

0           True
1           True
2           True
           ...  
2020719    False
2020720    False
2020721    False
Name: Year, Length: 2020722, dtype: bool

注意，对Series进行布尔比较会生成布尔Series。这几乎等同于：

`is_2020` `=` `[``]`
`for` `value` `in` `baby``[``'``Year``'``]``:`
    `is_2020``.``append``(``value` `==` `2020``)`

但布尔比较比for循环更容易编写，执行起来也更快。

现在告诉pandas仅保留评估为True的行：

`baby``.``loc``[``baby``[``'``Year``'``]` `==` `2020``,` `:``]`

	名称	性别	计数	年份
0	Liam	M	19659	2020
1	Noah	M	18252	2020
2	Oliver	M	14147	2020
...	...	...	...	...
31267	Zylynn	F	5	2020
31268	Zynique	F	5	2020
31269	Zynlee	F	5	2020

31270 rows × 4 columns

提示

将布尔Series传递到.loc中，仅保留Series具有True值的行。

筛选有简写方式。这将计算与前面代码段相同的表，但不使用.loc：

`baby``[``baby``[``'``Year``'``]` `==` `2020``]`

	名称	性别	计数	年份
0	Liam	M	19659	2020
1	Noah	M	18252	2020
2	Oliver	M	14147	2020
...	...	...	...	...
31267	Zylynn	F	5	2020
31268	Zynique	F	5	2020
31269	Zynlee	F	5	2020

31270 rows × 4 columns

最后，为了找到 2020 年最常见的名字，请按降序排列数据框中的Count。将较长的表达式用括号括起来，可以轻松添加换行符，使其更易读：

`(``baby``[``baby``[``'``Year``'``]` `==` `2020``]`
 `.``sort_values``(``'``Count``'``,` `ascending``=``False``)`
 `.``head``(``7``)` `# take the first seven rows`
`)`

	名称	性别	计数	年份
0	Liam	M	19659	2020
1	Noah	M	18252	2020
13911	Emma	F	15581	2020
2	Oliver	M	14147	2020
13912	Ava	F	13084	2020
3	Elijah	M	13034	2020
13913	Charlotte	F	13003	2020

我们看到，Liam、Noah 和 Emma 是 2020 年最受欢迎的婴儿名字。

示例：Luna 最近成为热门名字吗？

纽约时报的文章提到，Luna 这个名字在 2000 年之前几乎不存在，但此后已成为女孩们非常流行的名字。Luna 究竟在何时变得流行？我们可以使用切片和筛选来检查。在解决数据处理任务时，我们建议将问题分解为较小的步骤。例如，我们可以这样思考：

1. 筛选：仅保留Name列中含有'Luna'的行。

2. 筛选：仅保留Sex列中含有'F'的行。

3. 切片：保留Count和Year列。

现在只需要将每个步骤翻译成代码即可：

`luna` `=` `baby``[``baby``[``'``Name``'``]` `==` `'``Luna``'``]` `# [1]`
`luna` `=` `luna``[``luna``[``'``Sex``'``]` `==` `'``F``'``]`     `# [2]`
`luna` `=` `luna``[``[``'``Count``'``,` `'``Year``'``]``]`      `# [3]`
`luna`

	计数	年份
13923	7770	2020
45366	7772	2019
77393	6929	2018
...	...	...
2014083	17	1883
2018187	18	1881
2020223	15	1880

128 rows × 2 columns

在本书中，我们使用一个叫做plotly的库进行绘图。我们不会在这里深入讨论绘图，因为我们在第十一章中更多地谈论了它。现在，我们使用px.line()制作一个简单的折线图：

`px``.``line``(``luna``,` `x``=``'``Year``'``,` `y``=``'``Count``'``,` `width``=``350``,` `height``=``250``)`

就像文章所说的那样。Luna 在 2000 年左右几乎不流行。换句话说，即使没有关于他们的其他任何信息，如果有人告诉你他们的名字是 Luna，你也可以很好地猜到他们的年龄！

纯属娱乐，这里是同样的 Siri 名字的图表：

`siri` `=` `(``baby``.``query``(``'``Name ==` `"``Siri``"``'``)`
        `.``query``(``'``Sex ==` `"``F``"``'``)``)`
`px``.``line``(``siri``,` `x``=``'``Year``'``,` `y``=``'``Count``'``,` `width``=``350``,` `height``=``250``)`

小贴士

使用.query类似于使用带有布尔系列的.loc。query()在过滤上有更多的限制，但可以作为一种简写方便使用。

为什么在 2010 年后突然变得不那么受欢迎呢？嗯，Siri 恰好是苹果的语音助手的名字，于 2011 年推出。让我们在 2011 年划一条线，看看：

`fig` `=` `px``.``line``(``siri``,` `x``=``"``Year``"``,` `y``=``"``Count``"``,` `width``=``350``,` `height``=``250``)`
`fig``.``add_vline``(`
    `x``=``2011``,` `line_color``=``"``red``"``,` `line_dash``=``"``dashdot``"``,` `line_width``=``4``,` `opacity``=``0.7`
`)`

看起来家长们不希望其他人对他们的手机说“嘿 Siri”时让他们的孩子感到困惑。

在这一节中，我们介绍了pandas中的数据框。我们涵盖了数据科学家对数据框进行子集切片和使用布尔条件进行筛选的常见方式。在下一节中，我们将解释如何将行聚合在一起。

聚合

本节介绍了数据框中行聚合的操作。数据科学家将行聚合在一起以对数据进行摘要。例如，包含每日销售额的数据集可以聚合以显示月销售额。本节介绍了分组和透视，这两种常见的聚合数据操作。

我们使用上一节介绍的婴儿姓名数据：

`baby` `=` `pd``.``read_csv``(``'``babynames.csv``'``)`
`baby`

	名字	性别	计数	年份
0	Liam	男	19659	2020
1	Noah	男	18252	2020
2	Oliver	男	14147	2020
...	...	...	...	...
2020719	Verona	女	5	1880
2020720	Vertie	女	5	1880
2020721	Wilma	女	5	1880

2020722 rows × 4 columns

基本的分组聚合

假设我们想要找出记录在案的总出生婴儿数。这只是 Count 列的总和：

`baby``[``'``Count``'``]``.``sum``(``)`

352554503

汇总姓名计数是一种简单的聚合数据的方式——它将多行的数据合并。

但假设我们想回答一个更有趣的问题：美国出生率是否随时间而上升？为了回答这个问题，我们可以在每年内对 Count 列求和，而不是在整个数据集上进行求和。换句话说，我们根据 Year 将数据分组，然后在每个组内对 Count 值求和。这个过程在 图 6-3 中有所描述。

图 6-3. 描述了示例数据的分组和聚合过程

我们将这个操作称为分组，然后是聚合。在pandas中，我们这样写：

`baby``.``groupby``(``'``Year``'``)``[``'``Count``'``]``.``sum``(``)`

Year
1880     194419
1881     185772
1882     213385
         ...   
2018    3487193
2019    3437438
2020    3287724
Name: Count, Length: 141, dtype: int64

注意，代码几乎与未分组版本相同，只是以 .groupby('Year') 开始。

结果是一个 pd.Series，其中包含数据集中每年出生的总婴儿数。请注意，该系列的索引包含唯一的 Year 值。现在我们可以绘制随时间变化的计数：

`counts_by_year` `=` `baby``.``groupby``(``'``Year``'``)``[``'``Count``'``]``.``sum``(``)``.``reset_index``(``)`
`px``.``line``(``counts_by_year``,` `x``=``'``Year``'``,` `y``=``'``Count``'``,` `width``=``350``,` `height``=``250``)`

在这张图中我们看到了什么？首先，我们注意到 1920 年之前出生的婴儿似乎非常少。一个可能的解释是社会保障局成立于 1935 年，因此其之前的出生数据可能不够完整。

我们还注意到了 1939 年第二次世界大战爆发时的下降以及 1946 年至 1964 年间的战后婴儿潮时期。

这是在pandas中进行分组的基本步骤：

`(``baby`                `# the dataframe`
 `.``groupby``(``'``Year``'``)`    `# column(s) to group`
 `[``'``Count``'``]`           `# column(s) to aggregate`
 `.``sum``(``)`              `# how to aggregate`
`)`

示例：使用.value_counts()

数据框中更常见的任务之一是计算列中每个唯一项出现的次数。例如，我们可能对以下 classroom 数据框中每个姓名出现的次数感兴趣：

`classroom`

	姓名
0	Eden
1	Sachit
2	Eden
3	Sachit
4	Sachit
5	Luke

实现这一点的一种方式是使用我们的分组步骤和.size() 聚合函数：

`(``classroom`
 `.``groupby``(``'``name``'``)`
 `[``'``name``'``]`
 `.``size``(``)`
`)`

name
Eden      2
Luke      1
Sachit    3
Name: name, dtype: int64

这个操作非常常见，pandas提供了一个简写——.value_counts() 方法用于 pd.Series 对象：

`classroom``[``'``name``'``]``.``value_counts``(``)`

name
Sachit    3
Eden      2
Luke      1
Name: count, dtype: int64

默认情况下，.value_counts() 方法会对结果系列按照从最高到最低的顺序进行排序，方便查看最常见和最不常见的值。我们提到这个方法是因为在书的其他章节中经常使用它。

在多列上进行分组

我们可以将多个列作为列表传递给 .groupby 方法，以一次性对多列进行分组。当需要进一步细分我们的组时，这非常有用。例如，我们可以按年份和性别对数据进行分组，以查看随时间变化的男女婴儿出生数：

`counts_by_year_and_sex` `=` `(``baby`
 `.``groupby``(``[``'``Year``'``,` `'``Sex``'``]``)` `# Arg to groupby is a list of column names`
 `[``'``Count``'``]`
 `.``sum``(``)`
`)`
`counts_by_year_and_sex` 

Year  Sex
1880  F        83929
      M       110490
1881  F        85034
              ...   
2019  M      1785527
2020  F      1581301
      M      1706423
Name: Count, Length: 282, dtype: int64

注意代码如何紧随分组的步骤。

counts_by_year_and_sex 系列具有我们称之为多级索引的两个级别，一个用于每列进行的分组。如果我们将系列转换为数据框，则更容易看到结果只有一列：

`counts_by_year_and_sex``.``to_frame``(``)`

		数量
年份	性别
---	---	---
1880	F	83929
M	110490
1881	F	85034
...	...	...
2019	M	1785527
2020	F	1581301
M	1706423

282 rows × 1 columns

索引有两个级别，因为我们按两列进行了分组。多级索引可能有点棘手，所以我们可以重置索引，回到具有单个索引的 dataframe：

`counts_by_year_and_sex``.``reset_index``(``)`

	年份	性别	数量
0	1880	F	83929
1	1880	M	110490
2	1881	F	85034
...	...	...	...
279	2019	M	1785527
280	2020	F	1581301
281	2020	M	1706423

282 rows × 3 columns

自定义聚合函数

分组后，pandas 为我们提供了灵活的方法来聚合数据。到目前为止，我们已经看到了如何在分组后使用 .sum()：

`(``baby`
 `.``groupby``(``'``Year``'``)`
 `[``'``Count``'``]`
 `.``sum``(``)` `# aggregate by summing`
`)`

Year
1880     194419
1881     185772
1882     213385
         ...   
2018    3487193
2019    3437438
2020    3287724
Name: Count, Length: 141, dtype: int64

pandas 还提供其他聚合函数，如 .mean()、.size() 和 .first()。下面是使用 .max() 进行相同分组的示例：

`(``baby`
 `.``groupby``(``'``Year``'``)`
 `[``'``Count``'``]`
 `.``max``(``)` `# aggregate by taking the max within each group`
`)`

Year
1880     9655
1881     8769
1882     9557
        ...  
2018    19924
2019    20555
2020    19659
Name: Count, Length: 141, dtype: int64

但有时 pandas 并没有我们想要使用的确切聚合函数。在这些情况下，我们可以定义并使用自定义聚合函数。pandas 通过 .agg(fn) 让我们能够做到这一点，其中 fn 是我们定义的函数。

例如，如果我们想要找出每个组中最大值和最小值之间的差异（数据的范围），我们可以首先定义一个名为 data_range 的函数，然后将该函数传递给 .agg()。这个函数的输入是一个包含一列数据的 pd.Series 对象。它会对每个组调用一次：

`def` `data_range``(``counts``)``:`
    `return` `counts``.``max``(``)` `-` `counts``.``min``(``)`

`(``baby`
 `.``groupby``(``'``Year``'``)`
 `[``'``Count``'``]`
 `.``agg``(``data_range``)` `# aggregate using custom function`
`)`

Year
1880     9650
1881     8764
1882     9552
        ...  
2018    19919
2019    20550
2020    19654
Name: Count, Length: 141, dtype: int64

我们首先定义了一个 count_unique 函数，用于计算系列中唯一值的数量。然后我们将该函数传递给 .agg()。由于这个函数很短，我们可以使用 lambda 表达式代替：

`def` `count_unique``(``s``)``:`
    `return` `len``(``s``.``unique``(``)``)`

`unique_names_by_year` `=` `(``baby`
 `.``groupby``(``'``Year``'``)`
 `[``'``Name``'``]`
 `.``agg``(``count_unique``)` `# aggregate using the custom count_unique function`
`)`
`unique_names_by_year`

Year
1880     1889
1881     1829
1882     2012
        ...  
2018    29619
2019    29417
2020    28613
Name: Name, Length: 141, dtype: int64

`px``.``line``(``unique_names_by_year``.``reset_index``(``)``,`
        `x``=``'``Year``'``,` `y``=``'``Name``'``,`
        `labels``=``{``'``Name``'``:` `'``# unique names``'``}``,`
        `width``=``350``,` `height``=``250``)`

我们发现，尽管自 1960 年代以来每年出生的婴儿数量已经趋于稳定，但独特名称的数量总体上还是在增加的。

数据透视

数据透视基本上是在使用两列进行分组时，方便地安排分组和聚合结果的一种方法。在本节的前面，我们将婴儿姓名数据按年份和性别分组：

`counts_by_year_and_sex` `=` `(``baby`
 `.``groupby``(``[``'``Year``'``,` `'``Sex``'``]``)` 
 `[``'``Count``'``]`
 `.``sum``(``)`
`)`
`counts_by_year_and_sex``.``to_frame``(``)`

		数量
年份	性别
---	---	---
1880	F	83929
M	110490
1881	F	85034
...	...	...
2019	M	1785527
2020	F	1581301
M	1706423

282 rows × 1 columns

这将生成一个包含计数的 pd.Series。我们还可以想象相同的数据，将 Sex 索引级别“透视”到 dataframe 的列中。通过一个例子更容易理解：

`mf_pivot` `=` `pd``.``pivot_table``(`
    `baby``,`
    `index``=``'``Year``'``,`   `# Column to turn into new index`
    `columns``=``'``Sex``'``,`  `# Column to turn into new columns`
    `values``=``'``Count``'``,` `# Column to aggregate for values`
    `aggfunc``=``sum``)`    `# Aggregation function`
`mf_pivot`

性别	F	M
年份
---	---	---
1880	83929	110490
1881	85034	100738
1882	99699	113686
...	...	...
2018	1676884	1810309
2019	1651911	1785527
2020	1581301	1706423

141 rows × 2 columns

注意

如我们在mf_pivot表中所见，数据框的索引也可以命名。为了阅读输出，重要的是注意数据框有两列，M 和 F，存储在名为 Sex 的索引中。同样，数据框有 141 行，每行有自己的标签。这些标签存储在名为 Year 的索引中。这里，Sex 和 Year 是数据框索引的名称，不是行或列标签本身。

注意，在透视表和使用.groupby()生成的表中，数据值是相同的；只是排列方式不同。透视表可以使用两个属性快速汇总数据，通常出现在文章和论文中。

函数px.line()也能很好地与透视表配合使用，因为该函数在表中的每一列数据上绘制一条线：

`fig` `=` `px``.``line``(``mf_pivot``,` `width``=``350``,` `height``=``250``)`
`fig``.``update_traces``(``selector``=``1``,` `line_dash``=``'``dashdot``'``)`
`fig``.``update_yaxes``(``title``=``'``Value``'``)`

本节介绍了使用pandas中的.groupby()函数以及一个或多个列使用pd.pivot_table()函数对数据进行聚合的常见方法。在下一节中，我们将解释如何将数据框连接在一起。

连接

数据科学家经常希望连接两个或多个数据框，以便跨数据框连接数据值。例如，一个在线书店可能有一个数据框，其中包含每位用户订购的书籍，以及一个包含每本书流派的第二个数据框。通过将这两个数据框连接起来，数据科学家可以看到每位用户偏爱哪些流派。

我们将继续查看婴儿姓名数据。我们将使用连接来检查纽约时报关于婴儿姓名的文章中提到的一些趋势。文章讨论了某些类别的姓名如何随时间变得更受欢迎或不受欢迎。例如，它提到了神话般的姓名如 Julius 和 Cassius 变得流行，而婴儿潮时期的姓名如 Susan 和 Debbie 则变得不那么流行。这些类别的流行度随时间的变化如何？

我们将纽约时报文章中的名称和类别放入了一个小数据框中：

`nyt` `=` `pd``.``read_csv``(``'``nyt_names.csv``'``)`
`nyt`

	nyt_name	category
0	Lucifer	forbidden
1	Lilith	forbidden
2	Danger	forbidden
...	...	...
20	Venus	celestial
21	Celestia	celestial
22	Skye	celestial

23 rows × 2 columns

要查看名称类别的流行程度，我们将nyt数据框与baby数据框连接以从baby获取名称计数：

`baby` `=` `pd``.``read_csv``(``'``babynames.csv``'``)`
`baby`

	Name	Sex	Count	Year
0	Liam	M	19659	2020
1	Noah	M	18252	2020
2	Oliver	M	14147	2020
...	...	...	...	...
2020719	Verona	F	5	1880
2020720	Vertie	F	5	1880
2020721	Wilma	F	5	1880

2020722 rows × 4 columns

对于直觉，我们可以想象在 baby 的每一行中进行以下询问：这个名称是否在 nyt 表中？如果是，那么将 category 列的值添加到该行。这就是连接背后的基本思想。让我们首先看一些较小数据框的示例。

内连接

我们首先制作 baby 和 nyt 表的较小版本，这样在连接表格时更容易看到发生的情况：

`nyt_small`

	nyt_name	category
0	Karen	战斗者
1	Julius	神话
2	Freya	神话

`baby_small`

	名称	性别	计数	年份
0	Noah	M	18252	2020
1	Julius	M	960	2020
2	Karen	M	6	2020
3	Karen	F	325	2020
4	Noah	F	305	2020

要在 pandas 中连接表，我们将使用 .merge() 方法：

`baby_small``.``merge``(``nyt_small``,`
                 `left_on``=``'``Name``'``,`        `# column in left table to match`
                 `right_on``=``'``nyt_name``'``)`   `# column in right table to match`

	名称	性别	计数	年份	nyt_name	category
0	Julius	M	960	2020	Julius	神话
1	Karen	M	6	2020	Karen	boomer
**	2	Karen	F	325	2020	Karen	boomer

注意，新表具有 baby_small 和 nyt_small 表的列。名称为 Noah 的行已消失。其余行从 nyt_small 中获得了匹配的 category。

注意

读者还应注意，pandas 具有 .join() 方法用于将两个数据框连接在一起。然而，.merge() 方法在数据框连接方面更加灵活，因此我们专注于 .merge()。我们鼓励读者查阅 pandas 文档，了解这两者之间的确切区别。

当我们将两个表连接在一起时，我们告诉 pandas 我们要使用哪些列（left_on 和 right_on 参数）来进行连接。当连接列中的值匹配时，pandas 将行进行匹配，如图 6-4 所示。

图 6-4。要进行连接，pandas 使用 Name 和 nyt_name 列中的值进行行匹配，删除没有匹配值的行。

默认情况下，pandas 执行内连接。如果任一表中的行在另一表中没有匹配项，pandas 将从结果中删除这些行。在本例中，baby_small 中的 Noah 行在 nyt_small 中没有匹配项，因此被删除。同样，nyt_small 中的 Freya 行也没有在 baby_small 中找到匹配项，因此也被删除。只有在两个表中都有匹配项的行才会留在最终结果中。

左连接、右连接和外连接

有时我们希望保留没有匹配项的行，而不是完全删除它们。还有其他类型的连接——左连接、右连接和外连接——即使在没有匹配项时也会保留行。

在左连接中，保留左表中没有匹配项的行，如图 6-5 所示。

图 6-5。在左连接中，保留左表中没有匹配值的行

要在pandas中进行左连接，调用.merge()时使用how='left'参数：

`baby_small``.``merge``(``nyt_small``,`
                 `left_on``=``'``Name``'``,`        
                 `right_on``=``'``nyt_name``'``,`
                 `how``=``'``left``'``)`           `# left join instead of inner`

	Name	Sex	Count	Year	nyt_name	category
0	Noah	M	18252	2020	NaN	NaN
1	Julius	M	960	2020	Julius	mythology
2	Karen	M	6	2020	Karen	boomer
3	Karen	F	325	2020	Karen	boomer
4	Noah	F	305	2020	NaN	NaN

注意Noah行在最终表中被保留。由于这些行在nyt_small数据框中没有匹配，连接在nyt_name和category列中留下NaN值。同时注意，nyt_small中的Freya行仍然被丢弃。

右连接与左连接类似，但保留右表中没有匹配的行而不是左表：

`baby_small``.``merge``(``nyt_small``,`
                 `left_on``=``'``Name``'``,`        
                 `right_on``=``'``nyt_name``'``,`
                 `how``=``'``right``'``)`

	Name	Sex	Count	Year	nyt_name	category
0	Karen	M	6.0	2020.0	Karen	boomer
1	Karen	F	325.0	2020.0	Karen	boomer
2	Julius	M	960.0	2020.0	Julius	mythology
3	NaN	NaN	NaN	NaN	Freya	mythology

最后，外连接保留两个表中的行，即使它们没有匹配：

`baby_small``.``merge``(``nyt_small``,`
                 `left_on``=``'``Name``'``,`        
                 `right_on``=``'``nyt_name``'``,`
                 `how``=``'``outer``'``)`

	Name	Sex	Count	Year	nyt_name	category
0	Noah	M	18252.0	2020.0	NaN	NaN
1	Noah	F	305.0	2020.0	NaN	NaN
2	Julius	M	960.0	2020.0	Julius	mythology
3	Karen	M	6.0	2020.0	Karen	boomer
4	Karen	F	325.0	2020.0	Karen	boomer
5	NaN	NaN	NaN	NaN	Freya	mythology

示例：纽约时报名字类别的流行程度

现在让我们回到完整的数据框baby和nyt。.head()用于获取前几行数据，节省空间：

`baby``.``head``(``2``)`

	Name	Sex	Count	Year
0	Liam	M	19659	2020
1	Noah	M	18252	2020

`nyt``.``head``(``2``)`

	nyt_name	category
0	Lucifer	forbidden
1	Lilith	forbidden

我们想要了解nyt中名字类别的流行度随时间的变化。为了回答这个问题：

1. 用baby和nyt进行内连接。

2. 将表按category和Year分组。

3. 使用求和对计数进行聚合：

`cate_counts` `=` `(`
    `baby``.``merge``(``nyt``,` `left_on``=``'``Name``'``,` `right_on``=``'``nyt_name``'``)` `# [1]`
    `.``groupby``(``[``'``category``'``,` `'``Year``'``]``)`                       `# [2]`
    `[``'``Count``'``]`                                            `# [3]`
    `.``sum``(``)`                                               `# [3]`
    `.``reset_index``(``)`
`)`
`cate_counts`

	category	Year	Count
0	boomer	1880	292
1	boomer	1881	298
2	boomer	1882	326
...	...	...	...
647	mythology	2018	2944
648	mythology	2019	3320
649	mythology	2020	3489

650 rows × 3 columns

现在我们可以绘制boomer名字和mythology名字的流行度：

据纽约时报文章称，自 2000 年以来，婴儿潮一代的名字变得不太流行，而神话名字则变得更受欢迎。

我们还可以一次性绘制所有类别的流行度。查看下面的图表，看看它们是否支持纽约时报文章中的观点：

在本节中，我们介绍了数据框的连接操作。当将数据框连接在一起时，我们使用.merge()函数匹配行。在连接数据框时，考虑连接的类型（内部、左侧、右侧或外部）是很重要的。在下一节中，我们将解释如何转换数据框中的值。

转换中

当数据科学家需要以相同的方式更改特征中的每个值时，他们会转换数据框列。例如，如果一个特征包含以英尺表示的人的身高，数据科学家可能希望将身高转换为厘米。在本节中，我们将介绍apply，这是一种使用用户定义函数转换数据列的操作：

`baby` `=` `pd``.``read_csv``(``'``babynames.csv``'``)`
`baby`

	Name	Sex	Count	Year
0	Liam	M	19659	2020
1	Noah	M	18252	2020
2	Oliver	M	14147	2020
...	...	...	...	...
2020719	Verona	F	5	1880
2020720	Vertie	F	5	1880
2020721	Wilma	F	5	1880

2020722 rows × 4 columns

在纽约时报的婴儿名字文章中，帕梅拉提到，以字母 L 或 K 开头的名字在 2000 年后变得流行。另一方面，以字母 J 开头的名字在 1970 年代和 1980 年代达到了流行高峰，但自那以后流行度下降。我们可以使用baby数据集验证这些说法。

我们通过以下步骤解决这个问题：

1. 将Name列转换为一个新列，其中包含每个Name值的第一个字母。

2. 根据第一个字母和年份对数据框进行分组。

3. 通过求和聚合名称计数。

要完成第一步，我们将应用一个函数到Name列。

应用

pd.Series对象包含一个.apply()方法，该方法接受一个函数并将其应用于系列中的每个值。例如，要找出每个名称的长度，我们应用len函数：

`names` `=` `baby``[``'``Name``'``]`
`names``.``apply``(``len``)`

0          4
1          4
2          6
          ..
2020719    6
2020720    6
2020721    5
Name: Name, Length: 2020722, dtype: int64

要提取每个名称的第一个字母，我们定义一个自定义函数，并将其传递给.apply()。函数的参数是系列中的一个单独值：

`def` `first_letter``(``string``)``:`
    `return` `string``[``0``]`

`names``.``apply``(``first_letter``)`

0          L
1          N
2          O
          ..
2020719    V
2020720    V
2020721    W
Name: Name, Length: 2020722, dtype: object

使用.apply()类似于使用for循环。前面的代码大致相当于写成：

`result` `=` `[``]`
`for` `name` `in` `names``:`
    `result``.``append``(``first_letter``(``name``)``)`

现在我们可以将首字母分配给数据框中的新列：

`letters` `=` `baby``.``assign``(``Firsts``=``names``.``apply``(``first_letter``)``)`
`letters`

	Name	Sex	Count	Year	Firsts
0	Liam	M	19659	2020	L
1	Noah	M	18252	2020	N
2	Oliver	M	14147	2020	O
...	...	...	...	...	...
2020719	Verona	F	5	1880	V
2020720	Vertie	F	5	1880	V
2020721	Wilma	F	5	1880	W

2020722 rows × 5 columns

注意

要在数据框中创建一个新列，您可能还会遇到以下语法：

`baby``[``'``Firsts``'``]` `=` `names``.``apply``(``first_letter``)`

这会通过添加一个名为 Firsts 的新列来改变 baby 表。在前面的代码中，我们使用了 .assign()，它不会改变 baby 表本身，而是创建了一个新的数据框。改变数据框并不是错误的，但可能是错误的常见源泉。因此，在本书中，我们大多数情况下会使用 .assign()。

示例： “L” 开头名字的流行程度

现在我们可以使用 letters 数据框来查看首字母随时间变化的流行程度：

`letter_counts` `=` `(``letters`
 `.``groupby``(``[``'``Firsts``'``,` `'``Year``'``]``)`
 `[``'``Count``'``]`
 `.``sum``(``)`
 `.``reset_index``(``)`
`)`
`letter_counts`

	Firsts	Year	Count
0	A	1880	16740
1	A	1881	16257
2	A	1882	18790
...	...	...	...
3638	Z	2018	55996
3639	Z	2019	55293
3640	Z	2020	54011

3641 rows × 3 columns

`fig` `=` `px``.``line``(``letter_counts``.``loc``[``letter_counts``[``'``Firsts``'``]` `==` `'``L``'``]``,`
              `x``=``'``Year``'``,` `y``=``'``Count``'``,` `title``=``'``Popularity of` `"``L``"` `names``'``,`
              `width``=``350``,` `height``=``250``)`
`fig``.``update_layout``(``margin``=``dict``(``t``=``30``)``)`

这张图表显示，在 1960 年代，“L”开头的名字很流行，在此后的几十年里有所下降，但自 2000 年以来确实重新流行起来。

那么，“J” 开头的名字怎么样？

`fig` `=` `px``.``line``(``letter_counts``.``loc``[``letter_counts``[``'``Firsts``'``]` `==` `'``J``'``]``,`
              `x``=``'``Year``'``,` `y``=``'``Count``'``,` `title``=``'``Popularity of` `"``J``"` `names``'``,`
              `width``=``350``,` `height``=``250``)`
`fig``.``update_layout``(``margin``=``dict``(``t``=``30``)``)`

纽约时报 的文章指出，“J”开头的名字在 1970 年代和 80 年代很流行。图表也证实了这一点，并显示自 2000 年以来它们变得不那么流行。

.apply() 的代价

.apply() 的强大之处在于其灵活性——你可以用任何接受单个数据值并输出单个数据值的函数来调用它。

然而，它的灵活性也有一个代价。使用 .apply() 可能会很慢，因为 pandas 不能优化任意函数。例如，对于数值计算，使用 .apply() 比直接在 pd.Series 对象上使用向量化操作要慢得多：

`%``%``timeit`

`# Calculate the decade using vectorized operators`
`baby``[``'``Year``'``]` `/``/` `10` `*` `10`

9.66 ms ± 755 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

`%``%``timeit`

`def` `decade``(``yr``)``:`
    `return` `yr` `/``/` `10` `*` `10`

`# Calculate the decade using apply`
`baby``[``'``Year``'``]``.``apply``(``decade``)`

658 ms ± 49.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

使用 .apply() 的版本要慢 30 倍！特别是对于数值计算，我们建议直接在 pd.Series 对象上进行向量化操作。

在本节中，我们介绍了数据转换。为了在数据框中转换值，我们通常使用 .apply() 和 .assign() 函数。在下一节中，我们将比较数据框与其他表示和操作数据表的方法。

数据框与其他数据表示方式有何不同？

数据框只是表示表中存储数据的一种方式。在实践中，数据科学家会遇到许多其他类型的数据表，如电子表格、矩阵和关系表。在本节中，我们将比较数据框与其他表示方式，解释为什么数据框在数据分析中如此广泛使用。我们还将指出其他表示方式可能更合适的场景。

数据框和电子表格

电子表格是计算机应用程序，用户可以在网格中输入数据并使用公式进行计算。今天一个众所周知的例子是微软 Excel，尽管电子表格可以追溯到至少 1979 年的VisiCalc。电子表格使得直接查看和操作数据变得很容易，因为电子表格公式可以在数据更改时自动重新计算结果。相比之下，数据框代码通常需要在数据集更新时手动重新运行。这些特性使得电子表格非常受欢迎——根据2005 年的估计，有超过 5500 万的电子表格用户，而行业中只有 300 万专业程序员。

数据框比电子表格有几个关键优势。在类似 Jupyter 的计算笔记本中编写数据框代码会自然地产生数据谱系。打开笔记本的人可以看到笔记本的输入文件以及数据是如何更改的。电子表格不显示数据谱系；如果一个人手动编辑单元格中的数据值，将很难让未来的用户看到哪些值是手动编辑的以及如何编辑的。数据框可以处理比电子表格更大的数据集，用户还可以使用分布式编程工具来处理很难加载到电子表格中的大数据集。

数据框和矩阵

矩阵是用于线性代数操作的二维数据数组。在下一个例子中，X是一个具有三行两列的矩阵：

X = [ 1 0 0 4 0 0 ]

矩阵是由它们允许的运算符定义的数学对象。例如，矩阵可以相加或相乘。矩阵也有转置。这些运算符具有数据科学家在统计建模中依赖的非常有用的特性。

矩阵和数据框之间的一个重要区别是，当矩阵被视为数学对象时，它们只能包含数字。另一方面，数据框还可以包含文本等其他类型的数据。这使得数据框更适合加载和处理可能包含各种数据类型的原始数据。实际上，数据科学家经常将数据加载到数据框中，然后将数据处理成矩阵形式。在本书中，我们通常使用数据框进行探索性数据分析和数据清洗，然后将数据处理成矩阵形式用于机器学习模型。

注意

数据科学家将矩阵称为数学对象，也称为程序对象。例如，R 编程语言有一个矩阵对象，而在 Python 中，我们可以使用二维的numpy数组表示矩阵。在 Python 和 R 中实现的矩阵可以包含除了数字以外的其他数据类型，但这样做会丧失数学属性。这是领域可以用同一术语指称不同事物的另一个例子。

数据框和关系

关系是数据库系统中使用的数据表表示形式，特别是像 SQLite 和 PostgreSQL 这样的 SQL 系统。（我们在第七章中介绍了关系和 SQL。）关系与数据框架有许多相似之处；它们都使用行来表示记录，列来表示特征。两者都有列名，并且列内的数据具有相同的类型。

数据框架的一个关键优势在于，它们不需要行来表示记录，也不需要列来表示特征。许多时候，原始数据并不以直接放入关系中的方便格式出现。在这些情况下，数据科学家使用数据框架来加载和处理数据，因为数据框架在这方面更加灵活。通常，数据科学家会将原始数据加载到数据框架中，然后处理数据，使其能够轻松地存储在关系中。

关系型数据库系统（如PostgreSQL）比数据框架具有的一个关键优势在于，它们具有非常有用的数据存储和管理功能。考虑一家运营大型社交媒体网站的数据科学家。数据库可能包含的数据量太大，无法一次性读入pandas数据框架；因此，数据科学家使用 SQL 查询来子集化和聚合数据，因为数据库系统更能处理大型数据集。此外，网站用户通过发布帖子、上传图片和编辑个人资料不断更新其数据。在这种情况下，数据库系统让数据科学家能够重复使用现有的 SQL 查询更新他们的分析，而不是反复下载大型 CSV 文件。

摘要

在本章中，我们解释了数据框架是什么，它们为什么有用，以及如何使用pandas代码与它们一起工作。子集化、聚合、连接和转换几乎在每个数据分析中都是有用的。在本书的其余部分中，特别是第 8、第 9 和第十章中，我们将经常依赖这些操作。