隐马尔可夫模型(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再来看看这个算法的实现过程

参考资料：

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