(概率论与数理统计专业论文)mcmc方法研究.pdf_第1页
(概率论与数理统计专业论文)mcmc方法研究.pdf_第2页
(概率论与数理统计专业论文)mcmc方法研究.pdf_第3页
(概率论与数理统计专业论文)mcmc方法研究.pdf_第4页
(概率论与数理统计专业论文)mcmc方法研究.pdf_第5页
已阅读5页,还剩32页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

m c m c 方法研究 赵琪 山东大学数学与系统科学学院 济南2 5 0 1 0 0 z a h 0 2 3 1 2 0 0 1 6 3 c o m 摘要 本文丰要考虑了马氏链蒙特卡洛( m c m c ) 方法,主要是m e t r o p o l i s 和g i b b s 算法。 介绍了这些算法的过程,分析了马氏链的收敛性,讨论了m c m c 方法的误差,最后讨 论了m c m c 方法在贝叶斯模型中的应用。本文共分四部分:第一部分介绍了静态的 蒙特卡洛方法,主要是频率法与期望法,还讨论了如何提高这些方法的精确度。第二 部分丰要介绍了两种特殊的m c m c 方法,即m e t r o p l i s 和g i b b s 算法,包括这些算法的 详细的实施过程,这些算法的关系及如何在应用中改善这些算法。第三部分主要讨论 了m c m c 方法的三个重要方面,即收敛性、误差和策略问题。在收敛问题中讨论了三 种假设的马氏链的情况。在误差问题中讨论了如何降低误差。在这一部分的最后讨论 了如何在应用中选择一个合适的策略。第四部分主要介绍了m c m c 方法在贝叶斯模型 中的应用。 关键词:m c m c 方法,m e t r o p o u s 算法,g i b b s 算法,分层贝叶斯模型 s t u d yo fm c m cm e t h o d z h a oq i s c h o o lo fm a t l l e m a t i c sa n ds y s t e ms c i e n c e s ,s h a n d o n gu n i v j i n a n2 5 0 1 0 0 z a l l 0 2 3 1 2 0 0 1 6 3 c o i i l a b s t r a c t t h i sp a p e rc o n s i d e r st h em a r l c h a i nm o n t ec 舡l om e t h o d 8m a i n l yt h em e t h o - p l i 8a l g o r i t h ma n dt h eg i b b s 出g o r n h mi n t r o d u c e st h ep r o c e s so ft h e s ea k o r i t h m a n a l v z et h ec o i l v e r g e n c eo ft h em a r k o vc h 咖,d i s c u s st h ep r o b l e mo fa s 8 e s s i n gt h e e r r o ro fm c m cm e t h o d s 。i nt h ee n dw ed i s c u s st h ea p p l i c a t i o no fm c m cm e t h o d si n h i e r 舡c h i c a lb a v e s i a nn l o d e t h e 矗r s ts e c t i o nn l a i n l yi n t r o d u c e st h es t 砒i cs t a t em o n t e c 盯l om e t h o d sm a i l l l yi n t r o d u c e sf r e q u e n c yi i l e t h o d sa n de x p e c t a t i o i lm e t h o d sa l s o d i s c u s st h ep r o b l e m so fh o wt oi m p r d v et h ep r e c i s i o i li nt h e s em e t h o d s t h es e c o n d s e c t i o nm 出i l l yi n t r o d u c e st w oe x a m p l e so fm c m cm e t h o d si e m e t r o p l i sa l g o r i t h i n a i l dg i b b sa l g o r i t h m ,i n c l u d et h ed e t a i l e dp r o c e 8 so ft h e s ed g o r i t h ma i l dt 1 1 er e l a t i o n o ft h e s ea l g o r 扎h i n w ea l s od i s c u s gh o wt oi m p r o v et h e s ea l g o r i t h n li na p p l i c a 上i o nt h e t h i t ds e c t i o nm a i n l ym t d u c e st h r e ! ea s p e c to fm c m cm e t h o di e c o i e r g e n c e 、e r r o r 、a n ds t r a t e g 矿w ed i s c u s st h er e a l i z a t i o n so ft l l r e eh y p o t h e t i c a lm a r k o vc h a i ni n t h ep r o b l e mo fc o n v e r g e n c e w ea 1 8 0d i s c u s 8h o wt or e d u c et h ee r r o ro fm o n t ec a r l o e s t i m a t e s i nt h ee n do ft h i ss e c t i o nw ed i 8 c 1 】8 sh o wt oc h 商p eaa p p r o p r i a t es t r a t e g y i nas i m u l a t i o nr u n st h ef o u r t hs e c t i o nm a i n l yi n t r o d u c e st h ea p p l i c a t i o no fm c m c m e t h o d si nb a y e s i c a l lm o d e k e yw o r d s :m a r k o vc h a i nm o n t ec a r l o ,m e t r o p o l i sa l g o r i m m ,g i b b sa k o m h m ,h l - e r a 珀c h i c a lb a y e s 原创性声明 本人郑重声明:所呈交的学位论文,是本人在导师的指导下,独 立进行研究所取得的成果。除文中已经注明引用的内容外,本论文不 包含任何其他个人或集体已经发表或撰写过的科研成果。对本文的研 究作出重要贡献的个人和集体,均已在文中以明确方式标明。本声明 的法律责任由本人承担。 寸h 论文作者签名:壑! 鉴日期:地2 :丝丝 关于学位论文使用授权的声明 本人完全了解山东大学有关保留、使用学位论文的规定,同意学 校保留或向国家有关部门或机构送交论文的复印件和电子版,允许论 文被查阅和借阅;本人授权山东大学可以将本学位论文的全部或部分 内容编入有关数据库进行检索,可以采用影印、缩印或其他复制手段 保存论文和汇编本学位论文。 ( 保密论文在解密后应遵守此规定) 论文作者签名:煮丛导师签孝兰尹舭 期:塑! z 丝丝 第一章蒙特卡洛方法 1 1 引言 l 随机模拟方法是通过随机变量函数的概率模拟、统计试验或随机抽样求解数学物 理、工程技术、生产管理等各种不同问题的近似数值解的概率统计方法。随机模拟的模 型常常来自博采。因此物理学者给了它一个更新颖的名字l 蒙特卡洛( m o n t ec ”l o ) 方 法。1 7 7 7 年,法国数学家b u 肋n 的计算丌值的投针试验可视为近代蒙特卡洛方法的雏 形。由于随机试验费时费力,因此在科学问题中真正系统地运用蒙特卡洛方法是计算 机诞生以后的事情。有了计算机后就无需具体实施这些随机试验,只要在计算机上进 行模拟计算即可。如果模拟过程中不改变抽样分布就叫做静态蒙特卡洛方法,这只能 解决一些较简单的问题。在上世纪五十年代,统计物理学家m e t m p o l i s 等人引入了马 氏链蒙特卡洛方法( m c m c ) ,它是一种动态的蒙特卡洛方法,可用于处理许多复杂的 物理系统。此后几十年里有许多学者又发展了这种方法,使它应用到各种不同的领域。 特别是上世纪八十年代以后,计算机技术的蓬勃发展极大地推动了这种方法的广泛应 用。象组合最优化,缺失观测的似然计算、非参数统计推断、统计遗传学、贝叶斯模型 及其计算等。1 9 8 4 年,g e m a n ,s 和g e i n a n ,d 引入了g i b b s 算法,现已成为一种标准化 的统计计算t 具。本文主要介绍了m e t r 叩0 1 i s 算法与g i b b s 算法,在m e t r o p o h b 算法中 对不同的建议分布进行了对比分析。此外还对m c m c 方法的收敛性、误差,策略和应 用等问题进行了一些讨论。 1 2 静态蒙特卡洛方法 许多科学问题的基本部分是计算下面形式的积分 ,= ,( z ) 如( 1 ) 三 这儿的,0 ) 是感兴趣的目标函数。d 通常是高维空间中的一个区域。通常情况下如果问 题稍微复杂一些,这个积分就难以计算。因此就要另辟蹊径。如果,( 叫能被分解成一个 函数9 ( z ) 与一个概率密度甬数”( ) 的乘积,那么上述积分可看做是9 ( z ) 在密度0 ) 下 的期望。 ,= 厂m ) 如= 厂如) m ) 出= 晶始) 2 如果我们能得到7 r ( z ) 的样本x ( ”,x ( m 就获得珀q 近似值厶 厶= :喜) 这就是m o n t ec a r l 0 积分,当x ( ”x ( n ) 独立时,由大数定律有 1 i m 矗= ,( 以概率为l 地) 对于( 4 ) 式,它的收敛率能由中心极限定理来估计 、,伍( 厶一,) 一( o ,矿2 ) ( 依分布收敛)( 5 ) 这儿矿= 0 r g ( x ) ) ,因此m o n t ec a r l 0 近似值的“误差项”是d ( :;= ) ,不论x 的维数是多 少,这使得m o n t ec a r l o 方法比通常的计算方法有明显的优越性。 m o n t ec a r l o 方法的基本思路 ( 1 ) 针对实际问题建立一个简单且便于实现的概率统计模型,使所求的解恰好足 所建模型的概率分布或其某个数字特征,比如是某个事件的概率,或者是该模型的期 望值。 ( 2 ) 对模型中的随机变量建立抽样方法,在计算机上进行模拟试验,抽取足够的 随机数,并对有关的事件进行统计。 ( 3 ) 对模拟试验结果加以分析,给5 所求解的估计及其精度( 方差) 的估计。 ( 4 ) 改进模型以降低估计方差和减少试验费用,提高模拟计算的效率。 下面以估计简单积分,= e ,( z ) d z 为例说明m o n t ee a r 】o 方法的几种具体实旌办 法以及如何提高估计精确度。 l 、以频率估计概率的方法( 随机投点法) 设o ! ,( t ) 曼m 冷q = h 叫【0 ,m 】,是q 卜的均匀分布的二维随机变量,则f 的 密度函数为: p ( 砌) :而与“蜓蜓玩o s ) o ,( 对垤 o ,n o o ) ( 8 ) 即,是珀q 相合估计 具体实施步骤是: ( 1 ) 在计算机上独立产生2 n 个u ( 0 ,1 ) 随机数“i ,饥,t = 1 ,肌 ( 2 ) 计算戤= n + 啦( 6 一口) ,玑= m 和,( 戤) ( 3 ) 统计玑,( 现) 的个数n o ( 4 ) 用( 6 ) 式估计 下面看一下此方法的精确度( 方差) 由于n o 一6 ( n ,p ) ,( 其中p = 丽;= 鬲) ,所以有e j = e ,( z ) 如即,是珀q 无偏估计。 ”n r ( j ) = 掣”。r ( n 。) ,= ;【a ,( 6 一n ) 一j 1 ( 9 ) 由( 9 ) 式,若用估计的标准差来衡量其精确度,则,的精确度的阶为n 孚 2 、期望法 期望法的核心是把积分看成某个随机变量的期望 例如对积分,= r ,( z ) 如,设g ( 。) 是( n ,6 ) 上的一个密度函数,则有 ,= z 6 鬻咖肛e 【筹】 ( 1 0 ) 最简单地设x u ( 口,6 ) ,9 ( 。) 是( n ,6 ) 上均匀分布密度函数,即g ( z ) = 话,设z - ,册。是 来自矿f o ,6 ) 的随机数,则,的一个估计为 4 j 显然是,的无偏估计 ( 1 1 ) 酬肛( 埘竺笔盟- ( 6 刊2 知肛) ) 2 书肛) 2 ) _ ( 埘:p ) 2 圭出一:尸 ( 1 2 ) u r ( ,) : p d ) m 。一朋= 口n r ( ,) ( 此处仍设o , ) m ) ,可见期望法比频率法精度高。而且期望法不一定要 求,( z ) m 期望法的具体计算步骤为: ( 1 ) 独立地产生n 个c ,( o ,1 ) 随机数“1 , ( 2 ) i 一算z 。= 口+ ( 6 一) “。和“z i ) ,i = l ,n ( 3 ) 用( 1 1 ) 式估计, 频率法与期望法都可以估计,的值,但:者的精确度却不相同,这促使人们寻找 更好的m o n t ec a r l 0 方法。上述方法是使用了均匀随机数,实际上从( n ,6 ) 上取值的任 意随机数出发,都可以得到的估计量。从理论上看 o r ( j ) = :陋( 器) 2 一,2 ,若,( 。) o ,取9 扛) = 掣则有 n r ( j ) = o 。但由于,未知,这个方法实际上不能真正用于计算, 但是由此我们得到一个启示,即9 ( z ) 与,( z ) 越近似,估计的精确度越高。这就是重要 抽样法的基本思想。在实践中往往灵活地寻找常用的已知类型的密度g ,使它在峰值 附近与l ,( z ) i 较接近,以便达到降低估计方差的目的。从直观上看( 1 1 ) 式给出的j 是 采用均匀抽样,诸q 对j 的贡献是1 i 同的,( 戤) 大则贡献大。但在抽样时这种差别未 能体现出来,因此精确度小会很高。要达到同样的精确度需要的抽样次数多。而重要抽 样法则希望贡献率大的随机数出现的概率大,贡献率小的随机数出现的概率小,从而 提高精确度利用贡献率大小来降低方差的方法还有分层抽样法,它是把样本空间d 分 成一些小区间d l i d 。且诸d ,不相交,u d = d 。然后在各小区间内的抽样数由其 贡献大小决定,对旬日献大的眈抽样多,可提高抽样效率。另外提高抽样的效率的方 法还有关联抽样法、修正的重要抽样法等,不一一赘述。 第二章m 盯k a vc h a i nm o n t ec a r l o 方法 1 1m c m c 方法简介及马氏链基础知识 “ 。 宰 = 措 。 1 一n = 5 m o n t ec a r l o 方法的一个基本步骤是产生( 伪) 随机数,使之服从一个概率分布”( x ) 。 当x 是一维的情况时,这很容易做到。现在有许多计算机软件都可以得到这样的随机 数,前面介绍的例子都是这种简单情况。但x 常取值于舻,直接产生符合”的独立 样本通常是不可行的。往往是要么样本不独立,要么不符合7 r ,或者二者都有。以前 有很多人设计出许多方法去克服这一点。目前最常用的是m c m c 方法。m e t m p 0 1 i s 方 法( m e t r o p o l i s ,e t a 1 1 9 5 3 ) 与h t i n g s 的概括莫定了m c m c 方法的基石。此方法就是在 以万为平稳分布的马氏链上产生相互依赖的样本。换句话说,m c m c 方法本质上是一 个m o n t ec a r l o 综合程序,它的随机样本的产生与一条马氏链有关。基于条件分布的迭 代取样是另外一种重要的m c m c 方法,其中最著名的特殊情况就是g i b b 8 抽样,现在 已成为统计计算的标准工具,它晟吸引人的特征是其潜在的马氏链是通过分解一系列 条件分布建立起来。本文主要介绍这两种方法。下面先介绍一下有关马氏链( 离散) 的 基本知识 f :t 0 为一随机序列,随机序列所有可能取道的值组成的集合记为s ,称为状 态空间。如果对于0 及任意状态岛,s ,8 0 ,8 n l 都有 p ( 。= 勺i 五= 即,五一l = 8 “。= ) = p ( 五十1 = l 置= s j ) ( 1 3 ) 则称f 托:0 为一马氏链直观上看对于一个马氏过程,要预测将来的唯一有用 的信息就足过程当前的状态,而与以前的状态无关马氏链的性质完全由它的转移概 率( 转移核) 来决定,它足由状态以到状态s ,的一步转移概率,这里用p ( i ,) 表示。 即 p ( i ,j ) = p ( i ,j ) = p ( 五+ l = 彤i k = 8 )( 1 4 ) 用( t ) = p ( x = s ,) 表示马氏链在t 时刻处于状态s ,的概率,用7 r ( ) 表示在步状 态空间概率的行向量,初始向量用7 r ( 0 ) 表示,通常是 ( 0 ) 中只有一个分量为1 ,其余全 是o ,意味着马氏链从一个特定的状态开始,随着时间的变化,概率值扩散到整个状态 空间。在时刻+ 1 马氏链有状态值的概率由c h a p m a n k o l o m o g r o v 方程给出 几( + 1 ) = p ( x + _ 8 i ) = p ( 五+ l = s i 氍= s ) p ( 托一钆) = p ( 一t ) i ( t ) 一p ( 七,i kki ( 1 5 ) c h 印m a n - k o l o m o g r w 方程的连续迭代描述了马氏链的发展 如果我们定义一个矩阵p ,使它的第 行第j 列的元素为p ( i ,j ) 则p 叫做转移概率 矩阵,尸( t ,j ) 就是从状态z 转移到的概率p “一j ) ,这就意味着此矩阵每一行元索的和 为1 ,( ,p ( i ,j ) = ,p 一j ) = 1 ) 这样c h 印m a n - k o l o m o g r o v 方程可写成吏紧凑的 矩阵形式: 6 丌 + 1 ) = 丌( t ) p 利用矩阵形式,反复运用c h 印哪一k o l o m o g r o v 方程可得到 7 r ( t ) = 7 r 一1 ) p = ( 丌0 2 ) 尸) p = p 2 连续应用这种方式得到。 7 r ( ) = r ( 0 ) p f 定义n 步转移概率背,表示从状态i 出发经过n 步至n 状态j 的概率,即 p 9 专p ( x f + 乱= s ,i 五= 巩) ( 1 6 ) ( 1 7 ) 易看出蹬是矩阵p “的t 行j 列元素。 如果对所有的i ,j 都有正整数n 巧,使的砖,则称此马氏链为不可约的,也就 是所有状态互相之间都是有联系的,从一个状态总能转移到另一状态( 虽然可能不止 一步) 。如果从状态z 到状态”的步数都不是某个整数的倍数( 即无功约数) ,就称此链 是非周期的,换句话说,马氏链不会在某些状态之间以固定的长度陷入循环。如果马 氏链从状态t 出发能以概率1 地最终回到i ,则称状态 为常返态,所有状态全是常返态 的链称为常返链。若状态i 的平均回访时刻不是o 。,则称为正常返的,正常返的非周 期状态称为遍历态。 例假设状态空间是( 雨,晴,阴) ,并且天气变化过程是一个马氏过程,这样,明 天天气的概率仅依赖今天的天气,与以前的无关。假设已知今天是雨的情况下的转移 概率是 p ( 雨i 雨) = o 5 ,p ( 晴l 雨) = o 2 5 ,p ( 阴i 雨) = 0 2 5 这样转移概率矩阵的第一行为( 0 5 ,o ,2 5 ,0 2 5 ) 如此假设转移概率矩阵为 7 r ( 2 ) = 7 r ( o ) p 2 = ( o 3 7 5 ,o 2 5 ,o 3 7 5 ) 丌( 7 ) = ( o ) p 7 = ( o ,4 ,o 2 ,o 4 ) 如果假设今天是雨,即”( o ) = ( 1 ,o ,o ) ,则有 ( 1 8 ) ( 1 9 ) 7 7 r ( 2 ) = ( o 4 3 7 5 ,o 1 8 7 5 ,o 3 7 5 ) 和7 r ( 7 ) = ( o 4 ,o 2 ,o 4 ) 实际上,经过一个充分长的时间。期望的天气状况与初始状态无关了,换句话说, 马氏链达到了一个平稳状态,这儿的概率值与初始值是独立的。正如上面的例子所说 明的,一个马氏链可能达到一个平稳分布矿,这儿做为某个特别给定状态的概率向量 是与初始值无关的,平稳分布就是满足下式的分布 霄= 7 r 。p ( 2 0 ) 平稳分布的条件是马氏链是不可约非周期的,当马氏链是周期的,它会在状态之 间以一个确定的方式循环。是平稳分布的一个充分条件是下面的细致平衡条件成立 。 p 0 ,) 哼= p ( ,j ) ( 对所有i 和j ) ( 2 1 ) 如果对所有 ,都有方程( 2 1 ) 成立,就说马氏链是可逆的,因此( 2 1 ) 式有称可逆 条件。注意这个条件蕴含着”= 7 r p ,记( 7 r p ) ,为”p 的第j 个元素则 ( ”p ) ,= 死尸( t ,j ) = p o ,1 ) = 乃p o ) = 丌j ( 2 2 ) 最后一个等式利用了转移矩阵每行元素的和为l 离散状态的马氏链的基本思想可 以推广到连续情形,设其转移核为p ( z ,) 满足 厂p ( z ,) d ”= l ( 2 3 ) 则连续的c h a p m a i l - k o l o m o g r o v 方程变成 删= 厂( 妒( 刚) 妇 ( 2 4 ) 相应地平稳分布满足 。 州) = 厂州妒( 训) 白 ( 2 5 ) 2 1m e t r o p l 珏h t i n g 方法 前面已提到m c m c 方法就是通过建立一个平稳分布为7 r ( z ) 的马氏链来得到7 r ( z ) 的样本,基于这些样本做各种统计推断,概括起来分为以下三步: ( 1 ) 在d ,p 选一个”合适”的马氏链,使其转移核为尸( ,) ,”合适”的含义就是指7 r ( z ) 是其相应的平稳分布 8 ( 2 ) 由d 上的某一点x ( o ) 出发,用( 1 ) 中的马氏链产生点序列x ( ”,x ( 哪 ( 3 ) 对某个m 和足够大的n ,任一函数,和) 的期单估计如下 岛,:熹妻肛) ”一“t z 1 ( 2 6 ) 可以看出采用m c m c 方法时,构造转移核是至关重要的,不同的m c m c 方法往往 也就是转移核的构造方法不同。其中出现比较早的一种方法是m e t r o p l i 6 等人在1 9 5 3 年 提出的,h t i n g s 随后对其加以推广,叫做m e t r o p u 8 - h t i i l g s 方法,其做法如下:比如 说我们想建立一个以”( z ) 为平稳分布的马氏链,现任选一个不可约转移概率g ( ,) 以 及一个函数a ( ,) ,( o a 1 ) ,对任一组合( z ,g ) ,0 ) 定义: p ( z ,) = p ( z + ) = g ( z ,可) o ( o ,)( 2 7 ) 直观上看就是:如果链在时刻处于状态z ,即x ( ) = z ,首先南“jz ) 产生一个 潜在的转移z 一口,然后根据概率口( z ,) 来决定是否转移,也就是说在潜在转移点 找到后,以概率n ( z ,) 接受f 做为链在下一时刻的状态值,而以概率l o ( z ,) 拒绝 转移到,从而链在下一时刻仍处于状态z 。于是,在有了后,我们可从 o ,1 】上均匀 分布抽取。一个随机数。则 分布q ( i z ) 称为建议分布因为我们的目标是使7 r ( 动为平稳分布,因此在有了q ( - ,) 后,应选择一个a ( ,) 使相应的p ( z ,) 以7 r ( z ) 为其平稳分布,一个常用的选择( h a s t i n g ,1 9 7 0 ) 是: 此时,p ( 。,) 为: 出m 一叩,糕端) ( 2 8 ) i 口( 茁,可) ,7 r ( 可) q ( ,z ) 27 r ( z ) 口( z ,) p p 硝卜i 咖卅碧,咖舭) 7 r ( ) q ( ,z ) ,这种情况下 口( z ,) = 署薯器舞,n ( 玑z ) = 1 凶止b ( z ,咖( 司= 口( z ,) o ( ,) ”( z ) = g ( 为) 簧搿舞7 r ( z ) = 口( 掣,z ) 丌( 箩) = 口( ,z ) n ( 掣,z ) 丌( 可) = p ( 耖,z ) 丌( ) 3 、7 r 0 ) 口( o ,) ”( f ) g ( 玑z ) ,这种情况下 n ( z ,) = 1 ,o ( 玑z ) = 器舞籍 因此,p ( ,咖( ) = g ( ,z ) o ( ,。) 7 r ( y ) = 口( ,z ) ( 器辫崭) 7 r ( 口) = q ( 。,可) 丌( 卫) = 口( z ,暑,) o ( z ,掣) 7 r ( z ) = p ( z ,) 7 r ( $ ) 建议分布q ( z ,) 可取各种形式,下面是两种常用的选择:( 1 ) 对称的建议分布 。, ( m e t r o p o l i 8 ,1 9 5 3 ) ,即 q ( z ,可) = = q ( ,z ) ,比,可 此时o ( 为) 简化为n ,) = m 伽 1 ,籍) ( 2 ) 独立抽样 如果q ( 。,) 与当前状态z 无关,即口( 。,) = g ( ) ,则称此方法为独立抽样方法,此 时n ( z ,) 变为 n ( z ,9 ) = m 锄 1 ,器 其中u ) = 勰 一般,独立的效果可能很好,也可能很小好。通常,要使独立抽样有好的效果,口缸) 应 接近7 r ( z ) 。比较安伞的办法是使g ( z ) 的尾比7 r ( z ) 重。 例考虑逆x 2 分布7 r ( 8 ) = c p 一2 e z p ( 一蠢) 假设我们希望抽取符合7 r ( p ) 的样本,这儿设n = 5 ,n = 4 ,应用m e t r o p l i s 方法, 假设我们用( 0 ,1 0 0 ) 上的均匀分布做为建议分布,当然大于1 0 0 的概率是有的,但假 设此概率很小,以至于可以忽略。取= l 做为初始值,假设得到一个候选值( 来 自u ( 0 ,1 0 0 ) ) ,伊= 3 9 8 2 这样 n = m 饥( 端) ,1 ) = m 伽( 等萼鲁,1 ) = o 0 0 0 7 因为这个o 1 ,矿以0 0 0 0 7 的概率被接受,这样在( 0 ,1 ) 上均匀地抽取一个t ,如 果a 则接受矿( 把它做为下l 次抽样的初值,继续以上运算) ,否则就拒绝它,在建议 分布中抽取另一个候选值,继续以上运算。如此等等。下面图( 2 ) 是由得到的前5 0 0 个口 的值所画出的图像。下面图( 3 ) 是当建议分布改变为c ,( 0 ,1 0 ) 时得到的5 0 0 个口的值所 1 0 图2 : 画出的图像。对两组数值分别作如下处理:每2 0 个值为一组共分成2 5 组分别算出其均 值。发现剑最后面5 组所得数值已相差不大。由图( 2 ) 看出口取得l o 以上的值的情况还 是有的,但取得5 0 以上的值的情况己很少。所以把建议分布取为u ( o ,1 0 0 ) 已是比较保 险的了但足在迭代次数相同的情况下用矿( o ,l o o ) 比用( 厂( o ,1 0 ) 所用的计算时问多。 接受函数的选择 有这样一个问题,在m e t r o p o l i s 算法中哪,一种接受函数较好,比如说下面两种形式 出棚m 饥( 1 1 罴) ( 2 9 ) a ( 训) 2 揣 ( 3 0 ) 哪一个更好,对于分量只有两种状态的系统,后者的选择相当于g i b b s 取样( 此时7 r 仕) = 0 ,o ( 茁,口) = 1 ) ,这个问题看起来已被p e s k u n ( 1 9 7 3 ) 解决,他指i i 对任何固定的建议分 图3 l l 1 2 布( 2 9 ) 式的选择在以大类函数中是最优的,这个结论应用于m e t r o p o h s 算法。这些算 法能根据一个单独的建议分布表示出来。它意味着可改变状态的任何部分。这包括所 谓的”整体”m e t r o p o l i s 算法,也就是候选状态一般与所有分置都不同。也包括”局部” m e t r o d 0 l i s 算法,也就是被更新的分晕是被随机选择的。因为这些算法能根据一个建 议分布表述。p e 8 k l l i l 所用的判断标准是由马氏链产生的基于状态的一些函数的期望 的m o n t ec a r l 0 估计的”渐进方差”,此”渐进方差”的定义是n 次迭代后的估计的方差, 当n 趋于无限时的极限。p e s k u n 的结果还意味着最好用这样的一个建议分布,它对于 当前状态是零概率。 p k u n 的结果并没有真正解决这个问题,有两个原因。一是只有渴望得到非常精 确的估计的时候,对评价使用渐进方差才是一个适当的标准,而经常不是这种情况。 在许多问题中,即使在平稳分布中抽取一个单独的状态也将会有许多信息,并且在结 果中的许多不确定性应归于怀疑是否马氏链真正收敛与平稳分布,而不是在估计中运 用有限个样本产生的方差。可能更重要的是,这个结果没有应用到最普通的方法一分 量更新以固定顺序的局部m e t r o p o l i s 算法。确实,在这儿,伴随着一个简易分布的标 准p e 8 k u n 对当前状态是零概率能产生一个甚至是不遍历的马氏链。考虑一个二元变量 系统,所有状态都有相同特概率,假设转移函数且k ,被系统地应用,乳以概率l 给候选状 态,这儿第个变量的值是被更新的。因为所有状态有相等的概率,所有候选状态当 用m e t r 0 d o l i s 接受函数时将都被接受。一个单独的轨道将改变所有变量的值,下一个 轨道将把它们又都改变成原来的值,别的状态根本不会出现。一个另外的问题是在一 个j 维二元变量的系统中,除去( o ,1 ,o ) ,( 1 ,o ,1 ) 之外的所有状态都有相同的概率,它 的状态时任意的。对问题的一个现实的例子见( f l i e d b e r ga n dc m e r o n1 9 7 0 ) 对 大多数问题,利用系统更新的m e t r o p o l i s 接受函数产生一个遍历链,但是渐进方差要 比用( 3 0 ) 式的接收函数的同样方法的要大。在,s 伽g 系统中这种以经验为主的发现 是由c u n n i n g h a m ,m e u e r 和p e s k u n 得出的。这些问题的更多的理论结果是由n i g e s s id i s l e o ,h w a l l gs h e n ( 1 9 9 3 ) 得出的。虽然在大多数情况下,通常认为( 2 9 ) 式给出的接 受函数比较好,但仍有一些问题不太清楚。 2 3g i b b s 抽样方法( g e m a na n dg e m a n1 9 8 4 ) 前面已提到g i b b s 抽样方法是一种基于条件分布的迭代取样方法。这就要用到条 件分布,特别是满条件分布。所谓满条件分布就是形如”( 研i 。一r ) 的条件分布,其 中z r = 甄,i t ,z r = 钆,igt ) ,rc = 1 ,m ) 注意到,在卜述的条件分 布中,所有的变量全部出现了( 或 h 现在条件中,或出现在变元中) 。实际上g i b b s 抽 样方法可看做m e t m p o l 缸h a s t i n g 算法的特例。考虑置ix _ i i = l ,m 的条件分布, 1 3 选择一个转移核q ( 孔一iz 一) 固定翟= x _ = z 一不变。由吼( 甄一z :iz t ) 产生一 个可能的z :,然后以概率啦( 戤一iz 1 ) = m 锄( 1 ,磊等筹l 三i 訾等) 决定是否接受z 傲 为链的下一状态。这就是单元素m e t r o p o l 珏h a 8 t i n g 算法。若在上述算法中取q ( z 一 z ) 为”( 文iz i ) ,易得此时口0 一z ) = 1 这就是g i b b s 抽样方法。由于g i b b 8 抽样只涉及 单变量抽样,这使之最具吸引力。设给定了一个m 维联合分布7 r ( z l ,z 。) ,在g i b b s 抽 样中构造如下的转移核 只,垒p 0 ,) = n 墨1 口( 挑iy 1 ,挑一l ,茹k + 1 z 。) 其中z = ( 却z m ) ,静= ( 们,弘m ) ,戤d ,挑d 而”( 弧l 1 ,舭一“。 + 1 ,z 。) 是在除第k 个分量外,将第1 至k 1 个分量固定为l ,讥一1 ,并 将第k + 1 至第m 个分量固定为。k + 。,名。的条件下,第k 个分量在批处的条件分布。可 验证( b ,) 确实是一个转移阵,即。,= 1 ,”( i ) = 。絮案譬蠢手= 署券筹= 1 ”p z ,= 跏撕一i 【p 17 r ( li 2 ,z m ) 】丌( 啦i 可l ,z 3 z m ) 7 r ( 可ml l ,掣讯_ 1 ) = 、 1 可验证”( z 1 ,z 。) 是以p ( z ,) 为转移阵的马氏链的平稳分布 ( ) 7 r ( z 1 7 一) p ( ( z l i 。) ,( y l j ) ) 岛。m 一1 。,7 r ( z l ,$ m ) 7 r ( l z 2 ,一z 仇) 丌( 0 2l 掣l ,幻,z m ) 7 r ( 3 mi i ,一! h 。一i ) 。”( 轧z m ) 碧箫鬻踹”( z z 。) ”协一- ) = 。:”( l ,奶,z 。) 7 r ( 2ly l ,z 3 ,。) 7 r ( 跏i 1 ,一1 ) = = 。7 r ( “z m 瓦端篇鲁雨= 口( 轧) g i b b s 抽样具体步骤如下: 由马氏链p ) 的样本可按以下程序得到x h - p ) 的一个样本( 轧鼽) 。 ( 1 ) 先得到服从分布 7 r ( 1z 2 ,z 。) ,y l d 1 的随机变量j 乙“1 ( u ) ,( 注 意z 2 ,z 。来自k ( u ) ) ( 2 ) 在得到服从分布 7 r ( 珈i 。,3 ,z 。) ,耽d 2 的随机变量墨l + 1 2 ) 的一 个样本耽,依次下去,得到服从分布 7 r ( 玑i l ,讥一l ,z k + l ,z 。) ,批d 的随机变 量弱+ l ,i p ) 的一个样本批= 1 ,一m 一1 ) 最后得到服从分布 7 r ( if l ,一1 ) ,d 的随机变量x 。+ l ,。p ) 的一个样 本 定义垒( 妣蜘) 就是x ;+ l ( “j ) 的一个样本。 现任取一个初值j 而0 ) = ,( o ) ,按上面方法得x l ) 的一个样本g ( 1 ) ,对n 归 纳地得到( u ) ,j 厶) 的样本( ”,( ”) ,当n 充分大时,马氏链k ( u ) 分布近似 1 4 于7 r ( z 。) ,就可以认为( 曲是近似服从7 r ( z l ,) 的一个样本。 简单地说g i b b s 抽样就是在随机变量x = ( z l ,z 。) 中选一个坐标,比如z ;。然后 用一个新样本z :去更新它,由条件分布”( f 。卜吲) 抽取,这个等待更新的以如何选取 方法并不唯一。上面介绍的方法只是其中一种,即所谓的系统扫描。它是按下标顺序 从( 1 ,m ) 中逐个选取。如果是按照一个给定的概率向量,比如说,( 去,击) 随机地选 一个 ,这就是随机扫描,有研究者说明在马氏链收敛速度方面随机扫描胜过系统扫描。 显然,随机扫描的g i b b 8 抽样中每一条件更新步骤7 r 都是不变的。假设x ( 2 ) 一7 r ,x 恐、是 在7 r 下的边际分布,那么 ,r ( z p ”ix 陷7 r ( x 出】) = 7 r ( x 凶,z p ”) 而x 岛”= x 陷 这就意味着更新后,新结构仍服从分布7 r 当空间是离散的,l i u f l 9 9 6 c ) 构建了一种”o v e r - r e l a t i o n 策略来提高g i b b s 抽样的 效率。令甄有m ;个可能的值,是感兴趣的分布。在一般的随机扫描g i b b s 抽样中,一 个下标 是随机的第+ 一个选择,当前值甄被对应的满条件分布中抽取的值鼽所替代。这 儿考虑这个过程的一个修正,就是鼽以概率m i n 1 , 三! 捌 代替否则,施保持不 变,l i u ( 1 9 9 6 c ) 证明了这个修正策略比随机扫描的g i b b s 抽样更有效。 回顾g i b b s 抽样,还可以从以下角度去认识它:每一一g i b b 8 更新可看做对当前状态 根据r 个选择的随机再布置。例如,如果第个坐标被选择,那么这个”再布置”可表示 为( z l ,“,z 。) 一( z l ,+ r 靠) ,这儿r 是在一个适当的分布中随机抽取,易 说明如果r 在2 ) ( r ) 。( ”( z + r ,x 一m ) 中抽取,则7 r 在这种更新下是不变的。这种看法对 设计更有效的m c m c 抽样方法是很有用的。 例假设z ,的联合分布如下, p ( z ,) = t i :! 戋丽旷+ “1 ( 1 一) ”一。+ 口一1 ,其中z = o ,1 ,n ,o 1 注意到这儿z 是离散的,但是连续的,它们的联合密度很复杂,但条件密度很简 单。同顾二项随机变量z 有密度 p ( zl 口,n ) o ( 穹i ! 等,( o z 1 ) 这儿o q t 时马氏链收敛到平稳分布。我们一般试图通过 马氏链的一个特殊的实现去判断是否收敛,而状态一般是高维的,因此难以直接估计。 大多数判断收敛的是基于监测一个或一些状态的数量函数,比如这样的醋数,它的期 图7 5 01 0 0 望是我们在估计中感兴趣的。最通常的方法是随着时间的增加,简单地画出这些函数, 然后在直观上判断什么时候达到平稳。 上面给出了这种状态函数的三种假设的图形,对应着三种不同的马氏链的实现。 图7 说明经过大约2 5 次迭代后收敛。但是若没有进一步理论分析,不能保证这是真正的 收敛。图8 就说明这种町能性。如果每一轨道只模拟迭代5 0 次,可能就认为图8 表示的轨 道在大约2 5 次迭代后就收敛了,就可以用2 5 次到5 0 次之间的平均值估计期单了。这样 做的话就会出错,因为在得到1 0 0 次迭代的数据后就发现直到7 5 次迭代才收敛。图9 说 明了另外一种情况,由模拟数据无任何理由说明轨道已收敛( 注意我们同样不能确定 图8 表示的轨道不受敛,每一轨道都可能是从真实分布中抽样,只是状奎之间的转换很 慢,当只进行1 0 0 次迭代时,它们看上去是来自不同分布的样本) 。如果仅从图9 的一条 轨道来看可能认为5 0 次迭代甚至更少后就收敛了。这说明当有明显的证据说明从一条 轨道判研收敛小h j 靠时,多条轨道的作用。在i 圣1 8 中如果单从一条轨道也小能发现5 0 次 迭代后的变化,只能从多条轨道才能发现收敛是在7 5 次迭代以后( 或者说,用多条轨 道是一种比较保险的方法) 有人建议通过计算实现对时间的函数值的平均值,寻找那些平均值己稳定的点来 判断收敛。这种方法对图7 ,图8 两种情况很有效,但对图9 这种方法显然行不通。对于 平稳分布是多峰的,就会出现图9 这种情况,有人己对这种情况作过分析。 如果我们正在监控状态的两个函数n ( z ) 6 ( z ) ,可能h 现从n ( z ) 看是收敛了,而 从6 ( z ) 看还没达到平稳。这时我们可能用实现的值去估计e 【o ( z ) 】。然而,除非我们理 解为什么a ( z ) 收敛而6 ( z ) 没有收敛,否则这样做是不妥的。好象6 ( z ) 没有收敛,仍无 b 2 5如 1 0 0 图8 2 5 5 07 51 0 0 图9 把握确定访问的状态是

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论