全文下载链接：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<-1000burnin<-101phi<-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采样用于回归的贝叶斯估