（概率论与数理统计专业论文）基于gibbs抽样和em算法的生物保守序列motif识别.pdf

基于g i b b s 抽样和e m 算法的m o t i f 识_ 斟 摘要 基于g i b b s 抽样和e m 算法的生物保守序列m o t i f 识别 中文摘要 随着许多物种的基因测序工程的完成和生物技术的发展，人类拥有了大 量的生物数据本世纪一个具有挑战性的问题就是挖掘这些数据中的生物 信息，其中发现生物序列中的保守片断是一个重要的问题这些保守片断 被称为m o t i f g i b b s 抽样是在生物序列中的m o t i f 识别中应用最广泛，最成功的算法 以往的研究都把m o t i f 的长度视为固定的，而实际情况是事先并不知道m o t i f 的长度本文通过把m o t i f 的长度看作缺失数据，通过算法来确定这个长 度实验结果表明，这个算法是可行的 b a i l e y 和e l k a n 在1 9 9 4 年通过二元混合模型把e m 算法用于生物序列中 的m o t i f 识别这个方法首先把原来的生物序列截断，然后用二元混合模型 来拟合新的数据集注意到新的数据集中有很多数据并不能由这个二元混 合模型来生成本文通过引入多元混合模型来拟合这个数据集，从而使每 个数据都能由这个多元混合模型生成由于我们的模型能更准确地描述数 据，从而能够使参数更快更准确地收敛于真正的参数值 关键词：m o t i f , g i b b s 抽样，e m 算法 作者：陈晓林 指导老师：汪四水( 副教授) 基干g i b b s 抽样和e m 算法的m o t i f 识别 a b s t r a c t m o t i fd i s c o v e r yb a s e do nt h eg i b b ss a m p l i n ga n d e m a l g o r i t h m a b s t r a c t w i t ht h ec o m p l e t i o no fg e n o m e ss e q u e n c i n go fm a n y8 p e d 曙a n dt h ea d v a n c eo f m i c r o a r r a yt e c h n o l o g y , w eb e g i nt op o s s e s saw e a l t ho fv a l u a b l eb i o l o g i c a ld a t a ，o n eo f t h em o s tc h a l l e n g i n gp r o b l e m so ft h i sc e n t u r yi st od e c i p h e rt h ei n f o r m a t i o ne m c o d e d i nt h e s ed a t a ，i ti sa l li m p o r t a n tp r o b l e mt od 埴螂t h ec o n s e r v e d 圈l b s e q u a 醒j 篦i n b i o l o g i c a ls e q u e n c e s t h e s ec o n s e r v e d8 1 1 b s e q u e n c e sa r ec a l l e dm o t i f s g i b b ss a m p l i n gi st h em o s ts u c c e s s f u lm e t h o do fd i s c o v e r yo fm o t i fa n da p p l i e d m o s tb r o a d l y t h el e n 时ho ft h em o t i fi sl o o k e d c o n s t a n ti nt h ep r e v i o u sl i t e r a - t u r e b u ti nf a c tt h el g t hc a nn o tb ek n o w ni na d v a n c e i nt h i sa t t i c l et h el 印g t hb t r e a t e d m i s s i n gd a t aa n dd e t e r m i n e db yt h ea l g o r i t h m 。a n dt h er e s u l t so f e x p e d m e n t s h o wt h a to u ta l g o r i t h mi sf e a s i b l e b a r l e ya n de 虬c a ni n t r o d u c e dt h ee ma l g o r i t h mt h r o u g ht w o - c o m p o n e n t 吐】( t l l 弛 m o d a li n t ot h ef i e l do f d i s c o v e r yo f m o t i f i nb i o l o g i c a ls e q u e n c e sa t1 9 9 4 ，i nt h e i rm e t h o d ， t h es e q u e n c e sa r ec u td o w nt os h o r t8 1 1 6 6 e q u e l l 懒a tf i r s t t h e l lt h et w o - e n m p o n e n tm i x - t u r em o d e li su s e dt of i tt h en e wd a t u s e tu s i n ge ma l g o r i t h m n o t i c et h a tm a n yd a t ai n t h ed e rd a t a s e tc a nn o tb eg e n e r a t e db yt h et w o - c o m p o n e n tm i x t u r em o d e l h e r e 骶 u s eam o r e - c o m p o n e n tm i x t l l r em o d e lt of i tt h en e wd a t a s e t e v e r yd a t ac a nh eg e n - e r a t e db yt h en e wm o d e l b e c a u s eo i e rm o d e lc a l ld e s c r i b et h ed a t am o r ep r e c i s e l y , t h e p a r a m e t e r sc a nc o n v e r g et ot h er e a lv a l u em o r eq u i c k l ya n dp r e c i s e l y k e y w o r d s ：m o t i f , g i b b 8s a m p l i n g ，e ma l g o r i t h m w r i t t e nb yc h e nx i a o l i n s u p e r v i s e db ya s s o c i a t ep r o f w a n gs i s h u i i i 苏州大学学位论文独创性声明及使用授权声明 学位论文独创性声明 本人郑重声明：所提交的学位论文是本人在导师的指导下，独立进行 研究工作所取得的成果除文中已经注明引用的内容外，本论文不含其他 个人或集体已经发表或撰写过的研究成果，也不含为获得苏州大学或其它 教育机构的学位证书而使用过的材料对本文的研究作出重要贡献的个人 和集体，均已在文中以明确方式标明本人承担本声明的法律责任 研究生签名：睦啦捶日期：卑：旦堑盏 学位论文使用授权声明 苏州大学、中国科学技术信息研究所、国家图书馆、清华大学论文合 作部、中国社科院文献信息情报中心有权保留本人所送交学位论文的复印 件和电子文档，可以采用影印、缩印或其他复制手段保存论文本人电子 文档的内容和纸质论文的内容相一致，除在保密期内的保密论文外，允许 论文被查阅和借阅，可以公布( 包括刊登) 论文的全部或部分内容论文的 公布( 包括刊登) 授权苏州大学学位办办理 研究生签名：日期： 导师签名：速婴基日期： 第一章引言 1 1生物背景知识介绍 1 1 1 d n a 和r n a 生物体的遗传物质是核酸核酸分为脱氧核糖核酸( d e o x y r i b o n u l e i ca c i d ， d n a ) 和核糖核酸( r i b o n u l e i ca c i d ，r n a ) 这两类核酸分布于生物体的每一个 细胞中，其中d n a 是主要的遗传物质 核酸的基本结构单位是核苷酸，其组成方式为碱基戊糖一磷酸 碱基包括嘌呤和嘧啶两类d n a 中的嘧啶为胞嘧啶( c y t o s i n e ，c ) 和胸 腺嘧啶( t h y m i n e ，t ) ；嘌呤为腺嘌呤( a d e n i n e ，a ) 和鸟嘌呤( g u a n i n e ，g ) r n a 和 d n a 的不同在于，r n a 没有胸腺嘧啶而有尿嘧啶( u r a c i l ，u ) 戊糖又分为肛核糖和n 二脱氧核糖两类，前者存在于r n a 中，后者存 在于d n a 中 d n a 和r n a 的组成见表1 , 1 表1 1 d n a 和r n a 的化学组成 成分 d n ar n a 腺嘌呤a腺嘌呤a 碱基 鸟嘌呤g鸟嘌呤g 胞嘧啶c胞嘧啶c 胸腺嘧啶t尿嘧啶u 戊糖脱氧核糖核糖 磷酸磷酸磷酸 d n a 是由核苷酸组成的线性序列因为d n a 的核苷酸根据碱基类型的 不同可以分为4 类，所以我们可以简单地用a ，c ，g ，t 来表示4 种核苷酸 w a t s o n 和c r i c k 在1 9 5 3 年提出了d n a 的双螺旋结构d n a 分子是由两条 反向平行的核苷酸链绕着同一中心缠绕而成，形成右手螺旋两条多核苷 酸链靠碱基间的氢键相连配对的碱基又称互补碱基对( b a * p a i r ) ，其中a 与t 配对，c 与g 配对因为d n a 分子的一条链可以认为是核苷酸的线 性序列，所以它可以用一条字符串表示，其中的字母来自字母表( a ，c ，g ，t 基于g i b b s 抽样和e m 算法的m o t i f 识别 第一幸引言 r n a 和d n a 不同，r n a 是单链分子r n a 主要有三种，信使r n a ( m r n a ) 核糖体r n a ( r r n a ) 及转运r n a ( t r n a ) 1 1 2蛋白质 蛋白质是生物体内占有特殊地位的大分子，它是生物体的基本构件，也 是生命活动的重要物质基础，几乎一切生命现象都要通过蛋白质的结构和 功能体现出来 氮基酸是蛋白质的基本结构单位，参加蛋白质组成的氨基酸有2 0 多种 各种氨基酸由肽键连接形成一条较长的氨基酸链这条氨基酸链就称为蛋 白质氨基酸在蛋白质中通常称为氨基酸残基( r e s i d u e ) 蛋白质在空间中有着更加复杂的结构，但是它的线形结构是我们的研究 目的 1 1 3分子生物学中心法贝口 d n a 是遗传物质，是携带遗传信息的载体信息从基因的核苷酸序列 中被提取出来。用来指导蛋白质合成的过程对地球上的所有生物都是一样 的，分子生物学家称之为中心法则生物体的遗传信息以密码形式编码在 d n a 分子上，表现为特定的核昔酸排列顺序，并通过d n a 的复制使遗传信 息从亲代传向子代在后代的生长发育过程中，d n a 分子中的遗传信息转 录到r n a 分子中，再由r n a 翻译生成体内的各种蛋白质，行使特定的生物 功能，翻译是在核糖体上进行的这样，通过遗传信息从亲代传向子代，并 在子代表达，使得子代获得亲代的遗传性状b n a 也能通过复制过程合成 与其自身相同的分子此外，生物界还存在由r n a 指导下的d n a 合成过 程，即逆转录这过程发生在逆转录病毒中通过基因转录和翻译得到的 蛋白质分子可以反过来作用于d n a ，调控其它基因的表达 1 1 4转录调控 分子生物学的中心法则表明某些d n a 片断转录成r n a ，r n a 又翻译形 成蛋白质在一个生物体中，任何细胞都带有相同的基因，但是一个基因 在不同的组织、不同的发展阶段和不同的环境中的表现并不一样这个不 同是由于在不同的细胞和不同阶段，合成的蛋白质是不同的如果一个阶 段某个蛋白质被合成，那么它的编码d n a ( 基因) 就称为表达的因此在一 个特殊的生理阶段的细胞可以看成一个特殊的机械系统，在这个系统中不 同的基因要么表达要么不表达 2 基于g i b b s 抽样和e m 算法的m o t i f 识删第一章；l 吉 在许多生物体内，编码蛋白质的d n a 只占全部d n a 的一小部分例 如。人类基因组中，基因只占整个基因组的百分之十五那些非编码d n a 原来被认为是没有用处的事实上，这些d n a 的非编码部分包含着蛋白质 合成的控制机制基因的大部分控制序列位于上游的调控区域，这个区域也 被称为转录调控区域或启动子基因的转录不仅需要转录调控区域的d n a 序列，还需要许多特殊的蛋白质，这些蛋白质被称为转录因子转录因子 结合到转录调控区域的特殊d n a 形式，从而控制基因的转录 1 2生物序列中的m o t i f 识别 1 2 1m o t i f 识别的意义 随着许多物种基因测序工程的完成我们开始拥有大量有价值的生物数 据但是这些数据远远不能被我们所利用因为这些数据是完全没有任何 语法规则和空格的字符串挖掘这些数据中的生命信息是后基因时代的一 项重要任务识别多个生物序列中的相似子序列具有非常重要的意义因 为这些子序列往往暗示了一些重要的功能和结构例如：转录因子可以结 合到转录调控上游区域的特殊d n a 片段从而控制基因的转录一个转录因 子可以和许多不同的上游区域结合，从而对基因进行调控相同转录因子 的结合位点( 即m o t i f ) 具有明显的保守性识别这些结合位点可以帮助我们 更好地理解基因转录调控 1 2 2m o t i f 的定义及描述 生物序列中的m o t i f 是一组生物序列中具有保守性的短的子序列，也就 是这组生物序列中相似的子序列这一段子序列可能是一个连续的块，也 可能是不连续的块，我们把这种不连续的m o t i f 称为带空格的m o t i f 现在，主要有两种方法来描述生物序列中的m o t i f 一致性序列方法和矩 阵方法下面我们以d n a 序列为例来说明这两种方法 一致性序列方法把d n a 序列中的m o t i f 看成通配符表( 表1 2 ) 中的字母 组成的字符串它描述了m o t i f 的每个位置可能的碱基类型例如，c a a t 转录因子结合的核酸序列m o t i f 表示为5 - g c c a a t c t - 3 ，热休克因子结合的 核酸序列m o t i f 表示为5 - c n n g a a n n t c c n n g - 3 一致性序列是关于生物序 列m o t i f 的一种定性描述它只能说明m o t i f 序列的每一个位置可能出现的 字母类型，但不能说明每种字母类型出现的可能性的大小而矩阵方法能 3 基于g i b b s 抽样和e m 算法的m o t i f 识别第一幸引言 很好的解决这个问题 d n a 序列m o t i f 的矩阵表示可以分为3 种：频数矩阵，它的每一个元素 吩i 记录了m o t i f 的第j 个位置是第k 个字母的频数；频率矩阵，它的每一个 元素肛= 等，其中n 是m o t i f 发生的总数；加权矩阵，它的每一列中的每 个元素都表示在j 当对中这一位置的某个字母出现的权重 表1 2 ：通配符表 符号含义说明 g g 鸟嘌呤 aa腺嘌呤 tt胸腺嘧啶 cc胞嘧啶 rgo ra嘌呤 yto rc 嘧啶 ma0 rc 氨基 kgo rt酮基 sgo rc强氢键 wao rt 弱氢键 hao r c o r t非g bgo r c o r t非a v go r co r a非t ( 非u ) dgo r to r a 非c ngo r co r t o r a任意碱基 1 2 3m o t i f 识别的方法 m o t i f 识别的方法可以分为两类：形式驱动的方法和序列驱动的方法形 式驱动的方法一般来说先对m o t i f 的形式给一个定义，然后检验序列中是否 存在这样的m o t i f 一个普遍的方法是把所有可能的m o t i f 的形式列举出来， 然后在序列中检查所有匹配的子串，对所有的m o t i f 形式打分，根据得分大 小排序，得分越高越有可能是我们要找的m o t i f 序列驱动的方法采取另一 种不同的方法这种方法直接根据序列的数据集来进行m o t i f 识别它利用 m o t i f 矩阵迭代更新的方法寻找生物序列中的m o t i le m 算法和g i b b s 抽样 就是序列驱动的方法 4 基于g i b l 抽样和e m 算法的m o t i f 识别 1 3m o t i f 识别中两个常用的方法 第一章引言 关于生物序列中的m o t i f 识别算法很多，这里我们介绍两个最著名的从 统计角度进行m o t i f 识别的算法 1 3 1 m e m e m e m e ( m u l t i p l ee mf o rm o t i fe l i c i t a t i o n ) 是由b a i l e y 和e l k a n ( b a i l e ya n d e l k a n ，1 9 9 4 ，1 9 9 5 ) 提出的m o t i f 识别算法它是对l a w r e n c e 和r e i l y ( l a w r e n c ea n d r e i l y , 1 9 9 0 ) 提出的基于e m 算法进行m o t i f 识别的方法的一种改进l a w r e n c e 和r e i l y 的方法只能发现一个m o t i f 类型而m e m e 通过用e m 算法拟合一 个二元素的混合模型来发现一个或多个m o t i f 类型这个算法首先把序列集 中的每一个序列截断成一些固定长度的子序列，把这些子序列看成一个新 的样本集合样本集中的每个子序列可能是m o t i f 序列，也可能不是这样 每一个子序列就可以看成是一个服从二元混合分布的随机变量的独立实现 值而二元混合分布中的一个分量是m o t i f 模型，另一个是非m o t i f 模型， 即背景模型这时，e m 算法被用来估计二元混合分布中的参数估计出 模型的参数后，就可以根据贝叶斯最优分类方法把每一个子序列分为m o t i f 序列或背景序列 与g i b b s 抽样方法( 1 3 2 ) 不同的是m e m e 是一种确定性策略，所以算法 有可能被困在参数的一个局部最优值所以选择多个初始值多次运行m e m e 是一个合理的方法它能帮助我们找到参数的全局最优值从而正确的确 定m o t i f 序列 1 3 2g i b b s 抽样 g i b b s 抽样算法是一种特殊的马尔可夫链蒙特卡罗( m a r k o vc h a i nm o n t e c a r l o ，m c m c ) 方法，该算法最早由l a w r e n c e ( l a w r e n c ee t a 1 1 9 9 3 ) 等引入蛋白 质序列中的m o t i f 识别g i b b s 抽样算法是应用最广泛最成功的m o t i f 识别算 法l i u 等( l i ue t a l1 9 9 5 ) 把g i b b s 抽样和贝叶斯方法结合起来用于生物序 列的m o t i f 识别，得到了很好的结果后来又有很多人对这个算法进行了改 进例如，r o t h ，h u g h e s ，e s t e p 和c h u r c h ( 1 9 9 8 ) ，l i u ，b r u t l a g 和l i u ( 2 0 0 1 ) 等这些 算法已经形成了成熟的软件，如：g i b b sm o t i fs a m p l e r ，a l i g n a c e ，b i o p r o s p e c t o r 等 g i b b s 抽样算法识别生物序列中m o t i f 的基本原理是通过随机抽样不断 更新m o t i f 模型和序列中m o t i f 出现的位置来最优化后验分布密度当满足 5 基于g i b b s 抽样和e m 算法的m o t i f 识别 第一章引言 一定的迭代终止条件，我们就得到了最优的m o t i f 位置和m o t i f 模型g i b b s 抽样算法和m e m e 算法的一个重要不同就是g i b b s 抽样算法采用的是随机 化的策略每次迭代g i b b e 抽样随机地选择m o t i f 出现的位置，因此能够保 证找到参数的全局最优值 6 第二章g i b b s 抽样和e m 算法 2 1马尔科夫链蒙特卡罗( m c m c ) 方法 2 1 1蒙特卡罗( m o n t ec a r l o ) 方法 m o n t ec a r l o 方法也就是随机模拟的方法它的起源和发展始于2 0 世纪 4 0 年代但从方法特征的角度来说，可以一直追溯到1 8 世纪后半叶的b u f f o n 随机投针试验，即著名的b u f f o n 问题m o n t ec a r l o 方法的基本思路可以分 为以下几步： 1 ：针对实际问题建立一个简单而且便于实现的概率统计模型，使所求 的解恰好是所建立的概率分布或模型的某个数字特征，比如，是某个事件 的概率或者该模型的期望值 2 ：对模型中的随机变量建立抽样方法，在计算机上进行模拟试验，对 关心的事件进行统计计算 3 ；对模拟试验结果加以分析，给出所求解的估计及其精度的估计 4 。必要时还应改进模型以降低估计的方差和减少试验的费用，提高模 拟计算的效率 m o n t ec a r l o 方法的应用范围非常广泛，它即能求解确定性问题，又能解 决随机性问题m o n t ec a r l o 方法近年来应用广泛我们下一节介绍的马尔 科夫链蒙特卡罗( m c m c ) 方法是近年来发展最快的统计分支 2 1 2马尔科夫链蒙特卡罗( m c m c ) 方法 m c m c 方法是一种基于马尔科夫链( m a r k o vc h a i n ) 的m o n t ec a r l o 方法 m c m c 方法最早出现在统计物理中由于它的实用性，最近十几年它成为 统计学研究的一个热点统计学者们不仅对m c m c 方法的统计理论进行了 大量的研究，而且把它广泛地应用到实际中例如，m c m c 方法在在生物 统计和图像处理中得到了广泛应用并取得了很大成功 为了阐明m c m c 方法的基本思想，我们首先介绍一下有关马尔科夫链 的知识 假定随机变量序列 一，x - ，妒， 满足马尔科夫性。即：在任意时刻 t ( o ) ，序列中下一时刻t + 1 处的x 蚪- 由条件分布p ( x l x t ) 产生，它只依 赖于时刻t 处的当前状态而与时刻t 以前的历史状态似o ，x - ，x 2 ，x 一 无关这样的一个随机变量序列就称为马尔科夫链 7 基于g i b b s 抽样和e m 算法的m o t i f 识别 第二幸g i b b s 抽样和e m 算法 一般地，令 x 。) t 郅为省上的马尔科夫链，其一步转移概率为 p ( x ，z ) = p ( z 一一) = p ( + 1 = z j 一= z ) ( 离散) ， p ( z ，聊= p ( z b ) = 上p ( 2 j 2 ) d ( 连续) ，( ，) 称为该马尔科夫链的转移核，转移核p ( x ，一) 代表处于该状态。下的过 程下一步转移到状态z ，的概率若p ( ，) 与时间t 无关，则称该马尔科夫链 是时间齐次的t 步转移概率函数为 p t ( x ，蛋) = p ( 一+ 1 = 一i z l = 习 相应的，矿( z ，z ) 代表于该状态z 下的过程经过t 步转移到状态一的概率记 x o 的分布为p 0 ) = p ( x o = 动，经过t 步后x 2 的边际分布为p ( = z ) = m 扛) 如果”( z ) 满足 ， fp ( x ，一) 口( z ) d x = ( z ) y x z ， j 则称”0 ) 为转移核p ( ，) 的不变分布假如经过t 步的边际分布为”( z ) ， 则马尔科夫链的t 步以后任意时刻都服从7 r ( z ) 的分布 m c m c 方法的基本思想是通过建立一个不变分布为_ l r ( z ) 的马尔科夫链 来得到”扛) 的样本，并基于这些样本作统计推断 2 1 3m e t r o p l i s - h a s t i n g s 算法 最早出现的m c m c 方法是m e t r o p l i s 等在1 9 5 3 年提出的m e t r o p l i s 算法后 来的所有m c m c 方法的思想都是基于这个方法在过去的5 0 多年里m e t r o - p l t s 算法已经广泛应用于统计物理中这一节简单介绍一下这个算法 假定马尔科夫链第t 步的状态为一，”( z ) 是我们的目标分布m e t r o p l i s 算法通过以下两步得到第t + 1 步的状态： 1 ：从建议分布丁k 一) 中生成，这里r ( z ，一) 是对称的概率转移函数， b pt ( z ，z ) = t ( 一，$ ) 2 ：生成随机数f ，使u 服从( 0 ，i ) 上的均匀分布如果us 毫 ，则令 ，+ l = 一，否则令x t + l = z 在m e t r o p l i s 算法中，建议分布t ( 毛一) 必须是对称的h a s t i n g s 在1 9 7 0 年把这个算法推广到非对称的情形但要求r ( z ，一) 满足：丁( 毛一) 0 当且 仅当t ( x ，z ) 0 8 基于g i b b s 抽样和e m 算法的m o t i f 识删 第二章g i b b s 抽样和e m 算法 假定马尔科夫链第t 步的状态为矿，”( z ) 是我们的目标分布h a s t i n g s 算 法通过以下两步得到第t + 1 步的状态： 1 ：从建议分布t ( x ，一) 中生成一 2 ：生成随机数u ，使u 服从( o ，1 ) 上的均匀分布令r = m i l ，瓣) ， 如果c ，r ，则令+ l = 一，否则令z 蚪1 = z 下面我们证明对于在上述转移核下生成的马尔科夫链，”( z ) 是其不变 分布 引理2 1 1 如果马尔秤夫链的转移核 ( 毛z 。) 满足条件a ( x ，一) ”( z ) = a f ，z ) ”( z ) ， 则* ( ) 是该马尔科夫链的不变分布 证明：由 a 扛，一) ”( 功= a ( 一，z ) 订( 一) ， 知 a ( z ，一) 丌( z ) d 】【= ( 一，功”扛) d ) 【= 丌( 一) a ( 一，z ) d x = w ( 一) 故结论成立 引理2 1 2 由m e 亡， o p l i s - h a s t i n g s 算法得到的马尔科夫链的转移核a ( x ，一) 满 足a ( z ，一) 7 r ( ? ) = ( 一，z ) ( 一) 证明：当z 一时， ( ) = t ( 叩) m i n 1 ，删r x st z t , x ) 故 ，如一t 一) 血 t ，薷骞磊智) = m i n r ( 。) t ( ，一) ，”( ) t ( 。，z ) ) _ 附枷t 磊蓦雾舞，- 当z = 一时，显然有a ( x ，一) 7 r ( z ) = a ( 一，z ) ”【z ) 结论得证 定理2 1 1 在上面的假设下，7 r ( z ) 是由m e t r o p l i s h a s t i n g s 算法得到的马尔 科夫链的不变分布 证明：由引理2 1 1 和引理2 1 2 立得 关于m e t r o p l i s - h a s t i n g s 算法，还有很多变形，像随机游走的m c 北r o p l i s - h a s t i n g s 算法，多点m e t r o p l i s - h a s t i n g s 算法等 9 基于g i b b s 抽样和e m 鼻法的m o t i f 识别第二章g i b b s 抽样和e m 鼻法 2 2缺失数据 缺失数据方法是在对复杂数据建模及设计计算方法时常用的方法缺失 数据的统计分析首先是在复杂的抽样调查分析时使用的，最近三十多年得 到了很大的发展缺失数据的个重要问题是缺失机制的问题如果缺失数 据的缺失只依赖于观测数据和某个未知的参数f 而与缺失数据本身无关， 则称缺失数据是随机缺失的如果缺失数据是随机缺失的，并且f 与模型参 数0 不同，则称缺失机制是可忽落略的在缺失机制可忽落略的情况下， 我们可以很容易地应用e m 算法和m c m c 方法进行统计推断一些复杂的 统计问题，看上去并不涉及缺失数据但是当我们引入一些变量并把它们 看成缺失数据时，我们可以很容易地对这些统计问题求解例如，混合模型 的参数估计本来是很难处理的，我们很难得到它们的最大似然估计，但是 如果我们引入一个变量来表示样本来自混合分布的哪一个分量，并把它们 看成缺失数据，我们能很容易地通过e m 算法得到参数的最大似然估计的 近似值 2 3g i b b s 抽样 g i b b s 抽样是最简单应用最广泛的m c m c 方法，它是由g e m a n 和g e m a n 于1 9 8 4 年提出的g i b b s 抽样的一个最显著的特点就是通过条件分布来构 造马尔科夫链的转移核这时候满条件分布就起了很重要的作用当满条 件分布得到后，剩下的问题就只是如何从满条件分布中抽样了 2 3 1满条件分布 设随机变量x 可以分解为x = 旺l ，) ，令五一 l = ( x j ：j a c ) ，x a = x j ：j 舢，其中a 是 1 ，d ) 的任子集注意到在条件分布( x a x l 一棚) 中，所有的变量全部出现( 或出现在条件中，或出现在变量中) ，所以称这种 条件分布为满条件分布 注意到对任意的z 和a ，满条件分布 7 r ( z i z 卜，卅) = 了j i ! ； ：霄( z ) 因此，在计算满条件分布时，只需保留”( z ) 表达式中与z 一有关的项下 面，我们通过一个例子来看一下满条件分布具体如何推导 1 0 基于g i b b s 抽样和e m 算法的m o t i f 识别 g - - 幸g i b b s 抽样和e m 算法 例2 3 1 设啦，是来自( p ，r - 1 ) 的样本，并设p ，r 的先验分布分别为 正态分布( o ，1 ) 和r 分布g a ( 2 ，1 ) ，且卢和r 相互独立记y = ，) ， 则y p ，f 的联合密度函数为 烈v ) = ( 2 ”) 一半r 2 唧 一；( 弘一p ) 2 一等) r 唧f r ) “r 的联合后验分布为 枷= 矗鹄拓 故满条件分布为 刊吖，高等 e 计弘1 圳( p 一需) 2 _ ( 需 ( 1 圳- 1 ) ，r ( r l p ，y ) f e 印 一；( 玑一p ) 2 r e z p 一叶 一= l -n = r 1 + 2 唧卜t 1 + i 1 一曲2 】 一i = 1 = g a ( 2 + ；，1 + ；( 玑一| 【) 2 】) 一扛1 2 3 2g i b b s 抽样 和前面一样，假定随机变量x 可以分解为x = ( x i ，) 算法上我们 有两种不同的g i b b s 抽样方案：随机扫描的g i b b s 抽样和系统扫描的g i b b s 抽样具体的方法步骤如下 随机扫描的g i b b s 抽样假定第t 次迭代，我们有一= ( z i ，或) 第t + 1 次迭代按以下两步进行： 1 ：根据给定的概率向量v = ( n - ，n a ) 从 1 ，d ) 中随机地选择一个 i 2 ：从满条件分布”( h 卅) 中生成吒t “令畸；= z f 一日，一+ 1 = ( z 1 ，z t l l l + 1 ) 系统扫描的g i b b s 抽样假定第t 次迭代，我们有一= ( 吐，吐) ，第t + 1 次迭代按以下步骤进行： 基于g i b b e 抽样和e m 算法的m o t i f 识别 第二幸g i b b s 抽样和e m 算法 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ - _ _ _ - _ _ _ _ _ _ _ _ - - _ _ - _ _ _ _ - _ _ _ - - _ _ _ _ _ _ - _ _ _ _ _ - _ _ _ _ _ _ _ _ _ _ _ _ - _ - _ _ _ - _ _ - _ _ - _ _ - _ - _ - _ _ 一_ _ _ _ - 一 1 ：从满条件分布”( 旧- l 】) 中生成z p l i ：从满条件分布”( j 尊1 ，堆j ，吨。，吐) 中生成申1 d ：从满条件分布”( h t 一+ 田1 ) 中生成z 1 令1 = ( 斧1 ，封1 ) 下面我们证明由这两种方法得到的马尔科夫链都使得”为其不变分布 定理2 3 1 对于按随机扫描方法得到的马尔科夫链，”是其不变分布 证明：对于任意的z 和一，假定x l - , 1 = z ；一日，则随机扫描的g i b b s 抽样的 转移核为a ( z ，z ) = 啦”( z 抽一目) 为了证明”是按随机扫描方法得到的马尔 科夫链的不变分布，只需验证引理2 1 1 成立 显然 r ( z ) ( z ，一) = r ( z ) o q ，r ( ：i z l 司) = z ( x i l x l - , 1 ) 霄c x l 一日) 啦7 r ( i 。【叫1 ) = 啦”( ，x l - , ：1 ) - ( x # z h ) = 霄( 一) a ( z ，z ) 故引理2 1 1 成立，从而证明了定理 定理2 3 2 对于按系统扫描方法得到的马尔科夫链，”是其不变分布 证明：对于任意的z 和z ，由系统扫描的g i b b s 抽样的具体步骤可知转 移核为 a ( z ，z ) = ( 。：i z 2 ，j ) 7 r ( 。：i z ：，z 3 ，j 铀) ”( i ，一。，+ 1 ，) ，r ( z ：l 。：，一。) 1 2 基于g i b b s 抽样和e m 鼻法的m o t i f 识别第二章g i b b s 抽样和e m 鼻法 为t 证明绪 仑，只需验证j ”( z ) a ( 毛) 出2 口扛) 丌( z ) a ( 而) 出= 厶正- ，r ( 甸 ( f ，z ) 出t 出一 = 上上”( z - ，) ”( z ：k ，) ”( z ：l z ：，z 。，钆) = 厶厶“阮川“阮蚧， ”( i z ：，一) 二”仁- ，：z d ) d x - 如a = z 。z ，”( z ：胁，黝) ”( 连l i ，幻，) = o 小矗z ：， 黝) 丌( 近j 霉：， 跏) 口( z ：i t ，吐1 ) 锄d x d = w ( i ，z ：) 由上式知结论正确 2 3 3g i b b s 抽样中的分解与分组 在进行贝叶斯推断时，我们经常会把一些讨厌参数积掉这一技术同样 被统计学者应用于m o n t ec a r l o 方法中例如，重要抽样中的维数减少的技 巧就是这样做的在m o n t ec a r l o 方法中，我们把积掉讨厌参数的技术称为 分饵这一节我们介绍一下g i b b s 抽样中的分解与分组 和前面一样，我们假定随机变量x 可以分解为d 个分量，x = ( x ，x d ) 再假设马尔科夫链 x 。= ( 翟，弼) ：t = 0 ，l ，) 是由系统扫描的g i b b s 抽 样产生的由系统扫描的g i b b s 抽样的定义可知马尔科夫链的转移核为 d a ( x t ，a t + 1 ) = i i ”( 霹+ 1 1 群1 ，x t + l ，硌。，弼) i = 1 前面我们已经证明了在这个转移核下”是不变分布 假设x 的后两个分量托“扎可以一起进行抽样这样我们可以按分 组f = ( x l ，m ，一。) 进行g i b b s 抽样，其中墨一，= ( 1 ，) 这样的 一个技巧我们就称为分组进一步，如果施可以被积掉我们就可以关于 善于g i b b s 抽样和e m 算法的m o t i f 识别第二章g i b b s 抽样和e m 鼻法 z ( x 一) = f ，r ( x ) d x d 进行g i b b s 抽样，其中x 一= x a ，托一。) 我们称这个 技巧为分解 下面我们以d = 3 具体说明这两种方法，此时x = ( x 1 ，x 2 ，x 3 ) 分组的 g i b b s 抽样交替地在满条件分布7 r ( 恐，x 3 1 x - ) 和，r ( 五i 局，玛) 中抽样分解的 g i b b s 抽样则不考虑，交替地在条件分布”( x 2 l x ，) 和”( x 。i x 2 ) 中进行抽 样 l i u 在1 9 9 4 年证明了对于系统扫描的g i b b s 抽样，在一定条件下分组和 分解可以改善g i b b s 抽样产生的马尔科夫链的收敛速度，而且分解的g i b b s 抽样的收敛速度是最快的 2 4e m 算法 e m 算法是在有缺失数据的情况下求解参数模型中参数的最大似然估计 的种迭代算法虽然有关e m 算法的文献可以追溯到很早，但是在最一般 情况下的e m 算法的公式化是由d e m p s t e r ，l a i r d 和r u b i n 在1 9 7 7 年解决的 他们还给出了e m 算法具体应用的例子，并在相当一般的情形下建立了e m 算法的收敛性质和一些基本性质e m 算法给不完全数据的分析带来了革 命性的变化，它使得存在缺失数据时我们能有效地进行参数估计e m 算 法还可以用来处理一些困难的统计问题，像混合模型和分层模型等这些 统计问题可能并不涉及缺失数据，但是如果我们恰当地引入缺失数据，这 些问题就很容易处理 2 4 1e m 算法 令y = ( 。 赢) 表示完全数据，其中，r 幽和y 毛分别表示观测数据和 缺失数据再令口表示未知的参数向量完全数据y 的密度函数可以写为 p ( y l o ) = p ( y ，幽i p ( x 曲i k ，口) ， 其中p ( ，毛i 幺，0 ) 表示给定参数口和观测数据，k 后缺失数据施的预测 分布e m 算法就是通过p ( y 赢l y 赢，0 ) 抓住了缺失数据k 。和参数0 之间 的相互依赖关系，当把p ( y 毛l y k ，p ) 看成概率函数时，它就表明了对于现在 的0 有关缺失数据毛的信息；当把p ( i y k ，p ) 看成口的函数时，它就传 达了包含在缺失数据，赢中有关0 的信息有了完全数据y 的密度函数， 我们可以把它的对数似然函数写出来 l ( o i y ) = l ( 0 1 y k ) + l o g p ( y m 。i 。，口) ， 1 4 基于g i b b s 抽样和e m 算法的m o t i f 识别 第二章g i b b s 抽样和e m 算法 其中l ( o i y ) = l o g p ( r 1 8 ) ，2 ( 口l y ，幽) = l o g p ( y , “l o ) ，称l ( o i y o n , ) 为观测数据的对数似 然函数 e m 算法是一种迭代算法它的每一次迭代分为两步：e 步和m 步假 定经过t 步，参数估计值为伊第t + 1 次迭代步骤如下， b 步：对完全数据的对数似然关于缺失数据的预测分布p ( 】，赢l 赢，0 t ) 求期 望，即求 卯= 咿l y ) p ( y 岫i d = e i f ( 日i y ) i ，山，胡 m - 步：找使得e 步得到的q ( o l o ) 最大的参数值，即求舻- ，使 扩1 = a r g n 1 a ：x q c 0 1 e ) 将上述两步进行直到达到某个停止准则，例如两次迭代参数的变化很小的 时候 2 4 2e m 算法的收敛性 e m 算法主要用来求似然函数的最大似然估计和后验分