_马尔可夫链蒙特卡洛_MCMC_方法在估计IRT模型参数中的应用
ＩＲＴ自２０世纪６０年代出现以来，由于其理论模型的科学性和精确性见长，一开始就受到心理和教育测量学的研究者和实际工作者的关注和兴趣。至今已成为考试技术学研究领域中最有影响的一种现代测量理论。但理论的严谨性又导致了计算的复杂性，因而也影响了ＩＲＴ的普及和应用乃至它的考试研究２００６年１０月第２卷第４期ＥｘａｍｉｎａｔｉｏｎｓＲｅｓｅａｒｃｈＯｃｔ．２００６Ｖｏｌ．２，Ｎｏ．４
“马尔可夫链蒙特卡洛”(M CM C)方法在估计IRT
模型参数中的应用[1][2]
王权编译【摘要】本文介绍和阐述怎样运用“马尔可夫链蒙特卡洛”（ＭＣＭＣ）技术，并结合Ｂａｙｅｓ方法来估计ＩＲＴ的模型参数。首先简要地概述了ＭＣＭＣ方法估计模型参数的基本原理；其次介绍ＭＣＭＣ方法估计模型参数的一般方法，涉及Ｇｉｂｂｓ抽样、取舍抽样、Ｍｅｔｒｏｐｏｌｉｓ－Ｈａｓｔｉｎｇｓ算法等概念和方法；最后以ＩＲＴ的“二参数逻辑斯蒂”（２ＰＬ）模型为例，重点介绍了用“Ｇｉｂｂｓ范围内的Ｍ－Ｈ算法”估计项目参数（β１ｊβ２ｊ）的算法过程。结束本文时还解说了ＭＣＭＣ方法的特点。
阅读本文需具有随机过程、Ｍａｒｋｏｖ链、Ｂａｙｅｓ方法等概率论的基本知识。
【关键词】项目反应理论
马尔可夫链蒙特卡洛Ｇｉｂｂｓ抽样取舍抽样作者简介王权，教授，浙江大学教育系。浙江杭州，３１００２８。４５
《考试研究》第２卷第４期
发展速度。令我们欣喜的是在２０世纪９０年代，国外统计学家又推陈出新地提出了参数估计的新方法，使ＩＲＴ的应用和发展又迈出了新的一步。
模型参数的估计是ＩＲＴ的核心内容。以往的参数估计方法主要有“条件极大似然估计”（ＣＭＬＥ）、“联合极大似然估计”（ＪＭＬＥ）、“边际极大似然估计”
（ＭＭＬＥ）和“条件期望—极大化算法”（Ｅ－ＭＡｌｇｏｒｉｔｈｍ）等，大致上后一种算法均是前一种算法的改进［３］。Ｅ－Ｍ算法是由Ｒ．Ｄ．Ｂｏｃｋ和Ｍ．Ａｉｔｋｉｎ于１９８１年创立，它是以ＭＭＬＥ方法为基础发展而成。在Ｅ－Ｍ算法中，Ｅ步要涉及精确的数字积分计算，或者在Ｍ步要涉及偏导计算，当模型较复杂时，计算就十分困难。加之，它还难以将项目参数估计中的“不可靠性”（ｕｎｃｅｒｔａｉｎｔｙ）结合进能力参数估计时不可靠性的计算；反之亦然。
“马尔可夫链蒙特卡洛”（ＭａｒｋｏｖＣｈａｉｎＭｏｎｔｅＣａｒｌｏ，ＭＣＭＣ）方法是一种动态的计算机模拟技术，它是根据任一多元理论分布，特别是根据以贝叶斯（Ｂａｙｅｓ）推断为中心的多元后验分布来模拟随机样本的一种方法。它在估计ＩＲＴ模型参数的应用中，一方面继承了以往估计能力参数和项目参数时所采用的“分而治之”（ｄｉｖｉｄｅ－ａｎｄ－ｃｏｎｑｕｅｒ）的策略，采用能力参数与项目参数交替迭代计算的方法生成Ｍａｒｋｏｖ链；然后采取迥然不同于极大似然方法的思路，充分发挥计算机模拟技术的优势，采集充分大的状态样本，用初等的方法来估计模型参数，绕开了Ｅ－Ｍ算法中的复杂计算，从而提高了估计的成功率。
—“Ｇｉｂｂｓ采样１９９２年统计学家Ｊ．Ｈ．Ａｌｂｅｒｔ首先将一种特殊的ＭＣＭＣ方法——
法”应用于ＩＲＴ问题的研究。现在它已被推广应用于多种复杂的ＩＲＴ模型，在应用于大范围的教育测验评价中尤显它的长处。本文主要介绍ＭＣＭＣ方法的基本原理和基本方法，为说明方便，只列举应用于较为简单状况的二参数逻辑斯蒂模型，它是进一步推广应用的基础。
一、ＭＣＭＣ方法的基本原理
用ＭＣＭＣ方法估计ＩＲＴ的模型参数的基本思路是：首先定义一Ｍａｒｋｏｖ链，Ｍ０，Ｍ１，Ｍ２，…，Ｍｋ，…状态Ｍｋ＝（θｋ，βｋ），ｋ＝１，２，…其中θ为能力参数，β为项目参数，θ和β可以为多维；然后根据Ｍａｒｋｏｖ链模拟观测（即模拟状态）；最后用所得的模拟观测推断参数θ和β。在一定的规则条件下，随着ｋ的增长，状态Ｍｋ的４６
“马尔可夫链蒙特卡洛”（ＭＣＭＣ）方法在估计ＩＲＴ模型参数中的应用
分布收敛到如以下（１）式定义的链的平稳分布π（θ，β）。如果是用Ｂａｙｅｓ方法推
断参数，则需要用平稳的后验分布ｐ（θ，β｜Ｘ）来定义Ｍａｒｋｏｖ链。
Ｍａｒｋｏｖ链的行径是由它的转移核（ｔｒａｎｓｉｔｉｏｎｋｅｒｎｅｌ）ｔ［（θ０，β０），（θ１，β１）］＝Ｐ［Ｍｋ＋１＝（θ１，β１）｜Ｍｋ＝（θ０，β０）］决定的，即由链的现在状态（θ０，β０）转移到新的状态（θ１，
β１）的概率所决定。平稳分布π（θ，β）满足：
∫
θ，βｔ［（θ０，β０），（θ１，β１）］π（θ０，β０）ｄ（θ０，β０）＝π（θ１，β１）（１）
如果我们定义的转移核ｔ［（θ０，β０），（θ１，β１）］导致π（θ，β）＝ｐ（θ，β｜Ｘ），那么删去
为首的ｋ次观测后，留下的“好”观测（θ（１），β（１））＝Ｍｋ＋１，（θ（２），β（２））＝Ｍｋ＋２，…（θ（Ｌ），β（Ｌ））＝Ｍｋ＋Ｌ
就可以被用来推断有关的参数，因为它们的分布就像从Ｐ（θ，β｜Ｘ）抽取的观测
的分布。
二、ＭＣＭＣ的一般方法
（一）吉布斯抽样（ＧｉｂｂｓＳａｍｐｔｉｎｇ）［４］［５］
利用等式（１）经过简短的计算表明转移核
ｔＧ［（θ０，β０），（θ１，β１）］＝Ｐ（θ１｜β０，Ｘ）Ｐ（β１｜θ１，Ｘ）（２）
这是首先由Ｓ．Ｇｅｍａｎ和Ｄ．Ｇｅｍａｎ引入，π（θ，β）＝Ｐ（θ，β｜Ｘ）作为它的平稳分布。用这种方法构造转移核的Ｍａｒｋｏｖ链叫作吉布斯采样法（ＧｉｂｂｓＳａｍｐｔｅｒｓ）；
因子Ｐ（θ｜β，Ｘ）和Ｐ（β｜θ，Ｘ）叫作模型的完全条件分布。根据吉布斯采样法模拟观
测（θｋ，βｋ），即是反复地从完全条件分布抽样；由（θｋ－１，βｋ－１）到（θｋ，βｋ）采取以下两
个转移步骤：
１．抽取θｋ～Ｐ（θ｜Ｘ，βｋ－１）；
２．抽取βｋ～Ｐ（β｜Ｘ，θｋ）。
吉布斯采样仿效标准的ＩＲＴ在参数估计中施行“分而治之”的策略，即在
推断一组参数时，假定其他参数均被固定，而且已知。吉布斯采样通过对各
Ｋ＝１，２，３，…迭代这种“分而治之”的步骤，为其他参数的“不可靠性”（ｕｎｃｅｒ－
ｔａｉｎｔｙ）而调整一组参数的推断，直至要求的模拟大小。实际进行时，我们可对
各Ｋ的采样器（ｓａｍｐｌｅｒ）划分成多于两个转移步骤，在多个转移步骤中一次只
能从β或θ的１个或两个分量中抽样，并以所有其他分量的现在值为条件。
由条件概率的定义可推出：
４７
《考试研究》第２卷第４期
ｐ（θ｜Ｘ，β）＝ｐ（Ｘ｜!，"）ｐ（!，"）
∫ｐ（Ｘ｜!，"）ｐ（!，"）ｄ!，ｐ（β｜Ｘ，θ）＝ｐ（Ｘ｜!，"）ｐ（!，"）
∫ｐ（Ｘ｜!，"）ｐ（!，"）ｄ"
（３）
所以ｐ（θ｜Ｘ，β）和ｐ（β｜Ｘ，θ）都与联合分布ｐ（Ｘ，θ，β）＝ｐ（Ｘ｜!，"）ｐ（!，"）成比例。这种完全条件分布称作吉布斯采样器，它需要计算正规化常数∫ｐ（Ｘ｜θ，β）ｐ（θ，β）ｄθ和∫ｐ（Ｘ｜θ，β）ｐ（θ，β）ｄβ。一些其他的ＭＣＭＣ方法都是为了简化或围绕这些计算而设计的。
（二）数据扩张和吉布斯抽样
Ｍ．Ａ．Ｔａｎｎｅｒ和Ｗ．Ｈ．ｗｏｎｇ（１９８７）将简化复杂混合模型计算的一般方法加以公式化，即把难以运作的模型表示成一个缺失数据的分析性的简单模型的平均，这种方法称作“数据扩张”（Ｄａｔａａｕｇｍｅｎｔａｔｉｏｎ）。数据扩张的简化方式常被应用于简化吉布斯采样中的正规化常数的计算。
假定我们希望建立一个吉布斯采样器，但不能方便地计算完全条件分布ｐ（θ｜β，Ｘ）和ｐ（β｜θ，Ｘ）的正规化常数。然而却有可能将似然函数表示成缺失数据Ｗ时的一个平均（数学期望）：
ｐ（Ｘ｜θ，β）＝∫ｐ（Ｘ，Ｗ｜θ，β）ｄＷ
＝∫ｐ（Ｘ｜θ，β，Ｗ）Ｐ（Ｗ｜θ，β）ｄＷ（４）其中Ｗ假设为不可观测的潜在数据或缺失数据，将Ｘ视为可观测的不完全数据，（Ｘ，Ｗ）才被认为是完全数据。于是类似于等式（２）就可以建立如下的一个３因子转移核的吉布斯采样器：
ｔ［（θ０，β０，Ｗ０）（θ１，β１，Ｗ１）］＝ｐ（θ１｜β０，Ｘ，Ｗ０）Ｐ（β１｜θ１，Ｘ，Ｗ０）Ｐ（Ｗ１｜θ１，β１，Ｘ）
在等式（４）中似然函数的数据扩张表示式的巧妙结构使上式中的３个完全条件分布：
ｐ（θ｜Ｘ，β，Ｗ）＝ｐ（Ｘ｜!，"，Ｗ）ｐ（!，"）ｐ（Ｗ｜!，"）
∫ｐ（Ｘ｜!，"，Ｗ）ｐ（!，"）ｐ（Ｗ｜!，"）ｄ!
、
ｐ（β｜Ｘ，θ，Ｗ）＝ｐ（Ｘ｜!，"，Ｗ）ｐ（!，"）ｐ（Ｗ｜!，"）
∫ｐ（Ｘ｜!，"，Ｗ）ｐ（!，"）ｐ（Ｗ｜!，"）ｄ"
、
ｐ（Ｗ｜Ｘ，θ，β）＝ｐ（Ｘ｜!，"，Ｗ）ｐ（!，"）ｐ（Ｗ｜!，"）
∫ｐ（Ｘ｜!，"，Ｗ）ｐ（!，"）ｐ（Ｗ｜!，"）ｄＷ
所需要的正规化常数，实质上比等式（２）的常数容易计算。当在β和θ上进４８
“马尔可夫链蒙特卡洛”（ＭＣＭＣ）方法在估计ＩＲＴ模型参数中的应用
行推断时，Ｗ的ＭＣＭＣ输出被简单地忽略了。Ｊ．Ｈ．Ａｌｂｅｒｔ（１９９２）等人应用吉布斯采样法于正态卵形模型，其中Ｗ代表呈正态分布的连续型向量，作为构成离散反应的观测基础。
（三）Ｇｉｂｂｓ范围内的取舍抽样［４］［５］
等式（２）中的正规化常数的计算，可以利用“取舍抽样”（ａｃｃｅｐｔａｎｃｅ／Ｒｅｊｅｃ－ｔｉｏｎｓａｍｐｌｉｎｇ）法，即通过抽签的办法来完全加以避免。例如要从分布ｐ（θ｜Ｘ，β）
随机抽取观测值，首先从一“建议性分布”（Ｐｒｏｐｏｓａｌｄｉｓｔｒｉｂｕｔｉｏｎ）ｑ（θ）抽取θ＊，ｑ
（?）是指我们已经知道怎样从中抽签的任一方便分布。然后抛掷一枚头像概
率α＝ｃ?ｐ（θ＊｜Ｘ，β）／ｑ（θ＊）的硬币（即从ｐ＝α的贝努里分布中抽签），其中Ｃ是在约束
条件ｏ!α!１（对所有θ
＊而言）下可以自由选择的１个尽可能大的固定常数———等式（２）中的正规化常数就属于Ｃ。如果硬币出现头像，就接受θ＝θ＊；否则，再抽
取另一θ
＊，连续进行，直到最后出现头像止。类似的方法也可使用于ｐ（β｜Ｘ，θ）。所以“Ｇｉｂｂｓ范围内的取舍抽样”（ａｃｃｅｐｔａｎｃｅ／ＲｅｊｅｃｔｉｏｎｓａｍｐｌｉｎｇｗｉｔｈｉｎＧｉｂｂｓ）法就是指用取舍抽样法来替代吉布斯采样器的一种抽样法。
这种方法的计算速度很慢，特别是θ和β为高维度时，或者建议性分布ｑ与目标分布没有足够好的匹配，在这种情况下，采样时的各个步骤要接受一
个候选值θ＊前，往往需要多次抽签。
（四）Ｍ－Ｈ算法［４］［５］
用取舍抽样法也可以直接构造Ｍａｒｋｏｖ链，这里仍然假定平稳分布π（θ，β）＝ｐ（θ
，β｜Ｘ）。由Ｍｅｔｒｏｐｏｌｉｓ和Ｈａｓｔｉｎｇｓ提出的这种取舍抽样法通常称作“Ｍ－Ｈ算法”（Ｍ－ＨＡｌｇｏｒｉｔｈｍ），它先是根据一更方便的建议性转移核ｑ［（θ０，β０），（θ１，β１）］
产生Ｍａｒｋｏｖ链的一候选（ｃａｎｄｉｄａｔｅ）步骤（θ＊，β＊），随后取（θｋ，βｋ）＝（θ＊，β＊），其接受概
率为：
α［（θ０，β０），（θ＊，β＊）］＝ｍｉｎ!（!＊，"＊）ｑ（!＊，"＊），（!０，"０）!（!０，"０）ｑ（!０，"０），（!＊，"＊），"#
１（５）否则就设（θｋ，βｋ）＝（θｋ－１，βｋ－１）。这样导致的马氏链有转移核：ｔＭ－Ｈ［（θ０，β０），（θ１，β１）］＝ｑ［（θ０，β０），（θ１，β１）］α［（θ０，β０），（θ１，β１）］＋δθ０，β
０（θ１，β１）｛１－∫ｑ
［（θ０，β０），（θ１，β１）］α［（θ０，β０），（θ１，β１）］ｄ（θ１，β１）｝，４９
《考试研究》第２卷第４期
其中δ
θ０，β０
（θ１，β１）是Ｄｉｒａｃｄｅｌｔａ函数，
１－∫ｑ［（θ０，β０），（θ１，β１）］α［（θ０，β０），（θ１，β１）］ｄ（θ１，β１）＝ｐ是不转移至（θ＊，β＊）的概率。可以证明这种Ｍｅｔｒｏｐｏｌｉｓ－Ｈａｓｔｉｎｇｓ链具有平稳分布π（?）。
由于π（?）出现在α的分子和分母中，我们只需知道接近于某个比例常数。当π（?）是后验或完全条件分布时是特别方便的，并且因此对似然函数与先验分布的积成比例。还有，如果我们选择的ｑ是一对称函数，于是，ｑ［（θ＊，β＊），（θ０，β０）］＝ｑ［（θ０，β０），（θ＊，β＊）］，所以建议性分布ｑ退出接受概率α的计算。
Ｓ．ｃｈｉｂ和Ｅ．Ｇｒｅｅｎｂｅｒｇ（１９９５）观测到建议性分布ｑ的选择对Ｍａｒｋｏｖ链达到平稳的速度有很大的影响。有些建议性分布很少有可能生成候补值，这时α系统地低，对于许多步骤，链会趋向于陷入单一的状态。类似地，如果ｑ导致接受概率系统地接近１，链可以取相对少的步骤。
（五）Ｇｉｂｂｓ范围内的Ｍ－Ｈ算法
对于给定一个我们关注的分布π，Ｇｉｂｂｓ抽样法和Ｍ－Ｈ算法各自构建了不同的链，但它们的平稳分布都是π。若将它们的策略和方法结合起来，我们称之谓Ｇｉｂｂｓ范围内的Ｍ－Ｈ法，通常会方便得多，而且结果仍然有平稳分布π。在这种算法中，我们可以根据Ｇｉｂｂｓ算法反复地从完全条件分布中抽取观测值；但如果不是已知完全条件分布的形式和正规化常数，我们就要使用Ｍ－Ｈ算法的单一迭代。
类似于前面说明的直接的Ｇｉｂｂｓ采样法，Ｇｉｂｂｓ范围内的Ｍ－Ｈ算法对各个
转移步骤使用了两个分隔的建议性分布ｑ
θ（θ０，θ１）ｑβ（β０，β１）。
１．试图抽取θｋ～ｐ（θ｜βｋ－１，Ｘ）：
（１）抽取θ＊～ｑθ（θｋ－１，θ）；
（２）接受θｋ－１＝θ＊，具有概率
α（θｋ－１，θ＊）＝ｍｉｎ
ｐ（Ｘ｜"＊，#ｋ－１）ｐ（"＊，#ｋ－１）ｑ"（"＊，"ｋ－１）
ｐ（Ｘ｜"ｋ－１，#ｋ－１）ｐ（"ｋ－１，#ｋ－１）ｑ"（"ｋ－１，"＊）
，
"#１（６）
反之，置θｋ＝θｋ－１。
２．试图抽取β＊～ｐ（β｜θｋ，Ｘ）：
（１）抽取β＊～ｑβ（βｋ－１，β）；
５０
“马尔可夫链蒙特卡洛”（ＭＣＭＣ）方法在估计ＩＲＴ模型参数中的应用
（２）接受βｋ－１＝β＊，具有概率
α（βｋ－１，β＊
）＝ｍｉｎｐ（Ｘ｜#ｋ，!＊）ｐ（#ｋ，!＊）ｑ!（!＊，!ｋ－１）ｐ（Ｘ｜#ｋ，!ｋ－１）ｐ（#ｋ，!ｋ－１）ｑ!（!ｋ－１，!＊），!"１（７）
反之，置βｋ＝βｋ－１。如此得到的Ｍａｒｋｏｖ链有平稳分布π（θ，β）＝ｐ（θ，β｜Ｘ）∞ｐ（Ｘ｜θ，β）ｐ（θ，β）。如果建议性分布ｑθ或ｑβ在自变量上是对称的，那么它在相应接受概率α中的比率也就消除了。
三、
ＭＣＭＣ方法在ＩＲＴ模型参数估计中的应用要获得有效的Ｍａｒｋｏｖ链，即很快地收敛于平稳分布的ｐ（θ，β｜Ｘ）的链，我们必须确定怎样把参数组合进将被一起更新的“程序块”（ｂｌｏｃｋｓ）；怎样为Ｇｉｂｂｓ算法范围内的Ｍ－Ｈ步骤选择有效的建议性分布ｑ。在大多数的ＩＲＴ问题中，自
然地假定参数都是一个先验地独立于另一个，所以ｐ（θ，β）＝ｐ（θ）ｐ（β）＝ΠＩ１ｐｉ（θｉ）
ΠＪ１ｐｊ（βｊ）。这就简化了Ｍ－Ｈ算法中接受概率α
的计算，也简化了程序块的操作。（一）参数的组块（ｂｌｏｃｋｉｎｇ）
前面介绍的ＭＣＭＣ的基本算法使用了一组非常简单的条件分布：在给定项目参数和观测数据条件下所有考生参数向量的条件分布ｐ（θ｜β，Ｘ），及给定考生参数和观测数据条件下所有项目参数向量的条件分布ｐ（β｜θ，Ｘ）。然而这种表示式只是为了便于说明算法，当全参数向量θ和β的维度太高时就难以有效地应用于实际的ＩＲＴ问题。如果模型被合理地参数化，我们就可自由地将参数组合进便于使用的任何参数块。选择相对大的参数块的优点是可以更有效地模拟有相关结构的随机变量组。选择相对小的参数块的优点是较容易调节Ｇｉｂｂｓ范围内的Ｍ－Ｈ算法有合理的接受概率。
对于ＩＲＴ，将参数组合进相应的单个考生或单个项目的参数块比较方
便；在大多数的ＩＲＴ模型中，各θｉ和βｊ大致上是１至３维。
对于ＩＲＴ的模拟作业，我们建议的完全条件分布是ｐ（θｉ｜Ｘ，β，θ＜ｉ，θ＞ｉ）和ｐ（βｊ｜Ｘ，θ，β＜ｊ，β＞ｊ），其中θ＜ｉ＝｛θｉ＇，ｉ′＜ｉ｝，θ＞ｉ＝｛θｉ＇，ｉ′＞ｉ｝，β＜ｊ和β＞ｊ与前类似。因此ＩＲＴ的一般ＭＣＭＣ算法可按以下步骤进行：
１．１抽取θｋ１～ｐ（θ１｜Ｘ，βｋ－１，θｋ－１
＞１）。５１
《考试研究》第２卷第４期
１．２抽取θｋ２～ｐ（θ２｜Ｘ，βｋ－１，θｋ－１＜２，θｋ－１
＞２）。
１．Ｉ抽取θｋ１～ｐ（θＩ｜Ｘ，βｋ－１，θｋ
＜Ｉ）。
２．１抽取βｋ１～ｐ（β１｜Ｘ，θｋ，βｋ－１
＞１）。
２．２抽取βｋ２～ｐ（β２｜Ｘ，θｋ，βｋ＜２，βｋ－１
＞２）。
┇
２．Ｊ抽取βｋＪ～ｐ（βＪ｜Ｘ，θｋ，βｋ
＜Ｊ）。
但用上面的算法直接从完全条件分布抽样通常是比较困难的，所以我们使用Ｇｉｂｂｓ范围内的Ｍ－Ｈ算法，根据建议性分布来完成这些抽样。建议性分布ｑ的自变量是对称的。例如个体ｉ和项目ｊ的参数可计算如下：
１．决定θｋ
ｉ
（１）抽取θ＊ｉ～ｑθ（θｋ－１
ｉ，θｉ）；
（２）接受θｋｉ＝θ＊
ｉ，具有概率
α（θｋ－１ｉ，θ＊ｉ）＝ｍｉｎｐ（Ｘ｜"ｋ－１，!ｋ＜ｉ，!＊ｉ，!ｋ－１＞ｉ）ｐ（!$
ｉ）ｑ!（!$，!ｋ－１）
ｐ（Ｘ｜"ｋ－１，!ｋ－１，!ｋ－１）ｐ（!ｋ－１ｉ）ｑ!（!ｋ－１，!＊），"#
１＝ｍｉｎｐ（Ｘ｜"ｋ－１，!ｋ１，…，!ｋｉ－１，!＊ｉ，!ｋ－１ｉ＋１…，!ｋ－１Ｉ）ｐ（!$
ｉ）
ｐ（Ｘ｜"ｋ－１，!ｋ１，…，!ｋｉ－１，!ｋ－１ｉ，!ｋ－１ｉ＋１，…，!ｋ－１
Ｉ）ｐ（!ｋ－１ｉ），"$
１
＝ｍｉｎｐ（Ｘｉ｜"ｋ－１，!＊ｉ）ｐ（!$
ｉ）ｐ（Ｘｉ｜"ｋ－１，!ｋ－１ｉ）ｐ（!ｋ－１ｉ）
，
%$１（８）
根据因子化（ｆａｃｔｏｒｉｚａｔｉｏｎ）方法，假设上式中所有参数先验独立，ｑθ为对称函数，于是α中的比率就得到了很大的简化。
２．决定βｋ
ｊ
（１）抽取β＊ｊ～ｑβ（βｋ－
１
ｊ，βｊ）；
（２）接受βｋ
ｊ＝β＊ｊ，具有概率５２
“马尔可夫链蒙特卡洛”（ＭＣＭＣ）方法在估计ＩＲＴ模型参数中的应用
α（βｋ－１ｊ，β＊ｊ）＝ｍｉｎｐ（Ｘ｜#ｋ，"＊）ｐ（"＊）ｑ"（"＊，"ｋ－１）ｐ（Ｘ｜#ｋ，"ｋ－１）ｐ（"ｋ－１）ｑ"（"ｋ－１，"＊），!"
１＝ｍｉｎ
ｐ（Ｘ｜#ｋ，"ｋ１，…，"ｋｊ－１，"!ｊ，"ｋ－１ｊ＋１，…，"ｋ－１
Ｊ）ｐ（"!ｊ）ｐ（Ｘ｜#ｋ，"ｋ１，…，"ｋｊ－１，"ｋ－１ｊ，"ｋ－１ｊ＋１，…，"ｋ－１Ｊ）ｐ（"ｋ－１ｊ），!"
１＝ｍｉｎｐ（Ｘｊ｜#ｋ，"＊ｊ）ｐ（"!ｊ）ｐ（Ｘｊ｜#ｋ，"ｋ－１ｊ）ｐ（"ｋ－１ｊ），!"
１（９）若我们运用条件独立性假设、先验独立性假设和ｑβ的对称性，就可导致（９）式中的似然函数得到简约，使α的计算再次得到简化。
（二）建议性密度函数ｑθ和ｑβ的选择
为Ｇｉｂｂｓ范围内的Ｍ－Ｈ算法选择建议性密度函数ｑ，有很大的灵活性。但在一般运用时有两种方法，一种是“随机游动的Ｍ－Ｈ算法”（Ｒａｎｄｏｍ－ｗａｌｋＭ－ＨＡｌｇｏｒｉｔｈｍ），它是以现在的参数状态为中心的对称分布产生候选参数值。例
如我们可以根据正态密度函数Ｎ（θｉ｜μ，σ２）、Ａμ＝θｋ－１ｉ和σ２＝Ｃθ２
ｉ，对各θｉ生成建议性分布
ｑθ（θｋ－１ｉ，θｉ）＝Ｎ（θｉ｜θｋ－１ｉ，Ｃθ２
ｉ）。产生的方差Ｃθ２ｉ（如果θｉ是多维的，Ｃθ２
ｉ将是方差—协方差矩阵）将决定候选参数θ＊ｉ的接受概率。类似地参数β的建议性分布可以选择ｑβ（βｋ－１ｊ，βｊ）＝Ｎ（βｊ｜βｋ－１ｊ，Ｃβ２
ｊ）。两种选择导致等式（６）和（７）中的ｑ被消去，因此等式（８）和（９）的α的表示式得到简化。
产生候选参数的第二种一般方法叫作“独立的Ｍ－Ｈ算法（Ｉｎｄｅｐｅｎｄｅｎｃｅ
Ｍ－ＨＡｌｇｏｒｉｔｈｍ），
因为其中的ｑ不对称，一般不会消去ｑ。例如，我们可以首先对各能力参数作一粗糙的估计θｉ（例如使用标准化总分），
同时还作出相应的方差估计。然后将中心为θｉ、方差为ｃ２#ｉ
的正态分布作为θｉ的建议性分布：ｑθ（#ｋ－１ｉ，θｉ）＝Ｎ（θｉ｜θｉ，ｃ２#ｉ），其中的接受概率将由ｃ２#ｉ决定（当θｉ多维时，ｃ２#ｉ
将是一方差—协方差矩阵）。类似地也可选择ｑβ（"ｋ－１ｊ，βｊ）＝Ｎ（βｊ｜βｊ，ｃ２#ｊ
）。～～～～～～～～～５３
《考试研究》第２卷第４期
这里建议的正态分布可用任何分布来替代，只要在一个特殊问题中能方便地奏效；另外，随机游动和独立性步骤可以在单个链内结合进行。下面所举的二参数模型一例，使用了随机游动的Ｍ－Ｈ算法，在随机游动步骤中被试能力θ和项目难度的建议性分布是对称的正态分布，项目区分度的建议性分布是非对称的对数正态分布。对于估计拉希（Ｒａｓｃｈ）模型的参数，只要在下面介绍的一般算法中简单地去除相应的区分度参数的估计步骤即可。
Ａ．Ｇｅｌｍａｎ等人（１９９６）介绍，选择的方差ｃ２
!
ｉ和ｃ２
"ｊ
对单个参数的抽取要能给
出约５０％的接受比率，而对高维参数组的抽取则要能降至２５％的接受比率。
他们观测到在这样范围内的任何接受比率，Ｍａｒｋｏｖ链会合理地有效。
（三）ＭＣＭＣ方法在二参数逻辑斯蒂模型估计中的应用
根据以上的一般算法，下面介绍“二参数逻辑斯蒂”（ｔｗｏ－ｐａｒａｍｅｔｅｒＬｏｇｉｓｔｉｃ，２ＰＬ）模型Ｇｉｂｂｓ范围内的Ｍ－Ｈ算法，使用的数据抽自“美国教育进步评价”（ＮａｔｉｏｎａｔＡｓｓｅｓｓｍｅｎｌｏｆＥｄｕｃａｔｉｏｎａｔＰｒｏｇｒｅｓｓ，ＮＡＥＰ）阅读部分成绩的试验性州际评价。对ＮＡＥＰ中的短结构反应项目使用了２ＰＬ模型。能力参数θｉ是单维的，ｉ＝１，２，…Ｉ；项目参数是βｊ＝（β１ｊ，β２ｊ）二维的，ｊ＝１，２，…Ｊ；个体ｉ对项目ｊ的反应变量ｘｉｊ服从贝努里分布（即二级评分），其中的正确反应概率由下式给出：
ｐｉｊ（θｉ，βｊ）＝ｐ（ｘｉｊ＝１｜θｉ，β１ｊ，β２ｊ）＝１
１＋ｅｘｐ（－"２ｊ!ｉ＋"１ｊ）
（１０）其中β１ｊ是项目ｊ特征曲线的截距参数（难度），"２ｊ是斜率参数（区分度），（１０）式分母中的指数通常使用更一般的参数化指数－１．７αｊ（θｉ－ｂｊ）。似然函数ｐ（Ｘ｜θ，β）是Ｉ×Ｊ个贝努里项目的积：ｐ（Ｘ｜θ，β）＝ΠｉΠｊｐｉｊ（θｉ，βｊ）ｘｉｊ［１－ｐｉｊ（θｉ，βｊ）］（１－ｘｉｊ）这里使用了与ＮＡＥＰ惯例相一致的做法，指定了如下的独立的先验分布：
ｐ（θｉ）＝Ｎ（θｉ｜０，σ２
θ
）
ｐ（β１ｊ）＝Ｎ（β１ｊ｜０，σ２
β
１?
）
ｐ（β２ｊ）＝ｌｏｇｎｏｒｍａｌ（β２ｊ｜０，σ２
β
２?
）（１１）其中ｉ＝１，２，…Ｉ和ｊ＝１，２，…Ｊ。“ｐ（ｘ）＝ｌｏｇｎｏｒｍａｌ（ｘ｜μ，σ２）”意味着［ｌｏｇ（ｘ）－ｌｏｇ（μ）］５４
“马尔可夫链蒙特卡洛”（ＭＣＭＣ）方法在估计ＩＲＴ模型参数中的应用
／σ为标准正态分布Ｎ（?｜０，１）。ＭＣＭＣ技术的重要优点是容易调控先验规则的选择。现使用随机游动Ｍ－Ｈ算法，它的ｋ次迭代计算如下：
１．试图从ｐ（θ｜βｋ－１，Ｘ）抽取θｋ：
（１）对各ｉ＝１，２，…，Ｉ独立地抽取θ＊ｉ～Ｎ（θｉ｜θｋ－１ｉ，ｃ２θ
）（２）计算接受概率αｉ（θｋ－１ｉ，θ＊ｉ）＝ｍｉｎ｛
１，Ｒθ｝，其中Ｒθ＝ｐ（Ｘｉ｜"＊ｉ／#ｋ－１）ｐ（"＊ｉ）
ｐ（Ｘｉ｜"ｋ－１ｉ，
#ｋ－１）ｐ（"ｋ－１ｉ）＝［!ｊｐｉｊ（"＊ｉ，#ｋ－１ｊ）ｘｉｊ（１－ｐｉｊ（"＊ｉ，#ｋ－１ｊ））１－ｘｉｊ］ｅｘｐ－１２!"（"＊ｉ）２［!ｊｐｉｊ（"ｋ－１ｉ，#ｋ－１ｊ）ｘｉｊ（１－ｐｉｊ（"ｋ－１ｉ，#ｋ－１ｊ））１－ｘｉｊ］ｅｘｐ－１２!"（"ｋ－１ｉ）２
ｉ＝１，２，…，Ｉ。（由于对称，建议性正态分布消去）。
（３）接受各θｋｉ＝θ＊ｉ，具有接受概率αｉ；反之，令θｋｉ＝θｋ－１ｉ。
２．试图从ｐ（β｜θｋ，Ｘ）抽取βｋ：
（１）对各ｊ＝１，２，…，Ｊ独立地抽取β＊１ｊ～Ｎ（β１ｊ｜βｋ－１１ｊ，ｃ２β１?）和β＊２ｊ～ｌｏｇｎｏｒｍａｌ（β２ｊ｜βｋ－１２ｊ，ｃ２β２?
）。（２）计算接受概率αｊ（βｋ－１ｊ，β＊ｊ）＝ｍｉｎ｛
１，Ｒβ｝，其中Ｒβ＝ｐ（Ｘｊ｜θｋ，β＊ｊ）ｐ（β＊ｊ）ｌｏｇｎｏｒｍａｌ（β＊２ｊ｜βｋ－１２ｊ，ｃ２
β２）ｐ（Ｘｊ｜θｋ，βｋ－１ｊ）ｐ（βｋ－１ｊ）ｌｏｇｎｏｒｍａｌ（βｋ－１２ｊ｜β＊２ｊ，ｃ２
β２）＝［!ｉｐｉｊ（"ｋ，
#＊ｊ）ｘｉｊ（１－ｐｉｊ（"ｋ，#＊ｊ））１－ｘｉｊ］［!ｉｐｉｊ（"ｋ，#ｋ－１ｊ）ｘｉｊ（１－ｐｉｊ（"ｋ，#ｋ－１ｊ））１－ｘｉｊ
］×ｅｘｐ－１２σβ１?（β＊１ｊ）２１β＊２ｊｅｘｐ－１２σβ２?（ｌｏｇ（β＊２ｊ））２ｅｘｐ－１２σβ１?（βｋ－１１ｊ）２１βｋ－１２ｊｅｘｐ－１２σβ２?
（ｌｏｇ（βｋ－１２ｊ））２×βｋ－１２ｊβ＊２ｊ，ｊ＝１，２，…，Ｊ。
（３）接受各βｋｊ＝β＊ｊ，具有概率αｊ；反之，令βｋ
ｊ＝βｋ－１。５５
《考试研究》第２卷第４期
比率Ｒβ中的因子βｋ－１２ｊ／β＊
２ｊ是由建议性的对数正态密度函数缺乏对称性而
造成的结果。接受概率受建议性分布中的建议性方差ｃ２θ、ｃ２β１?和ｃ２
β２?所调控。如
果对步骤２加以修改以使β＊１ｊ和β＊
２ｊ联合地抽自二元分布，且该二元分布的相关
结构匹配于β的极大似然估计的Ｆｉｓｈｅｒ信息函数的逆，就可更有效地探索各参数组（β１ｊ，β２ｊ）的建议性分布。采用二元建议性分布比独立地依次抽取β的分量能更好地匹配于β的联合后验分布。
例：ＮＡＥＰ的６个短结构反应项目的参数估计
项目反应数据来自于ＮＡＥＰ的１９９２年试验性的州际评价中的一项目组，内容是四年级的阅读测验。利用这些数据获得的结果可与通常的项目参数估计，例如用ＢＩＬＯＧ软件估计的结果直接进行比较。但应指出：不像通常的ＩＲＴ拟合，ＭＣＭＣ的模拟输出还允许我们立即去估计各考生能力θ的后验分布；不要求对经验的贝叶斯分布追究到底或进行其他的计算。我们还会看到ＭＣＭＣ方法能方便地评价参数估计的标准误、参数估计间的相关等等。
使用加权抽样法从要求回答这组项目的学生中抽取大小为３０００的样本，该学生样本对６个项目的反应结果如表１所示。
（１１）式中的截距和斜率的先验标准差σβ１和σβ２分别采用２和１２，由于容量是３０００的大样本，参数估计对这种先验规定不太敏感。θ的先验标准差σθ采用１。建议性标准差ｃθ、ｃβ１和ｃβ２各取１．１、０．３和０．３。随着这些值的确立，在链的预备
性运行中，对于抽取单维θ产生接近于５０％的接受比率，而对于抽取二维β产生接近于２５％的接受比率。例如，首先试验ｃθ＝１，观测到小于５０％的θ＊的拒绝比率；增加ｃθ到１．１，拒绝比率就移动到较近于５０％。获得合理的拒绝比率并不困难，获得精确的、最优的拒绝比率也是不必要的。一般斜率的初始值采用１，截距的初始值采用０。
表２是ＭＣＭＣ方法的后验平均和标准差与使用ＢＩＬＯＧ软件得到的估计的比较。其中ＭＣＭＣ方法的后验估计是７４００次迭代的Ｍａｒｋｏｖ链的５次重复试验的估计的平均数，两者的结果相当一致。但ＭＣＭＣ方法优于ＢＩＬＯＧ方法的重５６
“马尔可夫链蒙特卡洛”（ＭＣＭＣ）方法在估计ＩＲＴ模型参数中的应用
５７
《考试研究》第２卷第４期
图１项目１的参数向量（β１１，β２１）的后验密度函数（左边）和若干等高线图（右边）要之处是能获得项目参数向量的精确的联合后验分布估计。例如图１．描绘了第１个项目的位置参数β１１和斜率参数β１２的边际联合后验分布的核密度估计。
β１１和β２１之间的后验相关－０．６７是明显的。
图１是第一个项目的截距参数和斜率参数向量（β１１，β２１），经ＭＣＭＣ方法估计后所得的结果，左图是（β１１，β２１）的后验密度函数，右图是左图中的部分等高线（ｃｏｎｔｏｕｒ）。从等高线的形状和坐标中的位置就可推知β１１和β２１的相关是负相关。
图２描绘了１０００条项目特征曲线样本，各条项目特征曲线的项目参数是５８
“马尔可夫链蒙特卡洛”（ＭＣＭＣ）方法在估计ＩＲＴ模型参数中的应用
图２从项目１的联合后验分布ｐ（β１１，β２１｜θ，Ｘ）抽取１０００对参数值（β（ｉ）１１，β（ｉ）
２１），ｉ＝１，２，…，１０００，据此观测作成的１０００条ＩＣＣ图。遮黑部分是项目１的真正ＩＣＣ的９５％后验置信概率的“点态包络”。从图１的联合后验分布抽取所得。也就是图２中的各条曲线代表一不同的抽
样参数对（β１１，β２１）。图中遮黑了真正的项目特征曲线的９５％后验概率的
“点态包络”（Ｐｏｉｎｔｗｉｓｅｅｎｖｅｌｏｐｅ）。通过对抽得的各对项目特征曲线作图，并计算这１０００条项目特征曲线在各θ处的概率值的第２．５％和第９７．５％的百分位数，就可获得这一分数带。
（四）评价收敛和蒙特卡洛标准误
用ＭＣＭＣ方法估计后验分布，重要的是要确定：（１）Ｍａｒｋｏｖ链达到了它的平稳分布；（２）与点估计关联的蒙特卡洛标准误是我们可以接受的小。
不论什么时候应用蒙特卡洛模拟方法估计统计模型中的量数，我们还应区分所获得的估计的“统计不可靠性”（ｓｔａｔｉｃｔｉｃａｌｕｎｃｅｒｔａｉｎｔｙ）和“蒙特卡洛不可靠性”（Ｍｏｎｔｅｃａｒｌｏｕｎｃｅｒｔａｉｎｔｙ）。统计不可靠性告诉我们多大的样本才是提供了被估计的数量的信息，一旦取定了样本和选定了模型，统计不可靠性就固定下来。对于ＭＬＥ统计不可靠性是用Ｆｉｓｈｅｒ信息函数的逆的平方根来度量；而对于贝叶斯估计，则用后验标准差来度量。蒙特卡洛不可靠性产生于模型特征的近似性，例如模型的期望值用蒙特卡洛样本平均数估计。这类不可靠性可用蒙特卡络样本计算合适的标准误来度量。要注意的是蒙特卡洛不可靠性不像统计不可靠性可以简单地通过增大样本来降低不可靠程度。５９
《考试研究》第２卷第４期
评价蒙特卡洛标准误和调查链的平稳性更要小心进行，这里检查了Ｍｏｒｋｏｖ链的几次独立重复试验和一长链的实现过程。短链是７４００步长，除去为首的４００次观测；长链是３７０００步长，除去为首的２０００次观测。表３和表４综合了重复试验的结果。
表３中标记为“短链”的５列，给出了５条短链的各截距参数β１１至β１６和斜率参数β２１至β２６的后验平均的ＭＣＭＣ估计，它们是各链保留下来的７０００个模拟参数值的简单的样本平均数。表４中的相应列，给出了从同样的短链得到的这些参数值的后验标准差的ＭＣＭＣ估计，它们是各短链保留下来的７０００个模拟参数值的简单的样本标准差。这些后验平均和后验标准差的估计是无偏的（ｕｎｂｉａｓｅｄ）和一致的（ｃｏｎｓｉｓｔｅｎｔ），假定已经达到平稳。各表的“ＭＣＭｅａｎ”列中的数值是５条短链的各后验平均数和标准差的平均数。这些平均数就是表２中报告的数值。各表中标记为“长链”的最后列数值是长链保留下来的３５０００个模拟参数值的后验平均数和标准差的样本估计，它们可与短链的结果进行比较。这些结果为我们提供了参数的点估计（后验平均数）和统计不可靠性（后验标准差）。
评价蒙特卡洛不可靠性的一种简单方法是根据所得的短链计算５个后验平均估计（表３）或后验标准差估计（表４）的样本标准差，这些数量报告在各表的“ＭＣＳＥ”列中。估计各平均数的标准误（ＭＣＳＥ）／５
!列在标有“ＭＣ６０
“马尔可夫链蒙特卡洛”（ＭＣＭＣ）方法在估计ＩＲＴ模型参数中的应用
ＳＥＭ”列中。例如β１１的后验平均的９９％置信概率的“蒙特卡洛区间”（Ｍｏｎｔｅ
ｃａｒｌｏｉｎｔｅｒｖａｌ）约是－１．８６８６±３（０．０００９），表示这一后验平均的蒙特卡洛（ＭＣ）估计
精确到３位小数，是好的。
由于链Ｍａｒｋｏｖ从平稳分布提供的是相依样本，所以链的个别观测的样
本标准差一般不能作为“蒙特卡洛标准误”（ＭＣＳＥ）。因此需要先从几个独立
的初始值开始各自模拟相应的几条链，然后再计算这些链的样本平均数的
样本标准差，这就是表３和表４中的“ＭＣＳＥ”的估计。根据相关观测计算“ＭＣ
ＳＥ”的另一种方法是分批处理（ｂａｔｃｈｉｎｇ）；将观测到的长链分隔成等长的相邻
批次（ｂａｔｃｈｅｓ），于是各批次的样本平均就近似地独立（可用自相关作图测量），
这些批次平均数的样本平均数就可作为后验平均的估计，将批次平均数的
样本标准误作为“ＭＣＳＥ”。
上面我们已经概要地介绍了ＭＣＭＣ方法应用于ＩＲＴ模型参数估计的一般
—“Ｇｉｂｂｓ范围内的
原理和方法，并且较为详细地阐述了它的一种特殊方法——
Ｍ－Ｈ算法”（Ｍ－ＨｗｉｔｈｉｎＧｉｂｂｓ），估计２ＰＬ模型参数的计算步骤。Ｍ－Ｈｗｉｔｈｉｎ
Ｇｉｂｂｓ大大降低了计算难度，对于估计模型参数，它的必需的主要知识仅是
ＩＲＴ的似然函数和先验分布等内容。它不像Ｅ－Ｍ算法，偏导计算和Ｎｅｗｔｏｎ－
Ｒａｐｈｓｏｎ迭代计算都是必需的内容，它也不像直接的Ｇｉｂｂｓ抽样，必须计算正
６１
《考试研究》第２卷第４期
规化常数。一般来说，不论用哪一种ＭＣＭＣ方法来估计模型参数，大致上可分作两个阶段，如上例所示，先是利用有关的后验分布产生随机数，从而生成Ｍａｒｋｏｖ链；随后再由生成的Ｍａｒｋｏｖ链反复产生所需的随机样本，由于大样本提供的信息相对来说信度较高，这就容许我们可以采用初等的方法对模型参数作出估计，因而避免了复杂的计算，这在大范围的教育测验数据的评价［６］中更为明显。
然而这种算法复杂性的降低却换来了增加完成计算时间的代价。本文上例中使用Ｓ－Ｐｌｕｓ算法语言运行７４００次迭代计算而获得短链约需５．７小时。其中一个重要原因是对算法编码时没有注意如何最优化；其次，Ｍ－ＨｗｉｔｈｉｎＧｉｂｂｓ中的取舍（ａｃｃｅｐｔａｎｃｅ／ｒｅｊｅｃｔｉｏｎ）过程耗费了较多时间；再有，求解的要求和性质也要影响完成的时间，在Ｅ－Ｍ算法中对各个参数只估计一两个值，而ＭＣＭＣ方法要估计的是全部参数的整个联合后验分布函数。所以提高ＭＣＭＣ方法的估计效率既有较大的空间，但又是一个不可忽视的十分重要的问题。
参考文献：
［１］ＲｉｃｈａｒｄＪ．ＰａｔｚａｎｄＢｒｉａｎＷ．Ｊｕｎｋｅｒ：“ＡＳｔｒａｉｇｈｔｆｏｒｗａｒｄＡｐｐｒｏａｃｈｔｏＭａｒｋｏｖＣｈａｉｎＭｏｎｔｅＣａｒｌｏＭｅｔｈｏｄｓｆｏｒＩｔｅｍＲｅｓｐｏｎｓｅＭｏｄｅｌｓ”，ＪｏｕｒｎａｌｏｆＥｄｕｃａｔｉｏｎａｌａｎｄＢｅｈａｖｉｏｒｉａｌＳｔａｔｉｓｔｉｃｓ，１９９９，Ｖｏｌ．２４，Ｎｏ．２，Ｐ．１４６－１７８．
［２］ＲｉｃｈａｒｄＪ．ＰａｔｚａｎｄＢｒｉａｎＷ．Ｊｕｎｋｅｒ：“ＡｐｐｌｉｃａｔｉｏｎｓａｎｄＥｘｔｅｎｓｉｏｎｓｏｆＭＣＭＣｉｎＩＲＴ：ＭｕｌｔｉｐｅＩｔｅｍＴｙｐｅｓ，ＭｉｓｓｉｎｇＤａｔａ，ａｎｄＲａｔｅｄＲｅｓｐｏｎｓｅｓ”，ＪｏｕｒｎａｌｏｆＥｄｕｃａｔｉｏｎａｌａｎｄＢｅｈａｖｉｏｒｉａｌＳｔａｔｉｓｔｉｃｓ，１９９９，
Ｖｏｌ．２４，Ｎｏ．４，Ｐ．３４２－３６６。
［３］漆书青：《现代测量理论在考试中的应用》，华中师范大学出版社，２００３年版，第６章，第２３３－２６８页。
［４］ＧｅｎｔｌｅＪ．Ｅ：《ＥｌｅｍｅｎｔｓｏｆＣｏｍｐｕｔａｔｉｏｎａｌＳｔａｔｉｓｔｉｃｓ》，科学出版社，２００６年版，ｃｈａｐｔｅｒ２，第３９－６６页。
［５］龚光鲁、钱敏平：《应用随机过程教程》，清华大学出版社，２００３年版，第８章，第１９１－２０２页。
［６］ＲｉｃｈａｒｄＪ．Ｐａｔｚ，ＢｒｉａｎＷ．Ｊｕｎｋｅｒ，ＭａｔｔｈｅｗＳ．Ｊｏｈｎｓｏｎ，ＬｏｕｉｓＴ．Ｍａｒｉａｎｏ：“ＴｈｅＨｉｅｒａｒｃｈｉｃａｌＲａｔｅｒＭｏｄｅｌｆｏｒＲａｔｅｄＴｅｓｔｉｔｅｍｓａｎｄｉｔｓＡｐｐｌｉｃａｔｉｏｎｔｏＬａｒｇｅ－ｓｃａｌＥｄｕｃａｔｉｏｎａｌＡｓｓｅｓｓｍｅｎｔＤａｔａ”，Ｊｏｕｒｎａｌｏｆ
ＥｄｕｃａｔｉｏｎａｌａｎｄＢｅｈａｖｉｏｒｉａｌＳｔａｔｉｓｔｉｃｓ，２００２，Ｖｏｌ．２７，Ｎｏ４，Ｐ．３４１－３８４。
６２
“马尔可夫链蒙特卡洛”（ＭＣＭＣ）方法在估计ＩＲＴ模型参数中的应用
ＡｐｐｌｉｃａｔｉｏｎｓｏｆＭａｒｋｏｖＣｈａｉｎＭｏｎｔｅ
Ｃａｒｌｏ（ＭＣＭＣ）Ｍｅｔｈｏｄｓ
ｉｎＥｓｔｉｍａｔｉｎｇＩＲＴＭｏｄｅｌＰａｒａｍｅｔｅｒｓ
ＷａｎｇＱｕａｎ
ＺｈｅｊａｎｇＵｎｉｖｅｒｓｉｔｙ，Ｚｈｅｊｉａｎｇ，Ｈａｎｇｚｈｏｕ，３１００２８
Ａｂｓｔｒａｃｔ：ＴｈｉｓｐａｐｅｒｄｅｍｏｎｓｔｒａｔｅｓｈｏｗｔｏｅｓｔｉｍａｔｅｔｈｅＩＲＴｍｏｄｅｌｐａｒａｍｅｔｅｒｓｂｙ
ｕｓｉｎｇＭＣＭＣｔｅｃｈｎｉｑｕｅｓ，ａｎｄｉｎｃｏｒｐｏｒａｔｉｎｇＢａｙｅｓｍｅｔｈｏｄｉｎｔｏｔｈｉｓｔｅｃｈｎｉｑｕｅｓ．
Ｔｈｅｆｉｒｓｔ，ｇｉｖｅｓａｂｒｉｅｆａｃｃｏｕｎｔｏｆｔｈｅｇｅｎｅｒａｌｉｔｉｅｓｔｈａｔｅｓｔｉｍａｔｅｓＩＲＴｍｏｄｅｌ
ｐａｒａｍｅｔｅｒｓｂｙｕｓｉｎｇＭＣＭＣｔｅｃｈｎｉｑｕｅｓ；Ｔｈｅｓｅｃｏｎｄ，ｄｅｍｏｎｓｔｒａｔｅｓｅｌｅｍｅｎｔａｒｙ
ｗａｙｓｏｆＭＣＭＣｔｅｃｈｎｉｑｕｅｓｔｈａｔｉｎｖｏｌｖｅｓＧｉｂｂｓｓａｍｐｌｉｎｇ，ａｃｃｅｐｔａｎｃｅ／ｒｅｊｅｃｔｉｏｎ
ｓａｍｐｌｉｎｇ，Ｍｅｔｒｏｐｏｌｉｓ－ＨａｓｔｉｎｇｓＡｌｇｏｒｉｔｈｍａｎｄｓｏｏｎ．Ｆｉｎａｌ，ｔａｋｅｔｈｅｔｗｏ
ｐａｒａｍｅｔｅｒｓＬｏｇｉｓｔｉｃｍｏｄｅｌａｓｅｘａｍｐｌｅ，ｅｍｐｈａｔｉｃａｌｌｙｄｅｍｏｎｓｔｒａｔｅｓｔｈｅａｌｇｏｒｉｔｈｍ
ｐｒｏｃｅｓｓｏｆｅｓｔｉｍａｔｉｎｇｉｔｅｍｐａｒａｍｅｔｅｒｓ（β１ｊ，β２ｊ）ｂｙｕｓｉｎｇＭｅｔｒｏｐｏｌｉｓ－Ｈａｓｔｉｎｇｓ
ａｌｇｏｒｉｔｈｍｗｉｔｈｉｎＧｉｂｂｓ．Ｔｈｅｅｎｄｏｆｔｈｉｓｐａｐｅｒｍａｋｅｃｏｍｍｅｎｔｓｏｎｃｈａｒａｃｔｅｒｏｆ
ＭＣＭＣｔｅｃｈｎｉｑｕｅｓ．
Ｋｅｙｗｏｒｄｓ：ＩｔｅｍＲｅｓｐｏｎｓｅＴｈｅｏｒｙＭａｒｋｏｖ／ＣｈａｉｎＭｏｎｔｅＣａｒｌｏ
Ｇｉｂｂｓｓａｍｐｌｉｎｇａｃｃｅｐｔａｎｃｅ／ｒｅｊｉｃｔｉｏｎｓａｍｐｌｉｎｇ
６３
基于马尔可夫排队模型的行程时间预测方法
第34卷 第4期吉林大学学报(工学版) Vol.34 No.4 2004年10月Journal of Jilin University(Engineering and Technology Edition) Oct.2004 文章编号:1671-5497(2004)04-0671-04 基于马尔可夫排队模型的行程时间预测方法 杨志宏1,杨兆升2,于德新2,陈 林2 (1.宝路集团,吉林长春 130022;2.吉林大学交通学院,吉林长春 130022) 摘 要:针对城市交通流诱导系统(U TF GS)亟待解决的综合路段行程时间预测这一关键问题,利用马尔可夫排队模型给出了车辆路段(含信号交叉口)实时行程时间预测的基本公式,并结合实际工程项目对公式中的一些参数进行了简化,提高了模型的实用性。人工调查数据验证表明该模型具有较高的精度。同时给出了相对误差图。 关键词:交通运输工程;城市交通流诱导系统(U TF GS);马尔可夫排队模型;排队等待时间;实时动态行程时间 中图分类号:U491.2 文献标识码:A T ravel time prediction method based on Malcov queuing model YAN G Zhihong1,YAN G Zhaosheng2,YU Dexin2,CHEN Lin2 (1.China B aolu Com pany,Changchun130022,China;2.College of T ransportation,Jilin U niversity,Changchun 130022,China) Abstract:Aiming at the key problem of synthetic Link travel time prediction in Urban Traffic Flow Guidance System(U TF GS).A Vehicle link travel time prediction algorithm based on Malcov Queuing model was presented.With a quantity of traffic measurement data,some model parameters were simplized and confirmed,thus getting a high precision and also making the model more become applicable. K ey w ords:traffic engineering;U TF GS;Malcov queuing model;queuing wait time;real2time dynamic travel time 0 引 言 交通流诱导以交通流预测和实时动态交通分配(D TA)为基础,应用现代通信技术、电子技术、计算机技术等为路网上的出行者提供必要的交通信息,为其指出当前的最佳行驶路线,从而避免盲目出行造成的交通阻塞,到达路网畅通、高效运行的目的[1,2]。交通流诱导的方式一般分为路边显示板式和车内显示屏式两种。前者主要适用于高速公路以及城市路网集体车辆诱导,后者主要适用于城市路网中的个体车辆诱导[2]。 为了准确、快速地给出路网的最佳行驶路线,需要估计路网中各路段的行程时间。路网中的路段均指含一个相邻的下游交叉口(有信号灯控制)的路段。当车辆进入路段后,其行程时间随交通流量的变 收稿日期:2004205219. 基金项目:“十五”国家智能交通重大科技攻关项目(2002BA404A22B). 作者简介:杨志宏(1971-),男,工程师.E2mail:yangzhihong0527@https://www.wendangku.net/doc/6b17821641.html, 通讯联系人:杨兆升(1938-),男,教授,博士生导师.E2mail:yangzs@https://www.wendangku.net/doc/6b17821641.html,
蒙特卡罗方法的应用【文献综述】
文献综述 信息与计算科学 蒙特卡罗方法的应用 在解决实际问题的时候, 为了模拟某一过程, 产生各种概率分布的随机变量和对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题, 我们应该怎么办? 蒙特·卡罗是一种十分有效的求出数值解的方法. 蒙特卡罗法（ monte-carlo method ）简称M -C 法 通过构造概率模型并对它进行随机试验来解算数学问题的方法. 以计算函数的定积分()()1 0I f x d x =?, ()01f x ≤≤为例, 首先构造一个概率模型: 取一个边长分别为和-的矩形, 并在矩形内随机投点M , 假设随机点均匀地落在整个矩形之内, 当点的掷点数N 充分大时, 则落在图中阴影区内的随机点数与投点总数N 之比M N 就近似等于积分值I . 蒙特卡罗法历史悠久. 1773年法国G.-L.L.von 布丰曾通过随机投针试验来确定圆周率π的近似值, 这就是应用这个方法的最早例子. 蒙特卡罗是摩纳哥著名赌城, 1945年 J.von 诺伊曼等人用它来命名此法, 沿用至今. 数字计算机的发展为大规模的随机试验提供了有效工具, 遂使蒙特卡罗法得到广泛应用. 在连续系统和离散事件系统的仿真中, 通常构造一个和系统特性相近似的概率模型, 并对它进行随机试验, 因此蒙特卡罗法也是系统仿真方法之一. 蒙特卡罗法的步骤是: 构造实际问题的概率模型; ②根据概率模型的特点, 设计和使用降低方差的各类方法, 加速试验的收敛; ③给出概率模型中各种不同分布随机变量的抽样方法; ④统计试验结果, 给出问题的解和精度估计. 概率模型用概率统计的方法对实际问题或系统作出的一种数学描述. 例如对离散事件系统中临时实体的到达时间、永久实体的服务时间的描述（见离散事件系统仿真方法）就是采用概率模型. 虽然由这些模型所确定的到达时间、服务时间可能与具体某一段时间内实际到达时间、服务时间有出入, 但它是通过多次统计获得的结果, 所以从概率分布的规律来说还是相符的. 概率模型不仅可用来描述本身就具有随机特性的问题或系统, 也可用来描述一个确定型问题. 例如参数寻优中的随机搜索法（见动力学系统参数寻优）就是将参数最优化问题构造为一个概率模型, 然后用随机投点、统计分析的方法来进行搜索.
马尔可夫链模型简介
马尔可夫链模型简介 设考察对象为一系统,若该系统在某一时刻可能出现的事件集合为,}{N N E E E E E E ??????,2,1,2,1,两两互斥,则陈i E 为状态。N i ???=,2,1。称该系统从一种状态i E 变化到另一状态j E 的过程称为状态转移，并把整个系统不断实现状态转移的过程称为马尔可夫过程。 定义1 具有下列两个性质的马尔可夫过程称为马尔可夫链： （1）无后效性，即系统的第n 次实验结果出现的状态，只与第1-n 次有关，而与它以前所处的状态无关； （2）具有稳定性，该过程逐渐趋于稳定状态，而与初始状态无关。 定义2 向量),,,(21n u u u u ???= 成为概率向量，如果u 满足： ?? ???=???=≥∑=n j j j u n j u 11,,2,10 定义3 如果方阵P 的每行都为概率向量，则称此方阵为概率矩阵。 如果矩阵A 和B 皆为概率矩阵，则AB ，k A ，k B 也都是概率矩阵（k 为正整数）。 定义4 系统由状态i E 经过一次转移到状态j E 的概率记为ij P ，称矩阵 ????????????????????????=32 12222111211N N N N N P P P P P P P P P P 为一次（或一步）转移矩阵。 转移矩阵必为概率矩阵，且具有以下两个性质： 1、P P P k k )1()(-=； 2、k k P P =)(
其中)(k P 为k 次转移矩阵。 定义5 对概率矩阵P ，若幂次方)(m P 的所有元素皆为正数，则矩阵P 称为正规概率矩阵。（此处2≥m ） 定理1 正规概率矩阵P 的幂次方序列P ，2P ，3P ，…趋近于某一方阵T ，T 的每一行均为同一概率向量t ，且满足t tP = 。 马尔可夫链模型如下： 设系统在0=k 时所处的初始状态 ),,() 0()0(2)0(1)0(N S S S S ???=为已知，经过k 次转移后的状态向量 ),,()()(2)(1)(k N k k k S S S S ???=),2,1(???=k ，则 ??????? ?????? ?????????????=NN N N N N k P P P P P P P P P S S 212222111211)0() ( 此式即为马尔可夫链预测模型。 由上式可以看出，系统在经过k 次转后所处的状态)(k S 取决与它的初始状态)0(S 和转移矩阵P 。 马尔可夫引例 例1：市场占有率预测 设有甲、乙、丙三家企业，生产同一种产品，共同供应1000家用户，各用户在各企业间自由选购，但不超出这三家企业，也无新的用户，假定在10月末经过市场调查得知，甲，乙，丙三家企业拥有的客户分别是：250户，300户，450户，而11月份用户可能的流动情况如下表所示：
基于马尔可夫模型的语言发展趋势预测
基于马尔可夫模型的语言发展趋势预测 发表时间：2019-03-14T15:24:06.727Z 来源：《知识-力量》2019年6月中作者：张浩1 姜晓丽1 朱英豪2 [导读] 为了预测世界语言发展趋势，将语言使用者分为两个部分来分别预测其数量。 （1.华北理工大学建筑工程学院，河北唐山 063210；2.华北理工大学以升教育创新基地，河北唐山 063210）摘要：为了预测世界语言发展趋势，将语言使用者分为两个部分来分别预测其数量。对于母语使用者，根据语言区域的自然增长率和净移民率计算出随时间变化的母语使用者的人数。对于第二或第三语言使用者，将影响使用者人数的三种因子归一化处理，利用层次分析法赋予相应的权重后得到各种语言的发展强度数值。建立马尔可夫预测模型模拟若干年后的第二或第三语言使用者数量，并模拟50年内排名前十四的语言的母语使用者数量的变化趋势。关键词：层次分析法；马尔可夫模型；聚类分析；语言使用者 人类不仅仅只掌握母语这一种语言，越来越多的人开始说第二语言甚至第三语言。在考虑某种语言的总使用人数时，需要在母语使用者人数的基础上加上第二或者第三语言使用者人数。根据可能影响语言的使用的因素，模拟各种语言的使用者随时间变化的分布。建立模型预测在未来50年里，英语的母语使用者的数量和语言的总使用者的数量的变化，并考虑它们是否会被另一种语言替代。 1.模型假设 ●忽略小概率灭绝事件，比如重大自然灾害的影响导致某一语言的灭绝等。 ●在几十年的时间里，各个语言区域都是稳定的发展，不会出现特别大的起伏的情况。 ●假设每个国家的移民一旦定居，他们的子孙都以此国家的官方语言为母语。 2.数量预测模型对于语言使用者数量的预测，我们需要将其分为母语使用者和其它的语言使用者（包括第二和第三语言使用者）两个方向来调查。 2.1母语使用者针对国家而言，母语使用者人数与该国家的居民人数直接相关。根据该国家的移民率，我们可以得到母语使用者人数随时间的变化为： 2.2 总使用者对于一种语言的总使用者人数，我们需要全面考虑它的变化，不仅仅考虑语言区域居民人数的增加或者减少，还需要考虑其它的语言使用者的变化。上文我们已经得知母语使用者的数量随时间的变化，下面我们将解决其它的语言使用者的预测问题。 2.2.1三种影响因子根据上文可得，我们将影响语言发展的因素分为区域的综合实力、商业往来和旅游业的发展状况三个部分。针对这三个部分，我们选取三个指标作为影响因子，分别是区域人均GDP、区域贸易对GDP的贡献度、区域国际游客数量。[1~2] 为进行统一，我们将十种语言的三种影响因子均除以该影响因子中的最大值。将得到的新结果运用层次分析法构造判断矩阵，得出三种影响因子的权重向量分别为0.545、0.272、0.183。我们可以得到关于语言发展强度的方程： 2.2.2马尔科夫模型以其亲代的第二语言作为他的初始状态，余下的九种语言是另外的九种状态，建立马尔科夫预测模型[3]。然后基于语言的发展强度，根据两种语言之间的强度比值来确定一个人的语言从一种状态转移到另一种状态的概率值。定义世界十大母语依次用数字0-9表示其语言状态，由此计算状态转移矩阵。 2.3 模型的应用 2. 3.1英语的语言使用者我们搜集到英语语言区域的平均自然增长率和平均净移民率[4]分别为1.04和0.0039，根据公式1我们可以求解得出英语的母语使用者在五十年以后的数量为：（4）
马尔可夫链蒙特卡罗在实践中的应用
2012年第12期 吉林省教育学院学报 No.12，2012 第28卷JOURNAL OF EDUCATIONAL INSTITUTE OF JILIN PROVINCE Vol .28（总300期） Total No .300 收稿日期：2012—11—14 作者简介：孟庆一（1989—），女，吉林长春人，新加坡籍华人，英国伦敦大学数学系，本科生，研究方向：MCMC 统计学。 浅议马尔可夫链蒙特卡罗在实践中的应用 孟庆一 (英国伦敦大学，英国伦敦) 摘要：本文概括地介绍了马尔可夫链蒙特卡罗(Markov chain Monte Carlo ———MCMC )，一种随机模拟贝叶斯推断的方法。主要的抽样方法包括吉布斯采样(Gibbs Sampling )和Metropolis －Hastings 算法。本文也对MCMC 主题和应用的拓展进行了讨论。 关键词：马尔可夫链;蒙特卡罗;Gibbs 抽样;Metropolis －Hastings 中图分类号：O29 文献标识码：A 文章编号：1671—1580(2012)12—0120—02 统计学中的贝叶斯推理在过去的几十年里有前 所未有的突破，统计学家们发现了一种非常简单，但又非常强大的模拟技术，统称为MCMC 。这种技术可以运用到各种复杂的贝叶斯范例和实际情况。 贝叶斯推理： 贝叶斯方法把所给的模型里所有的未知量的不确定性联系在一起。利用所知的信息，贝叶斯方法用联合概率分布把所有未观察到的数量综合起来，从而得出的推论。在这里，给定已知的未知分布被称为后验分布。有关未知量的推理被称为预测，它们的边缘分布称作为预测分布。 贝叶斯推理根据贝叶斯规则计算后验概率： P （H |E ）= P （E |H ）·P （H ） P （E ）然而，在大多数情况下，所给的模型的复杂性不允许我们运用这个简单的操作。因此，我们需要使用随机模拟， 或蒙地卡罗技术来代替。概述MCMC ： MCMC 采用未知量的高维分布，为难度极高的模拟复杂模型的问题提供了一个答案。 一个马尔可夫链是一个序列的随机变量X 1，X 2，X 3，．．．这个序列有马尔可夫的属性———给予目前的状态，未来和过去的状态是独立的。从数学公 式上看， Pr （X n +1=x |X 1=x 1，X 2=x 2，…，X n =x n ）=Pr （X n +1=x |X n =x n ）X i 的可能的值可数的集合S 称 为链的状态空间。 幸运的是，在马尔可夫链里，我们也有与大数定律和中心极限定理类似的定理。 另外一个问题存在于如何建立一个马尔可夫链的极限分布与所需的分配一模一样。一种可行的解决方案是Gibbs 抽样。它是基于一个马尔可夫链，其前身的依赖性是由模型中出现的条件分布所决定的。另一种可能性是Metropolis －Hastings 算法。它是基于一个马尔可夫链，其前身的依赖性是分裂成两个部分：一个是建议，另一个是接受这一建议。 Metropolis －Hastings 算法： Metropolis －Hastings 算法，可以从任何概率分布中抽取样品，只要求是可计算函数的密度成正比。在贝叶斯的应用程序中，归一化因子计算往往是非常困难的，所以，和其他常用的抽样算法一样，能够在不知道这个比例常数的情况下产生样本是Metropolis －Hastings 算法的重要特征。 该算法的总体思路是产生一系列在一个马尔可 夫链里的样品。在足够长的时间后，所生成的样品的分布与分布相匹配。 该算法基本上按如下方式工作（这是一个特殊 的例子，其建议密度是对称的情况下）：首先，选择一个任意的概率密度Q （x'|x t ），这表明一个新的采样值x'给定样本值x t 。对于简单的Metropolis 算法，这个建议密度必须是对称的Q （x'| 21
浅析蒙特卡洛方法原理及应用
浅析蒙特卡洛方法原理及应用 于希明 （英才学院1236103班测控技术与仪器专业6120110304） 摘要：本文概述了蒙特卡洛方法产生的历史及基本原理，介绍了蒙特卡洛方法的最初应用——蒲丰投针问题求圆周率，并介绍了蒙特卡洛方法在数学及生活中的一些简单应用，最后总结了蒙特卡洛方法的特点。 关键词：蒙特卡洛方法蒲丰投针生活应用 蒙特卡洛方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。它是以概率统计理论为基础, 依据大数定律( 样本均值代替总体均值) , 利用电子计算机数字模拟技术, 解决一些很难直接用数学运算求解或用其他方法不能解决的复杂问题的一种近似计算法。蒙特卡洛方法在金融工程学，宏观经济学，计算物理学（如粒子输运计算、量子热力学计算、空气动力学计算）等领域应用广泛。 一、蒙特卡洛方法的产生及原理 蒙特卡洛方法于20世纪40年代美国在第二次世界大战中研制原子弹的“曼哈顿计划”计划的成员S.M.乌拉姆和J.冯·诺伊曼首先提出。数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种方法，为它蒙上了一层神秘色彩。在这之前，蒙特卡洛方法就已经存在。1777年，法国数学家蒲丰（Georges Louis Leclere de Buffon，1707—1788）提出用投针实验的方法求圆周率π。这被认为是蒙特卡洛方法的起源。 其基本原理如下：由概率定义知，某事件的概率可以用大量试验中该事件发生的频率来估算，当样本容量足够大时，可以认为该事件的发生频率即为其概率。因此，可以先对影响其可靠度的随机变量进行大量的随机抽样，然后把这些抽样值一组一组地代入功能函数式，确定结构是否失效，最后从中求得结构的失效概率。蒙特卡洛法正是基于此思路进行分析的。 设有统计独立的随机变量Xi(i=1，2，3，…，k)，其对应的概率密度函数分别为fx1，fx2，…，fxk，功能函数式为Z=g(x1，x2，…，xk)。首先根据各随机变量的相应分布，产生N组随机数x1，x2，…，xk值，计算功能函数值Zi=g(x1，x2，…，xk)(i=1，2，…，N)，若其中有L组随机数对应的功能函数值Zi≤0，则当N→∞时，根据伯努利大数定理及正态随机变量的特性有：结构失效概率，可靠指标。 二、蒲丰投针问题 作为蒙特卡洛方法的最初应用, 是解决蒲丰投针问题。1777 年, 法国数学家蒲丰提出利用投针实验求解圆周率的问题。设平面上等距离( 如为2a) 画有一些平行线, 将一根长度为2l( l< a) 的针任意投掷到平面上, 针与任一平行线相交的频率为p 。针的位置可以用针的中心坐标x 和针与平行线的夹角θ来决定。任意方向投针, 便意味着x与θ可以任意取一值, 只是0≤x ≤a, 0≤θ≤π。那么, 投针与任意平行线相交的条件为x ≤ l sinθ。相交频率p 便可用下式求
数学建模之马尔可夫预测
马尔可夫预测 马尔可夫过程是一种常见的比较简单的随机过程。该过程是研究一个系统的 状况及其转移的理论。它通过对不同状态的初始概率以及状态之间的转移概率的研究,来确定状态的变化趋势,从而达到对未来进行预测的目的。 三大特点： （1）无后效性 一事物的将来是什么状态，其概率有多大，只取决于该事物现在所处的状态如何，而与以前的状态无关。也就是说，事物第n 期的状态，只与第n 期内的变化和第n-1期状态有关,而与第n-1期以前的状态无关。 （2）遍历性 不管事物现在所处的状态如何，在较长的时间内马尔可夫过程逐渐趋于稳定状态，而与初始状态无关。 （3）过程的随机性。 该系统内部从一个状态转移到另一个状态是，转变的可能性由系统内部的原先历史情况的概率值表示。 1.模型的应用， ①水文预测， ②气象预测， ③地震预测， ④基金投资绩效评估的实证分析， ⑤混合动力车工作情况预测， ⑥产品的市场占有情况预测。 2.步骤 ①确定系统状态 有的系统状态很确定。如：机床工作的状态可划分为正常和故障，动物繁殖后代可以划分为雄性和雌性两种状态等。但很多预测中，状态需要人为确定。如：根据某种产品的市场销售量划分成滞销、正常、畅销等状态。这些状态的划分是依据不同产品、生产能力的大小以及企业的经营策略来确定的，一般没有什么统一的标准。在天气预报中，可以把降水量划分为旱、正常和涝等状态。 ②计算初始概率()0i S 用i M 表示实验中状态i E 出现的总次数，则初始概率为 ()()0 1 1,2,i i i n i i M S F i n M =≈= =∑L ③计算一步转移概率矩阵
令由状态i E 转移到状态j E 的概率为()|ij j i P P E E =，则得到一步转移概率矩阵为： 1112121 2221 2n n n n nn p p p p p p P p p p ??????=??????L L M M M M L ④计算K 步转移概率矩阵 若系统的状态经过了多次转移，则就要计算K 步转移概率与K 步转移概率矩阵。 K 步转移概率矩阵为： 11121212221 2()k n n k n n nn p p p p p p P k p p p p ??????==??????L L M M M M L ⑤预测及分析 根据转移概率矩阵对系统未来所处状态进行预测，即： () ()111210212221 2K n K n n n nn p p p p p p S S p p p ??????=??????L L M M M M L 例题： 设某企业生产洗涤剂为A 型,市场除A 型外,还有B 型、C 型两种。为了生产经营管理上的需要,某企业要了解本厂生产的A 型洗涤剂在未来三年的市场占有倩况。为此,进行了两项工作,一是进行市场调查,二是利用模型进行预测。 市场调查首先全面了解各型洗涤剂在市场占有情况。年终调查结果:市场洗涤剂目前总容量为100万件,其中A 型占40万,B 型和C 型各占30万。 再者,要调杏顾客购买各型洗涤剂的变动情况。调查发现去年购买A 型产品的顾客,今年仍购A 型产品24万件,转购B 型和C 型产品备占8万件,去年购买B 型产品顾客,今年仍购B 型产品9万件,转购A 型15万件,转购C 型6万件,去年购买C 型产品的顾客,今年仍购C 型产品9万件,转购A 型15万件,转购B 型6万件。计算各型产品保留和转购变动率。 模型的建立： ①计算初始概率 用i M 表示i E 型产品出现的总次数，则初始概率为 ()()0 1 1,2,i i i n i i M S F i n M =≈= =∑L （1） ②计算各类产品保留和转购变动率
Markov的各种预测模型的原理与优缺点介绍
Markov的各种预测模型的原理与优缺点介绍 建立有效的用户浏览预测模型，对用户的浏览做出准确的预测，是导航工具实现对用户浏览提供有效帮助的关键。 在浏览预测模型方面，很多学者都进行了卓有成效的研究。AZER提出了基于概率模型的预取方法，根据网页被连续访问的概率来预测用户的访问请求。SARUKKAI运用马尔可夫链进行访问路径分析和链接预测，在此模型中，将用户访问的网页集作为状态集，根据用户访问记录，计算出网页间的转移概率，作为预测依据。SCHECHTER构造用户访问路径树，采用最长匹配方法，寻找与当前用户访问路径匹配的历史路径，预测用户的访问请求。XU Cheng Zhong等引入神经网络实现基于语义的网页预取。徐宝文等利用客户端浏览器缓冲区数据，挖掘其中蕴含的兴趣关联规则，预测用户可能选择的链接。朱培栋等人按语义对用户会话进行分类，根据会话所属类别的共同特征，预测用户可能访问的文档。在众多的浏览模型中，Markov模型是一种简单而有效的模型。Markov模型最早是ZUKERMAN等人于1999年提出的一种用途十分广泛的统计模型，它将用户的浏览过程抽象为一个特殊的随机过程——齐次离散Markov模型，用转移概率矩阵描述用户的浏览特征，并基于此对用户的浏览进行预测。之后，BOERGES等采用了多阶转移矩阵，进一步提高了模型的预测准确率。在此基础上，SARUKKAI建立了一个实验系统[9],实验表明，Markov预测模型很适合作为一个预测模型来预测用户在Web站点上的访问模式。 1 Markov模型 1.1 Markov模型 Markov预测模型对用户在Web上的浏览过程作了如下的假设。 假设1（用户浏览过程假设）：假设所有用户在Web上的浏览过程是一个特殊的随机过程——齐次的离散Markov模型。即设离散随机变量的值域为Web空间中的所有网页构成的集合，则一个用户在Web中的浏览过程就构成一个随机变量的取值序列，并且该序列满足Markov性。 一个离散的Markov预测模型可以被描述成三元组，S代表状态空间；A是转换矩阵，表
基于马尔可夫链蒙特卡洛方法的数据关联算法研究
第３１卷第６期２００７年１２月 武汉理工大学学报（鸯望霾差） ＪｏｕｒｎａｌｏｆＷｕｈａｎＵｎｉｖｅｒｓｉｔｙｏｆＴｅｃｈｎｏｌｏｇｙ （ＴｒａｎｓｐｏｒｔａｔｉｏｎＳｃｉｅｎｃｅ＆Ｅｎｇｉｎｅｅｒｉｎｇ） Ｖ０１．３ｌＮｏ．６ Ｄｅｃ．２００７ 基于马尔可夫链蒙特卡洛方法的 数据关联算法研究＊ 李景熹１’２’王树宗Ｄ王航宇３’ （海军工程大学海军兵器新技术应用研究所＂武汉４３００３３） （海军驻４２６厂军代室２’大连１１６００５）（海军工程大学电子工程学院∞武汉４３００３３） 摘要：数据关联是杂波环境下多目标跟踪问题的难点之一．文中提出了一种基于马尔可夫链蒙特卡洛（ＭＣＭＣ）方法的数据关联算法（ＭＣＭＣＤＡ），该算法通过在相应的关联事件空问中采样，可以有效地估计数据的边际关联概率，而且算法的估计精度可根据需要进行调节．仿真结果表明。在需要跟踪的目标数目较多．探测概率较低、杂波概率较高的情况下，ＪＰＤＡ算法因出现“组合爆炸”问题而难以在实际中应用；ＭＣＭＣＤＡ算法则能在保持较高估计精度的情况下降低计算负荷，从而能够较好地满足实时跟踪系统的要求． 关键词：数据关联；马尔可夫链蒙特卡洛；多目标跟踪；杂波 中图法分类号：ＴＰ３０１．６ ０引言 数据关联是杂波环境下多目标跟踪问题的难点之一．一直以来，联合概率数据关联算法（ＪＰＤＡ）是解决数据关联问题的经典算法Ⅱ砣］．但是，当目标数目增多、杂波概率提高、观测值数目增大时，目标与观测值之间的假设关联事件的数目将呈指数增长甚至出现“组合爆炸”现象，从而极大地加重了（ＪＰＤＡ）算法的计算负荷，使其无法满足实际工程应用的需要． 马尔可夫链蒙特卡洛方法（ＭＣＭＣ）是一种贝叶斯网络计算方法，广泛应用于统计学、经济计量学和计算科学领域，尤其适于处理高维和复杂概率分布问题［３。］．其基本思想是通过建立一个平稳分布为，ｒ（ｚ）的马尔可夫链，通过某种统计采样方法（ＭＨ，Ｇｉｂｂｓ）得到平稳分布万（ｚ）的样本，然后基于这些样本做出各种统计推断． 针对ＪＰＤＡ算法存在的问题，本文提出了一种基于马尔可夫链蒙特卡洛方法的数据关联算法（ＭＣＭＣＤＡ）．该算法以统计抽样思想近似估计关联概率，解决了ＪＰＤＡ算法以解析思想处理数据关联问题所引发的“组合爆炸”问题［５］．仿真算例表明，在保持较高跟踪性能的同时，ＭＣＭＣＤＡ算法大大降低了计算负荷，具有很高的工程实用价值． １问题描述 在杂波环境下多目标跟踪问题中，如果丁个目标的跟踪门出现交集，且交集内有候选回波，这就是典型的数据关联问题．观测值落入跟踪门相交区域，说明某些观测值在不同情况下可能来源于不同目标，数据关联算法的目的就是计算每一个观测值与其可能的各种源目标相关联的概率．设口（志）一｛Ｏｉ（五））筵。为ｋ时刻所有联合事件的集合．式中：巩为侈（奄）中元素的个数 坍（｛） 矾（屉）一ｎ只。（志）（１） 』皇０ 为第ｉ个联合事件．其中：以（志）为观测值７在第ｉ 收稿ａ期｛２００７—０５—２２ 李景熹：男，２８岁，博士生．主要研究领域为武器系统仿真试验技术、机动目标跟踪’国防顼研项目资助（批准号：４０１０８０４０４０１０１） 万方数据
马尔可夫链模型
马尔可夫链模型 马尔可夫链模型（Markov Chain Model） 目录 [隐藏] ? 1 马尔可夫链模型概述 ? 2 马尔可夫链模型的性质 ? 3 离散状态空间中的马尔可夫链 模型 ? 4 马尔可夫链模型的应用 o 4.1 科学中的应用 o 4.2 人力资源中的应用 ? 5 马尔可夫模型案例分析[1] o 5.1 马尔可夫模型的建 立 o 5.2 马尔可夫模型的应 用 ? 6 参考文献 [编辑] 马尔可夫链模型概述 马尔可夫链因安德烈·马尔可夫（Andrey Markov，1856－1922）得名，是数学中具有马尔可夫性质的离散时间随机过程。该过程中，在给定当前知识或信息的情况下，过去（即当期以前的历史状态）对于预测将来（即当期以后的未来状态）是无关的。 时间和状态都是离散的马尔可夫过程称为马尔可夫链, 简记为。 马尔可夫链是随机变量的一个数列。这些变量的范围，即他们所有可能 取值的集合，被称为“状态空间”，而Xn的值则是在时间n的状态。如果Xn + 1对于过去状态的条件概率分布仅是Xn的一个函数，则 这里x为过程中的某个状态。上面这个恒等式可以被看作是马尔可夫性质。
马尔可夫在1906年首先做出了这类过程。而将此一般化到可数无限状态空间是由柯尔莫果洛夫在1936年给出的。 马尔可夫链与布朗运动以及遍历假说这两个二十世纪初期物理学重要课题是相联系的，但马尔可夫寻求的似乎不仅于数学动机，名义上是对于纵属事件大数法则的扩张。 马尔可夫链是满足下面两个假设的一种随机过程： 1、t+l时刻系统状态的概率分布只与t时刻的状态有关，与t时刻以前的状态无关； 2、从t时刻到t+l时刻的状态转移与t的值无关。一个马尔可夫链模型可表示为=(S，P，Q)，其中各元的含义如下： 1）S是系统所有可能的状态所组成的非空的状态集，有时也称之为系统的状态空间，它可以是有限的、可列的集合或任意非空集。本文中假定S是可数集(即有限或可列)。用小写字母i,j(或S i,S j)等来表示状态。 2）是系统的状态转移概率矩阵，其中P ij表示系统在时刻t处于状态i，在下一时刻t+l处于状态i的概率，N是系统所有可能的状态的个数。对于任意i∈s，有 。 3）是系统的初始概率分布，q i是系统在初始时刻处于状态i的概率， 满足。 [编辑] 马尔可夫链模型的性质 马尔可夫链是由一个条件分布来表示的 P(X n + 1 | X n) 这被称为是随机过程中的“转移概率”。这有时也被称作是“一步转移概率”。二、三，以及更多步的转移概率可以导自一步转移概率和马尔可夫性质：
中天会计事务所马尔可夫模型例题(最完整的例题分析)
中天会计事务所马尔可夫模型例题一、问题分析 中天会计事务所由于公司业务日益繁忙，常造成公司事务工作应接不暇，解决该公司出现的这种问题的有效办法是要实施人力资源的供给预测技术。根据对该公司材料的深入分析，可采用马尔可夫模型这一供给预测方法对该事务所的人力资源状况进行预测。 马尔可夫分析法是一种统计方法，其方法的基本思想是：找出过去人力资源变动的规律，用以来推测未来人力变动的趋势。马尔可夫分析法适用于外在环境变化不大的情况下，如果外在环境变化较大的时候这种方法则难以用过去的经验情况预测未来。马尔可夫分析法的分析过程通常是分几个时期来收集数据，然后在得出平均值，利用这些数据代表每一种职位的人员变动频率，就可以推测出人员的变动情况。 二、项目策划 （一）第一步是编制人员变动概率矩阵表。 根据公司提供的内部资料：公司的各职位人员如下表1所示。 表1：各职位人员表 职位代号人数 合伙人P 40 经理M 80 高级会计师S 120 会计员 A 160 制作一个人员变动概率矩阵表，表中的每一个元素表示从一个时期到另一个时期（如从某一年到下一年）在两个工作之间调动的雇员数量的历年平均百分比（以小数表示）。（注：一般以3—5年为周期来估计年平均百分比。周期越长，根据过去人员变动所推测的未来人员变动就越准确。） 表2：历年平均百分比人员变动概率矩阵表 职位合伙人 P 经理M 高级会计师S 会计员A 职位年度离职升为 合伙 人 离职升为经 理 降为 会计 员 离职升为高级 会计师 离职 2005 0.20 0.08 0.13 0.07 0.05 0.11 0.12 0.11 2006 0.23 0.07 0.27 0.05 0.08 0.12 0.15 0.29 2007 0.17 0.13 0.20 0.08 0.03 0.10 0.17 0.20 2008 0.21 0.12 0.21 0.03 0.07 0.09 0.13 0.19 2009 0.19 0.10 0.19 0.02 0.02 0.08 0.18 0.21 平均0.20 0.10 0.20 0.05 0.05 0.10 0.15 0.20
用贝叶斯方法重建基因进化历史
实验3 用贝叶斯方法重建基因进化历史传统的系统进化学研究一般采用的要么是表型的数据，要么是化石的证据。化石的证据依赖于考古学的发现，而表型数据往往极难量化，所以往往会得到许多极具争议的结论。如今，现代分子生物学尤其是测序技术的发展为重建进化史提供了大量的数据，如多态性数据（如SNPs或微卫星）、基因序列、蛋白序列等等。常规的做法一般都是利用某一个或者几个基因来构建物种树（species tree），但是一个基因的进化史能不能完全代表所有被研究物种的进化史呢？这是非常值得讨论的问题，但这不是我们本次实验的重点，在这里就不多赘述了。所以，我们这里所指的进化树如非特别说明，指的都是基因树（gene tree）。 经典的研究系统进化的方法主要有距离法、最大简约法（maximum parsimony，MP）、最大似然法（maximum likelihood，ML）等等。这些方法各有各的优点，也分别有其局限性，例如距离法胜在简单快速、容易理解，但是其模糊化了状态变量，将其简化为距离，也就不可避免的丧失了许多序列本身所提供的信息。而最大简约法虽然用的是原始数据，但也只是原始数据的一小部分。特别是在信息位点比较小的情况下，其计算能力还不如距离法。相对来说，最大似然法虽然考虑问题更加全面，但带来的另一个结果是其计算量大大增加，因此常常需要采用启发式（heuristic）方法推断模型参数，重建进化模型。 本实验利用的是贝叶斯方法来重建基因进化史。 1.贝叶斯方法概述 不可免俗的，我们还是要来看看贝叶斯模型，并分别对模型内部的一系列内容一一进行简单的介绍。 Bayes模型将模型参数视作随机变量（r.v.），并在不考虑序列的同时为参数假设先验分布（prior distribution）。所谓先验分布，是对参数分布的初始化估计。根据Bayes定理，可以不断对参数进行改进： f(θ|D)=f(D|θ)f(θ)f(D)(1) 其中f(θ|D)为后验概率分布（posterior probability distribution），而f(θ)是先验概率分布（prior probability distribution），而f(D|θ)为似然值。此外 f(D)=∫f(D|θ)f(θ) Ωdθ (2)
蒙特卡罗方法及应用实验讲义2016资料
蒙特卡罗方法及应用 实验讲义 东华理工大学核工系 2016.8
实验一 蒙特卡罗方法基本思想 一、实验目的 1、了解蒙特卡罗方法方法的基本思想； 2、掌握蒙特卡罗方法计算面积、体积的方法； 3、掌握由已知分布的随机抽样方法。 二、实验原理 Monte Carlo 方法，又称统计模拟方法或计算机随机模拟方法，是一种基于“随机数”进行数值模拟的方法，一种采用统计抽样理论近似求解物理或数学问题的方法。 如待求量可以表述成某些特征量的期望值、某些事件出现的概率或两者的函数形式，那么可采用蒙特卡罗方法求解。在求解某些特征量的期望值或某些事件出现的概率时，必须构建合符实际的数学模型。例如采用蒙特卡罗方法计算某函数所围面积时，构建的数学模型是构造一已知面积的可均匀抽样区域，在该区域投点，由伯努利定理大数定理可知，进入待求区域投点的频率依概率1收敛于该事件出现的概率（面积之比）。 由已知分布的随机抽样方法指的是由已知分布的总体中抽取简单子样。具体方法很多，详见教材第三章。 三、实验内容 1、安装所需计算工具（MATLAB 、fortran 、C++等）； 2、学习使用rand(m,n)、unifrnd(a,b,m,n)函数 3、求解下列问题： 3.0、蒲丰氏投针求圆周率。 3.1、给定曲线y =2 – x 2 和曲线y 3 = x 2，曲线的交点为：P 1( – 1，1 )、P 2( 1，1 )。曲线围成平面有限区域，用蒙特卡罗方法计算区域面积； 3.2 、计算1z z ?≥??≤??所围体积 其中{(,,)|11,11,02}x y z x y z Ω=-≤≤-≤≤≤≤。 4、对以下已知分布进行随机抽样：
人力供给预测之马尔科夫模型
人力供给预测之马尔科夫模型 马尔科夫模型是根据历史数据，预测等时间间隔点上的各类人员分布状况。此方法的基本思想是根据过去人员变动的规律，推测未来人员变动的趋势。因此，运用马尔科夫模型时假设——未来的人员变动规律是过去变动规律的延续。既是说，转移率要么是一个固定比率，要么可以通过历史数据以某种方式推算出。 步骤： （1）根据历史数据推算各类人员的转移率，得出转移率的转移矩阵； （2）统计作为初始时刻点的各类人员分布状况； （3）建立马尔科夫模型，预测未来各类人员供给状况。 运用马尔科夫模型可以预测一个时间段后的人员分布，虽然这个时间段可以自由定义，但较为普遍的是以一年为一个时间段，因为这样最为实用。在确定转移率时，最粗略的方法就是以今年的转移率作为明年的转移率，这种方法认为最近时间段的变化规律将继续保持到下一时间段。虽然这样很简便，但实际上一年的数据过于单薄，很多因素没有考虑到，一个数据的误差可能非常大。因为以一年的数据得出的概率很难保证稳定，最好运用近几年的数据推算。在推算时，可以采用简单移动平均法、加权移动平均法、指数平滑法、趋势线外推法等，可以在试误的过程中发现哪种方法推算的转移率最准确。尝试用不同的方法计算转移率，然后用这个转移率和去年的数据来推算今年的实际情况，最后选择与实际情况最相符的计算方法。转移率是一类人员转移到另一类人员的比率，计算出所有的转移率后，可以得到人员转移率的转移矩阵。 转移出i类人员的数量 i类人员的转移率 = （3-1） i类人员原有总量 人员转移率的转移矩阵： P11 P12 (1) P21 P22 (2) P = P31 P32 (3) （3-2） ┇┇┇ P K1 P K2 ……P KK 一般是以现在的人员分布状况作为初始状况，所以只需统计当前的人员分布情况即可。这是企业的基本信息，人力资源部门可以很容易地找到这些数据。 建立模型前，要对员工的流动进行说明。流动包括外部到内部、内部之间、内部到外部的流动，内部之间的流动可以是提升、降职、平级调动等。由于推测的是整体情况，个别特殊调动不在考虑之内。马尔科夫模型的基本表达式为：
AI术语
人工智能专业重要词汇表 1、A开头的词汇： Artificial General Intelligence/AGI通用人工智能 Artificial Intelligence/AI人工智能 Association analysis关联分析 Attention mechanism注意力机制 Attribute conditional independence assumption属性条件独立性假设Attribute space属性空间 Attribute value属性值 Autoencoder自编码器 Automatic speech recognition自动语音识别 Automatic summarization自动摘要 Average gradient平均梯度 Average-Pooling平均池化 Accumulated error backpropagation累积误差逆传播 Activation Function激活函数 Adaptive Resonance Theory/ART自适应谐振理论 Addictive model加性学习 Adversarial Networks对抗网络 Affine Layer仿射层 Affinity matrix亲和矩阵 Agent代理/ 智能体 Algorithm算法 Alpha-beta pruningα-β剪枝 Anomaly detection异常检测 Approximation近似 Area Under ROC Curve／AUC R oc 曲线下面积 2、B开头的词汇 Backpropagation Through Time通过时间的反向传播Backpropagation/BP反向传播 Base learner基学习器 Base learning algorithm基学习算法 Batch Normalization/BN批量归一化 Bayes decision rule贝叶斯判定准则 Bayes Model Averaging／BMA贝叶斯模型平均 Bayes optimal classifier贝叶斯最优分类器 Bayesian decision theory贝叶斯决策论 Bayesian network贝叶斯网络
融合马尔科夫链_蒙特卡洛算法的改进通用似然不确定性估计方法在流域水文模型中的应用
2009年4月 水 利 学 报 SHUILI XUEBAO 第40卷 第4期 收稿日期:2008-04-23 基金项目:国家自然科学基金重点项目(40730632);教育部新世纪优秀人才支持计划(NCET-05-0624);霍英东青年教师基金资助 项目(101077) 作者简介:卫晓婧(1984-),女,山西阳泉人,硕士生,主要从事水文水资源方面的研究。E -mail:hellomuki@to https://www.wendangku.net/doc/6b17821641.html, 文章编号:0559-9350(2009)04-0464-10融合马尔科夫链-蒙特卡洛算法的改进通用 似然不确定性估计方法在流域水文模型中的应用 卫晓婧,熊立华,万 民,刘 攀 (武汉大学水资源与水电工程科学国家重点实验室,湖北武汉 430072) 摘要:本文在Blasone 研究工作的基础上,进一步提出了基于马尔科夫链-蒙特卡洛算法的改进通用似然不确定性估计方法(Markov Chain -Monte Carlo based Modified Generalized Likelihood Uncertainty Esti mation,MMGLUE)。该方法结合近年来被广泛用于推求参数后验分布的MC MC 方法,对基于Mon te Carlo 随机取样方法的传统GLUE 方法进行改进,并以预测区间性质最优为标准,对可行参数组阈值进行判断与选择,提出了衡量预测区间对称性的标准,并就预测区间性质与可行参数组个数的相关关系进行了探索。在汉江玉带河流域的实例研究证明,MMGLUE 方法较传统的GLUE 方法能够推求出性质更为优良的预测区间,从而更真实合理地反映水文模型的不确定性。 关键词:MC MC;GLUE;MMGLUE;预测区间;覆盖率;区间宽度;区间对称性 中图分类号:P333文献标识码:A 1 研究背景 近10年来,流域水文模型的不确定性研究逐渐成为当今水文界广泛研究的热点之一,各国的水文学家就此做了大量的工作[1]。Beven [2-3]于1992年率先提出了流域水文模型/异参同效性0的观点,并针对流域水文模型的不确定性研究问题提出了通用似然不确定性估计(Generalized Likelihood Uncertainty Estimation,GLUE)方法。该方法结合Monte Carlo 随机取样技术与Bayesian 框架,对由模型结构、参数冗余及相关性、输入输出误差等因素造成的不确定性进行综合分析。GLUE 方法原理简单,易于操作,但由于其自身理论结构的缺陷,越来越多的研究者就GLUE 方法提出了质疑[4-5],即并非经典的Bayesian 方法、主观判断参数可行域阈值和推求的参数后验概率分布不具有显著的统计特征。因此,基于不同假设的其他不确定性研究方法,如:基于经典Bayesian 理论的Ba RE(Bayesian Recursive Estimation)方法 [6],基于全局卡尔曼滤波理论的EnKF(Ense mble Kalman Filter )方法[7] ,多目标方法如MOSCE M (Mult-i objective Shuffled Complex Evolution Metropolis)方法[8]等被用于估计模型的不确定性工作中。然而,上述方法尽管 理论结构相对复杂,应用效果与GLUE 方法相比却并没有明显的提高。 同时期另一种基于经典Bayesian 理论的马尔科夫链-蒙特卡洛(Markov Chain Monte Carlo,MC MC)方法也被广泛应用于推求参数后验分布的研究中。特别是SCE M -UA (The Shuffled Complex E volution Metropolis Algorithm)方法[9]能够有效地探索参数空间,使Markov Chain 能够朝着高概率密度区进化,从而 推导出具有显著统计特征的水文模型参数的后验分布。 因此,Blasone [10]提出将两种方法结合起来,采用SCE M -UA 采样方法替代传统的GLUE 方法中的) 464)