标签：基本算法 process 简单 有趣 title 慢慢 tin name 频繁
原文链接：http://tecdat.cn/?p=12200
对于许多模型，例如物流模型，没有共轭先验。因此，吉布斯采样不适用。
这篇文章展示了我们如何使用Metropolis-Hastings（MH）从每次Gibbs迭代中的非共轭条件后验对象中进行采样–比网格方法更好的替代方法。
我将说明该算法，给出一些R代码结果，然后分析R代码以识别MH算法中的瓶颈。
此示例的模拟数据是包含
?患者的横截面数据集。有一个二元结果，?一个二元治疗变量，?一个因子age。年龄是具有3个等级的分类变量。我用贝叶斯逻辑回归建模：
?
对于大都市吉布斯来说，这是一个相当现实的示例：
?，因此其他系数在某种意义上是令人讨厌的参数。我不会说这是一个“高维”设置，但肯定会给采样器带来压力。让我们看一下该模型的（非标准化）条件后验。我不会进行推导，但是它遵循我以前的帖子中使用的相同过程。
?
此条件分布不是已知分布，因此我们不能简单地使用Gibbs从中进行采样。相反，在每个gibbs迭代中，我们需要另一个采样步骤来从该条件后验中提取。第二个采样器将是MH采样器。
目标是从中取样
?。请注意，这是4维密度。
MH采样器的工作方式如下：
?
?。然后，我们使用以提案为中心的提案分布重复步骤2-4 ，然后生成新提案。如果该比率小于1，则当前建议值的密度低于先前建议。因此，总是接受产生更高条件的后验评估的提议。但是，有时仅接受具有较低密度评估的提案-提案的相对密度评估越低，其接受的可能性就越低。
经过多次迭代，从后验的高密度区域开始的抽样被接受，并且被接受的序列“爬升”到高密度区域。一旦序列到达此高密度区域，它将趋于保持在那里。因此，这也类似于模拟退火。
这种表示法很容易扩展到我们的4维示例：提案分布现在是4维多元高斯模型。代替标量方差参数，我们有一个协方差矩阵。因此，我们的建议是系数的向量。从这个意义上讲，我们运行的是Gibbs –使用MH每次迭代绘制整个系数块。
?也很可能接近1。因此，我们会非常频繁地接受，但由于接受的值彼此之间非常接近，因此我们会攀升至较高在许多次迭代中慢慢降低密度区域。如果方差太大，则序列到达高密度区域后可能无法保留在该区域。这是我们感兴趣的4个参数的MCMC链。红线表示真实值。
profvis(expr = {
for(i in 2:gibbs_iter){
# sample from posterior of phi
gibbs_res[i,p+1] <- rcond_post_phi(gibbs_res[i-1,1:p],
alpha, gamma, lambda, p)
# sample from posterior of beta vector ( using 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‘)
# calculate posterior means and credible intervals
post_burn_trim<-gibbs_res[seq(1000,gibbs_iter,100),]
colMeans(post_burn_trim)
apply(post_burn_trim, 2, quantile, p=c(.025,.975))
?
有一些改进的空间：
关于贝叶斯范式的好处是，所有推断都是使用后验分布完成的。现在，系数估计值是对数刻度，但是如果我们需要比值比，则只需对后验取幂。如果我们想要对比值比进行区间估计，那么我们就可以获取指数后验平局的2.5％和97.5％。
下面是使用R分析，显示了这一点。for循环运行Gibbs迭代。在每个Gibbs迭代中，我都调用函数rcond_post_beta_mh（），该函数使用MH从参数向量的条件后验中得出图形。
?
深入研究rcond_post_beta_mh（），我们看到子例程log_cond_post_beta（）是MH运行中的瓶颈。此函数是beta载体的对数条件后验密度，将其评估两次。
?
使用R语言进行Metroplis-in-Gibbs采样和MCMC运行分析
标签：基本算法 process 简单 有趣 title 慢慢 tin name 频繁
原文地址：https://www.cnblogs.com/tecdat/p/12731381.html
