Metroplis-in-Gibbs采样

目标是从中取样。

MH采样器的工作方式如下：

这种表示法很容易扩展到我们的4维示例：提议分布现在是4维多元高斯模型。代替标量方差参数，我们有一个协方差矩阵。因此，我们的建议是系数的向量。从这个意义上讲，我们运行的是Gibbs –使用MH每次迭代绘制整个系数。

• 跳跃分布的方差是重要的参数。如果方差太小，则当前提议可能会非常接近最后一个值，因此也很可能接近1。因此，我们会非常频繁地接受，但由于接受的值彼此之间非常接近，因此我们会攀升至较高在许多次迭代中慢慢降低密度区域。如果方差太大，则序列到达高密度区域后可能无法保留在该区域。
• 许多“自适应” MH方法是此处描述的基本算法的变体，但包括调整周期以找到产生最佳接受率的跳跃分布方差。
• MH中计算量最大的部分是密度评估。对于每个Gibbs迭代，我们必须两次评估4维密度。
• 尽管很容易扩展到高维度，但性能本身在高维度上会变差。

结果

这是我们感兴趣的4个参数的MCMC链。红线表示真实值。

profvis(expr = {
for(i in 2:gibbs_iter){
  # 来自 phi 后验分布的样本
  gibbs_res[i,p+1] <- rcond_post_phi(gibbs_res[i-1,1:p], 
                                     alpha, gamma, lambda, p)
  #  来自beta后验分布的样本 （ 使用 Mh ）
  mh_draw <- rcond_post_beta_mh(gibbs_res[i-1,1:p], gibbs_res[i,p+1], 
                                lambda, X, Y, mh_trials=5, jump_v=.01)
}
})
par(mfrow=c(2,2))
plot(gibbs_res[,1],type='l',xlab='MCMC Iterations',ylab=c('Coefficient Draw'),
     main='Intercept')
abline(h=-1,col='red')
plot(gibbs_res[,2],type='l',xlab='MCMC Iterations',ylab=c('Coefficient Draw'),
     main='Age1')
abline(h=.7,col='red')
plot(gibbs_res[,3],type='l',xlab='MCMC Iterations',ylab=c('Coefficient Draw'),
     main='Age2')
abline(h=1.1,col='red')
plot(gibbs_res[,4],type='l',xlab='MCMC Iterations',ylab=c('Coefficient Draw'),
     main='Treatment')
abline(h=1.1,col='red')
# 计算后验分布和置信区间
post_burn_trim<-gibbs_res[seq(1000,gibbs_iter,100),]
colMeans(post_burn_trim)
apply(post_burn_trim, 2, quantile, p=c(.025,.975))

有一些改进的空间：

• 接受率只有18％，我本可以调整跳跃分布协方差矩阵来获得更好的接受率。
• 我认为更多的迭代肯定会在这里有所帮助。这些链看起来不错，但仍然是自相关的。

如果我们想要对比值进行区间估计，那么我们就可以获取指数后验的2.5％和97.5％。

下面是使用R分析，显示了这一点。for循环运行Gibbs迭代。在每个Gibbs迭代中，我都调用函数使用MH从参数向量的条件后验中得出图形。