关于数据挖掘:拓端tecdatR语言中Gibbs抽样的Bayesian贝叶斯简单线性回归

44次阅读

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

原文链接: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 模型实现

正文完
 0