贝叶斯分析助你成为优秀的调参侠：先验和似然对后验的影响

©PaperWeekly 原创 · 作者｜庞龙刚

学校｜华中师范大学

研究方向｜能核物理、人工智能

上一节 介绍了科学的研究方法，先验、似然、后验以及如何使用 MCMC 方法中的 Metropolis-Hastings 方法从模型参数 满足的非归一化后验分布抽样。

这一节以一个具体的例子介绍先验和似然对后验的影响。

学习内容

• 举例：扔硬币实验

• 扔硬币实验的先验

• 扔硬币实验的似然函数与后验分布

• 先验与测量如何影响后验分布

• 使用 MCMC 算法从后验分布抽样

举例：扔硬币实验

这里以扔硬币实验深入介绍先验、似然和后验。

如果一个人投币 10 次，发现 7 次正面朝上，

问：如果投币无数次，正面朝上的概率 是多少？

频率学派 的会说：概率为 ，其中 k 是 n 次投币中正面朝上的次数。

但只有当投币次数无穷多的情况下可以这样算，只投 10 次，总不能说 吧？

贝叶斯学派的会说，我们需要算个后验概率，

这是一个条件概率，表示 n 次投币观测到 m 次正面向上的情况下， 所满足的分布。

• m out of n 表示实验数据

• 表示模型参数，待定

根据前一节介绍，只需要一个非归一化的贝叶斯公式，

扔硬币实验的先验

上面式子中 表示先验（a prior)。

先验指一个人根据过去经验，脑海中对某事件发生的概率估计（或信仰 Belief）。

对于一枚硬币，估计大部分人会认为正面向上的概率为 50%， 即 或写成：

如果做了 10 次投币实验，7 次正面朝上，部分人可能会完全改变信仰，认为

另一部分人可能认为 在 0.5 和 0.7 之间，满足某种分布。

在先验概率未知时，可以用 Beta 分布函数来表示。忽略归一化常数，

是 Beta 分布的两个参数，改变这两个参数，Beta 分布会有很多不同的形状。

下面可视化一下 Beta 函数的表示能力，（注意代码接上一节）

from scipy.stats import beta    
import numpy as np
import matplotlib.pyplot as plt
# 需要安装 SciencePlots 库：
# pip install SciencePlots
plt.style.use(['science', 'ieee'])

def plot_beta_dist(axis, a=1, b=1):
    '''plot the Beta distribution function
    :a: parameter $\alpha$ in Beta distribution
    :b: parameter $\beta$ in Beta distribution'''
    x = np.linspace(0, 1, 10000)
    y = beta.pdf(x, a, b)
    axis.plot(x, y)
    axis.set_ylabel(r"$P(\theta)$")
    axis.set_xlabel(r"$\theta$")
    axis.set_title(r"$\alpha=%.1f, \beta=%.1f$"%(a, b))
# Let's investigate how does Beta distribution change with alpha and beta

# parameters in 2 rows, 3 columns
params = [[(1, 1), (0.6, 0.6), (10, 10)],
          [(0.3, 1.0), (10, 1), (1, 10)]]

fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(9, 6))
for i in range(2):
    for j in range(3):
        plot_beta_dist(ax[i][j], *params[i][j])

plt.tight_layout()

上面第一行最左边这张图对应 ，此时 Beta 分布是均匀分布，即某人相信硬币正面朝上的概率 都一样。

第一行中间这张图对应 ，此时 Beta 分布是 U 型，表示有另一人相信硬币被做了手脚，要么正面朝上的概率非常大 ，要么反面朝上的概率非常大， 。

第一行右边这张图对应 ，此时 Beta 分布是钟型，表示还有人相信硬币没大问题，正面朝上概率为 50%（ ）的可能最大，也有其它可能， 非 函数。

下面一行对应偏到一边的先验。

无论是哪种先验，经过多次实验，随着不确定性的减少，大脑会逐渐调整对这枚硬币正面朝上概率的估计。

后面会调整先验的参数 ，观察先验对后验分布的影响。

扔硬币实验的似然与后验

做 n 次实验有 m 次正面朝上，与正面朝上概率为 的模型的似然函数为，

即二项分布函数（Binomial distribution）。

二项分布很特殊，一般情况下，我们并不知道模型输出与实验数据的似然函数（likelihood）如何定义。一个比较常用的似然函数是，

其中 是归一化系数，在贝叶斯分析中不重要。

后面的 exp 函数是一个多变量的正态分布函数，方差 里保存模型或数据的不确定性。

选择了 Beta 分布函数做先验，知道了似然函数是二项分布，接下来就可以定义扔硬币实验的后验分布，

后验分布同时受到先验分布参数 的影响以及测量结果 的影响。

先验和测量如何影响后验分布

from ipywidgets import interact
import ipywidgets as widgets
print_to_pdf = False

def plot_posterior(m, n, a, b):
    theta = np.linspace(0.0001, 1, 100, endpoint=False)
    #poster = theta**(a + m - 1) * (1 - theta)**(b + n - m - 1)
    ln_poster = (a +m - 1) * np.log(theta) + (b + n - m - 1) * np.log(1 - theta)
    poster = np.exp(ln_poster)
    prior =  beta.pdf(theta, a, b)
    plt.axvline(m/n, color='c', label='frequency')
    plt.plot(theta, poster/poster.max(), label="posterior")
    plt.plot(theta, prior/prior.max(), 'k--', label="prior")
    plt.legend(loc='best', fontsize="small")
    plt.xlabel(r"$\theta$")

if print_to_pdf:
    plot_posterior(m=35, n=162, a=0.6, b=10)
else:
    interact(plot_posterior, 
         m=(10, 100), n=(20,200), 
         a = (0, 10.0, 0.2), b = (0, 10.0, 0.2))

接下来我们改变先验的参数 ，看看对后验分布的影响。

第一组，实验 110 次，55次正面向上。先验（prior）是虚线，中心在 0.5，但有很宽的分布，后验（posterior）分布是实线，宽度比较窄，表示测量减小了不确定度。频率（frequency）等于 55/110。

第二组，观测数据不变，先验参数为 。先验( Prior ) 偏向右边，而后验分布的中心不在频率（frequency）处，微微偏向先验。如果选 ，先验和后验都会偏向左边。

第三组，观测数据变少，做 20 次实验，10 次正面朝上。先验参数仍为 ，即初始信仰为正面朝上的概率很大。此时可以看到先验对后验的影响增大。后验向先验靠近更多。

总结一下：

• 后验分布的中心一般在先验和频率之间

• 测量次数少时，先验对后验的影响比较大

• 测量次数很多时，先验的影响可忽略不计

使用 MCMC 算法从后验分布抽样

接下来，就要对投硬币任务， 满足的后验分布函数进行抽样，找到对于给定的实验观测，最可能的 的值以及 的分布。

▲ 对投硬币任务，正面朝上概率满足的后验分布抽样。

代码如下：

def draw_posterior(m, n, a, b, nsamples=1000000):
    def poster(theta, m=m, n=n, a=a, b=b):
        ln_poster = (a +m - 1) * np.log(theta) \
                  + (b + n - m - 1) * np.log(1 - theta)
        poster = np.exp(ln_poster)
        return poster


    samples = mcmc(poster, np.array([[0.0001, 0.9999]]),
                   n=nsamples, warmup=0.1, thin=10)

    # compute the distribution of these samples using histogram
    hist, edges = np.histogram(samples.flatten(), bins=100)

    x = 0.5 * (edges[1:] + edges[:-1])
    theta = np.linspace(0.0001, 0.9999, 100)
    y = poster(theta)

    plt.plot(x, hist/hist.max(), '-', label="MCMC")
    plt.plot(theta, y/y.max(), '-', label = r"$P(\theta | \rm m\ out\ of\ n)$")

    plt.legend(loc='best', fontsize='small')
    plt.xlabel(r'$\theta$')
    plt.xlim(x[0], x[-1])

draw_posterior(m=100, n=150, a=10, b=10, nsamples=2000000)

总结

这一节以一个扔硬币任务，更加详细的介绍了何为先验，以及先验、测量对后验分布的影响。使用前一节引入的 MCMC 算法，对非归一化的后验分布抽样。推荐大家在一个具体的物理模型调参任务上走一遍流程，会对贝叶斯推理和 MCMC 有更深刻的认识。

如果物理模型运行速度很慢，可使用 Latin Cube 方法在参数空间取点，使用 Gauss Process Emulator 对参数空间插值，如果实验的可观测量比较多，还可以使用主成分分析将数据降维。LBNL 的 Ke WeiYao 博士曾到我们实验室给过一个贝叶斯分析的讲座，介绍了很多这方面的知识，这一节很多内容借鉴了他开源的 jupyter-notebook。

参考文献

[1] https://github.com/keweiyao/BayesExample

[2] http://lgpang.gitee.io/mlcphysics/htmls/bayes_analysis.slides.html#/

更多阅读

#投 稿 通 道#

 让你的论文被更多人看到 

如何才能让更多的优质内容以更短路径到达读者群体，缩短读者寻找优质内容的成本呢？答案就是：你不认识的人。

总有一些你不认识的人，知道你想知道的东西。PaperWeekly 或许可以成为一座桥梁，促使不同背景、不同方向的学者和学术灵感相互碰撞，迸发出更多的可能性。 

PaperWeekly 鼓励高校实验室或个人，在我们的平台上分享各类优质内容，可以是最新论文解读，也可以是学习心得或技术干货。我们的目的只有一个，让知识真正流动起来。

📝 来稿标准：

• 稿件确系个人原创作品，来稿需注明作者个人信息（姓名+学校/工作单位+学历/职位+研究方向） 

• 如果文章并非首发，请在投稿时提醒并附上所有已发布链接 

• PaperWeekly 默认每篇文章都是首发，均会添加“原创”标志

📬 投稿邮箱：

• 投稿邮箱：hr@paperweekly.site 

• 所有文章配图，请单独在附件中发送 

• 请留下即时联系方式（微信或手机），以便我们在编辑发布时和作者沟通

🔍

现在，在「知乎」也能找到我们了

进入知乎首页搜索「PaperWeekly」

点击「关注」订阅我们的专栏吧

关于PaperWeekly

PaperWeekly 是一个推荐、解读、讨论、报道人工智能前沿论文成果的学术平台。如果你研究或从事 AI 领域，欢迎在公众号后台点击「交流群」，小助手将把你带入 PaperWeekly 的交流群里。