首页
学习
活动
专区
工具
TVP
发布
精选内容/技术社群/优惠产品,尽在小程序
立即前往

使用R语言进行Metroplis-in-Gibbs采样和MCMC运行分析

全文链接:http://tecdat.cn/?p=12200

这篇文章展示了我们如何使用Metropolis-Hastings(MH)从每次Gibbs迭代中的非共轭条件后验对象中进行采样–比网格方法更好的替代方法。

相关视频

我将说明该算法,给出一些R代码结果,然后分析R代码以识别MH算法中的瓶颈。

模型

此示例的模拟数据是包含患者的横截面数据集。有一个二元因变量Y,一个二元处理变量A,一个因子变量age。年龄是具有3个等级的分类变量。我用贝叶斯逻辑回归建模:

对于Metroplis-in-Gibbs吉布斯采样来说,这是一个相当不错的示例:

我们有一个二进制结果,为此我们采用了非线性链接函数。

我们有一个需要调整的因素。

我们正在估计我们关心的更多参数,但肯定会给采样器带来压力。

非规范条件后验

让我们看一下该模型的(非标准化)条件后验。

此条件分布不是已知分布,因此我们不能简单地使用Gibbs从中进行采样。相反,在每个gibbs迭代中,我们需要另一个采样步骤来从该条件后验中提取。第二个采样器将是MH采样器。

Metroplis-in-Gibbs采样

目标是从中取样。

MH采样器的工作方式如下:

开始采样。

让我们假设将提议分布的方差设置为某个常数。

我们计算在上一次绘制时评估的非标准化密度与当前提议的比率: 

如果该比率大于1,则当前提议分布的密度高于先前值的密度。因此,我们“接受”了提议并确定了。然后,我们使用以提议为中心的提议分布重复步骤2-4  ,然后生成新提议。如果该比率小于1,则当前建议值的密度低于先前建议。

因此,总是接受产生更高条件的后验评估的提议。但是,有时仅接受具有较低密度评估的提议-提议的相对密度评估越低,其接受的可能性就越低。

经过多次迭代,从后验的高密度区域开始的抽样被接受,并且被接受的序列“爬升”到高密度区域。一旦序列到达此高密度区域,它将趋于保持在那里。因此,这也类似于模拟退火。

这种表示法很容易扩展到我们的4维示例:提议分布现在是4维多元高斯模型。代替标量方差参数,我们有一个协方差矩阵。因此,我们的建议是系数的向量。从这个意义上讲,我们运行的是Gibbs –使用MH每次迭代绘制整个系数。

跳跃分布的方差是重要的参数。如果方差太小,则当前提议可能会非常接近最后一个值,因此![R也很可能接近1。因此,我们会非常频繁地接受,但由于接受的值彼此之间非常接近,因此我们会攀升至较高在许多次迭代中慢慢降低密度区域。如果方差太大,则序列到达高密度区域后可能无法保留在该区域。

许多“自适应” MH方法是此处描述的基本算法的变体,但包括调整周期以找到产生最佳接受率的跳跃分布方差。

MH中计算量最大的部分是密度评估。对于每个Gibbs迭代,我们必须两次评估4维密度。

尽管很容易扩展到高维度,但性能本身在高维度上会变差。

结果

这是我们感兴趣的4个参数的MCMC链。红线表示真实值。

  • 发表于:
  • 原文链接https://kuaibao.qq.com/s/20221121A058MM00?refer=cp_1026
  • 腾讯「腾讯云开发者社区」是腾讯内容开放平台帐号(企鹅号)传播渠道之一,根据《腾讯内容开放平台服务协议》转载发布内容。
  • 如有侵权,请联系 cloudcommunity@tencent.com 删除。

相关快讯

扫码

添加站长 进交流群

领取专属 10元无门槛券

私享最新 技术干货

扫码加入开发者社群
领券