2008年1月 系统士程理论与实践 第1期 文章编号：1000-6788(2008)01-0100-09 Metropolis-Hastings自适应算法及其应用 陈平'，徐若曦2 (1.东南大学数学系，南京210096：2.俄亥俄州立大学统计系，美国43210-1247) 摘要：首先阐述Metropolis-Hastings算法实现的具体步骤，然后证明由此产生的Markov链满足细致平衡· 条件，从而以目标分布为不变分布接下来给出几个计算实例，以说明提议函数及其方差的选取对采样 结果的影响，并由此推出一种改进的自适应算法用以寻找合适的提议函数及其方差.最后，通过贝叶斯 Lgtc模型的例子说明M-H方法在贝叶斯分析中的应用，同时也检验M-H自适应算法的效果. 关簧词：MCMC;Metropolis-Hastings算法；马尔可夫链：R软件：贝叶斯分析 中图分类号：F830;TP183 文献标志码：A Metropolis-hastings adaptive algorithm and its application CHEN Ping',XU Ruo-xi2 (1.Department of Mathematics.Southeast Univerity,Nanjing 210096,China 2.Department of Statistics,The Ohio Suate University, Columbus,OH 43210-1247) Abstract:Markoy chain Monte Carlo(MCMC)methods is an important class of computer based simulation techniques.This paper investigates one MCMC method known as the Metropolis-Hastings algorithm.In this paper,we first introduce readers the proceedings of the Metropolis-Hastings algorithm.Then we prove the resulting chain satisfies detailed balance,and hence has the target distributionas the invariant distribution.Next,e providellustrative examples that show the influence of the propoeal function and its variance on the resulting chain,and developan adaptive method to find optimal propo for the random walk sampler.Finally,we discuss the relationship between M- H algorithm and Bayeian analysis.The Bayesian Logistic model is used to ilurative the application of M-H algorithm in Bayesian anaysis and to test the proposed adaptive method. Key words:MCMC:Metropolis-Hastings algorithm:markov chain:R software;bayesian analysis 1引言 MCMC(Markov chain Monte Carlo)方法又被称为动态Monte Carlo方法，它是以动态构造Markov链为基 础，通过遍历性约束来实现模拟目标分布的一类随机模拟方法.在过去的20年里，MCMC方法对统计学， 尤其是贝叶斯分析的发展产生了深远的影响.本文借助于R软件对MCMC方法中的Metropolis-Hastings(M- H)算法作进一步的探讨，并由此推出一种改进的自适应MH算法，从而可以提高算法的效率.同时还将借 助于贝叶斯分析的实例阐明MH算法与现代统计推断间的紧密联系， 假设我们要计算函数F(x)关于概率密度π(x)的平均值，即 E[F()]=F()x(s)dx 如果这个积分用分析的方法无法解决，那么随机模拟的方法就可以用来估计积分值.前提是我们能够从π (x)中产生独立同分布的随机数如果(X,,X)i.i.d.一π(x),那么依据大数定律，∑F(X)依概率 n 收敛到F(x)π(x)dx.于是问题转化为如何从概率分布x(x)中产生独立同分布的随机数.虽然在很多 收稿日期：2006-07-13 资助项目：国家自然科学基金(10671032) 作者葡介：陈平(1960-)，男，江苏溧阳人，教授，博士，研究方向：金融与工程时间序列分析、生存分析等，E-mail:pl8@ 263.net.cn. 万方数据
2008年1月 系统土程理论与实践 第1期 文章编号：1000．6788(2008)01-0100—09 Metropolis．Hastings自适应算法及其应用 陈 平1，徐若曦2 (1．东南大学数学系，南京210096；2．俄亥俄州立大学统计系，美国；43210-1247) 摘要：首先阐述Metropolis-Hastings算法实现的具体步骤，然后证明由此产生的Markov链满足细致平衡 条件，从而以目标分布为不变分布．接下来给出几个计算实例，以说明提议函数及其方差的选取对采样 结果的影响，并由此推出一种改进的自适应算法用以寻找合适的提议函数及其方差．最后，通过贝叶斯 Logistic模型的例子说明M．H方法在贝叶斯分析中的应用，同时也检验M．H自适应算法的效果． 关键词：MCMC；Metropofis-Hastings算法；马尔可夫链；R软件；贝叶斯分析 中图分类号：F830；TPl83 · 文献标志码：A Metropolis—hastings adaptive algorithm and its application CHEN Pin91，XU Ruo．xi2 (1．Department of Mathematics，Southeast University，Nanjing 210096，China；2．Department of Statistics，The Ohio State University， Columbus，OH 43210．1247) ’ Abstract： Markov chain Monte Carlo(MCMC)methods is an important class of computer based simulation techniques．This paper investigates one MCMC method known聃the Metropolis-Hastings algorithm．In this paper．we first introduce readel葛the proceedings of the Metropolis-Hastings algorithra．Then we prove the resulting chain satisfies detailed balance，and hence has the target distribution a8 the invariant distribution．Next，we provide some illustrative examples that show the influence of the proposal function and its variance on the陀吼millg chain，and develop an adaptive method to find optimal p∞p088l for the random walk sampler．Finally，we discuss the relationship between MH algorithm and Bayesian analysis．The Bayesian logistic model is used to illustrative the application of M-H algorithm in Bayesian analysis and to test the proposed adaptive method． Key words： MCMC；Metropolis-Hastings algorithm；mⅡkov chain；R software；bayesian analysis 1 引言 MCMC(Markov chain Monte Carlo)方法又被称为动态Monte Carlo方法，它是以动态构造Markov链为基 础，通过遍历性约束来实现模拟目标分布的一类随机模拟方法．在过去的20年里，MCMC方法对统计学， 尤其是贝叶斯分析的发展产生了深远的影响．本文借助于R软件对MCMC方法中的Metropolis．Hastings(M． H)算法作进一步的探讨，并由此推出一种改进的自适应M．H算法，从而可以提高算法的效率．同时还将借 助于贝叶斯分析的实例阐明M．H算法与现代统计推断间的紧密联系． 假设我们要计算函数F(并)关于概率密度玎(算)的平均值，即 E[F(茗)] =．r。，(茗)，r(髫)d省 如果这个积分用分析的方法无法解决，那么随机模拟的方法就可以用来估计积分值．前提是我们能够从丌 (髫)中产生独立同分布的随机数．如果(xl’．一，x。)i．i．d．一，r(茹)，那么依据大数定律，丢萃n，(置)依概率 收敛到．r。F(聋)，r(髫)d髫．于是问题转化为如何从概率分布，r(髫)中产生独立同分布的随机数．虽然在很多 收稿日期：2006．07．13 资助项目：国家自然科学基金(10671032) 作者简介：陈平(1960一)，男，江苏溧阳人，教授，博士，研究方向：金融与工程时间序列分析、生存分析等，E-mail：cpl8@ 263．net．c11． 万方数据
第1期 Metropoli:Hasting粉自适应算法及其应用 101 问题中，直接产生独立同分布的随机数并不容易，但是很多MCMC方法，如Gibbe采样法，Metropolis采样 法，却能够通过随机模拟产生收敛到x(x)的Markov链.例如，通过100次模拟产生100条长度为5000的 独立的Markov链，那么每条链的最后一个值X0就可以近似地看作来自极限分布π(x)的大小为100的 独立同分布的样本，当然，在实际计算中，更高效的做法是在Markov链运行了一段时间（比如500次）后， 间隔k个点采样一次，得到的样本仍可近似看作来自π(x)的独立同分布的样本.本文主要讨论MCMC方 法中很重要的一类算法：M-H算法.我们将分析M-H算法构造Markov链的方法，证明其具有不变分布（极 限分布)，探讨影响算法效率的因素，以及算法在贝叶斯分析中的一些应用. 具体安排是：第二节阐述Metropolis-Hastings算法的实现过程及一些常用的形式，给出算法可逆性的证 明；第三节给出一些计算实例，提出MH自适应算法及改进一般算法的途径并通过一些例子做比较分析. 第四节给出贝叶斯分析的有关理论，通过一个具体的贝叶斯模型例子来阐明M-H自适应算法在贝叶斯分 析中的应用，同时再次检验第三节中提出的算法的效果.最后给出简要结论， 2 Metropolis-Hastings算法及其性质 M-H算法的思想是构造一个以目标分布π(x)为不变分布的Markov链.为了实现这一目标，M-H算法 借助于一个辅助的概率密度函数9(x,y),通常被称为“提议函数”，q(x,y)通常要满足以下三个条件（其 中S表示状态空间)：(a)对于固定的x,9(x,·)是一个概率密度函数.(b)对于Vx,y∈S,g(x,y)的值要 能够计算出来.(c)对于固定的x,能够方便地从q(x,y)中产生随机数. 在满足上述三个条件的情况下，9(x,y)的选取是任意的.下面是算法的具体步骤，假设当前状态X, =,那么 1)从提议函数q(x,·)中产生一个新状态y. 2)计算接受概率 a,)m282} (1) 3)以概率a(x,y)置X.+1=y:以概率1-a(x,y)置X1=x, 我们将证明此算法构造的Markov链以目标分布π为不变分布. 从理论上讲，提议函数g(x,y)的选取是任意的，但在实际计算中，提议函数的选取对于算法的效率 的影响是相当大的.一般认为提议函数的形式与目标分布越接近，则模拟的效果越好，我们将在第三节里 仔细讨论提议函数的选取问题.在算法的第3步中，如果置X.+1=y,我们就称之为“接受提议”；否则就称 为“拒绝提议”.在实际计算中以概率a(x,y)接受y可以这样实现：从U(0,1)产生一个随机数“，如果4 0,那么在后续的计算中分母为 零的概率将是零，从而不会对实际计算造成麻烦. 如果M-H算法中的提议函数q(x,y)不仅满足对称性，而且只与￥-y有关，那么算法就演变为通常 意义下的随机游动采样法.最常见的一种随机游动采样法以正态分布，即9(x,y)=(x-y)为提议函数. 这里中是任意一种多维正态分布的密度函数，也就是说y~(x,Σ).（Σ是任意的正定矩阵，x是当前的 状态).在实际计算中，由于Σ要求是正定矩阵，所以常取为Σ=，·是一个参数.这时一个重要的问题是 ·应取多少才合适.如果取的过大，那么大多数提出的新状态将落人分布的尾部，从而被拒绝.如果取的过 小，那么采样的自相关系数就会很高，而且收敛到目标分布的速度也会很慢.不论哪种情况，我们都不能在 时间有限的采样过程中得到较好地服从目标分布的随机数.一般认为在模拟多维正态随机数时，如果调整 σ的值使得在整个模拟过程中提议被接受的比例大约为20%，将得到比较好的模拟效果. 在时齐且有条件转移密度的情形，我们简记 P(x,y)p()=p(x.) 本文所考虑的Markov链都是时齐的，并总假定条件转移密度存在.M-H算法构造的Markov链的条件转移 万方数据
第1期 Metropolis-Hastings白适应算法及其应用 101 问题中，直接产生独立同分布的随机数并不容易，但是很多MCMC方法，如Gibbs采样法，Metropolis采样 法，却能够通过随机模拟产生收敛到，r(石)的Markov链．例如，通过100次模拟产生100条长度为5000的 独立的Markov链，那么每条链的最后一个值x舯就可以近似地看作来自极限分布，r(菇)的大小为100的 独立同分布的样本．当然，在实际计算中，更高效的做法是在Markov链运行了一段时间(比如500次)后， 间隔k个点采样一次，得到的样本仍可近似看作来自，r(髫)的独立同分布的样本．本文主要讨论MCMC方 法中很重要的一类算法：M-H算法．我们将分析M．H算法构造Markov链的方法，证明其具有不变分布(极 限分布)，探讨影响算法效率的因素，以及算法在贝叶斯分析中的一些应用． 具体安排是：第二节阐述Metropolis．Hastings算法的实现过程及一些常用的形式，给出算法可逆性的证 明；第三节给出一些计算实例，提出M．H自适应算法及改进一般算法的途径并通过一些例子做比较分析． 第四节给出贝叶斯分析的有关理论，通过一个具体的贝叶斯模型例子来阐明M—H自适应算法在贝叶斯分 析中的应用，同时再次检验第三节中提出的算法的效果．最后给出简要结论． 2 Metropolis．Hastings算法及其性质 M．H算法的思想是构造一个以目标分布丌(戈)为不变分布的Markov链．为了实现这一目标，M．H算法 借助于一个辅助的概率密度函数q(髫，Y)，通常被称为“提议函数”．q(茗，，，)通常要满足以下三个条件(其 中S表示状态空间)：(a)对于固定的茗，q(髫，·)是一个概率密度函数．(b)对于V茗，Y∈S，q(茗，Y)的值要 能够计算出来．(c)对于固定的龙，能够方便地从口(石，y)中产生随机数． 在满足上述三个条件的情况下，q(善，Y)的选取是任意的．下面是算法的具体步骤，假设当前状态瓦 =戈，那么 1)从提议函数q(髫，·)中产生一个新状态Y． ． 2)计算接受概率 · 出’y)“n{-，揣) 3)以概率a(菇，y)置以+l-Y；以概率1一口(茗，，，)置置+I=茗． 我们将证明此算法构造的Markov链以目标分布，r为不变分布． 从理论上讲，提议函数q(茗，Y)的选取是任意的，但在实际计算中，提议函数的选取对于算法的效率 的影响是相当大的．一般认为提议函数的形式与目标分布越接近，则模拟的效果越好．我们将在第三节里 仔细讨论提议函数的选取问题．在算法的第3步中，如果置以+。=Y，我们就称之为“接受提议”；否则就称 为“拒绝提议”．在实际计算中以概率口(菇，y)接受Y可以这样实现：从v(0，1)产生一个随机数“，如果扯 0，那么在后续的计算中分母为 零的概率将是零，从而不会对实际计算造成麻烦． 如果M．H算法中的提议函数g(石，，，)不仅满足对称性，而且只与菇一Y有关，那么算法就演变为通常 意义下的随机游动采样法．最常见的一种随机游动采样法以正态分布，即口(茗，，，)=声(菇一，，)为提议函数． 这里声是任意一种多维正态分布的密度函数，也就是说Y一Ⅳ(茗，三)．(三是任意的正定矩阵，菇是当前的 状态)．在实际计算中，由于Z要求是正定矩阵，所以常取为三=al，口是一个参数．这时一个重要的问题是 口应取多少才合适．如果取的过大，那么大多数提出的新状态将落入分布的尾部，从而被拒绝．如果取的过 小，那么采样的自相关系数就会很高，而且收敛到目标分布的速度也会很慢．不论哪种情况，我们都不能在 时间有限的采样过程中得到较好地服从目标分布的随机数．一般认为在模拟多维正态随机数时，如果调整 盯的值使得在整个模拟过程中提议被接受的比例大约为20％，将得到比较好的模拟效果． 在时齐且有条件转移密度的情形，我们简记 P(茗，Y)暑P‘“。+u(茹，Y)=PW4’(茗，Y) 本文所考虑的Markov链都是时齐的，并总假定条件转移密度存在．M．H算法构造的Markov链的条件转移 万方数据
102 系统工程理论与实践 2008年1月 密度为： [q(x.r)a(x,y), Vy≠常 p(x,y)= 1-9(,y)a(x,yay,y=￥ (2) 下面的定理说明M-H算法构造的Markov链以目标分布π(x)为不变分布，并且是时间可逆的. 定理1若以P为条件转移密度的马尔科夫链具有可逆初分布π，那么π也一定是其不变分布， 证明 Pr(X...=y)=x(=)P(x.y)dx =(y)P(y.x)dz =x(y)P(y,x)dx=x(y） y的边缘分布即为π，这就证明了这个定理. 有了这个定理，我们下面所要做的就是验证M-H算法构造的Markov链满足细致平衡条件， 定理2设P(x,y)是由(2)式定义的条件转移密度.则在状态空间S上，对于任意的初始分布，以 P为条件转移密度的Markov链满足细致平衡条件（时间可逆的），即 x(x)P(x,y)=π(y)P(y,x),x,y∈S 证明如果y=,那么自然π(x)P(x,x)与其自身相等，对于y≠x的情况， (1)若 π(y）q(y,x) (9￥，y≤1 则 ）=名82 并且a(y,x)=1 从而有 x到P(,=g(）=xg,2号 =x(y)q(y,a)=(y)g(y,)a(y,x)=(y)p(y,x) (i)若 π(y)q(y,x) x()9￥，≥1 则 -8 ，并且a(x,y）=1 从而有 (y）P(y,)=x(yg(y,)a(,)=x()g(y,)g2 π（y)g(y,x) =π(x)q(,y）=r(x)q(x,y)a(x,y)=r(x)P(x,y） 不论哪种情况，均有π(x)P(x,y)=π(y)P(y,x),证毕. 3两种Metropolis--Hastings算法比较 R软件是一种基于S语言而开发的用于数据处理，数值计算以及图象显示的免费软件.虽然R软件不 仅仅用于统计计算，但多数人将其当作一种统计软件来使用，原因是R软件的工作环境中集成了很多经 典或现代的统计计算方法.其中一部分包含在基本工作环境中，更多的是以数据包的形式提供给使用者， 数据包是R软件集成环境的体现，任何使用者都可以编写数据包，供自己或他人使用.R软件的基本工作 环境自带约25个标准数据包，更多其他的数据包可以从任一CRAN(Comprehensive R Archive Network)站点 及其镜象站点获得. 万方数据
102 系统工程理论与实践 2008年1月 密度为： ’ 。 fq(茹，Y)a(石，Y)， V，，≠茗 P‘茗，，，’2 il—f口(菇，，，)口(菇，y)d，，， ，，：茗 ‘2’ 下面的定理说明M．H算法构造的MaIkov链以目标分布丌(茹)为不变分布，并且是时间可逆的． 定理1若以P为条件转移密度的马尔科夫链具有可逆初分布丌，那么，r也一定是其不变分布． 证明 Pr(Xm+l=Y)=I，r(茗)p(x，y)dx=I，r(，，)P(Y，x)dx ． 一 J =，r(，，)J．P(y，茹)d髫=，r(，，) y的边缘分布即为7f，这就证明了这个定理． 有了这个定理，我们下面所要做的就是验证M—H算法构造的Markov链满足细致平衡条件． 定理2‘设P(茗，，，)是由(2)式定义的条件转移密度．则在状态空间S上，对于任意的初始分布，r，以 P为条件转移密度的Markov链满足细致平衡条件(时间可逆的)，即 ，r(茗)P(茹，Y)=，r(，，)P(Y，茗)，V茹，y∈S 证明 如果Y=茹，那么自然，r(髫)P(髫，膏)与其自身相等．对于y≠菇的情况， (i)若 毫语q舞Y引 ， ，’ ，r(茗)(菇，)每1 则 出∽=老毁乌并且山㈠～ 从而有 · ，r(髫)尸(茁，)，)=丌(菇)口(菇，，，)口(露，，，)=，r(菇)口(茗，，，)罢{渊 =，r(，，)q(Y，茹)=，r(y)q(，，，髫)口(Y，戈)=7f(，，)P(y，髫) (ii)若 ．老{捌q≥t ，f(菇)(髫，，，)乎1 则 山，菇)=老燎高并且出∽～ 从而有 丌(，，)尸(，，，戈)=丌(y)g(，，，菇)口(，，，茹)=丌(，，)g(，，，茹)罢渊 =，r(菇)q(x，，Y)=，r(茗)q(x，，，)a(茁，Y)=7f(聋)P(x，Y) 不论哪种情况，均有，r(菇)P(菇，’，)=丌(Y)P(Y，髫)，证毕． 3 两种Metropolis．Hastings算法比较 R软件是一种基于S语言而开发的用于数据处理，数值计算以及图象显示的免费软件．虽然R软件不 仅仅用于统计计算，但多数人将其当作一种统计软件来使用，原因是R软件的工作环境中集成了很多经 典或现代的统计计算方法．其中一部分包含在基本工作环境中，更多的是以数据包的形式提供给使用者． 数据包是R软件集成环境的体现，任何使用者都可以编写数据包，供自己或他人使用．R软件的基本工作 环境自带约25个标准数据包，更多其他的数据包可以从任一CRAN(Comprehensive R Archive Network)站点 及其镜象站点获得． 万方数据
第1期 Metropolis-Hasting自适应算法及其应用 103 例1假设目标分布函数为N(3,5),提议函数采用随机游动，即新状态y~N(X,¤)，其中X.为当 前状态，a为标准差，取为2.M-H算法的步骤如下： 1)取初值X1=100. 2)从提议函数N(X。,)中产生一个新状态y, 3)计算接受概率 a(s)=m1,3}=mm1,e--2-67+6} 其中r(x)为正态分布N(3,5)的密度函数. 4)以概率a(x,y)置X.1=y;以概率1-a(x,y)置X1=x. 我们构造了一条长度为10000的Markov链，为了消除初值对采样的影响，从第500个样本开始取样分 析，从直方图上看采样结果令人满意，估计的密度函数与目标分布函数十分接近.但为了得到近似独立的 样本，往往采用间隔采样法.从样本的自相关系数图（略）中可看出在lg=40左右时，样本的自相关系数已 接近于零，故若我们每隔40个点取样一次，则取出的样本可近似看作独立的样本，但这样一来，我们从这 条链中最多只能取得250个近似独立的样本.因此，在实际工作中，总是希望1g能取的越小越好. 进一步用不同的g的值(。=1,2,3,10,15,40,100,200)重复上面的实验，得到的结果显示：当取c=10 或15时，用同样的算法得到的Markov链可以用不到20的间隔进行采样，从而所能得到的近似独立的样本 大小将两倍于前者，效率大大提高.这就是说。的选取确实对算法有很大的影响，但是。的选取和哪些因 素有关呢？是否存在普遍适用的值呢？为此，我们再用一个例子来说明这一点。 例2我们将目标分布换为Bta(3,5),提议函数仍然采用随机游动，即新状态y~N(X.,o),其中X。 为当前状态.同样，我们对不同的。值(a=0.05,0.1,0.4,0.6,1,3,10,15)进行重复实验，得到的结果如图 1所示， 从图1中易见较为合适的σ值为0.4,0.6,1.若如上例一样取10或15作为g的值，则无法得到较好 的采样结果、例1与例2说明。的选取要考虑目标分布.对应于不同的目标分布，合适的。值可能相差甚 远，在实际计算中必须具体问题具体分析， 随机游动采样法是MH算法中较为常用的一种形式，我们已经看到对于随机游动采样法而言，如何 选择提议函数的方差成为影响算法效率的主要因素，对于一维情况，就是要选择一个合适的。值，而对于 d>1的多维情况，我们需要一个合适的正定方差矩阵Σ.那么现在的问题是什么样的▣值或者方差矩阵 Σ才是合适的呢？一种想法是调整接受概率到合适的值.根据关于采用多维正态作为提议函数的随机游 动采样法的一些理论结果，一般认为如果目标分布是d维的且各分量之间独立，则调整接受概率到0.234 左右将最大限度地优化随机游动采样法的效率，并且接受概率0.234与目标分布的具体形式无关， 基于这样的想法，我们尝试构造一种自适应算法来寻找合适的提议函数的方差：大致的思想是让接受 比率落人一个可以接受的范围内，比如区间[a-e,a+e].以一维的情况为例，先给定一个初值o,用这 个值产生一条链，来估计接受比率P。.若P。在给定的区间中，则已符合要求，算法停止：若不然，则以。 为中心做一次随机游动，得到一个·，再计算P,是否满足要求，以此类推，直至找到合适的。值，具体的 实现步骤如下： 1)给定g一个初值o,置n=0. 2)以N(X-1,o,)为提议函数，用随机游动采样法产生一条长度为m的Markov链. 3)计算第2步中产生的Markov链接受新状态的比率P. 4)如果Pn∈[a-e,a+e],则置c=o.,退出；如果P.e[a-e,a+e],转人第5步. 5)如果n>2并且1P.-al<|P.-1-al,则置o.=o.-1· 6)从N(a。,G)中产生一个随机数，记为c+1,置n=n+1,转入第2步. 注第5步的判断是为了保证接受比率P.要优于P1,即。至少不会比0。-要坏；第6步中的的 选取没有固定的规律可循，主要是看对初始。的估计是否接近于合适的：值.若估计。与。较为接近， 万方数据
第1期 Metropolis-Hastings自适应算法及其应用 103 例1假设目标分布函数为N(3，5)，提议函数采用随机游动，即新状态Y．N(Xm，盯)，其中以为当 前状态，口为标准差，取为2．M．H算法的步骤如下： ‘1)取初值X。=100． 2)从提议函数Ⅳ(以，口)中产生一个新状态Y． 3)计算接受概率 出∽=min{·，器)=min{1’exp(一嘉(，，2一X2—6y+6x))) 其中，r(石)为正态分布N(3，5)的密度函数． 4)以概率口(菇，y)置兄+。=Y；以概率1一口(菇，，，)置X。+。=茹． 我们构造了一条长度为10000的Markov链，为了消除初值对采样的影响，从第500个样本开始取样分 析．从直方图上看采样结果令人满意，估计的密度函数与目标分布函数十分接近．但为了得到近似独立的 样本，往往采用间隔采样法．从样本的自相关系数图(略)中可看出在lag=40左右时，样本的自相关系数已 接近于零，故若我们每隔40个点取样一次，则取出的样本可近似看作独立的样本．但这样一来，我们从这 条链中最多只能取得250个近似独立的样本．因此，在实际工作中，总是希望lag能取的越小越好． 进一步用不同的d的值(盯=1，2，3，10，15，40，100，：ZOO)重复上面的实验，得到的结果显示：当取叮=10 或15时，用同样的算法得到的Markov链可以用不到20的间隔进行采样，从而所能得到的近似独立的样本 大小将两倍于前者，效率大大提高．这就是说盯的选取确实对算法有很大的影响，但是口的选取和哪些因 素有关呢?是否存在普遍适用的值呢?为此，我们再用一个例子来说明这一点． 例2我们将目标分布换为Beta(3，5)，提议函数仍然采用随机游动，即薪状态Y—N(Xn，盯)，其中以 为当前状态．同样，我们对不同的口值(口=0．05，0．1，0．4，0．6，1，3，10，15)进行重复实验，得到的结果如图 l所示． 从图1中易见较为合适的盯值为0．4，0．6，1．若如上例一样取10或15作为口的值，则无法得到较好 的采样结果、例1与例2说明盯的选取要考虑目标分布．对应于不同的目标分布，合适的口值可能相差甚 远，在实际计算中必须具体问题具体分析． 随机游动采样法是M．H算法中较为常用的一种形式，我们已经看到对于随机游动采样法而言，如何 选择提议函数的方差成为影响算法效率的主要因素．对于一维情况，就是要选择一个合适的口值，而对于 d>1的多维情况，我们需要一个合适的正定方差矩阵三．那么现在的问题是什么样的J值或者方差矩阵 三才是合适的呢?一种想法是调整接受概率到合适的值．根据关于采用多维正态作为提议函数的随机游 动采样法的一些理论结果．一般认为如果目标分布是d维的且各分量之间独立，则调整接受概率到O．234 左右将最大限度地优化随机游动采样法的效率，并且接受概率0．234与目标分布的具体形式无关． 基于这样的想法，我们尝试构造一种自适应算法来寻找合适的提议函数的方差：大致的思想是让接受 比率落入一个可以接受的范围内，比如区间[口一￡，a+e]．以一维的情况为例，先给定一个初值口。，用这 个值产生一条链，来估计接受比率P。．若Po在给定的区间中，则已符合要求，算法停止；若不然，则以口。 为中心做一次随机游动，得到一个口。，再计算P，是否满足要求，以此类推，直至找到合适的口值．具体的 实现步骤如下： 1)给定口一个初值％，置厅=0． 2)以Ⅳ(五一。，盯。)为提议函数，用随机游动采样法产生一条长度为m的Markov链． 3)计算第2步中产生的Markov链接受新状态的比率P。． 4)如果只∈[口一e，口+e]，则置口=口。，退出；如果P。薯[口一￡，口+￡]，转入第5步． 5)如果n>2并且I P。一口I<I PI—l一口l，则置盯。=仃。一1． 6)从N(口。，；)中产生一个随机数，记为口川，置n=n+1，转入第2步． 注第5步的判断是为了保证接受比率只要优于P。一。，即仃。至少不会比(Ylt-I要坏；第6步中的；的 选取没有固定的规律可循，主要是看对初始口。的估计是否接近于合适的口值．若估计口。与d较为接近， 万方数据
104 系统工程理论与实践 2008年1月 sigma-0.05 sigma-0.1 sigma-0.4 sigma-0.5 1.0 10f 1.0 1.0 08 0g、 08- 0.8 0.6 0.6- 0.6- 0.6 0.4- 总4 0.4- 0.2 02- 0.2 0.2- 0.0- 0.0- 0.0 0.0 0”0 III 020 60 02060 Lag Lag Lag Lag gm】 sigma-3 sigma=10 sigma=15 10 1.0 1.0 1.01 0.8 0.8- 0.8 0.8- 0.6 0.6- 0.6 0.6 0.4- 0.4- 04 0.2 02- 02- 02 0.0 0.0+ 0.0卡 00f.---- 十TTT 02060 020 60 02060 Lag Lag Lag 图1取。=0.05,0.1,0.4,0.6,1,3,10,15时，随机游动采样法得到的样本自相关系数图.横坐标为滞 后系数1，纵坐标为对应的自相关系数，两条虚线表示自相关系数为零的95%的置信区间， 则：可选择得小一些；若估计o。与。相差较远，则G可选择得大一些；对于d>1的多维情况，如果要确 定一个任意的正定矩阵三，则要确定d2个参量，显得有些麻烦更多情形下，我们取Σ=σl这样一来要确 定的参数就只有一个，同时上述算法也可以用于多维情形，只须在6步中“以N(X1,,)为提议函数”改 为“以N(X1,o,I)为提议函数”即可. 例3用例2的Bta(3,5)分布来检验上述自适应算法.取0。=10，m=2500,a=0.23,e=0.01,d=1, 计算得到符合要求的。值为0.82.图2为采样结果对比图，第一行为未经自适应算法改进得到的样本的自 相关系数图，轨迹图及依据样本估计的密度曲线；第二行为经过自适应算法改进得到的样本的自相关系数 图，轨迹图及依据样本估计的密度曲线. 从第一列自相关系数图的比较可看出经过自适应算法改进得到的样本要更符合实际计算需要，这表 现为其自相关系数图为一个很陡的下降曲线，并且在1g=20左右其自相关系数就已可以视为零了.相比 较而言，未经自适应算法改进得到的样本的自相关系数图为一个较平缓的下降曲线，且在1熙g=80时自相 关系数仍然很大，以致不能视为零，第二列两种算法轨迹图的比较显示经过自适应算法改进得到的Mkov 链移动的更快，从而能够更有效的搜索整个状态空间，从第三列估计的密度曲线图可明显看出，依据经过 自适应算法改进得到的样本估计出的密度曲线图要比依据未经过自适应算法改进得到的样本估计出的密 度曲线图要光滑. 4 Metropolis-Hastings算法与贝叶斯分析 贝叶斯估计方法就是把未知参数日视为一个具有已知分布π(x)的随机变量，从而将先验信息数字形 式化并加以利用的一种方法，通常称π(x)为先验分布，运用MH方法，我们可以从复杂的多维分布中采 样，从而能够对贝叶斯估计方法进行很好的模拟计算， 设总体多的分布密度为p(x,0),0∈⊙，0的先验分布为π(8)，由于0为随机变量并假定已知0的先 验分布，所以总体名的分布密度p(x,日)应看作给定日时多的条件分布密度，于是总体男的分布密度P (x,)需改用p(x|8)来表示.设X=(X1,…,X,)为取自总体多的一个样本，当给定样本值x=（1,…, 万方数据
104 系统工程理论与实践 2008年1月 1．O 0．8 函O·6 譬o．4 O·2 0．0 lignutffi0．05 O 20 60 。Lag 1．0- 1．01 O 20 60 Lag O 20 ∞ L略 图I取口=0．05，0．1，0．4，0．6，I，3，10，15时，随机游动采样法得到的样本自相关系数图．横坐标为滞 后系数189，纵坐标为对应的自相关系数，两条虚线表示自相关系数为零的95彩的置信区间． 则；可选择得小一些；若估计ar0与盯相差较远，则；可选择得大一些；对于d>1的多维情况，如果要确 定一个任意的正定矩阵：，则要确定d2个参量，显得有些麻烦．更多情形下，我们取三=01．这样一来要确 定的参数就只有一个，同时上述算法也可以用于多维情形，只须在6步中“以Ⅳ(五．。，吼)为提议函数”改 为“以Ⅳ(置一，，d。，)为提议函数”即可． 例3用例2的Beta(3；5)分布来检验上述自适应算法．取盯。=10，m=2500，口=0．23，s=0．01，孑=1， 计算得到符合要求的d值为0．82．图2为采样结果对比图，第一行为未经自适应算法改进得到的样本的自 相关系数图，轨迹图及依据样本估计的密度曲线；第二行为经过自适应算法改进得到的样本的自相关系数 图，轨迹图及依据样本估计的密度曲线． 从第一列自相关系数图的比较可看出经过自适应算法改进得到的样本要更符合实际计算需要，这表 现为其自相关系数图为一个很陡的下降曲线，并且在lag=20左右其自相关系数就已可以视为零了．相比 较而言，未经自适应算法改进得到的样本的自相关系数图为一个较平缓的下降曲线，且在lag：80时自相 关系数仍然很大，以致不能视为零．第二列两种算法轨迹图的比较显示经过自适应算法改进得到的Markov 链移动的更快，从而能够更有效的搜索整个状态空间．从第三列估计的密度曲线图可明显看出，依据经过 自适应算法改进得到的样本估计出的密度曲线图要比依据未经过自适应算法改进得到的样本估计出的密 度曲线图要光滑． 4 Metropolis．Hastings算法与贝叶斯分析 贝叶斯估计方法就是把未知参数口视为一个具有已知分布，r(茗)的随机变量，从而将先验信息数字形 式化并加以利用的一种方法，通常称，r(算)为先验分布．运用M．H方法，我们可以从复杂的多维分布中采 样，从而能够对贝叶斯估计方法进行很好的模拟计算． 设总体劈的分布密度为P(善，口)，口∈O，8的先验分布为丌(8)，由于8为随机变量并假定已知8的先 验分布，所以总体穷的分布密度p(聋，口)应看作给定口时劈的条件分布密度，于是总体劈的分布密度p (算，口)需改用P(茹l口)来表示．设x=(xI，．”，以)为取自总体留的一个样本，当给定样本值工=(菇l，．一， 万方数据
第1期 Metropoli-Hasting自适应算法及其应用 105 1.0 0. 2.5 0.8 0.b 2.0 0.6 15 04 0.4 1.0 02- 02 05 0.0- 0.0 020406080 0100020003000 0.00.40.8 (a)a-10 (b)trace (c)estimated density 1.0 0.8 08 2.0 0.6 0.6 15- 04A- 0.4 10- 02 0.2 0.5 0.0 0.0 0.0 020406080 0100020003000 0.00.40.8 (dg-082 (e)trace (f)estimated density 图2图(a)-(c)分别为采用一般算法(。=10)进行采样得到的样本的自相关系数图，轨迹图及依据 样本估计的密度曲线.图(d)-(0则分别为采用自适应算法(。=0.82)进行采样得到的样本的自相 关系数图，轨迹图及依据样本估计的密度曲线， )时，样本X=(X,,X.)的联合密度为：g(,,10)=亚p(10),或表示为：9(x10)= J. (18),从而，样本X和8的联合分布密度为：f(x,8)=q(xI8)π().由乘法公式知： f(x,8)=x(8)q(x1θ)=g(x)h(91x) 于是有： h(81x)=(9)g(x10) (x) 8∈g (3) 它是给定样本后日的条件分布.其中g(x)是(X,8)关于样本X的边缘分布.如果日是连续型随机变量， 则有：g(x)=J。(x10)x(9)d0:如果0是离散型随机变量，则有：8(x)=29(x1)x(9). 贝叶斯估计方法认为后验分布集中体现了样本和先验分布两者所提供的关于总体的信息，因而应在 后验分布的基础上来进行统计分析，在过去的实际应用中，由于g(x)不便于用分析的方法计算，使得贝 叶斯估计方法的应用大为受阻.为了方便计算而过度简化的模型常被用来凑合一下，MH方法的运用改变 了这一状况. 由于在贝叶斯估计方法中g(x)不依赖于日，在计算后验分布中仅起到一个正则化因子的作用，若把 g(x)省略，可将贝叶斯公式改写为如下等价形式： h(01x)xx(8)q(x1) (4) 其中符号“c”表示两边仅相差一个不依赖于日的常数因子.(4)式的右端虽不是正常的密度函数，但它是 后验分布h(91x)的主要部分.不难看出，M-H算法用于此类模型的一个主要优势在于只要(4)式右边已 知，即可从后验分布中采样，从而避免了计算g(x)的问题. 万方数据
第1期 Metropolis-Ha蛳自适应算法及其应用 105 I．o 。·s 。．6 室o．4 。·2 0．o 0 20 40 60 80 (吣口一lO 0 20 40 60 80 (d)口|o．82 0．0 0．4 0．8 (c)estimated demity 0．0 0．4 0．8 (O estimated demity 图2图(a)一(c)分别为采用一般算法(d=tO)进行采样得到的样本的自相关系数图，轨迹图及依据 样本估计的密度曲线．图(d)一(f)则分别为采用自适应算法(盯=O．82)进行采样得到的样本的自相 关系数图，轨迹图及依据样本估计的密度曲线． 髫。)时，样本X=(x1’．一，疋)的联合密度为：q(x。，…，毛f口)=Ⅱp(茗；l口)，或表示为：q(x I口)=ⅡP (毛I口)，从而，样本石和目的联合分布密度为：／(工，口)=q(引口)丌(口)．由乘法公式知： f(x，口)=石(口)q(x I口)=g(x)^(口I工) 于是有： ． ， h(O㈩=坐黔产，口∈@ (3) F～工， 它是给定样本后口的条件分布．其中g(x)是(X，∞关于样本石的边缘分布．如果口是连续型随机变量， 则有：g(x)=l q(x l口)，r(0)d0；如果口是离散型随机变量，则有：g(x)=三日q(x I口)7r(口)． 贝叶斯估计方法认为后验分布集中体现了样本和先验分布两者所提供的关于总体的信息，因而应在 后验分布的基础上来进行统计分析．在过去的实际应用中，由于g(z)不便于用分析的方法计算，使得贝 叶斯估计方法的应用大为受阻．为了方便计算而过度简化的模型常被用来凑合一下．M．H方法的运用改变 了这一状况． 由于在贝叶斯估计方法中g(x)不依赖于口，在计算后验分布中仅起到一个正则化因子的作用，若把 g(工)省略，可将贝叶斯公式改写为如下等价形式： 、 ．h(一I．It)％，r(口)q(工l a) (4) 其中符号“∞”表示两边仅相差一个不依赖于口的常数因子．(4)式的右端虽不是正常的密度函数，但它是 后验分布h(日I工)的主要部分．不难看出，M．H算法用于此类模型的一个主要优势在于只要(4)式右边已 知，即可从后验分布中采样，从而避免了计算g(x)的问题． 万方数据
106 系统工程理论与实践 2008年1月 在贝叶斯Logistic回归中假定样本观测值Y=(y1,y2,“,y,)服从参数为p的两点分布，即 ·P(y=1)=p(0),i=1,2,…,n 其中 p9）=ga(@)=,t=1,2,n (5) 并且， 月0+名1月++x13 (8,82,…,8.)r=8= 月。+xu月1+…+x2 人月。+x1月+…+xB 变换函数og(p)=log[p/(1-p)],logt-1为logit的反变换，X为n×k的矩阵，B为k×1的列向量.y:,i =1,2,…,n为取值为零或一的随机样本，x:为任意的实数， 例4失效的密封圈 1986年，“挑战者”号航天飞机在发射升空72秒后爆炸，造成7名宇航员罹难，事故原因怀疑是密封圈 失效，下表为航天飞机在不同温度下发射时密封圈失效的记录. 表1 航天飞机在不同温度下发射时密封圈失效的记录 飞行号 复 9 23 10 13 15 4 3 8 17 失效(=1) 1 1 1 0 0 0 0 0 0 0 0 温度（华氏） .53 57 58 63 66 67 67 67 68 69 70 70 飞行号 11 6 7 16 21 19 22 12 20 18 失效(=1) 1 .1 0 0 0 1 0 0 0 0 0 温度（华氏） 70 70 72 73 75 75 76 76 78 79 81 我们用贝叶斯L0gsc回归模型对表1中的数据做分析.设x表示发射时的温度，p(8)表示在不同温 度时密封圈失效的概率，y=1表示密封圈失效，y=0表示密封圈没有失效， 现在已有一组观测值(x,y:),i=1,2,…,23,按照模型假设 eo 为-B(1,p(8,)》并且p(0)=1+e可 先验分布取为：π(P)=1,似然函数为： L(BI Y) l÷+ (6) 后验分布正比于似然函数 我们先用类似于例1的随机游动采样法进行模拟，取提议函数为N(X,-1,1),构造一条长度为10000 的Markov链.模拟结果见图3. 图中显示估计A。的Markov链几乎没有平稳性可言.同时 8 0.95%的接受比率也使得序列的自相关系数非常高，因而，模拟结 0 果很难令人满意.观察发现回归系数B。与月1高度相关.这启发我 4 200040006000800010000 们如果对回归变量进行变换，可能会得到较好的结果，比如用x- 1.0- x作为新的回归变量设B。,P,为变换后的回归系数，则由：月。+ 0.5 0 1=B。+(x-*),可得： 03、 月。=。-，月= (7) 200040006000800010000 有了（门）式，我们就可以方便地转回到初始的回归系数值上.经过 图3上图为A的轨迹图， 变换后，模拟的结果有了很大的改进.从图象上看，P。,B,的轨迹 下图为B,的轨迹图 万方数据
106 系统工程理论与实践 2008年1月 其中 并且 在贝叶斯Logistic回归中假定样本观测值Y=(y。，Y：，…，Y。)服从参数为P的两点分布，即 。P(Y；=1)=P(0；)，i=l，2，…，厅 口． p(0；)=logi￡。1(仇)=南，}=1’2’2…，聘 (0，，0：，…，0。)7=邵= 风+菇ll卢l+…+髫l以 岛+髫2lpl+…+菇2以 风+髫。lpl+…+髫础& (5) 变换函数lo础(p)=log[p／(1一p)]，蛔咖。1为蛔咖的反变换，x为n×．|}的矩阵，卢为J}X 1的列向量．Y；，f =l，2，…，n为取值为零或一的随机样本，茗#为任意的实数． 例4失效的密封圈 1986年，“挑战者”号航天飞机在发射升空72秒后爆炸，造成7名宇航员罹难，事故原因怀疑是密封圈 失效．下表为航天飞机在不同温度下发射时密封圈失效的记录． 表1 航天飞机在不同温度下发射时密封圈失效的记录 飞行号 14 9 23 10 l 5 13 15 4 3 8 17 失效(=1) l l l 1 0 O O O O O O 0 温度(华氏) 53 57 58 63 66 67 67 67 68 69 70 ．70 飞行号 2 11 6 7 16 21 19 22 12 20 18 失效(：1) 1 ．1 0 O 0 l 0 0 O 0 0 温度(华氏) 70 70 72 73 75 75 76 76 78 79 8l 我们用贝叶斯Logistic回归模型对表1中的数据做分析．设菇表示发射时的温度，P(口)表示在不同温 度时密封圈失效的概率，Y=1表示密封圈失效，Y=0表示密封圈没有失效． 现在已有一组观测值(戤，Y；)，i：1，2，…，23，按照模型假设 Y；一B(1，p(B))并且p(∞=亡i≥ 先验分布取为：7t"(卢)毫1，似然函数为： 郴!叭熏(燕斋)^(d‰r ㈤ 后验分布正比于似然函数． 我们先用类似于例1的随机游动采样法进行模拟，取提议函数为Ⅳ(五一。，1)，构造一条长度为10000 的Markov链．模拟结果见图3． 图中显示估计岛的Markov链几乎没有平稳性可言．同时 0．95％的接受比率也使得序列的自相关系数非常高，因而，模拟结 果很难令人满意．观察发现回归系数凤与角高度相关．这启发我 们如果对回归变量进行变换，可能会得到较好的结果，比如用茗一 二作为新的回归变量．设∥。，∥。为变换后的回归系数，则由：风+ 筇l=∥o+(菇一菇)∥l可得： 母o=F q—xF l，pl=F l q、 有了(7)式，我们就可以方便地转回到初始的回归系数值上．经过 变换后，模拟的结果有了很大的改进．从图象上看，风，p。的轨迹 0 2000 4000 6000 8000 10000 图3上图为岛的轨迹图， 下图为风的轨迹图 o 4 o'r m帖o” 万方数据
第1期 Metropolis-Hastings自适应算法及其应用 107 都己表现出平稳的态势，其自相关系数也缩小到可以接受的范围以内.进一步研究发现，经过变换后，以 N(X-1,1)为提议函数进行模拟计算已经能够得到较为满意的结果，但在上面的模拟过程中，接受的比率 大约为10%，仍然偏低.为了进一步改进得到的样本的性质，我们用第3节中提到的自适应算法对这一问 题做模拟计算，从计算得到的A。,B,的轨迹对比图（图4）中易见自适应算法得到的样本轨迹比一般随机 游动算法得到的样本轨迹移动的更加迅速，也更加的平稳.自相关系数对比图（图5）则显示，不论是对B。 还是B,而言，自适应算法得到的样本的自相关系数都比一般随机游动算法得到的要小.这说明经过改进 的自适应算法更具有优势 10 40 0 03 20 00- 4料 0 50010001500200025003000 。30010001500202500300 .0 20 0、 -20- 0.0- 0- 0.5 -60- 050010001500200023003000 050010001500200025003000 图4一般随机游动与自适应算法得到的样本的轨迹对比图.第一行为一般随机游动算法得到的 风，月，的轨迹，第二行为自适应算法得到的A,月，的轨迹 0.8 08 0.4 0.044 4 66时 00 0 102030405060 0102030403060 08 0. 0.0ll4444sss5ss8939 0 102030405060 0 102030405060 图5一般随机游动与自适应算法得到的样本的自相关系数对比图.第一行为一般随机游动算法得 到的风，月，的自相关系数图，第二行为自适应算法得到的风，月的自相关系数图 值的一提的是，在运用上述自适应算法进行模拟时，取G=0.2的原因是由于初始值¤。得到的样本性 质已经比较好，所以我们估计较合适的。值与初始值。相差不大，另外，由于。值一定是正数，故在自适 应算法的第6步中以。，为中心进行随机游动时，若得到的是一个负数，则必须舍去， 基于以上得到的样本，我们可以对密封圈失效与温度之间的关系作出一些粗略的推断.图6为根据样 本估计得到的在不同温度下密封圈失效的概率密度.易见，温度为华氏50度时密封圈失效的概率（大部分 集中在“1”即失效附近)要明显大于温度为华氏0度时密封圈失效的概率，从而暗示导致密封圈失效的原 因很有可能是寒冷的天气. 万方数据
第1期 Meh叩oKs．m鲥n炉自适应算法及其应用 107 都己表现出平稳的态势，其自相关系数也缩小到可以接受的范围以内．进一步研究发现，经过变换后，以 Ⅳ(五一。，1)为提议函数进行模拟计算已经能够得到较为满意的结果，但在上面的模拟过程中，接受的比率 大约为10％，仍然偏低．为了进一步改进得到的样本的性质，我们用第3节中提到的自适应算法对这一问 题做模拟计算，从计算得到的风，卢，的轨迹对比图(图4)中易见自适应算法得到的样本轨迹比一般随机 游动算法得到的样本轨迹移动的更加迅速，也更加的平稳．自相关系数对比图(图5)则显示，不论是对风 还是卢，而言，自适应算法得到的样本的自相关系数都比一般随机游动算法得到的要小．这说明经过改进 的自适应算法更具有优势． 图4一般随机游动与自适应算法得到的样本的轨迹对比图．第一行为一般随机游动算法得到的 风，区的轨迹，第二行为自适应算法得到的卢o，卢。的轨迹 ’ 图5一般随机游动与自适应算法得到的样本的自相关系数对比图．第一行为一般随机游动算法得 到的岛，岛的自相关系数图，第二行为自适应算法得到的风，且的自相关系数图 值的一提的是，在运用上述自适应算法进行模拟时，取；=0．2的原因是由于初始值口。得到的样本性 质已经比较好，所以我们估计较合适的口值与初始值or。相差不大．另外，由于口值一定是正数，故在自适 应算法的第6步中以以为中心进行随机游动时，若得到的是一个负数，则必须舍去． 基于以上得到的样本，我们可以对密封圈失效与温度之间的关系作出一些粗略的推断．图6为根据样 本估计得到的在不同温度下密封圈失效的概率密度．易见，温度为华氏50度时密封圈失效的概率(大部分 集中在。1”即失效附近)要明显大于温度为华氏70度时密封圈失效的概率，从而暗示导致密封圈失效的原 因很有可能是寒冷的天气． 万方数据
108 系统工程理论与实践 2008年1月 157 10 5 o了h 0 0.00.20.40.60.8 0.002040.60.81.0 图6左图：温度为华氏0度时，根据样本估计得到的密封圈失效的概率密度；右图：温度为华氏50 度时，根据样本估计得到的密封圈失效的概率密度。 5结论 本文主要对MCMC方法中一类很重要的算法、即Metropolis-Hastings算法作了分析和改进.从理论上来 说，M-H算法构造的Markov链具有可逆性，从而保证了其唯一的极限分布为目标分布.从实际计算结果来 看，在提议函数选择合适的情况下，用文中提出的MH自适应算法模拟得到的样本能够很好地反映目标 分布的特性，虽然并非在任何情况下算法的效率都很高，但正如文中所述，只要我们能够对算法进行一些 实验分析，我们就有可能找到合适的提议函数使得算法的效率大大提高. MCMC方法对统计推断的发展的贡献是巨大的，尤其是在贝叶斯模型中，MCMC方法解决了长久以来 一直困扰贝叶斯模型求解的计算问题，本文以贝叶斯Logistic模型为例说明了M-H自适应算法在贝叶斯 模型中的应用.从分析结果上看，以后验分布为依据作出的统计推断与实际情况吻合的较好.，MCMC方法 是一种正在不断发展完善的计算方法，本文对MH算法做了比较分析及改进讨论，这将有助于在实践中 寻找合适高效的算法用于解决实际问题. 参考文献： [I Altaleb A,Chauveau D.Bayesian analysis of the Logit model and comparison of two Metropolis-Hastings strategics[]. Computational Statistics Data Analysis,2002,3(1):137-152. [2]Roberts GO,Rosenthal JS.Optimal scaling for various Metropolis-Hastings algorithms[J].Statistical Science,2001,6(4):351 -367. [3]Geweke J,Tanizaki H.Note on the sampling distribution for the Metropolis-Hastings algorithm[J].Communication in Statistics- Theory and Methods,2003,32(4):775-789. [4]Sawyer S.The Metropolitan-Hastings algorithm and extensions[J].Washington University,April 17,2004. [5]孟庆芳，张强，牟文英.混沌序列自适应多步预测及在股票中的应用[J].系统工程理论与实践，2005,25(12)：62- 68. Meng QF.Zhang Q,Mu WY.A novel adaptive multi-step-prediction met hod for chaotic time series and its applications in stock market[J].Systems Engineering-Theory Practice.2005,25(12):62-68. [6]Chen P.Some nonparametric estimator and their properties under the competing risks case[J].Sankhya:Indian Series A. 1998.60(2):293-304. [7]Chen P.Estimators and some behaviors for a partially linear model with censored data[J].Acta Mathematica Scientia,1999.19 (3):321-331. [8]Chen P.Yan F R,Wu Y Y,et al.Detection of outliers in ARMAX time series models[J].The 5th IIGSS Workshop.Wuhan.June, 2007,to appear [9]Chen P,Chen Y.The Identification of Outliers in ARMAX Models via Genetic Algorithm[J].The 5th IIGSS Workshop,Wuhan June,2007,to appear. [10]陈平，达庆利.运用SAS软件系统对我国农作物受灾及成灾面积的预测分折[J].系统工程理论与实践，2001,21(4)： 141-144. Chen P.Da QL.Prediction and analysis for the disaster area of crops in our country by SAS system[]].Systems Engineering Theory&Practice,2001,21(4):141-144. [11]陈平，王成震，周俊，等.运用SAS软件对上证指数月线数据的综合预测分析[J].系统工程理论与实践，2003,23 (6):36-41. Chen P.Wang CZ.Zhou J,et al.Prediction and analysis for the index of Shanghai stock exchange by SAS ystem on monthly data []]Systems Engineering-Theory Practice,2003,23(6):36-41. 万方数据
108 系统工程理论与实践 2008年1月 ● 窘3 童2 l O 15 鲁1。 善 5 0 图6左图：温度为华氏70度时，根据样本估计得到的密封圈失效的概率密度；右图：温度为华氏50 度时，根据样本估计得到的密封豳失效的概率密度． 5 结论 本文主要对MCMC方法中一类很重要的算法、即Metropolis．Hastings算法作了分析和改进．从理论上来 说，M．H算法构造的Markov链具有可逆性，从而保证了其唯一的极限分布为目标分布．从实际计算结果来 看，在提议函数选择合适的情况下，用文中提出的M．H自适应算法模拟得到的样本能够很好地反映目标 分布的特性．虽然并非在任何情况下算法的效率都很高，但正如文中所述，只要我们能够对算法进行一些 实验分析，我们就有可能找到合适的提议函数使得算法的效率大大提高． MCMC方法对统计推断的发展的贡献是巨大的，尤其是在贝叶斯模型中，MCMC方法解决了长久以来 一直困扰贝叶斯模型求解的计算问题．本文以贝叶斯logistic模型为例说明了M．H自适应算法在贝叶斯 模型中的应用．从分析结果上看，以后验分布为依据作出的统计推断与实际情况吻合的较好．MCMC方法 是一种正在不断发展完善的计算方法．本文对M．H算法做了比较分析及改进讨论，这将有助于在实践中 寻找合适高效的算法用于解决实际问题． 参考文献： ‘ [I] Ahaleb A，Chauveau D．Bayesian analysis of the IDgit model and comparison of two Meu-opolls-Hastings strategies[J]． Computational Statistics＆Data Analysis，2002，39(1)：137—152． [2]Roberts G 0；Rosenthal J S．Optimal sealing for various Meuupolis-l-lastings algorithms[J]．Statistical Science，2001，6(4)：351 —367． [3]Geweke J，Tanizaki H。Note on the sampling distribution for the Metropofis-Hastings algorithm[J]．Communication in StatisticsTheory and Methods，2003，32(4)：775—789． [4]Sawyer S．The Metropolitan-ttastings algorithm and extensions[J]．Washington University，April 17，2004． [5]盂庆芳，张强，牟文英．混沌序列自适应多步预测及在股票中的应用[J]．系统工程理论与实践，2005，25(12)：62— 68． Meng Q F，Z1mng Q，Mu W Y．A novel adaptive multi·step-prediction met hod for chaotic time series and its applications in stock market[J]．sy咖脚Engineering—Theory＆Practice，2005，25(12)：62—68．． 【6]Chert P．Some nonparametrie estimators and their properties under the咖peling risks c∞e[J]．S=khyi：Indian J Statist Series A， 1998，60(2)：293—304． [7] Chen P．Estimators and 80me behaviors for 8 partially linear model Wim censored data[J]．Acta Mathematica$cientia，1999，19 (3)：321—331． [8]ChenP，Yan F R，Wu Y Y，et a1．Detection ofoutliers in ARMAX time series models[J]．‰5th IIGSS Workshop，Wuhan．June， 2007．to appear． [9]Chen P，Chen Y．The Identification of Outliers in ARMAX Models via Genetic Algorithm[J]．11le 5th IIGSS Workshop，Wuhan June，2007，to appear． [10]陈平，达庆利．运用SAS软件系统对我国农作物受灾及成灾面积的预测分析[J]．系统工程理论与实践，2001，21(4)： 141—144． Chert P，Da Q L．Prediction and analysis for the disaster a脱of crops in哪country by SA$s，蜘[J]．Syst舢Engineering— Theory＆Practice，2001，21(4)：141—144． [11] 陈平，王成震，周俊，等．运用SAS软件对上证指数月线数据的综合预测分析[J]．系统工程理论与实践，2003，23 (6)i 36—41． Chen P，Wang C Z，Zhou J，et a1．Prediction and analysis for the index of Shanghai stock exchange by SAS system on monthly data [J]．sy咖娜Enoneermg—Theory＆Practice，2003，23(6)：36—41． 万方数据