第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 <a(x,y),则置X.1=y;否则置X+1=x.如果将一个使得π(x)≠0的点x作为整条Markov链的起始 点，那么(1)式中比值的分母将不会为零.原因是使得g(x,y)=0的新状态y被提出的概率为零，同时如 果π(y)=0,那么y被接受的概率为零，提议将被拒绝.所以只要π(X,)>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)产生一个随机数“，如果扯 <口(茗，Y)，则置疋+，=Y；否则置瓦+l=茹．如果将一个使得，r(茗)≠0的点茹作为整条Markov链的起始 点，那么(1)式中比值的分母将不会为零．原因是使得q(石，，，)=0的新状态Y被提出的概率为零，同时如 果巧(Y)=0，那么Y被接受的概率为零，提议将被拒绝．所以只要，r(x。)>0，那么在后续的计算中分母为 零的概率将是零，从而不会对实际计算造成麻烦． 如果M．H算法中的提议函数g(石，，，)不仅满足对称性，而且只与菇一Y有关，那么算法就演变为通常 意义下的随机游动采样法．最常见的一种随机游动采样法以正态分布，即口(茗，，，)=声(菇一，，)为提议函数． 这里声是任意一种多维正态分布的密度函数，也就是说Y一Ⅳ(茗，三)．(三是任意的正定矩阵，菇是当前的 状态)．在实际计算中，由于Z要求是正定矩阵，所以常取为三=al，口是一个参数．这时一个重要的问题是 口应取多少才合适．如果取的过大，那么大多数提出的新状态将落入分布的尾部，从而被拒绝．如果取的过 小，那么采样的自相关系数就会很高，而且收敛到目标分布的速度也会很慢．不论哪种情况，我们都不能在 时间有限的采样过程中得到较好地服从目标分布的随机数．一般认为在模拟多维正态随机数时，如果调整 盯的值使得在整个模拟过程中提议被接受的比例大约为20％，将得到比较好的模拟效果． 在时齐且有条件转移密度的情形，我们简记 P(茗，Y)暑P‘“。+u(茹，Y)=PW4’(茗，Y) 本文所考虑的Markov链都是时齐的，并总假定条件转移密度存在．M．H算法构造的Markov链的条件转移 万方数据