尽管你已是一个编程老手，但bug仍有可能在代码中存在。于是，在实现了一段特别难的算法之后，你决定先来一个简单的测试用例。这个用例通过了。接着你用了一个稍微复杂的测试用例。再次通过了。接下来更难的测试用例也通过了。这时，你开始觉得也许这段代码已经没有bug了。
如果你这样想，那么恭喜你：你已经在用贝叶斯的方式思考！简单地说，贝叶斯推断是通过新得到的证据不断地更新你的信念。贝叶斯推断很少会做出绝对的判断，但可以做出非常可信的判断。在上面的例子中，我们永远无法100%肯定我们的代码是无缺陷的，除非我们测试每一种可能出现的情形，这在实践中几乎不可能。但是，我们可以对代码进行大量的测试，如果每一次测试都通过了，我们更有把握觉得这段代码是没问题的。贝叶斯推断的工作方式就在这里：我们会随着新的证据不断更新之前的信念，但很少做出绝对的判断，除非所有其他的可能都被一一排除。
和更传统的统计推断不同，贝叶斯推断会保留不确定性。乍听起来，这像一门糟糕的统计方法，难道不是所有的统计都是期望从随机性里推断出确定性吗？要协调这些不一致，我们首先需要像贝叶斯派一样思考。 在贝叶斯派的世界观中，概率是被解释为我们对一件事情发生的相信程度，换句话说，这表明了我们对此事件发生的信心。事实上，我们一会儿就会发现，这就是概率的自然的解释。 为了更清晰地论述，让我们看看来自频率派关于概率的另一种解释。频率派是更古典的统计学派，他们认为概率是事件在长时间内发生的频率。例如，在频率派的哲学语境里，飞机事故的概率指的是长期来看，飞机事故的频率值。对许多事件来说，这样解释概率是有逻辑的，但对某些没有长期频率的事件来说，这样解释是难以理解的。试想一下：在总统选举时，我们经常提及某某候选人获选的概率，但选举本身只发生一次！频率论为了解决这个问题，提出了“替代现实”的说法，从而用在所有这些替代的“现实”中发生的频率定义了这个概率。
贝叶斯派，在另一方面，有更直观的方式。它把概率解释成是对事件发生的信心。简单地说，概率是观点的概述。某人把概率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，排除了之前可能的一些疾病。
在上述例子中，即便获得了新的证据，我们也并没有完全地放弃初始的信念，但我们重新调整了信念使之更符合目前的证据（也就是说，证据让我们对某些结果更有信心）。
通过引入先验的不确定性，我们事实上允许了我们的初始信念可能是错误的。在观察数据、证据或其他信息之后，我们不断更新我们的信念使得它错得不那么离谱。这和硬币预测正相反，我们通常希望预测得更准确。
如果频率推断和贝叶斯推断是一种编程函数，输入是各种统计问题，那么这两个函数返回的结果可能是不同的。频率推断返回一个估计值（通常是统计量，平均值是一个典型的例子），而贝叶斯推断则会返回概率值。
例如，在代码测试的例子中，如果你问频率函数：“我的代码通过了所有测试，它现在没有bug了吗？”频率函数会给出“yes”的回答。但如果你问贝叶斯函数：“通常我的代码有bug，现在我的代码通过了所有测试，它是不是没有bug了？”贝叶斯函数会给出非常不同的回答，它会给出“yes”和“no”的概率，例如“‘yes’的概率是80%，‘no’的概率是20%。”
这和频率函数返回的结果是非常不同的。注意到贝叶斯函数还有一个额外的信息——“通常的我的代码有bug”，这个参数就是先验信念。把这个参数加进去，贝叶斯函数会将我们的先验概率纳入考虑范围。通常这个参数是可省的，但我们将会发现缺省它会产生什么样的结果。
加入证据 当我们添加更多的证据，初始的信念会不断地被“洗刷”。这是符合期望的，例如如果你的先验是非常荒谬的信念类似“太阳今天会爆炸”，那么你每一天都会被打脸，这时候你希望某种统计推断的方法会纠正初始的信念，或者至少让初始的信念变得不那么荒谬。贝叶斯推断就是你想要的。
让N表示我们拥有的证据的数量。如果N趋于无穷大，那么贝叶斯的结果通常和频率派的结果一致。因此，对于足够大的N，统计推断多多少少都还是比较客观的。另一方面，对于较小的N，统计推断相对而言不稳定，而频率派的结果有更大的方差和置信区间。贝叶斯在这方面要胜出了。通过引入先验概率和返回概率结果（而不是某个固定值），我们保留了不确定性，这种不确定性正是小数据集的不稳定性的反映。
有一种观点认为，对于大的N来说，两种统计技术是无差别的，因为结果类似，并且频率派的计算要更为简单，因而倾向于频率派的方法。在采纳这种观点之前，也许应该先参考Andrew Gelman的论述：
“样本从来都不是足够大的。如果N太小不足以进行足够精确的估计，你需要获得更多的数据。但当N“足够大”，你可以开始通过划分数据研究更多的问题（例如在民意调查中，当你已经对全国的民意有了较好的估计，你可以开始分性别、地域、年龄进行更多的统计）。N从来无法做到足够大，因为当它一旦大了，你总是可以开始研究下一个问题从而需要更多的数据。” （参见http://andrewgelman.com/2005/07/31/n_is_never_largl.）
不。频率方法仍然非常有用，在很多领域可能都是最好的办法。例如最小方差线性回归、LASSO回归、EM算法，这些都是非常有效和快速的方法。贝叶斯方法能够在这些方法失效的场景发挥作用，或者是作为更有弹性的模型而补充。
出乎意料的是，通常解决大数据预测型问题的方法都是些相对简单的算法。因此，我们会认为大数据预测的难点不在于算法，而在于大规模数据的存储和计算。（让我们再次回顾Gelman的关于样本大小的名言，并且提问：“我真的有大数据吗？”）
中等规模或者更小规模的数据会使得分析变得更为困难。用类似Gelman的话说，如果大数据问题能够被大数据方法简单直接地解决，那么我们应该更关注不那么大的数据集上的问题。
我们感兴趣的估计，可以通过贝叶斯的思想被解释为概率。我们对事件A有一个先验估计——例如，在准备测试之前，我们对代码中的漏洞就有了一个先验的估计。
接下来，观察我们的证据。继续拿代码漏洞为例：如果我们的代码通过了X个测试，我们会相应地调整心里的估计。我们称这个调整过后的新估计为后验概率。调整这个估计值可以通过下面的公式完成，这个公式被称为贝叶斯定理，得名于它的创立者托马斯·贝叶斯。
∝P(X|A)P(A) (∝ 的意思是“与之成比例”)
上面的公式并不等同于贝叶斯推论，它是一个存在于贝叶斯推论之外的数学真理。在贝叶斯推论里它仅仅被用来连接先验概率P(A)和更新的后验概率P(A|X)。
几乎所有统计书籍都包含一个抛硬币的实例，那我也从这个开始着手吧。假设你不确定在一次抛硬币中得到正面的概率（剧透警告：它是50%），你认为这里肯定是存在某个比例的，称之为p，但是你事先并不清楚p大概会是多少。
我们开始抛硬币，并记录下每一次抛出的结果——正面或反面，这就是我们的观测数据。一个有趣的问题是：“随着收集到越来越多的数据，我们对p的推测是怎么变化的呢？”
说得更具体一些，当面对着很少量的数据或拥有大量数据时，我们的后验概率是怎么样的呢？下面，我们按照观测到的越来越多的数据（抛硬币数据），逐次更新我们的后验概率图。
在图中我们用曲线表示我们的后验概率，曲线越宽，我们的不确定性越大。如图2.1所示，当我们刚刚开始观测的时候，我们的后验概率的变化是不稳定的。但是最终，随着观测数据（抛硬币数据）越来越多，这个概率会越来越接近它的真实值p=0.5（图中用虚线标出）。
注意到图中的波峰不一定都出现在0.5那里，当然它也没有必要都这样。应该明白的是我们事前并不知道p会是多少。事实上，如果我们的观测十分的极端，比如说抛了8次只有1次结果是正面的，这种情况我们的分布会离0.5偏差很多（如果缺少先验的知识，当出现8次反面1次正面时，你真的会认为抛硬币结果是公平的吗？）。随着数据的累积，我们可以观察到，虽然不是每个时候都这样，但越来越多地，概率值会出现在p=0.5。
下面这个实例就简单地从数据角度演示一下贝叶斯推断。
图2.1 后验概率的贝叶斯变换
下面这个故事灵感来自于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，这样，
结合以上：
这个值并不算高，但是考虑到农民的数量比图书管理员的数量多这么多，这个结果也非常合理了。在图2.2中，对比了在史蒂文为农民和史蒂文为图书管理员时的先验和后验概率。
%matplotlib inline
from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
figsize(12.5, 4)
plt.rcParams['savefig.dpi'] = 300
plt.rcParams['figure.dpi'] = 300
colors = ["#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的观测值之后，史蒂文为图书管理员的概率增加了，虽然增加的不是很多，史蒂文为农民的可能性依旧是相当大的。
这是一个关于贝叶斯推断和贝叶斯法则的一个简单的实例。不幸的是，除了在人工结构的情况下，要执行更加复杂的贝叶斯推断所使用到的数学只会变得更加的复杂。在后面我们将看到执行这种复杂的属性分析并没有必要。首先，我们必须扩充我们的建模工具。
图2.2 史蒂文职业的先验及后验概率
首先定义以下希腊文字的发音：
α = alpha
β = beta
λ = lambda
µ = mu
σ = sigma
τ = tau
很好。接下来正式开始讲概率分布。首先快速地回忆一下什么是概率分布。设Z为一个随机变量，那么就存在一个跟Z相关的概率分布函数，给定Z任何取值，它都得到一个相应的概率值。
我们把随机变量分为3种类别。
如果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(λ)
Poisson分布的一个重要性质是：它的期望值等于它的参数。即：
E[Z|λ]=λ
这条性质以后经常会被用到，所以记住它是很有用的。在图3.1中，展示了不同λ取值下的概率质量分布。首先值得注意的是，增大λ的取值，k取大值的概率会增加。其次值得注意的是，虽然x轴在15的时候截止，但是分布并没有截止，它可以延伸到任意非负的整数。
figsize(12.5, 4)
import scipy.stats as stats
a = np.arange(16)
poi = stats.poisson
lambda_ = [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");
图3.1 不同λ取值情况下，Poisson随机变量的概率质量函数
对应于离散情况下的概率质量函数，连续情况下概率分布函数被称为概率密度函数。虽然在命名上作这样的区分看起来是没有必要的，但是概率质量函数和概率密度函数有着本质的不同。举一个连续型随机变量的例子：指数密度。指数随机变量的密度函数如下式：
fZ(z|λ)=λe-λz,z≥0
类似于Poisson随机变量，指数随机变量只可以取非负值。但是和Poisson分布不同的是，这里的指数可以取任意非负值，包括非整数，例如4.25或是5.612401。这个性质使得计数数据（必须为整数）在这里不太适用；而对于类似时间数据、温度数据（当然是以开氏温标计量）或其他任意对精度有要求的正数数据，它是一种不错的选择。图3.2展示了λ取不同值时的概率密度函数。
当随机变量Z拥有参数为λ的指数分布时，我们称Z服从于指数分布，并记为：
Z～Exp(λ)
对指定的参数λ，指数型随机变量的期望值为λ的逆，即
E[Z|λ]=1/λ
a = np.linspace(0, 4, 100)
expo = stats.expon
lambda_ = [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");
图3.2 不同λ取值情况下，指数分布的概率密度函数
这里值得注意的是，概率密度方程在某一点的值并不等于它在这一点的概率。这个将会在后面讲到，当然如果你对它感兴趣，可以加入下面的讨论：http://stats.stackexchange.com/questions/4220/a-probability-distribution-value-exceeding-1-is-ok。
这个问题可以理解为统计的动机是什么。在现实世界，我们并不知道λ的存在，我们只能直观地感受到变量Z，如果确定参数λ的话，那就一定要深入到整个事件的背景中去。这个问题其实很难，因为并不存在Z到λ的一一对应关系。对于λ的估计有很多的设计好的方法，但因为λ不是一个可以真实观察到的东西，谁都不能说哪种方式一定是最好的。
贝叶斯推断围绕着对λ取值的估计。与其不断猜测λ的精确取值，不如用一个概率分布来描述λ的可能取值。
起初这看起来或许有些奇怪。毕竟，λ是一个定值，它不一定是随机的！我们怎么可能对一个非随机变量值赋予一个概率呢？不，这样的考虑是老旧的频率派思考方式。如果我们把它们理解为估计值的话，在贝叶斯的哲学体系下，它们是可以被赋予概率的。因此对参数λ估计是完全可以接受的。
接下来模拟一个有趣的实例，它是一个有关用户发送和收到短信的例子。
你得到了系统中一个用户每天的短信条数数据，如图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);
图4.1 用户的短信使用行为是否随着时间发生改变?
在建模之前，仅仅从图4.1中你能猜到什么吗？你能说在这一段时间内用户行为出现了什么变化吗？
我们怎么模拟呢？像前文提到的那样，Possion随机变量能很好地模拟这种计数类型的数据。用Ci表示第i天的短信条数。
Ci～Poi(λ)
我们不能确定参数λ的真实取值，然而，在图4.1中，整个观察周期的后期收到短信的几率升高了，也可以说，λ在某些时段增加了（当λ取大值的时候更容易得到较大的结果值。在这里，也就是说一天收到短信条数比较多的概率增大了）。
怎么用数据表示这种观察呢？假设在观察期的某些天（称它为τ），参数λ的取值突然变得比较大。所以参数λ存在两个取值：在时段τ之前有一个，在其他时段有另外一个。在一些资料中，像这样的一个转变称之为转换点：
如果实际上不存在这样的情况，即λ1=λ2，那么λ的先验分布应该是均等的。
对于这些不知道的λ我们充满了兴趣。在贝叶斯推断下，我们需要对不同可能值的λ分配相应的先验概率。对参数λ1和λ2来说什么才是一个好的先验概率分布呢？前面提到过λ可以取任意正数。像我们前面见到的那样，指数分布对任意正数都存在一个连续密度函数，这或许对模拟λi来说是一个很好的选择。但也像前文提到的那样，指数分布也有它自己对应的参数，所以这个参数也需要包括在我们的模型里面。称它为参数α。
λ1～Exp(α)
λ2～Exp(α)
α被称为超参数或者是父变量。按字面意思理解，它是一个对其他参数有影响的参数。按照我们最初的设想，α应该对模型的影响不会太大，所以可以对它进行一些灵活的设定。在我们的模型中，我们不希望对这个参数赋予太多的主观色彩。但这里建议将其设定为样本中计算平均值的逆。为什么这样做呢？既然我们用指数分布模拟参数λ，那这里就可以使用指数函数的期望值公式得到：
使用这个值，我们是比较客观的，这样做的话可以减少超参数对模拟造成的影响。另外，我也非常建议大家对上面提到的不同时段的两个λi使用不同的先验。构建两个不同α值的指数分布反映出我们的先验估计，即在整个观测过程中，收到短信的条数出现了变化。
对于参数τ，由于受到噪声数据的影响，很难为它挑选合适的先验。我们假定每一天的先验估计都是一致的。用公式表达如下：
τ～DiscreteUniform(1,70)
＝﹥P(τ=k)=1/70
做了这么多了，那么未知变量的整体先验分布情况应该是什么样的呢？老实说，这并不重要。我们要知道的是，它并不好看，包括很多只有数学家才会喜欢的符号。而且我们的模型会因此变得更加复杂。不管这些啦，毕竟我们关注的只是先验分布而已。
下面会介绍PyMC，它是一个由数学奇才们建立起来的关于贝叶斯分析的Python库。
PyMC是一个做贝叶斯分析使用的Python库。它是一个运行速度快、维护得很好的库。它唯一的缺点是，它的说明文档在某些领域有所缺失，尤其是在一些能区分菜鸟和高手的领域。本书的主要目标就是解决问题，并展示PyMC库有多酷。
下面用PyMC模拟上面的问题。这种类型的编程被称之为概率编程，对此的误称包括随机产生代码，而且这个名字容易使得使用者误解或者让他们对这个领域产生恐惧。代码当然不是随机的，之所以名字里面包含概率是因为使用编译变量作为模型的组件创建了概率模型。在PyMC中模型组件为第一类原语。
Cronin对概率编程有一段激动人心的描述：
“换一种方式考虑这件事情，跟传统的编程仅仅向前运行不同的是，概率编程既可以向前也可以向后运行。它通过向前运行来计算其中包含的所有关于世界的假设结果（例如，它对模型空间的描述），但它也从数据中向后运行，以约束可能的解释。在实践中，许多概率编程系统将这些向前和向后的操作巧妙地交织在一起，以给出有效的最佳的解释。
由于“概率编程”一词会产生很多不必要的困惑，我会克制自己使用它。相反，我会简单地使用“编程”，因为它实际上就是编程。
PyMC代码很容易阅读。唯一的新东西应该是语法，我会截取代码来解释各个部分。只要记住我们模型的组件（τ，λ1，λ2）为变量：
import pymc as pm
alpha = 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)
在这段代码中，我们产生了对应于参数λ1和λ2的PyMC变量，并令他们为PyMC中的随机变量，之所以这样称呼它们是因为它们是由后台的随机数产生器生成的。为了表示这个过程，我们称它们由random（）方法构建。在整个训练阶段，我们会发现更好的tau值。
print "Random output:", tau.random(), tau.random(), tau.random()
[Output]:
Random output: 53 21 42
@pm.deterministic
def 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，但事实上我们可以把它想象成为一个随机变量——之前的随机参数λ。注意，因为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包含我们的数据countdata，它是由变量lambda用我们的数据产生机制得到。我们将observed设定为True来告诉PyMC这在我们的分析中是一个定值。最后，PyMC 希望我们收集所有变量信息并从中产生一个Model实例。当我们拿到结果时就会比较好处理了。
下面的代码将在第3章中解释，但在这里我们展示我们的结果是从哪里来的。可以把它想象成为一个不断学习的过程。这里使用的理论称为马尔科夫链蒙特卡洛（MCMC），在第3章中会给出进一步的解释。利用它可以得到参数λ1、λ2和τ先验分布中随机变量的阈值。我们对这些随机变量作直方图，观测他们的先验分布。接下来，将样本（在MCMC中我们称之为迹）放入直方图中。结果如图1.4.2所示。
# Mysteri-ous 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 sec
lambda_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 samples
ax = plt.subplot(311)
ax.set_autoscaley_on(False)
plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85,
la-bel="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,
la-bel="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,
la-bel=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");
图4.2 参数 λ1、λ2、τ的后验分布
回想一下，贝叶斯方法返回一个分布。因此，我们现在有分布描述未知的λ和τ。我们得到了什么？马上，我们可以看我们估计的不确定性：分布越广，我们的后验概率越小。我们也可以看到参数的合理值：λ1大概为18，λ2大概是23。这两个λ的后验分布明显是不同的，这表明用户接收短信的行为确实发生了变化。（
你还能做哪些其他的观测呢？再看看原始数据，你是否觉得这些结果合理呢？
还要注意到λ的后验分布并不像是指数分布，事实上，后验分布并不是我们从原始模型中可以辨别的任何分布。但这挺好的！这是用计算机处理的一个好处。如果我们不这么做而改用数学方式处理这个问题，将会非常的棘手和混乱。使用计算数学的方式可以让我们不用在方便数学处理上考虑太多。
我们的分析页返回了τ的一个分布。由于它是一个离散变量，所以它的后验看起来和其他两个参数有点不同，它不存在概率区间。我们可以看到，在45天左右，有50%的把握可以说用户的行为是有所改变的。没有发生变化，或者随着时间有了慢慢的变化，τ的后验分布会更加的分散，这也反映出在很多天中τ是不太好确定的。相比之下，从真实的结果中看到，仅仅有三到四天可以认为是潜在的转折点。
在这本书的其余部分，我们会面对这样一个问题，它是一个能带领我们得到强大结果的说明。现在，用另外一个实例结束本文。
我们会用后验样本回答下面的问题：在第t（0≤t≤70）天中，期望收到多少条短信呢？Poisson分布的期望值等于它的参数λ。因此问题相当于：在时间t中，参数λ的期望值是多少。
在下面的代码中，令i指示后验分布中的变量。给定某天t，我们对所有可能的λ求均值，如果 t < τi (也就是说，如果并没有发生什么变化)，令λi=λ1，i，否则我们令λi = λ2，i。
figsize(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 points
for 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()) / N
plt.plot(range(n_count_data), expected_texts_per_day, lw=4,
col-or="#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",
al-pha=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]
在图4.3中展示的结果，很明显地说明了转折点的重要影响。但是对此我们应该保持谨慎的态度，从这个观点看，这里并不存在我们非常希望看到的“短信的期望数量”这样的一条直线。我们的分析结果非常符合之前的估计——用户行为确实发生了改变（如不然，λ1和λ2的值应该比较接近），而且这是一个突然的变化，而非一种循序渐进的变化（τ的先验分布突然出现了峰值）。我们可以推测这种情况产生的原因是：短信费用的降低，天气提醒短信的订阅，或者是一段新的感情。
图4.3 实际收到的短信量和期望收到的量
这一章介绍了频率派和贝叶斯派对概率的解释的差别。同时我们也学到了两个重要的分布：Poisson分布和指数分布。这是今后我们构建更多贝叶斯模型的两块重要基石，就像我们在短信接收例子中所做的那样。在第2章中，我们会探讨更多的建模和PyMC策略。
在短信接收例子中，我们直观地观测了λ1和λ2的先验信息并认为它们是不同的。这很公平，毕竟先验的位置基本离得非常远。但如果这并不是真相，有一部分分布其实是重合的呢？我们怎么才能让上面的结论更加的正式呢？
一种方法就是计算出P(λ1<λ2|data)，即在获得观察数据的前提下，λ1的真实值比λ2小的概率。如果这个概率接近50%，那相当于抛硬币得到的结果，这样我们就不能确定它们是否真的不同。如果这个概率接近100%，那么我们就能确定地说这两个值是不同的。利用后验中的样本，这种计算非常简单——我们计算λ1后验中的样本比λ2后验中的样本小的次数占比：
print (lambda_1_samples < lambda_2_samples)
# Boolean array: True if lambda_1 is less than lambda_2.
[Output]:
[ True True True True ..., True True True True]
# How often does this happen?
print (lambda_1_samples < lambda_2_samples).sum()
# How many samples are there?
print lambda_1_samples.shape[0]
[Output]:
29994
30000
# The ratio is the probability. Or, we can just use .mean:
print (lambda_1_samples < lambda_2_samples).mean()
[Output]:
0.9998
结果很显然，有几乎100%的把握可以说这两个值是不等的。
我们也可以再问详细一点，比如：“两个值之间相差1、2、5、10的概率有多大？”
# The vector abs(lambda_1_samples - lambda_2_samples) > 1 is a boolean,
# True if the values are more than 1 apart, False otherwise.
# How often does this happen? Use .mean()
for d in [1,2,5,10]:
v = (abs(lambda_1_samples - lambda_2_samples) >= d).mean()
print "What is the probability the difference is larger than %d\
? %.2f"%(d,v)
[Output]:
What is the probability the difference is larger than 1? 1.00
What is the probability the difference is larger than 2? 1.00
What is the probability the difference is larger than 5? 0.49
What is the probability the difference is larger than 10? 0.00
读者们或许会对前面模型中转折点个数的扩充，即如果不止一个转折点会怎么样感兴趣，或者会对只有一个转折点的结论表示怀疑。下面我们把模型扩充至两个转折点（意味着会出现3个λi）。新模型跟之前的比较相像。
其中
λ1～Exp(α)
λ2～Exp(α)
λ3～Exp(α)
并且
τ1～DiscreteUniform(1,69)
τ2～DiscreteUniform(τ1,70)
我们把这个模型也编译成代码，跟前面的代码看上去差不多。
lambda_1 = pm.Exponential("lambda_1", alpha)
lambda_2 = pm.Exponential("lambda_2", alpha)
lambda_3 = pm.Exponential("lambda_3", alpha)
tau_1 = pm.DiscreteUniform("tau_1", lower=0, upper=n_count_data-1)
tau_2 = pm.DiscreteUniform("tau_2", lower=tau_1, upper=n_count_data)
@pm.deterministic
def lambda_(tau_1=tau_1, tau_2=tau_2,
lamb-da_1=lambda_1, lambda_2=lambda_2, lambda_3 = lambda_3):
out = np.zeros(n_count_data) # number of data points
out[:tau_1] = lambda_1 # lambda before tau is lambda_1
out[tau_1:tau_2] = lambda_2
out[tau_2:] = lambda_3 # lambda after (and including) tau
# is lambda_2
return out
observa-tion = pm.Poisson("obs", lambda_, value=count_data,observed=True)
mod-el = pm.Model([observation, lambda_1, lambda_2, lambda_3, tau_1,
tau_2])
mcmc = pm.MCMC(model)
mcmc.sample(40000, 10000)
图6.1展示了5个未知数的后验。我们可以看到模型的转折点大致在第45天和第47天的时候取得。对此你怎么认为呢？我们的模型是否对数据过拟合呢？
确实，我们都可能对数据中有多少个转折点抱有疑惑的态度。例如，我就认为一个转折点好过两个转折点，两个转折点好过三个转折点，以此类推。这意味着对于应该有多少个转折点可以设置一个先验分布并让模型自己做决定！在对模型进行调整之后，答案是肯定的，一个转折点确实比较适合。代码在本章就不展示了，这里我只是想介绍一种思想：用怀疑数据那样的眼光审视我们的模型。
图6.1 扩充后的短信模型中5个未知参数的后验分布
1．利用 lambda_1_samples 和 lambda_2_samples，怎么获得参数λ1和λ2后验分布的平均值 ?
2．计算短信频率提升比例的期望值是多少。提示：需要计算 (lambda_2_samples-lambda_1_samples)/lambda_1_samples的均值。注意这个结果和(lambda_2_ samples.mean()-lambda_1_samples.mean())/ lambda_1_samples.mean()计算出来的结果是有很大区别的。
3．在τ小于45的前提下，计算λ1的均值。也就是说，在我们被告知行为的变化发生在第45天之前时，对λ1的期望会是多少? （不需要重复PyMC 那部分，只需要考虑当 tau_samples < 45时所有可能的情况。）
1．计算后验的均值（即后验的期望值），我们只需要用到样本和a.mean函数。
print lambda_1_samples.mean()
print lambda_2_samples.mean()
2．给定两个数a 和 b，相对增长可以由 (a − b)/b给出。在我们的实例中，我们并不能确定λ1和λ2的值是多少。通过计算
(lambda_2_samples-lambda_1_samples)/lambda_1_samples
我们得到另外一个向量，它表示相对增长的后验，如图7.1所示。
relative_increase_samples = (lambda_2_samples-lambda_1_samples)
/lambda_1_samples
print relative_increase_samples
[Output]:
[ 0.263 0.263 0.263 0.263 ..., 0.1622 0.1898 0.1883 0.1883]
figsize(12.5,4)
plt.hist(relative_increase_samples, histtype='stepfilled',
bins=30, alpha=0.85, color="#7A68A6", normed=True,
label='posterior of relative increase')
plt.xlabel("Relative increase")
plt.ylabel("Density of relative increase")
plt.title("Posterior of relative increase")
plt.legend();
为了计算这个均值，需要用到新向量的均值：
print relative_increase_samples.mean()
[Output]:
0.280845247899
图7.1 相对增长的后验
3．如果已知 τ < 45，那么所有样本都需要考虑到这点：
ix = tau_samples < 45
print lambda_1_samples[ix].mean()
[Output]:
17.7484086925
内容简介
本书基于PyMC语言以及一系列常用的Python数据分析框架，如NumPy、SciPy和Matplotlib，通过概率编程的方式，讲解了贝叶斯推断的原理和实现方法。该方法常常可以在避免引入大量数学分析的前提下，有效地解决问题。书中使用的案例往往是工作中遇到的实际问题，有趣并且实用。作者的阐述也尽量避免冗长的数学分析，而让读者可以动手解决一个个的具体问题。通过对本书的学习，读者可以对贝叶斯思维、概率编程有较为深入的了解，为将来从事机器学习、数据分析相关的工作打下基础。本书适用于机器学习、贝叶斯推断、概率编程等相关领域的从业者和爱好者，也适合普通开发人员了解贝叶斯统计而使用。