metropolis算法 Stata软件贝叶斯统计应用：MCMC和Metropolis

本文由中国科学软件网翻译整理
在本文中，我从非技术的角度介绍马尔可夫链蒙特卡罗，通常简称MCMC。MCMC广泛应用于贝叶斯统计模型拟合。MCMC有很多不同的方法，这里我主要介绍Metropolis-Hastings算法。为了简单起见，我将忽略一些细节，我强烈建议您在练习使用MCMC之前阅读[BAYES]手册。
我们继续上一篇文章贝叶斯统计，介绍第1部分:基本概念中提到的投掷硬币的例子。我们感兴趣的是参数θ的后验分布，即抛硬币时“头”在上面的概率。我们的先验分布是平坦的，没有β分布参数1和1的信息。然后我们会用二项式似然函数来量化我们的实验数据，以及抛硬币10次和“抬头”4次的结果。我们可以使用MCMC的M-H算法来生成后验分布θ的样本。这个样本然后可以用来估计问题，如后验分布的平均值。
这项技术有三个基本部分:
1.蒙特卡洛
2.马尔可夫链
3.m-H算法
蒙特卡罗方法
蒙特卡洛是指产生伪随机数(简称伪随机数)的方法。图1说明了蒙特卡罗实验的一些特征。
图1:提案分布、跟踪图和密度图

我们不知道后验分布的函数形式，但是在这个例子中，我们可以用先验分布和似然函数的乘积来代替。比率的计算并不那么简单，但我们可以让事情变得更容易。
图4的步骤1示出了(θnew=0.380)和(θ t1 = 0.0.517)的后验概率比等于1.307。在步骤2中，我们计算接受概率α (θ new，θ t1)，它是r(θnew，θt -1)和1的最低比率。步骤2是必要的，因为比率必须介于0和1之间。
接受概率等于1，所以我们接受(θnew=0.380)，并建议在下一次迭代中，分布的平均值变为0.380。
图5:带M-H算法的MCMC:迭代2

文章图片

图5示出了根据平均值为(= 0.380)的后验分布绘制下一次迭代(= 0.286)。步骤1
计算出的比值r (θ new，θ t1)等于0.747且小于1。在步骤2中计算的接受概率α(θnew，θt1)等于0.747。
我们不会因为接受概率小于1而自动拒绝θnew。相反，我们可以从步骤3中画出均匀分布的随机数u(0，1)。如果u小于接受概率，我们接受建议值θnew，否则拒绝θnew，保留θ t1。
这里u=0.094小于接受概率(0.747)，所以我们接受(θnew=0.286)。在图5中，我用绿色的θnew和0.286表示它们已被接受。
图6:带M-H算法的MCMC:迭代3

文章图片

图6示出了根据平均值为(= 0.286)的后验分布绘制下一次迭代(= 0.088)。步骤1中计算的比值r (θ new，θ t1)等于0.039，小于1。步骤2中计算的接受概率α (θ new，θ t1)等于0.039。步骤3中计算的u值为0.247，大于接受概率。因此，在步骤4中，我们拒绝θnew=0.088，并保持θ t1 = 0.286。
让我们用M-H算法的MCMC画出10000个θ的随机值，看看它是如何工作的。
动画3:具有M-H算法的MCMC

文章图片

动画图3显示了几个问题。首先，建议分布会随着大多数迭代而改变。应该注意的是，如果θ新值被接受，它们将以绿色表示，如果被拒绝，它们将以红色表示。其次，当我们只使用MCMC时，轨迹图不显示随机游走模式，所有迭代中的变化都是相似的。最后，密度图看起来像一个有用的分布。
旋转密度图结果，仔细观察。
图7:一个来自后验分布的样本，它是用多模型预测控制算法生成的