作者：加 Cameron Davidson Pilon 卡梅隆 戴维森 皮隆
出版社：人民邮电出版社有限公司
格式: AZW3, DOCX, EPUB, MOBI, PDF, TXT
贝叶斯方法：概率编程与贝叶斯推断试读：
前言
贝叶斯方法是一种常用的推断方法，然而对读者来说它通常隐藏在乏味的数学分析章节背后。关于贝叶斯推断的书通常包含两到三章关于概率论的内容，然后才会阐述什么是贝叶斯推断。不幸的是，由于大多数贝叶斯模型在数学上难以处理，这些书只会为读者展示简单、人造的例子。这会导致贝叶斯推断给读者留下“那又如何？”的印象。实际上，这曾是我自己的先验观点。
最近贝叶斯方法在一些机器学习竞赛上取得了成功，让我决定再次研究这一主题。然而即便以我的数学功底，我也花了整整3天时间来阅读范例，并试图将它们汇总起来以便理解这一方法。那时并没有足够的文献将理论和实际结合起来。而让我产生理解偏差的正是由于没能将贝叶斯数学理论和概率编程实践结合起来。当然，如今读者已经无需再遭遇我当时的情景。本书就是为了填补这一空缺而编写的。
如果我们最终是要进行贝叶斯推断，那么一方面我们可以采用数学分析来实现这一目的，而另一方面，随着计算成本的下降，我们已经可以通过概率编程来完成这一任务。后一种方法更加有用，因为它避免了在每一步介入数学干预，而这也使得进行贝叶斯推断不再以通常很棘手的数学分析为前提。简而言之，后一种计算途径，是从问题起点经过小幅中间步骤到达问题终点，而前一种途径则大幅跃进，并通常最后远离目标。此外，如果没有深厚的数学功底，也根本无法完成前一种途径所需要的数学分析。
本书首先从计算和理解的角度，而后从数学分析的角度对贝叶斯推断进行了介绍。当然，作为一本入门书籍，本书将停留在入门阶段。对于受过数学训练的人来说，本书产生的疑问可通过其他偏重数学分析的书来解答。对于缺少数学背景的爱好者，或是仅对贝叶斯方法的实践而非数学理论感兴趣的读者来说，本书足以胜任且蕴含趣味。
选择PyMC作为概率编程语言有两方面原因。首先，在写本书之时，并没有集中的关于PyMC的说明和实例等资料。官方文档面向具有贝叶斯推断和概率编程背景知识的人。而我们希望本书可以鼓励各个层次的人了解PyMC。其次，随着近来用Python实现科学计算框架的流行及其核心进展，PyMC可能很快会成为核心组件之一。
PyMC的运行需要一些依赖库，包括NumPy以及可选的SciPy。为了不产生限制，本书的实例只依赖PyMC、NumPy、SciPy和Matplotlib。
本书内容安排如下。第1章介绍贝叶斯推断方法以及与其他推断方法的比较。我们会看到第一个贝叶斯模型，并对其进行建立和训练。第2章以实例为重点，讲述如何用PyMC构建模型。第3章介绍计算推断背后的一个强大算法——马尔科夫链蒙特卡洛，以及一些贝叶斯模型的调试技术。在第4章里，我们再次回到推断的样本量问题上，并解释为何样本量大小如此重要。第5章介绍强大的损失函数，它将在真实世界的问题与数学推断之间建立连接。我们将在第6章回顾贝叶斯先验，并通过启发式的方法找到先验的更优解。最后，我们在第7章探索如何将贝叶斯推断用于A/B测试。
本书用到的所有数据集都可以从这里获得：https:// github.com/CamDavidsonPilon/ Probabilistic-Programming-and- Bayesian-Methods-for-Hackers。致谢
谨将本书献给许多重要的人：我的父母、兄弟和我最好的朋友。此外，本书要献给开源社区，是他们每天都在为我们默默地做出贡献。
我要感谢参与到本书写作里的人们。首先要感谢的是为这本书的网络版做出贡献的人。这些作者里很多人提交的代码、思路或文章使本书得以完成。然后，我要感谢对本书进行审校的Robert Mauriello和Tobi Bosede，他们牺牲自己的时间来把一些难以理解的抽象概念变得浅显易懂，并缩减内容以便更好地阅读体验。最后，我要感谢我的朋友以及同事，他们在整个过程中给与我支持。关于作者
Cameron Davidson-Pilon，接触过数学在多个领域的应用——从基因和疾病的动态演化，到金融价格的随机模型。他对于开源社区最主要的贡献包括这本书以及lifelines项目。Cameron成长于加拿大的安大略省圭尔夫市，而就读于滑铁卢大学以及莫斯科独立大学。如今他住在安大略省渥太华市，并在电商领军者Shopify工作。关于译者
辛愿，浙江大学硕士毕业，腾讯公司基础研究高级工程师，舆情系统开发经理。曾在百度从事推荐系统、用户画像、数据采集等相关研究工作，拥有多项专利，组织过上海大数据技术沙龙。目前专注于文本挖掘、舆情分析、智能聊天机器人等相关领域。
钟黎，腾讯公司研究员。曾在中国科学院、微软亚洲研究院、IBM研究院（新加坡）从事图像处理、语音处理、机器学习等相关研究工作，拥有多项专利，目前聚焦在自然语言处理、深度学习和人工智能等相关领域。
欧阳婷，华南理工大学硕士毕业，腾讯公司后台策略工程师。在电信、互联网行业参与过推荐系统、资源优化、KPI预测、用户画像等相关项目，拥有多项专利，目前聚焦在欺诈检测、时序分析、业务安全等相关领域。关于审校者
余凯博士，地平线机器人技术创始人、CEO，国际著名机器学习专家，中组部国家“千人计划”专家，中国人工智能学会副秘书长。余博士是前百度研究院执行院长，创建了百度深度学习研究院。他在百度所领导的团队在广告变现、搜索排序、语音识别、计算机视觉等领域做出杰出贡献，创纪录地连续三次获得公司最高荣誉——“百度最高奖”。他还创建了中国公司第一个自动驾驶项目，后发展为百度自动驾驶事业部。回国前，余博士在德国和美国的工业界工作了12年，服务于西门子、微软、NEC硅谷实验室等机构。他发表的学术论文被国际同行引用超过12 000次，2011年在斯坦福大学计算机系主讲课程“CS121: Introduction to Artificial Intelligence”。余博士在南京大学获得学士和硕士学位，在德国慕尼黑大学获得计算机科学博士学位。
岳亚丁博士，腾讯公司专家研究员，腾讯技术职级评委会基础研究岗位的负责委员。岳博士拥有19年在金融、电信、互联网行业的数据挖掘经验，主导或参与过用户画像、在线广告、推荐系统、CRM、欺诈检测、KPI预测等多种项目。他曾在微软（加拿大）从事行为定向广告的模型研发，另有11年的工程结构、海洋水文气象的力学研究及应用的工作经历。岳博士在华中科技大学获得力学博士学位，在美国圣约瑟夫大学获得计算机科学硕士学位。第1章 贝叶斯推断的哲学1.1 引言
尽管你已是一个编程老手，但bug仍有可能在代码中存在。于是，在实现了一段特别难的算法之后，你决定先来一个简单的测试用例。这个用例通过了。接着你用了一个稍微复杂的测试用例。再次通过了。接下来更难的测试用例也通过了。这时，你开始觉得也许这段代码已经没有bug了。
如果你这样想，那么恭喜你：你已经在用贝叶斯的方式思考！简单地说，贝叶斯推断是通过新得到的证据不断地更新你的信念。贝叶斯推断很少会做出绝对的判断，但可以做出非常可信的判断。在上面的例子中，我们永远无法100%肯定我们的代码是无缺陷的，除非我们测试每一种可能出现的情形，这在实践中几乎不可能。但是，我们可以对代码进行大量的测试，如果每一次测试都通过了，我们更有把握觉得这段代码是没问题的。贝叶斯推断的工作方式就在这里：我们会随着新的证据不断更新之前的信念，但很少做出绝对的判断，除非所有其他的可能都被一一排除。1.1.1 贝叶斯思维
和更传统的统计推断不同，贝叶斯推断会保留不确定性。乍听起来，这像一门糟糕的统计方法，难道不是所有的统计都是期望从随机性里推断出确定性吗？要协调这些不一致，我们首先需要像贝叶斯派一样思考。
在贝叶斯派的世界观中，概率是被解释为我们对一件事情发生的相信程度，换句话说，这表明了我们对此事件发生的信心。事实上，我们一会儿就会发现，这就是概率的自然的解释。
为了更清晰地论述，让我们看看来自频率派关于概率的另一种解释。频率派是更古典的统计学派，他们认为概率是事件在长时间内发生的频率。例如，在频率派的哲学语境里，飞机事故的概率指的是长期来看，飞机事故的频率值。对
许多事件来说，这样解释概率是有逻辑的，但对某些没有长期频率的事件来说，这样解释是难以理解的。试想一下：在总统选举时，我们经常提及某某候选人获选的概率，但选举本身只发生一次！频率论为了解决这个问题，提出了“替代现实”的说法，从而用在所有这些替代的“现实”中发生的频率定义了这个概率。
贝叶斯派，在另一方面，有更直观的方式。它把概率解释成是对事件发生的信心。简单地说，概率是观点的概述。某人把概率0赋给某个事件的发生，表明他完全确定此事不会发生；相反，如果赋的概率值是1，则表明他十分肯定此事一定会发生。0和1之间的概率值可以表明此事可能发生的权重。这个概率定义可以和飞机事故概率一致。如果除去所有外部信息，一个人对飞机事故发生的信心应该等同于他了解到的飞机事故的频率。同样，贝叶斯概率的定义也能适用于总统选举，并且使得概率（信心）更加有意义：你对候选人A获胜的信心有多少？
请注意，在之前，我们提到每个人都可以给事件赋概率值，而不是存在某个唯一的概率值。这很有趣，因为这个定义为个人之间的差异留有余地。这正和现实天然契合：不同的人即便对同一事件发生的信心也可以有不同的值，因为他们拥有不同的信息。这些不同并不说明其他人是错误的。考虑下面的例子。
1．在抛硬币中我们同时猜测结果。假设硬币没有被做手脚，我和你应该都相信正反面的概率都是0.5。但假设我偷看了这个结果，不管是正面还是反面，我都赋予这个结果1的概率值。现在，你认为正面的概率是多少？很显然，我额外的信息（偷看）并不能改变硬币的结果，但使得我们对结果赋予的概率值不同了。
2．你的代码中也许有一个bug，但我们都无法确定它的存在，尽管对于它是否存在我们有不同的信心。
3．一个病人表现出x、y、z三种症状，有很多疾病都会表现出这三种症状，但病人只患了一种病。不同的医生对于到底是何种疾病导致了这些症状可能会有稍微不同的看法。
对我们而言，将对一个事件发生的信心等同于概率是十分自然的，这正是我们长期以来同世界打交道的方式。我们只能了解到部分的真相，但可以通过不断收集证据来完善我们对事物的观念。与此相对的是，你需要通过训练才能够以频率派的方式思考事物。
为了和传统的概率术语对齐，我们把对一个事件A发生的信念记为P(A)，这个值我们称为先验概率。
伟大的经济学家和思想者John Maynard Keynes曾经说过（也有可能是杜撰的）：“当事实改变，我的观念也跟着改变，你呢？”这句名言反映了贝叶斯派思考事物的方式，即随着证据而更新信念。甚至是，即便证据和初始的信念相反，也不能忽视了证据。我们用P(A|X)表示更新后的信念，其含义为在得到证据X后，A事件的概率。为了和先验概率相对，我们称更新后的信念为后验概率。考虑在观察到证据X后，以下例子的后验概率。
1．P(A)：硬币有50%的几率是正面。P(A|X)：你观察到硬币落地后是正面，把这个观察到的信息记为X，那么现在你会赋100%的概率给正面，0%的概率给负面。
2．P(A)：一段很长很复杂的代码可能含有bug。P(A|X):代码通过了所有的X个测试；现在代码可能仍然有bug，不过这个概率现在变得非常小了。
3．P(A)：病人可能有多种疾病。P(A|X)：做了血液测试之后，得到了证据X，排除了之前可能的一些疾病。
在上述例子中，即便获得了新的证据，我们也并没有完全地放弃初始的信念，但我们重新调整了信念使之更符合目前的证据（也就是说，证据让我们对某些结果更有信心）。
通过引入先验的不确定性，我们事实上允许了我们的初始信念可能是错误的。在观察数据、证据或其他信息之后，我们不断更新我们的信念使得它错得不那么离谱。这和硬币预测正相反，我们通常希望预测得更准确。1.1.2 贝叶斯推断在实践中的运用
如果频率推断和贝叶斯推断是一种编程函数，输入是各种统计问题，那么这两个函数返回的结果可能是不同的。频率推断返回一个估计值（通常是统计量，平均值是一个典型的例子），而贝叶斯推断则会返回概率值。
例如，在代码测试的例子中，如果你问频率函数：“我的代码通过了所有测试，它现在没有bug了吗？”频率函数会给出“yes”的回答。但如果你问贝叶斯函数：“通常我的代码有bug，现在我的代码通过了所有测试，它是不是没有bug了？”贝叶斯函数会给出非常不同的回答，它会给出“yes”和“no”的概率，例如“‘yes’的概率是80%，‘no’的概率是20%。”
这和频率函数返回的结果是非常不同的。注意到贝叶斯函数还有一个额外的信息——“通常的我的代码有bug”，这个参数就是先验信念。把这个参数加进去，贝叶斯函数会将我们的先验概率纳入考虑范围。通常这个参数是可省的，但我们将会发现缺省它会产生什么样的结果。
加入证据 当我们添加更多的证据，初始的信念会不断地被“洗刷”。这是符合期望的，例如如果你的先验是非常荒谬的信念类似“太阳今天会爆炸”，那么你每一天都会被打脸，这时候你希望某种统计推断的方法会纠正初始的信念，或者至少让初始的信念变得不那么荒谬。贝叶斯推断就是你想要的。
让N表示我们拥有的证据的数量。如果N趋于无穷大，那么贝叶斯的结果通常和频率派的结果一致。因此，对于足够大的N，统计推断多多少少都还是比较客观的。另一方面，对于较小的N，统计推断相对而言不稳定，而频率派的结果有更大的方差和置信区间。贝叶斯在这方面要胜出了。通过引入先验概率和返回概率结果（而不是某个固定值），我们保留了不确定性，这种不确定性正是小数据集的不稳定性的反映。
有一种观点认为，对于大的N来说，两种统计技术是无差别的，因为结果类似，并且频率派的计算要更为简单，因而倾向于频率派的方法。在采纳这种观点之前，也许应该先参考Andrew Gelman的论述：“样本从来都不是足够大的。如果N太小不足以进行足够精确的估计，你需要获得更多的数据。但当N“足够大”，你可以开始通过划分数据研究更多的问题（例如在民意调查中，当你已经对全国的民意有了较好的估计，你可以开始分性别、地域、年龄进行更多的统计）。N从来无法做到足够大，因为当它一旦大了，你总是可以开始研究下一个问题从而需要更多的数据。”1.1.3 频率派的模型是错误的吗？
不。频率方法仍然非常有用，在很多领域可能都是最好的办法。例如最小方差线性回归、LASSO回归、EM算法，这些都是非常有效和快速的方法。贝叶斯方法能够在这些方法失效的场景发挥作用，或者是作为更有弹性的模型而补充。1.1.4 关于大数据
出乎意料的是，通常解决大数据预测型问题的方法都是些相对简单的算法。因此，我们会认为大数据预测的难点不在于算法，而在于大规模数据的存储和计算。（让我们再次回顾Gelman的关于样本大小的名言，并且提问：“我真的有大数据吗？”）
中等规模或者更小规模的数据会使得分析变得更为困难。用类似Gelman的话说，如果大数据问题能够被大数据方法简单直接地解决，那么我们应该更关注不那么大的数据集上的问题。1.2 我们的贝叶斯框架
我们感兴趣的估计，可以通过贝叶斯的思想被解释为概率。我们对事件A有一个先验估计——例如，在准备测试之前，我们对代码中的漏洞就有了一个先验的估计。
接下来，观察我们的证据。继续拿代码漏洞为例：如果我们的代码通过了X个测试，我们会相应地调整心里的估计。我们称这个调整过后的新估计为后验概率。调整这个估计值可以通过下面的公式完成，这个公式被称为贝叶斯定理，得名于它的创立者托马斯·贝叶斯。∝P(X|A)P(A) (∝ 的意思是“与之成比例”)
上面的公式并不等同于贝叶斯推论，它是一个存在于贝叶斯推论之外的数学真理。在贝叶斯推论里它仅仅被用来连接先验概率P(A)和更新的后验概率P(A|X)。1.2.1 不得不讲的实例：抛硬币
几乎所有统计书籍都包含一个抛硬币的实例，那我也从这个开始着手吧。假设你不确定在一次抛硬币中得到正面的概率（剧透警告：它是50%），你认为这里肯定是存在某个比例的，称之为p，但是你事先并不清楚p大概会是多少。
我们开始抛硬币，并记录下每一次抛出的结果——正面或反面，这就是我们的观测数据。一个有趣的问题是：“随着收集到越来越多的数据，我们对p的推测是怎么变化的呢？”
说得更具体一些，当面对着很少量的数据或拥有大量数据时，我们的后验概率是怎么样的呢？下面，我们按照观测到的越来越多的数据（抛硬币数据），逐次更新我们的后验概率图。
在图中我们用曲线表示我们的后验概率，曲线越宽，我们的不确定性越大。如图1.2.1所示，当我们刚刚开始观测的时候，我们的后验概率的变化是不稳定的。但是最终，随着观测数据（抛硬币数据）越来越多，这个概率会越来越接近它的真实值p=0.5（图中用虚线标出）。
注意到图中的波峰不一定都出现在0.5那里，当然它也没有必要都这样。应该明白的是我们事前并不知道p会是多少。事实上，如果我们的观测十分的极端，比如说抛了8次只有1次结果是正面的，这种情况我们的分布会离0.5偏差很多（如果缺少先验的知识，当出现8次反面1次正面时，你真的会认为抛硬币结果是公平的吗？）。随着数据的累积，我们可以观察到，虽然不是每个时候都这样，但越来越多地，概率值会出现在p=0.5。
下面这个实例就简单地从数据角度演示一下贝叶斯推断。图1.2.1 后验概率的贝叶斯变换1.2.2 实例：图书管理员还是农民
下面这个故事灵感来自于Daniel Kahneman的《思考，快与慢》一书，史蒂文被描述为一个害羞的人，他乐于助人，但是他对其他人不太关注。他非常乐见事情处于合理的顺序，并对他的工作非常细心。你会认为史蒂文是一个图书管理员还是一个农民呢？从上面的描述来看大多数人都会认为史蒂文看上去更像是图书管理员，但是这里却忽略了一个关于图书管理员和农民的事实：男性图书管理员的人数只有男性农民的1/20。所以从统计学来看史蒂文更有可能是一个农民。
怎么正确地看待这个问题呢？史蒂文实际上更有可能是一个农民还是一个图书管理员呢？把问题简化，假设世上只有两种职业——图书管理员和农民，并且农民的数量确实是图书管理员的20倍。
设事件A为史蒂文是一个图书管理员。如果我们没有史蒂文的任何信息，那么P（A）=1/21=0.047。这是我们的先验。现在假设从史蒂文的邻居们那里我们获得了关于他的一些信息，我们称它们为X。我们想知道的就是P（A|X）。由贝叶斯定理得：
我们知道P（A）是什么意思，那P（X|A）是什么呢？它可以被定义为在史蒂文真的是一个图书管理员的情况下，史蒂文的邻居们给出的某种描述的概率，即如果史蒂文真的是一个图书管理员，他的邻居们将他描述为一个图书管理员的概率。这个值很可能接近于1。假设它为0.95。
P（X）可以解释为：任何人对史蒂文的描述和史蒂文邻居的描述一致的概率。现在这种形式有点难以理解，我们将其做一些逻辑上的改造：P(X)=P(X and A)+P(X and ～A)=P(X|A)P(A)+P(X|～A)P(～A)其中～A表示史蒂文不是一个图书管理员的事件，那么他一定是一个农民。现在我们知道P（X|A）和P(A)，另外也可知P(～A)=1-P(A)=20/21。现在我们只需要知道P(X|～A)，即在史蒂文为一个农民的情况下，史蒂文的邻居们给出的某种描述的概率即可。假设它为0.5，这样，。
结合以上：
这个值并不算高，但是考虑到农民的数量比图书管理员的数量多这么多，这个结果也非常合理了。在图1.2.2中，对比了在史蒂文为农民和史蒂文为图书管理员时的先验和后验概率。%matplotlib inlinefrom IPython.core.pylabtools import figsizeimport numpy as npfrom matplotlib import pyplot as pltfigsize(12.5, 4)plt.rcParams['savefig.dpi'] = 300plt.rcParams['figure.dpi'] = 300colors = ["#348ABD", "#A60628"]prior = [1/21., 20/21.]posterior = [0.087,1-0.087]plt.bar([0, .7], prior, alpha=0.70, width=0.25, color=colors[0], label="prior distribution", lw="3", edgecolor="#348ABD")plt.bar([0+0.25, .7+0.25], posterior, alpha=0.7, width=0.25, color=colors[1], label="posterior distribution", lw="3", edgecolor="#A60628")plt.xticks([0.20, 0.95], ["Librarian", "Farmer"])plt.title("Prior and posterior probabilities of Steve's\ occupation")plt.ylabel("Probability")plt.legend(loc="upper left");
在我们得到X的观测值之后，史蒂文为图书管理员的概率增加了，虽然增加的不是很多，史蒂文为农民的可能性依旧是相当大的。
这是一个关于贝叶斯推断和贝叶斯法则的一个简单的实例。不幸的是，除了在人工结构的情况下，要执行更加复杂的贝叶斯推断所使用到的数学只会变得更加的复杂。在后面我们将看到执行这种复杂的属性分析并没有必要。首先，我们必须扩充我们的建模工具。下一章的概率分布，如果你已经对它很熟悉了，可以选择跳过（或只是浏览一下），但是对于不熟悉的读者，下一章是很有必要的。图1.2.2 史蒂文职业的先验及后验概率1.3 概率分布
首先定义以下希腊文字的发音：α = alphaβ = betaλ = lambdaµ = muσ = sigmaτ = tau
很好。接下来正式开始讲概率分布。首先快速地回忆一下什么是概率分布。设Z为一个随机变量，那么就存在一个跟Z相关的概率分布函数，给定Z任何取值，它都得到一个相应的概率值。
我们把随机变量分为3种类别。● Z为离散的。离散随机变量的取值只能是在特定的列表中。像货
币、电影收视率、投票数都是离散随机变量。当我们将它与连续
型随机变量对比时，这个概念就比较清楚了。● Z为连续的。连续型随机变量的值可以是任意精确数值。像温度、
速度和时间都是连续型变量，因为对于这些数值你可以精确到任
何程度。● Z为混合的。混合型随机变量的值可以为以上两种形式，即结合
了以上两种随机变量的形式。1.3.1 离散情况
如果Z是离散的，那么它的分布为概率质量函数，它度量的是当Z取值为k时的概率，用P(Z=k)表示。可以注意到，概率质量函数完完全全地描述了随机变量Z，即如果知道Z的概率质量方程，那么Z会怎么表现都是可知的。下面介绍一些经常会碰到的概率质量方程，学习它们是十分有必要的。第一个要介绍的概率质量方程十分重要，设Z服从于Poisson分布：
λ 被称为此分布的一个参数，它决定了这个分布的形式。对于Poisson分布来说，λ 可以为任意正数。随着λ 的增大，得到大值的概率会增大；相反地，当λ 减小时，得到小值的概率会增大。λ 可以被称为Poisson分布的强度。
跟λ 可以为任意正数不同，值k可以为任意非负整数，即k必须为0、1、2之类的值。这个是十分重要的，因为如果你要模拟人口分布，你是不可以假设有4.25个或是5.612个人的。
如果一个变量Z存在一个Poisson质量分布，我们可以表示为：Z～Poi(l)
Poisson分布的一个重要性质是：它的期望值等于它的参数。即：E[Z|λ]=λ
这条性质以后经常会被用到，所以记住它是很有用的。在图1.3.1中，展示了不同λ 取值下的概率质量分布。首先值得注意的是，增大λ 的取值，k取大值的概率会增加。其次值得注意的是，虽然x轴在15的时候截止，但是分布并没有截止，它可以延伸到任意非负的整数。figsize(12.5, 4)import scipy.stats as statsa = np.arange(16)poi = stats.poissonlambda_ = [1.5, 4.25]colors = ["#348ABD", "#A60628"]plt.bar(a, poi.pmf(a, lambda_[0]), color=colors[0], label="$\lambda = %.1f$" % lambda_[0], alpha=0.60, edgecolor=colors[0], lw="3")plt.bar(a, poi.pmf(a, lambda_[1]), color=colors[1], label="$\lambda = %.1f$" % lambda_[1], alpha=0.60, edgecolor=colors[1], lw="3")plt.xticks(a + 0.4, a)plt.legend()plt.ylabel("Probability of $k$")plt.xlabel("$k$")plt.title("Probability mass function of a Poisson random variable,\ differing \$\lambda$ values");图1.3.1 不同λ取值情况下，Poisson随机变量的概率质量函数1.3.2 连续情况
对应于离散情况下的概率质量函数，连续情况下概率分布函数被称为概率密度函数。虽然在命名上作这样的区分看起来是没有必要的，但是概率质量函数和概率密度函数有着本质的不同。举一个连续型随机变量的例子：指数密度。指数随机变量的密度函数如下式：-λzf(z|λ)=λe,z≥0Z
类似于Poisson随机变量，指数随机变量只可以取非负值。但是和Poisson分布不同的是，这里的指数可以取任意非负值，包括非整数，例如4.25或是5.612401。这个性质使得计数数据（必须为整数）在这里不太适用；而对于类似时间数据、温度数据（当然是以开氏温标计量）或其他任意对精度有要求的正数数据，它是一种不错的选择。图1.3.2展示了λ取不同值时的概率密度函数。
当随机变量Z拥有参数为λ的指数分布时，我们称Z服从于指数分布，并记为：Z～Exp(λ)
对指定的参数λ，指数型随机变量的期望值为l的逆，即E[Z|λ]=1/λa = np.linspace(0, 4, 100)expo = stats.exponlambda_ = [0.5, 1]for l, c in zip(lambda_, colors): plt.plot(a, expo.pdf(a, scale=1./l), lw=3, color=c, label="$\lambda = %.1f$" % l) plt.fill_between(a, expo.pdf(a, scale=1./l), color=c, alpha=.33)plt.legend()plt.ylabel("Probability density function at $z$")plt.xlabel("$z$")plt.ylim(0,1.2)plt.title("Probability density function of an exponential random\ variable, differing $\lambda$ values");图1.3.2 不同λ取值情况下，指数分布的概率密度函数
这里值得注意的是，概率密度方程在某一点的值并不等于它在这一点的概率。这个将会在后面讲到，当然如果你对它感兴趣，可以加入stackexchange网站上面的讨论。1.3.3 什么是λ
这个问题可以理解为统计的动机是什么。在现实世界，我们并不知道λ的存在，我们只能直观地感受到变量Z，如果确定参数l的话，那就一定要深入到整个事件的背景中去。这个问题其实很难，因为并不存在Z到λ的一一对应关系。对于λ的估计有很多的设计好的方法，但因为λ不是一个可以真实观察到的东西，谁都不能说哪种方式一定是最好的。
贝叶斯推断围绕着对λ取值的估计。与其不断猜测λ的精确取值，不如用一个概率分布来描述λ的可能取值。
起初这看起来或许有些奇怪。毕竟，λ是一个定值，它不一定是随机的！我们怎么可能对一个非随机变量值赋予一个概率呢？不，这样的考虑是老旧的频率派思考方式。如果我们把它们理解为估计值的话，在贝叶斯的哲学体系下，它们是可以被赋予概率的。因此对参数λ估计是完全可以接受的。1.4 使用计算机执行贝叶斯推断
接下来模拟一个有趣的实例，它是一个有关用户发送和收到短信的例子。1.4.1 实例：从短信数据推断行为
你得到了系统中一个用户每天的短信条数数据，如图1.4.1中所示。你很好奇这个用户的短信使用行为是否随着时间有所改变，不管是循序渐进地还是突然地变化。怎么模拟呢？（这实际上是我自己的短信数据。尽情地判断我的受欢迎程度吧。）figsize(12.5, 3.5)count_data = np.loadtxt("data/txtdata.csv")n_count_data = len(count_data)plt.bar(np.arange(n_count_data), count_data, color="#348ABD")plt.xlabel("Time (days)")plt.ylabel("Text messages received")plt.title("Did the user's texting habits change over time?")plt.xlim(0, n_count_data);图1.4.1 用户的短信使用行为是否随着时间发生改变?
在建模之前，仅仅从图1.4.1中你能猜到什么吗？你能说在这一段时间内用户行为出现了什么变化吗？
我们怎么模拟呢？像前文提到的那样，Possion随机变量能很好地模拟这种计数类型的数据。用C表示第i天的短信条数。iC～Poi(λ)i
我们不能确定参数l的真实取值，然而，在图1.4.1中，整个观察周期的后期收到短信的几率升高了，也可以说，λ在某些时段增加了（在前文中有提到过，当l取大值的时候更容易得到较大的结果值。在这里，也就是说一天收到短信条数比较多的概率增大了）。
怎么用数据表示这种观察呢？假设在观察期的某些天（称它为τ），参数λ的取值突然变得比较大。所以参数λ存在两个取值：在时段τ之前有一个，在其他时段有另外一个。在一些资料中，像这样的一个转变称之为转换点：
如果实际上不存在这样的情况，即λ=λ，那么λ的先验分布应该12是均等的。
对于这些不知道的λ我们充满了兴趣。在贝叶斯推断下，我们需要对不同可能值的λ分配相应的先验概率。对参数λ和λ来说什么才λ2是一个好的先验概率分布呢？前面提到过λ可以取任意正数。像我们前面见到的那样，指数分布对任意正数都存在一个连续密度函数，这或许对模拟λ来说是一个很好的选择。但也像前文提到的那样，指数i分布也有它自己对应的参数，所以这个参数也需要包括在我们的模型里面。称它为参数α。λ～Exp(α)1λ～Exp(α)2
α被称为超参数或者是父变量。按字面意思理解，它是一个对其他参数有影响的参数。按照我们最初的设想，α应该对模型的影响不会太大，所以可以对它进行一些灵活的设定。在我们的模型中，我们不希望对这个参数赋予太多的主观色彩。但这里建议将其设定为样本中计算平均值的逆。为什么这样做呢？既然我们用指数分布模拟参数λ，那这里就可以使用指数函数的期望值公式得到：
使用这个值，我们是比较客观的，这样做的话可以减少超参数对模拟造成的影响。另外，我也非常建议大家对上面提到的不同时段的两个λ使用不同的先验。构建两个不同α值的指数分布反映出我们的i先验估计，即在整个观测过程中，收到短信的条数出现了变化。
对于参数τ，由于受到噪声数据的影响，很难为它挑选合适的先验。我们假定每一天的先验估计都是一致的。用公式表达如下：τ～DiscreteUniform(1,70)P(τ=k)=1/70
做了这么多了，那么未知变量的整体先验分布情况应该是什么样的呢？老实说，这并不重要。我们要知道的是，它并不好看，包括很多只有数学家才会喜欢的符号。而且我们的模型会因此变得更加复杂。不管这些啦，毕竟我们关注的只是先验分布而已。
下面会介绍PyMC，它是一个由数学奇才们建立起来的关于贝叶斯分析的Python库。1.4.2 介绍我们的第一板斧：PyMC
PyMC是一个做贝叶斯分析使用的Python库。它是一个运行速度快、维护得很好的库。它唯一的缺点是，它的说明文档在某些领域有所缺失，尤其是在一些能区分菜鸟和高手的领域。本书的主要目标就是解决问题，并展示PyMC库有多酷。
下面用PyMC模拟上面的问题。这种类型的编程被称之为概率编程，对此的误称包括随机产生代码，而且这个名字容易使得使用者误解或者让他们对这个领域产生恐惧。代码当然不是随机的，之所以名字里面包含概率是因为使用编译变量作为模型的组件创建了概率模型。在PyMC中模型组件为第一类原语。
Cronin对概率编程有一段激动人心的描述：“换一种方式考虑这件事情，跟传统的编程仅仅向前运行不同的是，概率编程既可以向前也可以向后运行。它通过向前运行来计算其中包含的所有关于世界的假设结果（例如，它对模型空间的描述），但它也从数据中向后运行，以约束可能的解释。在实践中，许多概率编程系统将这些向前和向后的操作巧妙地交织在一起，以给出有效的最佳的解释。
由于“概率编程”一词会产生很多不必要的困惑，我会克制自己使用它。相反，我会简单地使用“编程”，因为它实际上就是编程。
PyMC代码很容易阅读。唯一的新东西应该是语法，我会截取代码来解释各个部分。只要记住我们模型的组件（τ，λ，λ）为变12量：import pymc as pmalpha = 1.0/count_data.mean()# Recall that count_data is the # variable that holds our text counts.lambda_1 = pm.Exponential("lambda_1", alpha)lambda_2 = pm.Exponential("lambda_2", alpha)tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data)
在这段代码中，我们产生了对应于参数λ和λ的PyMC变量，并12令他们为PyMC中的随机变量，之所以这样称呼它们是因为它们是由后台的随机数产生器生成的。为了表示这个过程，我们称它们由random（）方法构建。在整个训练阶段，我们会发现更好的tau值。print "Random output:", tau.random(), tau.random(), tau.random()[Output]:Random output: 53 21 42@pm.deterministicdef lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2): out = np.zeros(n_count_data) # number of data points out[:tau] = lambda_1 # lambda before tau is lambda_1 out[tau:] = lambda_2 # lambda after (and including) tau is # lambda_2 return out
这段代码产生了一个新的函数lambda_，但事实上我们可以把它想象成为一个随机变量——之前的随机参数l。注意，因为lambda_1、lambda_2、tau是随机的，所以lambda_也会是随机的。目前我们还没有计算出任何变量。
@pm.deterministic是一个标识符，它可以告诉PyMC这是一个确定性函数，即如果函数的输入为确定的话（当然这里它们不是），那函数的结果也是确定的。observation = pm.Poisson("obs", lambda_, value=count_data, observed=True)model = pm.Model([observation, lambda_1, lambda_2, tau])
变量observation包含我们的数据count_data，它是由变量lambda_用我们的数据产生机制得到。我们将observed设定为True来告诉PyMC这在我们的分析中是一个定值。最后，PyMC 希望我们收集所有变量信息并从中产生一个Model实例。当我们拿到结果时就会比较好处理了。
下面的代码将在第3章中解释，但在这里我们展示我们的结果是从哪里来的。可以把它想象成为一个不断学习的过程。这里使用的理论称为马尔科夫链蒙特卡洛（MCMC），在第3章中会给出进一步的解释。利用它可以得到参数λ、λ和τ后验分布中随机变量的阈值。12我们对这些随机变量作直方图，观测他们的后验分布。接下来，将样本（在MCMC中我们称之为迹）放入直方图中。结果如图1.4.2所示。# Mysterious code to be explained in Chapter 3. Suffice it to say,# we will get# 30,000 (40,000 minus 10,000) samples back.mcmc = pm.MCMC(model)mcmc.sample(40000, 10000)[Output]:[-----------------100%-----------------] 40000 of 40000 complete in 9.6 seclambda_1_samples = mcmc.trace('lambda_1')[:]lambda_2_samples = mcmc.trace('lambda_2')[:]tau_samples = mcmc.trace('tau')[:]figsize(14.5, 10)# histogram of the samplesax = plt.subplot(311)ax.set_autoscaley_on(False)plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85, label="posterior of $\lambda_1$", color="#A60628", normed=True)plt.legend(loc="upper left")plt.title(r"""Posterior distributions of the parameters\ $\lambda_1,\;\lambda_2,\;\tau$""")plt.xlim([15, 30])plt.xlabel("$\lambda_1$ value")plt.ylabel("Density")ax = plt.subplot(312)ax.set_autoscaley_on(False)plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.85, label="posterior of $\lambda_2$", color="#7A68A6", normed=True)plt.legend(loc="upper left")plt.xlim([15, 30])plt.xlabel("$\lambda_2$ value")plt.ylabel("Density")plt.subplot(313)w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples)plt.hist(tau_samples, bins=n_count_data, alpha=1, label=r"posterior of $\tau$", color="#467821", weights=w, rwidth=2.)plt.xticks(np.arange(n_count_data))plt.legend(loc="upper left")plt.ylim([0, .75])plt.xlim([35, len(count_data)-20])plt.xlabel(r"$\tau$ (in days)")plt.ylabel("Probability");图1.4.2 参数 λ、λ、τ的后验分布121.4.3 说明
回想一下，贝叶斯方法返回一个分布。因此，我们现在有分布描述未知的λ和τ。我们得到了什么？马上，我们可以看我们估计的不确定性：分布越广，我们的后验概率越小。我们也可以看到参数的合理值：λ大概为18，λ大概是23。这两个λ的后验分布明显是不同的，12这表明用户接收短信的行为确实发生了变化。（请参阅1.6补充说明中的正式参数。）
你还能做哪些其他的观测呢？再看看原始数据，你是否觉得这些结果合理呢？
还要注意到λ的后验分布并不像是指数分布，事实上，后验分布并不是我们从原始模型中可以辨别的任何分布。但这挺好的！这是用计算机处理的一个好处。如果我们不这么做而改用数学方式处理这个问题，将会非常的棘手和混乱。使用计算数学的方式可以让我们不用在方便数学处理上考虑太多。
我们的分析页返回了τ的一个分布。由于它是一个离散变量，所以它的后验看起来和其他两个参数有点不同，它不存在概率区间。我们可以看到，在45天左右，有50%的把握可以说用户的行为是有所改变的。没有发生变化，或者随着时间有了慢慢的变化，τ的后验分布会更加的分散，这也反映出在很多天中τ是不太好确定的。相比之下，从真实的结果中看到，仅仅有三到四天可以认为是潜在的转折点。1.4.4 后验样本到底有什么用？
在这本书的其余部分，我们会面对这样一个问题，它是一个能带领我们得到强大结果的说明。现在，用另外一个实例结束这一章。
我们会用后验样本回答下面的问题：在第t（0≤t≤70）天中，期望收到多少条短信呢？Poisson分布的期望值等于它的参数λ。因此问题相当于：在时间t中，参数λ的期望值是多少。
在下面的代码中，令i指示后验分布中的变量。给定某天t，我们对所有可能的λ求均值，如果 t < τ (也就是说，如果并没有发生什么i变化)，令λ=λ，否则我们令λ = λ。i1，ii2，ifigsize(12.5, 5)# tau_samples, lambda_1_samples, lambda_2_samples contain# N samples from the corresponding posterior distribution.N = tau_samples.shape[0]expected_texts_per_day = np.zeros(n_count_data) # number of data pointsfor day in range(0, n_count_data):# ix is a bool index of all tau samples corresponding to# the switchpoint occurring prior to value of "day."ix = day < tau_samples# Each posterior sample corresponds to a value for tau.# For each day, that value of tau indicates whether we're# "before"# (in the lambda 1 "regime") or# "after" (in the lambda 2 "regime") the switchpoint.# By taking the posterior sample of lambda 1/2 accordingly,# we can average# over all samples to get an expected value for lambda on that day.# As explained, the "message count" random variable is# Poisson-distributed,# and therefore lambda (the Poisson parameter) is the expected# value of# "message count."expected_texts_per_day[day] = (lambda_1_samples[ix].sum()\ + lambda_2_samples[˜ix].sum()) / Nplt.plot(range(n_count_data), expected_texts_per_day, lw=4, color="#E24A33", label="Expected number of text messages\ received")plt.xlim(0, n_count_data)plt.xlabel("Day")plt.ylabel("Number of text messages")plt.title("Number of text messages received versus expected number\ received")plt.ylim(0, 60)plt.bar(np.arange(len(count_data)), count_data, color="#348ABD", alpha=0.65, label="Observed text messages per day")plt.legend(loc="upper left")print expected_texts_per_day[Output]:[ 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7707 17.7708 17.7708 17.7712 17.7717 17.7722 17.7726 17.7767 17.9207 18.4265 20.1932 22.7116 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117 22.7117]
在图1.4.3中展示的结果，很明显地说明了转折点的重要影响。但是对此我们应该保持谨慎的态度，从这个观点看，这里并不存在我们非常希望看到的“短信的期望数量”这样的一条直线。我们的分析结果非常符合之前的估计——用户行为确实发生了改变（如不然，λ1和λ的值应该比较接近），而且这是一个突然的变化，而非一种循序2渐进的变化（τ的先验分布突然出现了峰值）。我们可以推测这种情况产生的原因是：短信费用的降低，天气提醒短信的订阅，或者是一段新的感情。
试读结束[说明：试读内容隐藏了图片]