程序员如何用贝叶斯方式思考

1 引言

尽管你已是一个编程老手，但bug仍有可能在代码中存在。于是，在实现了一段特别难的算法之后，你决定先来一个简单的测试用例。这个用例通过了。接着你用了一个稍微复杂的测试用例。再次通过了。接下来更难的测试用例也通过了。这时，你开始觉得也许这段代码已经没有bug了。

如果你这样想，那么恭喜你：你已经在用贝叶斯的方式思考！简单地说，贝叶斯推断是通过新得到的证据不断地更新你的信念。贝叶斯推断很少会做出绝对的判断，但可以做出非常可信的判断。在上面的例子中，我们永远无法100%肯定我们的代码是无缺陷的，除非我们测试每一种可能出现的情形，这在实践中几乎不可能。但是，我们可以对代码进行大量的测试，如果每一次测试都通过了，我们更有把握觉得这段代码是没问题的。贝叶斯推断的工作方式就在这里：我们会随着新的证据不断更新之前的信念，但很少做出绝对的判断，除非所有其他的可能都被一一排除。

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.2 贝叶斯推断在实践中的运用

如果频率推断和贝叶斯推断是一种编程函数，输入是各种统计问题，那么这两个函数返回的结果可能是不同的。频率推断返回一个估计值（通常是统计量，平均值是一个典型的例子），而贝叶斯推断则会返回概率值。

例如，在代码测试的例子中，如果你问频率函数：“我的代码通过了所有测试，它现在没有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.）

1.3　频率派的模型是错误的吗？

不。频率方法仍然非常有用，在很多领域可能都是最好的办法。例如最小方差线性回归、LASSO回归、EM算法，这些都是非常有效和快速的方法。贝叶斯方法能够在这些方法失效的场景发挥作用，或者是作为更有弹性的模型而补充。

1.4　关于大数据

出乎意料的是，通常解决大数据预测型问题的方法都是些相对简单的算法。因此，我们会认为大数据预测的难点不在于算法，而在于大规模数据的存储和计算。（让我们再次回顾Gelman的关于样本大小的名言，并且提问：“我真的有大数据吗？”）

中等规模或者更小规模的数据会使得分析变得更为困难。用类似Gelman的话说，如果大数据问题能够被大数据方法简单直接地解决，那么我们应该更关注不那么大的数据集上的问题。

2 我们的贝叶斯框架

我们感兴趣的估计，可以通过贝叶斯的思想被解释为概率。我们对事件A有一个先验估计——例如，在准备测试之前，我们对代码中的漏洞就有了一个先验的估计。

接下来，观察我们的证据。继续拿代码漏洞为例：如果我们的代码通过了X个测试，我们会相应地调整心里的估计。我们称这个调整过后的新估计为后验概率。调整这个估计值可以通过下面的公式完成，这个公式被称为贝叶斯定理，得名于它的创立者托马斯·贝叶斯。

∝P(X|A)P(A) (∝ 的意思是“与之成比例”)

上面的公式并不等同于贝叶斯推论，它是一个存在于贝叶斯推论之外的数学真理。在贝叶斯推论里它仅仅被用来连接先验概率P(A)和更新的后验概率P(A|X)。

2.1　不得不讲的实例：抛硬币

几乎所有统计书籍都包含一个抛硬币的实例，那我也从这个开始着手吧。假设你不确定在一次抛硬币中得到正面的概率（剧透警告：它是50%），你认为这里肯定是存在某个比例的，称之为p，但是你事先并不清楚p大概会是多少。

我们开始抛硬币，并记录下每一次抛出的结果——正面或反面，这就是我们的观测数据。一个有趣的问题是：“随着收集到越来越多的数据，我们对p的推测是怎么变化的呢？”

说得更具体一些，当面对着很少量的数据或拥有大量数据时，我们的后验概率是怎么样的呢？下面，我们按照观测到的越来越多的数据（抛硬币数据），逐次更新我们的后验概率图。

在图中我们用曲线表示我们的后验概率，曲线越宽，我们的不确定性越大。如图2.1所示，当我们刚刚开始观测的时候，我们的后验概率的变化是不稳定的。但是最终，随着观测数据（抛硬币数据）越来越多，这个概率会越来越接近它的真实值p=0.5（图中用虚线标出）。

注意到图中的波峰不一定都出现在0.5那里，当然它也没有必要都这样。应该明白的是我们事前并不知道p会是多少。事实上，如果我们的观测十分的极端，比如说抛了8次只有1次结果是正面的，这种情况我们的分布会离0.5偏差很多（如果缺少先验的知识，当出现8次反面1次正面时，你真的会认为抛硬币结果是公平的吗？）。随着数据的累积，我们可以观察到，虽然不是每个时候都这样，但越来越多地，概率值会出现在p=0.5。

下面这个实例就简单地从数据角度演示一下贝叶斯推断。

图2.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，这样，

结合以上：

这个值并不算高，但是考虑到农民的数量比图书管理员的数量多这么多，这个结果也非常合理了。在图2.2中，对比了在史蒂文为农民和史蒂文为图书管理员时的先验和后验概率。

%matplotlib inlinefrom IPython.core.pylabtools import figsizeimport numpy as np
from matplotlib import pyplot as pltfigsize(12.5, 4)plt.rcParams['savefig.dpi'] = 300plt.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　史蒂文职业的先验及后验概率

3　概率分布

首先定义以下希腊文字的发音：

α = alpha

β = beta

λ = lambda

µ = mu

σ = sigma

τ = tau

很好。接下来正式开始讲概率分布。首先快速地回忆一下什么是概率分布。设Z为一个随机变量，那么就存在一个跟Z相关的概率分布函数，给定Z任何取值，它都得到一个相应的概率值。

我们把随机变量分为3种类别。

• Z为离散的。离散随机变量的取值只能是在特定的列表中。像货币、电影收视率、投票数都是离散随机变量。当我们将它与连续型随机变量对比时，这个概念就比较清楚了。

• Z为连续的。连续型随机变量的值可以是任意精确数值。像温度、速度和时间都是连续型变量，因为对于这些数值你可以精确到任何程度。

• Z为混合的。混合型随机变量的值可以为以上两种形式，即结合了以上两种随机变量的形式。

　离散情况

如果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。

3.3　什么是λ

这个问题可以理解为统计的动机是什么。在现实世界，我们并不知道λ的存在，我们只能直观地感受到变量Z，如果确定参数λ的话，那就一定要深入到整个事件的背景中去。这个问题其实很难，因为并不存在Z到λ的一一对应关系。对于λ的估计有很多的设计好的方法，但因为λ不是一个可以真实观察到的东西，谁都不能说哪种方式一定是最好的。

贝叶斯推断围绕着对λ取值的估计。与其不断猜测λ的精确取值，不如用一个概率分布来描述λ的可能取值。

起初这看起来或许有些奇怪。毕竟，λ是一个定值，它不一定是随机的！我们怎么可能对一个非随机变量值赋予一个概率呢？不，这样的考虑是老旧的频率派思考方式。如果我们把它们理解为估计值的话，在贝叶斯的哲学体系下，它们是可以被赋予概率的。因此对参数λ估计是完全可以接受的。

4　使用计算机执行贝叶斯推断

接下来模拟一个有趣的实例，它是一个有关用户发送和收到短信的例子。

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);

图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库。

4.2　介绍我们的第一板斧：PyMC

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,

其中

λ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个未知参数的后验分布

7习题

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时所有可能的情况。）

8答案

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 < 45print lambda_1_samples[ix].mean()

[Output]:
　17.7484086925

本文节选自《贝叶斯方法：概率编程与贝叶斯推断》