第四节 马尔可夫链蒙特卡洛方法（第1页）

第四节马尔可夫链蒙特卡洛方法

一、基本思想

关于马尔可夫链蒙特卡洛方法需要首先理解以下几个基本概念：随机过程、马尔可夫过程、马尔可夫链、蒙特卡洛方法。

随机过程研究的是随时间演变的随机现象。对于这种现象，一般来说，人们已不能用简单的随机变量或多维随机变量进行合理地表达，而需要用一族随机变量来描述。设T是一无限实数集，我们把依赖于参数t∈T的一族（无限多个）随机变量称为随机过程，记为{X（t），t∈T}，X（t）是一随机变量，T叫作参数集。我们常把t看作时间，称X（t）为时刻t时过程的状态，对于一切t∈T，X（t）所有可能取的一切值的全体称为随机过程的状态空间。对随机过程{X（t），t∈T}进行一次实际试验（即在T上进行一次全程观测），其结果是t的函数，记为x（t），t∈T，称它为随机过程的一个样本函数或样本曲线。当然，试验可以重复进行，所有不同的试验结果构成一族样本函数。因此，随机过程与其样本函数的关系就像数理统计中的总体与样本的关系一样。

关于马尔可夫过程，设随机过程{X（t），t∈T}的状态空间为I，如果对时间t的任意n个数值t1＜t2＜…＜tn，n≥3，ti∈T，在条件X（ti）=xi，xi∈I，i=1，2，…，n-1下，X（tn）的条件分布函数恰等于在条件X（tn-1）=xn-1下X（tn）的条件分布函数，即

则称随机过程{X（t），t∈T}具有马尔可夫性或无后效性，并称此过程为马尔可夫过程。通俗地说，马尔可夫过程是这样一类随机过程，就是在已经知道过程“现在”的条件下，其“将来”不依赖于“过去”。而时间和状态都是离散取值的马尔可夫过程称为马尔可夫链。蒙特卡洛方法是一种使用随机数（或伪随机数）来解决很多计算问题的方法。当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过实际试验的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某类数字特征，并将其作为问题的解。蒙特卡洛方法的求解过程可以归结为三个主要步骤：构造或描述概率过程；从已知概率分布进行抽样；建立各种估计量。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，蒙特卡洛方法是一种求出数值解的有效的方法。

马尔可夫链蒙特卡洛方法，简称MCMC，产生于19世纪50年代早期，是在贝叶斯理论框架下，通过计算机进行模拟的蒙特卡洛方法。该方法将马尔可夫（Markov）过程引入蒙特卡洛模拟中，实现抽样分布随模拟的进行而改变的动态模拟，弥补了传统的蒙特卡洛积分只能静态模拟的缺陷。马尔可夫链蒙特卡洛方法通过模拟的方式对高维积分进行计算，进而使原本异常复杂的高维积分计算问题迎刃而解，使贝叶斯方法仅适用于解决简单低维问题的状况大有改观，为贝叶斯方法的应用开辟了新的道路。

马尔可夫链蒙特卡洛方法的基本思想，首先是通过序列抽样，产生关于模型参数的随机变量序列{X（0），X（1），X（2），…，X（t-1），X（t）}，这个随机变量序列最终收敛到某个平稳分布（stationarydistribution）π（x），该随机变量序列称为马尔可夫链；然后通过该随机变量序列的样本来推论（估计）模型参数的各种特性。因此，该方法有两个关键步骤：一是建立马尔可夫链；二是运用马尔可夫链的样本来估计模型参数的各种特性。

为了便于叙述和理解，下文将把马尔可夫链中的各个值称为观测值。

二、马尔可夫链的建立方法

如何产生收敛于某个平稳分布π（x）的随机变量序列（马尔可夫链），是实现马尔可夫链蒙特卡洛方法的关键。根据马尔可夫链的建立原则，链中的各个观测值的产生过程应该是依次进行的，后一个观测值的产生只依赖于前一个观测值，而与更早产生的观测值无关。也就是说，在建立马尔可夫链过程中的任一时刻t（t≥0），序列中下一时刻t+1处的取值X（t+1）由条件分布P（x|X（t））决定，它只依赖于时刻t处的状态，而与时刻t以前的历史状态无关。因此，马尔可夫链的建立过程就是从一个时刻状态向下一个时刻状态转移的过程，其一步（次）转移概率函数为P（X（t+1）=x′|X（t）=x），该转移概率函数可以记为p（x，x′），称为该马尔可夫链的转移核，通过该转移核产生的且收敛的随机变量序列的分布即为平稳分布π（x）。因此，转移核的构造就成了实现马尔可夫链蒙特卡洛方法的关键，不同的马尔可夫链蒙特卡洛方法，其实质往往也就是构造转移核的方法不同。而在现在的实际应用研究中，大多数马尔可夫链蒙特卡洛方法在构造转移核时，关键就是要写出模型参数（变量）的满条件分布（fullals）。所谓满条件分布，就是在条件分布中，所有的模型参数（变量）应该全部出现，或出现在条件中，或出现在变元中。在关于项目反应理论模型的相关研究文献中一般会提供各相关参数的满条件分布。

由于模型参数的特性是通过建立的马尔可夫链的样本进行估计的，而该马尔可夫链应该收敛到某个平稳分布π（x）上，因此，理论上，该平稳分布应该要服从模型参数的后验分布（posteriordistribution）。

另外，任何一个马尔可夫链的建立都需要有一个起始点，即初始值X（0）。虽然经过研究发现非常不同的初始值最终也能够收敛到平稳分布，但选择一个更加合理的初始值应该能够提高收敛的效果。接着的另外一个问题是，什么时候才能认为抽样过程已经收敛了呢？这虽然没有统一的标准检验判断方法，但是，从经验出发，如果产生的观测值足够多，如6000或更多，那么，它应该是收敛的。另外，也可以通过以下经简化的方法来进行验证，在建立马尔可夫链的过程中，每隔一段计算该段所有观测值的平均数，当所有平均数趋于稳定时，可以认为抽样过程收敛。

马尔可夫链的具体建立方法主要有以下几种：Gibbs抽样方法、Metropolis-Hastings抽样方法、Gibbs结合Metropolis-Hastings方法。

在项目反应理论研究和应用中，Gibbs方法就是从模型参数的满条件分布中直接进行反复迭代抽样，这个过程与项目反应理论中的联合极大似然估计过程类似，即一个参数估计（产生）值作为估计（产生）另一个参数的已知值，如此循环进行迭代更新，因此，Gibbs抽样是单元素抽样，因为每一步只从单个参数的满条件分布中进行抽样。然而，该方法的实现通常涉及计算正态化常数，使样本难以抽取，这个步骤使得该方法显得不够简洁便利。

三、马尔可夫链的使用

建立了服从某个平稳分布的马尔可夫链之后，通常以这个马尔可夫链的样本平均数作为参数的估计值。当然也可以通过这个马尔可夫链的样本来估计参数估计的标准差等。在实现马尔可夫链的过程中，开头的部分观测值会显得非常不稳定，因此，在实践研究当中，通常是以删除了马尔可夫链开头的部分观测值（burn-in）之后剩余的观测值作为研究样本。

下面还是以DINA模型为例，简单介绍利用马尔可夫链蒙特卡洛方法估计项目参数和被试属性掌握模式的过程。这里需要估计的参数为项目的猜测参数s、失误参数g，以及被试的属性掌握模式α。这些参数的联合后验分布可以表示如下：

当给定余下相关的参数，这些参数的完全条件分布可以分别表示如下：

在马尔可夫链蒙特卡洛算法的第t次迭代时，按如下的过程进行抽样。

①对于s和g，从建议分布中抽取它们可能的值，并且按照如下的接受概率进行更新。

②对于α，同样从建议分布中抽取可能的值，并且按照如下的接受概率进行更新。

需要说明的是，虽然这里对于s和g是按照整体来进行抽样的，但并不代表必须这样做，对它们分别进行抽样也是可以的。建议分布的选取理论上是可以任意选择的，但是很多研究表明建议分布的选取会影响到算法的收敛速度。因此，为了算法更好更快地收敛，建议选择使样本相关较小的分布。关于这方面的内容，可以参考相关的文献。

本章小结

参数估计是获得认知诊断结果的技术。参数估计方法的合适性主要是看其实现过程的高效性和估计结果的稳定一致性。

不同的参数估计方法均是基于已有的信息，包括测试样本数据信息和过去的经验信息（先验信息），结合具体的模型，通过寻找极大似然估计值来推断被试特征的过程。测试样本数据信息和先验信息的质量是影响参数估计结果的主要因素。测试样本数据信息的质量受到被试样本质量和容量、测验题目质量和数量的影响，先验信息的质量主要受到我们关于被试特质和题目特征的认知深度的影响。对测试样本信息和先验信息的认识和使用方式也是影响参数估计结果的关键因素，并由此形成了不同的参数估计方法和实现策略。

思考题

1。什么是先验分布？什么是后验分布？

2。边际极大似然估计的一般过程是什么？

3。极大后验估计和期望后验估计方法的关键不同之处是什么？