HMM隐马尔科夫模型

马尔科夫过程

在概率论及统计学中，马尔可夫过程（英语：Markov process）是一个具备了马尔可夫性质的随机过程，因为俄国数学家安德雷·马尔可夫得名。马尔可夫过程是不具备记忆特质的。换言之，马尔可夫过程的条件概率仅仅与系统的当前状态相关，而与它的过去历史或未来状态，都是独立、不相关的。

一个马尔科夫过程是状态间的转移仅依赖于前n个状态的过程。这个过程被称之为n阶马尔科夫模型，其中n是影响下一个状态选择的（前）n个状态。最简单的马尔科夫过程是一阶模型，它的状态选择仅与前一个状态有关。这里要注意它与确定性系统并不相同，因为下一个状态的选择由相应的概率决定，并不是确定性的。

马尔可夫链（Markov Chain），描述了一种状态序列，其每个状态值取决于前面有限个状态。马尔可夫链是具有马尔可夫性质的随机变量的一个数列。这些变量的范围，即它们所有可能取值的集合，被称为“状态空间”，而Xn的值则是在时间n的状态。若Xn=i，就说过程在n时刻处于状态i，假设每当过程处于状态i，则过程在下一时
刻处于状态j的概率aij为一定值。即对任意n≥1，Xn+1对于过去状态的条件概率分布仅是Xn的一个函数

这样的随机过程称为Markov链。

状态转移矩阵

对于有N个状态的一阶马尔科夫模型，共有N*N个状态转移，因为任何一个状态都有可能是所有状态的下一个转移状态。每一个状态转移都有一个概率值，称为状态转移概率——这是从一个状态转移到另一个状态的概率。所有的N*N个概率可以用一个状态转移矩阵表示。

定义下面的矩阵A为状态转移概率矩阵：

满足条件：

   1=<i,j<=N.

注意这些概率并不随时间变化而不同——这是一个非常重要（但常常不符合实际）的假设。

隐马尔科夫模型HMM

隐马尔可夫模型（Hidden Markov Model，HMM）是统计模型，它用来描述一个含有隐含未知参数的马尔可夫过程。其难点是从可观察的参数中确定该过程的隐含参数。然后利用这些参数来作进一步的分析，例如模式识别。
在正常的马尔可夫模型中，状态对于观察者来说是直接可见的。这样状态的转换概率便是全部的参数。而在隐马尔可夫模型中,状态并不是直接可见的，但受状态影响的某些变量则是可见的。每一个状态在可能输出的符号上都有一概率分布。因此输出符号的序列能够透露出状态序列的一些信息。

HMM是一种用参数表示的用于描述随机过程统计特性的概率模型，它是一个双重随机过程。HMM 由两部分组成:马尔可夫链和一般随机过程。其中马尔可夫链用来描述状态的转移，用转移概率描述。一般随机过程用来描述状态与观察序列间的关系，用观察值概率描述。对于HMM模型，状态转换过程是不可观察的，因而称之为“隐”马尔可夫模型。

一个HMM的例子

设有N个罐，每个罐子中有很多彩色的球，球的颜色由一组概率分布矢量来描述。实验是这样进行的，根据某个初始概率分布，随机选取N个罐予中的一个，例如第i个罐子，再根据这个罐子中彩色球颜色的概率分布，随机的选择一个球，记下球的颜色，记为O1，再把球放回罐子中，又根据描述罐子的转移概率分布随机地选择下一个罐子，例如第j个罐子，再从罐子中随机选一个球，记下球的颜色，迎为O2，一直进行下去，可以得到一个描述球的颜色序列O1，O2，⋯，由于这是一些观测到的事件，因而称之为观测值序列。但罐子之间的转移以及每次选取的罐子都被隐藏起来了，并不能直接观测到。而且从每个罐子中选取球的颜色并不是与罐子一一对应，而是由该罐子中球的颜色的概率分布随机决定的。此外，每次选择哪一个罐子由一组转移概率所决定。

一个HMM可以由下列参数描述：

1.   N：模型中，Markov链的状态数目。记N个状态为θ1，θ2，⋯，θN，记t时刻Markov链所处的状态为qt，显然qt∈(θ1，θ2，⋯，θN)。在罐子和球的实验中罐子就相当于HMM中的状态。

2．M：每个状态对应的可能的观测值数目。记M个观测值为v1，v2，⋯，vM，记t时刻的观测值为Ot，其中Ot∈(v1，v2，⋯，vM)。在罐子和球的实验中所选球的颜色就是HMM模型中的观测值。

3． π：初始概率分布矢量 π =( π1， π2，⋯， πN)，其中
 πi=P(qt=θi)，1=<i<=N.
在罐子和球的实验中，指实验开始时选择的某个罐子的概率。

4．A：状态转移概率矩阵，A=(aij)N*N。其中，
aij=P(qt+1=θj／qt=θi)              1≤i，j≤N 
在罐子和球的实验中是指每次选取当前罐子的条件下，选取下一个罐子的概率。

5．B：观测值概率矩阵，B=(bjk)N*M，其中
bjk=P(Ot=vk／qt=θj)，1≤j≤N，1≤k≤M
在罐子和球的实验中，bjk就是第j个罐子中球的颜色k出现的概率。
这样，记HMM为：
 λ =(N，M，π，A，B) 
或简写为：
λ =(π，A，B) 

HMM可以分为两部分，一个是Markov链，由π、A来描述，产生的输出为状态序列：另一个随机过程，由B来述，产生的输出为观测值序列。

HMM模型中存在两个假设：一是输出观察值之间严格独立，二是状态的转移过程中当前状态只与前一状态有关(一阶马尔可夫模型)。

HMM基本算法

给定的一个HMM，必需要解决三个基本问题，围绕着三个基本问题，人们研究出了三个基本算法，这三个问题是：

问题1 HMM的概率计算问题
给定观测序列o={ol，o2，⋯，ot}和模型λ，怎样有效地计算观测变量序列o的在给定模型下的概率P(OIλ)

问题2 HMM的最优状态序列问题
给定观测序列O={o1，o2，⋯，ot)和模型λ，怎样选择一个相应的状态序列q={q1，q2，⋯，qt)，能够在某种意义上最优(例如，更好的解释“观测变量”)?

问题3：HMM的训练问题(参数估计问题)

给定观测序列O={o1，o2，⋯，ot，)和初始条件，如何调整模型参数λ=(π，A，B)使P(Olλ)为最大?

前向一后向算法

这个算法是上述第一个问题的解决方案。

前向算法

定义前向变量为：
αt(i)=P(o1，o2，…，ot，qt=θiIλ)，l≤t≤T
那么，有，
a)初始化
α1(i)=πi*bi(o1)
b)递归

c)终结

推导过程

后向算法

定义后向变量为：
βt(i)=P(ot+I，ot+2，…，oT|qt=θi，λ)，1≤t≤T-1
其中，βT(i)=1，后向算法的计算过程如下：
a)初始化
βT(i)=1，1≤i≤N
b)递归

c)终结

后向算法初始化时对于所有的状态i定义βT(i)=1。

Viterbi算法

问题2是寻求“最优”状态序列。‘‘最优’’的意义有很多种，由不同的定义可得到不同的结论。这里讨论的最优意义上的状态序列Q*，是指使P(Q|O，λ)最大时所确定的状态序列。这实际上是一个动态规划问题。这个过程可以用Viterbi算法来实现。Viterbi算法可以叙述如下：
定义δt(i)为t时刻沿一条路径ql，q2，⋯，qt，且qt=θi，产生出o1,o2,⋯,ot的最大概率，即有

那么求取最优状态序列Q*的过程为：
a)初始化

b)递归

C)终结

其中符号argmax符号定义为，如果i=I，f(i)达到最大值，那么

d)求最佳状态序列

Baum—Welch算法

这个算法解决了上面提到的问题三，即HMM训练或参数估计问题，即给定一个观测值序列O={o1，o2，⋯，oT}，该算法能够确定一个模型λ=(π，A，B)使P(Olλ)为最大。这是一个泛函极值问题，因而不存在一个最佳方案来估计λ。在这种情况下，Baum—Welch算法利用递归的思想，使P(Olλ)为局部最大，最后得到模型的参数。其实就是用了EM算法的思想。
定义ξt(i,j)为给定训练序列O和模型λ时，时刻t时Markov链处于θi状态和时刻t+1时处于θj状态时的概率，即

ξt(i,j)=P(O，qt=θi，qt+1=θj|λ)
根据前向变量和后向变量的定义可以导出：
ξt(i,j)=[αt(i)aijbj(ot+1)βt+1(j)]/P(Olλ)
那么t时刻Markov链处于θi状态的概率为：

于是．表示从θi状态转移出去的期望值数目，而表示从θi状态转移到θj状态的期望值数目。由此导出了Baum-Welch算法中著名的重估公式

那么HMM的参数λ=(π，A，B)的求取过程为；根据观测序列O和选取的初始模型λ0=(π，A，B)，由重估公式求得一组新的参数露亦即得到了一个新模型，可以证明。即有重估公式得到的比λ在表示观测值序列O的方面要好，那么重复这个过程，逐步改进模型参数，直到满足一定的收敛条件，即不再明显增大，此时的就是所求之模型，

HMM用于模式识别

HMM是一种动态模式识别工具，能够对一个时间跨度上的信息进行统计建模和分类。

在语音识别中,要首先建立一种对应关系,例如,使一个字对应一个HMM,这里的状态就是指这个音所包含的全部可能的音素(或其细分,或其组合)。对应于此字的一个观测样本,这些音素按照一定的先后顺序出现,这就形成了HMM中的状态序列,现实中不可观测到。相应的观测过程的现实就是每个字母所对应的声音信号的振幅。为了建立上述对应关系,首先应对
该字的一组观测样本(该字的若干个声音信号)进行学习,也就是说相应的状态序列缺失的情况下进行HMM的参数估计。用数理统计的语言来说,就是不完全数据的参数估计。

在学习了每个字的参数后,就可以用来识别。也就是对一组观测样本(一个字的声音信号),找到最大可能产生该观测样本的那个模型作为该字的代表。

基于方法的孤立词识别系统的基本思想:在训练阶段,用HMM的训练算法(如Baum一welch算法),建立系统词汇表中每个单词对应的HMM,记为λ。在识别阶段,使用前向一后向算法或Viterbi算法求出各个概率P(O/λi),其中O为待识别词的观测值序列。最后选取最大的P(O/λi))所对应的单词斌为O的识别结果。

下面是我的python实现

<span style="font-size:14px;"><span style="font-size:10px;">import numpy as np

#########author Marsharll#########
#############2015.11.9############
##########version 1.0#############



class DHMM:
    #HMM状态为N个，可以观察到的观测数为M
    #transition_matrix是N*N的状态转移矩阵
    #transition_matrix每一行的概率和为1，不一定对称
    #observation_matrix是观测概率矩阵，含义是某时刻t处于状态s的条件下生成观测v的概率，大小为N*M
    #observation_matrix每一行和为1
    #PI是初始状态概率分布，为N*1元组，含义是初始时刻处于某一状态的概率
    def __init__(self,state_num,visible_num):
        self.__state_num=state_transition_matrix.shape[0]
        self.__visible_num=observation_prob_matrix.shape[1]
        self.__transition_matrix=np.zeros(state_num,state_num)
        self.__PI=np.zeros(state_num)
        self.__observation_matrix=np.zeros(state_num,visible_num)
        self.__epsilon=0.0001
    
    
    #observations为t=1到t=T时刻的观测序列，T*1元组
    def __forward(self,observations):
        #实质是动态规划算法
        #__alpha为N*1矩阵
        T=len(observations)
        temp=self.__observation_matrix
        temp=temp.T
        temp1=temp[observations[0]]
        self.__alpha=np.zeros(T,self.__state_num)
        #1. 初始化：根据observations[0]计算t=0时刻所有状态的局部概率alpha
        self.__alpha[0]=np.multiply(self.__PI,temp1)#矩阵对应元素相乘
        #2. 归纳：递归计算每个时间点，t=1，… ，T-1时的局部概率
        for t in xrange(T-1):#0,1,...T-2
            for i in xrange(self.__state_num):
                sumnum=0
                for j in xrange(self.__state_num):
                    #弄清楚计算表达式的实际含义，i、j不要弄混淆，时间下标t不要弄错
                    sumnum+=self.__alpha[t][j]*self.__transition_matrix[j][i]
                self.__alpha[t+1][i]=sumnum*self.__observation_matrix[i][observations[t+1]]
        #3. 终止：观察序列的概率等于T时刻所有局部概率之和
        return np.sum(self.__alpha[T-1])        
    
    #后向算法
    def __backward(self,observations):
        T=len(observations)
        #1. 初始化
        self.__beta=np.zeros(T,self.__state_num)
        for i in xrange(self.__state_num):
            self.__beta[T-1][i]=1.0
        
        
        #2. 递归
        for t in xrange(T-2,-1,-1):#T-2,T-3,...,1,0
            for i in xrange(self.__state_num):
                sumnum=0
                for j in xrange(self.__state_num):
                    sumnum+=self.__transition_matrix[i][j]*self.__observation_matrix[j][observations[t+1]]*self.__beta[t+1][j]
                self.__beta[t][i]=sumnum
        
        
        #3. 终止        
        
        return np.sum(self.__beta[0])
            
    
    
    def viterbi(self,observations):
        #维特比算法，根据观测序列求取隐藏状态的最佳路径
        #实质是动态规划算法
        
        #delta[t][i]是t时刻到达状态i的所有序列概率中最大的概率
        delta=np.zeros(len(observations),self.__state_num)
        #psi记录路径
        psi=np.zeros(len(observations),self.__state_num)
        
        #初始化0时刻
        for x in xrange(self.__state_num):
            delta[0][x]=self.__PI[x]*self.__observation_matrix[x][observations[0]]
            psi[x][i]=0
        
        for t in xrange(1,len(observations)):#1,2,...,T-1
            for j in xrange(self.__state_num):
                for i in xrange(self.__state_num):
                    #弄清楚计算式的物理含义，i、j不要弄混淆，时间下标t不要用错
                    if (delta[t][j] < delta[t-1][i]*self.__transition_matrix[i][j]):
                        delta[t][j] = delta[t-1][i]*self.__transition_matrix[i][j]
                        #记录到达最大概率时的路径
                        psi[t][j] = i
                    delta[t][j] *= self.__observation_matrix[j][observations[t]]
                    
        # termination: find the maximum probability for the entire sequence (=highest prob path)
        p_max = 0 # max value in time T (max)
        path = np.zeros((len(observations)))
        for i in xrange(self.__state_num):
            if (p_max < delta[len(observations)-1][i]):
                p_max = delta[len(observations)-1][i]
                path[len(observations)-1] = i
                    
        # path backtracing
        # path = numpy.zeros((len(observations)),dtype=self.precision)
        for i in xrange(1, len(observations)):
            path[len(observations)-i-1] = psi[len(observations)-i][ path[len(observations)-i] ]
        return path                    
        
    
    def __calXi(self,observations):
        #计算估计参数时用到的中间变量
        xi=np.zeros((len(observations),self.__state_num,self.__state_num))
        for t in xrange(len(observations)-1):#0,1,...,T-2
            sumnum=0.0
            for i in xrange(self.__state_num):
                for j in xrange(self.__state_num):
                    temp=self.__alpha[t][i]*self.__transition_matrix[i][j]*self.__observation_matrix[j][observations[t+1]]*self.__beta[t+1][j]
                    sumnum+=temp
            for i in xrange(self.__state_num):
                for j in xrange(self.__state_num): 
                    numerator=self.__alpha[t][j]*self.__transition_matrix[i][j]*self.__observation_matrix[j][observations[t+1]]*self.__beta[t+1][j]
                    xi[t][i][j]=numerator/sumnum
        return xi
    
    def __calgamma(self,xi,T):
        #计算估计参数时用到的中间变量
        self._gamma=np.zeros(T,self.__state_num)
        for t in xrange(T):
            for i in xrange(self.__state_num):
                self._gamma[t][i]=sum(xi[t][i])
        return self._gamma
    
    
    
    def __update_trans_matrix(self,observations,xi):
        #更新状态转移矩阵
        new_trans_matrix=np.zeros(self.__state_num,self.__state_num)
        for i in xrange(self.__state_num):
            for j in xrange(self.__state_num):
                numerator=0.0
                denominator=0.0
                for t in xrange(len(observations)-1):
                    numerator+=self.__eta1(t,len(observations)-1)*xi[t][i][j]
                    denominator+=self.__eta1(t,len(observations)-1)*gamma[t][i]
                new_trans_matrix[i][j]=numerator/denominator
        return new_trans_matrix
    
    
    
    def __eta1(self,t,T):
        #控制时间序列中不同观测值的权重，这里都是1
        return 1.    
    
    
    def __update_B(self,observations):
        #更新状态-观测矩阵
        B_new=np.zeros(self.__state_num,self.__visible_num)
        for j in xrange(self.__state_num):
            for k in xrange(self.__visible_num):
                numer = 0.0
                denom = 0.0
                for t in xrange(len(observations)):
                    if observations[t] == k:
                        numer += self._gamma[t][j]
                    denom += self._gamma[t][j]
                B_new[j][k] = numer/denom
        
        return B_new        
        
       
    
    def __BaumWelch(self,observations):
        #非监督学习，估计参数，是训练过程
        # E-step
        self.__forward(observations)
        self.__backward(observations)
        xi=self.__calXi(observations)
        self.__calgamma(xi,len(observations))
        
        # M-step
        newmodel={}
        #newmodel["PI"]=self._gamma[0]/double(np.sum(self._gamma[0]))
        newmodel["PI"]=self._gamma[0]
        newmodel["transition_matrix"]=self.__update_trans_matrix(observations,xi)
        newmodel["B"]=self.__update_B(observations)
        return newmodel
    
   
    def __update_model(self,newmodel):
        #更新模型参数
        self.__PI=newmodel["PI"]
        self.__transition_matrix=newmodel["transition_matrix"]
        self.__observation_matrix=newmodel["B"]
    
    
    
    def train(self,observations,max_iteration=100):
        #训练HMM模型参数
        prob_old=100000
        iteration=0
        while(iteration<max_iteration):
            iteration+=1
            newmodel=self.__BaumWelch(observations)
            prob_new=np.log(sum(self.__alpha[-1]))
            self.__update_model(newmodel)
            #如果概率收敛，结束训练
            if(abs(prob_new-prob_old) < self.__epsilon):
                break
            #否则继续训练
            prob_old=prob_new</span></span>

HMM的介绍及实现 - 阿凡卢 - 博客园

几个常用机器学习算法 - 隐马尔可夫模型

如何用简单易懂的例子解释隐马尔可夫模型

浅析隐马尔可夫模型

https://github.com/tostq/Easy_HMM