再谈贝叶斯估计, 从MCMC和MH算法说起
可有偿投稿计量经济圈，计量相关则可
邮箱：econometrics666@sina.cn
所有计量经济圈方法论丛的do文件, 微观数据库和各种软件都放在社群里.欢迎到计量经济圈社群交流访问.这篇文章来自计量社群, 来源于网络但没找到源头.
推荐阅读：
1.实证研究中用到的135篇文章, 社科学者常用toolkit
2.1998-2016年中国地级市年均PM2.5数据release
今天，计量经济圈将为圈友引荐一篇文章，主要帮助圈友理解贝叶斯估计中运用到的算法。贝叶斯估计在可能当下的争议还比较大，因为先验信息的判断较为主观，这给了频率学派很大的舆论阵地去为主流估计理论辩护。但，如果我们准备把p值埋葬的话(美国宣布禁用p值，原来p值很危险，如何取代p值？; 科学家倡议P值需要0.005，显著性判断才成立)，那贝叶斯给出的判断准则可能就是模型的新评估标准，所以还是先对马尔可夫蒙特卡洛等概念有一些了解比较好。
在这篇文章中，我从非技术层面来介绍一下Markov chain Monte Carlo，通常简称为MCMC。MCMC被广泛应用于贝叶斯统计模型拟合。MCMC有很多不同的方法，这里我将主要介绍Metropolis–Hastings 算法。
我们对参数θ的后验分布感兴趣，也就是投掷硬币时，“头”在上面的概率。我们的先验分布是扁平的,无信息Beta分布参数1和1。然后我们将使用二项式似然函数对我们的实验数据进行量化，10次投掷硬币4次“头”朝上的结果。我们可以使用MCMC的M-H算法来生成后验分布θ的样本。然后可以使用这个样本来估算如后验分布的均值等问题。
这种技术有三个基本部分：
1． Monte Carlo
2． Markov chains
3. M–H algorithm
Monte Carlo methods
蒙特卡罗是指依靠伪随机数的生成方法（简称伪随机数）。图1说明了蒙特卡罗实验的一些特点。
Figure 1: Proposal distributions, trace plots, and density plots
在这个例子中，我将从正态分布的均值0.5和任意方差σ中生成一系列随机数。正态分布称为建议分布。
图中右上角部分称为轨迹图，它会根据绘图的顺序来显示θ值。它也显示建议分布顺时针旋转90度的情况，我每画一个θ值就会将它向右平移。绿点轨迹图显示了θ的当前值。
左上角的密度图是逆时针旋转90度的直方图样本。我旋转了坐标轴，因此θ轴线是平行的。同样的，绿色的点密度图显示了θ的当前值。
让我们绘制10,000个θ 随机值来看看是如何工作的。
Animation 1: A Monte Carlo experiment
动画图1说明了几个重要的特点。首先建议分布不会从一个迭代中更改到另一个中。其次，随着样本容量的增加，密度图看起来越来越像建议分布。第三，轨迹图是平稳的--所有迭代的变量，变化形式看起来都是一样的。
蒙特卡罗仿真生成的样本看起来像建议分布，很多时候也正是我们所需的。但是我们需要一个额外的工具从后验分布中生成一个样本。
Markov chain Monte Carlo methods
Markov chain是一个数字序列,每个数字都依赖于序列中的前一个数。例如，我们从一个正常的建议分布与平均等于θ的前一个值来绘制θ值。图2显示了 轨迹图和密度图，当前的θ值（θt=0.530）是从一个平均等于θ前值的建议分布 (θt−1=0.712)来绘制的。
Figure 2: A draw from a MCMC experiment
下面图3显示序列中的下一次迭代的轨迹图和密度图。建议密度的均值现在是θt−1=0.530，并且我们从新的建议分布来绘制出了随机数值θt=0.411。
Figure 3: The next iteration of the MCMC experiment
让我们使用MCMC绘制10,000个θ 随机值来看看是如何工作的。
Animation 2: A MCMC experiment
动画图2显示了Monte Carlo实验和MCMC实验直接的差异。首先，建议分布是随着每个迭代变化而变化的。这就生成一个随机游走模式的轨迹图：所有迭代的变化是不一样的。其次，产生的密度图看上去不像建议分布或其他任何有用的分布，肯定也不是一个后验分布。
MCMC从后验分布获取样本失败的原因是我们还没在生成样本的过程中输入过任何后验分布的信息。我们可以通过保持θ值来提高样本，θ值更可能是在后验分布下的值而不太可能是丢弃值。
但是基于后验分布，θ值接受和拒绝建议值最明显的困难是我们通常不知道后验分布的函数形式。如果我们知道它的函数形式，就可以直接计算概率，而不用生成一个随机样本。那么我们如何不用知道函数形式就可以接受或拒绝建议θ值呢? 答案就是M-H算法！
MCMC和M–H 算法
M-H算法可以用在我们不知道后验分布函数形式的情况下来判断是否接受θ建议值。让我们将M-H算法分成几步，然后经过几次迭代来看看它是怎么实现的。
Figure 4: MCMC with the M–H algorithm: Iteration 1
图4显示了等于(θt−1=0.0.517)的建议分布轨迹图和密度图。我们已经绘制了一个预估值θ(θnew=0.380)，我将这个值定义为θnew是因为还没有确定是否要用这个值。
我们从步骤1计算比率开始使用M-H算法
我们不知道后验分布的函数形式，但是在这个例子中，我们可以替代先验分布和似然函数的乘积。比率的计算并不是都那么简单，但是我们可以把事情变得简单些。
图4步骤1显示了(θnew=0.380) 和(θt−1=0.0.517) 的后验概率比等于1.307。步骤2我们计算接受概率α(θnew,θt−1)，也就是最低比值r(θnew, θt -1)和1。步骤2是必须的，因为比率肯定在0和1之间。
接受概率等于1，所以我们接受(θnew=0.380),建议分布的平均值在下个迭代中就变成了0.380.
Figure 5: MCMC with the M–H algorithm: Iteration 2
图5显示了下一个迭代(θnew=0.286)是根据平均值为 (θt−1=0.380)后验分布绘制的。步骤1
计算的比率r(θnew,θt−1)等于0.747，小于1。步骤2中计算的接受概率α(θnew,θt−1)等于0.747。
我们不会自动拒绝θnew，因为接受概率小于1。相反，我们可以绘制来自步骤3中的均匀分布（0，1）的随机数U。如果U小于接受概率，我们就接受建议值θnew，反之，我们拒绝θnew并保留 θt−1。
在这里u=0.094小于了接受概率(0.747)，所以我们就接受(θnew=0.286)。我在图5中 将θnew和0.286以绿色来表示它们已经被接受。
Figure 6: MCMC with the M–H algorithm: Iteration 3
图6显示了下一个迭代(θnew=0.088)是根据平均值为(θt−1=0.286)后验分布绘制的。步骤1计算的比率 r(θnew,θt−1)等于0.039，小于1. 步骤2中计算的接受概率α(θnew,θt−1)等于0.039.步骤3中计算的U 的值是0.247，大于接受概率。因此在步骤4中，我们拒绝θnew=0.088并保留θt−1=0.286 。
让我们使用M-H算法的MCMC绘制10,000个θ 随机值来看看是如何工作的。
Animation 3: MCMC with the M–H algorithm
动画图3表明了几个问题。首先，建议分布会随着大部分迭代变化而变化。需要注意的是, 如果θnew值被接受，它们会以绿色表示，如果被拒绝会以红色表示。其次，我们只使用MCMC的时候，轨迹图不显示随机游走模式，所有迭代中的变化是相似的。最后，密度图看起来像一个有用的分布。
旋转一下密度图结果，更仔细的看看。
Figure 7: A sample from the posterior distribution generated with MCMC with the M–H algorithm
图7显示了我们使用M-H算法MCMC生成的样本直方图。在这个例子中，我们知道后验分布的参数是Beta 5和7. 红色线条显示后验分布分析，蓝色线条是样本的核密度图。核密度图跟Beta（5,7）分布相似，这说明我们的样本是一个很好的近似的理论后验分布。
我们可以使用样本来计算后验分布均值或中位数，95%的置信区间的概率，或θ落在任意区间的概率。
Stata中的MCMC和M-H算法
让我们回到使用bayesmh的投掷硬币实例中。随着二项式的可能性，我们指定一个beta分布参数1和1。
Example 1: Using bayesmh with a Beta(1,1) prior
输出结果告诉我们bayesmh跑完了12,500个MCMC迭代。“Burn-in”告诉我们为了减少链中随机起始值的影响，迭代中的2500个被丢弃。下面那行告诉我们，最终的MCMC样本大小为10,000, 再下面一行表明我们的数据集有10个观测值。接受概率θ建议值比例包含在最终的MCMC样本中。我建议您阅读Stata Bayesian Analysis Reference Manual 并了解效率的定义，但是高效率表明低自相关，而低效率表明高自相关。 Monte Carlo standard error (MCSE)显示在系数表中，估算后验均值的近似误差。
检查chain的收敛性
“convergence”这个词在MCMC和极大似然中的表述有不同的含义。用于最大似然估计算法迭代直至收敛到最大。MCMC链不迭代直到确定最佳值。外链简单迭代直到达到要求的样本量大小，任何算法停止。事实上，外链停止运行并不表明后验分布中最佳样本的生成。我们必须检查样本以制止问题的出现。我们可以用bayesgraph diagnostics图形形式来检查样本。
图8显示包含MCMC样本的轨迹图、直方图和密度图的诊断图和correlegram。轨迹图有一个平稳模式，这也正是我们希望看到的。动画图2中的随机游走模式显示了chain存在的问题。直方图没有任何特殊特征如多模式。完整样本的核密度图，chain的前半部分，chain的后半部分都看起来很相似并没有任何特殊的特征如chain的前半部分和后半部分密度不同。使用Markov chain来生成样本并产生一个自相关，但是自相关会随着较大的滞后值迅速减小。这些图形并不能说明我们的样本有任何问题。
总结
这篇文章介绍了使用M-H算法MCMC背后的想法。请注意我省略了一些细节，忽略了一些假定条件以便让事情变得简单，顺着我们的感觉往下发展。Stata的bayesmh命令实际上执行了一个更难的算法，我们称之为带M-H算法的自适应MCMC。但是基本概念是一样的，希望能给大家一些启发。
所有计量经济圈方法论丛的do文件, 微观数据库和各种软件都放在社群里.欢迎到计量经济圈社群交流访问.
可以到计量经济圈社群进一步访问交流各种学术问题，这年头，我们不能强调一个人的英雄主义，需要多多汲取他人的经验教训来让自己少走弯路。
计量经济圈当前有几个阵地，他们分别是如下4个matrix：
①小鹅社群：数据软件书籍等所有资料(最多)，
②微信群：服务于计量经济圈社群群友(最活跃)，
③研究小组：因果推断, 空间计量, 面板数据(最专业)，
④QQ群：2000人大群服务于社群群友(最大)。
计量经济圈是中国计量第一大社区，我们致力于推动中国计量理论和实证技能的提升，圈子以海内外高校研究生和教师为主。计量经济圈六多精神：计量资料多，社会科学数据多，科研牛人多，名校人物多，热情互助多，前沿趋势多。如果你热爱计量并希望长见识，那欢迎你加入到咱们这个大家庭(戳这里)，要不然你只能去其他那些Open access圈子了。注意：进去之后一定要看小鹅社群“群公告”，不然接收不了群息，也不知道怎么进入咱们独一无二的微信群和QQ群。
只有进去之后才能够看见这个群公告