R语言使用Profviz进行Metropolis-in-Gibbs抽样和运行时间分析

系统管理员 2021-12-15 02:29 302阅读 0赞

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

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

**模型**

该示例的模拟数据是![N = 1000][N _ 1000]![点击并拖拽以移动][Image 1]患者的横截面数据集。有一个二元结果，![ÿ][latex.php_latex_Y_bg_ffffff_fg_000000_s_0]![点击并拖拽以移动][Image 1]一个二元治疗变量![一个][latex.php_latex_A_bg_ffffff_fg_000000_s_0]![点击并拖拽以移动][Image 1]，和一个混淆年龄。年龄是一个分类变量，有3个级别。我使用2个假人控制它（一个类别作为参考）。我用贝叶斯逻辑回归建模：

![logit（Y）= \\ beta\_0 + \\ beta\_1A + \\ beta\_2age\_1 + \\ beta\_3age\_2 = X \\ vec \\ beta][logit_Y_ _ beta_0 _ _ beta_1A _ _ beta_2age_1 _ _ beta_3age_2 _ X _ vec _ beta]![点击并拖拽以移动][Image 1]

![\\ beta\_0，\\ beta\_1，\\ beta\_2，\\ beta\_3 \\ sim N（\\ lambda，\\ phi）][beta_0_ beta_1_ beta_2_ beta_3 _ sim N_ lambda_ phi]![点击并拖拽以移动][Image 1]

以上，![\\拉姆达][latex.php_latex_5Clambda_bg_ffffff_fg_000000_s_0]![点击并拖拽以移动][Image 1]假设已知。注意我![X][]![点击并拖拽以移动][Image 1]用来表示1000×4模型矩阵并![\\ vec \\ beta][vec _ beta]![点击并拖拽以移动][Image 1]表示4×1参数向量。我还在![\\披][latex.php_latex_5Cphi_bg_ffffff_fg_000000_s_0]![点击并拖拽以移动][Image 1]已知的超参数之前放置了反伽马。对于Metropolis-in-Gibbs来说，这是一个相当现实的激励示例：

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

# 非标准化的后验条件  #

让我们来看看这个模型的（非标准化的）条件后验。我不会进行推导，但它遵循我以前的帖子中使用的相同程序。  
![log f（\\ vec \\ beta | \\ phi，X，Y）\\ propto \\ sum\_ \{i = 1\} ^ n y\_i \\ cdot log（p\_i）+（1  -  y\_i）\\ cdot log（1  -  p\_i） -  \\ frac \{1\} \{2 \\ phi\}（\\ vec \\ beta  -  \\ vec \\ lambda）'（\\ vec \\ beta  -  \\ vec \\ lambda）][log f_ vec _ beta _ _ phi_X_Y_ propto _ sum_ _i _ 1_ _ n y_i _ cdot log_p_i_1 - y_i_ cdot log_1 - p_i_ - _ frac _1_ _2 _ phi_ vec _ beta - _ vec _ lambda_ vec _ beta - _ vec _ lambda]![点击并拖拽以移动][Image 1]  
回想一下，我们正在建模![p\_i = expit（x\_i'\\ vec \\ beta）][p_i _ expit_x_i_ vec _ beta]![点击并拖拽以移动][Image 1]。这里![X\_I][X_I]![点击并拖拽以移动][Image 1]表示![I ^ \{\}第i个][I _ _i]![点击并拖拽以移动][Image 1]模型矩阵的行![X][X 1]![点击并拖拽以移动][Image 1]，并![一世][latex.php_latex_i_bg_ffffff_fg_000000_s_0]![点击并拖拽以移动][Image 1]指示患者。

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

旁注：条件后验![\\披][latex.php_latex_5Cphi_bg_ffffff_fg_000000_s_0]![点击并拖拽以移动][Image 1] *是*共轭的。因此，在每个Gibbs迭代中，我们可以使用标准函数从反伽马中进行采样。这里不需要第二个采样步骤。唷！

# MH采样 #

目标是从中抽样![log f（\\ vec \\ beta | \\ phi，X，Y）][log f_ vec _ beta _ _ phi_X_Y]![点击并拖拽以移动][Image 1]。请注意，这是一个四维密度。为了简化说明，假设我们只有一个![\\公测][latex.php_latex_5Cbeta_bg_ffffff_fg_000000_s_0]![点击并拖拽以移动][Image 1]，并且它的条件密度是单峰的。MH采样器的工作原理如下：

1.  制定一个初步的“提案”，![\\的β^ \{（0）\}][_0]![点击并拖拽以移动][Image 1]以便开始采样。
2.  从围绕中心的分布画出![\\的β^ \{（0）\}][_0]![点击并拖拽以移动][Image 1]。这称为提案分布，在标准情况下，必须是对称的  ![\\的β^ \{（0）\}][_0]![点击并拖拽以移动][Image 1]。我们可以用一个![N（\\ beta ^ \{（0）\}，\\ sigma ^ 2）][N_ beta _ _0_ sigma _ 2]![点击并拖拽以移动][Image 1]。现在让我们假设我们将提案分布的方差设置为某个常数。
3.  现在我们从提案分发中抽取新的数据![\\测试^ \*][latex.php_latex_5Cbeta_5E_2A_bg_ffffff_fg_000000_s_0]![点击并拖拽以移动][Image 1]。然后，我们计算在上一次抽奖中评估的非标准化密度![\\的β^ \{（0）\}][_0]![点击并拖拽以移动][Image 1]与当前建议的比率![\\测试^ \*][latex.php_latex_5Cbeta_5E_2A_bg_ffffff_fg_000000_s_0]![点击并拖拽以移动][Image 1]：   ![r = \\ frac \{f（\\ beta ^ \*）\} \{f（\\ beta ^ \{（0）\}）\}][r _ _ frac _f_ beta _ _ _f_ beta _ _0]![点击并拖拽以移动][Image 1]
4.  如果该比率大于1，则当前建议的密度高于先前值的密度。所以我们“接受”提案并设定![\\ beta ^ \{（1）\} = \\ beta ^ \*][beta _ _1_ _ _ beta _]![点击并拖拽以移动][Image 1]。然后，我们使用以中心为中心的提案分发重复步骤2-4  ![\\的β^ \{（1）\}][_1]![点击并拖拽以移动][Image 1]，然后生成新提案。如果该比率小于1，则当前建议值的密度低于先前的提议。在这种情况下，我们接受![\\测试^ \*][latex.php_latex_5Cbeta_5E_2A_bg_ffffff_fg_000000_s_0 1]![点击并拖拽以移动][Image 1]概率![\[R][R]![点击并拖拽以移动][Image 1]。

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

![\\西格马^ 2][2]![点击并拖拽以移动][Image 1]我们有一个协方差矩阵，而不是标量方差参数。因此，我们的提议是系数的向量。在这个意义上，我们运行一个被阻挡的Gibbs - 使用MH在每次迭代中绘制整个系数块。

*  跳跃分布的方差是一个重要参数。如果方差太小，当前的提议可能会非常接近最后一个值，因此![\[R][R]![点击并拖拽以移动][Image 1]将接近1.因此我们会非常频繁地接受，但因为接受的值彼此如此接近，我们会爬到高位在许多迭代中缓慢地密度区域。如果方差太高，一旦到达那里，序列可能无法保留在高密度区域中。文献表明，在高维度（超过5个参数）中，最佳接受率约为24％。
 *  许多“自适应”MH方法是这里描述的基本算法的变体，但包括调谐周期以找到产生最佳接受率的跳跃分布方差。
 *  MH中计算密集度最高的部分是密度评估。对于每个Gibbs迭代，我们必须两次评估4维密度：一次在当前提案中，一次在先前接受的提议中。
 *  虽然符号很容易扩展到高维度，但性能本身在更高的维度上会恶化。原因很简单，但非常有趣。[Michael Betancourt的][Michael Betancourt]这篇[论文][Michael Betancourt]  解释了Gibbs和MH在更高维度上的不足，并概述了哈密尔顿蒙特卡罗（HMC）如何克服这些困难。据我所知：在更高的尺寸，密度不等于体积。由于进入高容量区域确实是我们想要的，并且由于标准MH到达高*密度*区域，因此MH在高维度上无法有效地探索高容量区域。在低尺寸，密度很接近体积，所以这不是一个问题。

# 结果 #

![LogitMCMCChains.png][]![点击并拖拽以移动][Image 1]

还有一些改进的余地：

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

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

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

![屏幕截图2017-11-06 at 9.37.11 PM][2017-11-06 at 9.37.11 PM]![点击并拖拽以移动][Image 1]

潜入rcond\_post\_beta\_mh（），我们看到子程序log\_cond\_post\_beta（）是MH运行的瓶颈。该函数是β向量的对数条件后验密度，其被评估两次。

![屏幕截图2017-11-06 at 9.41.55 PM.png][2017-11-06 at 9.41.55 PM.png]![点击并拖拽以移动][Image 1]

## 还有问题吗？联系我们！ ##

[大数据部落][Link 1] \-中国专业的第三方数据服务提供商，提供定制化的一站式数据挖掘和统计分析咨询服务

统计分析和数据挖掘咨询服务：[y0.cn/teradat][Link 1]（咨询服务请联系[官网客服][Link 2]）

[![点击这里给我发消息][20181218152712945]][20181218152712945 1]![点击并拖拽以移动][Image 1]QQ：3025393450

![2019031315051532.png][]![点击并拖拽以移动][Image 1]

**【服务场景】  **

科研项目; 公司项目外包;线上线下一对一培训;数据采集;学术研究;报告撰写;市场调查。

**【大数据部落】**提供定制化的一站式数据挖掘和统计分析咨询服务

[![【大数据部落】大数据部落提供定制化的一站式数据挖掘和统计分析咨询服务][20181218152712961]][Link 1]![点击并拖拽以移动][Image 1]

转载于:https://www.cnblogs.com/tecdat/p/10751633.html

[N _ 1000]: https://s0.wp.com/latex.php?latex=N%3D1000&bg=ffffff&fg=000000&s=0
[Image 1]: 
[latex.php_latex_Y_bg_ffffff_fg_000000_s_0]: https://s0.wp.com/latex.php?latex=Y+&bg=ffffff&fg=000000&s=0
[latex.php_latex_A_bg_ffffff_fg_000000_s_0]: https://s0.wp.com/latex.php?latex=A&bg=ffffff&fg=000000&s=0
[logit_Y_ _ beta_0 _ _ beta_1A _ _ beta_2age_1 _ _ beta_3age_2 _ X _ vec _ beta]: https://s0.wp.com/latex.php?latex=logit%28Y%29+%3D+%5Cbeta_0+%2B+%5Cbeta_1A+%2B+%5Cbeta_2age_1+%2B+%5Cbeta_3age_2+%3D+X%5Cvec+%5Cbeta&bg=ffffff&fg=000000&s=0
[beta_0_ beta_1_ beta_2_ beta_3 _ sim N_ lambda_ phi]: https://s0.wp.com/latex.php?latex=%5Cbeta_0%2C+%5Cbeta_1%2C+%5Cbeta_2%2C+%5Cbeta_3+%5Csim+N%28%5Clambda%2C+%5Cphi%29+&bg=ffffff&fg=000000&s=0
[latex.php_latex_5Clambda_bg_ffffff_fg_000000_s_0]: https://s0.wp.com/latex.php?latex=%5Clambda&bg=ffffff&fg=000000&s=0
[X]: https://s0.wp.com/latex.php?latex=X+&bg=ffffff&fg=000000&s=0
[vec _ beta]: https://s0.wp.com/latex.php?latex=%5Cvec+%5Cbeta&bg=ffffff&fg=000000&s=0
[latex.php_latex_5Cphi_bg_ffffff_fg_000000_s_0]: https://s0.wp.com/latex.php?latex=%5Cphi&bg=ffffff&fg=000000&s=0
[beta_1]: https://s0.wp.com/latex.php?latex=%5Cbeta_1&bg=ffffff&fg=000000&s=0
[log f_ vec _ beta _ _ phi_X_Y_ propto _ sum_ _i _ 1_ _ n y_i _ cdot log_p_i_1 - y_i_ cdot log_1 - p_i_ - _ frac _1_ _2 _ phi_ vec _ beta - _ vec _ lambda_ vec _ beta - _ vec _ lambda]: https://s0.wp.com/latex.php?latex=log+f%28+%5Cvec+%5Cbeta+%7C+%5Cphi%2C+X%2C+Y%29+%5Cpropto+%5Csum_%7Bi%3D1%7D%5En+y_i%5Ccdot+log%28p_i%29+%2B+%281+-+y_i%29+%5Ccdot+log%281+-+p_i%29+-+%5Cfrac%7B1%7D%7B2%5Cphi%7D%28%5Cvec+%5Cbeta+-+%5Cvec+%5Clambda%29%27%28%5Cvec+%5Cbeta+-+%5Cvec+%5Clambda%29+&bg=ffffff&fg=000000&s=0
[p_i _ expit_x_i_ vec _ beta]: https://s0.wp.com/latex.php?latex=p_i+%3D+expit%28x_i%27+%5Cvec+%5Cbeta%29&bg=ffffff&fg=000000&s=0
[X_I]: https://s0.wp.com/latex.php?latex=x_i&bg=ffffff&fg=000000&s=0
[I _ _i]: https://s0.wp.com/latex.php?latex=i%5E%7Bith%7D&bg=ffffff&fg=000000&s=0
[X 1]: https://s0.wp.com/latex.php?latex=X&bg=ffffff&fg=000000&s=0
[latex.php_latex_i_bg_ffffff_fg_000000_s_0]: https://s0.wp.com/latex.php?latex=i&bg=ffffff&fg=000000&s=0
[log f_ vec _ beta _ _ phi_X_Y]: https://s0.wp.com/latex.php?latex=log+f%28%5Cvec+%5Cbeta+%7C+%5Cphi%2C+X%2C+Y%29&bg=ffffff&fg=000000&s=0
[latex.php_latex_5Cbeta_bg_ffffff_fg_000000_s_0]: https://s0.wp.com/latex.php?latex=%5Cbeta&bg=ffffff&fg=000000&s=0
[_0]: https://s0.wp.com/latex.php?latex=%5Cbeta%5E%7B%280%29%7D&bg=ffffff&fg=000000&s=0
[N_ beta _ _0_ sigma _ 2]: https://s0.wp.com/latex.php?latex=N%28%5Cbeta%5E%7B%280%29%7D%2C+%5Csigma%5E2%29+&bg=ffffff&fg=000000&s=0
[latex.php_latex_5Cbeta_5E_2A_bg_ffffff_fg_000000_s_0]: https://s0.wp.com/latex.php?latex=%5Cbeta%5E%2A+&bg=ffffff&fg=000000&s=0
[r _ _ frac _f_ beta _ _ _f_ beta _ _0]: https://s0.wp.com/latex.php?latex=r+%3D+%5Cfrac%7Bf%28%5Cbeta%5E%2A%29%7D%7Bf%28%5Cbeta%5E%7B%280%29%7D%29%7D&bg=ffffff&fg=000000&s=0
[beta _ _1_ _ _ beta _]: https://s0.wp.com/latex.php?latex=%5Cbeta%5E%7B%281%29%7D+%3D+%5Cbeta%5E%2A&bg=ffffff&fg=000000&s=0
[_1]: https://s0.wp.com/latex.php?latex=%5Cbeta%5E%7B%281%29%7D+&bg=ffffff&fg=000000&s=0
[latex.php_latex_5Cbeta_5E_2A_bg_ffffff_fg_000000_s_0 1]: https://s0.wp.com/latex.php?latex=%5Cbeta%5E%2A&bg=ffffff&fg=000000&s=0
[R]: https://s0.wp.com/latex.php?latex=r&bg=ffffff&fg=000000&s=0
[2]: https://s0.wp.com/latex.php?latex=%5Csigma%5E2&bg=ffffff&fg=000000&s=0
[Michael Betancourt]: https://arxiv.org/pdf/1701.02434.pdf
[LogitMCMCChains.png]: https://stablemarkets.files.wordpress.com/2017/11/logitmcmcchains.png?w=1100
[2017-11-06 at 9.37.11 PM]: https://stablemarkets.files.wordpress.com/2017/11/screen-shot-2017-11-06-at-9-37-11-pm.png?w=1100
[2017-11-06 at 9.41.55 PM.png]: https://stablemarkets.files.wordpress.com/2017/11/screen-shot-2017-11-06-at-9-41-55-pm.png?w=1100
[Link 1]: http://y0.cn/teradat
[Link 2]: http://dwz.cn/YntAj
[20181218152712945]: https://img-blog.csdnimg.cn/20181218152712945
[20181218152712945 1]: http://wpa.qq.com/msgrd?v=3&uin=3025393450&site=qq&menu=yes
[2019031315051532.png]: https://img-blog.csdnimg.cn/2019031315051532.png
[20181218152712961]: https://img-blog.csdnimg.cn/20181218152712961