基因识别及隐马模型

即每次叠代后非完整数据的似然性都会 增加，而似然性显然是有界的，因此叠 代必然收敛。
隐马模型的参数估计
利用EM算法进行隐马模型的参数估计， 又称为Baum-Welch算法。它要解决的主 要问题，是求
W arg max P(OW )
W
隐马模型复习
隐马模型：
– 存在隐序列H，由下列参数生成：
令上式等于0，由于 pm1(1 p)nm1 0，有：
m mpˆ npˆ mpˆ 0
pˆ m
n
即：发芽率p的极大似然估计为
pˆ m n
EM算法的一般概念
X：观测得到的数据，它是不完全的； Z=(X,Y)，是完全的；其中Y是缺失的或隐藏
的。它的联合分布为：
例：基因组编码区的隐马模型
设基因组由两种功能区域组成，即编码 区和非编码区。分别由字母c，n代表。 转移矩阵为同种字母延伸或变为另一种 字母的概率。初始状态概率为第一个字 母出现c或n的概率。明序列由A，C，G， T四个字母组成，生成概率分别为编码区 和非编码区四个字母出现的概率。
隐马氏模型的三种典型问题
随机过程：一族无穷多个、相互有关联的
随机变量。记为：X (t), t T
由于参数 t 经常代表时间，故称为随机过 程。T常为自然数，整数或区间。当参数 取值为整数时，也称为随机序列。
马尔可夫过程：取值为整数的随机过程， 若 t = i 时刻的取值只与时刻 i-1 取值有关， 则称为马尔可夫过程。
显然 iδT1((α))就是[m解{a码x}问 t题(的)t解。] eoi1
学习算法
HMM的学习算法有许多种，包括 EM （期望最大化）算法，不同形式的梯度 下降法，模拟退火算法等。下面我们主 要介绍 EM 算法。
EM算法的用途
EM算法（The expectation maximization algorithm）是根 据不完整数据作最大似然估计 的一般方法。这种不完整可以 是缺失了某些数据，也可以是 存在某些无法观测的隐藏变量。
P(O | w)
通过计算上述概率的极大值，可以得到 位置 i 处隐序列最可能的状态。 但在解码问题中，我们更关心的是最可 能的隐序列。它是一个整体，不能分解 开一位一位算。解决这个问题，需要使 用Viterbi算法。
Viterbi算法
定义 i ( ) max P(o1oi , hi | w) Hi 其中Hi 为一切长为i、以α结束的隐序列。 显然δi(α)对应的隐序列就是最可能生成 前i个字母明序列的隐序列。递推公式为：
{ ,T ,{ }}
– 观测序列O称为明序列，由下列参数生成：
{H , Ea ,{a}}
前向算法
定义在模型w下，隐序列第i位置状态为α， 明序列为 o1…oi时的概率为：
f i ( ) P(hi , o1 oi | w)
则有递推公式：
f i1( )
f i ( )teoi1
解：对于任意一条隐序列H，有：
T
P(O, H w) (t eoi ) i 1
其中α为隐序列H第i-1位置的状态， β为i 位置的状态。当i＝1时，tαβ应由πβ代替。
因此，所求概率为：
T
P(O w) P(O, H w)
(t eoi )
H
H i1
其中H为一切可能的隐序列。
log P(x ) log P(x t ) Q( t ) Q( t t )
当且仅当θt=θ，或P(y|x,θt)=P(y|x,θ)时等 式成立。
若取θt+1=argmaxQ(θ|θt)， 显然有 Q(θt+1|θt)>Q(θt|θt)，因此
log P(x t1) log P(x t )
bi ( )
0 arg max L( X )

为简化计算，实际工作中常使用的是对 数似然估计。即：
0 arg max ln L( X )

具体计算公式为：
d ln L( ) 0 d
如果参数不止一个，相应令偏导数为0即 可。
例题1 正态分布的最大似然估计
设x1 ,x2 ,… xn是取自正态总体 N(μ ,σ 2) 的简单随机子样，μ 与σ 2是未知参数，求
log P(x ) log P(x t ) Q( t ) Q( t t ) P( y x, t ) log P( y x, t )
y
P( y x, )
上式中最后一项是P(y|x,θt)与P(y|x,θ)之间 的Kullback-Leibler距离，它总是大于0 的，因此
可能性问题：给定模型参数，当观察到 一个明序列时，这一明序列确实由给定 模型生成的概率有多大？
解码问题：给定模型参数，当观察到一 个明序列时，这一明序列所对应的最可 能的隐序列是什么？
学习问题：观察到足够多明序列时，如 何估计转移概率和生成概率？
基本算法
可能性问题：已知参数为w的隐马模型和 长为T的序列O＝{o1o2…ot…oT}；求在给 定模型下观测到序列O的概率。
n
ln L( X ) ln P(xi ) i 1
最大似然估计
目标：根据已知观测数据估计总体参数。 似然函数可视为以未知参数为自变量的函
数，它的统计意义又是在未知参数下观测 到该组数据的概率，因此我们很自然地取 使似然函数达到最大值的未知参数为估计 值。这就是最大似然估计。即：
,) P(X )
对于这个联合分布函数，我们可以定义它的似 然函数： L( Z) L( X ,Y) P(X ,Y )
称为完整数据的似然函数。而X的似然函数
L( X ) P(X )
称为非完整数据的似然函数。
注意完整数据似然函数的自变量中，X和 θ为常数，而Y是随机变量，且服从于某种 由X和θ所决定的分布。参数X是观测数据， 是完全确定的；参数θ则是我们需要估计 的，在计算过程中它会不断调整。设θ(i-1) 为当前我们使用的θ估计值，则我们可以 定义完整数据对数似然的期望：
似然函数的概念
我们从某个未知分布得到了一组观察值 X={x1,x2,…xn}。此未知分布是由一组参数 Θ={θ1,θ2, … θm}决定。定义似然函数为：
n
由于对数L(函X数) 是 单i1 调P(递xi 增)的 P，( X因此) 不会改变
极值点的位置；而且它可以把连乘变成连加， 从而大大简化计算。因此实际工作中常使用 的是对数似然函数。即：
μ 和σ 的极大似然估计。
解：由于
f (x, , 2 )
1
2
exp

1
2
2
(x

)2

故有似然函数
L(, 2 )
(
1
2
)n
exp
1
2
2

n i1
(xi

)2

取对数，有：
ln L(, 2 ) n ln(2
2
Q(, (i1) ) E[ln P( X ,Y ) X , (i1) ]
ln P( X ,Y ) f ( y X , (i1) )dy yY
函数Q(θ,θ(i-1) )就是EM算法中，E步骤所 要计算的值。在上述函数中，θ(i-1)是上一 步计算的结果，是确定的值。而θ是接着 的M步骤中要调整的参数，调整它以便 使似然函数的值增加。即：
2
)

1
2
2
n
(xi
i 1
)2
∴ 似然方程为：
ln L(, 2 )



ln
L(,
2
)
2
0 n
2
1
2
2
(2)
n i1
( xi

)

0

2 2
2

1
2 4

n i1
( xi
)2

0
(1) (2)
由（1）解得：
对上式两边同时乘以P(y|x,θt)，并对所有
y求和得到
log P(x ) P( y x, t )log P(x, y ) P( y x, t )log P( y x, )
y
y
Q( t ) P( y x, t ) log P( y x, )
y
因此有：
(i) arg max Q(, (i1) )

在EM算法中，上述E步骤和M步骤不断 重复，每一次重复似然值都会增加，直 到收敛到似然函数的局部极大值。
EM算法收敛性的证明
由于 P(x, y) P(y x,)P(x) 对数似然函数可以改写成： log P(x) log P(x, y) log P(y x,)
生物信息学中的典型数学问题
基因识别 与隐马氏模型
基因结构及重要信号
基因识别模型
问题：从基因组序列出发，识别其中的 基因（基因组注释）。
输入：基因组序列 输出：各种基因元件（包括启动子，外
显子，内含子等）在基因组上的位置。 如果需要，也可输出翻译的蛋白质序列， 预测结果的可靠程度等。
基本概念
由于H的数量随T的增加以指数增长， 这种算法实际并不可行。
前向算法
定义在模型w下，隐序列第i位置状态为α， 明序列为 o1…oi时的概率为：
f i ( ) P(hi , o1 oi | w)
则有递推公式：
f i1()
f i ( )teoi1
{ }
初值为：
n
(xi ) 0
i 1
ˆ

1 n
n i 1
xi

x
代入（2），得：
1
2
n
(xi x)2
i 1
n
ˆ 2

1 n
n i 1
(xi

x)2

S
2 n
即：μ 和σ 2的极大似然估计分别
为
x
和
S
2 n
。
例题2 二项分布的最大似然估计
取n粒种子作发芽试验，其中有m粒发芽， 求发芽率p的最大似然估计。
f 1( ) eo1
则有：
P(O | w) f T () { }
后向算法
定义在模型w，隐序列第i位置状态为α条 件下，明序列为 oi＋1…oT时的概率为：
bi ( ) P(oi1 oT | hi , w)
则有递推公式：
bi ( ) bi1( )t eoi1 { }
{ }
初值为：
f 1( ) eo1
则有：
P(O | w) f T () { }
后向算法
定义在模型w，隐序列第i位置状态为α条 件下，明序列为 oi＋1…oT时的概率为：
bi ( ) P(oi1 oT | hi , w) 则有递推公式：
初值为：
bT ( ) 1
则有：
P(O | w) b1( ) eo1 { }
使用前向或后向算法，计算P(O|w)的复 杂度约为O(N2)级。
联合应用前向和后向算法，可以容易地 计算隐序列位置i处于状态α的概率：
pi () P(hi | O, w) f i ()bi ()
隐马氏模型（HMM）
模型的数学描述：
– 存在一个隐序列H，它是不可观测的，且由以 下参数生成： { ,T ,{ }}
其中πα为初始状态出现概率；Tαβ为转移概率， 即tαβ = P(hi= β | hi-1= α)； α, β属于{σ}； {σ}为 字符集，即隐序列由哪些字符组成。
– 观测的结果称为明序列O，它由隐序列按照生 成概率e αa生成。其中eαa =P(a | α)；{a}为明序 列字符集。
解：每粒种子发芽与否可视为两点分布： 发芽， 则 X=1, 其概率为 p 不发芽，则 X=0, 其概率为1-p 由似然函数的构造，有： L(p)=P(X=x1|p)·P(X=x2|p)…P(X=xn|p)
由于共有m粒发芽，(n-m)粒不发芽， ∴ L(p)=pm(1-p)n-m
dL( p) mpm1 (1 p)nm (n m) (1) p m (1 p)nm1 dp p m1 (1 p)nm1[m(1 p) (n m) p]