用 Python 进行贝叶斯模型建模(1)

本系列:

第1节:估计模型参数

在这一节,我们将讨论贝叶斯方法是如何思考数据的,我们怎样通过 MCMC 技术估计模型参数。

贝叶斯方法如何思考数据?

当我开始学习如何运用贝叶斯方法的时候,我发现理解贝叶斯方法如何思考数据是很有用的。设想下述的场景:

一个好奇的男孩每天观察从他家门前经过的汽车的数量。他每天努力地记录汽车的总数。一个星期过去后,他的笔记本上记录着下面的数字:12,33,20,29,20,30,18

从贝叶斯方法的角度看,这个数据是由随机过程产生的。但是,既然数据被观测,它便固定了并且不会改变。这个随机过程有些模型参数被固定了。然而,贝叶斯方法用概率分布来表示这些模型参数的不确定性。

由于这个男孩调查的是计数(非负整数),一个通常的做法是用泊松分布对数据(如随机过程)建模。泊松分布只有一个参数 μ,它既是数据的平均数,也是方差。下面是三个不同 μ 值的泊松分布。

 p(x \ | \ \mu) = \frac{e^{-\mu}\mu^{x}} {x!} \mbox{ for }x = 0, 1, 2, \cdots

\lambda = E(x) = Var(\mu)

代码:

在上一节中,我们引入我的 hangout 聊天数据集。特别地,我对我的回复时间(response_time)感兴趣。鉴于 response_time 是计数数据,我们可以用泊松分布对其建模,并估计参数 μ 。我们将用频率论方法和贝叶斯方法两种方法来估计。

频率论方法估计μ

在进入贝叶斯方法之前,让我们先看一下频率论方法估计泊松分布参数。我们将使用优化技术使似然函数最大。

下面的函数poisson_logprob()返回在给定泊松分布模型和参数值的条件下,观测数据总体的可能性。用方法opt.minimize_scalar找到在观测数据基础上参数值 μ 的最可信值(最大似然)。该方法的机理是,这个优化技术会自动迭代可能的mu值直到找到可能性最大的值。

所以,μ 的估计值是18.0413533867。优化技术没有对不确定度进行评估,它只返回一个点,效率很高。

下图描述的是我们优化的函数。对于每个μ值,图线显示给定数据和模型在μ处的似然度。优化器以登山模式工作——从曲线上随机一点开始,不停向上攀登直到达到最高点。

上述优化方法估计泊松分布的参数值为  18 .我们知道,对任意一个泊松分布,参数 μ 代表的是均值和方差。下图描述了这个分布。

上述泊松分布模型和 μ 的估计表明,观测小于10 或大于 30 的可能性很小,绝大多数的概率分布在 10 和 30 之间。但是,这不能反映我们观测到的介于 0 到 60 之间的数据。

贝叶斯方法估计 μ

如果你之前遇到过贝叶斯理论,下面的公式会看起来很熟悉。我从没认同过这个框架直到我读了John K. Kruschke的书“贝叶斯数据分析”,从他优美简洁的贝叶斯图表模型视角认识这个公式。
\overbrace{p(\mu \ |\ Data)}^{\text{posterior}} = \dfrac{\overbrace{p(Data \ | \ \mu)}^{\text{likelihood}} \cdot \overbrace{p(\mu)}^{\text{prior}}}{\underbrace{p(Data)}_{\text{marginal likelihood}}}

代码:

上述模式可以解读如下(从下到上):

  • 对每一个对话(i)观测计数数据(y)
  • 数据由一个随机过程产生,我们认为可以用泊松分布模拟
  • 泊松分布只有一个参数,介于 0 到 60 之间

由于我们没有任何依据对μ在这个范围内进行预估,就用均匀分布对 μ 建模

神奇的方法MCMC

下面的动画很好的演示了马尔科夫链蒙特卡罗方法(MCMC)的过程。MCMC采样器从先验分布中选取参数值,计算从这些参数值决定的某个分布中得到观测数据的可能性。
\overbrace{p(\mu \ |\ Data)}^{posterior} \varpropto \overbrace{p(Data \ | \ \mu)}^{likelihood} \cdot \overbrace{p(\mu)}^{prior}

这个计算可以作为MCMC采样器的导引。从参数的先验分布中选取值,然后计算给定数据条件下这些参数值可能性——导引采样器到更高概率的区域。

与上面讨论的频率论优化技术在概念上有相似之处,MCMC采样器向可能性最高的区域运动。但是,贝叶斯方法不关心寻找绝对最大值,而是遍历和收集概率最高区域附近的样本。所有收集到的样本都可认为是一个可信的参数。

上面的代码通过遍历 μ 的后验分布的高概率区域,收集了 200,000 个 μ 的可信样本。下图(左)显示这些值的分布,平均值与频率论的估值(红线)几乎一样。但是,我们同时也得到了不确定度,μ 值介于 17 到 19 之间都是可信的。我们稍后会看到这个不确定度很有价值。

丢弃早期样本(坏点)

你可能疑惑上面 MCMC 代码中pm.find.MAP()的作用。MAP 代表最大后验估计,它帮助 MCMC 采样器寻找合适的采样起始点。理想情况下,模型从高概率区域开始,但有时不是这样。因此,样本集中的早期样本(坏点)经常被丢弃。

模型的收敛性

样本轨迹

由于上面模型对 μ 的估计不能表明估计的效果很好,可以采用一些检验手段。首先,观察样本输出,样本的轨迹应该轻微上下浮动,就像一条毛虫。如果你看到轨迹上下跳跃或滞留在某点,这就是模型不收敛的标志,通过 MCMC 采样器得到的估计没有意义。

自相关图

也可采另一种测试,自相关测试(见下图),这是评估MCMC采样链中样本前后相关关系的一种方法。相关性弱的样本比相关性强的样本能够给参数估计提供更多的信息。

直观上看,自相关图应该快速衰减到0附近,然后在 0 附近上下波动。如果你的自相关图没有衰减,通常这是弱融合的标志,你需要重新审查你的模型选择(如似然度)和抽样方法(如Metropolis)

打赏支持我翻译更多好文章,谢谢!

打赏译者

打赏支持我翻译更多好文章,谢谢!

任选一种支付方式

1 4 收藏 评论

关于作者:JLee

Keep Learning , Keep Coding 个人主页 · 我的文章 · 13 ·   

相关文章

可能感兴趣的话题



直接登录
跳到底部
返回顶部