隐形马尔可夫链代码——python源代码，看完你就完全懂了

隐形马尔可夫链代码

隐形马尔可夫链在李航的统计学习方法里面已经介绍的很详细了(李航统计学习方法P193-P196)，我这里就不献丑了，这篇文章主要是分享隐形马尔可夫链的代码。

隐形马尔可夫有三个问题：
1、概率计算问题：
给定模型$\lambda =(A,B,\pi )$和观测序列$O=(o_1,o_2,...,o_T)$，计算在模型$\lambda$下观测序列$O$出现的概率$P(O|\lambda )$。

2、学习问题：
已知观测序列$O=(o_1,o_2,...,o_T)$,估计模型$\lambda =(A,B,\pi )$参数，使得在该模型下观测序列概率$P(O|\lambda )$最大，即用最大似然估计的方法估计参数。

3、预测问题
也成为解码。已知模型$\lambda =(A,B,\pi )$和观测序列$O=(o_1,o_2,...,o_T)$，求对给定观测序列条件概率$P(O|\lambda )$最大的状态序列$I=(i_1,i_2,...i_T)$

一、概率计算问题
该问题的计算方法有前向算法和后向算法，详情可见统计学习方法P198-P203。模型和代码一起看更容易消化呢。

class HMM1:
    def __init__(self):
        self.A=[[0,1,0,0],[0.4,0,0.6,0],[0,0.4,0,0.6],[0,0,0.5,0.5]]       ##转移概率
        self.B=[[0.5,0.5],[0.3,0.7],[0.6,0.4],[0.8,0.2]]                   ##观测概率
        self.W=[0.25,0.25,0.25,0.25]
        ##初始状态
    def cal_q(self,Q):
        """
        前向算法
        :param Q: 观测数据 1，0，1，0，1，0
        :return:
        """
        a=[[0 for i in range(len(self.B))] for j in range(len(Q))]
        for i in range(len(self.B)):
            a[0][i]=self.W[i]*self.B[i][Q[0]]
        for t in range(1,len(Q)):
            for i in range(len(self.B)):
                a[t][i]=0
                for j in range(len(self.B)):
                    a[t][i]+=a[t-1][j]*self.A[j][i]*self.B[i][Q[t]]
        return sum(a[-1])
    def cal_q2(self,Q):
        """
        后向算法
        :param Q: 观测数据 1，0，1，0，1，0
        :return:
        """
        b = [[0 for i in range(len(self.B))] for j in range(len(Q))]
        for i in range(len(self.B)):
            b[-1][i]=1
        for t in range(len(Q)-2,-1,-1):
            for i in range(len(self.B)):
                b[t][i]=0
                for j in range(len(self.B)):
                    b[t][i]+=b[t+1][j]*self.A[i][j]*self.B[j][Q[t+1]]
        q=0
        for i in range(len(self.B)):
            q+=b[0][i]*self.W[i]*self.B[i][Q[0]]
        return q
T=HMM1()
print(T.cal_q([0,1,1,0]))
print(T.cal_q2([0,1,1,0]))

二、学习问题
学习问题的算法有Baum-Welch算法，该算法采用的是EM算法。同样，模型推荐还是看书，李航的那本书将的是最清楚的，当然配合下面的代码会更加清晰。

import numpy as np
class Baum:
    def __init__(self,Q):
        self.l=len(Q)
        self.Q=Q
        self.n=len(set([i for j in Q for i in j]))  ##观测状态的个数
        self.m=4
        self.A=np.zeros((4,4))
        self.B=np.zeros((4,self.n))
        self.W=np.zeros(4)
        self.T=len(Q[0])
        self.initial()
    def initial(self):
        for i in range(4):
            for j in range(4):
                self.A[i,j]=0.25
        for i in range(4):
            for j in range(self.n):
                self.B[i,j]=1lf.n
        for i in range(4):
            self.W[i]=0.25

    def main(self):
        for iter in range(20):
            for l in range(len(self.Q)):
                gama,sigma=self.update(self.Q[l])
                for i in range(4):
                    for j in range(4):
                        self.A[i,j]=sum([sigma[t,i,j] for t in range(self.T-1)])/sum([gama[t,i] for t in range(self.T-1)])
                for i in range(4):
                    for j in range(self.n):
                        temp=0
                        for t in range(self.T):
                            if self.Q[l][t]==j:
                                temp+=gama[t,i]
                        self.B[i,j]=temp/sum([gama[t,i] for t in range(self.T)])
                for i in range(4):
                    self.W[i]=gama[0][i]
                
    def update(self,O):
        a = [[0 for i in range(len(self.B))] for j in range(len(O))]
        for i in range(len(self.B)):
            a[0][i]=self.W[i]*self.B[i][O[0]]

        for t in range(1,len(O)):
            for i in range(len(self.B)):
                a[t][i]=0
                for j in range(len(self.B)):
                    a[t][i]+=a[t-1][j]*self.A[j,i]*self.B[i,O[t]]

        b = [[0 for i in range(len(self.B))] for j in range(len(O))]
        for i in range(len(self.B)):
            b[-1][i]=1
        for t in range(len(O)-2,-1,-1):
            for i in range(len(self.B)):
                b[t][i]=0
                for j in range(len(self.B)):
                    b[t][i]+=b[t+1][j]*self.A[i,j]*self.B[j,O[t+1]]

        """
        gama t i 表示给定模型个观测数据的情况下在时刻t系统状态为i的概率
        """
        gama = np.zeros((self.T, 4))
        for t in range(self.T):
            for i in range(4):
                gama[t,i]=a[t][i]*b[t][i]/sum([a[t][k]*b[t][k] for k in range(4)])

        """
        sgma t i，j 表示给定模型个观测数据的情况下在时刻t系统状态为i,且t+1时刻状态为j的概率
        """
        sigma = np.zeros((self.T, 4, 4))
        for t in range(self.T-1):
            temp=0
            for i in range(4):
                for j in range(4):
                    temp+=a[t][i]*self.A[i,j]*self.B[i,O[t]]*b[t+1][j]
            for i in range(4):
                for j in range(4):
                    sigma[t,i,j]=a[t][i]*self.A[i,j]*self.B[i,O[t]]*b[t+1][j]/temp

        return gama,sigma

T=Baum([[0,1,1,0,2,3,0,1],[0,2,3,0,1,2,3,0],list(map(int,list("12300212")))])
T.main()

三、预测问题
预测问题采用的是维比特算法，该算法采用的是动态规划方法求最优路径。

import numpy as np
class 维比特算法:
    def __init__(self):
        self.A=[[0.5,0.2,0.3],[0.3,0.5,0.2],[0.2,0.3,0.5]]       ##转移概率
        self.B=[[0.5,0.5],[0.4,0.6],[0.7,0.3]]                   ##观测概率
        self.W=[0.2,0.4,0.4]
        self.O=[0,1,0]
        self.m=len(self.A) ## 状态的种数
        self.n=len(self.O) ## 观测的个数
        self.sigma=[[0 for _ in range(self.m)] for _ in range(self.n)]
        self.psi=[[0 for _ in range(self.m)] for _ in range(self.n)]


    def main(self):
        self.sigma[0][:]=[self.W[i]*self.B[i][self.O[0]] for i in range(self.m)]
        for i in range(1,self.n):
            for j in range(self.m):
                tem=[self.sigma[i-1][k]*self.A[k][j] for k in range(self.m)]
                tem_max=max(tem)
                self.sigma[i][j]=tem_max*self.B[j][self.O[i]]
                self.psi[i][j]=np.argmax(tem)

        res=[np.argmax(self.sigma[-1][:])]
        for i in range(self.n-1):
            res.insert(0,np.argmax(self.psi[res[-1]][:]))
        print(res)

T=维比特算法()
T.main()

想了解更多算法，数据分析，数据挖掘等方面的知识，欢迎关注ChallengeHub公众号，添加微信进入微信交流群，或者进入QQ交流群。关注公众号或者进群可以直接获得上面的代码，而且也有很多机器学习的资料可以获取。