0%

统计学习方法-隐马尔可夫模型(HMM)笔记

最近在看到CNV的时候,发现一些检测CNV的分析方法都是基于隐马尔可夫模型(HMM),而HMM在生物信息学中也应用广泛。如早期的DNA序列分析(对一家族序列建立HMM模型),还有预测DNA编码区域(对已知或者指定序列进行训练从而构建HMM模型),以及CNV分析(在Segmentation中利用HMM模型预估CNV数目,确定相应区域的拷贝数)等等

HMM是生物信息学中比较流行的机器学习和模式识别方法,它具有对模型中的一些隐性参数进行识别和优化的能力,使得这种模型具有很强的鲁棒性,并且可以随着训练深入进一步提高识别精度,同时这种自适应的模型能够有效地实现检测过程中的参数优化

因此整理下<<统计学习方法>>的隐马尔可夫模型和一些网上资料 基本概念:

  • 隐马尔可夫模型是一个关于时序的概率模型,可以用于根据一些已知的来推断未知的东西
  • 马尔可夫链是一个随机过程模型,服从马尔可夫性质:无记忆性,某一时刻的状态只受前一个时刻的影响
  • 状态序列是由马尔可夫链随机生成的,是隐藏的,不可观测的;每个状态又有其对应的观测结果,这些观测结果根据时序排序后即为观测序列;也可以说观测序列长度跟时刻长度是一致的

HMM基本假设:

  • 假设隐藏的马尔可夫链在任意时刻的状态只依赖于前一个时刻(t−1时)的状态,与其他时刻的状态及观测无关,也与时刻t无关
  • 假设任意时刻的观测只依赖于该时刻的马尔可夫链状态,与其它观测及状态无关

因此组成HMM有两个空间(状态空间Q和观测空间V),三组参数(状态转移矩阵A,观测概率矩阵以及状态初始概率分布π),有了这些参数,一个HMM就被确定了

HMM模型的三个基本问题:

  • 评估观测序列的概率,在给定HMM模型λ=(A,B,π)和观测序列O下,计算模型λ下观测序列的出现概率P(O|λ),用到的是前向-后向算法
  • 预测(解码)问题,在给定HMM模型λ=(A,B,π)和观测序列O下,求解观测序列的条件概率最大的条件下,最可能出现的状态序列,用到的是维特比算法
  • 模型参数学习问题,在给定观测序列O下,估计模型λ的参数,使得该模型下观测序列的条件概率P(O|λ)最大,用到的是EM算法的鲍姆-韦尔奇算法

向前/向后算法

根据<<统计学习方法>>中的第10章隐马尔可夫模型的例子,记录下向前/向后算法对于HMM第一个基本问题的实现过程,Python代码如下:

class Hidden_Markov_Model:
    def forward(self, A, B, PI, Q, V, O):
        Q_n = len(Q)
        O_n = len(O)
        alphas = np.zeros((O_n, Q_n))
        for t in range(O_n):
            index = V.index(O[t])
            if t == 0:
                alphas[t] = PI * B[:,index]
            else:
                alphas[t] = np.dot(alphas[t-1], A) * B[:,index]
        P = np.sum(alphas[O_n-1])
        print("P(O|lambda) =", P)
        return alphas

    def backward(self, A, B, PI, Q, V, O):
        Q_n = len(Q)
        O_n = len(O)
        betas = np.ones((O_n, Q_n))
        for t in range(O_n-2,-1,-1):
            index = V.index(O[t+1])
            betas[t] = np.dot(A * B[:,index], betas[t+1])
        P = np.dot(PI * B[:,V.index(O[0])], betas[0])
        print("P(O|lambda) =", P)
        return betas

对于10.1习题例子(主要根据前向和后向算法计算结果):

A = np.array([[0.5,0.2,0.3],[0.3,0.5,0.2],[0.2,0.3,0.5]])
B = np.array([[0.5,0.5],[0.4,0.6],[0.7,0.3]])
PI = np.array([0.2,0.4,0.4])
V = ["red","white"]
Q = [1,2,3]
O = ["red","white","red","white"]

HMM = Hidden_Markov_Model()
HMM = Hidden_Markov_Model()
print("P(O|lambda) =", HMM.forward(A,B,PI,Q,V,O))
>>前向 P(O|lambda) =  0.06009079999999999
print("P(O|lambda) =", HMM.backward(A,B,PI,Q,V,O))
>>后向 P(O|lambda) =  0.06009079999999999

对于10.2习题例子(根据公式10.24公式,加上前向和后向概率矩阵,从而计算结果):

A = np.array([[0.5,0.1,0.4],[0.3,0.5,0.2],[0.2,0.2,0.6]])
B = np.array([[0.5,0.5],[0.4,0.6],[0.7,0.3]])
PI = np.array([0.2,0.3,0.5])
V = ["red","white"]
Q = [1,2,3]
O = ["red","white","red","red","white","red","white","white"]

HMM = Hidden_Markov_Model()
alpha = HMM.forward(A,B,PI,Q,V,O)
beta = HMM.backward(A,B,PI,Q,V,O)
p = alpha[3,2]*beta[3,2]/sum(alpha[3] * beta[3])
print("P(i4=q3|O,lambda) = ", p)
>>P(i4=q3|O,lambda) =  0.5369518160647323

维特比算法

对于这个HMM预测问题,最容易想到的是枚举法,如果有3个状态和3个观测序列,那么枚举27次,计算其概率,找出概率最大的那个状态序列即可,但是这种方法运算过大(复杂度是指数增长),所以才有了维特比(Viterbi)算法

维特比算法实际上是用动态规划思路来解决隐马尔可夫的预测问题,相当于是寻找最优路径(概率最大的路径),这路径则对应一个时序状态序列;根据马尔可夫性质,当前状态的概率只跟前一次有关(第t+1时刻的状态概率仅由第t时刻的状态决定),因此依靠前向迭代获取达到每个状态的所有路径的最大概率,然后再反向搜索获取最优路径

因此其还是用前向算法计算了每个观测序列下每个状态的概率,而不是每步只取最好的(A*算法),比如在观测2下状态2是最大概率,那么在计算观测3下哪个状态的概率最大时,不是仅仅考虑之前观测序列的状态2,而是将状态1和3也加入考虑的,代码如下:

class Hidden_Markov_Model:
    def vertibi(self, A, B, PI, Q, V, O):
        Q_n = len(Q)
        O_n = len(O)
        deltas = np.zeros((O_n, Q_n))
        psis = np.zeros((O_n, Q_n))
        I = np.zeros(O_n)
        for t in range(O_n):
            index = V.index(O[t])
            for i in range(Q_n):
                if t == 0:
                    deltas[t,i] = PI[i] * B[i,index]
                    psis[t,i] = 0
                else:
                    deltas[t,i] = np.max(A[:,i] * deltas[t-1,:]) * B[i,index]
                    psis[t,i] = np.argmax(A[:,i] * deltas[t-1,:]) + 1
        I[O_n-1] = np.argmax(deltas[O_n-1,:]) + 1
        for t in range(O_n-2,-1,-1):
            I[t] = psis[t+1,int(I[t+1])-1]
        print("I =", I)
        return

以习题10-3为例:

A = np.array([[0.5,0.2,0.3],[0.3,0.5,0.2],[0.2,0.3,0.5]])
B = np.array([[0.5,0.5],[0.4,0.6],[0.7,0.3]])
PI = np.array([0.2,0.4,0.4])
V = ["red","white"]
Q = [1,2,3]
O = ["red","white","red","white"]

HMM = Hidden_Markov_Model()
HMM.vertibi(A,B,PI,Q,V,O)
>>I = [3. 2. 2. 2.]

除了上述两个HMM问题外,还有一个HMM训练(学习)问题;如果训练数据集既有观测序列又有对应的状态序列,那么则用监督学习,不然如果只有观测序列,那么则是非监督学习;监督学习比较简单,相当于以频数估计状态转移概率和观测概率以及初始状态概率;而非监督学习则是用Baum-Welch(鲍姆-韦尔奇)算法,等看完EM再来看看这个算法的实现过程


参考资料:

HMM隐马尔可夫模型
统计学习方法:习题笔记
第10章 隐马尔可夫模型(HMM) 隐马尔科夫模型HMM(一)HMM模型
隐马尔可夫模型HMM

本文出自于http://www.bioinfo-scrounger.com转载请注明出处