R语言中Gibbs抽样的Bayesian简单线性回归

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

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日以前的平局。

条件后验分布

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)))

仿真结果

现在我们可以从每个参数的条件后验进行采样，我们可以实现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方法依赖于每个参数的条件后验分布的顺序绘制。这是一个容易识别的已知的分布。对于斜率和截距项，我们决定用网格方法来规避代数。