Python贝叶斯分析（第2版）

版权信息

书名：Python贝叶斯分析（第2版）

ISBN：978-7-115-60089-9

本书由人民邮电出版社发行数字版。版权所有，侵权必究。

您购买的人民邮电出版社电子书仅供您个人使用，未经授权，不得以任何方式复制和传播本书内容。

我们愿意相信读者具有这样的良知和觉悟，与我们共同保护知识产权。

如果购买者有侵权行为，我们可能对该用户实施包括但不限于关闭该帐号等维权措施，并可能追究法律责任。

版　　权

著    [阿根廷] 奥斯瓦尔多·马丁（Osvaldo Martin）

译    张天旭　黄雪菊

责任编辑 胡俊英

人民邮电出版社出版发行　　北京市丰台区成寿寺路11号

邮编　100164 　电子邮件　315@ptpress.com.cn

网址　http://www.ptpress.com.cn

读者服务热线：(010)81055410

反盗版热线：(010)81055315

读者服务：

微信扫码关注【异步社区】微信公众号，回复“e60089”获取本书配套资源以及异步社区15天VIP会员卡，近千本电子书免费畅读。

内 容 提 要

本书是一本概率编程的入门书。本书使用概率编程库PyMC3以及可视化库ArviZ对贝叶斯统计分析的相关知识进行讲解，包括概率思维、概率编程、线性回归建模、广义线性模型、模型比较、混合模型、高斯过程以及推断引擎等知识。全书图文并茂，通俗易懂，适合具备一定Python基础的读者学习使用。学完本书，读者可以利用概率思维建立贝叶斯模型并解决自己的数据分析问题。

序　　言

概率编程是一类计算机编程框架，用于灵活地构建贝叶斯模型，一旦构建好模型，强大的推断算法便可以独立于特定的模型而工作，并通过模型拟合数据。把灵活的模型定义与自动推断结合在一起便得到了一个强大的工具，方便研究者们快速地构建、分析和迭代新的统计模型。这个迭代过程与以往用贝叶斯模型拟合数据的方式有很大的不同：以往的推断算法仅对某一特定的模型有效。这不仅导致人们在构建模型和设计推断算法的时候，需要具备很强的数据技巧，还降低了迭代速度：需要先修改模型，然后重新设计推断算法。而概率编程则将统计建模的过程大众化，大大降低了对使用者的数学水平的要求，缩短了构建新模型时所需花费的时间，同时还能增强使用者对数据内涵的洞察。

概率编程背后的思想并不新鲜。BUGS是最早的概率编程实践之一，于1989年首次发布。由于能够成功应用的模型非常有限，而且推断过程很慢，因此这些第一代语言不太实用。如今，人们已经开发出了许多概率编程语言，并在学术界和各大公司（如谷歌、微软、亚马逊）广泛地用于解决各种复杂问题。那现代的概率编程语言有哪些变化呢？最大的变化来自哈密顿蒙特卡洛采样算法，它相较以往的采样算法要高出若干个数量级，以往的算法只能用来解决一些玩具问题，而如今的算法可以用于解决非常复杂的大规模问题。尽管这些采样器起源于1987年，但得益于最近的一些概率编程系统，如Stan和PyMC3，它们才被广泛地使用起来。

本书将从务实的角度介绍概率编程这一强大而灵活的工具，它将影响你如何思考和解决复杂的分析问题。作为PyMC3的核心开发者之一，没有人比Osvaldo Martin更适合来写这本书了。Osvaldo非常擅长将复杂的问题拆解成容易理解和吸收的部分，他宝贵的实战经验将带领读者穿过这片复杂的领域。书中的图表和代码是非常有用的资源，这些都将增进读者对背后理论知识的直观理解。

此外，我还要称赞此刻拿起本书的亲爱的读者。如今这个时代，各大新闻头条都在鼓吹深度学习才是解决所有当前和未来分析问题的技术，为某个具体目的而特地构建一个定制化的模型似乎显得不那么吸引人。不过，通过概率编程，你将能够解决其他方法很难解决的复杂问题。

这并不是说深度学习不是那么令人兴奋的技术。事实上，概率编程并不局限于经典的统计模型。阅读最近一些机器学习文献，你就会发现，贝叶斯统计正作为一个强大的框架被用来表示和理解下一代深度神经网络。本书不仅会让你具备分析和解决复杂问题的技能，还会带你走进可能是人类智慧的最前沿：开发人工智能。尽情享受吧！

Tomas Wiecki 博士

Quantopian首席研究员

前　　言

贝叶斯统计距今已经有超过250年的历史，其间该方法既饱受赞誉又备受轻视。直到近几十年，得益于理论进步和计算能力的提升，贝叶斯统计才越来越多地受到来自统计学以及其他学科乃至学术圈以外工业界的重视。现代的贝叶斯统计主要是计算统计学，人们对模型的灵活性、透明性以及统计分析结果的可解释性的追求最终造就了该趋势。

本书将从实用的角度来介绍贝叶斯统计，不会太在意其他统计范式及其与贝叶斯统计的关系。本书的目的是介绍如何做贝叶斯数据分析，尽管与之相关的哲学讨论也很有趣，不过受限于篇幅，这部分内容并不在本书的讨论范围之内，有兴趣的读者可以通过其他方式深入了解。

这里我们采用建模的方式介绍统计学，介绍如何从概率模型的角度思考问题，并应用贝叶斯定理来推导模型和数据的逻辑结果。这种建模方式是可编程的，其中，模型部分会用PyMC3和ArviZ编写。PyMC3是用于贝叶斯统计的Python库，它为用户封装了大量的数学细节和计算过程。ArviZ是一个用于贝叶斯模型探索式分析的Python库。

贝叶斯方法在理论上以概率论为基础，很多介绍贝叶斯方法的书都充斥着复杂的数学公式。学习统计学方面的数学知识显然有利于构建更好的模型，同时还能让你对问题、模型和结果有更好的把握。不过类似PyMC3的库能够帮助你在有限的数学知识水平下学习并掌握贝叶斯统计。在阅读本书的过程中，你将亲自见证这一过程。

读者对象

如果你是一名学生、数据科学家、自然科学或者社会科学的研究人员或者开发人员，想要着手贝叶斯数据分析和概率编程，那么本书就是为你准备的。本书是入门级的，因此并不需要你提前掌握统计知识，当然，有一点儿Python和NumPy的使用经验就更好了。

本书结构

第1章，概率思维。本章介绍贝叶斯统计的基本概念及其在数据分析中的意义。本章包含本书其余章节中用到的基本思想。

第2章，概率编程。本章从编程的角度重新回顾第1章提到的概念，并介绍PyMC3和ArviZ。PyMC3是一个概率编程库，ArviZ是一个用于贝叶斯模型探索性分析的Python库。此外，本章还将用几个例子解释分层模型。

第3章，线性回归建模。线性回归是构建更复杂模型的基石，有着广泛的应用，本章将涵盖线性模型的基础部分。

第4章，广义线性模型。本章介绍如何将线性模型扩展到高斯分布之外的其他分布，为解决许多数据分析问题“打开大门”。

第5章，模型比较。本章讨论如何使用WAIC、LOO和贝叶斯因子来比较、选择和平均模型，同时讨论这些方法的一般注意事项。

第6章，混合模型。本章讨论如何通过混合简单的分布来构建更复杂的分布以增加模型的灵活性。本章还会介绍本书的第一个非参模型——狄利克雷过程。

第7章，高斯过程。本章涵盖高斯过程背后的基本思想，以及如何使用它们在函数上构建非参数模型，从而解决一系列的问题。

第8章，推断引擎。本章介绍数值逼近后验分布的方法，同时讨论一个从使用者的角度来说非常重要的主题——如何诊断后验估计的可靠度。

第9章，拓展学习。最后一章列出一系列继续深入学习的资料以及一段简短的演讲。

准备工作

本书的代码采用Python 3.6编写，要安装Python及各种Python库，我建议用Anaconda，这是一个用于科学计算的软件，可以从Anaconda官网下载。Anaconda会在你的系统上安装很多非常有用的Python库。此外，你还需要安装两个库。首先，用conda命令安装PyMC3：

conda install -c conda-forge pymc3

然后用下面的命令安装ArviZ：

pip install arviz

另一种方式是你可以安装本书所有必需的Python库：安装Anaconda之后可以从本书配套源码文件中找到一个bap.yml文件，然后用下面的命令安装所有必需的库：

conda env create -f bap.yml

本书用到的Python库列举如下。

● IPython 7.0。

● Jupyter 1.0 （或Jupyter-lab 0.35）。

● NumPy 1.14.2。

● SciPy 1.1。

● pandas 0.23.4。

● Matplotlib 3.0.2。

● Seaborn 0.9.0。

● ArviZ 0.3.1。

● PyMC3 3.6。

对于每章的代码部分，如果你安装了上面的这些库，那么与其从本书复制并粘贴代码运行，不如直接从异步社区网站的本书页面下载代码，然后用Jupyter Notebook或者Jupyter Lab运行。我会为PyMC3或ArviZ的新版本更新GitHub上面的代码。如果你在运行本书相关代码的过程中发现任何技术问题、语法错误或者任何其他错误，可以直接在上面的链接里提问，我会尽快解答。

本书的大多数图或表格都是用代码生成的，你会发现，一般一段代码后面都会跟着这段代码生成的图或表格，希望这种方式能够让熟悉Jupyter Notebook或者Jupyter Lab的读者感到亲切，但愿这样做不会让读者感到困惑。

格式约定

本书遵照一些文本编排的惯例。

代码中的单词、文件名或函数名在正文段落中仍以代码体呈现。

代码片段的格式一般如下：

μ = 0. 
σ = 1.
X = stats.norm(μ, σ)
x = X.rvs(3)

粗体：标明一个新的词汇，或者重要的词。

作者简介

奥斯瓦尔多·马丁（Osvaldo Martin）是阿根廷国家科学与技术研究理事会（CONICET）的一名研究员。他曾从事蛋白质、多糖及RNA分子等结构生物信息学方面的研究，此外，在应用马尔可夫链蒙特卡洛方法模拟分子动力学方向上有着丰富的经验，他喜欢用Python解决数据分析中的问题。

他曾讲授结构生物信息学、数据科学以及贝叶斯数据分析相关的课程，在2017年带头组建了阿根廷圣路易斯PyData委员会。同时，他也是PyMC3以及ArviZ两个项目的核心开发者之一。

英文版审校者简介

Eric J.Ma是诺华生物医学研究所的数据科学家，他负责生物医药方面的数据科学研究，关注使用贝叶斯方法为患者制药。他于2017年春天获得博士学位，并在2017年夏天成为Insight Health Data的研究员。

Eric J.Ma还是一名开源软件开发者，他曾主导开发了NetworkX的可视化包nxviz，以及简化Python清理数据的pyjanitor。此外他还参与并贡献了一系列开源工具，包括PyMC3、Matplotlib、brokeh以及CuPy。

Austin Rochford是Monetate Labs的首席数据科学家，他开发的产品用于帮助众多零售商进行个性化销售。他是一位训练有素的数学家，并且是贝叶斯方法的积极倡导者。

致　　谢

作者的话

我要感谢Romina长久以来的支持，还要感谢Walter Lapadula、Bill Engels、Eric J.Ma以及Austin Rochford给本书提出的宝贵意见和建议。此外，还要感谢PyMC3和ArviZ社区的核心开发者和贡献者，正是因为他们对这些库以及整个社区的热爱和付出才有了本书的问世。

译者的话

本书的第1版由田俊翻译，在此对田俊老师的工作表示感谢！第2版的翻译和校对由我在业余时间进行，断断续续持续了很长时间，期间烦琐的校对大都有赖于我的妻子黄雪菊的帮助才能完成，在此深表感谢！

在中译本初稿完成之后，我和妻子又对全书进行了多次勘核：对书中出现的英文术语多番查阅资料，挑选合适的中文用词；对书中生涩的隐喻、典故等，查询其出处，并加上了必要的注释。另外，我们还修正了原书中的一些错误。希望此番努力能减少读者朋友的阅读障碍，能对读者的实际工作和学习有所帮助。当然，由于我们水平有限，疏漏之处在所难免，希望读者批评指正。欢迎大家通过电子邮件（zhang_tianxu@sina.com）与我交流探讨。

服务与支持

本书由异步社区出品，社区（https://www.epubit.com）为您提供后续服务。

您还可以扫码右侧二维码， 关注【异步社区】微信公众号，回复“e60089”直接获取,同时可以获得异步社区15天VIP会员卡，近千本电子书免费畅读。

配套资源

本书提供配套资源，请在异步社区本书页面中点击，跳转到下载界面，按提示进行操作即可。

提交勘误信息息信息

作者、译者和编辑尽最大努力来确保书中内容的准确性，但难免会存在疏漏。欢迎您将发现的问题反馈给我们，帮助我们提升图书的质量。

当您发现错误时，请登录异步社区，按书名搜索，进入本书页面，单击“发表勘误”，输入错误信息，单击“提交勘误”按钮即可，如下图所示。本书的作者和编辑会对您提交的错误信息进行审核，确认并接受后，您将获赠异步社区的100积分。积分可用于在异步社区兑换优惠券、样书或奖品。

与我们联系

我们的联系邮箱是contact@epubit.com.cn。

如果您对本书有任何疑问或建议，请您发邮件给我们，并请在邮件标题中注明本书书名，以便我们更高效地做出反馈。

如果您有兴趣出版图书、录制教学视频，或者参与图书翻译、技术审校等工作，可以发邮件给我们；有意出版图书的作者也可以到异步社区投稿（直接访问www.epubit.com/contribute即可）。

如果您所在的学校、培训机构或企业想批量购买本书或异步社区出版的其他图书，也可以发邮件给我们。

如果您在网上发现有针对异步社区出品图书的各种形式的盗版行为，包括对图书全部或部分内容的非授权传播，请您将怀疑有侵权行为的链接通过邮件发送给我们。您的这一举动是对作者权益的保护，也是我们持续为您提供有价值的内容的动力之源。

关于异步社区和异步图书

“异步社区”是人民邮电出版社旗下IT专业图书社区，致力于出版精品IT图书和相关学习产品，为作译者提供优质出版服务。异步社区创办于2015年8月，提供大量精品IT图书和电子书，以及高品质技术文章和视频课程。更多详情请访问异步社区官网https://www.epubit.com。

“异步图书”是由异步社区编辑团队策划出版的精品IT专业图书的品牌，依托于人民邮电出版社的计算机图书出版积累和专业编辑团队，相关图书在封面上印有异步图书的LOGO。异步图书的出版领域包括软件开发、大数据、人工智能、测试、前端、网络技术等。

异步社区

微信服务号

第1章　概率思维

“归根结底，概率论只不过是把常识简化为计算。”

——西蒙·拉普拉斯（Simon Laplace）^[1]

[1]　西蒙·拉普拉斯是天体力学的主要奠基人，是天体演化学的创立者之一，是分析概率论的创始人，是应用数学的先驱。——译者注

本章将介绍贝叶斯统计中的核心概念，以及一些用于贝叶斯分析的基本工具。大部分是一些理论介绍，其中也会涉及一些 Python代码。本章中提及的大多数概念会在本书后文中反复提到。本章内容有点儿偏理论，对习惯代码的你来说可能会感到焦虑，不过学习这些理论知识，会让你在后面应用贝叶斯统计方法解决问题时更容易一些。本章包含以下主题。

● 统计建模。

● 概率与不确定性。

● 贝叶斯定理及统计推断。

● 单参数推断以及经典的抛硬币问题。

● 先验选择。

● 报告贝叶斯分析结果。

1.1　统计学、模型以及本书采用的方法

统计学主要是关于收集、组织、分析并解释数据的科学，统计学的基础知识对数据分析来说至关重要。在数据分析中，主要有以下两种统计学方法。

❏ 探索性数据分析（Exploratory Data Analysis，EDA）：数值统计，比如均值、众数、标准差以及四分位距等，这部分内容也称作描述性统计。此外EDA还涉及用一些你可能已经熟悉的工具（如直方图或散点图）对数据做可视化分析。

❏ 统计推断：主要是指在已有数据基础上做陈述。我们可能希望了解一些特定的现象，也可能是想对未来（或尚未观测到）的数据进行预测，又或者是希望从对观测值的多个解释中找出最合理的一个。统计推断为解决这类问题提供了一系列方法和工具。

提示：本书重点关注如何做贝叶斯统计推断，然后用EDA对贝叶斯推断的结果做总结、解释、检查和交流。

大多数统计学入门课程，至少对非统计学专业的人而言就像一份菜谱，这些菜谱或多或少是这样的：首先，到统计学的后厨拿一瓶罐头并打开，加点儿数据上去尝尝，然后不停搅拌直到得出一个稳定的值，该值最好低于0.05。这类课程的目的是教会你如何选择一瓶合适的罐头。我从来不喜欢这种方法，主要是因为最常见的结果是人们会很困惑，甚至连概念都无法掌握。本书采用的是另外一种方式：首先我们也需要点儿原料，不过这次是自己亲自做的而不是买来的罐头，然后学习如何把新鲜的食材混合在一起以适应不同的烹饪场景，更重要的是教会你如何把这些概念应用到本书例子之外的地方。

采用这种方式有两方面原因。

● 本体论：统计学是建立在概率论数学框架之下的一种统一的建模方式。概率论的方法能为一些看起来非常不一样的方法提供统一的视角，比如统计学和机器学习（Machine Learning，ML）从概率论的角度来看就非常相似。

● 技术：如PyMC3这样的现代软件允许实践者以相对简单的方式定义和解决问题。在几年前，这类问题可能是无法解决或者需要很高的数学水平和技术复杂度。

1.1.1　与数据打交道

数据是统计学十分基本的组成部分。数据有多种来源，比如实验、计算机模拟、调查以及实地观测等。假如我们是负责数据生成或收集的人，首先考虑的是要解决什么问题以及打算采用什么方法，然后再着手准备数据。事实上，统计学有一个叫作实验设计的分支，专门研究如何获取数据。在这个数据泛滥的年代，我们有时候会忘了获取数据并非总是很方便的。比如，一个大型强子对撞机（Large Hadron Collider，LHC）一天能产生上百TB的数据，但建造这个装置却要花费数年的人力和智力。

通常，可以认为生成数据的过程是随机的，这可能是事物本身、技术性因素又或者是认知的不确定性导致的。也就是说，系统本身具有不确定性，一些技术性问题会增加噪声或限制我们无法以任意精度观测数据。此外还有一些概念层面的理解局限导致我们难以揭示系统的细节。以上这些原因，使得我们需要在模型的背景之下来解释数据，包括心理模型和形式模型。数据不会说话，但能通过建模来表达。

本书假设我们已经收集到了数据，并且这些数据都是干净、整洁的（通常这在现实世界中很少见），这个假设能让我们把注意力放到本书的主题上来。我想强调的是，尽管本书并没有涵盖数据清洗这部分内容，但想要成功地与数据打交道，这些是你应该学习和实践的重要技能。

在数据分析中，掌握一门编程语言（比如Python）是非常实用的。考虑到我们生活在一个复杂的世界中，数据也是杂乱无章的，操作数据通常是必要的，而编程有助于我们完成这类任务。就算你的数据非常干净、整洁，编程仍然非常有用，因为现代贝叶斯统计主要是通过Python或R等编程语言完成的。

如果你想学习如何用Python清洗和操作数据，我推荐Jake VanderPlas写的Python Data Science Handbook一书。

1.1.2　贝叶斯建模

模型是对给定系统或过程的一种简化描述。这些描述只关注系统中某些重要的部分，因此，大多数模型的目的并不是解释整个系统。这也是为什么更复杂的模型并非总是更好的模型。

模型分为很多种，本书主要关注贝叶斯模型。贝叶斯建模过程可以总结为以下3步。

（1）给定一些数据以及这些数据是如何生成的假设，然后通过组合一些概率分布来设计模型。大多数情况下，这些模型是粗略的近似，不过正是我们所需要的。

（2）根据贝叶斯定理将数据添加到模型里，然后把数据和假设结合起来推导出逻辑结果，这就是根据数据调整模型。

（3）检查模型是否有意义可以依据不同的标准，包括数据、我们在这方面的专业知识，有时还通过比较几个模型来评价模型。

通常，实际的建模过程并非是严格按照这3个步骤的顺序进行的。我们会在一些特定点重复这些步骤：可能是犯了一个愚蠢的编程错误，也可能是找到了某种改进模型的方法，又或者是需要增加更多的数据或收集不同的数据集。

贝叶斯模型是基于概率构建的，因此也称作概率模型。为什么基于概率呢？因为概率这个数学工具能够很好地模拟不确定性，接下来让我们进一步了解概率这个工具。

1.2　概率论

本节的标题似乎有些“自命不凡”，不过这里我并不打算通过短短几页就把概率论讲清楚，而且这也不是我的目的。本节主要介绍概率论中一些普通而重要的概念，这些概念能让读者更好地理解贝叶斯方法，同时也为阅读本书后面的内容打下基础。如果有必要，我们会根据需要再展开或介绍一些与概率论相关的新概念。如果你想深入学习概率论，可以参考阅读Joseph K. Blitzstein和Jessica Hwang写的Introduction to Probability。另外一本非常有用的书是Sumio Watanabe写的Mathematical Theory of Bayesian Statistics，该书比第一本更偏向贝叶斯统计，也更侧重数学。

1.2.1　解释概率

尽管概率论是一个相当成熟和完善的数学分支，但关于概率的诠释不止一种。从贝叶斯的角度看，概率是衡量某一命题不确定性水平的方法。按照这种定义，提出以下问题都是自然且合理的：火星上有生命的概率，电子质量大约为千克，或者布宜诺斯艾利斯在1816年7月9日是晴天的概率。值得注意的是，类似于火星上是否有生命这种问题的答案是二值化的：要么有，要么没有。不过鉴于我们并不知道真实情况，更明智的做法是试图找出火星上存在生命的可能性有多大。该命题取决于我们当前所掌握的信息，而非客观的自然属性。由于这种概率的定义与我们大脑的认知有关，因此人们也常称之为概率的主观定义。不过需要注意的是，具备科学素养的人并不会用其个人信念来回答这类问题。相反，他们会用所有与火星相关的地理数据，以及相应条件下是否适合生命生存的生物知识等来回答这类问题。因此，贝叶斯概率，或扩展到贝叶斯统计，与我们已有的其他成熟的科学方法一样主观（或者客观）。

对于一个问题，如果我们没有相关信息，那么可以说每个可能的事件都有同等的可能性，从形式上讲，这相当于给每个可能的事件分配相同的概率。在没有任何信息的时候，事件的不确定性是最大的。假如某些事件的可能性更高，那么就可以给这些事件赋予更高的概率，相应的其他事件的概率就低一些。请注意，当我们在统计学中谈论事件时，并不局限于可能发生的事情，比如，小行星撞击地球或者我姨妈的60岁生日宴。事件只是变量可以采用的任何可能值（或值的子集），比如，你的年龄大于30岁，或者某品牌巧克力蛋糕的价钱，又或者是2017年全世界卖出的自行车的数量。

概率的概念也常常与逻辑相关。在亚里士多德学派，或者传统的逻辑学派看来，一个命题只能是真或者是假。而在贝叶斯学派所定义的概率中，这些不过是一个特例，一个真的命题所对应的概率是1，而一个假的命题对应的概率是0。当有足够数据表明，火星上存在能够生长、繁殖以及其他符合生命体征描述的生物时，我们才将有火星生命这一命题的概率赋值为1。另外，将一个命题的概率置为0要更为困难一些，因为我们可以猜测：或许某些火星人的足迹还未被勘测到，或者是我们所做的实验有些错误，又或者是某些其他原因导致我们错误地认为火星上不存在生命而实际上是存在的。与此相关的是克伦威尔准则（Cromwell’s Rule），其含义是在对逻辑上正确或错误的命题赋予概率值时，应当避免使用0或者1。有意思的是，Richard Cox在数学上证明了如果想在逻辑推理中加入不确定性，就必须使用概率论的知识。我们很快就会看到，贝叶斯定理只是概率规则的逻辑结果。从另一个角度来看，贝叶斯统计是对逻辑学处理不确定性问题的一种延伸，当然这里没有任何对主观推理的轻蔑。

总的来说，用概率来对不确定性建模，并不一定与自然界在本质上是确定性还是随机性的争论有关，也不一定与主观的自我认知有关。相反，这只是一种对不确定性建模的方法而已。我们认识到，大多数现象很难理解，因为我们面对的是不完整而且充满着干扰的数据，而且还受到由进化塑造的灵长类大脑的限制，以及一些其他原因。因此，我们要采用一种明确考虑了不确定性的建模方法。

提示：从实用的角度来看，本节主要想说明，贝叶斯学派使用概率作为量化不确定性的工具。

我们已经讨论了概率的贝叶斯解释，接下来让我们学习概率的一些数学性质。

1.2.2　定义概率

概率值介于0～１（包括0和1），其计算遵循一些法则，其中之一是乘法法则：

　　（1.1）

上述表达式中，A和B同时发生的概率值等于在B发生的条件下A也发生的概率值乘B发生的概率值，其中表示A和B的联合概率，表示条件概率，二者的现实意义是不同的。例如，路面是湿滑的概率跟下雨时路面湿滑的概率是不同的。条件概率可以大于、小于或等于无条件概率。如果B并不能提供任何关于A的信息，那么，也就是说，只有当A和B是相互独立的时候，这个等式才成立。如果事件B能够给出关于事件A的一些信息，那么根据事件B提供的信息不同，事件A可能发生的概率会变得更高或是更低。让我们来看一个简单的例子，掷一个六面骰子的时候，3朝上的概率（）是多少呢？答案是1/6，因为每个数字朝上的概率都是相同的。假设已知掷出的骰子是个奇数，那么它是3的概率（）是多少呢？答案是1/3，因为我们已经知道它是奇数，那么只有可能是中的某个数，而且每种的概率相同。最后，假设已知掷出的骰子数是个偶数，那么它是3的概率（）是多少呢？答案是0，因为我们已经知道这个数是偶数，那么只可能是中的数，因此3不可能出现。

从上面简单的例子可以看出，通过对观测数据进行调节，我们有效地改变了事件原有的概率，也改变了事件的不确定性。条件概率是统计学的核心，不论你要解决的问题是掷骰子还是自动驾驶。

概率分布是数学中的一个概念，用来描述不同事件发生的可能性，通常这些事件限定在一个集合内，比如{1,2,3,4,5,6}，代表了所有可能发生的事件。在统计学里可以这么理解：数据是从某种参数未知的概率分布中生成的。推断，就是根据从这些真实分布中得到的样本（也称作数据集）找出那些参数。通常我们没法直接获取真实的概率分布，只能退而求其次，设计出一个模型来逼近真实的分布。概率模型就是通过恰当地组合概率分布得到的。

提示：注意，通常我们并不知道模型是否正确，因此需要对模型做评估以便获得信心并说服他人我们的模型是适合所要探索或解决的问题的。

如果一个变量X可以用一个概率分布来描述，那么我们把X称为一个随机变量。通常，人们用大写的字母来表示随机变量（比如X），用小写的字母来表示随机变量的一个实例（比如x）。x可以是一个向量，因此可以包含很多独立的值，比如。下面用Python来看一个例子，我们的真实分布是一个正态分布（或者也叫高斯分布），均值为，标准差为。这两个参数准确无误地刻画了一个正态分布。使用SciPy，可以通过stats.norm(μ,σ)定义一个符合正态分布的随机变量X，然后可以通过rvs方法得到一个实例x，下面的例子中，得到了3个值。

μ = 0.
σ = 1.
X = stats.norm(μ, σ)
x = X.rvs(3)

你会发现，每次执行上面的代码，你都会得到不同的随机结果。一旦一个分布的参数确定了，那么x的取值概率也就确定了，唯一随机的是每次实验所得到的x是随机的。关于随机性的定义有一个常见的误解，即可以从一个随机变量中得到任意可能的值，或者所有的值都是等概率的。实际上，一个随机变量可以取到的值以及对应的概率是严格受到概率分布控制的，而随机性只是因为我们无法准确预测每一次试验的结果值。每次执行上面的代码，都会得到3个不同的值，但如果重复执行上千次之后，根据经验可以检查样本的均值会在0附近，并且95%的样本会落在[−1.96,+1.96]内。这里请不要这么果断地相信我的结论，用你的Python技巧实际验证一下。如果你学过与正态分布有关的数学知识，这里就会获得一致的结论。

统计学中，用来描述一个变量服从参数为和的正态分布的写法如下：

　　（1.2）

这里的波浪线～表示服从于某种分布。

提示：在大多数书中，正态分布用方差而不是标准差来表示，因此人们会写成。本书将使用标准差，首先是因为这样更容易解释，其次是因为PyMC3中也是这样表示的。

本书会涉及多个概率分布，每次介绍一个新的分布时，都会先花点时间介绍它。我们从正态分布开始介绍是因为它几乎是概率分布之源。一个变量如果服从正态分布，那么它的概率值可以用下面的表达式来描述：

　　（1.3）

这就是正态分布的概率密度函数，你不用记住这个表达式，给你看这个表达式只是想让你知道这些数字是从哪里来的。前面已经提到过，和是正态分布的两个参数，这两个参数的值一旦确定就完全定义了一个正态分布，可以看到，式（1.3）中的其他值都是常量。可以是任意实数，即，其决定了正态分布的均值（同时还有中位数和众数，它们都相等）。是标准差，其取值只能为正，描述了正态分布的分散程度，取值越大，分布越分散。由于和有无限种组合，因此高斯分布的实例也就有无限多种，所有这些分布属于同一个高斯分布族。

虽然数学公式这种表达形式简洁明了，甚至有人说它很美，不过得承认第一眼看到公式的时候还是不够直观，尤其是数学不太好的人。我们可以尝试用Python代码将公式的含义重新表示出来。先来看看高斯分布都长什么样。

mu_params = [-1, 0, 1]
sd_params = [0.5, 1, 1.5]
x = np.linspace(-7, 7, 200)
_, ax = plt.subplots(len(mu_params), len(sd_params), sharex=True,
                     sharey=True,
                     figsize=(9, 7), constrained_layout=True)
　
for i in range(3): 
    for j in range(3):
        mu = mu_params[i]
        sd = sd_params[j]
        y = stats.norm(mu, sd).pdf(x) 
        ax[i,j].plot(x, y)
       ax[i,j].plot([], label="μ = {:3.2f}\nσ = {:3.2f}".format(mu, sd),                          alpha=0)
        ax[i,j].legend(loc=1) 
ax[2,1].set_xlabel('x')
ax[1,0].set_ylabel('p(x)', rotation=0, labelpad=20)
ax[1,0].set_yticks([])

上面的代码主要是为了画图，其中与概率相关的部分是y = stats.norm(mu, sd).pdf(x)，x服从参数为mu和sd的正态分布，这行代码的作用是为一组x值计算其对应的概率密度。上面的代码会生成图1.1，每个子图中，用蓝色的线表示不同和值对应的正态分布。

图1.1

提示：本书大部分图是通过图前面的代码直接生成的，尽管有些时候并没有直接说明。对熟悉Jupyter Notebook或者Jupyter Lab的读者来说应该会很熟悉。

随机变量分为两种，连续变量和离散变量。连续随机变量可以从某个区间内取任意值（可以用Python中的浮点型数据来表示），而离散随机变量则只能取某些特定的值（可以用Python中的整型数据来描述）。正态分布属于连续分布。

需要注意的是，生成图1.1的代码中省略了yticks，这是一个功能而非Bug。省略的原因是，这些值并没有提供太多有用的信息，反而有可能使某些人产生困惑。让我来解释下，y轴的数据其实不太重要，重要的是它们的相对值。假设从变量x中取两个值，比方说x_i和x_j，然后发现（在图中两倍高），我就可以放心地说，x_i的概率值是x_j的两倍，这就是大多数人直观的理解（幸运的是，这也是正确的理解）。唯一棘手的地方在于，对于连续分布，y轴上的值并非表示概率值，而是概率密度。想要得到概率值，需要在给定区间内做积分，也就是说，需要计算（针对该区间的）曲线下方的面积。虽然概率不能大于1，但是概率密度却可以大于1，只不过其密度曲线下方的总面积限制为1。从数学的角度来看，理解概率和概率密度之间的差别非常重要。不过本书采用的是偏实践的方法，可以稍微简单些，毕竟只要你能根据相对值解释图表，这种差别就没有那么重要了。

很多模型假设随机变量的连续值都是从同一分布中采样的，并且这些值相互独立。在这种情况下，我们称这些变量为独立同分布变量。用数学符号来描述的话就是：如果两个随机变量x和y对于所有可能的取值都满足，那么称这两个变量相互独立。

非独立同分布变量的一个常见示例是时间序列，其中随机变量中的时间相关性是一个应该考虑的关键特性。下面例子中的数据记录了从1959年到1997年大气中二氧化碳的含量。可以用代码把图画出来。

data = np.genfromtxt('../data/mauna_loa_CO2.csv', delimiter=',')
plt.plot(data[:,0], data[:,1]) 
plt.xlabel('year') 
plt.ylabel('$CO_2$ (ppmv)')
plt.savefig('B11197_01_02.png', dpi=300)

图1.2

图1.2中的数据对应每个月测得的大气中二氧化碳的含量，可以很明显地看出数据的时间相关性。实际上，这里有两种趋势：一种是季节性的，这与植物周期性生长和衰败有关；另一种是整体性的，这表明大气中二氧化碳浓度在不断增加。

现在我们已经学习了概率论中的一些基本概念和术语，接下来就是激动人心的时刻了。话不多说，让我们以庄严的姿态思考贝叶斯定理：

　　（1.4）

这看起来稀松平常，似乎跟小学课本里的公式差不多？用Richard Feynman的话来说：“这就是关于贝叶斯统计的所有知识。”

学习贝叶斯定理的来源将有助于我们理解它的含义。根据式（1.1）有：

　　（1.5）

也可以写成：

　　（1.6）

假设式（1.5）和式（1.6）的左边是相等的，那么可以合并起来得到下面的式子：

　　（1.7）

重新把式（1.7）组织一下，就得到了式（1.4），即贝叶斯定理。

现在让我们看看式（1.4）的含义，理解它为什么重要。首先，它表明不一定和是一样的。这一点非常重要，在日常分析中，即使是系统学习过统计学和概率论的人也很容易忽略这点。我们举个简单例子来说明为什么二者不一定相同：一个阿根廷人成为教皇的概率和一个教皇是阿根廷人的概率并不相同。因为假设有44 000 000个阿根廷人，而只有一个是教皇，所以有p(教皇｜阿根廷人)≈1/44 000 000而p(阿根廷人｜教皇)=1。

如果用假设代替，用数据代替y，那么贝叶斯定理告诉我们的是，给定数据y，如何计算一个假设 的概率。你会在很多地方见到这种关于贝叶斯定理的解释，但是我们怎么把一个假设转换成可以放入贝叶斯定理中的东西呢？可以通过概率分布来实现。通常来讲，我们所说的假设是一个非常狭隘的假设。如果我们讨论为模型中的参数找到一个合适的值，即概率分布的参数，这会更精确。顺便说一点，不要试图把 设定为诸如“独角兽真的存在”等这类命题，除非你真的愿意构建一个关于独角兽真实存在的概率模型！

贝叶斯定理是贝叶斯统计的核心，我们将在第2章中看到，使用PyMC3等工具可以让我们在每次构建贝叶斯模型的时候都不必显式地编写贝叶斯定理。尽管如此，了解贝叶斯定理各个部分的名字还是非常重要的，因为后面它们会被反复提及。此外，理解各个部分的含义也有助于对模型进行概念化。

● ：先验。

● ：可能性、似然。

● ：后验。

● ：边缘似然、证据。

先验分布反映的是在观测到数据y之前对参数 的了解。如果我们对参数一无所知，那么可以用扁平先验（它不会传递太多信息）。读完本书后你会发现，通常我们能做到的要比扁平先验更好。为什么有些人会认为贝叶斯统计是主观的，原因就是使用了先验，先验不过是构建模型时的另一个假设，与任何其他假设是一样的主观（或客观）。

似然是指如何在分析中引入数据，它表达的是给定参数的数据合理性。在一些文章中，你会发现有些人也称之为采样模型、统计模型或者就叫模型。本书坚持使用似然这一名称。

后验分布是贝叶斯分析的结果，反映的是（在给定数据和模型的条件下）我们对问题的全部了解。后验指的是模型中参数的概率分布而不是单个值，这种分布是先验与似然之间的平衡。有这么个有名的笑话：“‘贝叶斯人’是这样的：他模糊地期待着一匹马（先验），并瞥见了一头驴（边缘似然），结果他坚信看到的是一头骡子（后验）。”对这个笑话的解释是如果似然和先验都是模糊的，那么也会得到一个模糊的后验。不管怎么说，我喜欢这个笑话，因为它讲出了这样一个道理，后验其实是对先验和似然的某种折中。从概念上讲，后验可以看作在观测到数据之后对先验的更新。事实上，一次分析中的后验，在收集到新的数据之后，也可以看作下一次分析中的先验。这使得贝叶斯分析特别适用于序列化的数据分析，比如实时处理来自气象站和卫星的数据以进行灾害预警，要了解更详细的内容，可以阅读机器学习方面的算法。

最后一个概念是边缘似然，也称作证据。正式地讲，边缘分布是在模型的参数取遍所有可能值的条件下得到指定观测值的概率的平均。不过，本书的大部分内容并不关心这个概念，我们可以简单地把它当作归一化系数。这么做没什么大问题，因为我们只关心参数的相对值而非绝对值。你可能还记得我们在1.2.2节讨论如何解释概率分布图时提到过这一点。把证据这一项忽略掉之后，贝叶斯定理可以表示成如下正比例形式：

　　（1.8）

理解贝叶斯定理中的每个概念可能需要点儿时间和更多的例子，本书剩余部分也将围绕这些内容展开。

1.3　单参数推断

在前面两节中，我们学习了几个重要的概念，其中有两个是贝叶斯统计的核心概念，这里我们用一句话再重新强调一遍。

提示：概率是用来衡量参数的不确定性的，而贝叶斯定理是用来在观测到新的数据时正确更新这些概率以期降低我们的不确定性。

现在我们已经知道什么是贝叶斯统计了，接下来就从一个简单的例子入手，通过推断单个未知参数来学习如何做贝叶斯统计。

抛硬币问题

抛硬币问题是统计学中的一个经典问题，描述如下：随机抛一枚硬币，重复一定次数，记录其正面朝上和反面朝上的次数，根据这些数据，回答诸如“这枚硬币是否”，以及“这枚硬币有多不公平”等问题。这些问题看起来似乎有点儿无聊，不过可别低估了它。抛硬币问题是一个学习贝叶斯统计非常好的例子，一方面是因为几乎人人都熟悉抛硬币这一过程，另一方面是因为这个模型很简单，我们可以很容易计算并解决这个问题。此外，很多真实问题都包含两个互斥的结果，例如0或者1、正或者负、奇数或者偶数、垃圾邮件或者正常邮件、安全或者不安全、健康或者不健康等。因此，即便我们讨论的是硬币，这个模型也同样适用于前面那些问题。

为了估计硬币的偏差，或者更广泛地说，想要用贝叶斯定理解决问题，我们需要数据和一个概率模型。对于抛硬币这个问题，假设我们已经验了一定次数并且记录了正面朝上的次数，也就是说数据部分已经准备好了，剩下的就是模型部分了。考虑到这是第一个模型，我们会列出所有必要的数学公式（别怕，我保证这部分不难），并且一步一步推导。第2章中，我们会回顾这个问题，并借用PyMC3和计算机来完成数学计算部分。

首先，要抽象出偏差的概念。如果一枚硬币总是正面朝上，那么我们说它的偏差就是1；如果总是反面朝上，那么它的偏差就是0；如果正面朝上和反面朝上的次数各占一半，那么它的偏差就是0.5。这里用参数来表示偏差，用y表示N次抛硬币实验中正面朝上的次数。根据贝叶斯定理，即式（1.4），首先需要指定先验和似然。让我们先从似然开始。

假设多次抛硬币的结果相互之间都没有影响，也就是说每次抛硬币都是相互独立的，同时还假设结果只有两种可能，正面朝上或者反面朝上。更进一步，我们假设所有硬币都服从相同的分布，因此，抛硬币过程中的随机变量是一个典型的独立同分布变量。但愿你能认同我们对这个问题做出的合理假设。基于这些假设，一个不错的似然候选是二项分布：

　　（1.9）

这是一个离散分布，表示的是N次抛硬币实验中y次正面朝上的概率（或者更一般的描述是，N次实验中，y次成功的概率）。

n_params = [1, 2, 4]  # 实验次数
p_params = [0.25, 0.5, 0.75]  # 成功概率
　
x = np.arange(0, max(n_params)+1)
f,ax = plt.subplots(len(n_params), len(p_params), sharex=True, 
                    sharey=True, 
                    figsize=(8, 7), constrained_layout=True)
　
for i in range(len(n_params)): 
    for j in range(len(p_params)):
        n = n_params[i] 
        p = p_params[j]
　
        y = stats.binom(n=n, p=p).pmf(x)
　
        ax[i,j].vlines(x, 0, y, colors='C0', lw=5)
        ax[i,j].set_ylim(0, 1)
        ax[i,j].plot(0, 0, label="N = {:3.2f}\nθ = 
                           {:3.2f}".format(n,p), alpha=0)
        ax[i,j].legend()
        ax[2,1].set_xlabel('y')
        ax[1,0].set_ylabel('p(y | θ, N)')
        ax[0,0].set_xticks(x)
        plt.savefig('B11197_01_03.png', dpi=300)

图1.3

图1.3展示了9个二项分布，每个子图中的标签显示了对应的参数。注意，在这个图中，我并没有省略y轴的值，这么做的目的是，读者能够自己动手确认每个子图的y值之和为1。也就是说，对离散分布而言，对应的y值就是概率值。

似然使用二项分布是一个合理选择，直观上讲，可以看作抛一次硬币时正面朝上的可能性（对于N=1，这很直观，不过对于任意的N值都适用），你可以将与图中y =1的柱状图对比。

假如知道了，那么就可以从二项分布得出硬币正面朝上的分布。关键是我们不知道！不过别灰心，在贝叶斯统计中，每当我们不知道某个参数的时候，就对它赋予一个先验，接下来继续选择先验。

这里我们选用贝叶斯统计中最常见的一个分布——贝塔分布，作为先验，其数学形式如下：

　　（1.10）

仔细观察式（1.10）可以看出，除了部分之外，贝塔分布和二项分布看起来很像。是希腊字母中大写的伽马，用来表示伽马函数。现在我们只需要知道，用分数表示的第一项是一个正则化常量，用来保证该分布的积分为1，此外和两个参数用来控制具体的分布形态。贝塔分布是我们到目前为止见到的第三个分布，利用下面的代码，可以深入了解其形态。

params = [0.5, 1, 2, 3]
x = np.linspace(0, 1, 100)
f, ax = plt.subplots(len(params), len(params), sharex=True, 
                     sharey=True, 
                     figsize=(8, 7), constrained_layout=True)
for i in range(4): 
    for j in range(4):
        a = params[i] 
        b = params[j]
        y = stats.beta(a, b).pdf(x)
        ax[i,j].plot(x, y)
        ax[i,j].plot(0, 0, label="α = {:2.1f}\nβ = {:2.1f}".format(a, 
                     b), alpha=0)
        ax[i,j].legend()
ax[1,0].set_yticks([])
ax[1,0].set_xticks([0, 0.5, 1])
f.text(0.5, 0.05, 'θ', ha='center')
f.text(0.07, 0.5, 'p(θ)', va='center', rotation=0)

图1.4

我很喜欢贝塔分布以及图1.4中的各种形状，但是我们为什么要在模型中用它呢？在抛硬币这个问题以及一些其他问题中使用贝塔分布的原因有很多，其中第一个原因是，贝塔分布的范围限制在0到1，这跟参数一样。通常，在需要对二值变量建模时，就会选用贝塔分布。第二个原因是其通用性，从图1.4可以看出，该分布可以有多种形状，包括均匀分布、类高斯分布、标准正态分布等。第三个原因是，贝塔分布是二项分布（之前我们用它描述似然）的共轭先验。似然的共轭先验是指，将这个先验分布与似然组合在一起之后，得到的后验分布与先验分布的表达式形式仍然是一样的。简单点儿说，就是每次用贝塔分布作为先验、二项分布作为似然时，会得到一个贝塔分布的后验。除贝塔分布之外还有很多其他共轭先验，例如高斯分布，它的共轭先验就是它自己。很多年来，贝叶斯分析都限定在共轭先验范围内，这主要是因为共轭能让后验在数学上变得更容易处理，要知道贝叶斯统计中一个常见问题的后验都很难从分析的角度去解决。在建立合适的计算方法来解决任意后验之前，这只是个折中的办法。从第2章开始，我们将学习使用现代的计算方法来解决贝叶斯问题而不必考虑是否使用共轭先验。

首先回忆一下，如贝叶斯定理，即式（1.4）所述，后验正比于似然乘先验。对于我们的问题，需要将二项分布乘贝塔分布：

　　（1.11）

现在，将上式简化。针对我们的实际问题，可以先把与不相关的项去掉而不影响结果，于是得到下式：

　　（1.12）

重新整理之后得到：

　　（1.13）

细心的读者可以看出式（1.13）和贝塔分布的形式很像（除了归一化部分），其对应的参数分别为和。也就是说，在抛硬币这个问题中，后验分布是如下贝塔分布：

　　（1.14）

现在已经有了后验的表达式，可以用Python对其计算并画出结果。下面的代码中，其实只有一行是用来计算后验结果的，其余的代码都是用来画图的。

plt.figure(figsize=(10, 8))
　
n_trials = [0, 1, 2, 3, 4, 8, 16, 32, 50, 150]
data = [0, 1, 1, 1, 1, 4, 6, 9, 13, 48]
theta_real = 0.35
　
beta_params = [(1, 1), (20, 20), (1, 4)] 
dist = stats.beta
x = np.linspace(0, 1, 200)
　
for idx, N in enumerate(n_trials): 
    if idx == 0:
        plt.subplot(4, 3, 2) 
        plt.xlabel('θ')
    else:
        plt.subplot(4, 3, idx+3) 
        plt.xticks([])
    y = data[idx]
    for (a_prior, b_prior) in beta_params:
        p_theta_given_y = dist.pdf(x, a_prior + y, b_prior + N - y) 
        plt.fill_between(x, 0, p_theta_given_y, alpha=0.7)
　
    plt.axvline(theta_real, ymax=0.3, color='k')
    plt.plot(0, 0, label=f'{N:4d} trials\n{y:4d} heads', alpha=0) 
    plt.xlim(0, 1)
    plt.ylim(0, 12) 
    plt.legend() 
    plt.yticks([])
plt.tight_layout()

在图1.5的第一个子图中，还没做任何试验，对应的3条曲线分别表示3种先验。

图1.5

● 蓝色表示均匀分布的先验，其含义是，偏差的所有可能取值都是等概率的。

● 橙色表示类高斯分布的先验，集中在值0.5附近，该先验反映了这样一种信息：通常硬币正面朝上和反面朝上的概率大致是差不多的。我们也可以说，该先验与“大多数硬币是公平的”这一信念是相符合的。虽然信念这个词在贝叶斯相关的讨论经常用到，但我们最好是讨论受数据影响的模型和参数。

● 绿色表示左偏的先验，更倾向于反面朝上。

其余子图描绘了后续实验的后验分布，每个子图都标注了实验（抛硬币）的次数和正面朝上的次数。此外每个子图中在横轴0.35附近还有一条黑色的竖线，表示的是的真实值。显然，在实际情况下，我们并不知道这个值，在这里标识出来只是方便理解。从图1.5中可以学到很多贝叶斯分析方面的知识，我们花点时间理解。

● 贝叶斯分析的结果是后验分布而不是某个值，该分布描述了不同数值的可能性，是根据给定数据和模型得到的。

● 后验最可能的值是由后验分布的众数决定的（也就是后验分布的峰值）。

● 后验分布的展开度与我们对参数的不确定性成正比例；分布越离散，不确定性越大。

● 直观上讲，观测到的数据越多，那么我们对某一结果就更有信心，尽管1/2=4/8=0.5，但是相比在4次试验中观测到2次正面朝上，在8次试验中观测到4次正面朝上，会让我们对偏差等于0.5这一描述更有信心。该直觉也同时反映在了后验分布上。你可以从图1.5中确认这一点，尤其注意第6个子图中的蓝色后验，尽管众数是一样的，但是第3个子图中曲线的离散程度（不确定性）相比第6个子图要更大。

● 在给定足够多的数据时，两个或多个不同先验的贝叶斯模型会趋近于收敛到相同的结果。在极限情况下，如果有无限多的数据，不论我们使用的是怎样的先验，最终都会得到相同的后验。注意这里说的无限多是指一个极限状态而非某个具体的数量，也就是说，从实际的角度来讲，通过有限且数量较少的数据点就可以得到近似的后验。

● 不同后验收敛到相同分布的速度取决于数据和模型。从图1.5中可以看出，蓝色的先验（均匀分布）和绿色的先验（左偏分布）在经过8次实验之后就很难看出区别了，而橙色的后验（类似高斯分布）即使进行了150次实验之后还能很容易地看出与另外两个后验的区别。

● 有一点从图1.5中不太容易看出来，如果我们一步一步地更新后验，最后得到的结果其实跟一次性计算得到的结果是一样的。换句话说，我们可以对后验进行150次计算，每次增加一个新的观测数据并将得到的后验作为下一次计算的先验，也可以在得到150次抛硬币的结果之后一次性计算出后验，而这两种计算方式得到的结果是完全一样的。这个特点非常有意义，当我们得到新的数据时，就可以很自然地更新估计值。这在很多数据分析问题中都很常见。

从前面的例子可以明显看出，先验对推理会有影响，这很正常，先验就是这样。一些贝叶斯分析的新手（以及一些诋毁该方法的人）会对如何选择先验感到茫然，因为他们不希望先验起到决定性作用，而是更希望数据本身替自己说话！有这样的想法很正常，不过我们得牢记，数据并不会真地“说话”，最多也只是发出“嗡嗡”声。数据只有在模型中才会有意义，包括数学模型和心理模型。面对同一主题下的同一份数据，不同人会有不同的看法，这类例子在科学史上有很多，即使你基于的是正式的模型，这种情况也还是会发生。

有些人青睐于使用无信息先验（也称作扁平先验^[2]、模糊先验或扩散先验），这类先验对分析过程的影响最小。尽管这么做是可行的，不过通常我们能做得比这更好。本书将遵循Gelman、McElreath和Kruschke^[3]这3位的建议，更倾向于使用弱信息先验（weakly-information prior）。在很多问题中，我们对参数可以取的值一般都会有些了解。比如，我们可能知道参数只能是正数，或者是知道参数近似的取值范围，又或者是希望该值接近0或大于/小于某值。这种情况下，我们可以给模型加入一些微弱的先验信息而不必担心它会掩盖数据本身的信息。由于这类先验会让后验近似位于某一合理的边界内，因此也被称作正则化先验。当然，如果已经有一些高质量的信息用于定义先验，那么使用带有较多信息量的强先验也是可行的。视具体的问题不同，有可能很容易或者很难找到这类先验，例如在我工作的领域（结构生物信息学），人们会尽可能地利用先验信息，通过贝叶斯或者非贝叶斯的方式来了解和预测蛋白质的结构。这样做是合理的，原因是我们在数十年间已经从上千次精心设计的实验中收集了数据，因而有大量可信的先验信息可供使用，如果不用的话也太荒唐了！所以，如果你有可信的先验信息，完全没有理由不去使用。试想一下，如果一个汽车工程师每次设计新车的时候，他都要重新发明内燃机、轮子乃至整个汽车！这显然不是正确的方式。

[2]　扁平先验（flat prior）用均匀分布来表示某个参数的先验分布，我们对参数的取值没有稳定的趋向，会认为参数取任何值的概率都相同，因此扁平先验也是无信息先验（non-information prior），即不提供任何信息的先验。

[3]　这3位分别是Bayesian Data Analysis、Statistical Rethinking: A Bayesian Course with Examples in R and Stan和Doing Bayesian Data Analysis的主要作者。——译者注

现在我们知道了先验有很多种，不过这并不能缓解我们选择先验时的焦虑。或许，最好是没有先验，这样事情就简单了。不过，不论是否基于贝叶斯，模型都在某种程度上拥有先验，即使这里的先验并没有明确表示出来。事实上，很多频率统计学方面的结果在某些情况下可以看作贝叶斯模型的特例，比如扁平先验。在频率学派中，一种常见的参数估计方法是最大似然法，该方法避免了设置先验，其原理就是找到合适的让似然最大化，该值通常用一个带尖的符号来表示，比如 或者。 是一个点估计（一个数值），而不是一个分布，对于前面抛硬币的问题，可以表示为：

　　（1.15）

如果你重新看图1.5，你可以亲自确认，蓝色后验（对应均匀/扁平先验）的众数是否与每个子图计算的 值相同。因此，至少对这个例子而言，尽管最大似然的方法并没有显式地引入任何先验，但仍然可以看作贝叶斯模型的一种（带均匀先验的）特例。

我们无法真正避免先验，不过如果你在分析中引入先验，得到的是合理值的分布而不只是最可能的一个值。明确引入先验的另一个好处是，我们会得到更透明的模型，这意味着这些模型更容易评判、调试（广义上）以及优化。构建模型是一个迭代的过程，有时候可能只需要几分钟，有时候则可能需要数年；有时候整个过程可能只涉及你本人，有时候则可能涉及你不认识的人。另外，模型复现很重要，而模型中透明的假设能有助于复现。在特定分析任务中，如果我们对某个先验或者似然不确定，可以使用多个先验或者似然进行尝试。模型构建过程中的一个环节就是质疑假设，而先验就是质疑的对象之一。不同的假设会得到不同的模型，根据数据和与问题相关的领域知识，我们可以对这些模型进行比较，本书第5章会深入讨论该部分内容。由于先验是贝叶斯统计中的一个核心内容，在接下来遇到新的问题时我们还会反复讨论它，因此如果你对前面讨论的内容感到有些疑惑，别太担心，先冷静冷静，要知道人们在这个问题上已经困惑了数十年并且相关的讨论一直在继续。

1.4　报告贝叶斯分析结果

创建报告并交流分析结果是统计学和数据科学中的核心环节。本节将简要介绍对采用贝叶斯模型得到的分析结果进行报告和交流的特别之处。在接下来的几章中，还会使用更多例子来介绍这部分内容。

1.4.1　模型表示和可视化

如果你想与人交流分析结果，那么同时你还需要与人交流你所使用的模型。以下是一种常用的概率模型表示法：

　　（1.16）

这是我们抛硬币例子里用到的模型。回忆一下，符号～表示左边随机变量的分布服从右边的分布形式，也就是说，这里服从于参数为和的贝塔（beta）分布，而y服从于参数为和的二项（binomial）分布。该模型还可以用Kruschke图表示，如图1.6所示。

图1.6

在第一层，根据先验生成了，然后通过似然生成最下面的数据。图中的箭头表示变量之间的依赖关系，符号～表示变量的随机性。本书中用到的Kruschke图都是由Rasmus Bååth提供的模板生成的，在这里我非常感谢他提供的这些模板。

1.4.2　总结后验

贝叶斯分析的结果是后验分布，关于给定数据集和模型的参数的所有信息都包含在后验分布中。因此，通过总结后验数据，可以总结模型和数据的逻辑结果。通常的做法是，报告每个参数的均值（或者众数、中位数），以了解分布的位置，以及一些度量值，如标准差，可以了解分布的离散程度从而了解我们估计中的不确定性。标准偏差适用于正态分布，不过对于其他类型的分布（如偏态分布）却可能得出误导性结论，因此，还可以采用以下度量方式。

最大后验密度（Highest-Posterior Density，HPD）区间常用于描述后验分布的分散程度。HPD区间是指包含一定比例概率密度的最小区间，最常见的比例是95%HPD或98%HPD，通常还伴随着一个50%HPD。如果我们说某个分析的HPD区间是[2, 5]，其含义是指，根据我们的模型和数据，参数位于2到5的概率是95%。

提示：选择95%还是50%或者其他值作为HPD区间的概率密度比例并没有特殊的地方，这些不过是常用的值罢了。如果愿意，我们完全可以选用比例为91.37%的HPD区间。如果你选的是95%，这完全没问题，只是要记住这只是个默认值，究竟选择多大比例仍然需要具体问题具体分析。

ArviZ是一个用于对贝叶斯模型做探索性数据分析的Python库，提供了很多方便的函数用于对后验做总结，比如，az.plot_posterior可以用来生成一个带有均值和HPD区间的分布图。下面用一个从贝塔分布中生成（并非实际分析中的后验分布）的随机样本为例来说明。

np.random.seed(1)
az.plot_posterior({'θ':stats.beta.rvs(5, 11, size=1000)})

图1.7

注意，图1.7中报告的区间是94% HPD，这是为了善意地提醒你，HPD可以是除了95%之外的任意值。每当ArviZ计算并报告HPD时，默认都会使用0.94（对应94%区间），你可以修改参数credible_interval的值来改变这一默认行为。

提示：如果你了解频率学派，那么需要注意，HPD区间和置信区间并不一样。HPD区间理解起来很直观，以至于人们往往容易将置信区间与HPD区间弄混。贝叶斯分析能让我们考虑参数取不同值的概率，而这对于频率学派是不可能的，因为参数在设计的时候就是固定的，频率学派的置信区间的含义是，某一区间是否包含参数的值的概率。

1.5　后验预测检查

贝叶斯方法的一个优势是，一旦得到了后验分布，就可以基于数据集y和估计到的参数，使用后验分布生成待预测的数据。后验预测分布如下：

　　（1.17）

可以看到，后验预测分布就是对的后验分布进行条件预测的平均。从概念上（以及计算的角度）讲，可以通过以下两步来估计式（1.17）中的积分。

（1）从后验分布中采样出。

（2）将“喂”给似然（或者叫采样分布），然后得到一个数据点。

提示：注意，这个过程组合了两种不确定性：一种是参数的不确定性，即通过后验来刻画的；另一种是采样的不确定性，即通过似然来刻画的。

生成的数据可以在需要预测数据的时候使用。此外还可以用预测数据和观测数据y做对比，找出两组数据之间的差异，从而对模型进行评价，这就是所谓的后验预测检查，其目的是检查自动一致性。生成的数据应该与观测数据看起来差不多，否则就说明建模的过程中出了问题，或者是向模型提供数据时出了问题。不过即使没有出错，二者也有可能不同。尝试理解其中的偏差有助于我们改进模型，或者至少能知道模型的极限。即使我们并不知道如何改进模型，也可以知道模型捕捉到了问题或数据的哪些方面以及没能捕捉到哪些方面，也许模型能够很好地捕捉到数据中的均值但没法预测出罕见值，这可能是个问题，不过如果我们只关心均值，那这个模型也还是可用的。通常我们的目的不是去断言一个模型是错误的，我们只想知道模型的哪个部分是值得信任的，并测试它是否在特定方面符合我们的预期。不同学科对模型的信任程度显然是不同的，物理学中研究的系统是在高可控条件下依据高级理论运行的，因而模型可以看作对现实的良好描述，而在一些其他学科（如社会学和生物学）中，研究的是错综复杂的孤立系统，因而模型对系统的认知较弱。尽管如此，不论你研究的是哪一门学科，都需要对模型做检查，利用后验预测和本章学到的探索性数据分析中的方法去检查模型是个不错的习惯。

1.6　总结

本章围绕统计建模、概率论和贝叶斯定理开启了我们的贝叶斯之旅，然后用抛硬币的例子介绍了贝叶斯建模和数据分析，借用这个经典例子传达了贝叶斯统计中的一些最重要的思想，比如通过概率分布构建模型并用它来表示不确定性。我们试图揭开先验的神秘面纱，并将其与建模过程中的其他元素（如可能性）或更多元问题（如我们为什么要首先解决某个特定问题）放在同等的地位。本章的最后讨论了如何解释和报告贝叶斯分析的结果。

图1.8来自Sumio Watanabe，很好地总结了本章描述的贝叶斯建模过程。

图1.8

我们假设存在一个真实分布，该分布通常是未知的（原则上也是不可知的），我们可以通过实验、调查、观察或模拟从中获得有限样本。为了从真实分布中学到些东西，假设我们只观测到了一个样本，并构建了一个概率模型。一个概率模型通常由先验和似然两部分构成。用模型和样本做贝叶斯推断并得到一个后验分布；在给定的模型和数据下，该分布囊括有关某个具体问题的全部信息。从贝叶斯的角度来看，后验分布及其衍生的其他信息（包括后验预测分布）是我们主要感兴趣的部分。由于后验分布（及其衍生值）是由模型和数据产生的，因此贝叶斯推断的作用受限于模型和数据的质量。一种评价模型的方法是，将后验预测分布与有限的观测数据做对比。这里注意区分，后验分布是模型中参数的分布（以观测到的样本为条件），而后验预测分布则是预测样本的分布（在后验分布上平均）。模型验证的过程非常重要，因为我们希望确保得到的模型是正确的。不过我们也知道，永远没有所谓的正确的模型。之所以做模型检查是为了确定我们得到的模型在某个具体的场景是否足够有用，如果不是，那就深入了解如何改善它们。

本章简要总结了贝叶斯数据分析的几个主要方面，本书其余章节中，我们以它为基础去理解更高阶的概念，并反复回顾这些概念，最终达到真正掌握它们的目的。第2章会介绍PyMC3和ArviZ这两个库，PyMC3是一个用于贝叶斯建模和概率机器学习的Python库，ArviZ是用于探索性分析贝叶斯模型的Python库。

1.7　练习

我们尚且不清楚大脑是如何运作的，是按照贝叶斯方式？还是类似贝叶斯的某种方式？又或者是进化过程中形成的某种启发式的方式？不管如何，我们至少知道自己是通过数据、例子和练习来学习的。你可能对此有异议，但我仍然强烈建议你完成每章最后的练习。

（1）下面的表达式中，哪一个与“1816年7月9日是晴天的概率”这一描述相符？

● p(晴天)

● p(晴天｜7月)

● p(晴天｜1816年7月9日)

● p(1816年7月9日｜晴天)

● p(晴天,1816年7月9日)/p(1816年7月9日)

（2）证明“随机选一个人是教皇”的概率与“教皇是人类”的概率不同。在《飞出个未来》这部动画片中，教皇是一个爬行动物，这会改变你之前的计算吗？

（3）指出下面定义的概率模型中的先验和似然：

（4）在上面的模型中，后验有多少个参数？把它和抛硬币问题中的模型做对比。

（5）请为练习（3）中的模型写出贝叶斯定理。

（6）假设有两枚硬币，抛第一个硬币的时候，一半概率正面朝上，另一半概率反面朝上。另外一枚硬币则总是正面朝上。如果随机抽取两枚硬币中的一个抛出并观测到其正面朝上，那么这枚硬币是第二枚硬币的概率是多少？

（7）修改生成图1.5的代码，给每个子图添加一条虚的竖线，用来表示观测到的正面朝上出现的比例，将其与每个子图中后验分布的众数进行比较。

（8）尝试用一些其他的先验（比如贝塔分布）和一些其他的数据重绘图1.5。

（9）换一些参数重新绘制高斯分布、二项分布和贝塔分布（对应图1.1、图1.3以及图1.4），当然你也可以只画一个图而不必像我那样画出多个子图的网格。

（10）了解克伦威尔准则的相关内容。

（11）了解荷兰赌（Dutch Book）定理的相关内容。

读者服务：

微信扫码关注【异步社区】微信公众号，回复“e60089”获取本书配套资源以及异步社区15天VIP会员卡，近千本电子书免费畅读。