拓端tecdat|R语言中Gibbs抽样的Bayesian贝叶斯简单线性回归
原文链接:http://tecdat.cn/?p=4612
原文出处:拓端数据部落公众号
贝叶斯分析的许多介绍都使用了相对简单的教学实例(例如,根据伯努利数据给出成功概率的推理)。虽然这很好地介绍了贝叶斯原理,但是这些原则的扩展并不是直截了当的。
这篇文章将概述这些原理如何扩展到简单的线性回归。我将导出感兴趣参数的后验条件分布,给出用于实现Gibbs采样器的R代码,并提出所谓的网格点方法。
贝叶斯模型
假设我们观察数据
对于我们的模型是
有兴趣的是作出推论
如果我们在方差项之前放置正态前向系数和反伽马,那么这个数据的完整贝叶斯模型可以写成:
假设超参数
是已知的,后面可以写成一个常数的比例,
括号中的术语是数据或可能性的联合分布。其他条款包括参数的联合先验分布(因为我们隐含地假设独立前,联合先验因素)。
伴随的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抽样算法。吉布斯采样是一种迭代算法,从每个感兴趣的参数的后验分布产生样本。它通过按照以下方式从每个参数的条件后面依次绘制:
可以看出,剩下的1,000个抽签是从后验分布中抽取的。这些样本不是独立的。绘制顺序是随机游走在后空间,空间中的每一步取决于前一个位置。通常还会使用间隔期(这里不做)。这个想法是,每一个平局可能依赖于以前的平局,但不能作为依赖于10日以前的平局。
条件后验分布
要使用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()中,我实际上评估了派生的条件后验分布的对数。然后,我通过从所有评估的最大值减去每个评估之前归一化,然后从对数刻度取回。
我们不需要使用网格方法来从条件的后面绘制。
因为它来自已知的分布
请注意,这种网格方法有一些缺点。
首先,这在计算上是复杂的。通过代数,希望得到一个已知的后验分布,从而在计算上更有效率。
其次,网格方法需要指定网格点的区域。如果条件后验在我们指定的[-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)
总结一下,我们首先推导了一个表达式,用于参数的联合分布。然后我们概述了从后面抽取样本的Gibbs算法。在这个过程中,我们认识到Gibbs方法依赖于每个参数的条件后验分布的顺序绘制。 这是一个容易识别的已知的分布。对于斜率和截距项,我们决定用网格方法来规避代数。
参考文献
1.matlab使用贝叶斯优化的深度学习
2.matlab贝叶斯隐马尔可夫hmm模型实现
3.R语言Gibbs抽样的贝叶斯简单线性回归仿真
4.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归
5.R语言中的Stan概率编程MCMC采样的贝叶斯模型
6.Python用PyMC3实现贝叶斯线性回归模型
7.R语言使用贝叶斯 层次模型进行空间数据分析
8.R语言随机搜索变量选择SSVS估计贝叶斯向量自回归(BVAR)模型
9.matlab贝叶斯隐马尔可夫hmm模型实现
拓端tecdat|R语言中Gibbs抽样的Bayesian贝叶斯简单线性回归相关推荐
- 拓端tecdat|R语言用LOESS(局部加权回归)季节趋势分解(STL)进行时间序列异常检测
最近我们被客户要求撰写关于LOESS(局部加权回归)的研究报告,包括一些图形和统计输出. 这篇文章描述了一种对涉及季节性和趋势成分的时间序列的中点进行建模的方法.我们将对一种叫做STL的算法进行研究, ...
- 拓端tecdat|R语言向量误差修正模型 (VECMs)分析长期利率和通胀率影响关系
最近我们被客户要求撰写关于向量误差修正模型的研究报告,包括一些图形和统计输出. 向量自回归模型估计的先决条件之一是被分析的时间序列是平稳的.但是,经济理论认为,经济变量之间在水平上存在着均衡关系,可以 ...
- 拓端tecdat|R语言线性回归和时间序列分析北京房价影响因素可视化案例
最近我们被客户要求撰写关于北京房价影响因素的研究报告,包括一些图形和统计输出. 目的 房价有关的数据可能反映了中国近年来的变化: 人们得到更多的资源(薪水),期望有更好的房子 人口众多 独生子女政策: ...
- 拓端tecdat|R语言逻辑回归(Logistic回归)模型分类预测病人冠心病风险
最近我们被客户要求撰写关于冠心病风险的研究报告,包括一些图形和统计输出. 相关视频:R语言逻辑回归(Logistic回归)模型分类预测病人冠心病风险 逻辑回归Logistic模型原理和R语言分类预测冠 ...
- R语言随机搜索变量选择SSVS估计贝叶斯向量自回归(BVAR)模型
介绍 最近我们被客户要求撰写关于向量自回归的研究报告,包括一些图形和统计输出.向量自回归(VAR)模型的一般缺点是,估计系数的数量与滞后的数量成比例地增加.因此,随着滞后次数的增加,每个参数可用的信息 ...
- 1071svm函数 r语言_如何利用R语言中的rpart函数建立决策树模型
决策树是根据若干输入变量的值构造出一个适合的模型,以此来预测输出变量的值,并用树形结构展示出来.决策树主要有两个类别:分类树和回归树.分类树主要针对离散的目标变量,回归树则针对连续的目标变量.R语言中 ...
- r语言中c函数错误,R语言中c()函数与paste()函数的区别说明
c()函数:将括号中的元素连接起来,并不创建向量 paste()函数:连接括号中的元素 例如 c(1, 2:4),结果为1 2 3 4 paste(1, 2:4),结果为"1 2" ...
- r语言中的shiny教程_如何使用Shiny在R中编写Web应用程序
r语言中的shiny教程 新年快乐! 这个月我忙于撰写一些较大的文章,因此请在接下来的几周内查找这些文章. 对于本月的Nooks和Crannies,我想简要指出一个我一直在用它进行自我教育的出色R库. ...
- R语言中if语句使用方法之超详细教程
在R语言中,if属于一种分支结构,即根据某个条件执行相关的语句.R中的if语句与else配合主要有3种结构. 单个if语句 if(cond) {expr} 其它语句 即当括弧中的cond条件为TRUE ...
- R语言中GCC编译的问题(续)
这篇文章承接R语言中GCC编译的问题,这篇文章主要解决我在Linux系统上安装"expm"出现的问题. 出现的问题 这个问题非常的有趣,因为我在两台服务器分别安装同一个包,其中一台 ...
最新文章
- 【通知】有三AI GPU平台上线新功能,GPU/CPU可灵活选择
- android P精简教程,华为EMUI 9.0发布:基于Android P打造 设置项精简10%
- python numeric_Python pandas.to_numeric函数方法的使用
- 使用Servlet实现用户注册
- python部分 + 数据库 + 网络编程
- ubuntu18.04如何安装mysql
- 中间滑动 头部底部固定_固定抗震弹性支座报价技术参数
- 基础知识—表达式与语句-语句
- obj是什么数据类型 python_Python入门级第一天
- JQuery插件,轻量级表单模型验证
- 真机调试 —— An unknown error occurred.
- 手机图片怎么免费转换成PDF格式?教程来了
- # Idea,2.5,软件安装,提示If you already have a 64-bit JDK installed ,defined a JAVA_HOME variable in Compu
- c# 爬网教程_Python Web爬网教程
- shopEx数据库错误,无法连接
- 北京大学肖臻老师《区块链技术与应用》公开课笔记-BTC
- C++ delete释放内存的本质
- java 画立体图形
- 一文读懂“个人经营收款码”和“个人收款码”的区别
- 如何求水平渐近线(例题讲解)