隐马尔可夫模型（Hidden Markov Model,HMM）—有监督学习方法、概率模型、生成模型

标签：状态概率模型模型 len HMM range line PI sum

定义

隐马尔可夫模型是关于时序的概率模型,描述由一个隐藏的马尔可夫链随机生成不可观测的状态随机序列,再由各个状态生成一个观测从而产生观测随机序列的过程。隐藏的马尔可夫链随机生成的状态的序列,称为状态序列(State Sequence);每个状态生成一个观测,而由此产生的观测的随机序列,称为观测序列(Observation Sequence)。序列的每一位置又可以看作是一个时刻。

基本问题

(1)概率计算问题。给定模型 λ = ( A , B , π ) \lambda=(A,B,\pi) λ=(A,B,π)和观测序列 O = ( o 1 , o 2 , ⋯ , o T ) O=(o_1,o_2,\cdots,o_T) O=(o1,o2,⋯,oT),计算在模型 λ \lambda λ下观测序列 O O O出现的概率 P ( O ∣ λ ) P(O|\lambda) P(O∣λ)。
(2)学习问题。已知观测序列 O = ( o 1 , o 2 , ⋯ , o T ) O=(o_1,o_2,\cdots,o_T) O=(o1,o2,⋯,oT),估计模型 λ = ( A , B , π ) \lambda=(A,B,\pi) λ=(A,B,π)参数,使得在该模型下观测序列概率 P ( O ∣ λ ) P(O|\lambda) P(O∣λ)最大。即用极大似然估计的方法估计参数。
(3)预测问题,也称为解码(decoding)。已知模型 λ = ( A , B , π ) \lambda=(A,B,\pi) λ=(A,B,π)和观测序列 O = ( o 1 , o 2 , ⋯ , o T ) O=(o_1,o_2,\cdots,o_T) O=(o1,o2,⋯,oT)，求对给定观测序列条件概率 P ( I ∣ O ) P(I|O) P(I∣O)最大的状态序列 I = ( i 1 , i 2 , ⋯ , i T ) I=(i_1,i_2,\cdots,i_T) I=(i1,i2,⋯,iT)。即给定观测序列，求最有可能的对应的状态序列。

输入空间

O = ( o 1 , o 2 , ⋯ , o T ) O=(o_1,o_2,\cdots,o_T) O=(o1,o2,⋯,oT)

特征空间（Feature Space）

O = ( o 1 , o 2 , ⋯ , o T ) O=(o_1,o_2,\cdots,o_T) O=(o1,o2,⋯,oT)

统计学习方法

模型

I ∗ = f ( λ , O ) I^* = f(\lambda,O) I∗=f(λ,O)

策略

P ∗ = m a x 1 ≤ i ≤ N δ T ( i ) P^* = \mathop{max}\limits_{1\leq i \leq N} \delta_T(i) P∗=1≤i≤NmaxδT(i)
i T ∗ = a r g ∗ m a x 1 ≤ i ≤ N [ δ T ( i ) ] i_T^* = arg * \mathop{max}\limits_{1\leq i \leq N} \big[ \delta_T(i) \big] iT∗=arg∗1≤i≤Nmax[δT(i)]

算法

α t + 1 ( i ) = [ ∑ j = 1 N α t ( j ) a j i ] b i ( o t + 1 ) , i = 1 , 2 , ⋯ , N \alpha_{t+1}(i)=\bigg[ \sum_{j=1}^N \alpha_t(j) a_{ji} \bigg]b_i(o_{t+1}),i=1,2,\cdots,N αt+1(i)=[j=1∑Nαt(j)aji]bi(ot+1),i=1,2,⋯,N
β t ( i ) = ∑ j = 1 N a i j b j ( o t + 1 ) β t + 1 ( j ) , 1 = 1 , 2 , ⋯ , N \beta_{t}(i)=\sum_{j=1}^Na_{ij}b_j(o_{t+1})\beta_{t+1}(j),1=1,2,\cdots,N βt(i)=j=1∑Naijbj(ot+1)βt+1(j),1=1,2,⋯,N
a i j ( n + 1 ) = ∑ t = 1 T − 1 ζ t ( i , j ) ∑ t = 1 T − 1 γ t ( i ) , 其中 ζ t ( i , j ) = α t ( i ) a i j b j ( o t + 1 ) β t + 1 ( j ) ∑ t = 1 N ∑ j = 1 N α t ( i ) a i j b j ( o t + 1 ) β t + 1 ( j ) , γ t ( i ) = α t ( i ) β t ( i ) ∑ j = 1 N α t ( j ) β t ( j ) a_{ij}^{(n+1)}=\frac{\sum_{t=1}^{T-1}\zeta_t(i,j)}{\sum_{t=1}^{T-1}\gamma_t(i)},其中\zeta_t(i,j)=\frac{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{\sum_{t=1}^N\sum_{j=1}^N\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)},\gamma_t(i)=\frac{\alpha_t(i)\beta_t(i)}{\sum_{j=1}^N\alpha_t(j)\beta_t(j)} aij(n+1)=∑t=1T−1γt(i)∑t=1T−1ζt(i,j),其中ζt(i,j)=∑t=1N∑j=1Nαt(i)aijbj(ot+1)βt+1(j)αt(i)aijbj(ot+1)βt+1(j),γt(i)=∑j=1Nαt(j)βt(j)αt(i)βt(i)
b j ( k ) ( n + 1 ) = ∑ t = 1 , o t = v k T γ t ( j ) ∑ t = 1 T γ t ( j ) b_j(k)^{(n+1)}=\frac{\sum_{t=1,o_t=v_k}^T\gamma_t(j)}{\sum_{t=1}^T\gamma_t(j)} bj(k)(n+1)=∑t=1Tγt(j)∑t=1,ot=vkTγt(j)
π i ( n + 1 ) = γ 1 ( i ) \pi_i^{(n+1)} = \gamma_1(i) πi(n+1)=γ1(i)

import numpy as np
import time

def trainParameter(fileName):
    '''
    数据集 下载地址：https://download.csdn.net/download/nanxiaotao/89738993
    依据训练文本统计PI、A、B
    :param fileName: 训练文本
    :return: 三个参数
    '''
    #定义一个查询字典，用于映射四种标记在数组中对应的位置，方便查询
    # B：词语的开头
    # M：一个词语的中间词
    # E：一个词语的结果
    # S：非词语，单个词
    statuDict = {'B':0, 'M':1, 'E':2, 'S':3}

    #每个字只有四种状态，所以下方的各类初始化中大小的参数均为4
    #初始化PI的一维数组，因为对应四种状态，大小为4
    PI = np.zeros(4)
    #初始化状态转移矩阵A，涉及到四种状态各自到四种状态的转移，因为大小为4x4
    A = np.zeros((4, 4))
    #初始化观测概率矩阵，分别为四种状态到每个字的发射概率
    #因为是中文分词，使用ord(汉字)即可找到其对应编码，这里用一个65536的空间来保证对于所有的汉字都能
    #找到对应的位置来存储
    B = np.zeros((4, 65536))
    #去读训练文本
    fr = open(fileName, encoding='utf-8')

    #文本中的每一行认为是一个训练样本
    #注：并没有使用Baum-Welch算法，只是借助了其内部的三个参数生成公式，其实
    #公式并不是Baum-Welch特有的，只是在那一节正好有描述
    for line in fr.readlines():
        #对单行句子按空格进行切割
        curLine = line.strip().split()
        #对词性的标记放在该列表中
        wordLabel = []
        #对每一个单词进行遍历
        for i in range(len(curLine)):
            #如果长度为1，则直接将该字标记为S，即单个词
            if len(curLine[i]) == 1:
                label = 'S'
            else:
                #如果长度不为1，开头为B，最后为E，中间添加长度-2个M
                #如果长度刚好为2，长度-2=0也就不添加了，反之添加对应个数的M
                label = 'B' + 'M' * (len(curLine[i]) - 2) + 'E'

            #如果是单行开头第一个字，PI中对应位置加1,
            if i == 0: PI[statuDict[label[0]]] += 1

            #对于该单词中的每一个字，在生成的状态链中统计B
            for j in range(len(label)):
                #遍历状态链中每一个状态，并找到对应的中文汉字，在B中
                #对应位置加1
                B[statuDict[label[j]]][ord(curLine[i][j])] += 1

            #在整行的状态链中添加该单词的状态链
            #注意：extend表直接在原先元素的后方添加，
            wordLabel.extend(label)

        #单行所有单词都结束后，统计A信息
        #因为A涉及到前一个状态，因此需要等整条状态链都生成了才能开始统计
        for i in range(1, len(wordLabel)):
            #统计t时刻状态和t-1时刻状态的所有状态组合的出现次数
            A[statuDict[wordLabel[i - 1]]][statuDict[wordLabel[i]]] += 1

    #对PI求和，概率生成中的分母
    sum = np.sum(PI)
    #遍历PI中每一个元素，元素出现的次数/总次数即为概率
    for i in range(len(PI)):
        if PI[i] == 0:  PI[i] = -3.14e+100
        #如果不为0，则计算概率，再对概率求log
        else: PI[i] = np.log(PI[i] / sum)

    #与上方PI思路一样，求得A的概率对数
    for i in range(len(A)):
        sum = np.sum(A[i])
        for j in range(len(A[i])):
            if A[i][j] == 0: A[i][j] = -3.14e+100
            else: A[i][j] = np.log(A[i][j] / sum)

    #与上方PI思路一样，求得B的概率对数
    for i in range(len(B)):
        sum = np.sum(B[i])
        for j in range(len(B[i])):
            if B[i][j] == 0: B[i][j] = -3.14e+100
            else:B[i][j] = np.log(B[i][j] / sum)

    #返回统计得到的三个参数
    return PI, A, B

#依据现有训练集统计PI、A、B
PI, A, B = trainParameter('HMMTrainSet.txt')

def loadArticle(fileName):
    '''
    下载地址：https://download.csdn.net/download/nanxiaotao/89738993
    加载文章
    :param fileName:文件路径
    :return: 文章内容
    '''
    #初始化文章列表
    artical = []
    #打开文件
    fr = open(fileName, encoding='utf-8')
    #按行读取文件
    for line in fr.readlines():
        #读到的每行最后都有一个\n，使用strip将最后的回车符去掉
        line = line.strip()
        #将该行放入文章列表中
        artical.append(line)
    #将文章返回
    return artical

(1)初始化
δ 1 ( i ) = π i b i ( o 1 ) , i = 1 , 2 , ⋯ , N \delta_1(i) = \pi_ib_i(o_1),i=1,2,\cdots,N δ1(i)=πibi(o1),i=1,2,⋯,N
ψ 1 ( i ) = 0 , i = 1 , 2 , ⋯ , N \psi_1(i)=0,i=1,2,\cdots,N ψ1(i)=0,i=1,2,⋯,N
(2)递推
δ t ( i ) = m a x 1 ≤ i ≤ N [ δ t − 1 ( j ) a ( j i ) ] b i ( o t ) , i = 1 , 2 , ⋯ , N \delta_t(i) = \mathop{max}\limits_{1\leq i \leq N}[\delta_{t-1}(j)a_(ji)]b_i(o_t),i=1,2,\cdots,N δt(i)=1≤i≤Nmax[δt−1(j)a(ji)]bi(ot),i=1,2,⋯,N
ψ t ( i ) = a r g ∗ m a x 1 ≤ i ≤ N [ δ t − 1 ( j ) a j i ] , i = 1 , 2 , ⋯ , N \psi_t(i) = arg * \mathop{max}\limits_{1\leq i \leq N}[\delta_{t-1}(j)a_{ji}],i = 1,2,\cdots,N ψt(i)=arg∗1≤i≤Nmax[δt−1(j)aji],i=1,2,⋯,N
(3)终止
P ∗ = m a x 1 ≤ i ≤ N δ T ( i ) P^* = \mathop{max}\limits_{1\leq i \leq N} \delta_T(i) P∗=1≤i≤NmaxδT(i)
i T ∗ = a r g ∗ m a x 1 ≤ i ≤ N [ δ T ( i ) ] i_T^* = arg * \mathop{max}\limits_{1\leq i \leq N} \big[ \delta_T(i) \big] iT∗=arg∗1≤i≤Nmax[δT(i)]

def participle(artical, PI, A, B):
    '''
    分词”
    :param artical:要分词的文章
    :param PI: 初始状态概率向量PI
    :param A: 状态转移矩阵
    :param B: 观测概率矩阵
    :return: 分词后的文章
    '''
    #初始化分词后的文章列表
    retArtical = []

    #对文章按行读取
    for line in artical:
        #四种概率值，因此长度时该行的长度
        delta = [[0 for i in range(4)] for i in range(len(line))]
        for i in range(4):
            #初始化δ状态链中第一个状态的四种状态概率
            delta[0][i] = PI[i] + B[i][ord(line[0])]
        #初始化ψ，初始时为0
        psi = [[0 for i in range(4)] for i in range(len(line))]

        #依次处理整条链
        for t in range(1, len(line)):
            #对于链中的米格状态，求四种状态概率
            for i in range(4):
                #初始化一个临时列表，用于存放四种概率
                tmpDelta = [0] * 4
                for j in range(4):
                    tmpDelta[j] = delta[t - 1][j] + A[j][i]

                #找到最大的那个δ * a，
                maxDelta = max(tmpDelta)
                #记录最大值对应的状态
                maxDeltaIndex = tmpDelta.index(maxDelta)

                #将找到的最大值乘以b放入，
                #注意：这里同样因为log变成了加法
                delta[t][i] = maxDelta + B[i][ord(line[t])]
                #在ψ中记录对应的最大状态索引
                psi[t][i] = maxDeltaIndex

        #建立一个状态链列表，开始生成状态链
        sequence = []
        #在上面for循环全部结束后，很明显就到了第三步了
        #获取最后一个状态的最大状态概率对应的索引
        i_opt = delta[len(line) - 1].index(max(delta[len(line) - 1]))
        #在状态链中添加索引
        #注：状态链应该是B、M、E、S，这里图方便用了0、1、2、3，其实一样的
        sequence.append(i_opt)
        #算法10.5 第四步：最优路径回溯
        #从后往前遍历整条链
        for t in range(len(line) - 1, 0, -1):
            #不断地从当前时刻t的ψ列表中读取到t-1的最优状态
            i_opt = psi[t][i_opt]
            #将状态放入列表中
            sequence.append(i_opt)
        #因为是从后往前将状态放入的列表，所以这里需要翻转一下，变成了从前往后
        sequence.reverse()

        #开始对该行分词
        curLine = ''
        #遍历该行每一个字
        for i in range(len(line)):
            #在列表中放入该字
            curLine += line[i]
            #如果该字是3：S->单个词  或  2:E->结尾词 ，则在该字后面加上分隔符 |
            #此外如果改行的最后一个字了，也就不需要加 |
            if (sequence[i] == 3 or sequence[i] == 2) and i != (len(line) - 1):
                curLine += '|'
        #在返回列表中添加分词后的该行
        retArtical.append(curLine)
    #返回分词后的文章
    return retArtical

#读取测试文章
artical = loadArticle('testArtical.txt')

#打印原文
print('-------------------原文----------------------')
for line in artical:
    print(line)

#进行分词
partiArtical = participle(artical, PI, A, B)

#打印分词结果
print('-------------------分词后----------------------')
for line in partiArtical:
    print(line)

假设空间（Hypothesis Space）

{ f ∣ f ( x ) = ( i 1 ∗ , i 2 ∗ , ⋯ , i T ∗ ) } \left\{f|f(x) = (i_1^*,i_2^*,\cdots,i_T^*) \right\} {f∣f(x)=(i1∗,i2∗,⋯,iT∗)}

输出空间

( i 1 ∗ , i 2 ∗ , ⋯ , i T ∗ ) (i_1^*,i_2^*,\cdots,i_T^*) (i1∗,i2∗,⋯,iT∗)

标签：状态,概率模型,模型,len,HMM,range,line,PI,sum
From： https://blog.csdn.net/nanxiaotao/article/details/142145961