HMM的学习

爱被打了一巴掌 2023-06-13 06:20 32阅读 0赞

### 隐马尔科夫模型 ###

*  隐马尔科夫模型的基本概念
 *  HMM的三个问题
 *   *  概率计算问题
     *   *  前向算法
     *  学习问题
     *  预测问题
     *   *  维特比算法

# 隐马尔科夫模型的基本概念 #

定义10.1（隐马尔可夫模型） ==隐马尔可夫模型==是关于时序的概率模型， 描述由一个
    隐藏的马尔可夫链随机生成不可观测的状态随机序列， 再由各个状态生成一个观测而产生
    观测随机序列的过程。 隐藏的马尔可夫链随机生成的状态的序列， 称为状态序列（state
    sequence） ； 每个状态生成一个观测， 而由此产生的观测的随机序列， 称为观测序列
    （observation sequence） 。 序列的每一个位置又可以看作是一个时刻。

隐马尔可夫模型由初始概率分布、 状态转移概率分布以及观测概率分布确定。 隐马尔  
可夫模型的形式定义如下：  
首先：马尔科夫要求状态和观测的集合是有限的。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70]![在这里插入图片描述][20191114194903505.png]  
其中隐马尔科夫模型有两个假设：  
1、齐次马尔科夫行假设，即假设隐藏的马尔科夫链在任意时刻t的状态只依赖于其前一时刻的状态，于其它时刻的状态及观测无关，也与时刻t无关。  
![在这里插入图片描述][20191114195330577.png]  
2、观测独立性假设，当前时刻的观测只依赖当前时刻的状态。

# HMM的三个问题 #

## 概率计算问题 ##

给定模型λ和观测序列O，计算在模型λ的基础下观测序列O出现的概率。  
需要的算法是前向算法与后向算法。

### 前向算法 ###

![在这里插入图片描述][20191114200831425.png]  
首先直接计算算法，即求出所有的符合观测序列O的状态序列I=（i1,i2,i3…it），求和计算该观测序列在给出模型下的概率。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70 1]  
缺点：计算量太大，存在重复计算问题。  
接下来是前向算法：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70 2]  
![在这里插入图片描述][20191114201515441.png]  
总结：递归计算，依次按照观测序列O中的观测值计算。  
看一个例题就明白了：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70 3]  
后向算法就是反着去推，前向算法从观测序列O的第一个往后去递归，后向算法相反。

## 学习问题 ##

已知观测序列O，学习模型λ的参数，使得在该模型下观测序列概率最大。  
需要的算法是Baum-Welch算法（EM算法）。

## 预测问题 ##

已知模型λ和观测序列O，求最有可能对应的状态序列I。  
主要学习维特比算法（Viterbi algorithm）。

### 维特比算法 ###

维特比算法实际是用动态规划解隐马尔可夫模型预测问题， 即用动态规划（dynamic
    programming）求概率最大路径（最优路径），这时一条路径对应着一个状态序列。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70 4]

例子：

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70 5]  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70 6]  
代码：

import numpy as np
    class HiddenMarkov:
        def forward(self, Q, V, A, B, O, PI):  # 使用前向算法
            N = len(Q)  # 状态序列的大小
            M = len(O)  # 观测序列的大小
            alphas = np.zeros((N, M))  # alpha值
            T = M  # 有几个时刻，有几个观测序列，就有几个时刻
            for t in range(T):  # 遍历每一时刻，算出alpha值
                indexOfO = V.index(O[t])  # 找出序列对应的索引
                for i in range(N):
                    if t == 0:  # 计算初值
                        alphas[i][t] = PI[t][i] * B[i][indexOfO]  # P176（10.15）
                        print('alpha1(%d)=p%db%db(o1)=%f' % (i, i, i, alphas[i][t]))
                    else:
                        alphas[i][t] = np.dot([alpha[t - 1] for alpha in alphas], [a[i] for a in A]) * B[i][
                            indexOfO]  # 对应P176（10.16）
                        print('alpha%d(%d)=[sigma alpha%d(i)ai%d]b%d(o%d)=%f' % (t, i, t - 1, i, i, t, alphas[i][t]))
                        # print(alphas)
            P = np.sum([alpha[M - 1] for alpha in alphas])  # P176(10.17)
            # alpha11 = pi[0][0] * B[0][0]    #代表a1(1)
            # alpha12 = pi[0][1] * B[1][0]    #代表a1(2)
            # alpha13 = pi[0][2] * B[2][0]    #代表a1(3)
    
        def backward(self, Q, V, A, B, O, PI):  # 后向算法
            N = len(Q)  # 状态序列的大小
            M = len(O)  # 观测序列的大小
            betas = np.ones((N, M))  # beta
            for i in range(N):
                print('beta%d(%d)=1' % (M, i))
            for t in range(M - 2, -1, -1):
                indexOfO = V.index(O[t + 1])  # 找出序列对应的索引
                for i in range(N):
                    betas[i][t] = np.dot(np.multiply(A[i], [b[indexOfO] for b in B]), [beta[t + 1] for beta in betas])
                    realT = t + 1
                    realI = i + 1
                    print('beta%d(%d)=[sigma a%djbj(o%d)]beta%d(j)=(' % (realT, realI, realI, realT + 1, realT + 1),
                          end='')
                    for j in range(N):
                        print("%.2f*%.2f*%.2f+" % (A[i][j], B[j][indexOfO], betas[j][t + 1]), end='')
                    print("0)=%.3f" % betas[i][t])
            # print(betas)
            indexOfO = V.index(O[0])
            P = np.dot(np.multiply(PI, [b[indexOfO] for b in B]), [beta[0] for beta in betas])
            print("P(O|lambda)=", end="")
            for i in range(N):
                print("%.1f*%.1f*%.5f+" % (PI[0][i], B[i][indexOfO], betas[i][0]), end="")
            print("0=%f" % P)
    
        def viterbi(self, Q, V, A, B, O, PI):
            N = len(Q)  # 状态序列的大小
            M = len(O)  # 观测序列的大小
            deltas = np.zeros((N, M))
            psis = np.zeros((N, M))
            I = np.zeros((1, M))
            for t in range(M):
                realT = t+1
                indexOfO = V.index(O[t])  # 找出序列对应的索引
                for i in range(N):
                    realI = i+1
                    if t == 0:
                        deltas[i][t] = PI[0][i] * B[i][indexOfO]
                        psis[i][t] = 0
                        print('delta1(%d)=pi%d * b%d(o1)=%.2f * %.2f=%.2f'%(realI, realI, realI, PI[0][i], B[i][indexOfO], deltas[i][t]))
                        print('psis1(%d)=0' % (realI))
                    else:
                        deltas[i][t] = np.max(np.multiply([delta[t-1] for delta in deltas], [a[i] for a in A])) * B[i][indexOfO]
                        print('delta%d(%d)=max[delta%d(j)aj%d]b%d(o%d)=%.2f*%.2f=%.5f'%(realT, realI, realT-1, realI, realI, realT, np.max(np.multiply([delta[t-1] for delta in deltas], [a[i] for a in A])), B[i][indexOfO], deltas[i][t]))
                        psis[i][t] = np.argmax(np.multiply([delta[t-1] for delta in deltas], [a[i] for a in A]))
                        print('psis%d(%d)=argmax[delta%d(j)aj%d]=%d' % (realT, realI, realT-1, realI, psis[i][t]))
            print(deltas)
            print(psis)
            I[0][M-1] = np.argmax([delta[M-1] for delta in deltas])
            print('i%d=argmax[deltaT(i)]=%d' % (M, I[0][M-1]+1))
            for t in range(M-2, -1, -1):
                I[0][t] = psis[int(I[0][t+1])][t+1]
                print('i%d=psis%d(i%d)=%d' % (t+1, t+2, t+2, I[0][t]+1))
            print(I)
    
    if __name__ == '__main__':
        Q = [1, 2, 3]
        V = ['红', '白']
        A = [[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]
        B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]
        # O = ['红', '白', '红', '红', '白', '红', '白', '白']
        O = ['红', '白', '红', '白']    #习题10.1的例子
        PI = [[0.2, 0.4, 0.4]]
        HMM = HiddenMarkov()
        # HMM.forward(Q, V, A, B, O, PI)
        # HMM.backward(Q, V, A, B, O, PI)
        HMM.viterbi(Q, V, A, B, O, PI)

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70]: https://img-blog.csdnimg.cn/20191114194750795.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg==,size_16,color_FFFFFF,t_70
[20191114194903505.png]: https://img-blog.csdnimg.cn/20191114194903505.png
[20191114195330577.png]: https://img-blog.csdnimg.cn/20191114195330577.png
[20191114200831425.png]: https://img-blog.csdnimg.cn/20191114200831425.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70 1]: https://img-blog.csdnimg.cn/20191114201223863.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg==,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70 2]: https://img-blog.csdnimg.cn/20191114201505786.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg==,size_16,color_FFFFFF,t_70
[20191114201515441.png]: https://img-blog.csdnimg.cn/20191114201515441.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70 3]: https://img-blog.csdnimg.cn/20191114201621121.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg==,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70 4]: https://img-blog.csdnimg.cn/20191114205417158.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg==,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70 5]: https://img-blog.csdnimg.cn/20191114205504106.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg==,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg_size_16_color_FFFFFF_t_70 6]: https://img-blog.csdnimg.cn/20191114205516707.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80Mzg0ODYwMg==,size_16,color_FFFFFF,t_70