0
点赞
收藏
分享

微信扫一扫

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析

 

对于许多模型,如物流模型,没有共轭先验 - 所以Gibbs不适用。正如我们在第一篇文章中看到的那样,蛮力网格方法太慢而无法扩展到真实环境。

这篇文章展示了我们如何使用Metropolis-Hastings(MH)从每个被阻挡的Gibbs迭代中的非共轭条件后验中进行采样 - 这是一种比网格方法更好的替代方案。

 模型

该示例的模拟数据是

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_03

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_05

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_07

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_09

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_11

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_13

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_15

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_17

  1. 我们有一个二元结果,我们使用非线性链接函数。
  2. 我们需要调整一个混淆器。
  3. 我们正在估算更多我们关心的参数。在这种情况下,我们真正关心治疗效果的估计​,因此其他系数在某种意义上是讨厌的参数。我不会说这是一个“高维”设置,但它肯定会让采样器紧张。

非标准化的后验条件 

让我们来看看这个模型的(非标准化的)条件后验。我不会进行推导,但它遵循我以前的帖子中使用的相同程序。

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_19

回想一下,我们正在建模

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_21

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_23

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_25

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_27

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_29

这里没有共轭!这种条件分布不是已知的分布,因此我们不能简单地使用Gibbs从中进行采样。相反,在每个gibbs迭代中,我们需要另一个采样步骤来从这个条件后验中抽取。第二个采样器将是MH采样器。如果您需要更新Gibbs,请参阅上面链接的上一篇文章。

旁注:条件后验

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_31

MH采样

目标是从中抽样

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_33

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_35

  1. 制定一个初步的“提案”,
  2. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_37

  3. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_38以便开始采样。
  4. 从围绕中心的分布画出
  5. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_39

  6. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_40。这称为提案分布,在标准情况下,必须是对称的  
  7. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_41

  8. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_42。我们可以用一个
  9. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_43

  10. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_44。现在让我们假设我们将提案分布的方差设置为某个常数。
  11. 现在我们从提案分发中抽取新的数据
  12. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_45

  13. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_46。然后,我们计算在上一次抽奖中评估的非标准化密度
  14. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_47

  15. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_48与当前建议的比率
  16. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_49

  17. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_50:   
  18. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_51

  19. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_52
  20. 如果该比率大于1,则当前建议的密度高于先前值的密度。所以我们“接受”提案并设定
  21. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_53

  22. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_54。然后,我们使用以中心为中心的提案分发重复步骤2-4  
  23. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_55

  24. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_56,然后生成新提案。如果该比率小于1,则当前建议值的密度低于先前的提议。在这种情况下,我们接受
  25. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_57

  26. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_58概率
  27. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_59

  28. 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_60

因此,总是接受产生更高条件后验评估的提议。然而,有时只接受密度评估较低的建议 - 提案的相对密度评估越低,其接受概率越低(直观!)。经过多次迭代,接受来自后部高密度区域的抽取,并且接受的建议序列“爬上”到高密度区域。一旦序列到达这个高密度区域,它往往会留在那里。因此,在许多方面,您可以将MH视为随机“爬山”算法。 

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_61

 

  • 跳跃分布的方差是一个重要参数。如果方差太小,当前的提议可能会非常接近最后一个值,因此
  • 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_63

  • 拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_64将接近1.因此我们会非常频繁地接受,但因为接受的值彼此如此接近,我们会爬到高位在许多迭代中缓慢地密度区域。如果方差太高,一旦到达那里,序列可能无法保留在高密度区域中。文献表明,在高维度(超过5个参数)中,最佳接受率约为24%。
  • 许多“自适应”MH方法是这里描述的基本算法的变体,但包括调谐周期以找到产生最佳接受率的跳跃分布方差。
  • MH中计算密集度最高的部分是密度评估。对于每个Gibbs迭代,我们必须两次评估4维密度:一次在当前提案中,一次在先前接受的提议中。
  • 虽然符号很容易扩展到高维度,但性能本身在更高的维度上会恶化。原因很简单,但非常有趣。​​Michael Betancourt的​​这篇​​论文​​  解释了Gibbs和MH在更高维度上的不足,并概述了哈密尔顿蒙特卡罗(HMC)如何克服这些困难。据我所知:在更高的尺寸,密度不等于体积。由于进入高容量区域确实是我们想要的,并且由于标准MH到达高密度区域,因此MH在高维度上无法有效地探索高容量区域。在低尺寸,密度很接近体积,所以这不是一个问题。
结果

 

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_方差_65

还有一些改进的余地:

  • 接受率仅为18%,我可以调整跳跃分布协方差矩阵以获得更优的速率。
  • 我认为更多的迭代肯定会有所帮助。这些链看起来还不错,但仍然相当自相关。

关于贝叶斯范式的好处是所有推理都是使用后验分布完成的。系数估计现在是对数尺度,但如果我们想要优势比,我们只是对后抽取进行取幂。如果我们想要比值比的区间估计,那么我们可以抓住指数后抽取的2.5和97.5百分位数。就这么简单 - 没有delta-method垃圾。

之前我曾提到MH是昂贵的,因为每次Gibbs迭代必须对log后验进行两次评估。下面是使用R包profviz进行的配置文件分析。for循环运行Gibbs迭代。在每个Gibbs迭代中,我调用函数rcond_post_beta_mh(),它使用MH从参数向量的条件后验产生绘制。请注意,它占用了大部分运行时。

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_迭代_67

潜入rcond_post_beta_mh(),我们看到子程序log_cond_post_beta()是MH运行的瓶颈。该函数是β向量的对数条件后验密度,其被评估两次。

拓端tecdat|R语言代写辅导使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析_采样器_69

 












举报

相关推荐

0 条评论