R语言BUGS/JAGS贝叶斯分析：马尔科夫链蒙特卡洛方法（MCMC）采样

全文链接：http://tecdat.cn/?p=17884

相关视频：马尔可夫链原理可视化解释与R语言区制转换Markov regime switching实例

马尔可夫链原理可视化解释与R语言区制转换Markov regime switching实例

相关视频

马尔可夫链蒙特卡罗方法MCMC原理与R语言实现

，时长08:47

马尔科夫链蒙特卡洛方法

MCMC的关键如下：

跳跃概率的比例与后验概率的比例成正比。

跳跃概率可以表征为：

概率（跳跃）*概率（接受）

从长远来看，该链将花费大量时间在参数空间的高概率部分，从而实质上捕获了后验分布。有了足够的跳跃，长期分布将与联合后验概率分布匹配。

MCMC本质上是一种特殊类型的随机数生成器，旨在从难以描述（例如，多元，分层）的概率分布中采样。在许多/大多数情况下，后验分布是很难描述的概率分布。

MCMC使您可以从实际上不可能完全定义的概率分布中进行采样！

令人惊讶的是，MCMC的核心并不难于描述或实施。让我们看一个简单的MCMC算法。

Metropolis-Hastings算法

该算法与模拟退火算法非常相似。

MH算法可以表示为：

Prob（acceptB | A）= min（1，Posterior（B）Posterior（A）⋅Prob（ba）Prob（ab））

请注意，从本质上讲，这与“ Metropolis”模拟退火算法相同，后验概率代替了概率，并且k参数设置为1。

二元正态例子

请记住，MCMC采样器只是随机数生成器的一种。我们可以使用Metropolis-Hastings采样器来开发自己的随机数生成器，生成进行简单的已知分布。在此示例中，我们使用MH采样器从标准双变量正态概率分布生成随机数。

对于这个简单的示例，我们不需要MCMC采样器。一种实现方法是使用以下代码，该代码从具有相关参数ρ的双变量标准正态分布中绘制并可视化任意数量的独立样本。

然后，我们可以使用MH采样器从该已知分布中获取随机样本…

让我们尝试解决一个问题。

MCMC对粘液瘤病进行调查

选择使用Gamma分布。这是经验分布：

我们需要估算最适合此经验分布的伽马速率和形状参数。这是适合此分布的Gamma的一个示例：

二维（对数）似然面：

这是MH算法的实现，用于找到后验分布！

首先，我们需要一个似然函数–这次，我们将返回真实概率–而不是对数转换的概率

然后，我们需要预先分配参数！在这种情况下，我们分配gamma（shape = 0.01，scale = 100）和gamma（shape = 0.1，scale = 10）的分布（均值为1且方差略低）：

我们假设形状和比例 在先验中是独立的（联合先验的乘法概率）。然而，并没有对后验参数相关性提出相同的假设，因为概率可以反映在后验分布中。

然后，我们需要一个函数，该函数可以计算参数空间中任何给定跳转的后验概率比率。因为我们正在处理 后验概率的比率，所以我们不需要计算归一化常数。

无需归一化常数，我们只需要计算加权似然比（即先验加权的似然比）

然后，我们需要一个函数进行新的推测或在参数空间中跳转：

现在，我们准备实现Metropolis-Hastings MCMC算法：

我们需要一个初始点：

测试函数

现在让我们看一下Metropolis-Hastings：

我们运行更长的时间

更长的时间

看起来更好！搜索算法可以很好地找到参数空间的高似然部分！

现在，让我们看一下“ shape”参数的链

对于比例参数

我们可以说这些链已经收敛于形状参数的后验分布吗？

首先，链的起点“记住”起始值，因此不是固定分布。我们需要删除链的第一部分。

但这看起来还不是很好。让我们运行更长的时间，看看是否得到的东西看起来更像是随机数生成器（白噪声）

让我们首先删除前5000个样本作为预烧期

现在，让我们再次看一下链条

在评估这些迹线图时，我们希望看到看起来像白噪声的“平稳分布”。该轨迹图看起来可能具有一些自相关。解决此问题的一种方法是稀疏MCMC样本：

现在我们可以检查我们的后验分布！

我们可以像以前一样可视化。

可以修改Metropolis-Hastings MCMC方法来拟合任意模型的任意数量的自由参数。但是，MH算法本身不一定是最有效和灵活的。在实验中，我们使用吉布斯采样，大多采用建模语言 BUGS 。

注意：BUGS实现（例如JAGS）实际上倾向于结合使用MH和Gibbs采样，MH和Gibbs采样器并不是唯一的MCMC例程。例如，“ stan”使用MH采样的一种改进形式，称为“ Hamiltonian Monte Carlo”。

吉布斯Gibbs采样器

Gibbs采样器非常简单有效。基本上，该算法从完整的条件概率分布（即， 在模型中所有其他参数的已知值作为条件的条件下，对任意参数i的后验分布）中进行 连续采样 。

在很多情况下，我们不能直接制定出我们的模型后验分布，但我们 可以 分析出条件后验分布。尽管如此，即使它在分析上不易处理，我们也可以使用单变量MH程序作为最后方法。

问：为什么Gibbs采样器通常比纯MH采样器效率更高？

二元正态例子

MCMC采样器只是随机数生成器的一种。我们可以使用Gibbs采样器来开发自己的随机数生成器，以实现相当简单的已知分布。在此示例中，我们使用Gibbs采样器从标准双变量正态概率分布生成随机数。注意，吉布斯采样器在许多方面都比MH算法更简单明了。

然后，我们可以使用Gibbs采样器从该已知分布中获取随机样本…

在这里，马尔可夫链的样本中有很多明显的自相关。Gibbs采样器经常有此问题。

示例

BUGS语言

最后，让我们为我们最喜欢的粘瘤病示例创建一个Gibbs采样器，为此，我们将使用BUGS语言（在JAGS中实现）来帮助我们！

JAGS，全称是Just another Gibbs sampler，是基于BUGS语言开发的利用MCMC来进行贝叶斯建模的软件包。它没有提供建模所用的GUI以及MCMC抽样的后处理，这些要在其它的程序软件上来处理，比如说利用R包（rjags）来调用JAGS并后处理MCMC的输出。JAGS相对于WinBUGS/OpenBUGS的主要优点在于平台的独立性，可以应用于各种操作系统，而WinBUGS/OpenBUGS只能应用于windows系统；JAGS也可以在64-bit平台上以64-bit应用来进行编译。

BUGS语言看起来与R类似，但是有几个主要区别：

首先，BUGS是一种编译语言，因此代码中的操作顺序并不重要

BUGS不是矢量化的-您需要使用FOR循环

在BUGS中，几个概率分布的参数差异很大。值得注意的是，正态分布通过均值和精度（1 / Variance ）进行参数化。

这是用BUGS语言实现的粘液病示例：

我们可以使用R中的“ cat”函数将此模型写到您的工作目录中的文本文件中：

现在我们已经将BUGS模型打包为文本文件，我们将数据捆绑到一个列表对象中，该列表对象包含BUGS代码中引用的所有相关数据：

然后，我们需要为所有参数定义初始值。将其定义为一个函数很方便，因此可以使用不同的起始值来初始化每个MCMC链。

现在我们可以调用JAGS！

评估收敛

第一步是视觉检查-我们寻找以下内容来评估收敛性：

当视为“轨迹图”时，每个参数的链应看起来像白噪声（模糊的毛毛虫）或类似的噪声

多个具有不同起始条件的链条看起来应该相同

我们可能在这里可以做得更好的一种方法是使链条运行更长的时间，并丢弃初始样本我们还可以。

通过细化链来尝试减少自相关，我们每20个样本中仅保留1个。

从视觉上看，这看起来更好。现在我们可以使用更多的定量收敛指标。

Gelman-Rubin诊断

一种简单而直观的收敛诊断程序是Gelman-Rubin诊断程序，该诊断程序基于简单的蒙特卡洛误差来评估链之间是否比应有的链距更大：

通常，1.1或更高的值被认为收敛不佳。为模型中的所有可用参数计算GR诊断。如果测试失败，则应尝试运行更长的链！

所以这个模型看起来不错！