因此： f(\theta|y)\propto \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}\frac{n!}{y!(n-y)!}\theta^y(1-\theta)^{n-y}

而针对选定的先验，参数 \alpha 和 \beta 是已知的，针对一组已知的观测数据，抛掷的次数 n 和正面向上的次数 y 也是已知的，因此 \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}​和 \frac{n!}{y!(n-y)!} ​都是与未知参数 \theta 无关的项，他们可以被合并到归一化项当中去，因此我们可以进一步化简：

最终实际的后验 f(\theta|y) 就是在 \theta^{\alpha+y-1}(1-\theta)^{\beta+n-y-1} 的基础上加上一个归一化的项，使之在整个 \theta 的取值域上积分为1。

此时你可能会问，后面我们再该怎么处理？我们从后验分布的概率密度函数的正比表达式： \theta^{\alpha+y-1}(1-\theta)^{\beta+n-y-1} 中惊奇的发现，后验分布同样也是一个 beta 分布，只不过这个 beta 分布的参数变了：

那么我们很容易就能根据先验分布和观测数据得到后验分布了。

5.模拟实验验证¶

首先我们写一段模拟抛硬币的程序，我们手上的这枚硬币，正面向上的实际概率为0.62，我们模拟随机抛掷1000次硬币的试验，并且记录抛掷过程中，抛掷次数为 [5, 10, 20, 100, 500, 1000] 时，正面出现的次数。这个0.62就是我们需要根据观测数据来估计的未知参数，并且我们可以用它来和估计值进行对比。

代码片段：

import random

def bernoulli_trial(p):
    u = random.uniform(0, 1)
    if u <= p:
        return 1
    else:
        return 0

def coin_experiments(n_array, p):
    y = 0
    n_max = max(n_array)
    results = []
    for n in range(1, n_max+1):
        y = y + bernoulli_trial(p)
        if n in n_array:
            results.append((y, n))

    return results

print(coin_experiments([5, 10, 20, 100, 500, 1000], 0.62))

运行结果：

[(2, 5), (4, 10), (11, 20), (60, 100), (306, 500), (614, 1000)]

这个结果中，每一个元组表示的含义为(正面的次数，试验的次数)，记录了完成1000次试验的过程中，抛掷到5次，10次，20次，100次，500次和1000次的时候，相应的正面向上的次数。

例如 (311, 500) 表示当试验进行了500次时，正面出现的次数为311次。

请大家放心的是，这个试验生成的结果完全是按照伯努利试验随机生成的，你可以再重复运行5次该实验，一定会得到完全不同的5组试验结果：

代码片段：

for i in range(5):
    print(coin_experiments([5, 10, 20, 100, 1000, 10000], 0.62))

运行结果：

[(3, 5), (7, 10), (13, 20), (55, 100), (306, 500), (622, 1000)]
[(4, 5), (6, 10), (12, 20), (68, 100), (319, 500), (632, 1000)]
[(2, 5), (6, 10), (12, 20), (57, 100), (319, 500), (645, 1000)]
[(4, 5), (7, 10), (14, 20), (64, 100), (309, 500), (626, 1000)]
[(4, 5), (8, 10), (15, 20), (56, 100), (309, 500), (607, 1000)]

最后，我们利用三组 (\alpha,\beta) 分别为 (0.25,0.25),(1,1),(10,10) 的 beta 分布作为先验分布：

代码片段：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta
import seaborn
seaborn.set()

params = [0.25, 1, 10]
x = np.linspace(0, 1, 100)

plt.plot(x, beta(0.25, 0.25).pdf(x), color='b', label='\\alpha=0.25,\\beta=0.25')
plt.fill_between(x, 0, beta(0.25, 0.25).pdf(x), color='b', alpha=0.25)

plt.plot(x, beta(1, 1).pdf(x), color='g',label='\\alpha=1,\\beta=1')
plt.fill_between(x, 0, beta(1, 1).pdf(x), color='g', alpha=0.25)

plt.plot(x, beta(10, 10).pdf(x), color='r',label='\\alpha=10,\\beta=10')
plt.fill_between(x, 0, beta(10, 10).pdf(x), color='r', alpha=0.25)

plt.gca().axes.set_ylim(0,10)
plt.gca().axes.set_xlabel('\\theta')
plt.gca().axes.set_ylabel('p(\\theta)')
plt.legend()
plt.show()

运行结果：

我们选择：

[(2, 5), (4, 10), (11, 20), (60, 100), (306, 500), (614, 1000)] 作为观测数据，利用贝叶斯推断的方法，来得到后验分布

代码片段：

from scipy.stats import beta
import matplotlib.pyplot as plt
import numpy as np
import seaborn

seaborn.set()

theta_real = 0.62
n_array = [5, 10, 20, 100, 500, 1000]
y_array = [2, 4, 11, 60, 306, 614]

beta_params = [(0.25, 0.25), (1, 1), (10, 10)]
x = np.linspace(0, 1, 100)

fig, ax = plt.subplots(2, 3, sharex=True, sharey=True)

for i in range(2):
    for j in range(3):
        n = n_array[3 * i + j]
        y = y_array[3 * i + j]
        for (a_prior, b_prior), c in zip(beta_params, ('b', 'r', 'g')):
            a_post = a_prior + y
            b_post = b_prior + n - y
            p_theta_given_y = beta.pdf(x, a_post, b_post)
            ax[i, j].plot(x, p_theta_given_y, c)
            ax[i, j].fill_between(x, 0, p_theta_given_y, color=c, alpha=0.25)

        ax[i, j].axvline(theta_real, ymax=0.5, color='k')
        ax[i, j].set_xticks([0, 0.2, 0.4, 0.6, 0.8, 1])
        ax[i, j].set_title('n={},y={}'.format(n, y))

ax[0, 0].set_ylabel('p(\\theta|y)')
ax[1, 0].set_ylabel('p(\\theta|y)')
ax[1, 1].set_xlabel('\\theta')
plt.show()

运行结果：

6.实验结果分析¶

我们来分析一下试验结果，首先我们设置了三种典型的先验分布，均匀分布、 U 型分布和类似正态分布，分别代表了我们对于未知参数 \theta 不同的认识，比如如果我们选择均匀分布，意味着我们认为， \theta 在 [0,1] 范围内的取值是等概率的； U 型分布代表了在靠近0和1的两头， \theta 的取值概率越高，类正态分布表示我们认为越靠近中心0.5的位置，取值概率越高。这三者的认识是截然不同的。

在整个实验的每一个阶段，通过先验分布与观测数据的综合，我们得到了每个阶段的后验结果，我们从图中可以看出，贝叶斯推断得到的是一个后验分布，而不像极大似然估计中得到的是一个具体值。他表示了在给定观测数据的情况下，推断得到的未知参数 \theta 的分布情况。

一般而言我们都会选择后验分布概率密度函数曲线的峰值作为我们最终对于未知参数的估计值。这就是贝叶斯推断中的最大后验概率（ MAP ）准则，即选择在一个给定数据下，具有最大后验概率的值。

从这个结果图中，我们可以看出很多的结论：

首先，随着观测数据的不断增多，后验分布会越来越集中，分布越集中表示对于参数的确定性越高，这很显然，观测数据的增多意味着有更多的数据、更多的信息来更新和支撑我们对于参数的认识。

其次，当观测数据的量足够多的时候，不同的先验分布对应的后验分布都会收敛到一个相同的结果，数据越多，通过最大后验概率准则得到的估计值就会与参数的实际值（黑色竖线）越接近。

7.关于共轭先验的问题¶

那么这里有一个重要的点要强调，当然也是我们当初没有谈到的，先验选择 beta 分布的第三个好处，那就是：

先验分布是 beta 分布，观测数据服从二项分布，得到的后验仍然是 beta 分布，也就是说 beta 分布是二项分布的共轭先验：即将先验 beta 分布与二项分布组合在一起之后，得到的后验分布与先验分布的表达式形式仍然是一样的。除此之外，正态分布也是自身的共轭先验。

这种分布的共轭特性，极大的简化了我们求解后验分布的计算复杂性。但是，共轭分布的情形并不普遍，因此如果不是在共轭先验的条件下去解决贝叶斯推断问题，那么在计算后验分布的过程中将会遇到非常大的困难，后验分布绝大多数情况下就不再是一个标准分布，甚至没有解析解，在这种情况下想了解后验分布的形态将遇到巨大的挑战，基于他去做后续的统计分析将难上加难。

先验分布与后验分布¶

什么是贝叶斯统计？¶

贝叶斯统计的重要特点在于，我们在建模前需要给出模型参数 \theta 的先验分布（prior distribution） p(\theta) 。也就是说，在得到任何数据，或者将任何数据对模型进行拟合之前，我们需要先给定模型参数服从的分布。例如，对于抛掷一枚硬币一次可能出现的结果，我们可以构建一个参数为 \theta 的伯努利分布概率模型 X\sim bernoulli(\theta) 。那么在使用数据估计参数 \theta 之前，我们需要给这个参数设定一个分布（注意，先验分布是关于模型参数的分布，而不是我们建模的对象本身，下文中将要介绍的后验分布，也是关于模型参数的分布）。给出参数的先验分布，是贝叶斯统计的核心部分之一，也是对于长期接触频率学派思想的贝叶斯初学者来说最容易困惑的部分：没有数据的支持，在贝叶斯统计中的先验分布到底如何设定？

实际上，先验分布的存在正是贝叶斯学派相对于频率学派更贴近实际问题的一点。虽然理论上我们在利用现有数据之前，不存在任何数据支持我们对 \theta 做出假设，但在现实中，我们对我们要建模的目标总会存在一些经验判断。由此，我们可以在利用数据来估计 \theta 之前，将我们的实际经验作为有价值的信息加入到模型参数的估计过程当中。

继续上面的例子，反复投掷同一枚硬币，我们通过常识可以知道其正面朝上的概率 \theta 在 0.5 左右浮动。概率的浮动是由这一枚硬币的实际情况导致（密度均匀程度，以及外形是否完美等等）。由于先验分布可以是任意的，我们因此可以依据经验，假设先验分布则

p(\theta)=4xI(0\leq x<0.5)+4(1-x)I(0.5\leq x\leq 1). \

这是我随意举例的一个先验分布。理论上，我们可以使用任何的分布作为先验分布。但一个好的先验分布需要囊括我们对参数的经验判断信息，而上文构造的分布则成功包含了关于 \theta 分布应有的部分特征：

1. 它当 \theta=0.5 时概率密度最大，符合上文所描述的对硬币的认知。
2. 概率密度函数从 \theta=0.5 处向两端递减，这同样具有现实意义。我们几乎不可能找到正面朝上概率为 0.99 的硬币，甚至无法想象。

在得到了模型参数的先验分布之后，我们便可以通过样本数据 x_1,...,x_n 来进一步更新模型参数的分布。这便是贝叶斯学派与频率学派的另一个重要的不同之处：在贝叶斯统计中，数据被用于更新参数的分布，而非作为参数极大似然估计的组成部分。事实上，我们已经能够通过研究模型参数的先验分布 p(\theta) ，了解 \theta 的性质。而利用数据更新模型参数的分布，能够使我们更加充分地利用已知信息进行统计推断。更新模型参数分布的过程，本质上是在给定 x_1,...,x_{10} 下得到的关于 \theta 的新的分布，即 p(\theta|x_1,...,x_{10})=p(\theta|X) ，其中 p(\theta|X) 被称作 \theta 的后验分布（posterior distribution）。根据贝叶斯定理，我们有

下面是一个特殊的例子，能够帮助读者深化对于贝叶斯思想的印象：

Example¶

当我非常肯定，所有的硬币投掷正面朝上概率为0.5时，我可以设先验分布 p(\theta=0.5)=1 , 除此之外 p(\theta)=0 。由于 \theta 的后验分布表示如下：

小结¶

事实上，本文所举的考察二项分布的参数 \theta 是较为简单的应用情形，贝叶斯学派的这一套思想可以进一步用来研究泊松分布、指数分布、正态分布乃至线性回归模型；除此之外，本文所构造的先验分布也是相对随意的。在实际问题中，我们往往会构造共轭先验分布（conjugate prior distribution），保证我们在不断使用数据更新参数后验分布时，该分布与先验分布属于同一已知分布族当中。我们将在下一篇文章中介绍共轭先验分布，敬请期待。

附录¶

1. 贝叶斯公式

我们回忆贝叶斯公式的形式：