详情
本书首先介绍贝叶斯框架的关键概念，而后探讨广义线性模型的强大功能和灵活性，以及如何把它们用到一系列广泛的问题中，包括回归和分类。同时，本书也探讨了一些高级课题，包括混合模型和高斯过程等。无论数据科学的新手，还是有经验的专业人士，都可以从本书受益。
书名：Python贝叶斯分析
ISBN：978-7-115-47617-3
本书由人民邮电出版社发行数字版。版权所有，侵权必究。
您购买的人民邮电出版社电子书仅供您个人使用，未经授权，不得以任何方式复制和传播本书内容。
我们愿意相信读者具有这样的良知和觉悟，与我们共同保护知识产权。
如果购买者有侵权行为，我们可能对该用户实施包括但不限于关闭该帐号等维权措施，并可能追究法律责任。
• 著 [阿根廷]奥斯瓦尔多•马丁（Osvaldo Martin）
译 田 俊
责任编辑 王峰松
• 人民邮电出版社出版发行 北京市丰台区成寿寺路11号
邮编 100164 电子邮件 315@ptpress.com.cn
• 读者服务热线：(010)81055410
反盗版热线：(010)81055315
Copyright © 2016 Packt Publishing. First published in the English language under the title Bayesian Analysis with Python, ISBN 978-1-78588-380-4. All rights reserved.
本书中文简体字版由Packt Publishing公司授权人民邮电出版社出版。未经出版者书面许可，对本书的任何部分不得以任何方式或任何手段复制和传播。
版权所有，侵权必究。
本书从务实和编程的角度讲解了贝叶斯统计中的主要概念，并介绍了如何使用流行的PyMC3来构建概率模型。阅读本书，读者将掌握实现、检查和扩展贝叶斯统计模型，从而提升解决一系列数据分析问题的能力。本书不要求读者有任何统计学方面的基础，但需要读者有使用Python编程方面的经验。
Osvaldo Martin是阿根廷国家科学与技术理事会（CONICET）的一名研究员。该理事会是负责阿根廷科技进步的主要组织。Osvaldo Martin曾从事结构生物信息学和计算生物学方面的研究，特别是在如何验证结构蛋白质的模型方面有着深入研究；此外，在应用马尔科夫—蒙特卡洛方法模拟分子方向上有着丰富的经验，尤其喜欢用Python解决数据分析问题。他曾讲授结构生物信息学、Python编程等课程，最近还开设了贝叶斯数据分析的课程。Python和贝叶斯统计改变了他对科学的认知和对问题的思考方式。他写本书的最大动力是希望借助Python帮助大家理解概率模型，同时，他也是PyMOL社区（一个基于C/Python的分子可视化社区）的活跃成员，最近他也对PyMC3社区做了一些贡献。
我要感谢我的妻子Romina在我写作本书的过程中对我的支持，以及对我所有“疯狂”项目的支持。我还要感谢Walter Lapadula、Juan Manuel Alonso和Romina Torres-Astorga，他们为本书提供了宝贵的反馈和建议。
此外，我还要特别感谢PyMC3的核心开发者们，本书正是因为他们的辛勤奉献才得以成为可能，希望本书能够促进PyMC3的传播和使用。
Austin Rochford是Monetate Labs的首席数据科学家，他开发的产品用于帮助零售商在每年数十亿的活动中进行个性化销售。同时他还是一位积极倡导贝叶斯方法的数学家。
田俊，计算机专业硕士。2016年毕业于中国科学院自动化研究所，主要研究方向为自然语言处理中的短文本分类，毕业后曾在滴滴出行担任算法工程师，目前在微软从事自然语言处理方面的工作。
劳俊鹏，心理学博士，PyMC团队成员。2014年毕业于英国格拉斯哥大学，主要研究认知神经心理学。2013年至今在瑞士弗里堡大学从事心理学研究，专攻数据建模分析和神经计算模型。
贝叶斯统计距今已经有超过250年的历史，其间该方法既饱受赞誉又备受轻视。直到近几十年，得益于理论进步和计算能力的提升，贝叶斯统计才越来越多地受到来自统计学以及其他学科乃至学术圈以外工业界的重视。现代的贝叶斯统计主要是计算统计学，人们对模型的灵活性、透明性以及统计分析结果的可解释性的追求最终造就了该趋势。
本书将从实用的角度来学习贝叶斯统计，不会过多地考虑统计学范例及其与贝叶斯统计之间的关系。本书的目的是借助Python做贝叶斯数据分析，尽管与之相关的哲学讨论也很有趣，不过受限于篇幅，这部分内容并不在本书的讨论范围之内，有兴趣的读者可以通过其他方式深入了解。
这里我们采用建模的方式学习统计学，学习如何从概率模型的角度思考问题，结合模型和数据利用贝叶斯理论推导出逻辑结论。这种建模方式使用的是一种数值计算的方法，其中，模型部分会由基于PyMC3的代码构建。PyMC3是用于贝叶斯统计的Python库，它为用户封装了大量的数学细节和计算过程。贝叶斯方法在理论上源于概率论，这也是为什么许多讲贝叶斯方法的书中都充斥着需要一定数学基础的公式。学习统计学方面的数学知识显然有利于构建更好的模型，同时还能让你对问题、模型和结果有更好的直觉。不过类似PyMC3的库能够帮助你在有限的数学知识水平下学习并掌握贝叶斯统计。在阅读本书的过程中，你将亲自见证这一过程。
第1章，概率思维——贝叶斯推断指南，介绍了贝叶斯理论及其在数据分析中的意义，并进一步阐述了贝叶斯思维方式的定义以及为什么和如何使用概率来处理不确定性。本章还包含本书其余章节中的一些基本概念。
第2章，概率编程——PyMC3编程指南，从计算的角度重新回顾了前一章提到的概念。这一章中我们将引入PyMC3，并学习如何用它来构建概率模型、对后验进行采样、判断采样是否正确以及分析和解释贝叶斯结果。
第3章，多参和分层模型，介绍了贝叶斯模型中最基础的内容，并在此基础上加入了一些更复杂的内容。我们将学习如何利用多个参数构建并分析模型，以及如何利用分层模型的优势往模型中添加结构。
第4章，利用线性回归模型理解并预测数据，介绍了线性回归模型的广泛应用以及如何将其应用于更复杂的模型。在本章中，我们将学习如何利用线性模型解决回归问题，以及如何处理异常值和多变量的问题。
第5章，利用逻辑回归对结果进行分类，在前一章的基础上对线性模型做了进一步推广，将其应用于解决多输入/多输出分类问题。
第6章，模型比较，讨论了统计和机器学习中一些常见的模型比较难点。我们将学习一些信息测准和贝叶斯因子方面的理论知识，以及如何将其应用于模型比较。
第7章，混合模型，讨论了如何将简单的模型混合在一起构建出更复杂的模型，这种方法将引导我们认识新的模型，并从混合模型的角度重新回顾前面几章中学到的模型。此外，本章还讨论了如何进行数据聚类和如何处理计数类型数据等问题。
第8章，高斯过程，简要讨论了一些非参数统计方面的高级概念作为本书的结束，包括什么是核函数、如何使用线性核回归以及如何将高斯过程用于回归问题。
本书代码部分使用的是Python 3.5以上版本，因此，建议你使用Python 3的最新版，尽管本书的大部分代码都能在更早的版本上运行（包括Python 2.7，不过可能需要稍微修改）。
安装Python和Python库最简单的方法是使用Anaconda（一个用于科学计算的软件），你可以通过网站https://www.continuum.io/downloads了解和下载Anaconda。安装好Anaconda之后，可以使用conda install库的名称来安装Python库。
本书会用到以下Python库：
Ipython 5.0；
NumPy 1.11.1；
SciPy 0.18.1；
Pandas 0.18.1；
Matplotlib 1.5.3；
Seaborn 0.7.1；
PyMC3 3.0。
本书的阅读对象为不熟悉贝叶斯统计方法，同时又希望学习如何进行贝叶斯数据分析的本科生、研究生或数据科学家。本书不要求读者有统计学或贝叶斯分析方面的背景，书中尽可能地减少了数学公式的使用，只在我们认为有利于读者理解相关概念的地方用到。此外，所有的概念都通过代码、图表以及文字进行了详细描述。本书假设读者知道如何使用Python进行编程，最好熟悉NumPy、Matplotlib或者Pandas。
在阅读本书过程中，你会看到一些不同的排版方式用于区分不同的信息，以下是这些排版的例子及其解释。
文字中的代码单词、数据表的名字、文件夹名、文件名、文件扩展名、路径、链接、用户输入都按以下方式排版：“为了准确计算HPD，我们将使用plot_post函数”。
代码片段的排版方式如下：
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)
for i in range(3):
for j in range(3):
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='b', lw=5)
ax[i,j].set_ylim(0, 1)
ax[i,j].plot(0, 0, label="n = {:3.2f}\np = {:3.2f}".format(n, p), alpha=0)
ax[i,j].legend(fontsize=12)
ax[2,1].set_xlabel('$\\theta$', fontsize=14)
ax[1,0].set_ylabel('$p(y|\\theta)$', fontsize=14)
ax[0,0].set_xticks(x)所有命令行的输入或输出都按以下方式排版：
conda install NamePackage
欢迎读者对本书的反馈，让我们知道你关于这本书的想法——喜欢什么，不喜欢什么。读者反馈对于我们很重要，它可以帮助我们开发读者真正需要的话题。想给我们发送反馈，只需要发送电子邮件至feedback@packtpub.com，并在邮件主题中告知书名。如果你是某个话题的专家，并且有兴趣编写书籍或者给予贡献，请查看我们的作者指导： www.packtpub.com/authors 。
你现在已经是Packt书籍的荣誉所有者。你还拥有以下权利。
你可以从http://www.packtpub.com的个人账户下载本书的示例代码文件。如果是从别的地方购买的本书，可以访问http://www.packtpub.com/support，注册后，代码文件会直接通过电子邮件发送给你。你也可以通过下列步骤下载代码文件。
（1）使用你的邮箱和密码在我们的网站登录并注册。
（2）在顶部的SUPPORT标签上悬停光标。
（3）单击Code Downloads & Errata按钮。
（4）在Search文本框中输入书名。
（5）选取代码文件所在的书籍。
（6）选择购书途径的下拉菜单。
（7）单击Code Download按钮。
你也可以在本书网站的页面单击Code Files按钮下载代码文件。可以通过在Search文本框中输入书名找到本书的网页。你需要登录自己的Packt账户。
文件下载完成后，确保使用下列软件的最新版解压或抽取文件：
WinRAR / 7-Zip for Windows；
Zipeg / iZip / UnRarX for Mac；
7-Zip / PeaZip for Linux。
本书涉及的代码也可以在GitHub上下载https://github.com/ PacktPublishing/Bayesian-Analysis-with-Python。此外，在https://github. com/PacktPublishing/中还可以查看一系列其他书籍和视频的代码。
我们还提供了本书中所有彩色图表的PDF文件，从而帮助你理解印刷造成的差异。你可以从https://www.packtpub. com/sites/default/files/downloads/BayesianAnalysiswithPython_ColorImages.pdf下载。
尽管我们已经非常细心地努力保证内容的正确性，但是错误还是不可避免。如果你在本书中发现错误，不管是文本错误或是代码错误，请告诉我们。你的善举会省去其他用户的烦恼，并帮助我们改进本书的后续版本。如果你发现了任何错误，请访问http://www.packtpub.com/submit-errata报告给我们。你只需选取书名，单击Errata Submission Form链接，输入勘误的具体信息。一旦勘误确定之后，我们会接受你的提交。勘误会上传到我们的网站，或者添加到书籍勘误部分已有的勘误列表下。要查看以前提交的勘误，访问https://www.packtpub.com/books/content/support，在搜索框中输入书名，所需信息便会出现在Errata部分下。此外，你也可以在以下地址查看和提交勘误信息：https://github.com/aloctavodia/BAP。
互联网上版权资料的盗版问题一直是所有媒介关心的问题。在Packt，我们一直严肃对待版权和许可的保护问题。如果你在互联网上遇到任何形式的我社出版物的非法副本，请立即把具体地址或者网站名称提供给我们，请联系copyright@packtpub.com，附上可疑的盗版材料的链接。我们非常感谢你在保护作者方面所做的努力，我们会注重提升自我能力，给你带来更有价值的内容。
如果你对本书有任何疑问，可以联系我们：questions@packtpub.com，我们会尽全力解决你的问题。
归根到底，概率论不过是把常识化作计算而已。
——皮埃尔—西蒙•拉普拉斯
本章我们将学习贝叶斯统计中的核心概念以及一些用于贝叶斯分析的基本工具。大部分内容都是一些理论介绍，其中会涉及一些Python代码，绝大多数概念会在本书其余章节中反复提到。尽管本章内容有点偏理论，可能会让习惯代码的你感到有点不安，不过这会让你在后面应用贝叶斯统计方法解决问题时容易一些。
本章包含以下主题：
统计学主要是收集、组织、分析并解释数据，因此，统计学的基础知识对数据分析来说至关重要。分析数据时一个非常有用的技巧是知道如何运用某种编程语言（如Python）编写代码。真实世界里充斥着复杂而杂乱的数据，因此对数据做一些预处理操作必不可少。即便你的数据已经是整理好的了，掌握一定的编程技巧仍然会给你带来很大帮助，因为如今的贝叶斯统计绝大多数都是计算统计学。
大多数统计学导论课程（对非统计学专业的人而言）一般就像展示一本菜谱书，每一种统计方法就是一个菜谱：首先，到统计学的后厨取出一个罐头打开，放点数据上去尝尝，然后不停搅拌直到得出一个稳定的p值，该值最好低于0.05（如果你不知道什么是p值，别担心，本书不会涉及这些概念）。这类课程的目的是教会你如何选择一个合适的罐头。本书采用的是另外一种方式：首先我们也需要点原料，不过这次是自己亲自做的而不是买来的罐头，然后学习如何把新鲜的食材混合在一起以适应不同的烹饪场景。在正式烹饪之前，我们先学点统计学的术语和概念。
数据是统计学最基本的组成部分。数据的来源多，比如实验、计算机模拟、调查以及观测等。假如我们是数据生成或收集人员，首先要考虑的是要解决什么样的问题以及打算采用什么方法，然后再去着手准备数据。事实上，统计学有一个叫做实验设计的分支专门研究如何获取数据。在这个数据泛滥的年代，我们有时候会忘了获取数据并非总是很便宜。例如，尽管大型强子碰撞加速装置一天能产生上百TB的数据，但其建造却要花费数年的人力和智力。本书假设我们已经获取了数据并且数据是整理好的（这在现实中通常很少见），以便关注到本书的主题上来。如果你想学习如何用Python做数据清洗和分析并进一步学习机器学习，你可以阅读Jake VanderPlas写的《Python Data Science Handbook》一书。
假设我们已经有了数据集，通常的做法是先对其探索并可视化，这样我们就能对手头的数据有个直观的认识。可以通过如下两步完成所谓的探索式数据分析过程：
其中，描述性统计是指如何用一些指标或统计值来定量地总结或刻画数据，例如你已经知道了如何用均值、众数、标准差、四分位差等指标来描述数据。数 据可视化是指用生动形象的方式表述数据，你大概对直方图、散点图等表现形式比较熟悉。乍看起来，探索式数据分析似乎是在复杂分析之前的一些准备工作，或者是作为一些复杂分析方法的替代品，不过探索式数据分析在理解、解释、检查、总结及交流贝叶斯分析结果等过程中依然有用。
有时候，画画图、对数据做些简单的计算（比如求均值）就够了。另外一些时候，我们希望从数据中挖掘出一些更一般性的结论。我们可能希望了解数据是怎么生成的，也可能是想对未来还未观测到的数据做出预测，又或者是希望从多个对观测值的解释中找出最合理的一个，这些正是统计推断所做的事情。模型分为许多种，统计推断依赖的是概率模型，许多科学研究（以及我们对真实世界的认识）也都是基于模型的， 大脑不过是对现实进行建模的一台机器，可以观看相关的TED演讲了解大脑是如何对现实进行建模的，网址为http://www.tedxriodelaplata.org/videos/m%C3%A1quina-construye-realidad。
什么是模型？模型是对给定系统或过程的一种简化描述。这些描述只关注系统中某些重要的部分，因此，大多数模型的目的并不是解释整个系统。此外，假如我们有两个模型能用来解释同一份数据并且效果差不多，其中一个简单点，另一个复杂一些，通常我们倾向于更简单的模型，这称作奥卡姆剃刀，我们会在第6章模型比较部分讨论贝叶斯分析与其之间的联系。
不管你打算构建哪种模型，模型构建都遵循一些相似的基本准则，我们把贝叶斯模型的构建过程总结为如下3步。
（1）给定一些数据以及这些数据是如何生成的假设，然后构建模型。通常，这里的模型都是一些很粗略的近似，不过大多时候也够用了。
（2）利用贝叶斯理论将数据和模型结合起来，根据数据和假设推导出逻辑结论，我们称之为经数据拟合后的模型。
（3）根据多种标准，包括真实数据和对研究问题的先验知识，判断模型拟合得是否合理。
通常，我们会发现实际的建模过程并非严格按照该顺序进行的，有时候我们有可能跳到其中任何一步，原因可能是编写的程序出错了，也可能是找到了某种改进模型的方式，又或者是我们需要增加更多的数据。
贝叶斯模型是基于概率构建的，因此也称作概率模型。为什么基于概率呢？因为概率这个数学工具能够很好地描述数据中的不确定性，接下来我们将对其进行深入了解。
尽管概率论是数学中一个相当成熟和完善的分支，但关于概率的诠释仍然有不止一种。对于贝叶斯派而言，概率是对某一命题不确定性的衡量。假设我们对硬币一无所知，同时没有与抛硬币相关的任何数据，那么可以认为正面朝上的概率介于0到1之间，也就是说，在缺少信息的情况下，所有情况都是有可能发生的，此时不确定性也最大。假设现在我们知道硬币是公平的，那么我们可以认为正面朝上的概率是0.5或者是0.5附近的某个值（假如硬币不是绝对公平的话），如果此时收集数据，我们可以根据观测值进一步更新前面的先验假设，从而降低对该硬币偏差的不确定性。按照这种定义，提出以下问题都是自然而且合理的：火星上有多大可能存在生命？电子的质量是9.1×10−31kg的概率是多大？1816年7月9号是晴天的概率是多少？值得注意的是，类似火星上是否有生命这种问题的答案是二值化的，但我们关心的是，基于现有数据以及我们对火星物理条件和生物条件的了解，火星上存在生命的概率有多大。该命题取决于我们当前所掌握的信息而非客观的自然属性。我们使用概率是因为我们对事件不确定，而不代表事件本身是不确定的。由于这种概率的定义取决于我们的认知水平，有时也被称为概率的主观定义，这也就解释了为什么贝叶斯派总被称作主观统计。然而，这种定义并不是说所有命题都是同等有意义的，仅是承认我们对世界的理解是基于现有的数据和模型，是不完备的。不基于模型或理论去理解世界是不可能的。即使我们能脱离社会预设，最终也会受到生物上的限制：受进化过程影响，我们的大脑已经与世界上的模型建立了联系。我们注定要以人类的方式去思考，永远不可能像蝙蝠或其他动物那样思考。而且，宇宙是不确定的，总的来说我们能做的就是对其进行概率性描述。需要注意的是，不管世界的本原是确定的还是随机的，我们都将概率作为衡量不确定性的工具。
逻辑是指如何有效地推论，在亚里士多德学派或者经典的逻辑学中，一个命题只能是对或者错，而在贝叶斯学派定义的概率中，确定性只是概率上的一种特殊情况：一个正确命题的概率值为1，而一个错误命题的概率值为0。只有在我们拥有充分的数据表明火星上存在生物生长和繁殖时，我们才认为“火星上存在生命”这一命题为真的概率值为1。不过，需要注意的是，认定一个命题的概率值为0通常是比较困难的，这是因为火星上可能存在某些地方还没被探测到，或者是我们的实验有问题，又或者是某些其他原因导致我们错误地认为火星上没有生命而实际上有。与此相关的是克伦威尔准则（Cromwell’s Rule），其含义是指在对逻辑上正确或错误的命题赋予概率值时，应当避免使用0或者1。有意思的是，考克斯（Cox）在数学上证明了如果想在逻辑推理中加入不确定性，我们就必须使用概率论的知识，由此来说，贝叶斯定理可以视为概率论逻辑上的结果。因此，从另一个角度来看，贝叶斯统计是对逻辑学处理不确定性问题的一种扩充，当然这里毫无对主观推理的轻蔑。现在我们已经熟悉了贝叶斯学派对概率的理解，接下来就了解下概率相关的数学特性吧。如果你想深入学习概率论，可以参考阅读Joseph K Blitzstein和Jessica Hwang写的《概率导论》（Introduction to Probability）。
概率值介于[0,1]之间（包括0和1），其计算遵循一些法则，其中之一是乘法法则：
上式中，A和B同时发生的概率值等于B发生的概率值乘以在B发生的条件下A也发生的概率值，其中，
表示A和B的联合概率， 
表示条件概率，二者的现实意义是不同的，例如，路面是湿的概率跟下雨时候路面是湿的概率是不同的。条件概率可能比原来的概率高，也可能低。如果B并不能提供任何关于A的信息，那么，
，也就是说，A和B是相互独立的。相反，如果事件B能够给出关于事件A的一些信息，那么根据事件B提供的信息不同，事件A可能发生的概率会变得更高或是更低。
条件概率是统计学中的一个核心概念，接下来我们将看到，理解条件概率对于理解贝叶斯理论至关重要。这里我们换个角度来看条件概率，假如我们把前面的公式调整下顺序，就可以得到下面的公式：
需要注意的是，我们不对0概率事件计算条件概率，因此，分母的取值范围是(0,1)，从而可以看出条件概率大于或等于联合概率。为什么要除以p(B)呢？因为在已经知道事件B发生的条件下，我们考虑的可能性空间就缩小到了事件B发生的范围内，然后将该范围内A发生的可能性除以B发生的可能性便得到了条件概率
。需要强调的是，所有的概率本质上都是条件概率，世间并没有绝对的概率。不管我们是否留意，在谈到概率时总是存在这样或那样的模型、假设或条件。比如，当我们讨论明天下雨的概率时，在地球上、火星上甚至宇宙中其他地方明天下雨的概率是不同的。同样，硬币正面朝上的概率取决于我们对硬币有偏性的假设。在理解概率的含义之后，接下来我们讨论另一个话题：概率分布。
概率分布是数学中的一个概念，用来描述不同事件发生的可能性，通常这些事件限定在一个集合内，代表了所有可能发生的事件。在统计学里可以这么理解：数据是从某种参数未知的概率分布中生成的。由于并不知道具体的参数，我们只能借用贝叶斯定理从仅有的数据中反推参数。概率分布是构建贝叶斯模型的基础，不同分布组合在一起之后可以得到一些很有用的复杂模型。
本书会介绍一些概率分布，在第一次介绍某个概率分布时，我们会先花点时间理解它。最常见的一种概率分布是高斯分布，又叫正态分布，其数学公式描述如下：
上式中，μ和σ是高斯分布的两个参数。第1个参数μ是该分布的均值（同时也是中位数和众数），其取值范围是任意实数，即
；第2个参数σ是标准差，用来衡量分布的离散程度，其取值只能为正。由于μ和σ的取值范围无穷大，因此高斯分布的实例也有无穷多。虽然数学公式这一表达形式简洁明了，也有人称之有美感，不过得承认公式还是有些不够直观，我们可以尝试用Python代码将公式的含义重新表示出来。首先看看高斯分布都长什么样：
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
import seaborn as sns
mu_params = [-1, 0, 1]
sd_params = [0.5, 1, 1.5]
x = np.linspace(-7, 7, 100)
f, ax = plt.subplots(len(mu_params), len(sd_params), sharex=True, sharey=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(0, 0,
label="$\\mu$ = {:3.2f}\n$\\sigma$ = {:3.2f}".format (mu, sd), alpha=0)
ax[i,j].legend(fontsize=12)
ax[2,1].set_xlabel('$x$', fontsize=16)
ax[1,0].set_ylabel('$pdf(x)$', fontsize=16)
plt.tight_layout()上面代码的输出结果如下：
由概率分布生成的变量（例如x）称作随机变量，当然这并不是说该变量可以取任意值，相反，我们观测到该变量的数值受到概率分布的约束，而其随机性源于我们只知道变量的分布却无法准确预测该变量的值。通常，如果一个随机变量服从在参数μ和σ下的高斯分布，我们可以这样表示该变量：
其中，符号～读作服从于某种分布。
随机变量分为两种：连续变量和离散变量。连续随机变量可以从某个区间内取任意值（我们可以用Python中的浮点型数据来表示），而离散随机变量只能取某些特定的值（我们可以用Python中的整型数据来表示）。
许多模型都假设，如果对服从于同一个分布的多个随机变量进行连续采样，那么各个变量的采样值之间相互独立，我们称这些随机变量是独立同分布的。用数学语言描述就是，如果两个随机变量x和y对于所有可能的取值都满足
，那么称这两个变量相互独立。
时间序列是不满足独立同分布的一个典型例子。在时间序列中，需要对时间维度的变量多加留心。下面的例子是从http://cdiac.esd.ornl.gov中获取的数据。这份数据记录了从1959年到1997年大气中二氧化碳的含量。我们用以下代码将它写出来：
data = np.genfromtxt('mauna_loa_CO2.csv', delimiter=',')
plt.plot(data[:,0], data[:,1])
plt.xlabel('$year$', fontsize=16)
plt.ylabel('$CO_2 (ppmv)$', fontsize=16)图中每个点表示每个月空气中二氧化碳含量的测量值，可以看到测量值是与时间相关的。我们可以观察到两个趋势：一个是季节性的波动趋势（这与植物周期性生长和衰败有关）；另一个是二氧化碳含量整体性的上升趋势。
到目前为止，我们已经学习了一些统计学中的基本概念和词汇，接下来让我们首先看看神奇的贝叶斯定理：
看起来稀松平常，似乎跟小学课本里的公式差不多，不过这就是关于贝叶斯统计你所需要掌握的全部。首先看看贝叶斯定理是怎么来的，这对我们理解它会很有帮助。事实上，我们已经掌握了如何推导它所需要的全部概率论知识。
现在，让我们看看这个式子的含义及其重要性。首先，上式表明
和
并不一定相等，这一点非常重要，日常分析中即使系统学习过统计学和概率论的人也很容易忽略这点。我们举个简单例子来说明为什么二者不一定相等：有两条腿的动物就是人的概率和人有两条腿的概率显然是不同的。几乎所有人都有两条腿（除了某些人因为先天性原因或者意外导致没有两条腿），但是有两条腿的动物中很多都不是人类，比如鸟类。
在前面的式子中，如果我们将H理解为假设，D理解为数据，那么贝叶斯定理告诉我们的就是，在给定数据的条件下如何计算假设成立的概率。不过，如何把假设融入贝叶斯定理中去呢？答案是概率分布。换句话说，H是一种狭义上的假设，我们所做的实际上是寻找模型的参数（更准确地说是参数的分布）。因此，与其称H为假设，不如称之为模型，这样能避免歧义。
贝叶斯定理非常重要，后面会反复用到，这里我们先熟悉下其各个部分的名称：
p(H )：先验；
p(D | H )：似然；
p(H | D)：后验；
p(D)：证据。
先验分布反映的是在观测到数据之前我们对参数的了解，如果我们对参数一无所知（就跟《权力的游戏》中的雪诺一样），那么可以用一个不包含太多信息的均匀分布来表示。由于引入了先验，有些人会认为贝叶斯统计是偏主观的，然而，这些先验不过是构建模型时的一些假设罢了，其主观性跟似然差不多。
似然是指如何在实验分析中引入观测数据，反映的是在给定参数下得到某组观测数据的可信度。
后验分布是贝叶斯分析的结果，反映的是在给定数据和模型的条件下我们对问题的全部认知。需要注意，后验指的是我们模型中参数的概率分布而不是某个值，该分布正比于先验乘以似然。有这么个笑话：贝叶斯学派就像是这样一类人，心里隐约期待着一匹马，偶然间瞥见了一头驴，结果坚信他看到的是一头骡子。当然了，如果要刻意纠正这个笑话的话，在先验和似然都比较含糊的情况下，我们会得到一个（模糊的）“骡子”后验。不过，这个笑话也讲出了这样一个道理，后验其实是对先验和似然的某种折中。从概念上讲，后验可以看做是在观测到数据之后对先验的更新。事实上，一次分析中的后验，在收集到新的数据之后，也可以看做是下一次分析中的先验。这使得贝叶斯分析特别适合于序列化的数据分析，比如通过实时处理来自气象站和卫星的数据从而提前预警灾害，更详细的内容可以阅读在线机器学习方面的算法。
最后一个概念是证据，也称作边缘似然。正式地讲，证据是在模型的参数取遍所有可能值的条件下得到指定观测值的概率的平均。不过，本书的大部分内容并不关心这个概念，我们可以简单地把它当作归一化系数。我们只关心参数的相对值而非绝对值。把证据这一项忽略掉之后，贝叶斯定理可以表示成如下正比例形式：
理解其中的每个概念可能需要时间和更多的例子，本书也将围绕这些内容展开。
前面，我们学习了几个重要概念，其中有两个是贝叶斯统计的核心概念，这里我们用一句话再重新强调下：概率是用来衡量参数不确定性的，贝叶斯定理就是用来在观测到新的数据时正确更新这些概率以期降低我们的不确定性。
现在我们已经知道什么是贝叶斯统计了，接下来就从一个简单的例子入手，通过推断单个未知参数来学习如何进行贝叶斯统计。
抛硬币是统计学中的一个经典问题，其描述如下：我们随机抛一枚硬币，重复一定次数，记录正面朝上和反面朝上的次数，根据这些数据，我们需要回答诸如这枚硬币是否公平，以及更进一步这枚硬币有多不公平等问题。抛硬币是一个学习贝叶斯统计非常好的例子，一方面是因为几乎人人都熟悉抛硬币这一过程，另一方面是因为这个模型很简单，我们可以很容易计算并解决这个问题。此外，许多真实问题都包含两个互斥的结果，例如0或者1、正或者负、奇数或者偶数、垃圾邮件或者正常邮件、安全或者不安全、健康或者不健康等。因此，即便我们讨论的是硬币，该模型也同样适用于前面这些问题。
为了估计硬币的偏差，或者更广泛地说，想要用贝叶斯学派理论解决问题，我们需要数据和一个概率模型。对于抛硬币这个问题，假设我们已试验了一定次数并且记录了正面朝上的次数，也就是说数据部分已经准备好了，剩下的就是模型部分了。考虑到这是第一个模型，我们会列出所有必要的数学公式，并且一步一步推导。下一章中，我们会重新回顾这个问题，并借用PyMC3从数值上解决它（也就是说那部分不需要手动推导，而是利用PyMC3和计算机来完成）。
首先，我们要抽象出偏差的概念。我们称，如果一枚硬币总是正面朝上，那么它的偏差就是1，反之，如果总是反面朝上，那么它的偏差就是0，如果正面朝上和反面朝上的次数各占一半，那么它的偏差就是0.5。这里用参数θ来表示偏差，用y表示N次抛硬币实验中正面朝上的次数。根据贝叶斯定理，我们有如下公式：
这里需要指定我们将要使用的先验
和似然
分别是什么。让我们首先从似然开始。
假设多次抛硬币的结果相互之间没有影响，也就是说每次抛硬币都是相互独立的，同时还假设结果只有两种可能：正面朝上或者反面朝上。但愿你能认同我们对这个问题做出的合理假设。基于这些假设，一个不错的似然候选是二项分布：
这是一个离散分布，表示N次抛硬币实验中y次正面朝上的概率（或者更通俗地描述是，N次实验中，y次成功的概率）。下面的代码生成了9个二项分布，每个子图中的标签显示了对应的参数：
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)
for i in range(3):
for j in range(3):
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='b', lw=5)
ax[i,j].set_ylim(0, 1)
ax[i,j].plot(0, 0, label="n = {:3.2f}\np =
{:3.2f}".format(n, p), alpha=0)
ax[i,j].legend(fontsize=12)
ax[2,1].set_xlabel('$\\theta$', fontsize=14)
ax[1,0].set_ylabel('$p(y|\\theta)$', fontsize=14)
ax[0,0].set_xticks(x)二项分布是似然的一个合理选择，直观上讲，θ可以看作抛一次硬币时正面朝上的可能性，并且该过程发生了y次。类似地，我们可以把“1−θ”看作抛一次硬币时反面朝上的概率，并且该过程发生了“N−y”次。
假如我们知道了θ，那么就可以从二项分布得出硬币正面朝上的分布。如果我们不知道θ，也别灰心，在贝叶斯统计中，当我们不知道某个参数的时候，就对其赋予一个先验。接下来继续选择先验。
这里我们选用贝叶斯统计中最常见的beta分布，作为先验，其数学形式如下：
仔细观察上面的式子可以看出，除了Γ部分之外，beta分布和二项分布看起来很像。Γ是希腊字母中大写的伽马，用来表示伽马函数。现在我们只需要知道，用分数表示的第一项是一个正则化常量，用来保证该分布的积分为1，此外，
和
两个参数用来控制具体的分布形态。beta分布是我们到目前为止见到的第3个分布，利用下面的代码，我们可以深入了解其形态：
params = [0.5, 1, 2, 3]
x = np.linspace(0, 1, 100)
f, ax = plt.subplots(len(params), len(params), sharex=True, sharey=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="$\\alpha$ = {:3.2f}\n$\\beta$ = {:3.2f}".format(a, b), alpha=0)
ax[i,j].legend(fontsize=12)
ax[3,0].set_xlabel('$\\theta$', fontsize=14)
ax[0,0].set_ylabel('$p(\\theta)$', fontsize=14)
plt.savefig('B04958_01_04.png', dpi=300, figsize=(5.5, 5.5))为什么要在模型中使用beta分布呢？在抛硬币以及一些其他问题中使用beta分布的原因之一是：beta分布的范围限制在0到1之间，这跟我们的参数
一样；另一个原因是其通用性，从前面的图可以看出，该分布可以有多种形状，包括均匀分布、类高斯分布、U型分布等。第3个原因是：beta分布是二项分布（前面我们使用了该分布描述似然）的共轭先验。似然的共轭先验是指，将该先验分布与似然组合在一起之后，得到的后验分布与先验分布的表达式形式仍然是一样的。简单说，就是每次用beta分布作为先验、二项分布作为似然时，我们会得到一个beta分布的后验。除beta分布之外还有许多其他共轭先验，例如高斯分布，其共轭先验就是自己。关于共轭先验更详细的内容可以查看https://en.wikipedia.org/wiki/Conjugate_prior。许多年来，贝叶斯分析都限制在共轭先验范围内，这主要是因为共轭能让后验在数学上变得更容易处理，要知道贝叶斯统计中一个常见问题的后验都很难从分析的角度去解决。在建立合适的计算方法来解决任意后验之前，这只是个折中的办法。从下一章开始，我们将学习如何使用现代的计算方法来解决贝叶斯问题而不必考虑是否使用共轭先验。
首先回忆一下贝叶斯定理：后验正比于似然乘以先验。
对于我们的问题而言，需要将二项分布乘以beta分布：
现在，对上式进行简化。针对我们的实际问题，可以先把与θ不相关的项去掉而不影响结果，于是得到下式：
重新整理之后得到：
可以看出，上式和beta分布的形式很像（除了归一化部分），其对应的参数分别为
、
。也就是说，在抛硬币这个问题中，后验分布是如下beta分布：
计算后验并画图 现在已经有了后验的表达式，我们可以用Python对其计算并画出结果。下面的代码中，其实只有一行是用来计算后验结果的，其余的代码都是用来画图的：
theta_real = 0.35
trials = [0, 1, 2, 3, 4, 8, 16, 32, 50, 150]
data = [0, 1, 1, 1, 1, 4, 6, 9, 13, 48]
beta_params = [(1, 1), (0.5, 0.5), (20, 20)]
dist = stats.beta
x = np.linspace(0, 1, 100)
for idx, N in enumerate(trials):
if idx == 0:
plt.subplot(4,3, 2)
else:
plt.subplot(4,3, idx+3)
y = data[idx]
for (a_prior, b_prior), c in zip(beta_params, ('b', 'r', 'g')):
p_theta_given_y = dist.pdf(x, a_prior + y, b_prior + N - y)
plt.plot(x, p_theta_given_y, c)
plt.fill_between(x, 0, p_theta_given_y, color=c, alpha=0.6)
plt.axvline(theta_real, ymax=0.3, color='k')
plt.plot(0, 0, label="{:d} experiments\n{:d} heads".format(N, y), alpha=0)
plt.xlim(0,1)
plt.ylim(0,12)
plt.xlabel(r"$\theta$")
plt.legend()
plt.gca().axes.get_yaxis().set_visible(False)
plt.tight_layout()
plt.savefig('B04958_01_05.png', dpi=300, figsize=(5.5, 5.5))在上图的第一行中，实验的次数为0，因此第一个图中的曲线描绘的是先验分布，其中有3条曲线，每条曲线分别表示一种先验。
剩余的子图描绘了后续实验的后验分布
，回想一下，后验可以看做是在给定数据之后更新了的先验。实验（抛硬币）的次数和正面朝上的次数分别标注在每个子图中。此外每个子图中在横轴0.35附近还有一个黑色的竖线，表示的是真实的θ，显然，在真实情况下，我们并不知道该值，在这里标识出来只是为了方便理解。从这幅图中可以学到很多贝叶斯分析方面的知识。
从前面的例子可以看出，先验对分析的结果会有影响。一些贝叶斯分析的新手（以及一些诋毁该方法的人）会对如何选择先验感到茫然，因为他们不希望先验起到决定性作用，而是更希望数据本身替自己说话！有这样的想法很正常，不过我们得牢记，数据并不会真的“说话”，只有在模型中才会有意义，包括数学上的和脑海中的模型。面对同一主题下的同一份数据，不同人会有不同的看法，这类例子在科学史上有许多，查看以下链接可以了解《纽约时报》最近一次实验的例子：http://www.nytimes.com/interactive/2016/09/20/upshot/the-error-the-polling-world-rarely-talks-about.html?_r=0。
有些人青睐于使用没有信息量的先验（也称作均匀的、含糊的或者发散的先验），这类先验对分析过程的影响最小。本书将遵循Gelman、McElreath和Kruschke 3人的建议[1]，更倾向于使用带有较弱信息量的先验。在许多问题中，我们对参数可以取的值一般都会有些了解，比如，参数只能是正数，或者知道参数近似的取值范围，又或者是希望该值接近0或大于/小于某值。这种情况下，我们可以给模型加入一些微弱的先验信息而不必担心该先验会掩盖数据本身的信息。由于这类先验会让后验近似位于某一合理的边界内，因此也被称作正则化先验。
当然，使用带有较多信息量的强先验也是可行的。视具体的问题不同，有可能很容易或者很难找到这类先验，例如在我工作的领域（结构生物信息学），人们会尽可能地利用先验信息，通过贝叶斯或者非贝叶斯的方式来了解和预测蛋白质的结构。这样做是合理的，原因是我们在数十年间已经从上千次精心设计的实验中收集了数据，因而有大量可信的先验信息可供使用。如果你有可信的先验信息，完全没有理由不去使用。试想一下，如果一个汽车工程师每次设计新车的时候，他都要重新发明内燃机、轮子乃至整个汽车，这显然不是正确的方式。
现在我们知道了先验有许多种，不过这并不能缓解我们选择先验时的焦虑。或许，最好是没有先验，这样事情就简单了。不过，不论是否基于贝叶斯，模型都在某种程度上拥有先验，即使这里的先验并没有明确表示出来。事实上，许多频率统计学方面的结果可以看做是贝叶斯模型在一定条件下的特例，比如均匀先验。让我们再仔细看看前面那幅图，可以看到蓝色后验分布的峰值与频率学分析中θ的期望值是一致的：
注意，这里
是点估计而不是后验分布（或者其他类型的分布）。由此看出，你没办法完全避免先验，不过如果你在分析中引入先验，得到的会是概率分布分布而不只是最可能的一个值。明确引入先验的另一个好处是，我们会得到更透明的模型，这意味着更容易评判、（广义上的）调试以及优化。构建模型是一个循序渐进的过程，有时候可能只需要几分钟，有时候则可能需要数年；有时候整个过程可能只有你自己，有时候则可能包含你不认识的人。而且，模型复现很重要，而模型中透明的假设能有助于复现。
在特定分析任务中，如果我们对某个先验或者似然不确定，可以自由使用多个先验或者似然进行尝试。模型构建过程中的一个环节就是质疑假设，而先验就是质疑的对象之一。不同的假设会得到不同的模型，根据数据和与问题相关的领域知识，我们可以对这些模型进行比较，本书第6章模型比较部分会深入讨论该内容。
由于先验是贝叶斯统计中的一个核心内容，在接下来遇到新的问题时我们还会反复讨论它，因此，如果你对前面的讨论内容感到有些疑惑，别太担心，要知道人们在这个问题上已经困惑了数十年并且相关的讨论一直在继续。
现在我们已经有了后验，相关的分析也就结束了。下面我们可能还需要对分析结果进行总结，将分析结果与别人分享，或者记录下来以备日后使用。
根据受众不同，你可能在交流分析结果的同时还需要交流模型。以下是一种简单表示概率模型的常见方式：
这是我们抛硬币例子中用到的模型。符号～表示左边随机变量的分布服从右边的分布形式，也就是说，这里θ服从于参数为α和β的Beta分布，而y服从于参数为n = 1和p = θ的二项分布。该模型还可以用类似Kruschke书中的图表示成如下形式：
在第一层，根据先验生成了θ，然后通过似然生成最下面的数据。图中的箭头表示变量之间的依赖关系，符号～表示变量的随机性。
本书中用到的类似Kruschke中的图都是由Rasmus Bååth（http://www.sumsar.net/blog/2013/10/diy-kruschke-style-diagrams/）提供的模板生成的。
贝叶斯分析的结果是后验分布，该分布包含了有关参数在给定数据和模型下的所有信息。如果可能的话，我们只需要将后验分布展示给观众即可。通常，一个不错的做法是：同时给出后验分布的均值（或者众数、中位数），这样能让我们了解该分布的中心，此外还可以给出一些描述该分布的衡量指标，如标准差，这样人们能对我们估计的离散程度和不确定性有一个大致的了解。拿标准差衡量类似正态分布的后验分布很合适，不过对于一些其他类型的分布（如偏态分布）却可能得出误导性结论，因此，我们还可以采用下面的方式衡量。
一个经常用来描述后验分布分散程度的概念是最大后验密度（ Highest Posterior Density，HPD）区间。一个HPD区间是指包含一定比例概率密度的最小区间，最常见的比例是95%HPD或98%HPD，通常还伴随着一个50%HPD。如果我们说某个分析的HPD区间是[2, 5]，其含义是指：根据我们的模型和数据，参数位于2～5的概率是0.95。这是一个非常直观的解释，以至于人们经常会将频率学中的置信区间与贝叶斯方法中的可信区间弄混淆。如果你对频率学的范式比较熟悉，请注意这两种区间的区别。贝叶斯学派的分析告诉我们的是参数取值的概率，这在频率学的框架中是不可能的，因为频率学中的参数是固定值，频率学中的置信区间只能包含或不包含参数的真实值。在继续深入之前，有一点需要注意：选择95%还是50%或者其他什么值作为HPD区间的概率密度比例并没有什么特殊的地方，这些不过是经常使用的值罢了。比如，我们完全可以选用比例为91.37%的HPD区间。如果你选的是95%，这完全没问题，只是要记住这只是个默认值，究竟选择多大比例仍然需要具体问题具体分析。
对单峰分布计算95%HPD很简单，只需要计算出2.5%和97.5%处的值即可：
def naive_hpd(post):
sns.kdeplot(post)
HPD = np.percentile(post, [2.5, 97.5])
plt.plot(HPD, [0, 0], label='HPD {:.2f} {:.2f}'.format(*HPD),
linewidth=8, color='k')
plt.legend(fontsize=16);
plt.xlabel(r'$\theta$', fontsize=14)
plt.gca().axes.get_yaxis().set_ticks([])
np.random.seed(1)
post = stats.beta.rvs(5, 11, size=1000)
naive_hpd(post)
plt.xlim(0, 1)对于多峰分布而言，计算HPD要稍微复杂些。如果把对HPD的原始定义应用到混合高斯分布上，我们可以得到：
np.random.seed(1)
gauss_a = stats.norm.rvs(loc=4, scale=0.9, size=3000)
gauss_b = stats.norm.rvs(loc=-2, scale=1, size=2000)
mix_norm = np.concatenate((gauss_a, gauss_b))
naive_hpd(mix_norm)
plt.savefig('B04958_01_08.png', dpi=300, figsize=(5.5, 5.5))从上图可以看出，通过原始HPD定义计算出的可信区间包含了一部分概率较低的区间，位于[0,2]。为了正确计算出HPD，这里我们使用了plot_post函数，你可以从本书附带的代码中下载对应的源码：
from plot_post import plot_post
plot_post(mix_norm, roundto=2, alpha=0.05)
plt.legend(loc=0, fontsize=16)
plt.xlabel(r"$\theta$", fontsize=14)从上图可以看出，95%HPD包含两个区间，同时plot_post函数也返回了两个众数。
贝叶斯方法的一个优势是：一旦得到了后验分布，我们可以根据该后验生成未来的数据y，即用来做预测。后验预测检查主要是对观测数据和预测数据进行比较从而发现这两个集合的不同之处，其目的是进行一致性检查。生成的数据和观测的数据应该看起来差不多，否则有可能是建模出现了问题或者输入数据到模型时出了问题，不过就算我们没有出错，两个集合仍然有可能出现不同。尝试去理解其中的偏差有助于我们改进模型，或者至少能知道我们模型的极限。即使我们并不知道如何去改进模型，但是理解模型捕捉到了问题或数据的哪些方面以及没能捕捉到哪些方面也是非常有用的信息。也许模型能够很好地捕捉到数据中的均值但却没法预测出罕见值，这可能是个问题，不过如果我们只关心均值，这个模型对我们而言也还是可用的。通常我们的目的不是去声称一个模型是错误的，我们更愿意遵循George Box的建议，即所有模型都是错的，但某些是有用的。我们只想知道模型的哪个部分是值得信任的，并测试该模型是否在特定方面符合我们的预期。不同学科对模型的信任程度显然是不同的，物理学中研究的系统是在高可控条件下依据高深理论运行的，因而模型可以看做是对现实的不错描述，而在一些其他学科如社会学和生物学中，研究的是错综复杂的孤立系统，因而模型对系统的认知较弱。尽管如此，不论你研究的是哪一门学科，都需要对模型进行检查，利用后验预测和本章学到的探索式数据分析中的方法去检查模型。
本书用到的代码是用Python 3.5写的，建议使用Python 3的最新版本运行，尽管大多数代码也能在更老的Python版本（包括Python 2.7）上运行，不过可能会需要做些微调。
本书建议使用Anaconda安装Python及相关库，Anaconda是一个用于科学计算的软件分发，你可以从以下链接下载并了解更多：https://www.continuum.io/downloads。在系统上装好Anaconda之后，就可以通过以下方式安装Python库了：
conda install NamePackage
我们会用到以下Python库：
在命令行中执行以下命令即可安装最新稳定版的PyMC3：
pip install pymc3 我们的贝叶斯之旅首先围绕统计建模、概率论和贝叶斯定理做了一些简短讨论，然后用抛硬币的例子介绍了贝叶斯建模和数据分析，借用这个经典例子传达了贝叶斯统计中的一些最重要的思想，比如用概率分布构建模型并用概率分布来表示不确定性。此外我们尝试揭示了如何选择先验，并将其与数据分析中的一些其他问题置于同等地位（怎么选择似然，为什么要解决该问题等）。本章的最后讨论了如何解释和报告贝叶斯分析的结果。本章我们对贝叶斯分析的一些主要方面做了简要总结，后面还会重新回顾这些内容，从而充分理解和吸收，并为后面理解更高级的内容打下基础。下一章我们会重点关注一些构建和分析更复杂模型的技巧，此外，还会介绍PyMC3并将其用于实现和分析贝叶斯模型。
我们尚不清楚大脑是如何运作的，是按照贝叶斯方式？还是类似贝叶斯的某种方式？又或许是进化过程中形成的某种启发式的方式？不管如何，我们至少知道自己是通过数据、例子和练习来学习的，尽管你可能对此有不同的意见，不过我仍然强烈建议你完成以下练习。
（1）修改生成本章的第3个图的代码，在图中增加一条竖线用来表示观测到的正面朝上的比例（正面朝上的次数/抛硬币的次数），将该竖线的位置与每个子图中后验的众数进行比较。
（2）尝试用不同的先验参数（beta值）和不同的实验数据（正面朝上的次数和实验次数）重新绘制本章的第3个图。
（3）阅读维基百科上有关Cromwell准则的内容：https://en.wikipedia.org/wiki/Cromwell%27s_rule。
（4）探索不同参数下高斯分布、二项分布和beta分布的图像，你可以为每个分布单独画一个图而不是全都画在一个网格中。
[1] 该3人分别是《Bayesian Data Analysis》《Statistical Rethinking: A Bayesian Course with Examples in R and Stan》和《Doing Bayesian Data Analysis》的主要作者。——译者注
前面我们已经对贝叶斯统计有了初步了解，接下来将学习如何利用一些计算工具构建概率模型，还将学习概率编程。主要目的是利用代码描述模型并用其进行推断，当然这并非因为我们太懒了所以才没有采用数学的方法去推断，也不是因为我们是最棒的极客（做梦都在写代码……）。该决定主要基于这样一个原因：许多模型都没有一个封闭的解析后验，也就是说我们只能使用数值的方法计算后验。学习概率编程的另一个原因是：现代的贝叶斯统计主要是通过编程实现的，概率编程提供了一个有效的方式去构建复杂模型，它让我们更加关注模型设计、评估和解释，而不用过多地考虑数学或计算细节。
本章我们将学习贝叶斯统计中的数值计算方法以及如何使用PyMC3。PyMC3是概率编程中一个非常灵活的库，了解PyMC3有助于我们从更实用的角度学习更高级的贝叶斯概念。
本章涵盖以下主题：
贝叶斯统计的概念很简单，我们有一些固定的数据（固定的意思是指我们无法改变观测值），和一些感兴趣的参数，剩下要做的就是探索这些参数可能的取值，其中所有的不确定性都通过概率进行建模。在别的统计学范式中，未知量有多种不同的表示方式，而在贝叶斯的框架中，所有未知量都是同等对待的，如果不知道某个变量，我们就给其赋予一个概率分布。因此，贝叶斯定理就是将先验概率分布（在观测到数据之前我们对问题的理解）转化成后验分布
（观测到数据之后所得到的信息），换句话说，贝叶斯统计就是一种机器学习的过程。
尽管概念上很简单，纯粹的概率模型得到的表达式通常分析起来很棘手。许多年来，这都是一个问题，大概也是影响贝叶斯方法广泛应用的最大原因之一。随着计算时代的到来，数值化方法的发展使得计算几乎任意模型的后验成为了可能，这极大地促进了贝叶斯方法的应用。我们可以把这些数值化方法当作通用引擎，按照PyMC3的核心开发者之一Thomas Wiecki的说法，只需要按下按钮，推断部分就可以自动完成了。
自动化推断促进了概率编程语言的发展，从而使得模型构建和推断相分离。在概率编程语言的框架中，用户只需要寥寥数行代码描述概率模型，后面的推断过程就能自动完成了。概率编程使得人们能够更快速地构建复杂的概率模型并减少出错的可能，可以预见，这将给数据科学和其他学科带来极大的影响。
我认为，编程语言对科学计算的影响可以与60多年前Fortran语言的问世相对比。尽管如今Fortran语言风光不再，不过在当时Fortran语言被认为是相当革命性的。科学家们第一次从计算的细节中解放出来，开始用一种更自然的方式关注构建数值化的方法、模型和仿真系统。类似地，概率编程将处理概率和推断的过程对用户隐藏起来，从而使得用户更多的去关注模型构建和结果分析。
即便某些情况下不太可能从分析的角度得到后验，我们也有办法将后验计算出来，其中一些方法列举如下。
（1）非马尔科夫方法：
（2）马尔科夫方法：
如今，贝叶斯分析主要是通过马尔科夫链蒙特卡洛（Markov Chain Monte Carlo，MCMC）方法进行，同时变分方法也越来越流行，特别是在一些较大的数据集上。学习贝叶斯分析并不需要完全掌握这些方法，不过从概念层面上对它们的工作原理有一定了解会很有帮助，比如调试模型。
我们首先讨论基于非马尔科夫方法的推断引擎。非马尔科夫方法在低维的问题上要比马尔科夫方法更快，对于某些问题，这类方法非常有用，而对另外一些问题，这类方法只能提供真实后验的粗略近似，不过没关系，它们可以为马尔科夫方法提供一个不错的初始点，后面我们会介绍马尔科夫的含义。
网格计算是一种暴力穷举的方法。即便你无法计算出整个后验，你也可以根据一些点计算出先验和似然。假设我们要计算某个单参数模型的后验，网格近似可以按照如下方式进行：
视情况，我们可能会对计算结果进行归一化（把每个点的计算结果除以所有点的计算结果之和）。
很容易看出，选的点越多（网格越密）近似的结果就越好。事实上，如果使用无限多的点，我们可以得到准确的后验。网格计算的方法不能很好地适用于多参数（又或者称多维度）的场景，随着参数的增加，采样空间相比后验空间会急剧增加，换言之，我们花费了大量时间计算后验值，但对于估计后验却几乎没有帮助，因而使得该方法对于大多数统计学和数据科学的问题都不太实用[1]。
下面的代码用网格计算的方法解决第一章中的抛硬币问题：
def posterior_grid(grid_points=100, heads=6, tosses=9):
"""
A grid implementation for the coin-flip problem
"""
grid = np.linspace(0, 1, grid_points)
prior = np.repeat(5, grid_points)
likelihood = stats.binom.pmf(heads, tosses, grid)
unstd_posterior = likelihood * prior
posterior = unstd_posterior / unstd_posterior.sum()
return grid, posterior假设4次抛硬币实验中只有一次观测到正面朝上，那么就有：
points = 15
h, n = 1, 4
grid, posterior = posterior_grid(points, h, n)
plt.plot(grid, posterior, 'o-', label='heads = {}\ntosses = {}'.format(h, n))
plt.xlabel(r'$\theta$')
plt.legend(loc=0)二次近似，也称拉普拉斯方法或者正态近似，利用的是高斯分布来近似后验。该方法通常很有效，原因是后验分布众数附近的区域通常是服从正态分布的，事实上很多情况下就是高斯分布。二次近似的计算过程分为两步。首先，找出后验分布的峰值，该值可以通过一些最优化算法得到（有许多开箱即用的算法用来求解函数的最大值或最小值），这样我们得到了用来近似的高斯分布的均值；然后估计函数在众数附近的曲率，根据该曲率可以得出近似高斯分布的标准差。在介绍完PyMC3之后，我们会讲解如何使用该方法。
现代的贝叶斯统计学大多都采用马尔科夫方法（下一节会介绍），不过该方法对某些问题求解很慢而且不能很好地并行计算。一种简单的做法是同时运行多个马尔科夫链，然后将结果合并，不过对大多数问题而言这并不是一个合适的解决方案，如何找到有效的并行计算方式是当前的一个研究热点。
对于较大的数据集或是某些计算量很大的似然而言，变分方法是一个更好的选择。此外，这类方法能快速得到后验的近似，为MCMC方法提供初始点。
变分方法的基本思想是用一个更简单的分布去近似后验，这听起来有点像拉普拉斯近似，不过深入细节你会发现二者有很大不同。变分方法的最大缺点是我们必须对每个模型设计一个特定的算法，因而变分方法并不是一个通用的推断引擎，而是与模型相关的。
当然，许多人都在尝试将变分方法自动化。最近提出的一个方法是自动差分变分推断（Automatic Differentiation Variational Inference，ADVI）（查看http://arxiv.org/abs/1603.00788了解更多）。从概念层面上讲，ADVI是按下面步骤工作的。
（1）对参数进行变换从而使得它们位于整个实轴上。例如，假设一个参数限定为正数，对其求log后得到的值位于无界区间
内。
（2）用高斯分布对无界参数近似。需要注意，转换后参数空间里的高斯分布在原来参数空间内是非高斯的，这也是与拉普拉斯近似方法的不同点。
（3）采用某种优化算法使得高斯近似的结果尽可能接近后验，该过程通过最大化证据下界（Evidence Lower Bound，ELBO）实现。如何衡量两个分布的相似性以及证据下界的具体含义都是一些数学细节。
ADVI算法已经在PyMC3中实现了，本书后面会使用到。
有许多相关的方法都称MCMC方法。对于网格计算近似方法而言，我们需要根据给定的点计算出似然和先验，并且近似出整个后验分布。MCMC方法要比网格近似更好，这是因为其设计思想是将更多的时间用于高概率区域的计算，而在较低概率区间花费较少时间。事实上，MCMC方法会根据相对概率访问参数空间中的不同区间。如果区间A的概率是区间B的两倍，那么我们从区间A采样的次数可能是从区间B采样次数的两倍。因此，即使我们无法从分析的角度得到整个后验，我们也可以采用MCMC的方法从中采样，并且采样数越多，效果越好。
要理解MCMC方法，我们先将其拆分成两个MC部分，即蒙特卡洛部分和马尔科夫链部分。
蒙特卡洛这部分可以用随机数的应用来解释。蒙特卡洛方法是一系列应用非常广泛的算法，其思想是通过随机采样来计算或模拟给定过程。蒙特卡洛是位于摩纳哥公国的一个非常有名的城市，蒙特卡洛方法的开发者之一是Stanislaw Ulam。Stan正是基于这一核心思想——尽管很多问题都难以求解甚至无法准确用公式表达，但我们可以通过采样或者模拟来有效地研究。据传，起因是为了计算单轮游戏中的概率，解决该问题的方法之一是组合分析法；另一种是Stan声称的，进行多次单轮游戏，最后计算其中有多少次是我们感兴趣的。这听起来似乎是显而易见的，或者至少是相当合理的，比如，你可能已经用重采样的方法来解决统计问题。不过这个实验是早在70年前进行的，当时，第一台计算机才刚开始研发。该方法首次应用于一个核物理问题。如今，个人计算机都已经足够强大，用蒙特卡洛方法可以解决许多有趣的问题，因此，该方法也被广泛应用到科学、工程、工业以及艺术等各个方面。
在使用蒙特卡洛方法计算数值的例子中，一个教科书上非常经典的是估计π。实际使用中有更好的方法来计算π，不过这个例子仍然具有教学意义。我们可以通过以下过程估计π：
（1）在边长为2R的正方形内随机撒N个点。
（2）在正方形内画一个半径为R的圆，计算在圆内的点的个数。
（3）得出
的估计值
。
需要注意：如果一个点满足以下条件，我们认为该点位于圆内：
正方形的面积是(2R)2，圆的面积是πR2，因此二者的面积之比是4/π，而圆和正方形的面积分别正比于圆内的点的个数和总的点数N。我们可以通过几行简单的代码来模拟该蒙特卡洛过程计算π值，同时计算出估计值与实际值之间的相对误差：
N = 10000
x, y = np.random.uniform(-1, 1, size=(2, N))
inside = (x**2 + y**2) <= 1
pi = inside.sum()*4/N
error = abs((pi - np.pi)/pi)* 100
outside = np.invert(inside)
plt.plot(x[inside], y[inside], 'b.')
plt.plot(x[outside], y[outside], 'r.')
plt.plot(0, 0, label='$\hat \pi$ = {:4.3f}\nerror = {:4.3f}%'.
format(pi, error), alpha=0)
plt.axis('square')
plt.legend(frameon=True, framealpha=0.9, fontsize=16);上面的代码中，outside是用来绘图的，在计算 的过程中没有用到。另外一点需要澄清的是，由于这里用的是单位圆，因此在判断一个点是否在圆内时没有计算平方根。
马尔科夫链是一个数学对象，包含一系列状态以及状态之间的转移概率，如果每个状态转移到其他状态的概率只与当前状态相关，那么这个状态链就称为马尔科夫链。有这样一个马尔科夫链之后，我们可以任取一个初始点，然后根据状态转移概率进行随机游走。假设能够找到这样一个马尔科夫链，其状态转移概率正比于我们想要采样的分布（如贝叶斯分析中的后验分布），采样过程就变成了简单地在该状态链上移动的过程。那么，如何在不知道后验分布的情况下找到这样的状态链呢？有一个概念叫做细节平衡条件（Detailed Balance Condition），直观上讲，这个条件是说，我们需要采用一种可逆的方式移动（可逆过程是物理学中的一个常见的近似）。也就是说，从状态i转移到状态j的概率必须和状态j转移到状态i的概率相等。
总的来说就是，如果我们能够找到满足细节平衡条件的马尔科夫链，就可以保证从中采样得到的样本来自正确的分布。保证细节平衡的最流行的算法是Metropolis-Hasting算法。
为了更形象地理解这个算法，我们用下面这个例子来类比。假设我们想知道某个湖的水容量以及这个湖中最深的点，湖水很浑浊以至于没法通过肉眼来估计深度，而且这个湖相当大，网格近似法显然不是个好办法。为了找到一个采样策略，我们请来了两个好朋友小马和小萌。经过讨论之后想出了如下办法，我们需要一个船（当然，也可以是竹筏）和一个很长的棍子，这比声呐可便宜多了，而且我们已经把有限的钱都花在了船上。
（1）随机选一个点，然后将船开过去。
（2）用棍子测量湖的深度。
（3）将船移到另一个地点并重新测量。
（4）按如下方式比较两点的测量结果。
如何决定接受还是拒绝新的测量值呢？这里的一个技巧便是使用Metropolis-Hastings准则，即接受新的测量值的概率正比于新旧两点的测量值之比。
按照以上过程迭代下去，我们不仅可以得到整个湖的水容量和最深的点，而且可以得到整个湖底的近似曲率。你也许已经猜到了，在这个类比中，湖底的曲率其实就是后验的分布，而最深的点就是后验的众数。根据小马的说法，迭代的次数越多，近似的效果越好。
事实上，理论保证了在这种情形下，如果我们能采样无数次，最终能得到完整的后验。幸运地是，实际上对于很多问题而言，我们只需要相对较少地采样就可以得到一个相当准确的近似。
现在让我们从更正式的角度来看看该算法。对于很多分布而言（如高斯分布），我们有相当高效的算法对其采样，但对于一些其他分布，情况就变了。Metropolis-Hastings算法使得我们能够从任意分布中以概率p(x)得到采样值，只要我们能算出某个与p(x)成比例的值。这一点很有用，因为在类似贝叶斯统计的许多问题中，最难的部分是计算归一化因子，也就是贝叶斯定理中的分母。Metropolis-Hastings算法的步骤如下。
（1）给参数xi赋一个初始值，通常是随机初始化或者使用某些经验值。
（2）从某个简单的分布
中选一个新的值
，如高斯分布或者均匀分布。这一步可以看做是对状态
的扰动。
（3）根据Metropolis-Hastings准则计算接受一个新的参数值的概率：
（4）从位于区间[0,1]内的均匀分布中随机选一个值，如果第（3）步中得到的概率值比该值大，那么就接受新的值，否则仍保持原来的值。
（5）回到第（2）步重新迭代，直到我们有足够多的样本，稍后会解释什么叫足够多。
有几点需要注意。
如果选取的分布是对称的，那么可以得到
，这通常称为Metropolis准则。
步骤（3）和步骤（4）表明：我们总是会转移到一个比当前状态（或参数）概率更大的状态（或参数），对于概率更小的，则会以xi+1与xi之比的概率接受。该准则中的接受步骤使得采样过程相比网格近似方法更高效，同时保证了采样的准确性。
目标分布（贝叶斯统计中的后验分布）是通过记录下来的采样值来近似的。如果我们接受转移到新的状态xi+1，那么我们就记录该采样值xi+1。如果我们拒绝转移到xi+1，那么我们就记录xi。
最后，我们会得到一连串记录值，有时候也称采样链或者迹。如果一切都正常进行，那么这些采样值就是后验的近似。在采样链中出现次数最多的值就是对应后验中最可能的值。该过程的一个优点是：后验分析很简单，我们把对后验求积分的过程转化成了对采样链所构成的向量求和的过程。
下面的代码展示了Metropolis算法的一个基本实现。这段代码并不是为了解决什么实际问题，只是在这里用来演示，如果我们知道怎么计算给定点的函数值，我们就可能得到该函数的采样。需要注意下面的代码中不包含贝叶斯相关的部分，既没有先验也没有数据。要知道，MCMC是一类能够用于解决很多问题的通用方法。例如，在一个（非贝叶斯的）分子模型中，我们可能需要一个函数来计算在某个状态x下系统的能量而不是简单地调用func.pdf(x)函数。
metropolis函数的第一个参数是一个SciPy的分布，假设我们不知道如何从中直接采样。
def metropolis(func, steps=10000):
"""A very simple Metropolis implementation"""
samples = np.zeros(steps)
old_x = func.mean()
old_prob = func.pdf(old_x)
for i in range(steps):
new_x = old_x + np.random.normal(0, 0.5)
new_prob = func.pdf(new_x)
acceptance = new_prob/old_prob
if acceptance >= np.random.random():
samples[i] = new_x
old_x = new_x
old_prob = new_prob
else:
samples[i] = old_x
return samples下面的例子中，我们把func定义成一个beta函数，因为beta函数可以通过改变参数调整分布的形状。我们把metropolis函数的结果用直方图画出来，同时用红色的线表示真实的分布。
func = stats.beta(0.4, 2)
samples = metropolis(func=func)
x = np.linspace(0.01, .99, 100)
y = func.pdf(x)
plt.xlim(0, 1)
plt.plot(x, y, 'r-', lw=3, label='True distribution')
plt.hist(samples, bins=30, normed=True, label='Estimated distribution')
plt.xlabel('$x$', fontsize=14)
plt.ylabel('$pdf(x)$', fontsize=14)
plt.legend(fontsize=14)现在你应该从概念上掌握了Metropolis-Hastings算法。也许你需要回过头去重新阅读前面几页才能完全消化。此外，我还强烈建议阅读PyMC3核心作者之一写的博文http://twiecki.github.io/blog/2015/11/10/mcmc-sampling/[2]。他用一个简单的例子实现了metropolis方法，并将其用于求解后验分布，文中用非常好看的图展示了采样的过程，同时简单讨论了最初选取的步长是如何影响结果的。
MCMC方法，包括Metropolis-Hastings，都在理论上保证如果采样次数足够多，最终会得到后验分布的准确近似。不过，实际中想要采样足够多次可能需要相当长的时间，因此，人们提出了一些Metropolis-Hastings算法的替代方案。这些替代方案，包括Metropolis-Hastings算法本身，最初都是用来解决统计力学中的问题。统计力学是物理学的一个分支，主要研究原子和分子系统的特性。汉密尔顿蒙特卡洛方法，又称混合蒙特卡洛（Hybrid Monte Carlo，HMC），是这类改进方案之一。简单来说，汉密尔顿这个词描述的是物理系统的总能量，而另外一个名称中的“混合”是指将Metropolis-Hastings算法与分子力学（分子系统中广泛应用的一种仿真技巧）相结合。HMC方法本质上和Metropolis-Hastings是一样的，改进的地方在于：原来是随机放置小船，现在有了一个更聪明的办法，将小船沿着湖底方向放置。为什么这个做法更聪明？因为这样做避免了Metropolis-Hastings算法的主要问题之一：探索得太慢而且采样结果自相关（因为大多数采样结果都被拒绝了）。
那么，如何才能不必深入其数学细节而理解汉密尔顿蒙特卡洛方法呢？假设我们还是在湖面上坐着船，为了决定下一步将要去哪，我们从当前位置往湖底扔了一个球，受“球状奶牛”的启发[3]，我们假设球面是理想的，没有摩擦，因而不会被泥巴和水减速。扔下球之后，让它滚一小会儿，然后把船划到球所在的位置。现在利用Metropolis-Hastings算法中提到的Metropolis准则来选择接受或者拒绝，重复整个过程一定次数。改进后的过程有更高的概率接受新的位置，即使它们的位置相比前一位置距离较远。
现在跳出我们的思维实验，回到现实中来。基于汉密尔顿的方法需要计算函数的梯度。梯度是在多个维度上导数的推广。我们可以用梯度信息来模拟球在曲面上移动的过程。因此，我们面临一个权衡；HMC计算过程要比Metropolis-Hastings更复杂，但是被接受概率更高。对于一些复杂的问题，HMC方法更合适一些。HMC方法的另一个缺点是：想要得到很好的采样需要指定一些参数。如果手动指定，需要反复尝试，这要求使用者有一定的经验。幸运地是，PyMC3中有一个相对较新的不掉向采样算法，该方法被证实可以有效提升HMC方法的采样效率，同时不必手动调整参数。
MCMC的方法很多，而且人们还在不断提出新的方法。如果你认为你能提升采样效率，会有很多人对你的想法感兴趣。讨论所有这类方法及其优缺点显然超出了本书的范围。不过，有些方法仍值得一提，因为你可能经常会听到人们讨论它们，所以最好是了解下他们都在讨论些什么。
在分子系统模拟中广泛应用的另外一种采样方法是副本交换（Replica Exchange），也称并行退火（Parallel Tempering）或者Metropolis Coupled MCMC（MC3；好多MC……）。该方法的基本思想是并行模拟多个不同的副本，每个副本都按照Metropolis-Hastings算法执行。多个副本之间的唯一不同是一个叫做温度的参数（又受到物理学的启发！），该参数用来控制接受低概率采样点的可能性。每隔一段时间，该方法都尝试在多个副本之间进行切换，切换过程同时遵循Metropolis-Hastings 准则来接受或拒绝，只不过现在考虑的是不同副本之间的温度。切换过程可以在状态链上进行随机切换，不过，通常更倾向于在相邻副本之间切换，也就是说，具有相似温度的副本会有更高的接受概率。该方法的直观理解是：如果我们提高温度，接受新位置的概率也会提升，反之则会降低。温度更高的副本可以更自由地探索系统，因为这些副本的表面会变得相当平坦从而更容易探索。对于温度无限高的副本，所有状态都是等概率的。副本之间的切换避免了较低温度的副本陷在局部最低点，因而该方法很适合探索有多个最低点的系统。
PyMC3是一个用于概率编程的Python库，当前最新的版本号是2016年10月4号发布的3.0.rc2。PyMC3提供了一套非常简洁直观的语法，非常接近统计学中描述概率模型的语法，可读性很高。PyMC3是用Python写的，其中的核心计算部分基于NumPy和Theano。Theano是一个用于深度学习的Python库，可以高效地定义、优化和求解多维数组的数学表达式。PyMC3使用Theano的主要原因是某些采样算法（如NUTS）需要计算梯度，而Theano可以很方便地进行自动求导。而且，Theano将Python代码转化成了C代码，因而PyMC3的速度相当快。关于Theano就需要了解这些，如果你想深入学习更多可以阅读Theano官网上的教程http://deeplearning.net/software/theano/tutorial/ index.html#tutorial。
让我们重新回顾下抛硬币问题，这次我们使用PyMC3。首先我们需要获取数据，这里我们使用手动构造的数据。由于数据是我们自己生成，所以知道真实的参数，以下代码中用theta_real变量表示。显然，在真实数据中，我们并不知道参数的真实值，而是要将其估计出来。
np.random.seed(123)
n_experiments = 4
theta_real = 0.35
data = stats.bernoulli.rvs(p=theta_real, size=n_experiments)
print(data)
array([1, 0, 0, 0])现在有了数据，需要再指定模型。回想一下，模型可以通过指定似然和先验的概率分布完成。对于似然，我们可以用参数分别为
和
的二项分布来描述，对于先验，我们可以用参数为
的beta分布描述。这个beta分布与[0,1]区间内的均匀分布是一样的。我们可以用数学表达式描述如下：
这个统计模型与PyMC3的语法几乎一一对应。第1行代码先构建了一个模型的容器，PyMC3使用with语法将所有位于该语法块内的代码都指向同一个模型，你可以把它看作是简化模型描述的“语法糖”，这里将模型命名为our_first_model。第2行代码指定了先验，可以看到，语法与数学表示很接近。我们把随机变量命名为，需要注意的是，这里变量名与Beta函数的第1个参数名一样；保持相同的名字是个好习惯，这样能避免混淆。然后，我们通过变量名从后验采样中提取信息。这里变量是一个随机变量，我们可以将该变量看做是从某个分布（在这里是beta分布）中生成数值的方法而不是某个具体的值。第3行代码用跟先验相同的语法描述了似然，唯一不同的是我们用observed变量传递了观测到的数据，这样就告诉了PyMC3我们的似然。其中，data可以是一个Python列表或者Numpy数组或者Pandas的DataFrame。这样我们就完成了模型的描述。
with pm.Model() as our_first_model:
theta = pm.Beta('theta', alpha=1, beta=1)
y = pm.Bernoulli('y', p=theta, observed=data)对于抛硬币这个问题，后验分布既可以从分析的角度计算出来，也可以通过PyMC3用几行代码从后验的采样中得到。代码中的第1行，调用了find_MAP函数，该函数调用SciPy中常用的优化函数尝试返回最大后验（Maximum a Posteriori，MAP）。调用find_MAP是可选的，有时候其返回值能够为采样方法提供一个不错的初始点，不过有时候却并没有多大用，因此大多数时候会避免使用它。然后，下一行定义了采样方法。这里用的是Metropolis-Hastings算法，函数名取了简写Metropolis。PyMC3可以让我们将不同的采样器赋给不同的随机变量；眼下我们的模型只有一个参数，不过后面我们会有多个参数。我们也可以省略该行，PyMC3会根据不同参数的特性自动地赋予一个采样器，例如，NUTS算法只对连续变量有用，因而不能用于离散的变量，Metropolis算法能够处理离散的变量，而另外一些类型的变量有专门的采样方法。总的来说，我们可以让PyMC3为我们选一个采样方法。最后一行是执行推断，其中第1个参数是采样次数，第2个和第3个参数分别是采样方法和初始点，可以看到这两个参数是可选的。
start = pm.find_MAP()
step = pm.Metropolis()
trace = pm.sample(1000, step=step, start=start)这样，只需要几行代码我们就完成了整个模型的描述和推断。感谢PyMC3的开发者们为我们提供了这么棒的库。
现在我们根据有限数量的样本对后验做出了近似，接下来要做的第一件事就是检查我们的近似是否合理。我们可以做一些测试，有些是可视化的，有些是定量的。这些测试会尝试从样本中发现问题，但并不能证明我们得到的分布是正确的，它们只能提供证据证明样本看起来是合理的。如果我们通过样本发现了问题，解决办法有如下几种。
本书剩余部分将会详细讲解这些方案。
通常，我们要做的第一件事就是查看结果长什么样，traceplot函数非常适合该任务：
burnin = 100
chain = trace[burnin:]
pm.traceplot(chain, lines={'theta':theta_real});对于未观测到的变量，我们得到了两幅图。左图是一个核密度估计（Kernel Density Estimation，KDE）图，可以看做是平滑之后的直方图。右图描绘的是每一步采样过程中得到的采样值。注意图中红色的线表示变量theta_real的值。
在得到这些图之后，我们需要观察什么呢？首先，KDE图看起来应该是光滑的曲线。通常，随着数据的增加，根据中心极限定理[4]，参数的分布会趋近于高斯分布。当然，这并不一定是正确的。右侧的图看起来应该像白噪声，也就是说有很好的混合度（mixing），我们看不到任何可以识别的模式，也看不到向上或者向下的曲线，相反，我们希望看到曲线在某个值附近震荡。对于多峰分布或者离散分布，我们希望曲线不要在某个值或区域停留过多时间，我们希望看到采样值在多个区间自由移动。此外，我们希望迹表现出稳定的相似性，也就是说，前10%看起来跟后50%或者10%差不多。再次强调，我们不希望看到规律的模式，相反我们期望看到的是噪声。下图展示了一些迹呈现较好混合度（右侧）与较差混合度（左侧）的对比。
如果迹的前面部分跟其他部分看起来不太一样，那就意味着需要进行老化处理，如果其他部分没有呈现稳定的相似性或者可以看到某种模式，那这意味着需要更多的采样，或者需要更换采样方法或者参数化方法。对于一些复杂的模型，我们可能需要结合使用前面所有的策略。
PyMC3可以实现并行地运行一个模型多次，因而对同一个参数可以得到多条并行的迹。这可以通过在采样函数中指定njobs参数实现。此时使用traceplot函数，便可在同一幅图中得到同一个参数的所有迹。由于每组迹都是相互独立的，所有的迹看起来都应该差不多。除了检查收敛性之外，这些并行的迹也可以用于推断，我们可以将这些并行的迹组合起来提升采样大小而不是扔掉多余的迹：
with our_first_model:
step = pm.Metropolis()
multi_trace = pm.sample(1000, step=step, njobs=4)
burnin = 0
multi_chain = multi_trace[burnin:]
pm.traceplot(multi_chain, lines={'theta':theta_real});一种定量地检测收敛性的方法是Gelman-Rubin检验。该检验的思想是比较不同迹之间的差异和迹内部的差异，因此，需要多组迹来进行该检验。理想状态下，我们希望得到
。根据经验，我们认为如果得到的值低于1.1，那么可以认为是收敛的了，更高的值则意味着没有收敛：
pm.gelman_rubin(multi_chain)
{'theta': 1.0074579751170656, 'theta_logodds': 1.009770031607315}我们还可以用forestplot函数将
和每个参数的均值、50%HPD和95%HPD可视化地表示出来：
pm.forestplot(multi_chain, varnames=['theta']);函数summary提供了对后验的文字描述，它可以提供后验的均值、标准差和HPD区间：
pm.summary(multi_chain)
theta:
Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
0.339 0.173 0.006 [0.037, 0.659]
Posterior quantiles:
2.5 25 50
75 97.5
|--------------|==============|==============|--------------|
0.063 0.206 0.318
0.455 0.709此外，df_summary函数会返回类似的结果，不过类型是Pandas中的DataFrame：
pm.df_summary(multi_chain)| mean |
sd |
mc_error |
hpd_2.5 |
hpd_97.5 |
|
| theta |
0.33883 |
0.17305 |
0.00592 |
0.03681 |
0.65916 |
其中，返回值之一是mc_error，这是对采样引入误差的估计值，该值考虑的是所有的采样值并非真的彼此独立。mc_error是迹中不同块的均值的标准差，每一块是迹中的一部分：
该误差值显然低于我们结果的准确度。由于采样方法是随机的，每次重跑的时候，summary或者df_summary返回的值都会不同，不过没关系，mc_error的值应该是相似的，如果返回的值有很大不同，则说明我们可能需要更多的样本。
最理想的采样应该不会是自相关的，也就是说，某一点的值应该与其他点的值是相互独立的。在实际中，从MCMC方法（特别是Metropolis-Hastings）中得到的采样值是自相关的。由于参数之间的相互依赖关系，可能模型会导致更多的自相关采样。PyMC3有一个很方便的函数用来描述自相关。
pm.autocorrplot(chain)该图显示了采样值与相邻连续点（最多100个）之间的平均相关性。理想状态下，我们不会看到自相关性，实际中我们希望看到自相关性降低到较低水平。参数越自相关，要达到指定精度的采样次数就需要越多，也就是说，自相关性不利于降低采样次数。
一个有自相关性的采样要比没有自相关性的采样所包含的信息量更少，因此，给定采样大小和采样的自相关性之后，我们可以尝试估计出该采样的大小为多少时，该采样没有自相关性而且包含的信息量不变，该值称为有效采样大小。理想情况下，两个值是一模一样的；二者越接近，我们的采样效率越高。有效采样大小可以作为我们的一个参考，如果我们想要估计出一个分布的均值，我们需要的最小采样数至少为100；如果想要估计出依赖于尾部分布的量，比如可信区间的边界，那么我们可能需要1000到10000次采样。
pm.effective_n(multi_chain)['theta']
667显然，提高采样效率的一个方法是换一个更好的采样方法；另一个办法是转换数据或者对模型重新设计参数，此外，还有一个常用的办法是对采样链压缩。所谓压缩其实就是每隔k个观测值取一个，在Python中我们称为切片。压缩会降低自相关性，但代价是同时降低了样本量。因此，实际使用中通常更倾向于增加样本量而不是切片。不过有时候，压缩会很有用，比如降低存储空间。如果仍不能避免高自相关性，我们就只能算出更长的采样链，模型中的参数很多的话，存储量会是个问题。而且，我们可能还会对后验做一些计算量很大的后处理，此时在自相关性尽可能小的前提下，采样数量的大小就显得尤为重要。
目前为止，所有的诊断测试都是经验性而非绝对的。实际使用中，我们会先运行一些测试，如果看起来没什么问题，我们就继续往下分析。如果发现了一些问题，就需要回过头解决它们，这也是建模过程的一部分。要知道，进行收敛性检查并非贝叶斯理论的一部分，由于我们是用数值的方式在计算后验，因而这只是贝叶斯实践过程中的一部分。
我们已经知道，贝叶斯分析的结果是后验分布，其包含了在已有数据和模型下，参数的所有信息。我们可以使用PyMC3中的plot_posterior函数对后验分布进行可视化总结，这个函数的核心参数是一个PyMC3的迹和或者一个NumPy的数组，默认情况下，该函数会画出参数的直方图以及分布的均值，此外图像的底端还有一个黑色的粗线用来表示95%HPD区间。可以通过设置alpha_level参数来改变HPD区间。我们将这类图称为Kruschke图，这是因为John K. Kruschke第一次在他的《Doing Bayesian Data Analysis》一书中引入了这种图。
pm.plot_posterior(chain, kde_plot=True)有时候，仅仅描述后验还不够，我们还需要根据推断结果做决策，即将连续的估计值收敛到一个二值化结果上：是或不是、受污染了还是安全的等等。回到抛硬币问题上，我们需要回答硬币是不是公平的。一枚公平的硬币是指的值为0.5，严格来说，出现这种情况的概率是0，因而，实际中我们会对定义稍稍放松，假如一枚硬币的值在0.5左右，我们就认为这枚硬币是公平的。这里“左右”的具体含义依赖于具体的问题，并没有一个满足所有问题的普适准则。因此决策也是偏主观的，我们的任务就是根据我们的目标做出最可能的决策。
直观上，一个明智的做法是将HPD区间与我们感兴趣的值进行比较，在我们的例子中，该值是0.5。前面的图中可以看出HPD的范围是0.06～0.71，包含0.5这个值，不过根据后验分布来看，硬币似乎倾向于反面朝上，我们无法就此裁定一个硬币是公平的。或许，我们需要收集更多的数据来降低数据的分散程度，从而得到一个更确定的决策；又或者是因为我们漏掉了某些关键信息，以至我们没能找到更完备的先验。
基于后验做决策的一种方案是实用等价区间（Region Of Practical Equivalence，ROPE），其实就是在感兴趣值附近的一个区间，例如我们可以说[0.45,0.55]是0.5的一个实用性等价区间。同样，ROPE是根据实际情况决定的。接下来我们可以将ROPE与HPD对比，结果至少有以下3种情况。
当然，如果选择区间[0,1]作为ROPE，那么不管结果怎样我们都会说这枚硬币是公平的，不过恐怕没人会同意我们对ROPE的定义。plot_posterior函数可以用来画ROPE。从图中可以看到，ROPE是一段较宽的半透明的红色线段，同时上面有两个数值表示ROPE的两个端点。
pm.plot_posterior(chain, kde_plot=True, rope=[0.45, 0.55])我们还可以给plot_posterior传递一个参考值，例如0.5，用来和后验进行对比。从图中可以看出我们会得到一个绿色的垂直线以及大于该值和小于该值的后验比例。
pm.plot_posterior(chain, kde_plot=True, ref_val=0.5)关于如何使用ROPE的更多细节，你可以阅读John Kruschke写的《Doing Bayesian Data Analysis》一书的第12章。这一章还讨论了在贝叶斯框架下如何做假设检验，以及一些（贝叶斯或者非贝叶斯的）假设检验方面的警告。
如果你觉得ROPE准则有些杂乱，想要更正式一些的形式，那么损失函数就是你想要的。做好决策很重要的一点是，参数的估计值要有很高的精度，同时还要考虑到预测出错的代价。对于收益/损失的权衡可以从数学形式上用代价函数来描述，有时候也称为损失函数。损失函数刻画的是当Y（硬币是不公平的）成立时预测X（硬币是公平的）成立时的代价。在许多问题中，决策的代价函数是不对称的，例如，在决定5岁以下的儿童是否应该接种某种疫苗这件事上，决定接种或者不接种可能造成完全不同的影响，一旦做出错误的决策，可能会导致上千人死亡，而假如能决定接种某种非常安全同时又相对便宜的疫苗，则可能避免一场健康危机。人们对如何在有限信息下做出决策的问题已经研究许多年了，该学科被称为决策理论。
本章，我们学习了概率编程，同时体会到了推断引擎的强大力量。我们从概念上讨论了MCMC方法的核心思想及其在现代贝叶斯数据分析中的地位。此外，我们第一次见证了PyMC3的强大和易用性。我们还重新回顾了前一章中的抛硬币问题，这次，我们使用PyMC3重新定义并解决了这个问题，同时还做了建模过程中非常重要的模型检查和模型诊断。
下一章会继续学习贝叶斯分析的一些技巧，我们将学习如何处理多个参数的模型，以及如何协调多个参数之间的层次关系。
（1）用其他先验尝试网格算法。例如，尝试将代码改成prior = (grid <= 0.5).astype(int)或者prior = abs(grid - 0.5)，或者你可以将其定义成你所想的任意形式。换一些其他数据重复实验，例如增加总的样本数或者将样本数根据你看到的正面朝上的次数适当调整。
（2）在我们估计p值的代码中，将N固定，重复多次实验。注意由于我们使用了随机数，每次的结果都会不一样，不过检查一下可以发现误差应该差不多。尝试把N调大然后再重跑，你能猜出N与误差之间的关系吗？你可能需要修改下代码，将N作为一个参数传给计算误差的函数，这样方便估计N与误差之间的关系。对于相同的N值，可以重复跑多次后计算误差的平均值和这些误差的标准差，然后使用matplotlib中的errobar()函数将它们画出来。可以尝试一些类似100、1000、10000的数作为N的值（每次提高一个数量级）。
（3）修改传给metropolis函数的参数。尝试使用第1章中用到的先验。将这段代码与网格计算进行比较，看一下哪部分需要修改后才能用于解决贝叶斯推断问题。
（4）将前一题的答案与下面链接中Thomas Wiechi的代码进行比较：http://twiecki.github.io/blog/2015/11/10/mcmc-sampling/。
（5）修改beta先验分布的参数以匹配前1章中的分布，试比较2者的结果。将beta分布替换成区间为[0,1]的均匀分布，其结果是否与beta(a=1,b=1)相等？采样速度相比是更快？还是更慢？或者是相等？如果换成更大的区间，比如[-1,2]呢？模型是否能正常运行？你得到的错误是什么意思？如果你不调用find_MAP()函数的话呢？（记得将sample函数的第1个参数也去掉，如果你是在Jupyter/IPython记事本中运行代码，尤其需要注意这点。）
（6）修改退化的数量。尝试将退化的数量修改为0或者500，还可以尝试不用find_MAP()函数，结果的差别有多大？提示：这里采样的模型非常简单，记住这个练习，后面在其他章节中我们还会在更复杂的模型上进行尝试。
（7）使用自己的数据。用你自己感兴趣的数据重新跑一边本章的内容，这个练习适用于本书剩余的所有章节。
（8）阅读PyMC3文档中的煤矿灾变模型：http://pymc-devs.github.io/pymc3/notebooks/getting_started.html#Case-study-2:-Coal-mining-disasters，试着自己实现并运行该模型。
除了每章最后的练习之外，你还可以将已经学到的内容应用到你感兴趣的问题上。也许，你需要重新定义你的问题，或者需要扩展或者修改你已经学到的模型。试着修改模型，如果你觉得这个任务已经超出了你实际掌握的部分，那么先将问题记下来，等读完本书其余章节后再回过头来重新思考这些问题。如果最后本书仍然无法解决你的问题，那么你可以查看下PyMC3的例子（http://pymc-devs.github.io/pymc3/examples.html），或者在PyMC3的论坛（https://discourse.pymc.io）上提问。
[1] 可阅读https://zh.wikipedia.org/wiki/维数灾难，了解更多信息。——译者注
[2] 该文章notebook形式的翻译见https://github.com/findmyway/Bayesian-Analysis-with-Python/blob/master/MCMC-sampling-for-dummies.ipynb。——译者注
[3] “球状奶牛”起源于一个物理学家的笑话，具体含义可自行搜索维基百科。——译者注
[4] 此处已根据原书的勘误更正。——译者注
[5] 第2版也有人转成了Python，https://github.com/JWarmenhoven/DBDA-python。——译者注
