R语言Gibbs抽样的贝叶斯简单线性回归仿真分析

原文连接：http://tecdat.cn/?p=4612

贝叶斯分析的许多介绍都使用了相对简单的教学实例（例如根据伯努利数据给出成功几率的推理）。虽然这能够很好地介绍贝叶斯原理，可是将这些原理扩展到回归并非直接的。

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

贝叶斯模型

假设咱们观察数据​

对于咱们的模型​是

​

有趣的是推断​和​。若是咱们将正态先验放在系数上，将反伽玛先验放在方差项上，则此数据的完整贝叶斯模型能够写为：

​

假设超参数​是已知的，后验能够被写入到一个比例常数，

​

括号中的项是是数据的联合分布或几率。其余项包括参数的联合先验分布。

​

R代码从指定的“真实”参数模型生成数据。咱们稍后将用这个数据估计贝叶斯回归模型来检查是否能够恢复这些真实的参数。

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个样本是从后验分布中抽取的。这些样本不是独立的。每一个步骤都取决于先前的位置。

条件后验分布

要使用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采样器。

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模型实现