你是一名经验丰富的程序员，但是bug仍然暗藏在你的代码中。实现一个极其困难的算法后，你决定在一个简单的例子上测试自己的代码。过了。然后在一个稍稍困难的问题上进行测试，还是过了。接着这样下去，更加复杂的问题，都过了！你开始相信自己的代码莫有问题了～

如果你这样子进行思考，那么祝贺你，你是在如贝叶斯主义者那样进行思考！贝叶斯推断只是简单地在考虑了新的证据后，更新你的信念。贝叶斯主义者很少对于一个结果很肯定，但是他们可以对某件事有一定的信心。就像上面的这个例子，我们不会100%地肯定我们的代码是无bug的除非我们在每个可能的情形下进行测试；这个在现实情形下是不可能的。相反，我们可以在大量的问题上对这代码进行测试，如果测试结果都OK，那么我们可以认为代码本身没有什么问题，对正确性相当地肯定，但仍会有所怀疑。贝叶斯推断的思想是完全一致的：我们更新自己对结果的信念；几乎不会绝对滴相信一件事的发生，除非我们已经排除了所有的可能性。

思维的贝叶斯状态

贝叶斯推断和传统的统计推断的不同在于它保留了不确定性。初看，这是一种不太好的统计技术。难道统计不是关于从随机性中推导出确定性的么？为了弄清楚这一点，我们需要开始像一位贝叶斯学家一样地考虑问题。

贝叶斯世界观下，概率被当作事件的可信度的度量，即我们对一个事件发生有多么的确定。事实上，我们一会儿会看到这就是概率的自然的解释。

为了更加清晰地说明这点，我们考虑另一个概率的解释：频率主义者，这是更加经典的统计观点的持有者，他们假设概率是一个事件的长期运行所获得的频率。例如，在频率主义的观点下，飞机失事的概率就是长期过程中飞机失事的频率。这种观点对于很多的情形的事件的概率是很合理的。然而你想想：我们经常对总统选举的结果赋值概率，但是这样的事件就只会进行一次！频率主义者通过启用另外的现实绕过这个问题，然后对所有这些现实进行汇总，然后出现的频率定义了这个概率。

而贝叶斯学派，有一种更加符合直觉的观点。他们将概率解释成信念（belief）的度量，或者一个事件发生的信心（confidence）。简单说，概率就是一个观点的总结。一个将信念为0赋给了事件，则表示不相信这个事件会发生；相反，如果对该事件赋值为1那么这个人是绝对相信事件会发生。而介于0和1之间的那些值则对应于其他的结果。这样定义方式，在飞机事件的例子中的概率是相符的，已经观察到了飞机事件的频率后，个人的信念应当与那个频率相同，而不包含其他外部的信息。类似的，在这个定义下，概率和信念相同，这样子在总统选举的例子中，我们对于选举结果的概率（信念）的描述就是很合理的：你对候选人A获得胜利有多确信？

注意在上面的例子中，我们把信念（概率）度量赋给了一个个体，而不是给大自然。这是非常有趣的，因为这个定义带来了个体之间冲突的信念的可能。不过这样也是更加合理的：不同的个体对事件的发生有不同的信念，因为他们拥有对这个世界的不一样的信息。这种不同信念的存在并不是说某些人就是错误的。考虑下面的场景，描述了个人信念和概率之间的关系：

• 我丢一个硬币，我们同时来猜测结果。假设硬币是公平的，我们可能都同意出现正面的概率是1/2。接着，假设我偷偷看了一眼。现在我已经知道硬币的正反了：我置概率1.0给正面或者反（这里随便给）。现在问你对硬币正反的猜测的信念？我对结果的知识并不影响硬币的结果。因此，我们对结果赋给了不同的概率。
• 你的代码有一个bug或者没有，但是我们都不能确信哪个是正确的，尽管我们都持有这样的信念——有或者没有bug。
• 一个病人出现了3个症状——x、y和z。有若干种疾病可能会导致所有这些症状，但是仅有一个疾病出现。一名医生对某个疾病持有信念，但是不同的医生却会有略微不同的信念。

将概率看作信念的哲学观点对人类来说是自然的。我们将会不断地使用这个观点，因为我们和世界交互的过程中只能看到部分的真相（partial truth），但是可以搜集更多的证据来形成信念。而对立着的是，我们需要被训练成频率主义者。

为了熟悉传统的概率形式，我们将对事件A的信念表示为P(A)。我们将这个数量称为先验概率。

John Maynard Keynes——一位著名的经济学家和思想家，曾说：“当事实发生变化，我改变自己的思维。你呢，先生？”这段话反映了一名贝叶斯主义者在看到证据后的信念的更新。甚至，尤其是，如果证据和刚开始时候的信念完全相违的时候，这个证据是不可以被忽略的。我们将更新后的信念表示为P(A|X)，解释为在给定证据X后的A事件发生的概率。我们称更新后的这个信念为后验概率以示区别。例如，考虑上面那些例子在观察到一些证据X后的后验概率：

• P(A): 硬币有50%的可能性出现正面。P(A|X): 你看了硬币，观察到正面出现，表示该信息为X，然后直接把正面的概率赋为1.0而反面的概率为　0.0
• P(A): 这大段复杂的代码可能有一个bug。P(A|X): 代码过了所有X的测试；当然仍会有可能出现bug，不过这个出现的概率已经相当小了
• P(A): 病人有任意数量的疾病。P(A|X): 执行了血液测试产生了证据X，排除了部分病因。

在每个例子中，我们都没有十足的把握在看到新的证据X能够完全放弃先验的信念，但是我们对先验信念进行的重新的权衡从而整合进新出现的那些证据。（我们给一些信念新的权值或者信任度）。

通过引入事件的先验可能性，我们已经认可了做出的任何猜测都是潜在地错误的。在观测了数据、证据或者其他信息，我们更新自己的信念，然后我们的猜测会变得略微正确一些。这是预测硬币的另一面，一般来说我们都朝着更加啊正确的方向前行。

实践中的贝叶斯推断

如果频率主义者和贝叶斯推断都是程序函数，其输入是统计问题，那么这两者在返回给用户的方式上是很不同的。频率主义的推断函数可能是返回一个数字，表示一个估计（典型的就是诸如样本均值这样的汇总数据），而贝叶斯函数会返回概率。

例如，在上面我们的调试问题中，调用频率主义函数的方法就是：“我的代码过了所有X个测试；它是不是不含有bug？”会返回一个YES。另一方面，对贝叶斯函数的问题：“我的代码通常会有bug。我的代码过了所有X个测试；我的代码是不是不含有bug？”将会返回非常不同的东西：YES和NO的不同。这个函数会返回：

YES的概率是0.8；NO的概率是0.2

这就和频率主义的观点得到的结果非常不同。请注意贝叶斯函数接受了“我的代码通常会有bug”这个额外的参数。这是先验的知识。通过包含先验的参数，我们告诉贝叶斯函数来加入我们关于处境的信念。技术上看，贝叶斯函数中的参数是可选的，但是我们将看到去除这点，会带来自身的后果。

引入证据

在我们获得了越来越多的证据后，我们的先验信念将被新的证据冲刷。例如，若你先验的信念是某些荒唐的事情（我认为今天太阳会爆炸）然后每天都在证明你是错误的，你可能会希望任何的推断会纠正你的观点，或者至少在某种程度上改变你的信念。贝叶斯推断会纠正这个信念。

用N来表示我们已有证据的例子的数目。当我们搜集了一个无穷多的证据时，不妨设N -> infinity，我们的贝叶斯结果通常和频率主义者相似。因此，对很大的N，统计推断多少是客观的。另一方面，当N比较小的时候，推断更加不稳定：频率主义的估计有更大的方差和置信区间。这个在贝叶斯分析超越。通过引入先验概率，和返回的概率，我们保留了非确定性，这就说明了小的数据集所具备的推断能力。

有人会觉得，对于大的N，我们可以忽略这两种技术之间的区别，并且会倒向计算代价小的频率主义方法。持这样观点的人们在做出上面的决定之前建议看下面Andrew Gelman所讲的这段话：

样本大小从来不会很大。如果N太小以至于难以获得一个充分精确的估计，你需要获得更多的数据（或者做出更多的假设）。但是，一旦N足够的大了，你可以开始把数据进行分解来学习更多的经验（例如，在一个公共观点池中，一旦你有了对于整个国家的好的估计，你可以估计男人和女人，北方和南方，不同的年龄群组等等）。N是从来都不足够的，因为，如果他如果是足够了后，你将需要开一个新的需要更多数据的任务。

频率主义方法就错了么？

不是

频率主义方法仍然在很多的领域是相对有效的。诸如最小方差线性回归，LASSO回归，和期望最大算法还是快速且强大的。贝叶斯方法填补了这些方法留下的空白领域，可以得到更多的灵活的模型。

大数据

让人费解的是，大数据的预测分析问题往往是使用相对简单的算法解决的。因此我们可以这样说，大数据预测的难点不是算法本身，而是数据的存储和计算执行的难度。（读者影带考虑一个Gelman的话并且问“我是真的拥有大数据么？”）

更加困难的分析问题包含了中等量的数据，而最为困难的是相当小的数据。使用一个类似于Gelman上述的观点，如果大数据足够大可以立即解决的问题，那么我们对于那些不是过大的数据集更有兴趣。

我们的贝叶斯框架

我们对信念感兴趣，这个可以解释为按照贝叶斯观点下的概率。我们对于一个事件A有一个先验概率，信念由先前的信息形成，例如我们在执行测试之前相信了bug可能会出现。

另外，我们观测我们的证据。继续我们的bug代码例子：如果我们的代码过了X个测试，我们希望更新增加了这个认识之后的信念。我们称这个新的信念为后验概率。通过下面的公式来更新我们的信念，也就是传说中的贝叶斯定理（以发现者Thomas Bayes命名）：

P(A|X) = P(X|A)P(A) / P(X) \propto P(X|A)P(A)

上面的公式并不是专属于贝叶斯推断的：这是一个数学事实，在贝叶斯推荐之外也经常这样使用。在这里，BI仅仅通过它来联结先验概率P(A)和后验概率P(A|X)。

例子：强制抛硬币例子

每本统计课本肯定会包含抛硬币的例子。我也会使用它。假设，你不确定抛掷一枚硬币出现正面的概率（拆台者说：这不是50%么！！！）。你相信在比例之后暗藏着某种真理，就叫p吧，但是并没有任何关于p究竟是多少的先验的认识。
我们开始抛硬币，记录下来观察值：正面或者反面。这就是我们的观察数据。一个有趣的问题就是我们的推断在观测发生变化后如何改变？具体讲，就是我们的后验概率变成什么样子，在我们数据很少和很多的时候分别是什么样子。
下面我们画出来在我们观测出更多的数据（抛掷硬币的结果）更新后验概率的序列。

from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
figsize(11, 9)import scipy.stats as statsdist = stats.beta
n_trials = [0, 1, 2, 3, 4, 5, 8, 15, 50, 500]
data = stats.bernoulli.rvs(0.5, size=n_trials[-1])
x = np.linspace(0, 1, 100)# For the already prepared, I'm using Binomial's conj. prior.
for k, N in enumerate(n_trials):sx = plt.subplot(len(n_trials) / 2, 2, k + 1)plt.xlabel("$p$, probability of heads") \if k in [0, len(n_trials) - 1] else Noneplt.setp(sx.get_yticklabels(), visible=False)heads = data[:N].sum()y = dist.pdf(x, 1 + heads, 1 + N - heads)plt.plot(x, y, label="observe %d tosses,\n %d heads" % (N, heads))plt.fill_between(x, 0, y, color="#348ABD", alpha=0.4)plt.vlines(0.5, 0, 4, color="k", linestyles="--", lw=1)leg = plt.legend()leg.get_frame().set_alpha(0.4)plt.autoscale(tight=True)plt.suptitle("Bayesian updating of posterior probabilities",y=1.02,fontsize=14)plt.tight_layout()

贝叶斯后验概率的变化

后验概率由曲线表示，我们的非确定性跟曲线的宽度相关的。如上图所示，当我们开始观察数据时，我们的后验概率开始偏移。最终，当我们看到越来越多的数据时，概率紧紧滴集中在p = 0.5的位置附近。

注意看到途中的峰点不是落在0.5。但是继续增加数据后，肯定会落在那个地方。下一个例子会展示贝叶斯推断的数学。

例子：bug还是仅仅甜点，不能预估的未来？

令A表示我们代码没有bug的事件。让X表示代码过了所有测试用例。现在，我们将没有bug的先验概率当作一个变量，即P(A) = p。

我们对P(A|X)感兴趣，就是给定了测试用例X后的没有bug的概率。为使用上面的公式，我们需要计算一些值。

什么是P(X|A)即已知没有bug的情况下过了X个测试的概率？当我傻啊，就是1啊，没有bug的代码肯定会过了所有的测试用例。

P(X)则有一点儿的tricker：事件X可以被分成两种情况：事件X出现，我们的代码是有bug的（~A），或者X出现，代码是没有bug的（A）。公式如下：

P(X) = P(X and A) + P(X and ~A) = P(X|A)P(A) + P(X|~A)P(~A) = P(X|A) p + P(X|~A) (1-p)

我们上面已经计算了P(X|A)。另外，P(X|~A)是主观的；我们的代码很可能在有bug的情况下没有被测试测出来。注意这个值是和已完成的测试的数量及测试的复杂程度等等因素是相关的。让我们保守一些，给出P(X|~A) = 0.5。那么

P(A|X) = 1*p / (1 * p + 0.5 * (1 - p))= 2p / 1 + p

这就是后验概率。这是我们的先验概率p的函数的参数，p \in [0, 1]？

figsize(12.5, 4)
p = np.linspace(0, 1, 50)
plt.plot(p, 2 * p / (1 + p), color="#348ABD", lw=3)
# plt.fill_between(p, 2*p/(1+p), alpha=.5, facecolor=["#A60628"])
plt.scatter(0.2, 2 * (0.2) / 1.2, s=140, c="#348ABD")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel("Prior, $P(A) = p
figsize(12.5, 4)
p = np.linspace(0, 1, 50)
plt.plot(p, 2 * p / (1 + p), color="#348ABD", lw=3)
# plt.fill_between(p, 2*p/(1+p), alpha=.5, facecolor=["#A60628"])
plt.scatter(0.2, 2 * (0.2) / 1.2, s=140, c="#348ABD")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel("Prior, $P(A) = p
figsize(12.5, 4)
p = np.linspace(0, 1, 50)
plt.plot(p, 2 * p / (1 + p), color="#348ABD", lw=3)
# plt.fill_between(p, 2*p/(1+p), alpha=.5, facecolor=["#A60628"])
plt.scatter(0.2, 2 * (0.2) / 1.2, s=140, c="#348ABD")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel("Prior, $P(A) = p
figsize(12.5, 4)
p = np.linspace(0, 1, 50)
plt.plot(p, 2 * p / (1 + p), color="#348ABD", lw=3)
# plt.fill_between(p, 2*p/(1+p), alpha=.5, facecolor=["#A60628"])
plt.scatter(0.2, 2 * (0.2) / 1.2, s=140, c="#348ABD")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel("Prior, $P(A) = p$")
plt.ylabel("Posterior, $P(A|X)$, with $P(A) = p$")
plt.title("Are there bugs in my code?")
quot;)
plt.ylabel("Posterior, $P(A|X)$, with $P(A) = p
figsize(12.5, 4)
p = np.linspace(0, 1, 50)
plt.plot(p, 2 * p / (1 + p), color="#348ABD", lw=3)
# plt.fill_between(p, 2*p/(1+p), alpha=.5, facecolor=["#A60628"])
plt.scatter(0.2, 2 * (0.2) / 1.2, s=140, c="#348ABD")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel("Prior, $P(A) = p
figsize(12.5, 4)
p = np.linspace(0, 1, 50)
plt.plot(p, 2 * p / (1 + p), color="#348ABD", lw=3)
# plt.fill_between(p, 2*p/(1+p), alpha=.5, facecolor=["#A60628"])
plt.scatter(0.2, 2 * (0.2) / 1.2, s=140, c="#348ABD")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel("Prior, $P(A) = p$")
plt.ylabel("Posterior, $P(A|X)$, with $P(A) = p$")
plt.title("Are there bugs in my code?")
quot;)
plt.title("Are there bugs in my code?")quot;)
plt.ylabel("Posterior, $P(A|X)$, with $P(A) = pfigsize(12.5, 4)
p = np.linspace(0, 1, 50)
plt.plot(p, 2 * p / (1 + p), color="#348ABD", lw=3)
# plt.fill_between(p, 2*p/(1+p), alpha=.5, facecolor=["#A60628"])
plt.scatter(0.2, 2 * (0.2) / 1.2, s=140, c="#348ABD")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel("Prior, $P(A) = p
figsize(12.5, 4)
p = np.linspace(0, 1, 50)
plt.plot(p, 2 * p / (1 + p), color="#348ABD", lw=3)
# plt.fill_between(p, 2*p/(1+p), alpha=.5, facecolor=["#A60628"])
plt.scatter(0.2, 2 * (0.2) / 1.2, s=140, c="#348ABD")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel("Prior, $P(A) = p$")
plt.ylabel("Posterior, $P(A|X)$, with $P(A) = p$")
plt.title("Are there bugs in my code?")
quot;)
plt.ylabel("Posterior, $P(A|X)$, with $P(A) = p
figsize(12.5, 4)
p = np.linspace(0, 1, 50)
plt.plot(p, 2 * p / (1 + p), color="#348ABD", lw=3)
# plt.fill_between(p, 2*p/(1+p), alpha=.5, facecolor=["#A60628"])
plt.scatter(0.2, 2 * (0.2) / 1.2, s=140, c="#348ABD")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel("Prior, $P(A) = p$")
plt.ylabel("Posterior, $P(A|X)$, with $P(A) = p$")
plt.title("Are there bugs in my code?")
quot;)
plt.title("Are there bugs in my code?")quot;)
plt.title("Are there bugs in my code?")

figure_2.png

如果我们观测了X测试用例过了当先验概率p很低的时候，我们可以看到最大的增益。如果我自认为是一个牛逼程序员，所有我给自己一个先验的概率为0.20，有20％的概率我的代码是没有bug的。更加实际地，先验概率应该和代码的复杂程度和规模相关，但是就让我这么设置吧。那么我更新后的无bug的信念是0.33。
回忆先验时概率：p是没有任何bug的先验概率，所以1-p是有bug的先验概率。

类似的，我们的后验同样是概率，给定所有测试都通过的情形下没有bug的概率是P(A|X)，因此1 - P(A|X)是给定所有测试用例都通过的情形下出现bug的概率。那么我们的后验概率是什么样的？下面是先验和后验概率的图。

figsize(12.5, 4)
colours = ["#348ABD", "#A60628"]prior = [0.20, 0.80]
posterior = [1. / 3, 2. / 3]
plt.bar([0, .7], prior, alpha=0.70, width=0.25,color=colours[0], label="prior distribution",lw="3", edgecolor=colours[0])plt.bar([0 + 0.25, .7 + 0.25], posterior, alpha=0.7,width=0.25, color=colours[1],label="posterior distribution",lw="3", edgecolor=colours[1])plt.xticks([0.20, .95], ["Bugs Absent", "Bugs Present"])
plt.title("Prior and Posterior probability of bugs present")
plt.ylabel("Probability")
plt.legend(loc="upper left");

figure_3.png

注意到，在我们观测到X发生后，不会出现bug的概率有所增加。通过增加测试的数量，我们确信不会有bug出现。

这是贝叶斯推断和贝叶斯规则的简单例子。不幸的是，进行更加复杂的贝叶斯推断，我们需要更多的数学知识。后面我们会看到这样的数学分析事实上是不必须的。首先我们需要扩展建模工具，下面介绍概率分布。如果你已经相当的熟悉，那就潇洒地跳过吧。

概率分布

让我们快速回忆一下什么是概率分布：令Z为某个随机变量。与Z相关的概率分布函数给Z的不同的值赋予了概率。图形上看起来，概率分布是一个曲线，其某个结果的概率曲线的高度成比例。你可以在本节的第一张图片中看到。
我们可以讲随机变量分成三类：

• 离散：离散随机变量可能仅仅假设值存在于一个指定的列表中。就像人口，电影评分和投票数目这些都是离散的随机变量。
• 连续：连续随机变量可以取任何一个准确的值。例如，温度、速度、事件、颜色等等都可以被建模成连续的值因为你总可以让值越来越精确。
• 混合：混合的随机变量给离散和连续的随机变量赋予概率，他是上面两种的组合。

离散情形

如果Z是离散的，那么他的分布被称为概率质量函数，这个东西度量了Z取值为k的概率，用P(Z=k)表示。注意这个概率质量函数完全地刻画了随机变量Z，如果我们知道了概率质量函数，我们就能了解Z的所有行为。有很多流行的概率质量函数：我们将根据需要引出这些质量函数，但现在我们隆重介绍第一个相当重要的概率质量函数。我们称Z是泊松分布，当：

P(Z=k) = \lambda^k e^{-\lambda} / k!, k = 0, 1, 2, ...

\lambda是这个分布的参数，它控制了分布图形的形状。对泊松分布，\lambda可以是任何正数。通过增加\lambda，我们增加了更多的概率给更大的值，而相反，通过降低\lambda我们增加了更大的概率给了较小的值。也常称\lambda为泊松分布的强度（intensity）

引入我们的第一个工具：PyMC

PyMC是一个用于贝叶斯分析的python库。他是一个快速的处于较好维护状态下的库。但是不太爽的是它本身缺乏某些领域的文档介绍，尤其是衔接初学者和黑客的桥梁的那部分。这本书本身的目标就是来解决这个问题，同样也展示了为何Py MC这样棒！
我们将会使用PyMC对这个问题进行模型。这种类型的程序设计被称为概率程序设计（probabilistic programming），这是一个令人困惑的术语，因为引入了随机产生的代码，可能已经吓住了不少用户了。代码本身不是随机的；其概率意味只是在我们在创建概率模型的时候使用的程序变量作为了模型的一部分。模型组成部分在PyMC框架内是一等的原语（primitives）。
B.Cronin有一段关于概率程序设计的有感染力的描述：

用另外的一种方式看：不同于传统的程序，他们只是运行在一个向前的方向上，概率程序在正向和反向上都有运行。向前运行就是计算了那些它包含的对这个世界的假设产生的后果，但是同样也可以通过数据来限制可能的解释而后向运行。在实践中，很多概率程序设计系统会聪明地交织这些前向和后向操作来高效地达到最佳的解释。
由于术语概率程序设计带来的困惑，我就将其说成编程（programming），因为它本来如是。

PyMC的代码读起来很简单。最为新鲜的就是它的语法了，我会将代码分解开来介绍每个部分的。只需要记住，我们表示模型的成分(\tau, \lambda_1, \lambda_2)为变量：

import pymc as pmalpha = 1.0 / count_data.mean()  # Recall count_data is the# variable that holds our txt 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)

原地址