你是一名程序员，今天你实现了一个非常复杂的算法。你不知道这段代码有没有 bug。根据你以往的经验，写出 bug 是在所难免的，毕竟你不是 Jeff Dean。你先拿一些简单的测试用例 run 了一下，通过了。你又试了一些复杂的测试用例，也通过了。接下来，你连续试了一大堆用例，都通过了。于是你松了一口气，开始觉得这段程序可能没有 bug 了。
如果你这样思考问题，那么恭喜你，你是一个贝叶斯主义者。所谓贝叶斯推断，说白了，就是在观察到新的证据之后，更新你的信念。贝叶斯主义者往往不会对某件事确信无疑，但他可以对某件事很有信心。拿上面的例子来说，我们无法保证代码 100% 没有 bug（除非穷尽所有可能的情况）。但当代码通过了大量的测试用例之后，我们对它更有信心了。
我们经常用到“概率”这个词，但对什么是“概率”，其实是有不同看法的。
- 在贝叶斯学派看来，一个事件的概率是人们根据已有的知识和经验对该事件发生的可能性所给出的个人信念[1]，这种信念用区间
中的一个数来表示，可能性大的对应较大的数。
- 与之相对的，在频率学派看来，一个事件在独立重复实验中发生的频率会趋于一个极限，这个极限就是该事件的概率。
乍一看，频率学派的观点好像更严谨。举例来说，你想知道抛一枚硬币国徽朝上的概率，你就不停地抛这枚硬币，统计国徽朝上的次数的占比，随着你抛的次数的增多，你会发现这个占比渐渐趋于稳定，这个极限值就是国徽朝上的概率。问题是，有些情形下，你是无法做独立重复实验的。比如老板问你目前负责的某个项目有多大的概率能拿到业务成果。这个项目只能做一次，要么成功，要么失败，怎么计算成功的频率呢？对此，频率学派可能会告诉你，假设有很多个平行世界，这个项目在一些平行世界里成功了，而在另一些平行世界里失败了。这个项目成功的概率就是在平行世界里成功的频率的极限。但作为贝叶斯主义者，我们不需要诉诸虚无缥缈的平行世界。基于以往的工作经验和对这个项目现状的分析，我们依然能够给出我们的信念，“老板，我有 7 成把握，这个项目能成。”
我们把对事件 发生的可能性的信念记为
，称作事件
的先验概率（prior probability）。当观察到证据
之后，我们更新了对事件
的信念。更新后的信念记为
，称作事件
的后验概率（posterior probability）[2]。那么，在观察到证据
之后，我们怎么才能更新先验概率来得到后验概率呢？很简单，我们使用贝叶斯定理[3]：
是当事件
发生时，我们观察到证据
的概率。它是关于证据
的函数，称作证据
的似然函数（likelihood）。
例 1：2020 年 12 月 23 日，你看完勇士 VS 篮网的比赛，库里三分球 10 投 2 中。你想预测一下库里在本赛季的三分球命中率。因为这是库里本赛季的第一场比赛，直接用
显然是不合理的。那么该怎么做呢？
查阅资料，我发现库里在最近的两个赛季三分球一共是 859 投 366 中。用 表示库里本赛季三分球的命中率。根据历史数据，我认为
服从参数为
和
的 Beta 分布。即我对库里本赛季三分球命中率
的信念为
Beta 分布的期望为
所以，在开始的时候，我觉得库里本赛季的三分球命中率会在 附近波动。
已知库里的三分命中率为 ，那么他投
次三分球命中的次数
应该服从参数为
和
的二项分布，即似然为
在看到本场比赛库里 10 投 2 中的结果后，我就可以更新我的信念了。我们有
其中
又有
因此
于是，我的信念更新为
也就是说，在观察到库里首场比赛三分球 10 投 2 中之后，我更新了我的信念，认为他本赛季的三分球命中率将会在 附近波动。
看到这些计算过程，有些读者可能已经崩溃了。然而这已经算是最简单的情形之一了。但是别着急，我们并不一定要把后验分布计算出来！我们所关心的无非是后验分布的期望、方差、偏度、峰度之类的统计量，如果我们能够设法从后验分布中采样，就可以利用这些样本来估计这些统计量。那么，有什么方法可以从一个复杂的概率分布中采样呢？非常幸运的是，马尔可夫链蒙特卡洛（Markov chain Monte Carlo, MCMC）方法可以做到这一点。关于这个方法，后续我会专门介绍。这里我们只讲一下实际操作。
在 Python 下，我们可以使用 PyMC3 轻松实现 MCMC 采样。拿上面的例子来说，
import pymc3 as pm
import seaborn as sns
with pm.Model() as model:
# theta 的先验分布为 Beta 分布
theta = pm.Beta(
name='theta',
alpha=366,
beta=493
)
# 似然函数为二项分布
likelihood = pm.Binomial(
name='likelihood',
p=theta,
n=10,
observed=[2]
)
# 利用 MCMC 方法从后验分布中采样
posterior = pm.sample()
# 利用后验分布的样本估计后验分布的期望
print(posterior.theta.mean())
# 0.4234174176263356
# 样本的直方图
sns.distplot(posterior.theta)
例 2： 下图展示了某商品在一段时间内每天的销量。请问在这段时间内销量的平均水平是否发生过变化？
用 表示第
天的销量，这种计数类型的随机变量适合用 Poisson 分布来刻画，即似然为
假设在第 天，销量的水平发生了变化，则
我们只需要给出 、
和
的先验分布，然后利用 MCMC 方法从他们的后验分布中采样即可。
关于 和
，我们只知道他们一定是正数，且和这段时间的平均销量应该相差不多。因此，我们假设他们的先验分布为指数分布
其中 为这段时间的平均销量。
关于 ，我们所知甚少，只能假设它的先验分布是离散均匀分布
同样地，我们使用 PyMC3 从后验分布中采样
with pm.Model() as model:
alpha = 1.0 / sales.mean()
lambda_1 = pm.Exponential('lambda_1', alpha)
lambda_2 = pm.Exponential('lambda_2', alpha)
tau = pm.DiscreteUniform('tau', lower=0, upper=len(sales)-1)
t = np.arange(len(sales))
lambda_ = pm.math.switch(t < tau, lambda_1, lambda_2)
likelihood = pm.Poisson('likelihood', lambda_, observed=sales)
posterior = pm.sample(20000, tune=8000, step=pm.Metropolis())
结果如下：
根据采样的结果，
这几乎就是正确答案了，因为所谓的“销量”实际上是我用下面的代码随机生成的
from scipy.stats import poisson
import numpy as np
sales = []
for _ in range(52):
sales.append(poisson.rvs(mu=18))
for _ in range(38):
sales.append(poisson.rvs(mu=23))
sales = np.array(sales)
参考资料
- Probabilistic Programming & Bayesian Methods for Hackers
- Beta distribution - Wikipedia
- Binomial distribution - Wikipedia
- 斯蒂芬-库里NBA球员虎扑篮球
