（概率论与数理统计专业论文）马链monte+carlo算法中的gibbs采样.pdf

中文摘要 摘要 上世纪9 0 年代以来，很多应用问题都存在着分析对象比较复杂与正确识别 模型结构的困难，而传统的经验方法由于不能逼近真实的过程，很难得到满意 的结果这其中涉及到很多计算问题，而这些计算问题却可以化成某个随机变 量函数的数学期望来解决但此数学期望往往难于计算，因此我们可以通过随 机抽样的方法，用样本均值来估计期望但在很多情况下，随机抽样在计算机 上也很难实现对于这样的随机变量，可以构造- - m a r k o v 链，此链以所给随机 变量的分布为平稳分布，且此链的抽样便于在计算机上实现，当链满足某种遍 历性时，就用链的样本( 用v o nn e u m a n 取舍原则进行取舍) 作为该随机变量的 近似这便足m a r k o v 链m o n t ec a r l o 方法，m c m c 的研究对建立可实际应用的统 计模型开辟了广阔的前景 本文主要介绍几种常见的m a r k o v 链m o n t ec a r l o 方法及其应用( 特别 是m e t r o p o l i s h a s t i n g s 算法与g i b b s 采样) ，以及学习这些方法的心得体会 关键词： c o nn e u m a n 取舍原则；m e t r o p o l i s - h a s t i n g s 算法；g i b b s 采样；l s i n g 模 型 湖北大学硕士学位论文 a bs t r a c t s i n c e19 9 0 s ，al o to f a p p l i e dp r o b l e m se x i s tt h ed i f f i c u l t yt h a tt h ea n a l y t i c a lo b j e c t a r ec o m p l e xa n di d e n t i f y i n gt h em o d e ls t r u c t u r ei sa l s od i f f i c u l t ，b u tb e c a u s et h et r a d i t i o n a le m p i r i c a lm e t h o dc a n ta p p r o a c ht ot h ea c t u a lp r o c e s s ，i ti sd i f f i c u l tt og a i n i n gs a t - i s f l e dr e s u l t s t h i si n v o l v e sal o to fc o m p u t i n gp r o b l e m ，b u tt h e s ec a l c u l a t i o np r o b l e m s c a nb et r a n s f o r m e di n t om a t h e m a t i c a le x p e c t a t i o no f s o m ef u n c t i o no f r a n d o mv a r i a b l e ， b u tt h i sm a t h e m a t i c a le x p e c t a t i o na r eu s u a l l ym o r ed i f f i c u l tt ob ed e a l tw i t h t h e r e f o r e b yt h er a n d o ms a m p l i n g ，w ec a nu s es a m p l ea v e r a g et oe s t i m a t et h ee x p e c t a t i o n b u t u n d e rm a n yc o n d i t i o n s ，t h er a n d o ms a m p l i n ga r ea l s ov e r yd i f f i c u l tt or e a l i z a t i o n f o r t h i sr a n d o mv a r i a b l e ，w ec a nc o n s t r u c to n em a r k o vc h a i n ，t h i sc h a i n ss t a t i o n a r yd i s t r i b u t i o ni st h eg i v e nr a n d o mv a r i a b l e ，a n dt h es a m p l i n go ft h i sc h a i ni se a s yt ob e c o m p u t e d ，w h e nt h ec h a i ns a t i s f i e ss o m ee r g o d i cp r o p e r t y , w ec a l la p p r o x i m a t e l yu s et h e s a m p l e s t or e p l a c et h er a n d o mv a r i a b l e t h i si ss o - c a l l e dm a r k o vc h a i nm o n t ec a r l oa l - g o r i t h m s ，t h em c m c r e s e a r c hd e v e l o p p e dv a s tf o r e g r o u n dt ob u i l d i n gu pt h ea p p l i e d m o d e l t h i st e x tm a i n l yi n t r o d u c e ss o m ef a m i l i a rm a r k o vc h a i nm o n t ec a r l oa l g o r i t h m s a n dt h e i ra p p l i c a t i o n ( g i b b ss a m p l i n ge s p e c i a l l y ) ，a n ds o m ee x p e r i e n c eo f m yo w n k e yw o r d s ：r e j e t i o ns a m p l i n go fv o nn e u m a n ；m e t r o p o l i s h a s t i n g sa l g o r i t h m ； g i b b ss a m p l i n g ；l s i n gm o d e l 湖北大学学位论文原创性声明和使用授权说明 原创性声明 本人郑重声明：所呈交的论文是本人在导师的指导下独立进行研究所取得 的研究成果。除了文中特别加以标注引用的内容外，本论文不包含任何其他个 人或集体已经发表或撰写的成果作品。对本文的研究做出重要贡献的个人和集 体，均已在文中以明确方式标明。本人完全意识到本声明的法律后果由本人承 担。 论文作者签名：专湛从 签名r j 胡：一彦年6 月厶卜日 学位论文使用授权说明 本学位论文作者完全了解学校有关保留、使用学位论文的规定，即： 按照学校要求提交学位论文的印刷本和电子版本；学校有权保存并向国家 有关部门或机构送交论文的复印件和电子版，并提供f i 录检索与阅览服务；学 校可以允许采用影印、缩印、数字化或其它复制手段保存学位论文；在不以赢 利为目的的前提下，学校可以公开学位论文的部分或全部内容。( 保密论文在 解密后遵守此规定) 论文作者签名：东：、：，扯 导师签名：了龟钮：乙 签名日期：汕。辟6 月钿 签名日期：叼年6 月乒日 第一章引言 1 1 方法简介 第一章引言 m o n t ec a r l o 方法亦称为随机模拟( r a n d o ms i m u l a t i o n ) 方法，有时也称作 随机抽样( r a n d o ms a m p l i n g ) 技术或统计试验( s t a t i s t i c a lt e s t i n g ) 方法【l 】1 2 】 它的基本思想是：为了求解数学、物理、工程技术以及生产管理等方面的问 题，首先建立一个概率模型或随机过程，使它的参数等于问题的解；然后通过 对模型或过程的观察或抽样实验来计算所求参数的统计特征，最后给出所求 解的近似值，而解的精确度可用估计值的标准误差来表示因为在概率论初 期发展时，随机模拟方法的原型常常来自博采，于是人们就以博采之都m o n t e c a r l o 作为随机模拟方法的别称久而久之，m o n t ec a r l o 方法作为名称倒比随机 模拟方法更为广泛的应用了【3 】【4 】 m c m c ( m a k o vc h a i nm o n t ec a r l o ) 方法最早应用于计算物理领域 ( m e t r o p o l i s 等1 9 5 3 ) ，h a s t i n g s ( 1 9 7 0 ) 的工作使其更为一般化，但主要应用于统 计问题g e l f a n d 和s m i t h ( 1 9 9 0 ) 的研究显示出m c m c 方法在贝叶斯计算上有着巨 大的潜力，他们将马尔科夫链的方法与t a n n e r t g l w o n g ( 1 9 8 7 ) 的增参数方法相结 合，证明了对含有潜在动态变量的计算模型的估计是极为成功 为要模拟服从给定分布丌的随机变量，可以生成一个易于实现的不可约遍 历链，使其平稳分布为丌，用它作为随机抽样，基于这些样本就可以作各种 统计推断，这种方法称为马氏链蒙特卡罗方法( m c m c ) 这是动态的m o n t e c a r l o 方法由于这种方法的问世，使随机模拟在很多领域的计算中，相对于决 定性的算法，显示出它的巨大优越性而有时随机模拟与决定性算法的结合使 用，会显示出更多的长处g i b b s 抽样和m e t r o p o l i s h a s t i n g s 算法足两种最为常用 的m c m c 方法 g i b b s 采样足一种m c m c 算法最早由g e m a n 和g e r m a n ( 1 9 8 4 ) 在一篇讨论图 象恢复的文章中提出在统计上，t a n n e r 弄c i w a n g ( 1 9 8 9 ) 将该算法用于遗失数据 问题接着，g e l f a n d 平t l s m i t h ( 1 9 9 0 ) 指出了g i b b s 抽样在贝叶斯计算中的应用 目前，g i b b s 抽样算法处理复杂统计模型的能力正逐渐为人们所认识【5 】 g i b b s 采样的原始模型是i s i n g 模型，用于解决格点阵列上能量函数的求 解问题，但在维数比较大时，直接求解计算量足惊人的于是我们可以通 湖北大学硕士学位论文 过此i s i n g 模型中的q 矩阵来构造一个m a r k o v 链，由于对不可约马链有华一 仇，( t _ 0 ) ，于是当a i t t d 时有，p o ( a t ) = q j a t 6 1 由此我们就可以得到一系 列的样本，再通过样本计算频率，就可以作出估计此方法的特点足：每一次 只在一个格点上发生改变在一般情形下的g i b b s 采样，可以采用条件分布族来 得到以7 r 为不变分布的m a r k o v 链，一次固定其它分景只改变一个坐标，以此可 以大大的减少计算量在第三章我们将它i 孑一般抽样方法的计算量进行比较 1 2m c m c 的作用 就求解数学和物理问题而言，m c m c 方法比过去的随机试验方法增添了许 多新的内容，至少有如下三点是值得提到的： ( 1 ) 应用概率模型模拟物理现象或数学方程，扩大了随机抽样的应用范 围 ( 2 ) 整个随机试验足在计算机上进行，这就使随机试验方法有可能更接近 真正的物理试验或其他数学物理方程解 ( 3 ) 在计算机上用程序产生伪随机数代替真正的随机数，不仅使用方便而 且提高了计算速度 随着计算机使用范围的日益广泛，它向各个学科的渗透也越来越深入，近 几年兴起的计算物理学【7 j 、计算概率统计学f 8 】f 9 】、计算机科学与统计学的界面 学【10 】等边缘学科，它们都署u m c m c 方法有着密切的关系 1 3m c m c 可以解决的问题 ( 1 ) 用于生成较复杂的随机数 实现对高维分布( 或高维格点分布) 7 r 的取样，得到丌随机数 是实现重要度采样的一种方法，对i ，( z ) i 的重要度采样，就是取得 丌一( x ) a e f = 摆施随机数 对于厂( z ) 0 ，作以丌= 丌( z ) 型刊箍为极限分布的m a r k 0 v 链，利用遍历 定理可以由这个m a r k o v 链的一条轨道，得到分布密度丌( z ) 的估计，记为亓( z ) ( 2 ) 实现高维积分( 或者项数极多的求和) 的数值计算( 典型的是g i b b s 分 布的各种泛函的甲均值的计算) 用籍作为积分j = r ，( z ) 出的估计 2 第一章引言 ( 3 ) 用模拟方法估计最可几轨道例如，如果模拟了1 0 0 条轨道，那么就 能以人概率推断，最可几轨道就在这些轨道的附近当统计量的分布未知时， 可以用模拟方法从频率估计置信限 ( 4 ) 用被估计参数的b a y e s 分布的取样，来估计参数【1 1 1 ( 5 ) 求复杂样本空间上函数的极值( 模拟退火【1 2 ) 以上问题用m c 方法很难解决 3 湖北大学硕- 上学位论文 2 1 随机样本的生成法 第二章预备知识 弟一早 耿亩大u 以 随机变量( 或随机向量) 的样本简称为随机数由于在统计中常用的是独 立样本列，不妨假设随机数之间都是独立的生成随机数的方法，也称为随机 数的取样法 2 1 1 一维均匀随机变量的计算机模拟 定义在【o ，1 1 上均匀分布的随机变量的独立样本称为均匀随机数( v o ，1 】随 机数) 在计算机上产生的称之为“伪随机数【1 3 】”的数列，是一一种具有非常长周期 的，且能通过数理统计中的独立性与均匀性假设检验的数列实践证明，伪 随机数足均匀随机数的一种可行的近似这种伪随机数虽然不足独立同分布 的v o ，1 】随机数随机变量的样本，而是在【o ，1 】中取值的周期数列，但足由于它 可以像均匀随机数一样地通过数理统计中的独立性与均匀性假设检验，而且它 的周期非常长，以至于在计算机实际运算过程中不会出重复，所以在实际计算 中它能很好的代替均匀随机数 最常用以产生伪随机数的方法是同余法典型的例子如下： + l = 5 1 3 y n ( m o d 2 3 6 ) ，y o = 1 ，z 。= 2 - 3 6 ( 周期约为2 1 0 1 0 ) + 1 = 5 1 7 y , , ( m o d 2 4 2 ) ，y o = 1 ，z 。= 2 - 4 2 蜘( 周期约为2 1 0 1 2 ) 鲰+ 2 = y n + 1 + y n ( m o d 2 4 4 ) ，y o = 1 ，y l = 1 ，z n = 2 一“鲰( 周期约为 1 0 1 4 ) 2 1 2 多维随机数 ( 1 ) 连续型多维随机数 对于已知分布密度，可以利用条件密度，把生成多维随机数归结为生成一 系列一维随机数：设随机向量( x 1 ，托) 的密度x j f ( x ，x 2 ，x d ) ，则有表达 式 ，( z l ，x 2 ，x d ) = 厶l ( x 1 ) f ( x 2 x 1 ) f ( x d l x l ，x d 一1 ) ， 第二章预备知谚! 其中丘，( x 1 ) 为x l 的边缘密度，( z k l z l ，- ，巩一1 ) 为在已知x l = z l ，x 女一1 = x k - i 条件下，x k 的条件密度于是可以先取一个厶。随机数z ；然后，在 固定z 。的情形下生成。个，( i z 。) 随机数；再在z 1 ，x 2 固定的情形下，生成一 个，( i x l ，x 2 ) 随机数；最后，在固定z 1 ，z d l 的情形下，生成一个随机数x d ， 这样得到的就足随机向量( x l 一，) 的一个随机数当然在生成各个条件密度 的随机数时，仍然可以使用n o nv e u m a n n 取舍原则注意在用取舍原则于一维 分布取样时，可以忽略一个因子 然而，在d 较大时，更常用的是g i b b s 取样法 在很多实际情形中，多维密度常常并不知道，更无法知道各个条件密度 这时最自然而粗糙的想法是用条件频率来近似条件密度 ( 2 ) 离散型多维随机数 只要用概率函数p ( z l ，z 2 ，x d ) a = e tp ( x 1 = z l 一，x d = 黝) 代替密度函 数，就可以采用连续型的方法 2 2 直接抽样方法 2 2 1 直接抽样方法的一般形式 对于给定的分布函数f ( z ) ，直接抽样方法的一般形式如下【1 4 】： 矗2f 豫。，肛, 2 , - - - ， 其中r l ，吃，为随机数序列 下面证明上式所确定的随机变量序列。，已，矗是由总体分布f ( z ) 中产生 的简单子样，只须证明它们具有同分布f ( z ) 且相互独立即可 p ( n z ) = p ( f 【t i n ) 2 f ，。 q ( y ，z ) 7 r ( ) ，则口( z ，3 f ) = 渊，口( 3 ，z ) = 1 ，所以 p ( x ，y ) 丌( z ) = q ( x ，y ) 口( z ，秒) 丌( z ) = g c z ，y ，； 兰黝丌c z ， = q ( y ，z ) 7 r ( y ) = q ( y ，x ) q ( y ，z ) o ( 剪，z ) 丌( y ) = p ( v ，z ) 丌( 3 ) 湖北大学硕士学位论文 ( 3 ) 若g ( z ，可) 丌( z ) q ( v ，z ) 丌( y ) ，贝u a ( y ，z ) = 渊，n ( y ，z ) = 1 ，所以 p ( y ，z ) 丌( 秒) = q ( y ，z ) q ( 可，z ) 7 r ( 可) = 口( y ，z ) 詈踹丌( 箩) = q ( x ，! ，) 7 r ( z ) = 口( z ，) q ( z ，可) 7 r ( z ) = p ( x ，可) 7 r ( z ) 3 1 2 应用举例 例考虑i 眨x 2 - 分布【1 9 】 p ( o ) = c 口一号e 一， 我们想从扎= 5 ，n = 4 的分布抽样 选取参考密度x j u o ，1 0 0 】，选取o o = 1 作为初始值，从u o ，x o o 中得到参考 值0 + = 3 9 8 2 ，则 n n ( 褊，) 一n 卜3 9 8 2 - 。z _ _ 钆5 e _ e 矗s 者8 2 ：，- ) = 0 - 蝴 由于a 1 0 1 8 咖) 个元素) ，这就使得分母中的求和无法实际 完成 m c m c 方法就是以通过构造一个以这个g i b b s 分布飞= 气昱；为不变分 嵇 布的离散时间的m a r k o v 链x n ，当n 充分大时，的一个样本可以近似地认为 取i ! i g i b b s 分布的一个样本，即按此m a r k o v 链沿任意一条轨道充分发展，就得 至u g i b b s 分布的近似取样 问题2g i b b s 分布的归一化常数( 称为配分常数) e 一卢h ( 叩) 是一个巨大的 刀s 求和，即一个“离散”的积分用随机模拟法计算这个”离散的”积分的最佳 随机数正服从g i b b s 分布( 即重要度采样) 对于g i b b s 分布的取样，用通常的 取舍原则常常并不可行 例如，参考密度为组态空间上的均匀分布时，取值常常小得超出计算精 度，而无法实际计算所以需要用m a r k o v 链m o n t ec a r l o 方法用m c m c 方法生 成了以g i b b s 分布为极限分布的m a r k o v 链以后，由遍历定理用m a r k o v 链的一条 轨道，可给出极限分布( g i b b s 分布) 的估计【2 2 】：对于充分大的n ，可令 = 丙1 ( 酬x + 1 ) + + 比) ( 恐) ) ， 还可以用e 一口日5 除以g i b b s 分布在处的估计值兰二警，作为配分函数的估 计在理论上这个估计应该与的取法无关但足，在实际计算中对多个不同的 组态，分别估计此和数后，再作平均常常能降低方差 用以完成m c m c 采样操作的m a r k o v 链，可构造如下： 2 0 第三章m a r k o v 链m o n t ec a r l o 方法 对于d 维台限格点上，由具有两个自旋的组态空问上的能量函数 日( 叩) = 一去7 7 ( z ) 叩( y ) 叩( z ) ， o ，暑，相邻 z 可构造如下的转移概率速率 c ( x ，f ) 的两种取法各为 fc ( z ，7 7 ) ， 叩= 专z ， 必，72t o ，其它情形 淝) ：e 一口( 日( r ) 一删：e 一烈。懒们+ 驴 或 1 口( z ，) 2 再i 疬丽 它们决定的g l a u b e r z 力力学，分别对的两个不同的连续时间的m a r k o v 链，它们都 以日( 7 7 ) 为能晕函数的g i b b s 分布7 r 为可逆不变分布，而且7 r 还是转移矩阵的极限 分布：p ( t ) ! 驾f r 丌 如果考虑问隔为a t 的采样时间，其中充分小，我们还可以进行如下的近 似计算 ( 1 ) 假定初始组态为( ，在时刻t 它以概率c ( x ，o a t 转移到组态r ，而以概、 率1 一t ec ( x ，( ) 停留在原来的组态，这就近似地得至l j t m a r k o v 链在时n 刻所处的组态( 1 ； ( 2 ) 再以类似的方式得至l j m a r k o v 链在2 a t 时刻所处的组态( 2 ； ( 3 ) 依次下去，当采样数n 充分大时，组态叩在这段中 臼，白) 出现在频 率，就用作的估计南 2 1 湖北大学硕士学位论文 第四章g i b b s 采样法的几个问题 4 1 利用可逆链防止多峰情形的偏取 在g i b b s 样本生成过程中，只有当m a r k o v 链 k ) 发展的“时间”r l 充分大 时，它的分布才近似于骨在实践中怎样才能判断n 足否已经足够大了呢? 通常 的判断标准一般足看在计算中的 k ) 的分布是否比较稳定了但是，在事实上 这种判断方法不尽正确下面的示意图给出了7 r = 7 r ( z ) 具有三峰的情形 从不同初值出发的轨道，可能分布稳定在e l ，易与马三个区域，而事实 上 还远未达到“近似于”不变分布的要求也就是说，按上一段的方法得 到的m a r k o v 链，从不同的初值出发，在计算过到看起来稳定时，它们还远不能 近似的样本原因在于这时的平稳分布丌( z ) 有三个峰，且在任意两峰中间转移 的概率都非常小，从而在一条轨道上从其中一个区域到达另外一个区域就要花 十分长的时问，有时甚至在现实生活中可能不会在有意义的时间内到达发现 与使用g i b b s 样本生成法的先驱们往往希望通过选取多个初始值，得到多条轨道 以解决这一问题，但是这并不能完全解决问题正如在当前的例子中的三条轨 道各自只局限在不同有区域易，e 2 和岛中，它们的分布，根本反映不出分别限 制在易，易和马上的7 r = 7 r ( z ) 足如何拼接的，即以多大的比例c 1 ：c 2 ：c 3 ，使 ( g + q + g ) 丌g i r i s l + c 2 7 r l s 2 + c 3 7 r l e , 只有找到了这个比例c 1 ：q ：( 7 3 ，才能从分别限制在局，易和岛上 的= i s 。，r i l e 。和丌l 岛，拼接出整体的丌( z ) 为了找出比例c 1 ：q ：c 3 ，我们利用可逆m a r k o v 链来生成g i b b s 样本，并从 可逆m a r k o v 链的特性得到不变分布在它的各个峰值上的比例 下面我们以三峰时三条轨道的情形，来阐述确定与估计这个比例的方法 设以p ( z ，秒) 为转移概率的可逆m a r k o v 链，有分别稳定在e 1 ，易和晶上 的三条轨道的样本分布对应于分布丌l ( z ) ，7 1 2 ( z ) 和丌3 ( z ) 我们近似地 令7 r ( z ) 在局，易和局外为0 于是由可逆m a r k o v 链的性质可知，这样应该分 一2 2 第四章g i b b s 采样法的几个问题 e i玛砀 图4 1 7 r = 7 r ) 具有三峰 别与 丌( z ) ：z 易) ， 7 r ( z ) ：z 五2 及 7 r ( z ) ：z e 3 觥t 例，即 那么 嘶) 怒，从而嵩= 鬻( v 训最， 2 3 ) 所以我们只需要求出几( z ) 的估计就可以了 具体确定g 可如下进行：对z ，y 晟，有 于足有 互嘶) _ ，觜( 可逆性) z e iz e t i 、1 。 哟，互端 2 几( y ) 邑【丽1 而】 1 2 3 湖北大学硕士学位论文 响) ( 岛【高】) 只要有x l 的n 次取样x i ，x ：( 忍中) ，我们可以得到几的估计 锄h 专娄南广 以上说明只要构造了以丌为可逆分布的m a r k o v 链k 转移概率p ( x ，) ，就可 以近似得到c 1 ：岛：岛和几，从而能近似地拼接出整体的丌这