抛硬币问题

抛硬币问题是统计学中的一个经典问题，描述如下：随机抛一枚硬币，重复一定次数，记录其正面朝上和反面朝上的次数，根据这些数据，回答诸如“这枚硬币是否”，以及“这枚硬币有多不公平”等问题。这些问题看起来似乎有点儿无聊，不过可别低估了它。抛硬币问题是一个学习贝叶斯统计非常好的例子，一方面是因为几乎人人都熟悉抛硬币这一过程，另一方面是因为这个模型很简单，我们可以很容易计算并解决这个问题。此外，很多真实问题都包含两个互斥的结果，例如0或者1、正或者负、奇数或者偶数、垃圾邮件或者正常邮件、安全或者不安全、健康或者不健康等。因此，即便我们讨论的是硬币，这个模型也同样适用于前面那些问题。

为了估计硬币的偏差，或者更广泛地说，想要用贝叶斯定理解决问题，我们需要数据和一个概率模型。对于抛硬币这个问题，假设我们已经验了一定次数并且记录了正面朝上的次数，也就是说数据部分已经准备好了，剩下的就是模型部分了。考虑到这是第一个模型，我们会列出所有必要的数学公式（别怕，我保证这部分不难），并且一步一步推导。第2章中，我们会回顾这个问题，并借用PyMC3和计算机来完成数学计算部分。

1．通用模型

首先，要抽象出偏差的概念。如果一枚硬币总是正面朝上，那么我们说它的偏差就是1；如果总是反面朝上，那么它的偏差就是0；如果正面朝上和反面朝上的次数各占一半，那么它的偏差就是0.5。这里用参数来表示偏差，用y表示N次抛硬币实验中正面朝上的次数。根据贝叶斯定理，即式（1.4），首先需要指定先验和似然。让我们先从似然开始。

2．选择似然

假设多次抛硬币的结果相互之间都没有影响，也就是说每次抛硬币都是相互独立的，同时还假设结果只有两种可能，正面朝上或者反面朝上。更进一步，我们假设所有硬币都服从相同的分布，因此，抛硬币过程中的随机变量是一个典型的独立同分布变量。但愿你能认同我们对这个问题做出的合理假设。基于这些假设，一个不错的似然候选是二项分布：

　　（1.9）

这是一个离散分布，表示的是N次抛硬币实验中y次正面朝上的概率（或者更一般的描述是，N次实验中，y次成功的概率）。

n_params = [1, 2, 4]  # 实验次数
p_params = [0.25, 0.5, 0.75]  # 成功概率
　
x = np.arange(0, max(n_params)+1)
f,ax = plt.subplots(len(n_params), len(p_params), sharex=True, 
                    sharey=True, 
                    figsize=(8, 7), constrained_layout=True)
　
for i in range(len(n_params)): 
    for j in range(len(p_params)):
        n = n_params[i] 
        p = p_params[j]
　
        y = stats.binom(n=n, p=p).pmf(x)
　
        ax[i,j].vlines(x, 0, y, colors='C0', lw=5)
        ax[i,j].set_ylim(0, 1)
        ax[i,j].plot(0, 0, label="N = {:3.2f}\nθ = 
                           {:3.2f}".format(n,p), alpha=0)
        ax[i,j].legend()
        ax[2,1].set_xlabel('y')
        ax[1,0].set_ylabel('p(y | θ, N)')
        ax[0,0].set_xticks(x)
        plt.savefig('B11197_01_03.jpg', dpi=300)

图1.3

图1.3展示了9个二项分布，每个子图中的标签显示了对应的参数。注意，在这个图中，我并没有省略y轴的值，这么做的目的是，读者能够自己动手确认每个子图的y值之和为1。也就是说，对离散分布而言，对应的y值就是概率值。

似然使用二项分布是一个合理选择，直观上讲，可以看作抛一次硬币时正面朝上的可能性（对于N=1，这很直观，不过对于任意的N值都适用），你可以将与图中y =1的柱状图对比。

假如知道了，那么就可以从二项分布得出硬币正面朝上的分布。关键是我们不知道！不过别灰心，在贝叶斯统计中，每当我们不知道某个参数的时候，就对它赋予一个先验，接下来继续选择先验。

3．选择先验

这里我们选用贝叶斯统计中最常见的一个分布——贝塔分布，作为先验，其数学形式如下：

　　（1.10）

仔细观察式（1.10）可以看出，除了部分之外，贝塔分布和二项分布看起来很像。是希腊字母中大写的伽马，用来表示伽马函数。现在我们只需要知道，用分数表示的第一项是一个正则化常量，用来保证该分布的积分为1，此外和两个参数用来控制具体的分布形态。贝塔分布是我们到目前为止见到的第三个分布，利用下面的代码，可以深入了解其形态。

params = [0.5, 1, 2, 3]
x = np.linspace(0, 1, 100)
f, ax = plt.subplots(len(params), len(params), sharex=True, 
                     sharey=True, 
                     figsize=(8, 7), constrained_layout=True)
for i in range(4): 
    for j in range(4):
        a = params[i] 
        b = params[j]
        y = stats.beta(a, b).pdf(x)
        ax[i,j].plot(x, y)
        ax[i,j].plot(0, 0, label="α = {:2.1f}\nβ = {:2.1f}".format(a, 
                     b), alpha=0)
        ax[i,j].legend()
ax[1,0].set_yticks([])
ax[1,0].set_xticks([0, 0.5, 1])
f.text(0.5, 0.05, 'θ', ha='center')
f.text(0.07, 0.5, 'p(θ)', va='center', rotation=0)

图1.4

我很喜欢贝塔分布以及图1.4中的各种形状，但是我们为什么要在模型中用它呢？在抛硬币这个问题以及一些其他问题中使用贝塔分布的原因有很多，其中第一个原因是，贝塔分布的范围限制在0到1，这跟参数一样。通常，在需要对二值变量建模时，就会选用贝塔分布。第二个原因是其通用性，从图1.4可以看出，该分布可以有多种形状，包括均匀分布、类高斯分布、标准正态分布等。第三个原因是，贝塔分布是二项分布（之前我们用它描述似然）的共轭先验。似然的共轭先验是指，将这个先验分布与似然组合在一起之后，得到的后验分布与先验分布的表达式形式仍然是一样的。简单点儿说，就是每次用贝塔分布作为先验、二项分布作为似然时，会得到一个贝塔分布的后验。除贝塔分布之外还有很多其他共轭先验，例如高斯分布，它的共轭先验就是它自己。很多年来，贝叶斯分析都限定在共轭先验范围内，这主要是因为共轭能让后验在数学上变得更容易处理，要知道贝叶斯统计中一个常见问题的后验都很难从分析的角度去解决。在建立合适的计算方法来解决任意后验之前，这只是个折中的办法。从第2章开始，我们将学习使用现代的计算方法来解决贝叶斯问题而不必考虑是否使用共轭先验。

4．推导后验

首先回忆一下，如贝叶斯定理，即式（1.4）所述，后验正比于似然乘先验。对于我们的问题，需要将二项分布乘贝塔分布：

　　（1.11）

现在，将上式简化。针对我们的实际问题，可以先把与不相关的项去掉而不影响结果，于是得到下式：

　　（1.12）

重新整理之后得到：

　　（1.13）

细心的读者可以看出式（1.13）和贝塔分布的形式很像（除了归一化部分），其对应的参数分别为和。也就是说，在抛硬币这个问题中，后验分布是如下贝塔分布：

　　（1.14）

5．计算后验并画图

现在已经有了后验的表达式，可以用Python对其计算并画出结果。下面的代码中，其实只有一行是用来计算后验结果的，其余的代码都是用来画图的。

plt.figure(figsize=(10, 8))
　
n_trials = [0, 1, 2, 3, 4, 8, 16, 32, 50, 150]
data = [0, 1, 1, 1, 1, 4, 6, 9, 13, 48]
theta_real = 0.35
　
beta_params = [(1, 1), (20, 20), (1, 4)] 
dist = stats.beta
x = np.linspace(0, 1, 200)
　
for idx, N in enumerate(n_trials): 
    if idx == 0:
        plt.subplot(4, 3, 2) 
        plt.xlabel('θ')
    else:
        plt.subplot(4, 3, idx+3) 
        plt.xticks([])
    y = data[idx]
    for (a_prior, b_prior) in beta_params:
        p_theta_given_y = dist.pdf(x, a_prior + y, b_prior + N - y) 
        plt.fill_between(x, 0, p_theta_given_y, alpha=0.7)
　
    plt.axvline(theta_real, ymax=0.3, color='k')
    plt.plot(0, 0, label=f'{N:4d} trials\n{y:4d} heads', alpha=0) 
    plt.xlim(0, 1)
    plt.ylim(0, 12) 
    plt.legend() 
    plt.yticks([])
plt.tight_layout()

在图1.5的第一个子图中，还没做任何试验，对应的3条曲线分别表示3种先验。

图1.5

● 蓝色表示均匀分布的先验，其含义是，偏差的所有可能取值都是等概率的。

● 橙色表示类高斯分布的先验，集中在值0.5附近，该先验反映了这样一种信息：通常硬币正面朝上和反面朝上的概率大致是差不多的。我们也可以说，该先验与“大多数硬币是公平的”这一信念是相符合的。虽然信念这个词在贝叶斯相关的讨论经常用到，但我们最好是讨论受数据影响的模型和参数。

● 绿色表示左偏的先验，更倾向于反面朝上。

其余子图描绘了后续实验的后验分布，每个子图都标注了实验（抛硬币）的次数和正面朝上的次数。此外每个子图中在横轴0.35附近还有一条黑色的竖线，表示的是的真实值。显然，在实际情况下，我们并不知道这个值，在这里标识出来只是方便理解。从图1.5中可以学到很多贝叶斯分析方面的知识，我们花点时间理解。

● 贝叶斯分析的结果是后验分布而不是某个值，该分布描述了不同数值的可能性，是根据给定数据和模型得到的。

● 后验最可能的值是由后验分布的众数决定的（也就是后验分布的峰值）。

● 后验分布的展开度与我们对参数的不确定性成正比例；分布越离散，不确定性越大。

● 直观上讲，观测到的数据越多，那么我们对某一结果就更有信心，尽管1/2=4/8=0.5，但是相比在4次试验中观测到2次正面朝上，在8次试验中观测到4次正面朝上，会让我们对偏差等于0.5这一描述更有信心。该直觉也同时反映在了后验分布上。你可以从图1.5中确认这一点，尤其注意第6个子图中的蓝色后验，尽管众数是一样的，但是第3个子图中曲线的离散程度（不确定性）相比第6个子图要更大。

● 在给定足够多的数据时，两个或多个不同先验的贝叶斯模型会趋近于收敛到相同的结果。在极限情况下，如果有无限多的数据，不论我们使用的是怎样的先验，最终都会得到相同的后验。注意这里说的无限多是指一个极限状态而非某个具体的数量，也就是说，从实际的角度来讲，通过有限且数量较少的数据点就可以得到近似的后验。

● 不同后验收敛到相同分布的速度取决于数据和模型。从图1.5中可以看出，蓝色的先验（均匀分布）和绿色的先验（左偏分布）在经过8次实验之后就很难看出区别了，而橙色的后验（类似高斯分布）即使进行了150次实验之后还能很容易地看出与另外两个后验的区别。

● 有一点从图1.5中不太容易看出来，如果我们一步一步地更新后验，最后得到的结果其实跟一次性计算得到的结果是一样的。换句话说，我们可以对后验进行150次计算，每次增加一个新的观测数据并将得到的后验作为下一次计算的先验，也可以在得到150次抛硬币的结果之后一次性计算出后验，而这两种计算方式得到的结果是完全一样的。这个特点非常有意义，当我们得到新的数据时，就可以很自然地更新估计值。这在很多数据分析问题中都很常见。

6．先验的影响以及如何选择合适的先验

从前面的例子可以明显看出，先验对推理会有影响，这很正常，先验就是这样。一些贝叶斯分析的新手（以及一些诋毁该方法的人）会对如何选择先验感到茫然，因为他们不希望先验起到决定性作用，而是更希望数据本身替自己说话！有这样的想法很正常，不过我们得牢记，数据并不会真地“说话”，最多也只是发出“嗡嗡”声。数据只有在模型中才会有意义，包括数学模型和心理模型。面对同一主题下的同一份数据，不同人会有不同的看法，这类例子在科学史上有很多，即使你基于的是正式的模型，这种情况也还是会发生。

有些人青睐于使用无信息先验（也称作扁平先验[2]、模糊先验或扩散先验），这类先验对分析过程的影响最小。尽管这么做是可行的，不过通常我们能做得比这更好。本书将遵循Gelman、McElreath和Kruschke[3]这3位的建议，更倾向于使用弱信息先验（weakly-information prior）。在很多问题中，我们对参数可以取的值一般都会有些了解。比如，我们可能知道参数只能是正数，或者是知道参数近似的取值范围，又或者是希望该值接近0或大于/小于某值。这种情况下，我们可以给模型加入一些微弱的先验信息而不必担心它会掩盖数据本身的信息。由于这类先验会让后验近似位于某一合理的边界内，因此也被称作正则化先验。当然，如果已经有一些高质量的信息用于定义先验，那么使用带有较多信息量的强先验也是可行的。视具体的问题不同，有可能很容易或者很难找到这类先验，例如在我工作的领域（结构生物信息学），人们会尽可能地利用先验信息，通过贝叶斯或者非贝叶斯的方式来了解和预测蛋白质的结构。这样做是合理的，原因是我们在数十年间已经从上千次精心设计的实验中收集了数据，因而有大量可信的先验信息可供使用，如果不用的话也太荒唐了！所以，如果你有可信的先验信息，完全没有理由不去使用。试想一下，如果一个汽车工程师每次设计新车的时候，他都要重新发明内燃机、轮子乃至整个汽车！这显然不是正确的方式。

---

[2]　扁平先验（flat prior）用均匀分布来表示某个参数的先验分布，我们对参数的取值没有稳定的趋向，会认为参数取任何值的概率都相同，因此扁平先验也是无信息先验（non-information prior），即不提供任何信息的先验。

[3]　这3位分别是Bayesian Data Analysis、Statistical Rethinking: A Bayesian Course with Examples in R and Stan和Doing Bayesian Data Analysis的主要作者。——译者注

现在我们知道了先验有很多种，不过这并不能缓解我们选择先验时的焦虑。或许，最好是没有先验，这样事情就简单了。不过，不论是否基于贝叶斯，模型都在某种程度上拥有先验，即使这里的先验并没有明确表示出来。事实上，很多频率统计学方面的结果在某些情况下可以看作贝叶斯模型的特例，比如扁平先验。在频率学派中，一种常见的参数估计方法是最大似然法，该方法避免了设置先验，其原理就是找到合适的让似然最大化，该值通常用一个带尖的符号来表示，比如 或者。 是一个点估计（一个数值），而不是一个分布，对于前面抛硬币的问题，可以表示为：

　　（1.15）

如果你重新看图1.5，你可以亲自确认，蓝色后验（对应均匀/扁平先验）的众数是否与每个子图计算的 值相同。因此，至少对这个例子而言，尽管最大似然的方法并没有显式地引入任何先验，但仍然可以看作贝叶斯模型的一种（带均匀先验的）特例。

我们无法真正避免先验，不过如果你在分析中引入先验，得到的是合理值的分布而不只是最可能的一个值。明确引入先验的另一个好处是，我们会得到更透明的模型，这意味着这些模型更容易评判、调试（广义上）以及优化。构建模型是一个迭代的过程，有时候可能只需要几分钟，有时候则可能需要数年；有时候整个过程可能只涉及你本人，有时候则可能涉及你不认识的人。另外，模型复现很重要，而模型中透明的假设能有助于复现。在特定分析任务中，如果我们对某个先验或者似然不确定，可以使用多个先验或者似然进行尝试。模型构建过程中的一个环节就是质疑假设，而先验就是质疑的对象之一。不同的假设会得到不同的模型，根据数据和与问题相关的领域知识，我们可以对这些模型进行比较，本书第5章会深入讨论该部分内容。由于先验是贝叶斯统计中的一个核心内容，在接下来遇到新的问题时我们还会反复讨论它，因此如果你对前面讨论的内容感到有些疑惑，别太担心，先冷静冷静，要知道人们在这个问题上已经困惑了数十年并且相关的讨论一直在继续。