Bayes统计学与MCMC方法——Metropolis-Hastings(M-H)算法的Matlab

第３５卷第１期　２０１８年２月　华东交通大学学报　Ｊｏｕｒｎａｌ　ｏｆ　Ｅａｓｔ　Ｃｈｉｎａ　Ｊｉａｏｔｏｎｇ　Ｕｎｉｖｅｒｓｉｔｙ　Ｖｏ１．３５　Ｎｏ．１　Ｆｅｂ．，２０１８　文章编号：１００５—０５２３（２０１８）０１—０００１—０８　Ｂａｙｅｓ统计学与ＭＣＭＣ方法　——Ｍｅｔｒ０ｐｏｌｉｓ—Ｈａｓｔｉｎｇｓ（Ｍ—Ｈ）算法的Ｍａｔｌａｂ程序实现　陈梦成，方　苇，杨　超，谢　力　（华东交通大学土木建筑学院，江西南昌３３００１３）　摘要：Ｂａｙｅｓ统计学能够从空中楼阁的理论广泛地落地于自然科学、经济学和社会学等领域，得益于计算机技术和马尔可夫链　蒙特卡洛（Ｍａｒｋ０ｖ　ｃｈａｌｎ　Ｍｏｎｔｅ　Ｃａｒｌｏ．简称ＭＣＭＣ）法的发展。文章介绍了ＭＣＭＣ方法在Ｂａｙｅｓ推断中的应用，主要讨论了　ＭＣＭＣ方法中的独立抽样和随机游走抽样的Ｍｅｔｒｏｐｏｌｉｓ—Ｈａｓｔｉｎｇｓ（Ｍ—Ｈ）算法，利用可读性较强的Ｍａｔｌａｂ程序来实现两种抽样算　法．并给出了详细的程序实施过程，分析了两种抽样的优缺点。模拟分析结果表明：独立抽样Ｍ—Ｈ算法比较容易实施，但是要　求建议分布和后验分布的吻合度较高，否则计算效率低下，而且模拟效果不理想；随机游走抽样的Ｍ—Ｈ算法不需要建议分布　接近后验分布。模拟效果也比较好，因此，克服了独立抽样算法的不足，适用范围更广。　关键词：Ｂａｙｅｓ；ＭＣＭＣ；Ｍ－Ｈ法；Ｍａｔｌａｂ程序　中图分类号：Ｃ８２９．２９　文献标志码：Ａ　笔者最初对Ｂａｙｅｓ统计学推断法产生兴趣，可以追溯到２０１４年对加拿大西安大略大学的访学。在访学　期间。师从国际著名结构工程专家Ｈｏｎｇ　Ｈａｎｐｉｎｇ教授和Ｚｈｏｕ　Ｗｅｎｘｉｎ教授，从事Ｂａｙｅｓ推断法在油气管道　失效评估中的应用研究．并深刻感悟到Ｂａｙｅｓ推断法的卓越性和趣味性。回国以后，笔者一直尝试将Ｂａｙｅｓ　推断法用于全寿命周期内土木工程结构性能退化与剩余寿命评估方面，阅读文献的过程中发现这方面的国　内外研究成果非常有　；因此，决定写一篇有关Ｂａｙｅｓ推断统计学中的热点话题——ＭＣＭＣ方法的文章。　众所周知．ＭＣＭＣ方法是一种利用计算机强大计算能力与仿真能力的Ｂａｙｅｓ推断方法　。然而，即使计算机　能力再强大．由于程序设计有优有劣，用ＭＣＭＣ方法进行统计推断时也会出现有时成功，有时失败的结果。　至今为止．虽然已经有一些统计方面的软件（如ＷｉｎＢＵＧＳ和ＯｐｅｎＢＵＧＳ）包含ＭＣＭＣ方法　”，由于看不到源　码，很难根据实际问题进行二次开发；即使网上公开的一些ＭＣＭＣ源程序，很多情况下也得不到正确解。基　于这个理由，本文在忠实追求ＭＣＭＣ理论框架的同时，用Ｍａｔｌａｂ语言编制了一个可读性较强的源程序，并　对该程序的运行结果进行分析和讨论。　１贝叶斯推断　本节拟对贝叶斯推断法做一个简单介绍。为此，我们先从传统的非Ｂａｙｅｓ统计学开始介绍。以硬币投掷　为例，考虑硬币投掷１次后出现正面的概率０的问题，用随机变量　来表示投掷硬币正反面出现的结果，也　就是考虑当硬币出现正面（取　＝１）和出现反面（Ｘ＝０）的概率问题。由于Ｘ＝Ｉ的概率是　，Ｘ＝Ｏ概率是（１－ｏ），　当随机变量　取　（＝０或１）时的概率可以用下式表示为　Ｐ＝０　（１＿０）１－，　收稿日期：２０１７—１２—２０　（１）　基金项目：国家自然科学基金（５１３７８２０６）　作者简介：陈梦成（１９６２一），男，教授，博士，博士生导师。主持完成国家科技部９７３计划前期研究专项课题，国家自然科学基金　项目．江西省重点自然基金项目等４０余项。获江西省自然科学二等奖２项，科技进步奖１项，江西省高校优秀科技成　果一等奖１项．三等奖ｌ项。出版专著１部，参与编著１部，参与编制江西地方规程１部，获得国家发明专利４项，实　用新型７项。在力学学报，Ｉｎｔｅｒ　Ｊ　Ｎｕｍｅｒ　Ｍｅｔａｏｄｓ和Ｉｎｔｅｒ　Ｊ　Ｓｏｌｉｄｓ　Ｓｔｒｕｃｔ等国内外学术刊物和国内外学术会议上发表论　文２００余篇。ＳＣＩ，ＥＩ，ＩＳＴＰ收录８０余篇。主要研究方向为计算与实验断裂力学，工程结构材料耐久性，组合结构。　

２　华东交通大学学报　２０１８正　从上式可以看出，当　给定时，式（１）为０的函数，又称似然函数。在传统的统计学推断理论中，未知参数０　由似然函数的极值法（极大似然函数法）确定。下面来观察／＇ｂ次硬币投掷的过程。假定ｎ次硬币投掷后，共　有　回正面出现。众所周知，　取观测值　（＝０，１，…，　）时的概率分布服从二项分布，即　Ｐ：　Ｏｘ（１一ｏ）ｏ－￣　（２）　／ｘ／　上式和前述一样。如果　给定的话，是　的函数，是似然函数。根据似然函数的极值法，要使似然函数取得极　值，就是要使式（２）对参数０取得极值，据此，当Ｏ＝ｘｌｎ（参数０的点推断）时似然函数取得最大值。推断值　＝　ｘ／ｎ实际上表示硬币投掷凡次后正面出现的比率，直观上就能得出的这个结论。　下面来考察一下０的点推断值具有什么样的特征？谈到特征，由于ｘｌｎ是一个常数，因此，ｘ／ｎ不应该会　有什么特征。然而，话题不能到此为止。由于　原本是随机变量　的观测值，所以我们应该讨论　ｘ／ｎ的特征。众所周知，　，而不是　是由独立同分布函数生成的样本数据的平均值，根据中心极限定理，当ｎ取得足　够大（ｎ＝３０即可）的时候，可以认为Ｘ／ｎ服从均值为０，方差为０（－ｏ）／ｎ的正态分布。这样，如果　作为０　的点推断值，则可以认为其期望值就是　（无偏估计值），并能构建显著性检验（如：检验０＝１／２是否正确）。　将这个问题作更一般叙述。当０给定时，Ｘ＝ｘ的概率用Ｐ（ｘｌＯ）表示。当　给定时，Ｐ（ｘｌＯ）是似然函数。似　然函数取得最大值时有ｏ＝ｂ，是　的函数，记为　（　）；另外，由于　原本是一个随机变量Ｘ，则　（　）可改写为　（Ｘ）。也是一个随机变量，所以　（　）具有随机变量的各种特征。利用这些特征可以进行推断和检验。　一言以蔽之，在传统统计学中，未知参数０是个常数，并不是随机变量，不存在着分布。然而，在贝叶斯　统计学中是以未知参数０为随机变量作为讨论的出发点，具有不确定性，参数空间　中每个未知参数０值　对应着一个分布。这是传统统计学和贝叶斯统计学的核心区别【１２１，也是传统统计学和贝叶斯统计学之间在哲　学上迄今争论不休的原因。本文不打算卷人他们这种纷争。姑且认为０是个随机变量，它与数据　的联合分　布，使用条件分布后，可以表示成以下形式　Ｐ（　，Ｏ）－－ｐ（￣ｌｏ）ｐ（口）＝ｐ（ＯＩｘ）ｐ（　）　（３）　由上式可以直接得到下式　Ｐ（ＯＩｘ）－－ｐ（ｘｌＯ）ｐ（　）／ｐ（　）∞ｐ（ｘｌＯ）ｐ（０）　（４）　式（４）即为贝叶斯定理，贝叶斯推断中利用贝叶斯定理来连接先验分布和后验分布。其中：Ｐ（０）为参数０的　先验分布，Ｐ（ＯＩｘ）为参数０的后验分布，ｐ（　Ｉ　）为似然函数，ｐ（　）为　的边缘分布。在Ｂａｙｅｓ统计学中，从０的　先验分布出发，用观测数据对先验分布进行修正，修正后的分布表现形式即为后验分布。值得注意的是，此　时　不再是随机变量（后验分布是　为给定时的０条件分布）。另外，由于边缘分布是个定值，所以后验分布　与先验分布和似然函数的乘积成比例。在下述的后验分布抽样中可以不考虑边缘的似然性。后验分布是普　通的概率分布．因此可以直接用这个分布对０进行各种推断（譬如用后验分布的模式和平均值作为０的推　断值）。由此可见，Ｂａｙｅｓ统计学中的后验分布和传统统计学中的　（　）分布属于两个体系，但相对而言，可以　认为它们的概念属地相同。顺便提一下，当先验分布不确定时，可以假定其服从均匀分布。从这个意义上讲，　笔者认为Ｂａｙｅｓ统计学包含了传统统计学内容，而且，在经济学等社会科学中其先验分布的假定通常是按　主观意识进行的，可以很好地将Ｂａｙｅｓ统计学应用于社会科学之中。　综上所述，在Ｂａｙｅｓ统计学中，先验分布是主观的，可以凭借经验与历史事实进行假定，主要认为是依　据先验分布和样本数据确定。但事实上，在很多情况下要用解析方法很难获得后验分布。此时最可行的方法　就是采用程序编制的马尔可夫链蒙特卡洛（ＭＣＭＣ）数值模拟法。　２基于Ｍｅｔｒｏｐｏｌｉｓ—Ｈａｓｔｉｎｇｓ（Ｍ—Ｈ）法的贝叶斯推断　在具体介绍ＭＣＭＣ方法的内容之前，先对贝叶斯统计学和传统统计学中可以共同利用的随机数（样本　信息）的推断方法进行阐述，并且探讨一下如何推断硬币投掷ｎ次后正面出现的概率０问题。如前所述，假　

第１期　陈梦成，等：Ｂａｙｅｓ统计学与ＭＣＭＣ方法　如取正面出现的次数为　，则０的推断值为　＝　。由此可见，　（ｘ）的概率分布是这个推断问题的核心。　前文已经指出，依据中心极限定理，　的分布服从正态分布。然而，对那些不了解中心极限定理的人来　说，是否还有其它方法可以确定　的分布呢？在计算机出现以前，这个问题的回答是否定的，但如今并非如　此，利用ＭＣＭＣ数值模拟方法，即使是一个统计学的门外汉也可以得到　（　）的大概分布形式。其实，求解原　理很简单，即假设真实参数０为适当值。观测硬币投掷　次后的实际结果。用　记录正面出现的次数，由此　可得到一个值Ｘ／ｎ；接着，再投掷１１，次，可得到第二个Ｘ／ｎ值；如此反复投掷硬币可得到无数个Ｘ／ｎ。根据这　些无数个Ｘ／ｎ值，绘出直方图，可以认为这个图大概反映了Ｘ／ｎ的概率分布。当然，随机数的生成实际上没　必要去依赖投掷硬币，而可以全部由计算机来实现。这种模拟的本质就是，通过符合　（　）分布的抽样（计算　硬币实际投掷ｎ次后正面出现的比例）达到了解　（　）的概率分布目的。因此，假如能够实现这种抽样的话，　那么根据这些样本数据就很容易绘出直方图。　与传统统计学中　（　）的概率分布相对应的贝叶斯统计学中　（　）的分布称为后验分布。因此，在贝叶斯　统计学中若能从后验分布抽样，那么依据后验分布进行推断，在理论上是完全可行的。当然，如果能直接用　解析法获得后验分布形式，就没有必要使用数值模拟，这就如同使用中心极限定理获得Ｏ（ｘ）的概率分布后　而没有必要进行数值模拟一样。然而，由于先验分布的假定多源自于分析者的自由构想，所以一般利用解析　法导出后验分布是不可能的，而是求助于从后验分布进行抽样进行数值模拟。当下流行的手法当数ＭＣＭＣ　数值模拟方法。　马尔可夫链的最终状态平稳分布常常是存在的，但是不一定是唯一的。我们利用马尔可夫链希望产生　唯一的平稳分布，这需要产生具有遍历性的马尔可夫链。用ＭＣＭＣ方法进行数值模拟时，首先需要构造一　条马尔可夫链，使其平稳分布为参数　的后验分布，然后通过构造出的这条马尔可夫链生成样本序列，截取　其中已经收敛的有效样本计算Ｂａｙｅｓ推断中的参数期望值。Ｍ—Ｈ算法是ＭＣＭＣ抽样算法的核心，本文通过　Ｍａｔｌａｂ实现了独立抽样和随机游走的Ｍ—Ｈ算法，并对其执行结果进行介绍。　２．１独立抽样的Ｍ—Ｈ算法　首先，从均值为　，方差为１的正态分布中生成１Ｏ个独立的样本数据，其均值为５．０９８　１，然后，采用　ＭＣＭＣ方法对相关参数　进行推断。很显然，　的似然函数服从均值为５．０９８　１，方差为１／１０的正态分布，　的先验分布则假定服从自由度为５的ｔ分布，因此，　的后验分布可以认为是正态分布与ｔ分布之积。于　是，参数　的推断过程可以用三个Ｍａｔｌａｂ函数程序来完成，它们简洁地给出如下：　％　的似然函数　ｆｕｎｃｔｉｏｎ　ｙ＝ｌｉｅｋｌｉｈｏｏｄ（ｍｙｕ，ｄａｔａ）；　ｙ：（１／ｓｑｒｔ（２　ｐｉ　（１／ｄａｔａ（２））））　ｅｘｐ（一（１／２）　（ｍｙｕ—ｄａｔａ（１））＾：！／（１／ｄａｔａｔ（２）））；　ｅｎｄ　％关于参数ｘ的后验分布　ｆｕｎｃｔｉｏｎ　ｙ＝ｐｏｓｔｅｒｉｏｒ（ｘ，ｄａｔａ）；　ｙ：ｐｒｉｏｒ（ｘ）　ｌｉｅｋｌｉｈｏｏｄ（ｘ，ｄａｔａ）；　ｅｎｄ　

４　华东交通大学学报　２０１８焦　需要注意的是在似然函数中参数是由参数“ｍｙｕ”和观测数据“ｄａｔａ”构成的，这正是似然函数完全可以　用参数和观测数据来描述的原因，它也与后验分布（ｐｏｓｔｅｒｉｏｒ函数）相似。此外，程序中的ｄａｔａ（１）存放样本数　据均值５．０９８　１，ｄａｔａ（２）存放样本数据个数１０。　如果能从以上获得的后验分布进行直接抽样的话，问题就简单得多。然而事实并非如此，需要采用Ｍ—　Ｈ抽样算法。Ｍ—Ｈ算法简述如下：首先，采用一个尽可能与后验分布相近的建议分布进行抽样，其次，通过　其它方式给定接受概率，最后，利用这个概率选择接受或者拒绝当前的抽样样本。正如大家所知的那样，这　样得到的样本数据序列是一条经过一定时问（燃烧期）、服从后验分布的样本数据链。本算例拟采用均值　５．０９８　１，方差为１／１０的正态分布作为建议分布。建议分布的抽样程序（函数ｐｒｃｐｏｓａｌ＿ｒａｎ）￣建议分布的密度　函数（函数ｐｒｏｐｏｓａｌ＿ｐｄｆ）程序给出如下：　％建议分布，依据　的先验分布（自由度为５的ｔ分布）与似然函数之积构成的后验分布确定。　％该建议分布与ｍｙｕＯ无关．是一条独立链。　ｆｕｎｃｔｉｏｎ　ｙ＝ｐｒｏｐｏｓａｌｒａｎ（ｍｙｕＯ）；　＿ｇｌｏｂａｌ　ｄａｌａ；　ｙ＝ｄａｔａ（１）＋ｓｑｒｔ（１／ｄａｔａ（２））＊ｒａｎｄｎ（１）；　ｎｄ　％建议分布的密度函数，该建议分布与此前　样本值ｍｙｕＯ的定义是彼此独立的（独立链）。　ｆｕｎｃｔｉｏｎ　ｙ：ｐｒ０ｐｏｓａ１ｐｄｆ（ｍｙｕ０，Ｉｒ【ｙｕ—ｐｒｉｍｅ）；　—ｇｌｏｂａｌ　ｄａｌａ；　ｙ＝ｌ／ｓｑｒｔ（２＇　ｐｉ／ｄａｔａ（２））　：ｅｘｐ（－（１／２）　（ｍｙｕ＿ｐｒｉｍｅ－ｄａｔａ（１））＾２／（１／ｄａｔａ（２）））；　ｅｎｄ　建议分布抽样程序给出了参数“ｍｙｕＯ”，表示上一次从建议分布中抽出的样本值。在很多情况下一般均　设定从建议分布中当前抽出的样本值依赖于上一次抽出的样本值。然而，正如本问题的程序没有使用这个　“ｍｙｕＯ”那样，从建议分布中抽出的样本值不依赖于过去的抽样样本，这就是所谓的独立链。下面给出建议分　布密度函数ｐｒｏｐｏｓａｌ＿ｐｄｆ，其参数是“ｍｙｕＯ”和“ｍｙｕ＿ｐｒｉｍｅ”。“ｍｙｕＯ”前面已经介绍过，“ｍｙｕ＿ｐｒｉｍｅ”表示当前　从建议分布中抽出的样本值。用以上两个参数定义密度函数ｖ。需要注意的是，由于建议分布中使用了独立　链，所以密度函数的表达式中没有出现“ｍｙｕＯ”。　依据某一给定的接受概率选取由建议分布生成的样本程序给出如下：　％ＭＨ算法样本选取程序　ｆｕｎｃｔｉｏｎ　ｙ＝ｓｅｌｅｃｔｉｏｎ（ｍｙｕＯ，ｍｙｕ＿ｐｒｉｍｅ，ｄａｔａ）；　ｕ＝ｒａｎｄ（１，１）；％生成一行一列的（Ｏ，１）之间的均匀随机数　ｎｕｍｅｒａｔｏｒ－＝ｐｏｓｔｅｒｉｏｒ（ｍｙｕ　ｐｒｉｍｅ，ｄａｍ）　ｐｒ０ｐ０ｓａｌ＿ｐｄｆ（ｍｙｕＪ）ｒｉｍｅ，ｍｙｕ０）；　ｄｅｎｏｍｉｎａｔｏｒ＝ｐｏｓｔｅｎｏｒ（ｍｙｕｏ，ｄａｔａ）　ｐｒｏｐｏｓｌａ　ｐｄｆ（ｍｙｕＯ，ｍｙｕ—ｐｒｉｍｅ）；　ｓｅｌｅｃｔｉｏｎｐ＝ｍｉｎ（１，ｎｕｍｅｒａｔｏｒ／ｄｅｎｏｍｉｎａｔｏｒ）；　＿ｉｆ　ｕ＜ｓｅｌｅｅｔｉｏｎｐ　＿ｙ　ｍｙｕ　ｒｉｍｅ；　ｅｌｓｅ　ｙ＝ｍｙｕｏ；　ｅｎｄ　ｅｎｄ　ｓｅｌｅｃｔｉｏｎ是以建议分布的上一次抽样样本“ｍｙｕＯ”、当前抽样样本“ｍｖｕ—ｐｒｉｍｅ”以及观测数据作为参数　的函数。样本数据的本质就是通过接受概率使用后验分布抽样的关系，也就是，用接受概率作为是否接受当　前抽样样本的标准，如果是接受则输出当前抽样样本，否则输出上一次抽样样本。另外，Ｍａｔｌａｂ程序中给｝｝｛　了接受概率的定义，ｓｅｌｅｃｔｉｏｎ＿ｐ函数表示。　

第１期　陈梦成等：Ｂａｙｅｓ统汁学与ＭＣＭＣ办法　，５　ＩＬｌｍｅｒａｌｏｒ＝＝：ｐｏｓｔｅｒｉｏｒ（ｍｙｕ．一＿ｐｒｉｍｅ，ｄａｔａ）　＇ｐｒｏｐｏｓａｌｐｄｆ（ｍｙｔ￣一ｐｒｉｍｅ．，ｍｙｕ０）；　一．ｄｅｎｏｍｉｎａｔｏｒ＝ｐｏｓｔｅｒｉｏｒ（ｍｙｕｏ，ｄａｔａ）　ｐｒｏｐｏｓａｌｐｄｆ（ｍｙｕＯ，￣ｎ　ｙｕｐｒｊｎ１ｅ）；　ｓｅｌｅｃｔｉｏｎｐ＝ｍｉｎ（１，，ｎｕｉｎ￣ｒ廿ｔ０ｒ／ｄｅｎｏｍｉｎａｔｏｒ）：　使用以上程序进行１万次抽样，去掉其　ｌ１初始抽样的１　０００次（燃烧期）后得到的余下样本序列的时间　序列图和直方图分别如图１（ａ）和网１（ｂ）所　马尔町夫链的收敛情况为初始抽样的１０％　奉算例『ｆＩ，图ｌ（ａ）建议分布使用ｒ止态分布，从得到的后验分布的抽样样小数据血方罔　ｌ（ｂ）来　看，我们认为，该图明显地再现了正态分布，冈此，就Ｍ—Ｈ算法而言，模拟结果令人感到非常满意　、但足需要　指　的足，独立抽样的Ｍ—Ｈ算法中要求建议分布　后验分布相近，当它们ｌ叶Ｊ现较大差异时，模拟结果小大　理想，如图ｌ（ｃ）和图１（ｄ）所示。冈此，独立抽样的Ｍ—Ｈ算法作为单独的算法很少使用　，　整个问题的Ｍａｔｌｅｂ主函数程序参见附录１。对于马尔可大链　的燃烧期的选择，在学术上日前还没有一个统一的标准，本示例中燃烧期的选取根据时间序列『皋］ｒ　所示的　ｌ０　３００　９　８　７　２５Ｏ　２００　６　ｓ　４　幽　ｉｊｉｈ　溘　出　ｊ　料　！毫５　１５０　哪　孵　×１０　Ｉ　ｌ　　ｌＩ　ｌ　　ｌ３　２　１（）０　５０　ｌ　０　２　３　４　５　６　７　８　９　１０　３　４　５　６　７　８　迭代次数Ⅳ（　Ｅ念分布）／次　（ａ）Ｍ一－Ｈ法数据时间序列　（独立抽样算法：建议分布为正态分布）　１０　一、　参数值　（ｂ）去除燃烧期后的后验直疗罔　　９　８　７　６　５　４　３　２　１　一　×１０　　ＩＩ　Ｉ　　Ｉｌ　　ｌＩ　ｌ　Ｉ　Ｏ　２　３　４　５　６　７　８　９　１０　迭代次数，ｖ（均匀分布），次　（ｃ）Ｍ—Ｈ法数据时问序列　（独立抽样算法：建议分布为均匀分布）　参数值　（ｄ）去除燃烧期后的后验卣厅图　图１　独立抽样的Ｍ—Ｈ算法的模拟结果（先验分布为ｔ分布）　Ｆｉｇ．１　Ｓｉｍｕｌａｌｉｏｎ　ｅｔ＂Ｍ－ｔＩ　ａｌｇｏｒｉｔｈｍ　ｆ０ｒ　ｉｎｄｅｐｅｎｄｅｎｔ　ｓａｍｐｌｉｎｇ（ｐｒｉｏｒ　ｄｉｓｔｒｉｂｕｔｉｏｎ、ｖｉｔｈ…Ｔ’ｄｉｓｔｒｉｂｕｔｉｏｎ）　２．２随机游走Ｍｅｔｒｏｐｏｌｉｓ算法　Ｍ—Ｈ算法抽样中，除了有独立抽样，还彳丁随机游走抽样。夺节拟讨沦随机游走的Ｍ—ｌ１算法，＿并给ｆｆＪ　祠ｌＪ应的模拟箅例。由于只改变了建议分布，所以似然函数、先验分布和后验分布的Ｍａｔｌａｂ程序无需改变，它　

６　华东交通大学学报　２０１８正　们的Ｍａｔｌａｂ程序均与前述２．１节给出的相同。下面使用随机游走抽样的建议分布中，随机游走的误差项服　从均值为０，方差ｔａｕ—ｓｑｕａｒｅ＝０．５的正态分布（函数ｐｒｏｐｏｓａｌ　ｒａｎ　ｒｗ）：　％建议分布，新的样本是由游走过程生成的。　ｆｕｎｃｔｉｏｎ　ｙ＝ｐｒｏｐｏｓａｌｒａｎｒｗ（ｍｙｕ０）；　＿ｇｌｏｂａｌ　ｔａｕｓｑｕａｒｅ；　＿ｙ：ｍｙｕ０十ｓｑｎ（ｔａｕ—ｓｑｕａｒｅ）母ｒａｎｄｎ（１）；　ｅｎｄ　与之相对应的建议分布的密度函数（函数ｐｍｐｏｓａ１　ｄｆ＿　）如下：　％该建议分布与此前　样本值ｍｙｕ０相关，遵循随机游走模型（游走链）。　ｆｕｎｃｔｉｏｎ　ｙ＝ｐｒｏｐｏｓａｌｐｄｆｒｗ（ｍｙｕ０，ｍｙ＿ｐｒｉｍｅ）；　＿＿ｇｌｏｂａｌ　ｔａｕｓｑｕａｒｅ；　＿ｙ＝ｌ／ｓｑｒｔ（２　ｐｉＡａｕ—ｓｑｕａｒｅ）　ｅｘｐ（一（１／２）　（ｍｙｕ＿ｐｒｉｍｅ—ｍｙｕ０）＇２／ｔａｕ—ｓｑｕａｒｅ）；　ｅｎｄ　依据某一给定的接受概率选取随机游走抽样建议分布生成的抽样样本程序可以完全与上一节的抽样　样本程序一样。但是，由于随机游走误差项的方差ｔａｕ—ｓｑｕａｒｅ为全局变量，所以程序中使用了“ｇｌｏｂａｌ　ｔａｕ—ｓｑｕａｒｅ”语句。利用该程序抽样１万次，去掉初始抽样（燃烧期）１　０００次后的抽样样本序列的时间序列图　和直方图如图２（ａ）和图２（ｂ）所示。　ｌ０　９　８　７　６　也　Ｉ　＾ｔ山Ａ　＿　ｍ　ｋｊ幽　ｉｉ。　ｋ　Ｂ如ｋ　ｉ止　　’Ｉ　－　　。　ｌｓ　呷帑　唧　４　１『｛　’哪Ｐ？　佟　３　２　１　一ｘｌ０：　Ｊ　ｌ　Ｉ　Ｊ　ｌ　　ｌｌ　Ｉ　ｌ　０　迭代次数Ⅳ（正态分布），次　（ａ）Ｍ—Ｈ法数据时间序列　（随机移动的Ｍ—Ｈ算法：误差项正态分布）　参数值　（ｂ）去除燃烧期后的后验直方图（误差项正态分布）　图２随机游走（误差项正态分布）Ｍ－Ｈ算法的模拟结果（先验分布为ｔ分布）　Ｆｉｇ．２　Ｓｉｍｕｌａｔｉｏｎ　ｏｆ　Ｍ—Ｈ　ａｌｇｏｒｉｔｈｍ　ｆｏｒ　ｒａｎｄｅｍ　ｗａｌｋ　ｓａｍｐ￣ｎｇ（ｅｒｒｏｒ　ｔｅｒｍ　ｍｏｒｍａｌ　ｄｉｓｔｒｉｂｕｔｉｏｎ）　注意到尽管图２中的建议分布使用的是随机游走分布，与图１中使用的正态分布完全不同，但得到的　抽样结果直方图几乎完全相同。虽说用的是随机游走。但最终得到的建议概率分布也不过是服从以上一次　样本数据为均值，０．５为方差的正态分布。从这个意义上讲，难道建议分布摆脱不了正态分布的魔咒？为了防　止读者因此而对该抽样方法感到疑惑，下面再给出误差项服从样本区间（一０．５，０．５）的均匀分布的随机游走　建议分布的模拟算例。与上一算例的不同之处在于建议分布的抽样程序和建议分布密度函数程序，它们分　别列出　：下：　％误　项服从均勺分布的随机游走链的建议分布，样本生成来自于随机游走过程Ｌ腿机游走链）。　ｆｕｎｃｔｉｏｎ　ｙ＝ｐｒｏｐｏｓａｌｒａｎ—＿ｒｗ＿ｕｎｉｆ（ｍｙｕＯ）；　ｙ＝ｍｙｕ０＋ｒａｎｄ（１，１）－０．５；　ｅｎｄ　

理　０　９　８　７　６　５　４　３　２　ｌ　Ｏ　第１期　陈梦成．等：Ｂａｙｅｓ统计学与ＭＣＭＣ方法　７　与其对应的建议分布的密度函数为常数（函数ｐｒｏｐｏｓａｌｐｄｆｒｗ），程序如下：　＿＿ｆｕｎｃｔｉｏｎ　ｙ＝ｐｒｏｐｏｓａｌｐｄｆｒｗ＿＿＿ｕｎｉｆ（ｍｙｕＯ，ｍｙ＿ｐｒｉｍｅ）；　ｙ＝１；　ｅｎｄ　建议分布的密度函数“ｐｒｏｐｏｓａｌ＿ｐｄｆ　ｒｗ　ｕｎｉｆ函数”反映了建议分布的均匀随机性，是个常数。这个程序　运行结果如下图３（ａ）和图３（ｂ）所示。图３完美再现了后验分布，再一次显示出Ｍ—Ｈ法的强大能力。　；　｝　矗　池＾　ｊ　．ｈ　ｉｊｉ｝ｊ埘　１孵　ｄ妇　＆捌　血ｉｉ矗娃ｊ　Ｉ　研哪　－、　１　｝　１　甜　一　×１Ｃ　　ｌＩ　ｌ　Ｉ　ｌ　Ｉ　Ｉ　ｌ　１　２　３　４　５　６　７　８　９　１０　迭代次数Ⅳ（均匀分布）／次　（ａ）Ｍ—Ｈ法数据时间序列　（随机移动的Ｍｅｔｒｏｐｏｌｉｓ算法：误差项均匀分布）　参数值　（ｂ）去除燃烧期后的后验直方图（误差项均匀分布）　图３随机游走（误差项均匀分布）Ｍ—Ｈ算法的模拟结果（先验分布为ｔ分布）　Ｆｉ窖．３　Ｓｉｍｕｌａｔｉｏｎ　ｏｆ　Ｍ—Ｈ　ａｌｇｏｒｉｔｈｍ　ｆｏｒ　ｒａｎｄｅｎｌ　ｗａｌｋ　ｓａｍｐｌｉｎｇ（ｅｒｒｏｒ　ｔｅｒｍ：ｕｎｉｆｏｒｍ　ｄｉｓｔｒｉｂｕｔｉｏｎ）　３结论　１）通过可读性强的Ｍａｔｌａｂ程序来实现了大家公认的、难度比较高的Ｍ—Ｈ（Ｍｅｓｔｒｏｐｌｉｓ　Ｈａｓｔｉｎｇｓ）算法。在　编制的源程序中，通过明确引入实参数和数据，计算什么，还需要什么，都一目了然。源程序包括了Ｍ—Ｈ法　中的独立抽样（独立链）和随机游走抽样算法。　２）在独立抽样Ｍ—Ｈ算法中，根据建议分布得到的前后两次样本值是相互独立的，因此马尔可夫链的　计算效率相对较低；当建议分布与后验分布接近的时候，马尔可夫链收敛速度快，可以得到很准确的计算结　果；但是实际应用中，后验分布的形式往往并不明确，若建议分布与后验分布差别很大时，独立抽样算法的　表现较差．甚至无法收敛。　３）在随机游走Ｍ—Ｈ算法中，根据建议分布得到的当前样本值与前一次样本值是相关的，克服了独立　抽样Ｍ—Ｈ算法计算效率较低的缺陷。文中分别给出了随机游走误差项服从正态分布和均匀分布作为最终　建议分布的样Ｍ—Ｈ算法的源程序，并给出了相应算例，均得到了计算效率高、收敛速度快的马尔可夫链和　精确的计算结果。该算法可用于贝叶斯推断过程中后验分布有日寸元　明确的情况。　参考文献：　［１］ＥＥＲＩＧＨＴ　Ｍ　Ｐ．Ｓｅｒｖｉｃｅ－ｌｉｆｅ　ｐｒｅｄｉｃｔｉｏｎ　ｏｆ　ｄｅｔｅｒｉｏｒａｔｉｎｇ　ｃｏｎｃｒｅｔｅ　ｂｒｉｄｇｅｓ［Ｊ］．Ｊｏｕｒｎａｌ　ｏｆ　Ｓｔｒｕｃｔｕｒａｌ　Ｅｎｇｉｎｅｅｒｉｎｇ，１９９８，１２４（３）：３０９—　３ｌ７．　【２］ＣＨＥＵＮＧ　Ｓ　Ｈ，ＢＥＣＫ　Ｊ　Ｌ．Ｂａｙｅｓｉａｎ　ｍｏｄｅｌ　ｕｐｄａｔｉｎｇ　ｕｓｉｎｇ　ｈｙｂｒｉｄ　ｍｏｎｔｅ　ｃａｒｌｏ　ｓｉｍｕｌａｔｉｏｎ　ｗｉｔｈ　ａｐｐｌｉｃａｔｉｏｎ　ｔｏ　ｓｔｕｃｔｒｕｒａｌ　ｄｙｎａｍｉｃ　ｍｏｄｅｌｓ　ｗｉｔｈ　ｍａｎｙ　ｕｎｃｅｒｔａｉｎ　ｐａｒｍｅｔａｅｒｓ［］Ｊ．Ｊｏｕｒｎａｌ　ｏｆＥｎｇｉｎｅｅｒｉｎｇ　Ｍｅｃｈａｎｉｃｓ，２００９，１３５（４）：２４３－２５５．　【３】ＷＡＮ　ＨＵＡＰＩＮＧ，ＲＥＮ　ＷＥＩＸＩＮ．Ｐａｒａｍｅｔｅｒ　ｓｅｌｅｃｔｉｏｎ　ｉｎ　ｉｆｎｔｅ　ｅｌｅｍｅｎｔ　ｍｏｄｅｌ　ｕｐｄａｔｉｎｇ　ｂｙ　ｇｌｏｂａｌ　ｓｅｎｓｉｔｉｖｉｔｙ　ａｎａｌｙｓｉｓ　ｕｓｉｎｇ　ｇａｕｓｓｉａｎ　ｐｒｏｃｅｓｓ　ｍｅｔａｍｏｄｅ［Ｊ］．Ｊｏｕｒｎａｌ　ｏｆ　Ｓｔｒｕｃｔｕｒｌａ　Ｅｎｇｉｎｅｅｒｉｎｇ，２０１５，１４１（６）：１－１０．　

８　华东交通大学学报　２０１８正　［４］李硕豪，张军．贝叶斯网络结构学习综述【Ｊ１．计算机应用研究，２０１５，３２（３）：６４１—６４６．　【５　ＣＨＩ５ＪＢ　Ｓ．Ｃｈａｐｔｅｒ　５７　ｍａｒｋｏｖ　ｃｈａｉｎ　ｍｏｎｔｅ　ｃａｒｌｏ　ｍｅｔｈｏｄｓ：ｃｏｍｐｕｔａｔｉｏｎ　ａｎｄ　ｉｎｆｅｒｅｎｃｅ［Ｊ］．Ｈａｎｄｂｏｏｋ　ｏｆ　Ｅｃｏｎｏｍｅｔｒｉｃｓ，２００１，５（１）：　３５６９—３６４９．　【６］伊庭幸人，獯村正美，大森裕浩，ｅｔ　ａ１．妍算统趼ｌＩ７，Ｌ，３７连锁乇、／于力，Ｌ／口法匕奄　周边【Ｍ】．东京：岩波害店，２００５．　Ｌ　王双成，高瑞，杜瑞杰．具有超父结点时间序列贝叶斯网络集成回归模型［Ｊ１．计算机学报，２０１７（４０）：１－１６．　【８】ＰＣ　ＬＡＭＢＥＲＴ，ＡＪ　ＳＵＴＦＯＮ，ＰＲ　ＢＵＲＴＯＮ　ｅｔ　ａ１．Ｈｏｗ　ｖａｇｕｅ　ｉｓ　ｖａｇｕｅ？Ａ　ｓｉｍｕｌａｔｉｏｎ　ｓｔｕｄｙ　ｏｆ　ｔｈｅ　ｉｍｐａｃｔ　ｏｆ　ｔｈｅ　ｕｓｅ　ｏｆ　ｖａｇｕｅ　ｐｒｉｏｒ　ｄｉｓｔｉｒｂｕｔｉｏｎｓ　ｉｎ　ＭＣＭＣ　ｕｓｉｎｇ　ｗｉｎＢＵＧＳ［Ｊ］．Ｓｔａｔｉｓｔｉｃｓ　ｉｎ　Ｍｅｄｉｃｉｎｅ，２００５，２４（１５）：２４０１—２４２８．　［９】余芳．基于ＭＣＭＣ方法的ＷｉｎＢＵＧＳ软件应用［Ｊ】．经营管理者，２０１０（１５）：３０３—３０３．　［１０】ＳＵＲＨＯＮＥ　Ｌ　Ｍ，ＴＥＮＮＯＥ　Ｍ　Ｔ，ＨＥＮＳＳＯＮＯＷ　Ｓ　Ｆ．ＯｐｅｎＢＵＧＳ［Ｍ］．Ｂｅｔａｓｃｒｉｐｔ　Ｐｕｂｌｉｓｈｉｎｇ，２０１０．　【ｌｌ】ＫＵＭＡＲ　Ｒ，ＫＵＭＡＲ　ＳＲＩＶＡＳＴＡＶＡ　Ａ，ＫＵＭＡＲ　Ｖ．Ｅｘｐｏｎｅｎｔｉａｔｅｄ　ｇｕｍｂｅｌ　ｍｏｄｅｌ　ｆｏｒ　ｓｏｆｔｗａｒｅ　ｒｅｌｉａｂｉｌｉｔｙ　ｄａｔａ　ａｎａｌｙ　ｓｉｓ　ｕｓｉｎｇ　ＭＣＭＣ　ｓｉｍｕｌａｔｉｏｎ　ｍｅｈｏｄ［ｔＪ］．Ｉｎｔｅｒｎａｔｉｏｎａｌ　Ｊｏｕｒｎａｌ　ｏｆＣｏｍｐｕｔｅｒＡｐｐｌｉｃａｔｉｏｎｓ，２０１３，６２（２０）：２４—３２．　［１２］韦来生，张伟平．贝叶斯分析［Ｍ］．合肥：中国科学技术大学出版社，２０１３．　Ｂａｙｅｓｉａｎ　Ｓｔａｔｉｓｔｉｃｓ　ａｎｄ　ＭＣＭＣ　Ｍｅｔｈｏｄ　Ｍａｔｌａｂ　Ｐｒｏｇｒａｍｍｉｎｇ　ｆｏｒ　Ｍｅｔｒｏｐｏｌｉｓ－Ｈａｓｔｉｎｇｓ（Ｍ－Ｈ）Ａｌｇｏｒｉｔｈｍ　Ｃｈｅｎ　Ｍｅｎｇｃｈｅｎｇ，Ｆａｎｇ　Ｗｅｉ，Ｙａｎｇ　Ｃｈａｏ，Ｘｉｅ　Ｌｉ　（Ｓｃｈｏｏｌ　ｏｆ　Ｃｉｖｉｌ　Ｅｎｇｉｎｅｅｒｉｎｇ　ｎｄ　ａＡｒｃｈｉｔｅｃｔｕｒｅ，Ｅａｓｔ　Ｃｈｉｎａ　Ｊｉａｏｔｏｎｇ　Ｕｎｉｖｅｒｓｉｔｙ，Ｎａｎｃｈａｎｇ　３３００１３，Ｃｈｉｎａ）　Ａｂｓｔｒａｃｔ：Ｂａｙｅｓ　ｓｔａｔｉｓｔｉｃｓ　ｉｎ　ａ　ｖａｃｕｕｍ　ｃａｎ　ｂｅ　ｂｒｏａｄｌｙ　ａｐｐｌｉｅｄ　ｔｏ　ｔｈｅ　ｉｅｌｆｄｓ　ｏｆ　ｎａｔｕｒａｌ　ｓｃｉｅｎｃｅ，ｅｃｏｎｏｍｉｃｓ　ａｎｄ　ＳＯ－　ｃｉｌ　ｓｃｉａｅｎｃｅ，ｗｈｉｃｈ　ｂｅｎｅｆｉｔｓ　ｆｒｏｍ　ｔｈｅ　ｄｅｖｅｌｏｐｍｅｎｔ　ｏｆ　ｃｏｍｐｕｔｅｒ　ｓｃｉｅｎｃｅ　ａｎｄ　ｔｅｃｈｎｏｌｏｇｙ　ａｎｄ　Ｍａｒｋｏｖ　ｃｈａｉｎ　Ｍｏｎｔｅ　Ｃａｒｌｏ（ＭＣＭＣ）ｍｅｔｈｏｄ．Ｉｎ　ｔｈｉｓ　ｐａｐｅｒ　ｔｈｅ　ａｐｐｌｉｃａｔｉｏｎ　ｏｆ　ＭＣＭＣ　ｍｅｔｈｏｄ　ｉｎｔｏ　Ｂａｙｅｓ　ｉｎｆｅｒｅｎｃｅ　ｗａｓ　ｉｎｔｒｏｄｕｃｅｄ，ａｎｄ　ｉｎｄｅｐｅｎｄｅｎｔ　ｓａｍｐｌｉｎｇ　ａｎｄ　ｒａｎｄｏｍ　ｗａｌｋ　ｓａｍｐｌｉｎｇ　ｏｆ　Ｍｅｔｒｏｐｏｌｉｓ－Ｈａｓｔｉｎｇｓ（Ｍ－Ｈ）ａｌｇｏｒｉｔｈｍ　ｗｅｒｅ　ｍａｉｎｌｙ　ｄｉｓｃｕｓｓｅｄ．　Ｔｈｅ　ｔｗｏ　ｓａｍｐｌｉｎｇ　ｍｏｄｅｓ　ｗｅｒｅ　ｒｅａｄａｂｌｙ　ｐｒｏｇｒａｍｍｅｄ　ｗｉｔｈ　Ｍａｔｌａｂ，ａｎｄ　ｔｈｅｉｒ　ｄｅｔａｉｌｅｄ　ｉｍｐｌｅｍｅｎｔａｔｉｏｎ　ｐｒｏｃｅｓｓｅｓ，　ｍｅｒｉｔｓ　ａｎｄ　ｄｅｍｅｒｉｔｓ　ｗｅｒｅ　ｔａｌｋｅｄ　ａｂｏｕｔ．Ｉｔ　ｗａｓ　ｓｈｏｗｎ　ｂｙ　ｐｒｅｓｅｎｔ　ｓｉｍｕｌａｔｉｏｎ　ｔｈａｔ，ｔｈｅ　ｉｎｄｅｐｅｎｄｅｎｔ　ｓａｍｐｌｉｎｇ　ｍｏｄｅ　ｉｓ　ｒｅｌａｔｉｖｅｌｙ　ｅａｓｙ　ｔｏ　ｉｍｐｌｅｍｅｎｔ，ｂｕｔ　ｎｅｅｄ　ｔｈｅ　ｐｒｏｐｏｓａｌ　ｄｉｓｔｒｉｂｕｔｉｏｎ　ｔｏ　ｂｅ　Ｃｌｏｓｅ　ｔｏ　ｔｈｅ　ｐｏｓｔｅｒｉｏｒ　ｄｉｓｔｉｂｕｔｒｉｏｎ；ｏｔｈｅｒ．　ｗｉｓｅ，ｔｈｅ　ｃａｌｃｕｌａｔｉｏｎ　ｅｆｉｃｉｆｅｎｃｙ　ｉｓ　ｌｏｗ　ａｎｄ　ｔｈｅ　ｓｉｍｕｌａｔｉｏｎ　ｅｆｆｅｃｔ　ｉｓ　ｕｎｓａｔｉｓｆａｃｔｏｒｙ．Ｔｈｅ　ｒａｎｄｏｍ　ｗａｌｋ　ｓａｍｐｌｉｎｇ　ｍｏｄｅ　ｄｏｎ’ｔ　ｎｅｅｄ　ｔｈｅ　ｐｒｏｐｏｓａｌ　ｄｉｓｔｉｒｂｕｔｉｏｎ　ｔｏ　ａｐｐｒｏａｃｈ　ｔｈｅ　ｐｏｓｔｅｒｉｏｒ　ｄｉｓｔｉｂｕｔｒｉｏｎ　ａｎｄ　ｉｔｓ　ｓｉｍｕｌａｔｉｏｎ　ｒｅｓｕｌｔｓ　ａｒｅ　ｓａｔｉｓｆａｃｔｏｒｙ．Ｔｈｅｒｅｆｏｒｅ，ｉｔ　ｏｖｅｒｃｏｍｅｓ　ｔｈｅ　ｌｉｍｉｔａｔｉｏｎｓ　ｏｆ　ｔｈｅ　ｉｎｄｅｐｅｎｄｅｎｔ　ｓａｍｐｌｉｎｇ　ｍｏｄｅ　ｈａｓ　ｍｏｒｅ　ｗｉｄｅｓｐｒｅａｄ　ａｐ—　ｐｌｉｃａｔｉｏｎ．　Ｋｅｙ　ｗｏｒｄｓ：Ｂａｙｅｓ；ＭＣＭＣ　ｍｅｔｈｏｄ；Ｍ－Ｈ　ａｌｇｏｒｉｔｈｍ；Ｍａｔｌａｂ　ｐｒｏｇｒａｍ