R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

ByteVoyant
• 阅读 591

全文下载链接:http://tecdat.cn/?p=4612

最近我们被客户要求撰写关于贝叶斯简单线性回归的研究报告,包括一些图形和统计输出。

贝叶斯分析的许多介绍都使用了相对简单的教学实例(例如,根据伯努利数据给出成功概率的推理)。虽然这很好地介绍了贝叶斯原理,但是这些原则的扩展并不是直截了当的

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

贝叶斯模型

假设我们观察数据

对于R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据我们的模型是

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

有兴趣的是作出推论

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

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

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

假设超参数

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据是已知的,后面可以写成一个常数的比例,

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

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

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

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

tphi<-rinvgamma(1, shape=a, rate=g)
tb0<-rnorm(1, m0, sqrt(t0) )
tb1<-rnorm(1, m1, sqrt(t1) )
tphi; tb0; tb1;

y<-rnorm(n, tb0 + tb1*x, sqrt(tphi))

吉布斯采样器

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

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

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


点击标题查阅往期内容

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

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

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

左右滑动查看更多

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

01

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

02

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

03

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

04

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

条件后验分布

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

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

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

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

条件后验:

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

同样的,

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

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

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

条件后验R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据不那么容易识别。但是如果我们愿意使用网格方法,我们并不需要经过任何代数。

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据考虑网格方法。网格方法是非常暴力的方式(在我看来)从其条件后验分布进行抽样。这个条件分布只是一个函数R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据。所以我们可以评估一定的密度R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据值。在R表示法中,这可以是grid = seq(-10,10,by = .001)。这个序列是点的“网格”。

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

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


  for(i in 1:length(p) ){
    p[i]<- (-(1/(2*phi))*sum( (y - (grid[i]+b1*x))^2 ))  + ( -(1/(2*t0))*(grid[i] - m0)^2)
  }
  

  draw<-sample(grid, size = 1, prob = exp(1-p/max(p)))

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

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

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

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

因为它来自已知的分布R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

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

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

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

仿真结果

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

iter<-1000
burnin<-101
phi<-b0<-b1<-numeric(iter)
phi[1]<-b0[1]<-b1[1]<-6

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

z <- kde2d(b0, b1, n=50)
plot(b0,b1, pch=19, cex=.4)
contour(z, drawlabels=FALSE, nlevels=10, col='red', add=TRUE)

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

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


R语言Gibbs抽样的贝叶斯简单线性回归仿真分析|附代码数据

点击文末 “阅读原文”

获取全文完整资料。

本文选自《R语言Gibbs抽样的贝叶斯简单线性回归仿真分析》。

点击标题查阅往期内容

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采样用于回归的贝叶斯估计

点赞
收藏
评论区
推荐文章
blmius blmius
3年前
MySQL:[Err] 1292 - Incorrect datetime value: ‘0000-00-00 00:00:00‘ for column ‘CREATE_TIME‘ at row 1
文章目录问题用navicat导入数据时,报错:原因这是因为当前的MySQL不支持datetime为0的情况。解决修改sql\mode:sql\mode:SQLMode定义了MySQL应支持的SQL语法、数据校验等,这样可以更容易地在不同的环境中使用MySQL。全局s
皕杰报表之UUID
​在我们用皕杰报表工具设计填报报表时,如何在新增行里自动增加id呢?能新增整数排序id吗?目前可以在新增行里自动增加id,但只能用uuid函数增加UUID编码,不能新增整数排序id。uuid函数说明:获取一个UUID,可以在填报表中用来创建数据ID语法:uuid()或uuid(sep)参数说明:sep布尔值,生成的uuid中是否包含分隔符'',缺省为
Wesley13 Wesley13
3年前
PPDB:今晚老齐直播
【今晚老齐直播】今晚(本周三晚)20:0021:00小白开始“用”飞桨(https://www.oschina.net/action/visit/ad?id1185)由PPDE(飞桨(https://www.oschina.net/action/visit/ad?id1185)开发者专家计划)成员老齐,为深度学习小白指点迷津。
Jacquelyn38 Jacquelyn38
4年前
2020年前端实用代码段,为你的工作保驾护航
有空的时候,自己总结了几个代码段,在开发中也经常使用,谢谢。1、使用解构获取json数据let jsonData  id: 1,status: "OK",data: 'a', 'b';let  id, status, data: number   jsonData;console.log(id, status, number )
Stella981 Stella981
3年前
PhoneGap设置Icon
参考:http://cordova.apache.org/docs/en/latest/config\_ref/images.html通过config.xml中的<icon标签来设置Icon<iconsrc"res/ios/icon.png"platform"ios"width"57"height"57"densi
Wesley13 Wesley13
3年前
mysql设置时区
mysql设置时区mysql\_query("SETtime\_zone'8:00'")ordie('时区设置失败,请联系管理员!');中国在东8区所以加8方法二:selectcount(user\_id)asdevice,CONVERT\_TZ(FROM\_UNIXTIME(reg\_time),'08:00','0
Wesley13 Wesley13
3年前
00:Java简单了解
浅谈Java之概述Java是SUN(StanfordUniversityNetwork),斯坦福大学网络公司)1995年推出的一门高级编程语言。Java是一种面向Internet的编程语言。随着Java技术在web方面的不断成熟,已经成为Web应用程序的首选开发语言。Java是简单易学,完全面向对象,安全可靠,与平台无关的编程语言。
Stella981 Stella981
3年前
Django中Admin中的一些参数配置
设置在列表中显示的字段,id为django模型默认的主键list_display('id','name','sex','profession','email','qq','phone','status','create_time')设置在列表可编辑字段list_editable
Wesley13 Wesley13
3年前
MySQL部分从库上面因为大量的临时表tmp_table造成慢查询
背景描述Time:20190124T00:08:14.70572408:00User@Host:@Id:Schema:sentrymetaLast_errno:0Killed:0Query_time:0.315758Lock_
Python进阶者 Python进阶者
1年前
Excel中这日期老是出来00:00:00,怎么用Pandas把这个去除
大家好,我是皮皮。一、前言前几天在Python白银交流群【上海新年人】问了一个Pandas数据筛选的问题。问题如下:这日期老是出来00:00:00,怎么把这个去除。二、实现过程后来【论草莓如何成为冻干莓】给了一个思路和代码如下:pd.toexcel之前把这
美凌格栋栋酱 美凌格栋栋酱
4个月前
Oracle 分组与拼接字符串同时使用
SELECTT.,ROWNUMIDFROM(SELECTT.EMPLID,T.NAME,T.BU,T.REALDEPART,T.FORMATDATE,SUM(T.S0)S0,MAX(UPDATETIME)CREATETIME,LISTAGG(TOCHAR(
ByteVoyant
ByteVoyant
Lv1
家在梦中何日到,春来江上几人还?
文章
3
粉丝
0
获赞
0