By

文章目录
  1. 1. 隐马尔可夫定义
  2. 2. 代码功能
  3. 3. 代码
  4. 4. 小记

隐马尔可夫定义

隐马尔可夫模型是关于时序的概率模型，描述由一个隐藏的隐马尔可夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测而产生观测随机序列的过程。

代码功能

下面的代码用python实现了隐马尔科夫模型的概率计算和预测部分，主要是前向后向算法和维特比算法。输入包括状态转移矩阵A和观测矩阵B,初始状态概率向量pi以及观测序列，实现用前向算法和后向算法分别打印前向概率矩阵和后向概率矩阵，打印观测序列的概率，实现用维特比算法打印在观测序列下，概率最大的状态序列，并打印概率和状态序列，输入变量在类的init函数中指定

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
from numpy import *
class HMM:
    def __init__(self):
        self.A=array([(0.5,0.2,0.3),(0.3,0.5,0.2),(0.2,0.3,0.5)])
        self.B=array([(0.5,0.5),(0.4,0.6),(0.7,0.3)])
        self.pi=array([(0.2),(0.4),(0.4)])
        self.o=[0,1,0]
        self.t=len(self.o)#观测序列长度
        self.m=len(self.A)#状态集合个数
        self.n=len(self.B[0])#观测集合个数
    def qianxiang(self):
        #t时刻部分观测序列为o1,o2,ot,状态为qi的概率用矩阵x表示，
        #则矩阵大小行数为观测序列数，列数为状态个数
        self.x=array(zeros((self.t,self.m)))
        #先计算出时刻1时，观测为o1,状态为qi的概率
        for i in range(self.m):
            self.x[0][i]=self.pi[i]*self.B[i][self.o[0]]
        for j in range(1,self.t):
            for i in range(self.m):
                #前一时刻所有状态的概率乘以转移概率得到i状态概率
                #i状态的概率*i状态到j观测的概率
                temp=0
                for k in range(self.m):
                    temp=temp+self.x[j-1][k]*self.A[k][i]
                self.x[j][i]=temp*self.B[i][self.o[j]]
        result=0
        for i in range(self.m):
            result=result+self.x[self.t-1][i]
        print u"前向概率矩阵及当前观测序列概率如下："
        print self.x
        print result
    def houxiang(self):
        #t时刻状态为qi,从t+1到T观测为ot+1,ot+2,oT的概率用矩阵y表示
        #则矩阵大小行数为观测序列数，列数为状态个数
        self.y=array(zeros((self.t,self.m)))
        #下面为对最终时刻的所有状态，接下来的观测序列概率初始化为1
        #(可以理解为接下来没有观测所有为1)
        for i in range(self.m):
            self.y[self.t-1][i]=1
        j=self.t-2
        #j时刻为i状态，转移到k状态，k状态观测为oj+1,
        #再乘以j+1时刻状态为k的B矩阵的值，对k遍历相加，
        #即为j时刻i状态前提下，后面满足观测序列的概率
        while(j>=0):
            for i in range(self.m):
                for k in range(self.m):
                    self.y[j][i]+=self.A[i][k]*self.B[k][self.o[j+1]]*self.y[j+1][k]
            j=j-1
        #第一个状态任意，观测为o1,所以对所有第一个状态概率相加
        result=0
        for i in range(self.m):
            result=result+self.pi[i]*self.B[i][self.o[0]]*self.y[0][i]
        print u'后向概率矩阵及当前观测序列概率如下：'
        print self.y
        print result
    def get_stateprobability(self,t,p):
        #打印在观测为self.o的前提下，t时刻，处于状态p的概率,
        #self.x[t][p]表示到t时刻观测为o1,o2,ot,状态为p的概率
        #self.y[t][p]表示在t时刻状态为p的前提下，接下来观测为ot+1,ot+2,oT的概率
        #self.x[t][p]*self.y[t][p]即表示观测为self.o，且t时刻处于状态p的概率,
        #利用贝叶斯公式，除以观测为self.o的概率即为所求
        if(t>self.t or p>self.m):
            print u'输入数据超过范围'
            return
        print u'在时刻'+str(t)+u'处于状态'+str(p)+u'的概率是：'
        temp=self.x[t-1][p-1]*self.y[t-1][p-1]
        total=0
        for i in range(self.m):
            total=total+self.x[t-1][i]*self.y[t-1][i]
        print temp/total
    def viterbi(self):
        #利用模型和观测序列找出最优的状态序列
        #时刻t时，很多路径可以到达状态i,且观测为self.o,
        #每个路径都有自己的概率，最大的概率用矩阵z记录,前一个状态用d矩阵记录
        self.z=array(zeros((self.t,self.m)))
        self.d=array(zeros((self.t,self.m)))
        for i in range(self.m):
            self.z[0][i]=self.pi[i]*self.B[i][self.o[0]]
            self.d[0][i]=0
        for j in range(1,self.t):
            for i in range(self.m):
                maxnum=self.z[j-1][0]*self.A[0][i]
                node=1
                for k in range(1,self.m):
                    temp=self.z[j-1][k]*self.A[k][i]
                    if(maxnum<temp):
                        maxnum=temp
                        node=k+1
                self.z[j][i]=maxnum*self.B[i][self.o[j]]
                self.d[j][i]=node
        #找到T时刻概率最大的路径
        max_probability=self.z[self.t-1][0]
        last_node=[1]
        temp=0
        for i in range(1,self.m):
            if(max_probability<self.z[self.t-1][i]):
                max_probability=self.z[self.t-1][i]
                last_node[0]=i+1
                temp=i
        i=self.t-1
        #self.d[t][p],t时刻状态为p的时候，t-1时刻的状态
        while(i>=1):
            last_node.append(self.d[i][temp])
            i=i-1
        temp=['o']
        temp[0]=int(last_node[len(last_node)-1])
        j=len(last_node)-2
        while j>=0:
            temp.append(int(last_node[j]))
            j=j-1
        print u'路径节点分别为'
        print temp
        print u'该路径概率为'+str(max_probability)
if __name__ == '__main__':
    test=HMM()
    test.qianxiang()
    test.houxiang()
    test.get_stateprobability(3,3)
    test.viterbi()

小记

注释比较多，代码也不是很规范，但是简单易懂，说明的功能都有实现，参考资料及测试数据来源于李航老师的统计学习方法一书。

文章目录
  1. 1. 隐马尔可夫定义
  2. 2. 代码功能
  3. 3. 代码
  4. 4. 小记