拓端tecdat|R语言使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析
原文链接：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每次迭代绘制整个系数。
- 跳跃分布的方差是重要的参数。如果方差太小，则当前提议可能会非常接近最后一个值，因此也很可能接近1。因此，我们会非常频繁地接受，但由于接受的值彼此之间非常接近，因此我们会攀升至较高在许多次迭代中慢慢降低密度区域。如果方差太大，则序列到达高密度区域后可能无法保留在该区域。
- 许多“自适应” MH方法是此处描述的基本算法的变体，但包括调整周期以找到产生最佳接受率的跳跃分布方差。
- MH中计算量最大的部分是密度评估。对于每个Gibbs迭代，我们必须两次评估4维密度。
- 尽管很容易扩展到高维度，但性能本身在高维度上会变差。
结果
这是我们感兴趣的4个参数的MCMC链。红线表示真实值。
profvis(expr = {
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')
abline(h=.7,col='red')
plot(gibbs_res[,3],type='l',xlab='MCMC Iterations',ylab=c('Coefficient Draw'),main='Age2')
abline(h=1.1,col='red')
plot(gibbs_res[,4],type='l',xlab='MCMC Iterations',ylab=c('Coefficient Draw'),main='Treatment')
abline(h=1.1,col='red')# 计算后验分布和置信区间
post_burn_trim<-gibbs_res[seq(1000,gibbs_iter,100),]
colMeans(post_burn_trim)
apply(post_burn_trim, 2, quantile, p=c(.025,.975))
有一些改进的空间：
- 接受率只有18％，我本可以调整跳跃分布协方差矩阵来获得更好的接受率。
- 我认为更多的迭代肯定会在这里有所帮助。这些链看起来不错，但仍然是自相关的。
关于贝叶斯范式的好处是，所有推断都是使用后验分布完成的。现在，系数估计值是对数化，但是如果我们需要比值，则只需对后验取幂。如果我们想要对比值进行区间估计，那么我们就可以获取指数后验的2.5％和97.5％。
下面是使用R分析，显示了这一点。for循环运行Gibbs迭代。在每个Gibbs迭代中，我都调用函数使用MH从参数向量的条件后验中得出图形。
我们看到子例程log_cond（）是MH运行中的瓶颈。此函数是beta的对数条件后验密度。
最受欢迎的见解
1.使用R语言进行METROPLIS-IN-GIBBS采样和MCMC运行
2.R语言中的Stan概率编程MCMC采样的贝叶斯模型
3.R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样
4.R语言BUGS JAGS贝叶斯分析 马尔科夫链蒙特卡洛方法（MCMC）采样
5.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归
6.R语言Gibbs抽样的贝叶斯简单线性回归仿真分析
7.R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数
8.R语言使用Metropolis- Hasting抽样算法进行逻辑回归
9.R语言中基于混合数据抽样(MIDAS)回归的HAR-RV模型预测GDP增长
拓端tecdat|R语言使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析相关推荐
- 拓端tecdat|R语言逻辑回归（Logistic回归）模型分类预测病人冠心病风险
最近我们被客户要求撰写关于冠心病风险的研究报告,包括一些图形和统计输出. 相关视频:R语言逻辑回归(Logistic回归)模型分类预测病人冠心病风险 逻辑回归Logistic模型原理和R语言分类预测冠 ...
- 拓端tecdat|R语言用LOESS(局部加权回归)季节趋势分解（STL）进行时间序列异常检测
最近我们被客户要求撰写关于LOESS(局部加权回归)的研究报告,包括一些图形和统计输出. 这篇文章描述了一种对涉及季节性和趋势成分的时间序列的中点进行建模的方法.我们将对一种叫做STL的算法进行研究, ...
- 拓端tecdat|R语言向量误差修正模型 (VECMs)分析长期利率和通胀率影响关系
最近我们被客户要求撰写关于向量误差修正模型的研究报告,包括一些图形和统计输出. 向量自回归模型估计的先决条件之一是被分析的时间序列是平稳的.但是,经济理论认为,经济变量之间在水平上存在着均衡关系,可以 ...
- 拓端tecdat|R语言线性回归和时间序列分析北京房价影响因素可视化案例
最近我们被客户要求撰写关于北京房价影响因素的研究报告,包括一些图形和统计输出. 目的 房价有关的数据可能反映了中国近年来的变化: 人们得到更多的资源(薪水),期望有更好的房子 人口众多 独生子女政策: ...
- 拓端tecdat荣获掘金社区入驻新人奖
2021年7月,由掘金发起了"入驻成长礼"颁奖活动.本次活动邀请到知名开发者.服务机构代表等业界人士. 据了解,掘金社区"新入驻创作者礼"主要对已经积累了一定历 ...
- 拓端tecdat荣获2022年度51CTO博主之星
相信技术,传递价值,这是51CTO每一个技术创作者的动力与信念,2022 年度,拓端tecdat 作为新锐的数据分析咨询公司,在51CTO平台上,不断的输出优质的技术文章,分享前沿创新技术,输出最佳生 ...
- 拓端tecdat|bilibili视频流量数据潜望镜
最近我们被客户要求撰写关于bilibili视频流量的研究报告,包括一些图形和统计输出. 最新研究表明,中国有超过7亿人在观看在线视频内容.Bilibili,被称为哔哩哔哩或简称为B站,是中国大陆第二个 ...
- R语言Bootstrap、百分位Bootstrap法抽样参数估计置信区间分析通勤时间和学生锻炼数据
最近我们被客户要求撰写关于Bootstrap的研究报告,包括一些图形和统计输出. 本文展示了如何使用 R 构建Bootstrap自举置信区间的示例.还强调了 R 包 ggplot2 用于图形的用途.但 ...
- R语言——（三）、随机数与抽样模拟
文章目录 前序r d p q 一.一元随机数的产生 1.均匀分布随机数runif 2.正态分布随机数的产生rnorm 3.指数分布随机数产生rexp 4.二项分布随机数的产生rbinom 二.多元随机 ...
- R语言广义加性模型GAMs分析温度、臭氧环境数据绘制偏回归图与偏残差图
最近我们被客户要求撰写关于广义加性模型的研究报告,包括一些图形和统计输出. 视频:R语言广义相加模型(GAM)在电力负荷预测中的应用 拓端tecdat:R语言广义相加模型(GAM)在电力负荷预测中的应 ...
最新文章
- Unix 和 Linux
- git的使用学习（三）时光机穿梭
- Zabbix linux agent 安装
- mysql某个表被行锁了_一文搞懂MySQL行锁、表锁、间隙锁详解
- [译] 如何用 Android vitals 解决应用程序的质量问题
- android flash游戏播放器代码,新flashplayer游戏播放器apk-新Flash游戏播放器(NEW Server Mode)预约V4.2安卓版-西西软件园...
- 优盘安装红帽linux系统,RedHat Linux系统U盘安装图文教程
- 用友凭证打印故障解决
- The TARGETDIR variable must be provided when invoking this installer的解决方案
- VS2019 后面有“::”的名称一定是类名或命名空间名 解决办法汇总
- MySQL 为什么需要 redo log？
- GPS公司内部管理系统
- 图片裁剪_如何裁剪某张图片的某一部分
- enumerate()说明
- python画图，坐标轴问题
- Android开发之监听手机来电
- raid管理软件MegaRAID Storage Manager(MSM)的安装和使用
- 常用的计算机技术博客,计算机常用命令！
- Python与简单网络爬虫的编写
- Arduino TFT_eSPI库来驱动SPI接口的LCD显示文字详解