（概率论与数理统计专业论文）利用贝叶斯统计建立生物信息学模型.pdf

摘要 在生物信息学的研究过程中，往往会伴随着大量的d n a 和蛋白质数 据，对于这些数据如何处理就变得十分必要，为此目的人们发展了很多的 方法。在f 4 1 文中指出，如果假设数据具有某种随机分布，可以通过贝叶 斯( b a y e s ) 统计方法利用获得的数据建立关于d n a 和蛋白质中核苷酸和 氨基酸分布的概率模型，通过一些经验约束条件谨慎选择先验分布并计算 在其条件下所获得数据的似然度，我们就可以求得模型参数的最大后验 概率估计( m a p ) 或者其它估计比如最大似然估计和石验概率均值估计。 最后利用这样的建模方法讨论了两个较为典型的模型：骰子模型和隐马 氏( h m m ) 模型。本文主要的工作就是运用贝叶斯决策方法详细讨论了建 模的主要过程和步骤，仔细分析了两个典型模型，做了大量的计算和证 明，对f 4 】文作了重要补充和改进。 c 贝叶斯统计，数据似然度，先验分布后验分布，最大后验 概搴估计，骰子模型，跨马氏模型 a b s t r a c t i nt h er e s e a r c hp r o c e s so fb i o i n f o r m a t i e s t h e r ea r el o t so fd n ad a t a a n dp r o t e i nd a t a ，s oh o wt od e a lw i t ht h e s ed a t ai sv e r yn e c e s s a r ya n dw e h a v ed e v e l o p e dm a n ym e a n sf o rt h i sp u r p o s e a r t i c l e 4 t e l lu 8 i fs u p p o s e d a t ai ss t o c h a s t i cs e q u e n c e ，u s i n g b a y e s - s t a t i s t i c sw e c a ne s t a b l i s hp r o b a b i b i t ym o d e l f o rr i b i t o l a m i n oa c i dd a t ao fd n a p r o t e i n u s i n gs o m ee x p e r i e n - t i a lo b l i g a t i o nw es e l e c tp r i o rd i s t r i b u t i o nc a r e f u l l ya n dc a l c u l a t el i k e l i h o o d o ft h e s ed a t ai nt h i sc o n d i t i o n ，s ow e g a i nm a p e s t i m a t i o no fm o d e l sp a r a m - e t e ra n do t h e re s t i m a t es u c ha sm l ( m a x i m u ml i k e l i h o o de s t i m a t i o n ) a n d m p ( m e a np o s t e r i o re s t i m a t i o n ) f i n a l l y , b yt h i sm e a n sw ed i s c u s st w or e p r e s e n t a t i v em o d e l ：d e v i l sb o n e sm o d e la n dh i d d e nm a r k o vm o d e l ( h m m ) i n t h i sp a p e r ，id i s c u s s e dc o u r s ea n ds t e po fm o d e l l i n gp a r t i c u l a r l yb yu s i n g b a y e s i a n - d e c i s i o nt h e o r y , a n a l y z e dt w or e p r e s e n t a t i v em o d e l sc a r e f u l l y , m a d e o u tl o t so f c a l c u l a t i n ga n dp r o o f s ，p u tu pi m p o r t a n tr e c r u i ta n dm e n d e d f o r a r t i c l e 4 。 k e y w o r d s ：b a y c ss t a t i s t i c s ，l i k e l i h o o d ，p r i o rd i s t r i b u t i o n ，p o s t e r i o rd i s t r i b u t i o n ，m a pe s t i m a t i o n ，d e v i l sb o n e sm o d e l ，h i d d e nm a r k o vm o d e l ( h m m ) 2 1 前言 生物信息学( b i o i n f o r m a t i c s ) 是一门相当年轻的学科，它是伴随着八 九十年代计算机技术的迅猛发展，同时得以获得自身的发展的。无论从理 论上来讲还是从现实情况来看，生物信息学的实质就是利用计算机科学和 网络技术来解决生物学问题。二十世纪尤其是末期，生物科学技术的迅猛 发展，无论从数量上还是从质量上，都极大地丰富了生物科学的数据资 源，数据资源的急剧膨胀首先迫使我们不得不考虑寻求一种强有力的工具 去组织他们，以利于对已知生物学知识的储存和进一步加工利用。大量多 样化的生物学数据资源中必然蕴含着大量重要的生物学规律，这些规律是 我们解决许多生命问题的关键所在，然而继续沿用传统手段以人脑来分析 如此庞杂的数据实在是太困难了1 人们同样需要寻求一种强有力的工具去 协助人脑完成这些分析工作。可以说，伴随着二十一世纪的到来，生物科 学的重点和潜在的突破点已经由二十世纪的试验分析和数据积累转移到数 据分析及其指导下的试验验证上来，生物科学也正在经历着一个从分析还 原思维到系统整合思维的转变。因此，我们所寻求的那种强有力的数据处 理分析工具就成为未来生物科学的关键所在：与此同时，伴随着生物科学 这一需求的加剧，以数据处理分析为本质的计算机科学技术和网络技术同 样获得了突飞猛进的进展，自然就成为生物科学家的必然选择，计算机科学 技术和网络技术日益渗透到生物科学的方方面面，- - 1 7 崭新的、拥有巨大 发展潜力的生物信息学也就发展和成熟起来了! 而生物信息学一生物科学 与计算科学的融合体一就成为下一代生物科学研究的重要工具。 随着生物信息学发展到今天，我们已经掌握了大量的关于生物微观 结构层面的信息，特别是d n a 序列和蛋白质建立了相应的数据库， 如g e n e b a n k 和p d b ，但关于这一领域还缺乏一套在分子水平上理解生命 组织的系统理论。于是退而求其次，利用庞大的数据库通过推理、模型匹 配或者样本学习在数据中自动学习理论的机器学习方法应运而生，机器学 习方法正是适合这类数据量大、含有噪声模式并且缺乏统一理论的领域， 关于机器学习方法的统一理论体系就是用于建模和推断的贝叶斯概率体 系。在接下来的叙述中，本文先是讨论了原始数据的处理问题，简单说明 在通过一些譬如对数据组赋权重和减少数据冗余度的工作后，得到新的数 据组，这种新的数据组是更加“自然的”和更有代表性的。接着运用统计 学中的贝叶斯决策的理论，首先讨论了建立生物模型的理论框架，分析了 在建立模型的过程中起主要作用的先验概率的选择问题和数据似然度的计 算问题，并分析了比较典型的先验概率如高斯( g a u s s ) 先验概率和狄利克 3 莱( d i r i c h e l e t ) 先验概率，证明高斯先验概率满足最大熵原则，而狄利克莱 先验概率是多项分布的共轭分布，正是这些特点使得他们更加有用。在做 了详细的讨论和分析之后，理论框架已经完备，接着分析了两个典型的模 型；骰予模型和隐马氏模型，骰子模型主要是假设序列中字符出现是随机 的和独立的，利用这样的假设在选择狄利克莱先验分布的情况下，运用序 列数据通过贝叶斯推断得出了模型参数的估计，包括最大后验概率估计和 后验概率均值估计，并分析了与最大似然估计的不同。隐马氏模型则假设 序列的生成和它上面每一个字符的出现都是随机的和独立的，并且这两个 随机序列本身也是独立的。我们同样假设模型参数的先验分布是狄利克莱 形式的，利用这样的假设建立的隐马氏模型可以模拟数据序列的生成，不 同的模型参数代表不同的具体模型，利用不同的具体模型生成特定序列的 概率的不同，在对已知数据组进行比较后，可以比较出最好的具体模型， 得到参数的最优估计。 2 数据的处理 与生物功能和进化相关的链状分子具有一种基本特性，即他们能够以 数字化符号序列的形式表示，d n a 、r n a 和蛋白质中的核苷酸和氨基酸单 体是确定的，这一点使得它明显不同于其他的科学实验数据，虽然如此， 在大规模d n a 测序的基因组项目中，或者进行蛋白质直接测序时，研究目 的、信息检索能力、承担项目的机构、伦理和经济等因素都会影响数据的 质量。另外还必须考虑生物序列的噪声，这种噪声表现为序列片断的多样 性，这种多样性部分地被来自进化和突变放大了的随机事件，特定序列在 不同物种中可能具有不同的形式，甚至在具有多态性的情况下在相同物种 中也具有不同的形式，在对一组序列进行综合研究时，由于具有特定功能 和结构的d n a 和蛋白质序列存在一些不确定的差异，这样序列模型也必 然是基于概率理论的模型。特定序列的数据的躁声来源于两个方面，一个 是对实验结果的错误解释和实验手段本身，另一个是公共数据库中的数据 处理存储不当产生的，并且由于生物序列数据都是用电子化的方式来存储 的，维护数据库的人员组成也极为不同，数据库中的数据更是由不同的生 物学家和生物信息学家提供的，因此后续信息的处理引入的误差更是远大 于初始的实验误差，正是由于这点，使缛采用概率理论建立模型成为必要 的。 另岁 一个困扰研究者的问题是数据冗余。蛋自质或者基因数据库中的 许多纪录代表一些蛋白质和基因家族的不同成员，或者是在不同生物体中 4 发现的同源基因的不同版本。一些研究小组可能分别提供了相同的序列数 据，因此相应的项即使不是完全相同的也是密切相关的。在测序项目 中，典型的数据冗余来自不同的实验手段本身，由于存入数据库中的数据 是以广泛多样的手段和渠道获得的，数据库中的不同记录可能表示同一基 因，虽然这些记录之间会存在不同程度的差异。在蛋白质数据库中，特定 基因也可能表示为某种氨基酸序列，但该序列并不对应于原始核苷酸序列 的直接翻译，比如蛋白质序列经常需要通过微小的修饰，以便获得更好的 结晶，利于采用x 射线晶体衍射法测定结构。氨基酸的删除和替换也是 产生数据冗余的常见原因。使用具有冗余的数据集至少会导致三种误差： 第一，如果数据集中的氨基酸或核苷酸序列包含很大的密切相关的序列家 族，统计分析将偏向这些家族并侧重描绘他们的特性；第二，序列不同位 置之间表面上的相关性可能是对序列数据进行有偏倚的采样所导致的认为 特征；第三，在我们使用数据集对某一特征进行预测或用于选择标定预测 方法时，如果用于训练和标定预测方法的训练集的数据与用于测试的序列 相关性过于密切，显然会过高估计预测方法的性能。因此，首先对数据集 进行整理，使低显现度( u n d e r r e p r e s e n t a t i o n ) 序列得到均等的机会，会获 得较好的效果。在大多数情况下，可以采用一个数据的相似度阈值，或者 根据、传统的蛋白质和基因家族列表，在每一个家族中选择一个成员从而 构造出一个代表数据集；另一种替代策略是保留数据集中的全部序列，根 据序列的奇异度赋予他们不同的权重，对于密切相关序列的预测将得到较 低的权值，丽相关距离较远的序列则赋予较大的权值，这一方法的风险在 于即使错误数据也赋予了较大的权值，从而放大了错误。序列谱方法可以 有效的利用数据冗余，它通过描述多重序列比对组织起来的每个位置上的 氨基酸的变化，可以设计为机器学习方法中输入数据的表示方法。经过处 理后的数据库数据，可以认为是从自然界中随机抽取的，这就为利用统计 方法建立模型提供了基础。 3贝叶斯统计决策 设经过处理后的数据集记为d ，我们试图通过它来推导出一个参数 化的模型m ) 是所有参数组成的向量，一个复杂模型可以很容易 的简化为一个二值命题：使用模型m 似) 解释数据d ，误差为。理 论上讲，设所有可供选择的参数模型m ( u ) 构成一个行动空间a ，数 据d 的分布为f 如，m ) ，五( 尬d ( o ) ) 表示模型是m 时选用j ( z ) a 引起的损失函数，设所有模型组成一个广义空间廊，日庙是府的一 5 切b o r e l 子集构成的口一域，若我们在b 时上引进测度u ，定义b a y e s 风 险r ，( 6 ( z ) ) = e ( l ( m ，6 ( z ) ) ) ，使扁，( 6 ( z ) ) 达到最小的5 0 ( x ) 就是我们要求 的贝叶斯决策。但是实际上，我们无法给出b 时的具体构造，也就谈不上 它上面的测度u ，我们也难以给出损失函数l ( m ，6 ( z ) ) ，因此这种讨论只 具有理论上的意义，在实际工作中，我们一般限定一类模型，而把问题简 化为估计参数，我们将利用贝时斯方法来估计它们的参数值。 我们假设选用的模型为m ( ”) ，则由贝叶斯定理 p ( 1 1 d ) _ 署- p ( m ) 鬻( 1 ) 先验概率p ( m ) 表示在没有得到任何数据的前提下所估计的模型m ( w ) 的 概率，后验概率p ( m i d ) 表示在得到观测数据d 之后重新选用模型m ( ”) 为真的概率，p ( d m ) 是数据似然度，p ( d ) 是常数参数。出于技术上的原 因，这些数值般都非常小，因此我们用相应的对数来处理 l o gp ( m l d ) = l o gp ( d i m ) 4 - l o gp ( m ) 一l o gp ( d ) 在将这一公式应用于任何种类的模型前我们要计算先验概率p ( m 1 和数 据似然度p ( d i m ) ，一旦这两项清楚了，初始的建模工作也就完成了，剩 下的只是进行计算。对于依次获得的数据d 1 ，d 2 ，d 。，由贝叶斯定理有 p ( m i d l - d ) = p ( m i d l ，d 一1 ) ! 箬盖兽叁暑嵩( 3 ) 这样，修正前的后验概率p ( m i d l ，d t 一1 ) 可以作为新的先验概率，利 用新获得的数据重新给出修正后的后验概率，从而可以计算引入新数据 后的更合适的参数模型。 3 1先验概率的讨论 使用先验概率是贝叶斯方法的一个优点，因为它允许将先验知识和约 束条件导入建模过程，但是先验概率带有主观性，这也是它易受攻击的地 方。针对这些反对意见，贝叶斯学派提供了几种解答： o ) 通常说来，随着数据量的增加，先验概率的作用逐渐减少。形式上 看，负对数似然度一l o g p ( d i m ) 随着d 中数据量的增加呈线性增加，而 先验概率一1 0 9 p ( m ) 保持不变。 b 在通常有些情况下，可以使用一些客观准则，如最大熵原则、群不 变性等来确定无信息的先验概率值。 6 c 在没有明显提到先验概率时，它们也被隐含的使用了，贝叶斯方法 在解决问题时不一定需要揭示出隐含的先验概率问题，这也是贝叶斯方法 的妙处之。 d 最后，也是最重要的，与不同模型及模型类别一样，不同的先验概 率的影响可以在贝叶斯体系内通过比较相应的概率来评估。 由于选择某种先验分布以及相应的数值结果直隐含在整个概率计算 过程中，圜此我们最好是抱着尝试的态度来选择先验分布，但在某些场 合，可以依靠最大熵原则和群论的知识来选择先验分布。最大熵原则规定 先验概率应该是与所有先验知识或约束条件相一致的概率分布中熵值最大 的那个概率。因此可以说先验分布是“最少假设的”，“最大模糊的”或 具有“最大不确定性的”。根据拉普拉斯( l a p l a c e ) 无差别原则，缺乏先验 约束将导致均匀分布。因此，若是除了参数”的范围以外没有其他关予 t o 的可用信息对，在此范围上的先验概率为均匀分布是自然的选择。在实际 建模中最大熵可以由分布f 或者相应的直方图来确定：若是先验分布的某 些约束条件可以用群论的语言来表达，则我们可以利用群不变性来限制先 验概率的种类甚至决定先验概率。在先验分布不是均匀分布时，经常应用 的是高斯先验分布和伽玛先验分布，单参数高鲰先验分布为： 脚0 2 ，= 去e 卵 _ 譬) 性质：当关于钾的信息只有均值p 和方差一2 时，w 的满足最大熵原则的 先验分布就是高斯先验分布( 协，口2 ) 证明：事实上，设”的密度函数为p ) 芦的取值范围为谚，则它的 微分熵为口( ”) = 一岛p ( t l ，) l o g p ( ) d w 利用约束条件岛p ( w ) w d w = p 和岛p ( ) 似一p ) 2 d w = 口2 以及归一化条件岛p ( ) d w = 1 构造拉格朗 日乘子 l = 一p ( ”) l o g p ( ”) d w 一 1l p ( t ，) ( 1 1 3 - - 1 2 , ) 2 d w wi j w一。2 jjj 吨阮p c 咖d w p 卜s 跏”) d w 一， 取工对p ( w ) 的偏导数，则有 瓣a l2 厶卜l o g p ( ) - 1 - ) q ( w - p ) 2 咄 d w 我们令a = 0 ，并考虑到谚的任意性，则有 一l o g p ( w ) 一1 一, k 1 ( w p ) 2 一a 2 ”一 3 = 0( 4 】 从而由上式可得 p ( ) 2 三币i 了1 丽e x p 一 l ( ”一p ) 2 1 2 w ) ( 5 ) 式中z ( a 1 ，a 2 ，a 3 ) 是归一化因子，由( 5 ) 式的表达式可知满足约束条件且 具有最大熵的先验分布是高斯先验分布，从而= ( 一0 0 ，+ 。) ，我们可具 体求得z ( a l ，a 2 ，扎) 及a 1 和a 2 。由j 茗p ( w ) d w = l 得 r + 1 _ 。丽轰两e 印，( w - - l z ) 2 _ a 2 w ) d w 1 考虑到z ( a l ， 2 ， 3 ) 是常数，我们积分该式得 堡! 墨二型：， 1 z ( a 1 ，入2 ， 3 ) ( 6 ) 由j ! ：p ( w ) w d w = p 把( 5 ) 式代入计算得 0 一杀) 瓮摅导：p 比较以上两式可得沁= 0 从而有 p ( ”) 2 郅拳丽唧 呐( ”一p ) 2 ) ( 7 ) 把( 6 ) 式代入j 茗p ( ) 一p ) 2 d w = 口2 既有 丽赫甜 ( 8 ) 2 z ( a 1 ，a 2 ，a 3 ) ；7 。 由( 6 ) 式( 8 ) 式联立得a l = j 1 ，z q l ， 2 ，a 3 ) = 、；所以 p c 妒志e 印 _ 譬) 这样我们就完成了该性质的证明 我们还可求得它的微分熵 丑( w ) = 一7p ( ”) l o g p ( w ) d = ；l o g ( 2 ，r e o 2 ) 对于n 维空问的随机向量x ，若是所知的信息只有均值向量p 和方差矩 阵c ，则符合最大熵原则的先验分布是r l , 维高斯先验分布n ( x i p ，g ) 具 有参数n 和a 的伽玛分布密度函数，其形式为 忡圳= 焉扩 以：；拦 8 其中r ( n ) = j ”z ”1 e d x ，通过调整a 和a 并改变w ，伽玛密度函数 可以是一个很大范围的先验分布，特别是当参数w 是单边有界时，伽玛先 验分布有着非常广泛的应用。 另外一个常用的先验分布是狄剥克莱( d i r c h l e t ) 分布，具有参数。和 向量q = ( q 1 ，q 2 ，口n ) 的随机向量p = ( v l ，p 2 ，h ) 的狄利克莱分布 的密度函数为 。 巩。( p ) = 鼎酚一 其中c t ，a ，吼0 俄= 1 ，肼= l ，正是由于p i = 1 ，我们可以改写 上式为 三k 口c p ，= f z 五；j i p ? 吼一1 ( - 一蓦n ) 。4 “一1 定理：设p = ( p l ，h ) 一d 。o ( p ) 则它的混合原点矩为 e ( p ：1 r - - - 1 ) = 絮券器等等雩掣o o ， 证明：令s = p ，= ( p ”一，p n 一1 ) ，p l ，p n - 1 0 p i 1 ，则 e o ：1t p n - i ) = z f i 苫i 蹁p ? q l + r l - - 1 一。a q a 。- l + r n 一，一l n 一1 o 如一1 l 一a d p r 蛾一z 一 厂r ( a )f ( a q l + r 1 ) r ( a q 。一1 + r n - 1 ) 1 1 ( a 口n ) j s r ( a q l ) r ( a q )r ( a + r 1 + - + 一1 ) 一 一 r ( 口q l 端r l 甚r ( 杀a q 1 等鲁。而r ( a q 黔1 1 叶r i - l + ) 一+ r ，)，“ n 一1 、“口n l 1 一弘) d p l d p n 一1 ：! ：f 婴! 韭：! ! 墼= ! 垒= 四堕2 1 幽 r ( a q z ) r ( a q n ) r ( q + r l + ，+ 一1 ) j 厂s f ( a q 。蔫高等， n - - 1 亘1 1 p p 。 + r 1 ) r ( o q 。一l + ) r ( ) ：， n 一1 o 如一1 l e p , 劫1 蕊一l 9 ：! ! ! ! ! ! 12 ：! 生! 堑= ! 垒= 12 1 ( 竺堑2 1 1 1 2 f ( o t q l ) r ( a q 。) r ( a + r l + + r 。一1 ) z d n ，“p ) d p - d p n 一- = 墅r ( 。等q l 崭r ( a q 堕。) r ( 未a 等等掣 ) ) -+ r l + - - + r n 一1 ) 卜17 我1 门令0 ，= o + r l + + r n l ，o = ( 学，塑= 芦，笔 ) ，则有 。o ，掣c p ，= f 再7 磊蘑茁1 ( ，一喜肼) “吐1 所以由归一化条件有b d 掣( p ) 咖l 。如h 一1 = 1 ，在上面的证明中，我们 就是利用这个结果得到了该定理。利用该定理，我们可以求出 e 眩) = 舔， l ，n w 慨) = b 簖) 一e 2 慨) = 丽( a q i + 1 ) o t q i 一毋= 警斧 以及当 j 时， c 鲫眩，巧) = e 慨乃) 一吼劬= 涨一吼劬= 一鬻 由此我们可以看出q 为该分布的均值。则决定了该分布曲线在均值附近的 光滑程度。狄利克莱之所以重要应用广泛，就是因为它是多项分布的共轭 分布，共轭分布就是说参数是狄利克莱分布，若样本是多项分布则可以引 入更多的样本来估计参数，参数的后验密度函数还是狄利克莱分布。这一 点缀容易说明，事实上若样本x = 扛l ，$ 。) 是参数为p = ( p l ，肌) 的多项分布，我们把它写为对参数p 的条件分布 肛= 志垂砖= 志黔( - 一斟k ， 其中l = 1 1 + + 2 。为样本总量，出现第i 种情形的样本量为z i ，若参数p 的先验分布为狄利克莱分布d 。o ( p ) ，其中p l + + 肌= 1 ，则我们可以 计算后验密度函数 h ( p x ) = c - d a q ( p ) i ( x i p ) = 耥静( t 一斟= lk 巍p p 俨- g - t ! r ( a 1 2 可_ 礤碌石f 确 l o 孽p ? + n “一1 ( t 一薹鼽) k + 。和一1 其中c 是归一化常数，由上式我们可以明显看出h ( p i x ) 是狄利克莱分布， 令s = ( p l ，p 。) p l ，p 。o ，p l + + p 。= 1 ) 我们利用归一化条 件如h ( p i x ) d p = 1 ，求得 。r ( a + 1 ) l l ! - f n ! f ( a q l ) r ( a q n ) u2 可面呵丽干百_ 币丽 从而h ( p x 1 的准确结果为 n c p i x ，= i 1 ：i i i j j ；i ；j 嚣p ? “+ “一1 ( - 一喜肌) 。4 “+ ”一1 ( 1 3 ) 至此我们完成了狄利克莱分布是多项分布共轭分布的证明。 3 2 数据似然度 要想定义数据似然度p ( d i m ) ，必须掌握如何从一个模型m 产生不同 的观测集d ，因为在用贝叶斯理论建立模型时，所生成的序列模型是概率 形式的，如果没有给定似然度，序列模型的科学描述它们如何拟合数据 以及它们之间怎样进行比较是不可能的。似然度问题与变异以及噪声问 题有关，生物序列本身就存在噪声，而进化使随机事件的效应放大并最终 导致变异。特定个体序列与一个家族中的“平均”序列之间不匹配和存在 差异是必然的，这一点需要进行量化。因为同样的d n a 序列或氨基酸序列 即使在同一物种的不同个体之间也是有差异的，在不同的物种之间差异更 大，所以我们用概率的观点来考察模型，那么从这一模型产生不同的数据 集d ，也是自然的了。似然度的计算当然是与模型相关的，因此我们无法就 其进行一般性的讨论，然而存在一些准则来推导模型，在利用这些准则之 后。似然度的计算并不困难。但是不论是用什么样的准则来衡量模型和数 据之间的差异和误差，这些准则都必须是来自一个明确给出的基本的概率 模型，并使之符合贝叶斯分析的检验。 如果具有参数向量 的模型m = m ( w ) 可以由某个被称为误差函数 的f ( w ，d ) 0 来评价，其值越大则说明其效果越差，那么相关的似然度可 以定义为 p ( d i m ( w ) ) = 百一 (14f(w) 。一，d ) z j 其中z = c 。一( ”，d ) d w 是保证似然度积分为1 的归一化常数因子。这 样，最小化误差函数与最大似然度等价，符合最大似然估计的参数估 计就是所求的估计。虽然对于一般的数据似然度没有确定的计算方法， 不过对于确定的模型则很容易给出，例如若某随机向量x 的密度函数 是f ( x ，日) ，口是模型的参数，则简单随机样本( x l ，恐，弱) 的似然度 为p ( ( x l ，x 。) l f ( x ，口) ) = n ，( 置，口) ，因此对于具体的含参数模型， = 1 n 我们可以直接给出样本的似然度，令l ( o ) = 兀，( x i ，口) ，我们可以利 用变分法来求得使l ( o ) 达到最大的参数日为了比较两个模型类m 1 ( w ) 和i 2 ( w ) ，设采得的数据集为d ，我们可以定义某一类模型的似然度 ，厂 p ( d i m ( w ) ) = p ( d ， i m ) d w = p ( d w ，m ) p ( w i m ) d w ( 1 5 ) ，”cj 伽g 其中d 为参数”的定义域。( 1 5 ) 式可以看作确定参数模型情形下数据的 似然度p ( d 1 w ，m ) 的数学期望；当似然度p ( 口i u ，m ) 在其最大值附近形 成形状较为尖锐的峰时，期望值可以用最大概率值近似，对于一般情形 只有借助计算量很大的蒙特卡罗( m o n t e - c a r l o ) 方法来模拟计算期望值， 若p ( d i m l ) p ( d i m 2 ) 则我们称尬优于尬。 矗就是要选用的，我们 还可进一步进行参数估计来确定参数。 在讨论了先验概率p ( m ) 和数据似然度p ( d i m ) 之后，我们就可以通 过公式( 2 ) 来比较任意的两个模型，从而选用最优的模型。我们把公式( 2 ) 在写一遍 l o gp ( m i d ) = l o gp ( d i m ) + l o g p ( m ) 一l o gp ( d ) 我们把问题分为两部分：一是首先选择模型类，二是对于选定的模型类确 定模型中的参数。在讨论先验概率的时候我们已经讨论过模型类的选择， 除了那里说的根据一些经验原则对模型的限制以外还有一个经常被采用的 原则一奥卡姆( o k a m ) 剃刀原则，即在各种因素都相同的情况下，人们应 尽量选用最简单的模型。这一原则有时候自动体现在计算过程中，首先人 们很容易选择对复杂模型进行惩罚的先验概率，即复杂模型的先验概率值 较小；即使不考虑先验概率，参数化的复杂模型也必然倾向于数据量大的 样本情况。其次是由于数据似然度p ( d i m ) 在数据空间上的和为1 ，即很 显然的f p ( d i m ) d d = 1 。如果p ( d i m ) 覆盖了数据空间的体积较大，那 么这个数据集的似然度的平均值自然较小，由于复杂模型产生的数据空间 体积更大，因此对于固定的观测数据，复杂模型得到较小的似然度。由于 正值更易于处理，我们把公式( 2 ) 改写为 一l o gp ( m i d ) = 一l o g v ( d i m ) 一l o g p ( m ) + l o g p ( d ) 1 2 从优化的观点看，先验概率的对数起到了正则因子的作用，这相当于一个 附加的惩罚项，它可以用来体现附加的约束条件如平滑性等。这在先验概 率的讨论中有过说明。p ( d ) 起到归。化常数的作用，它不依赖于w ，因此 与优化问题无关，这样问题就简化为优化一l o g p ( d i m ) 一1 0 9 p ( m ) ，我们 在根据经验选择出几类模型后，首先最小化一l o g p ( dj m ) 一l o g p ( m ) 得 到模型类，这些需要优化的函数因为选择模型的不同可能变得十分复杂， 因此发展了一系列相应的算法来专门针对这一问题；在得到模型类后，确 定参数就变得相应容易了。接下来我们简单讨论两个模型 4 序列数据的骰子模型 数据集d 由d n a 链构成，每条链由字符集a = ( q c ，g ，) 中的字符 组合而成，这样每条d n a 链就是一组字符串，我们假设字符串是通过独 立投掷同一个四面体骰子得到的，这样我们可以把多个字符串看成一条较 长的字符串，由于投掷的独立性，它并不影响我们的分析。我们假设这个 观察序列即投掷字符串长度为n ，则可以记为d = x 1 x 2 x “) ，x 。 a 。我们考虑一个四个参数的骰子模型肘，其参数为r ，p c ，p 口，b 且满 足r + r + p 口+ 最= 1 ，这里它们分别表示对应各个字符出现的概率，则 显然在模型m 下出现数据集d 的似然度为 p ( d i m ) = p ( x 1 x 2 x i m ) = p ( x 1 i m ) p ( x 2 i m ) p ( x i m ) = p ：。p ! c 式。式t 这里x 表示字符x 在序列d 中出现的次数。我们把上式代入 一l o gp ( m d ) = 一l o g p ( d l m ) 一l o g p ( m ) + l o gp ( m ) 因此负对数后验概率为- l o g p ( mj d ) = 一n x l o g 取一l o g p ( m ) + x e a l o g p ( d ) 事实上由于我们是对固定的模型估计参数，我们可以改写上式为 一l o g p ( w l d ，m ) = 一文l o g p x l o g p ( w j m ) + l o g p ( d )( 1 6 ) x e a 我们先假设这四个参数具有均匀分布，则最大后验概率参数估计就是最大 似然估计，也就是最小化 一l o g p ( w i d ，m ) = 一n x l o g p x 0 7 ) x e a 1 3 我们构造拉个朗日( l a g r a n g e ) 乘子 l = 一n x l o g p x a ( 1 一既) ( 1 8 ) 令o l o p x = 0 ，得取= n x ) , ，由p x = l 则得到a = n ，最后我 们得到尸y 的最后的最大似然估计 砖= 等 ( 1 9 ) 当很大时，用观测频率n x ，估计段是很自然的。大数定律也 保证当值足够大时，观测频率将与真实的段很接近。但当n 很小 时，比如在一个长度为2 0 的序列中只观测到一次字符o ，我们用频率t 2 0 来估计只就是不准确的，因为我们没有理由认为骰予是高度偏倚的，也就 是说我们的先验知识并不认为模型的参数r 这么小。在这种情况下合理的 参数先验分布并不是均匀分布，而应该是关于参数向量p = ( r ，p c ，b ，p t ) 的狄利克莱分布d a q ( p ) ，其中n 和q = ( q 。，q e ，q g ，q ) 是待定的参数。我 们写出它的表达式如下： 玩舻卜蚤尚飘璎”1 x o 此时数据组的似然度仍然是p ( w i d ，m ) = n 世。先验概率的对数为 x 6 a l o g p ( w l m ) = l o g d 。口( p ) = ( a q x 1 ) l o g p x l o g z x e a 其中上式中k 唔z = 一l o g r ( a ) + l o g r ( a q x ) 是狄利克莱分布的归一化 x 6 a 常数与参数p x 无关。我们把它代入( 1 6 ) 式，化简整理后得到负对数后验 概率为 一l o g p ( w l d ，m ) = 一( n x + a 驻一1 ) l o g p x + l o g z + l 。g p ( d ) ( 2 0 ) x 6 a 由( 2 0 ) 式我们与( 1 7 ) 式比较看出除了由n x + a q x 一1 来代替n x 外，最 大后验概率估计闯题与( 1 7 ) 相同，通过构造拉格朗日乘子经过相同的计算 我们得到对于任意的x a 最= 面n x + i o t q x 可- 1 = 等等 ( 2 1 ) 舷5 可了而2 订i = r l 趴j 1 4 这样我们看到，采用狄利克莱先验分布的作用相当于在观测到的次数上 增加一个虚计数项，通过适当的选择参数q 和a 我们可以更好的来估 计取特别对于字符出现很不均衡的序列更是如此。通过( 2 0 ) 式我们可以 得到 p ( w l d ，m ) = 寿面e 璎对鲫”1 如果我们令卢= o + n ，r = ( r 。，r 。，勺，r t ) 其中v x a ，7 。x2 ( x + a ) ，( o + ) 。可以验证后验概率分布p ( w i d ，m ) 就是d p r ( p ) ，由我们 在上节证明的定理得e ( p x ) = r x ，这不同于( 2 2 ) 式中的估计，我们可以 使用r x 作为p x 的另一种估计称为后验分布均值估计 p ：n x + a q x( 2 2 ) 单一狄利克莱先验分布的骰子模型是一个十分简单的模型，我们可以利用 解析方法进行更高层次的贝叶斯推断，例如计算p ( d ) ，事实上 p ( d )= p ( d i 研，a 彳) p ( 叩i m ) d w ：型d，n(p)dpr(f，)l-ir ( a q x ) x f 2 庐z ：。x ：尸x = l r ( 口) nr ( 砂x ) r ( z ) y if ( a q x ) r ( 卢) e ar ( f j r x ) d p x ( 2 3 ) e a 根明显，这实际上是先验分布与后验分布归一化常数的比值 如果我们假设先验分布是狄利克莱分布的混合，即r ( p ) = ) , i d a 。硪( p ) 时，v x a 我们利用相同的方法同样可以求得最大后验概率估计尸未和后 验分布均值估计p # 。结果如下： n x + a i q x 一1 厮5 可了曩而 1 5 擎 严 p 鬻弦揣删也州蹦zr船鳃揣 n x + o t i q x 蹿。可索 而后验分布为群( p ) = k d 风几( p ) 其中屈= n + o t i ，n x = ( n x 十 啦q x ) ( n + a ，) 。事实上我们可以选用多项分布的任何共轭分布的混合来作 为先验分布，由此得出的后验分布就是相应的共轭分布的混合，只不过参 数改变了而已。 对于这个彀子模型，我们还可以考虑数据集由字符出现的次数d = 蜥) 组成的情形。在这种情况下，数据的似然度为 p ( 。嗍2 戛n 舷! r 始1 p n n x ( 2 4 ) x “ 其中n = n x 。这是显然的，前面的因子项表示由字符集实现的长 度为的字符串中字符出现次数为d = 蜥) 的字符串的所有可能的 方式，这就是排列数，后面就是出现一个确定字符串的似然度。如果参 数向量尸具有狄利克莱分布d a q ( p ) ，通过与确定字符串相似的计算我 们得到参数向量的后验概率分布也是狄利克莱分布记为d 口冗( p ) ，其中参 数卢= o t + n ，r 的分量啊= ( n x + a q x ) 陋+ ) ，v x a 同样我们 也可以求得相应的最大后验概率估计尸；和后验分布均值估计尸。但实 际上。先验狄利克莱分布不是满足虽大熵原则的先验分布，我j f 乖j 用斯特 令( s t i r l i n g ) 公式竹! v 西茅磊( ：) “对( 2 4 ) 式取对敷，得 l o g p ( d 1 w ，m ) = l o g n ! 一l o g n x ! + n x l o g p x l o g n ! 一+ l o g 、孬丽+ 坛l o g p x 一( n x l 。g n x n x + l o g 厕) = l o g n ! 一n x l o g 警+ g ( 2 5 ) 其中g 是仅依赖于n x 的常量，由于n x n = 1 ，则 n x n 是经 验概率分布从而定义它与既的相对熵如下 n引二删n。，xh(nx - d - l o g 麓? t r i o ，段) = 盐 这样，( 2 5 ) 式就可以改写为 l o gp ( d 1 w ，m ) = g 7 一h ( n x n , 既)( 2 6 ) 1 6 式中c ，是另一一个只与n x 有关的常量。在参数p 的分布是均匀分布的情 况下很明显u c x l n ，尸x ) 与经验分布的熵h ( n x n ) 只差一个负号和一 个附加的常数，这里经验分布 n x n ) 的熵定义为 日( 蜥) = 一三百n x l 锯万n x 因此这时候似然度可以表示为 。h ( n x n ) - p ( d 1 w ，m ) = 厂一 ( 2 7 ) 刖 其中z 是归一化常数，这被称作熵分布。也就是说，由均匀分布p 可以导 出一个关于次数蜥的、在所有可能直方图构成的空间之上的熵分布，这 正是最大熵原则的主要初衷之一，相当于采用先验熵分布。如果采用先验 熵分布 。 p ( i m ) e - h ( 尸) e x p ( 一p l o g p ) t 2 万一 我们计算得到的最大后验概率估计仍然是p ；= n x y ，后验概率分布则 十分复杂，我们不再作讨论。 尽管简单的骰子模型十分粗糙，但它正是我们计算l 阶统计量时所使 用的模型例如一个给定的序列集合( 如外显子、内含子或蛋白质家族) 中每个字符出现的比例。简单彀子模型还可以通过将其每一面的字符换成 固定的字符串而得到更为一般的推广，这等价于扩展字符集，比如可以使 用一个6 4 面的骰子构造d n a 密码子模型。我们还可以把简单骰子模型 稍加推广得到复合骰子模型，在这样的模型中，我们假设数据组由耳个 序列组成，每个序列的长度为，例如我们可以考虑一个耳个序列的多 重序列比对，其中间隙符号“一”可以看作字符集中的一个字符，在这样 的复合骰子模型中我们假定有个独立的骰予，每个骰子对应与序列的 一个位置，这样每个序列都是这个骰子按照一定的顺序投掷得到的结 果。令p 表示第i 个骰子掷出字符x 的概率， 表示在第i 个位置 上出现字符x 的次数，则显然有 = k 。由于假定骰子和序列都是独 n 7 i 立的，则数据集的似然度为p ( d 1 w ，m ) = h 兀p 釜x ，若所有的骰子具 = 1 正亘一 有相同的均匀分布，则最大后验概率估计就是最大似然估计，对v x a 有p 荨= 譬g 若所有骰子均具有相同的狄利克莱分布d 。o ( p ) ，则最大 后验概率估计为毁= ( 峨+ a 似一1 ) ( k + o t i a i ) 后验概率均值估计 为p p = ( + a q ) ( k + a ) ，同样我们也可以假设各个骰子具有不同 的先验狄利克莱分布的情形或者是不同的混合情形，其计算方法都是相似 1 7 的，骰子模型虽然简单，但它可以看作是一个迭代建模过程的第一步，因 此有其意义。 5隐马氏模型( h m m ) 简介 在最近的几十年中，随着各类基因组、c d n a 和其它测序计划的发 展，尤其是测序过程中产生的表达序列标签( e x p r e s s e ds e q u e n c et a g ，e s t ) 的积累，大规模片断数据库变得越来越实用。对这些片断进行分类和识 别，并从中挖掘更多的有用信息就变得很有意义。利用多重序列比对提取