R语言Gibbs抽样的贝叶斯简单线性回归仿真分析

贝叶斯分析的许多介绍都使用了相对简单的教学实例(例如,根据伯努利数据给出成功概率的推理)。

由Kaizong Ye,Liao Bao撰写

这篇文章将概述这些原理如何扩展到简单的线性回归。我将导出感兴趣参数的后验条件分布,给出用于实现Gibbs采样器的R代码,并提出所谓的网格点方法

贝叶斯模型

假设我们观察数据

对于我们的模型是

×

Gibbs采样是一种特殊的马尔可夫链算法,常被用于解决包括矩阵分解、张量分解等在内的一系列问题,也被称为交替条件采样(alternating conditional sampling),其中,“交替”一词是指Gibbs采样是一种迭代算法,并且相应的变量会在迭代的过程中交替使用,除此之外,加上“条件”一词是因为Gibbs采样的核心是贝叶斯理论,围绕先验知识和观测数据,以观测值作为条件从而推断出后验分布。

简单来说,已知观测值为 [公式] ,给定一个带有参数的向量 [公式] ,若[公式] 表示第 [公式] 个参数 [公式] 在第 [公式] 次迭代的采样值,则该采样值随机地取自概率分布 [公式] .

[示例] 多元正态分布

问题描述:假设观测值 [公式] 服从均值(未知)为 [公式] ,协方差(已知)为 [公式] ,且变量 [公式] 服从均匀分布,如何利用Gibbs采样估计 [公式] ?

在贝叶斯分析中,后验分布 [公式] 通常被定义为

[公式]

其中, [公式] 是先验分布,[公式] 是关于观测值的似然函数,符号 [公式] 表示“正比于”。

因此,这里的后验分布为

[公式]

尽管 [公式] 可以直接通过这里的联合后验分布得到,但需要注意的是,Gibbs采样要求参数 [公式] 和 [公式] 在迭代过程中交替采样,即

[公式]

[公式]

[示例] Gibbs采样过程

假设协方差矩阵中的 [公式] ,观测值为 [公式] ,并给定四个独立变量 [公式] ,相应的初始值为 [公式] , [公式] , [公式] , [公式] ,如图1.

图1 参数在前10代的采样值

迭代500次后,参数 [公式] 趋于收敛,如图2.

图2 参数迭代500次的所有采样值


有兴趣的是作出推论

如果我们在方差项之前放置正态前向系数和反伽马,那么这个数据的完整贝叶斯模型可以写成:

假设超参数

是已知的,后面可以写成一个常数的比例,


热门课程

R语言数据分析挖掘必知必会

面对扑面而来的数据浪潮,包含Google、Facebook等国际企业,都已采用R语言进行数据分析

探索课程

括号中的术语是数据或可能性的联合分布。其他条款包括参数的联合先验分布(因为我们隐含地假设独立前,联合先验因素)。

伴随的R代码的第0部分为该指定的“真实”参数从该模型生成数据。我们稍后将用这个数据估计一个贝叶斯回归模型来检查我们是否可以恢复这些真实的参数。


吉布斯采样器

为了从这个后验分布中得出,我们可以使用Gibbs抽样算法。吉布斯采样是一种迭代算法,从每个感兴趣的参数的后验分布产生样本。它通过按照以下方式从每个参数的条件后面依次绘制:

可以看出,剩下的1,000个抽签是从后验分布中抽取的。这些样本不是独立的。绘制顺序是随机游走在后空间,空间中的每一步取决于前一个位置。通常还会使用间隔期(这里不做)。这个想法是,每一个平局可能依赖于以前的平局,但不能作为依赖于10日以前的平局。

条件后验分布

要使用Gibbs,我们需要确定每个参数的条件后验。

它有助于从完全非标准化的后验开始:

为了找到参数的条件后验,我们简单地删除不包含该参数的关节后验的所有项。例如,常数项

条件后验:

同样的,

条件后验可以被认为是另一个逆伽马分布,有一些代数操作。

条件后验不那么容易识别。但是如果我们愿意使用网格方法,我们并不需要经过任何代数。

考虑网格方法。网格方法是非常暴力的方式(在我看来)从其条件后验分布进行抽样。这个条件分布只是一个函数。所以我们可以评估一定的密度值。在R表示法中,这可以是grid = seq(-10,10,by = .001)。这个序列是点的“网格”。

那么在每个网格点评估的条件后验分布告诉我们这个抽取的相对可能性。

然后,我们可以使用R中的sample()函数从这些网格点中抽取,抽样概率与网格点处的密度评估成比例。

这在R代码的第一部分的函数rb0cond()和rb1cond()中实现。

使用网格方法时遇到数值问题是很常见的。由于我们正在评估网格中未标准化的后验,因此结果可能会变得相当大或很小。这可能会在R中产生Inf和-Inf值。

例如,在函数rb0cond()和rb1cond()中,我实际上评估了派生的条件后验分布的对数。然后,我通过从所有评估的最大值减去每个评估之前归一化,然后从对数刻度取回。

我们不需要使用网格方法来从条件的后面绘制。

因为它来自已知的分布

请注意,这种网格方法有一些缺点。

首先,这在计算上是复杂的。通过代数,希望得到一个已知的后验分布,从而在计算上更有效率。

其次,网格方法需要指定网格点的区域。如果条件后验在我们指定的[-10,10]的网格间隔之外具有显着的密度?在这种情况下,我们不会从条件后验得到准确的样本。记住这一点非常重要,并且需要广泛的网格间隔进行实验。所以,我们需要聪明地处理数字问题,例如在R中接近Inf和-Inf值的数字。

仿真结果

现在我们可以从每个参数的条件后验进行采样,我们可以实现Gibbs采样器。这是在附带的R代码的第2部分中完成的。它编码上面在R中概述的相同的算法。

结果很好。下图显示了1000个吉布斯(Gibbs)样品的序列。红线表示我们模拟数据的真实参数值。第四幅图显示了截距和斜率项的后面联合,红线表示轮廓。

总结一下,我们首先推导了一个表达式,用于参数的联合分布。然后我们概述了从后面抽取样本的Gibbs算法。在这个过程中,我们认识到Gibbs方法依赖于每个参数的条件后验分布的顺序绘制。 这是一个容易识别的已知的分布。对于斜率和截距项,我们决定用网格方法来规避代数。


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498