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

61次阅读

共计 3573 个字符,预计需要花费 9 分钟才能阅读完成。

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

最近咱们被客户要求撰写对于 Gibbs 抽样的钻研报告,包含一些图形和统计输入。
贝叶斯剖析的许多介绍都应用了绝对简略的教学实例(例如,依据伯努利数据给出胜利概率的推理)。尽管这很好地介绍了贝叶斯原理,然而这些准则的扩大并不是含糊其辞的

这篇文章将概述这些原理如何扩大到简略的线性回归。我将导出感兴趣参数的后验条件散布,给出用于实现 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 日以前的平局。


点击题目查阅往期内容

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

左右滑动查看更多

01

02

03

04

条件后验散布

要应用 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 办法依赖于每个参数的条件后验散布的程序绘制。这是一个容易辨认的已知的散布。对于斜率和截距项,咱们决定用网格办法来躲避代数。


点击文末 “浏览原文”

获取全文残缺材料。

本文选自《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 采样用于回归的贝叶斯估

正文完
 0