_马尔可夫链蒙特卡洛_MCMC_方法在估计IRT模型参数中的应用

ＩＲＴ自２０世纪６０年代出现以来，由于其理论模型的科学性和精确性见长，一开始就受到心理和教育测量学的研究者和实际工作者的关注和兴趣。

至今已成为考试技术学研究领域中最有影响的一种现代测量理论。

但理论的严谨性又导致了计算的复杂性，因而也影响了ＩＲＴ的普及和应用乃至它的考试研究２００６年１０月第２卷第４期ＥｘａｍｉｎａｔｉｏｎｓＲｅｓｅａｒｃｈＯｃｔ．２００６Ｖｏｌ．２，Ｎｏ．４“马尔可夫链蒙特卡洛”(M CM C)方法在估计IRT模型参数中的应用[1][2]王权编译【摘要】本文介绍和阐述怎样运用“马尔可夫链蒙特卡洛”（ＭＣＭＣ）技术，并结合Ｂａｙｅｓ方法来估计ＩＲＴ的模型参数。

首先简要地概述了ＭＣＭＣ方法估计模型参数的基本原理；其次介绍ＭＣＭＣ方法估计模型参数的一般方法，涉及Ｇｉｂｂｓ抽样、取舍抽样、Ｍｅｔｒｏｐｏｌｉｓ－Ｈａｓｔｉｎｇｓ算法等概念和方法；最后以ＩＲＴ的“二参数逻辑斯蒂”（２ＰＬ）模型为例，重点介绍了用“Ｇｉｂｂｓ范围内的Ｍ－Ｈ算法”估计项目参数（β１ｊβ２ｊ）的算法过程。

结束本文时还解说了ＭＣＭＣ方法的特点。

阅读本文需具有随机过程、Ｍａｒｋｏｖ链、Ｂａｙｅｓ方法等概率论的基本知识。

【关键词】项目反应理论马尔可夫链蒙特卡洛Ｇｉｂｂｓ抽样取舍抽样作者简介王权，教授，浙江大学教育系。

浙江杭州，３１００２８。

４５《考试研究》第２卷第４期发展速度。

令我们欣喜的是在２０世纪９０年代，国外统计学家又推陈出新地提出了参数估计的新方法，使ＩＲＴ的应用和发展又迈出了新的一步。

模型参数的估计是ＩＲＴ的核心内容。

以往的参数估计方法主要有“条件极大似然估计”（ＣＭＬＥ）、“联合极大似然估计”（ＪＭＬＥ）、“边际极大似然估计”（ＭＭＬＥ）和“条件期望—极大化算法”（Ｅ－ＭＡｌｇｏｒｉｔｈｍ）等，大致上后一种算法均是前一种算法的改进［３］。

Ｅ－Ｍ算法是由Ｒ．Ｄ．Ｂｏｃｋ和Ｍ．Ａｉｔｋｉｎ于１９８１年创立，它是以ＭＭＬＥ方法为基础发展而成。

在Ｅ－Ｍ算法中，Ｅ步要涉及精确的数字积分计算，或者在Ｍ步要涉及偏导计算，当模型较复杂时，计算就十分困难。

加之，它还难以将项目参数估计中的“不可靠性”（ｕｎｃｅｒｔａｉｎｔｙ）结合进能力参数估计时不可靠性的计算；反之亦然。

“马尔可夫链蒙特卡洛”（ＭａｒｋｏｖＣｈａｉｎＭｏｎｔｅＣａｒｌｏ，ＭＣＭＣ）方法是一种动态的计算机模拟技术，它是根据任一多元理论分布，特别是根据以贝叶斯（Ｂａｙｅｓ）推断为中心的多元后验分布来模拟随机样本的一种方法。

它在估计ＩＲＴ模型参数的应用中，一方面继承了以往估计能力参数和项目参数时所采用的“分而治之”（ｄｉｖｉｄｅ－ａｎｄ－ｃｏｎｑｕｅｒ）的策略，采用能力参数与项目参数交替迭代计算的方法生成Ｍａｒｋｏｖ链；然后采取迥然不同于极大似然方法的思路，充分发挥计算机模拟技术的优势，采集充分大的状态样本，用初等的方法来估计模型参数，绕开了Ｅ－Ｍ算法中的复杂计算，从而提高了估计的成功率。

—“Ｇｉｂｂｓ采样１９９２年统计学家Ｊ．Ｈ．Ａｌｂｅｒｔ首先将一种特殊的ＭＣＭＣ方法——法”应用于ＩＲＴ问题的研究。

现在它已被推广应用于多种复杂的ＩＲＴ模型，在应用于大范围的教育测验评价中尤显它的长处。

本文主要介绍ＭＣＭＣ方法的基本原理和基本方法，为说明方便，只列举应用于较为简单状况的二参数逻辑斯蒂模型，它是进一步推广应用的基础。

一、ＭＣＭＣ方法的基本原理用ＭＣＭＣ方法估计ＩＲＴ的模型参数的基本思路是：首先定义一Ｍａｒｋｏｖ链，Ｍ０，Ｍ１，Ｍ２，…，Ｍｋ，…状态Ｍｋ＝（θｋ，βｋ），ｋ＝１，２，…其中θ为能力参数，β为项目参数，θ和β可以为多维；然后根据Ｍａｒｋｏｖ链模拟观测（即模拟状态）；最后用所得的模拟观测推断参数θ和β。

在一定的规则条件下，随着ｋ的增长，状态Ｍｋ的４６“马尔可夫链蒙特卡洛”（ＭＣＭＣ）方法在估计ＩＲＴ模型参数中的应用分布收敛到如以下（１）式定义的链的平稳分布π（θ，β）。

如果是用Ｂａｙｅｓ方法推断参数，则需要用平稳的后验分布ｐ（θ，β｜Ｘ）来定义Ｍａｒｋｏｖ链。

Ｍａｒｋｏｖ链的行径是由它的转移核（ｔｒａｎｓｉｔｉｏｎｋｅｒｎｅｌ）ｔ［（θ０，β０），（θ１，β１）］＝Ｐ［Ｍｋ＋１＝（θ１，β１）｜Ｍｋ＝（θ０，β０）］决定的，即由链的现在状态（θ０，β０）转移到新的状态（θ１，β１）的概率所决定。

平稳分布π（θ，β）满足：∫θ，βｔ［（θ０，β０），（θ１，β１）］π（θ０，β０）ｄ（θ０，β０）＝π（θ１，β１）（１）如果我们定义的转移核ｔ［（θ０，β０），（θ１，β１）］导致π（θ，β）＝ｐ（θ，β｜Ｘ），那么删去为首的ｋ次观测后，留下的“好”观测（θ（１），β（１））＝Ｍｋ＋１，（θ（２），β（２））＝Ｍｋ＋２，…（θ（Ｌ），β（Ｌ））＝Ｍｋ＋Ｌ就可以被用来推断有关的参数，因为它们的分布就像从Ｐ（θ，β｜Ｘ）抽取的观测的分布。

二、ＭＣＭＣ的一般方法（一）吉布斯抽样（ＧｉｂｂｓＳａｍｐｔｉｎｇ）［４］［５］利用等式（１）经过简短的计算表明转移核ｔＧ［（θ０，β０），（θ１，β１）］＝Ｐ（θ１｜β０，Ｘ）Ｐ（β１｜θ１，Ｘ）（２）这是首先由Ｓ．Ｇｅｍａｎ和Ｄ．Ｇｅｍａｎ引入，π（θ，β）＝Ｐ（θ，β｜Ｘ）作为它的平稳分布。

用这种方法构造转移核的Ｍａｒｋｏｖ链叫作吉布斯采样法（ＧｉｂｂｓＳａｍｐｔｅｒｓ）；因子Ｐ（θ｜β，Ｘ）和Ｐ（β｜θ，Ｘ）叫作模型的完全条件分布。

根据吉布斯采样法模拟观测（θｋ，βｋ），即是反复地从完全条件分布抽样；由（θｋ－１，βｋ－１）到（θｋ，βｋ）采取以下两个转移步骤：１．抽取θｋ～Ｐ（θ｜Ｘ，βｋ－１）；２．抽取βｋ～Ｐ（β｜Ｘ，θｋ）。

吉布斯采样仿效标准的ＩＲＴ在参数估计中施行“分而治之”的策略，即在推断一组参数时，假定其他参数均被固定，而且已知。

吉布斯采样通过对各Ｋ＝１，２，３，…迭代这种“分而治之”的步骤，为其他参数的“不可靠性”（ｕｎｃｅｒ－ｔａｉｎｔｙ）而调整一组参数的推断，直至要求的模拟大小。

实际进行时，我们可对各Ｋ的采样器（ｓａｍｐｌｅｒ）划分成多于两个转移步骤，在多个转移步骤中一次只能从β或θ的１个或两个分量中抽样，并以所有其他分量的现在值为条件。

由条件概率的定义可推出：４７《考试研究》第２卷第４期ｐ（θ｜Ｘ，β）＝ｐ（Ｘ｜!，"）ｐ（!，"）∫ｐ（Ｘ｜!，"）ｐ（!，"）ｄ!，ｐ（β｜Ｘ，θ）＝ｐ（Ｘ｜!，"）ｐ（!，"）∫ｐ（Ｘ｜!，"）ｐ（!，"）ｄ"（３）所以ｐ（θ｜Ｘ，β）和ｐ（β｜Ｘ，θ）都与联合分布ｐ（Ｘ，θ，β）＝ｐ（Ｘ｜!，"）ｐ（!，"）成比例。

这种完全条件分布称作吉布斯采样器，它需要计算正规化常数∫ｐ（Ｘ｜θ，β）ｐ（θ，β）ｄθ和∫ｐ（Ｘ｜θ，β）ｐ（θ，β）ｄβ。

一些其他的ＭＣＭＣ方法都是为了简化或围绕这些计算而设计的。

（二）数据扩张和吉布斯抽样Ｍ．Ａ．Ｔａｎｎｅｒ和Ｗ．Ｈ．ｗｏｎｇ（１９８７）将简化复杂混合模型计算的一般方法加以公式化，即把难以运作的模型表示成一个缺失数据的分析性的简单模型的平均，这种方法称作“数据扩张”（Ｄａｔａａｕｇｍｅｎｔａｔｉｏｎ）。

数据扩张的简化方式常被应用于简化吉布斯采样中的正规化常数的计算。

假定我们希望建立一个吉布斯采样器，但不能方便地计算完全条件分布ｐ（θ｜β，Ｘ）和ｐ（β｜θ，Ｘ）的正规化常数。

然而却有可能将似然函数表示成缺失数据Ｗ时的一个平均（数学期望）：ｐ（Ｘ｜θ，β）＝∫ｐ（Ｘ，Ｗ｜θ，β）ｄＷ＝∫ｐ（Ｘ｜θ，β，Ｗ）Ｐ（Ｗ｜θ，β）ｄＷ（４）其中Ｗ假设为不可观测的潜在数据或缺失数据，将Ｘ视为可观测的不完全数据，（Ｘ，Ｗ）才被认为是完全数据。

于是类似于等式（２）就可以建立如下的一个３因子转移核的吉布斯采样器：ｔ［（θ０，β０，Ｗ０）（θ１，β１，Ｗ１）］＝ｐ（θ１｜β０，Ｘ，Ｗ０）Ｐ（β１｜θ１，Ｘ，Ｗ０）Ｐ（Ｗ１｜θ１，β１，Ｘ）在等式（４）中似然函数的数据扩张表示式的巧妙结构使上式中的３个完全条件分布：ｐ（θ｜Ｘ，β，Ｗ）＝ｐ（Ｘ｜!，"，Ｗ）ｐ（!，"）ｐ（Ｗ｜!，"）∫ｐ（Ｘ｜!，"，Ｗ）ｐ（!，"）ｐ（Ｗ｜!，"）ｄ!、ｐ（β｜Ｘ，θ，Ｗ）＝ｐ（Ｘ｜!，"，Ｗ）ｐ（!，"）ｐ（Ｗ｜!，"）∫ｐ（Ｘ｜!，"，Ｗ）ｐ（!，"）ｐ（Ｗ｜!，"）ｄ"、ｐ（Ｗ｜Ｘ，θ，β）＝ｐ（Ｘ｜!，"，Ｗ）ｐ（!，"）ｐ（Ｗ｜!，"）∫ｐ（Ｘ｜!，"，Ｗ）ｐ（!，"）ｐ（Ｗ｜!，"）ｄＷ所需要的正规化常数，实质上比等式（２）的常数容易计算。

当在β和θ上进４８“马尔可夫链蒙特卡洛”（ＭＣＭＣ）方法在估计ＩＲＴ模型参数中的应用行推断时，Ｗ的ＭＣＭＣ输出被简单地忽略了。

Ｊ．Ｈ．Ａｌｂｅｒｔ（１９９２）等人应用吉布斯采样法于正态卵形模型，其中Ｗ代表呈正态分布的连续型向量，作为构成离散反应的观测基础。

（三）Ｇｉｂｂｓ范围内的取舍抽样［４］［５］等式（２）中的正规化常数的计算，可以利用“取舍抽样”（ａｃｃｅｐｔａｎｃｅ／Ｒｅｊｅｃ－ｔｉｏｎｓａｍｐｌｉｎｇ）法，即通过抽签的办法来完全加以避免。

例如要从分布ｐ（θ｜Ｘ，β）随机抽取观测值，首先从一“建议性分布”（Ｐｒｏｐｏｓａｌｄｉｓｔｒｉｂｕｔｉｏｎ）ｑ（θ）抽取θ＊，ｑ（・）是指我们已经知道怎样从中抽签的任一方便分布。

然后抛掷一枚头像概率α＝ｃ・ｐ（θ＊｜Ｘ，β）／ｑ（θ＊）的硬币（即从ｐ＝α的贝努里分布中抽签），其中Ｃ是在约束条件ｏ!α!１（对所有θ＊而言）下可以自由选择的１个尽可能大的固定常数———等式（２）中的正规化常数就属于Ｃ。

如果硬币出现头像，就接受θ＝θ＊；否则，再抽取另一θ＊，连续进行，直到最后出现头像止。

类似的方法也可使用于ｐ（β｜Ｘ，θ）。

所以“Ｇｉｂｂｓ范围内的取舍抽样”（ａｃｃｅｐｔａｎｃｅ／ＲｅｊｅｃｔｉｏｎｓａｍｐｌｉｎｇｗｉｔｈｉｎＧｉｂｂｓ）法就是指用取舍抽样法来替代吉布斯采样器的一种抽样法。

这种方法的计算速度很慢，特别是θ和β为高维度时，或者建议性分布ｑ与目标分布没有足够好的匹配，在这种情况下，采样时的各个步骤要接受一个候选值θ＊前，往往需要多次抽签。

（四）Ｍ－Ｈ算法［４］［５］用取舍抽样法也可以直接构造Ｍａｒｋｏｖ链，这里仍然假定平稳分布π（θ，β）＝ｐ（θ，β｜Ｘ）。

由Ｍｅｔｒｏｐｏｌｉｓ和Ｈａｓｔｉｎｇｓ提出的这种取舍抽样法通常称作“Ｍ－Ｈ算法”（Ｍ－ＨＡｌｇｏｒｉｔｈｍ），它先是根据一更方便的建议性转移核ｑ［（θ０，β０），（θ１，β１）］产生Ｍａｒｋｏｖ链的一候选（ｃａｎｄｉｄａｔｅ）步骤（θ＊，β＊），随后取（θｋ，βｋ）＝（θ＊，β＊），其接受概率为：α［（θ０，β０），（θ＊，β＊）］＝ｍｉｎ!（!＊，"＊）ｑ（!＊，"＊），（!０，"０）!（!０，"０）ｑ（!０，"０），（!＊，"＊），"#１（５）否则就设（θｋ，βｋ）＝（θｋ－１，βｋ－１）。