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

于 2022-11-21 16:16:02 发布
阅读量859 收藏 19
点赞数 2
文章标签： r语言 开发语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/tecdat/article/details/128090481
版权

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

对于许多模型，例如逻辑模型，没有共轭先验分布。因此，吉布斯采样不适用（点击文末“阅读原文”获取完整代码数据）。

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

相关视频

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

模型

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

4f48dac67fb4306634db966ebf309492.png

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

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

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

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

非规范条件后验

让我们看一下该模型的（非标准化）条件后验。
75c0481293d4c4d3942955dc325b0a92.png
此条件分布不是已知分布，因此我们不能简单地使用Gibbs从中进行采样。相反，在每个gibbs迭代中，我们需要另一个采样步骤来从该条件后验中提取。第二个采样器将是MH采样器。

Metroplis-in-Gibbs采样

目标是从中取样cc9aa2c6bf4c847c2ffe27b5428319fc.png。

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

  1. 开始采样。

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

  3. 我们计算在上一次绘制时评估的非标准化密度与当前提议的比率：   4891ac4bc1261042b7466e6ee5c6065a.png

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

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

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

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

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

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

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

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

结果

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

for(i in 2:gibbs_iter){
  # 来自 phi 后验分布的样本
  gibbs_res[i,p+1] <- rcond_post_phi(gibbs_res[i-1,1:p], 
                                     alpha, gamma, lambda, p)
  #  来自beta后验分布的样本 （ 使用 Mh ）
  mh_draw <- rcond_post_beta_mh(gibbs_res[i-1,1:p], gibbs_res[i,p+1], 
                                lambda, X, Y, mh_trials=5, jump_v=.01)
}
par(mfrow=c(2,2))
plot(gibbs_res[,1],type='l',xlab='MCMC Iterations',ylab=c('Coefficient Draw'),
     main='Intercept')
abline(h=-1,col='red')
plot(gibbs_res[,2],type='l',xlab='MCMC Iterations',ylab=c('Coefficient Draw'),
     main='Age1')
# 计算后验分布和置信区间
post_burn_trim<-gibbs_res[seq(1000,gibbs_iter,100),]
colMeans(post_burn_trim)
apply(post_burn_trim, 2, quantile, p=c(.025,.975))

d13d5b3220636f874d219606de7b2ca1.png

点击标题查阅往期内容

d2232dd0c958cfd2156daa83ae75371a.jpeg

R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

outside_default.png

左右滑动查看更多

outside_default.png

01

df72dfddc2a41aef5d1e154ad14f05d3.png

02

44f1b16070348e52c216585acdf1d3d7.png

03

e9154cefdf30722371972aae8f62bba7.png

04

47a670c750adb7845a47e4fc73895f34.png

有一些改进的空间：

  • 接受率只有18％，我本可以调整跳跃分布协方差矩阵来获得更好的接受率。

  • 我认为更多的迭代肯定会在这里有所帮助。这些链看起来不错，但仍然是自相关的。

关于贝叶斯范式的好处是，所有推断都是使用后验分布完成的。现在，系数估计值是对数化，但是如果我们需要比值，则只需对后验取幂。如果我们想要对比值进行区间估计，那么我们就可以获取指数后验的2.5％和97.5％。

下面是使用R分析，显示了这一点。for循环运行Gibbs迭代。在每个Gibbs迭代中，我都调用函数使用MH从参数向量的条件后验中得出图形。

e1198bdc64b71307b46a7320225e727c.png

我们看到子例程log_cond（）是MH运行中的瓶颈。此函数是beta的对数条件后验密度。

8ed96c75a9db198b37df742c3331cbae.png

341fbf95376d5a522989a901130812fa.jpeg

点击文末“阅读原文”

获取全文完整代码数据资料。

本文选自《使用R语言进行Metroplis-in-Gibbs采样和MCMC运行分析》。

点击标题查阅往期内容

R语言BUGS/JAGS贝叶斯分析: 马尔科夫链蒙特卡洛方法（MCMC）采样

R语言MCMC:Metropolis-Hastings采样用于回归的贝叶斯估计

Python用MCMC马尔科夫链蒙特卡洛、拒绝抽样和Metropolis-Hastings采样算法

R语言贝叶斯METROPOLIS-HASTINGS GIBBS 吉布斯采样器估计变点指数分布分析泊松过程车站等待时间

R语言马尔可夫MCMC中的METROPOLIS HASTINGS，MH算法抽样（采样）法可视化实例

python贝叶斯随机过程：马尔可夫链Markov-Chain，MC和Metropolis-Hastings，MH采样算法可视化

Python贝叶斯推断Metropolis-Hastings（M-H）MCMC采样算法的实现

Metropolis Hastings采样和贝叶斯泊松回归Poisson模型

Matlab用BUGS马尔可夫区制转换Markov switching随机波动率模型、序列蒙特卡罗SMC、M H采样分析时间序列

R语言RSTAN MCMC：NUTS采样算法用LASSO 构建贝叶斯线性回归模型分析职业声望数据

R语言BUGS序列蒙特卡罗SMC、马尔可夫转换随机波动率SV模型、粒子滤波、Metropolis Hasting采样时间序列分析

R语言Metropolis Hastings采样和贝叶斯泊松回归Poisson模型

R语言贝叶斯MCMC：用rstan建立线性回归模型分析汽车数据和可视化诊断

R语言贝叶斯MCMC：GLM逻辑回归、Rstan线性回归、Metropolis Hastings与Gibbs采样算法实例

R语言贝叶斯Poisson泊松-正态分布模型分析职业足球比赛进球数

R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数

R语言逻辑回归、Naive Bayes贝叶斯、决策树、随机森林算法预测心脏病

R语言中贝叶斯网络（BN）、动态贝叶斯网络、线性模型分析错颌畸形数据

R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

Python贝叶斯回归分析住房负担能力数据集

R语言实现贝叶斯分位数回归、lasso和自适应lasso贝叶斯分位数回归分析

Python用PyMC3实现贝叶斯线性回归模型

R语言用WinBUGS 软件对学术能力测验建立层次（分层）贝叶斯模型

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

R语言和STAN,JAGS：用RSTAN,RJAG建立贝叶斯多元线性回归预测选举数据

R语言基于copula的贝叶斯分层混合模型的诊断准确性研究

R语言贝叶斯线性回归和多元线性回归构建工资预测模型

R语言贝叶斯推断与MCMC：实现Metropolis-Hastings 采样算法示例

R语言stan进行基于贝叶斯推断的回归模型

R语言中RStan贝叶斯层次模型分析示例

R语言使用Metropolis-Hastings采样算法自适应贝叶斯估计与可视化

R语言随机搜索变量选择SSVS估计贝叶斯向量自回归（BVAR）模型

WinBUGS对多元随机波动率模型：贝叶斯估计与模型比较

R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样

R语言贝叶斯推断与MCMC：实现Metropolis-Hastings 采样算法示例

R语言使用Metropolis-Hastings采样算法自适应贝叶斯估计与可视化

视频：R语言中的Stan概率编程MCMC采样的贝叶斯模型

R语言MCMC:Metropolis-Hastings采样用于回归的贝叶斯估计

8374564e19ed348ca03194e38495eb53.png

221afb568f66eaa0f5e8d329a77969b4.jpeg

03fd88d002402da5cfa8fe8d7c6c8876.png

拓端研究室TRL
关注
  • 2
    点赞
  • 踩
  • 19
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
博客
【专题】2024中国零售业发展报告-拥抱“性价比时代”报告合集PDF分享（附原数据表）...
04-26 9
原文链接：https://tecdat.cn/?p=36050在消费升级与降级交织的复杂市场背景下，中国零售市场正迈向性价比主导新阶段，高价值追求成核心竞争力。国家统计局数据显示，一季度GDP和居民收入增长，但社零和零售增速放缓，消费信心待提振（查看文末了解专题报告合集PDF版本免费获取方式）。相关视频2023年快消品零售市场销售额微降，但内容电商和便利店等新兴渠道表现亮眼，特别是内容电商平台通过...
博客
【专题】中国美妆行业白皮书报告合集PDF分享（附原数据表）
04-25 2
原文链接：https://tecdat.cn/?p=36047在过去的十年里，中国美妆市场持续增长，展现出巨大的活力和潜力。2023年，护肤与彩妆两大品类市场规模分别突破1700.9亿元和1193.5亿元，引领市场潮流（查看文末了解专题报告合集PDF版本免费获取方式）。国货美妆品牌在面膜、精华等领域表现出色，彩妆品类亦迅猛增长，底妆增速尤为亮眼，高达22.4%。这些品牌凭借高品质、中国特色成分及对...
博客
【专题】2023-2024年游客满意度调查报告合集PDF分享（附原数据表）
04-24 15
原文链接：https://tecdat.cn/?p=36043游客不仅是旅游服务质量的最终评判者，还是旅游业的界定者、城市旅游的评价者，更是塑造国家级和世界级旅游城市的关键力量（查看文末了解专题报告合集PDF版本免费获取方式）。他们的满意度直接揭示了旅游发展的动力与方向，为旅游城市的建设指明路径。2011年，我国荣获联合国世界旅游组织政策管理创新一等奖，开创先河。阅读原文，获取专题报告合集全文，解...
博客
【专题】2023-2024年二手车市场消费需求洞察报告合集PDF分享（附原数据表）
04-23 16
原文链接 :https://tecdat.cn/?p=360212023年，乘用车二手车市场展现出了蓬勃的发展态势，交易量攀升至1478万辆，同比增长高达15%，创下了近五年的新纪录（查看文末了解专题报告合集PDF版本免费获取方式）。相关视频在这一繁荣景象的背后，以旧换新政策的助力功不可没，它极大地激发了消费者的换购热情，为二手车市场注入了源源不断的新活力，推动了市场的持续繁荣与发展。阅读原文，获...
博客
【视频】N-Gram、逻辑回归反欺诈模型文本分析招聘网站欺诈可视化讲解|附数据代码...
04-23 676
原文链接：https://tecdat.cn/?p=36028分析师：Zhuo Chen随着互联网的快速发展，招聘网站已成为求职者与雇主之间的重要桥梁。然而，随之而来的欺诈行为也日益猖獗，给求职者带来了极大的困扰和风险（点击文末“阅读原文”获取完整代码数据）。视频因此，如何帮助客户有效地识别和防范招聘网站上的欺诈行为，已成为一个亟待解决的问题。逻辑回归模型作为一种强大的分类工具，在识别欺诈行为方面...
博客
Python遗传算法GA对长短期记忆LSTM深度学习模型超参数调优分析司机数据|附数据代码...
04-22 13
全文链接：https://tecdat.cn/?p=36004随着大数据时代的来临，深度学习技术在各个领域中得到了广泛的应用。长短期记忆（LSTM）网络作为深度学习领域中的一种重要模型，因其对序列数据的强大处理能力，在自然语言处理、时间序列预测等领域中取得了显著的成果（点击文末“阅读原文”获取完整代码数据）。相关视频然而，LSTM模型的性能在很大程度上取决于其超参数的设置。因此，如何有效地帮助客户...
博客
【专题】2024年中国人工智能产业研究报告合集PDF分享（附原数据表）
04-22 26
原文链接:https://tecdat.cn/?p=36000纵观人类进化史，制造和使用工具无疑是我们发展进程中的一大里程碑。而在当今时代，如何更好地利用AI工具已成为社会各个层面——从产业应用到日常生活、学习工作——热议的焦点（查看文末了解专题报告合集PDF版本免费获取方式）。相关视频特别是随着大模型和生成式AI技术的崛起，它们强大的数据处理、学习泛化及内容生成能力，极大地推动了各行业智能化进程...
博客
文本挖掘分析多元应用：新能源汽车股市、英国封锁、疫情旅游与舆情分析
04-19 536
全文链接：https://tecdat.cn/?p=35975分析师：Yaning Wu，Zexi Jin，Chang Zhang，Wei Zuo在当今信息爆炸的时代，文本分析作为一种重要的数据处理方法，已经广泛应用于各个领域的研究中（点击文末“阅读原文”了解更多）。相关视频本文旨在通过运用文本分析技术，帮助客户深入探究新能源汽车股市涨跌影响、英国全国性封锁对零售配送产业的影响、疫情下的旅游微博数...
博客
【专题】2024业务视角的中国企业AI+学习发展报告合集PDF分享（附原数据表）
04-19 64
原文链接：https://tecdat.cn/?p=35970学习型组织的数字化转型已成为企业竞相追逐的焦点。在数字战略的推动下，人才升级显得尤为关键，而数字化学习则成为了加速人才培养的强有力工具（查看文末了解专题报告合集PDF版本免费获取方式）。相关视频AI技术的广泛运用，为企业的学习型组织插上了数字化转型的翅膀，助力其腾飞。这份报告合集深入解析了AI时代企业的学习发展，全面剖析了AI在业务赋能...
博客
【专题】2024新能源及储能参与电力市场交易白皮书报告合集PDF分享（附原数据表）...
04-18 49
原文链接：https://tecdat.cn/?p=359582019至2023年，我国新能源市场化交易电量持续增长，2023年更是达到6,845亿千瓦时，占新能源发电总量的47.3%（查看文末了解专题报告合集PDF版本免费获取方式）。相关视频同年，国家电网公司绿电结算电量跃升至576亿千瓦时，绿证交易也激增15倍，达到2,364万张。阅读原文，获取专题报告合集全文，解锁文末204份电力相关行业研...
博客
【专题】2024年中国跨境进口电商B2C行业消费趋势白皮书报告合集PDF分享（附原数据表）...
04-17 11
原文链接：https://tecdat.cn/?p=35940报告合集指出，随着中国经济开放程度的加深和全球化的推进，跨境进口电商行业规模持续扩大，抖音电商全球购等平台成为国际品牌进入中国市场的重要跳板（查看文末了解专题报告合集PDF版本免费获取方式）。相关视频消费者对品质生活的追求日益提升，对产品的专业性和真实性要求更加严格，更倾向于兴趣导向的购物决策。阅读原文，获取专题报告合集全文，解锁文末2...
博客
【专题】2024中国生成式AI开发者洞察报告合集PDF分享（附原数据表）
04-16 84
原文链接：https://tecdat.cn/?p=35908生成式人工智能正成为产业创新的关键力量，与实体经济的深度融合将推动新质产业的蓬勃发展（查看文末了解专题报告合集PDF版本免费获取方式）。相关视频本报告合集重点关注那些在各行各业致力于生成式AI技术研发或利用AI工具高效工作，并具备相应专业技术能力的开发者。这些开发者不仅包括传统的AI程序员，还涵盖了数据分析师、运营人员、市场专员及测试工...
博客
【视频】R语言支持向量回归SVR预测水位实例讲解|附代码数据
04-16 711
全文链接：https://tecdat.cn/?p=35914分析师：Miaoqiao Wang当我们面对样本需要建立相应模型时，使用传统统计方法建立模型需要大量的样本数据，只有在样本量足够大时，该模型才具有一定的可靠性，而实际实验中，不一定每次实验都拥有足够大的样本，甚至是小样本，这时使用传统统计方法来建立出的模型，在可靠性方面就存在一定的局限，难以达到理想的效果（点击文末“阅读原文”获取完整代...
博客
R语言Stan贝叶斯回归置信区间后验分布可视化模型检验|附数据代码
04-15 16
全文链接：https://tecdat.cn/?p=35899贝叶斯回归是一种统计方法，它使用贝叶斯定理来估计回归模型的参数。与传统的频率派回归方法不同，贝叶斯回归提供了参数的后验分布，而不仅仅是点估计。这意味着我们可以得到参数的不确定性度量，而不仅仅是单一的估计值（点击文末“阅读原文”获取完整代码数据）。相关视频R语言用于拟合贝叶斯广义线性模型。我们可以方便地帮助客户拟合贝叶斯线性回归、逻辑回归...
博客
【专题】中国纯电新能源汽车-市场发展和用车报告2024年报告合集PDF分享（附原数据表）...
04-15 55
原文链接：https://tecdat.cn/?p=358782023年，纯电车型在新能源市场中占据高达七成的市场份额，尽管技术挑战仍然存在。插混及增程车型在续航里程和驾驶体验上取得显著进步，但纯电车型仍占据主导地位（查看文末了解专题报告合集PDF版本免费获取方式）。相关视频预计未来几年，插混及增程车型市场份额将持续攀升，为市场增添更多竞争与选择。阅读原文，获取专题报告合集全文，解锁文末336份新...
博客
数据分享|R语言机器学习预测案例合集：众筹平台、机票折扣、糖尿病患者、员工满意度...
04-12 670
全文链接：https://tecdat.cn/?p=35835分析师：Yujie Zhou，Weixin Hu，Xiaoyan Xu，Xu Tian在数据驱动的当代社会，机器学习已成为揭示复杂现象、预测未来趋势的重要工具。特别是在商业决策、健康管理、交通出行等多个领域，机器学习技术的应用日益广泛（点击文末“阅读原文”获取完整数据）。相关视频本文集合了多个基于R语言、python的机器学习预测案例，...
博客
【专题】2024休闲食品白皮书报告合集PDF分享（附原数据表）
04-12 47
原文链接：https://tecdat.cn/?p=35814根据艾媒咨询发布的报告，2023年中国零食市场规模已达到惊人的11654亿元，整体市场展现出稳健的增长态势（查看文末了解专题报告合集PDF版本免费获取方式）。在这一庞大的万亿级消费市场中，休闲食品的消费习惯以及销售渠道成为我们深入剖析的两大核心要素。阅读原文，获取专题报告合集全文，解锁文末110份休闲食品相关行业研究报告。同时，根据企查...
博客
R语言、WEKA从数据角度看员工为什么离职？
04-11 11
原文链接：http://tecdat.cn/?p=3986马云说：“员工离职的原因总是只有两个：钱，没有到位;心委屈了。”（点击文末“阅读原文”获取完整代码数据）。相关视频现在很多老板都抱怨说，年轻人的流动率太高了，员工觉得老板的钱太少了，最后还是多指责。为什么我们最好和最有经验的员工过早离职？到底如何解决这个困境？拓端数据tecdat使用数据分析员工离开的原因，希望能从中找到线索。我们搜集的数据...
博客
【专题】2024TikTok出海达人营销白皮书报告合集PDF分享（附原数据表）
04-11 12
原文链接：https://tecdat.cn/?p=35781自互联网时代崛起以来，达人的影响力已经日益增强，遍及全球各地（查看文末了解专题报告合集PDF版本免费获取方式）。相关视频事实上，不论是在世界的哪个角落，也不管文化背景的迥异，每一个消费者对于达人的需求和喜爱都是相似的。这是因为全世界的人们都怀揣着相同的愿景，那就是追求更加美好的生活，并期望在不同的生活场景中，都能有达人为他们提供丰富多样...
博客
【专题】2024年中国游戏营销趋势报告合集PDF分享（附原数据表）
04-10 51
原文链接：https://tecdat.cn/?p=357452023年，全球游戏行业表现卓越，不仅用户规模扩大至33.81亿，行业营收也攀升至1.35万亿人民币，呈现出强劲的增长态势（查看文末了解专题报告合集PDF版本免费获取方式）。相关视频然而，与此同时，全球游戏创业公司在风险投资上的大幅缩减也揭示了行业面临的某些挑战。阅读原文，获取专题报告合集全文，解锁文末104份游戏相关行业研究报告。在另...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值