隐马尔模型----五：python实现

实现一

本节主要介绍隐马尔的实现，在以下部分我们实现四部分内容并利用课本上的例子进行了验证：

1、前向计算法；
2、后巷计算法
3、在给定模型且观测为self.o的前提下，t时刻，处于状态p的概率
4、维特比预测算法

# -*- encoding:utf-8 -*-

import numpy as np

class HMM:
    def __init__(self):
        self.A=np.array([(0.5,0.2,0.3),(0.3,0.5,0.2),(0.2,0.3,0.5)])
        self.B=np.array([(0.5,0.5),(0.4,0.6),(0.7,0.3)])
        self.pi=np.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 foreward(self):
        '''
        t时刻部分观测序列为o1,o2,...,ot，状态为qi的概率用矩阵X表示，
        矩阵大小为：行数为观测序列数，列为状态个数
        :return:
        '''
        self.X=np.array(np.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 '前向概率矩阵如下：'
        print self.X
        print '当前观测序列概率如下：'
        print result

    def backward(self):
        #t时刻状态为qi，从t+1到T观测为ot+1,ot+2,..,oT的概率矩阵用y来表示
        #则矩阵大小为：行数为观测序列，列数为状态个数
        self.y=np.array(np.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-=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 '后向概率矩阵如下：'
        print self.y
        print '当前观测序列概率如下：'
        print result


    def get_stateprobability(self,t,p):
        '''

        :param t: t时刻
        :param p: 状态p
        :return:
        打印在给定模型且观测为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 '输入数据超过范围'
            return
        print '在时刻'+str(t)+'处于状态'+str(p)+'的概率是：'
        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矩阵记录
        :return:
        '''
        self.z=np.array(np.zeros((self.t,self.m)))
        self.d=np.array(np.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-=1
        print '路径节点分别为：'
        print temp
        print '该路径的概率为：'+str(max_probability)


if __name__ == '__main__':
    test=HMM()
    test.foreward()
    test.backward()
    test.get_stateprobability(3,3)
    test.viterbi()

实验结果为(和课本中的结果是一致的)：

前向概率矩阵如下：
[[ 0.1       0.16      0.28    ]
 [ 0.077     0.1104    0.0606  ]
 [ 0.04187   0.035512  0.052836]]
当前观测序列概率如下：
0.130218
后向概率矩阵如下：
[[ 0.2451  0.2622  0.2277]
 [ 0.54    0.49    0.57  ]
 [ 1.      1.      1.    ]]
当前观测序列概率如下：
0.130218
在时刻3处于状态3的概率是：
0.405750357093
路径节点分别为：
[3, 3, 3]
该路径的概率为：0.0147

实现二

在本部分我们主要接受利用开源包来调用hmm，我们依赖的包为：hmmlearn.
pip install hmmlearn该命令可以安装我们需要的包。

hmmlearn概述

　　hmmlearn实现了三种HMM模型类：GaussianHMM、GMMHMM以及MultinomialHMM。按照观测状态是连续状态还是离散状态，可以分为两类：GaussianHMM和GMMHMM是连续观测状态的HMM模型，而MultinomialHMM是离散观测状态的模型。
　　对于MultinomialHMM的模型，使用比较简单，”startprob_”参数对应我们的隐藏状态初始分布Π, “transmat_”对应我们的状态转移矩阵A, “emissionprob_”对应我们的观测状态概率矩阵B。
　　对于连续观测状态的HMM模型，GaussianHMM类假设观测状态符合高斯分布，而GMMHMM类则假设观测状态符合混合高斯分布。一般情况下我们使用GaussianHMM即高斯分布的观测状态即可。以下对于连续观测状态的HMM模型，我们只讨论GaussianHMM类。
　　在GaussianHMM类中，”startprob_”参数对应我们的隐藏状态初始分布Π, “transmat_”对应我们的状态转移矩阵A, 比较特殊的是观测状态概率的表示方法，此时由于观测状态是连续值，我们无法像MultinomialHMM一样直接给出矩阵B。而是采用给出各个隐藏状态对应的观测状态高斯分布的概率密度函数的参数。
　　如果观测序列是一维的，则观测状态的概率密度函数是一维的普通高斯分布。如果观测序列是N维的，则隐藏状态对应的观测状态的概率密度函数是N维高斯分布。高斯分布的概率密度函数参数可以用μ表示高斯分布的期望向量，Σ表示高斯分布的协方差矩阵。在GaussianHMM类中，“means”用来表示各个隐藏状态对应的高斯分布期望向量μ形成的矩阵，而“covars”用来表示各个隐藏状态对应的高斯分布协方差矩阵Σ形成的三维张量。

使用

现在我们调用hmmlearn来跑一下课本中的例子（维特比预测方法）

# -*- encoding:utf-8 -*-

#参考链接：http://www.cnblogs.com/pinard/p/7001397.html

import numpy as np
from hmmlearn import hmm
import math

states=['box1','box2','box3']
n_states=len(states)
observations=['red','white']
n_observations=len(observations)

start_probability=np.array([0.2,0.4,0.4])
transition_probability=np.array([[0.5, 0.2, 0.3],
  [0.3, 0.5, 0.2],
  [0.2, 0.3, 0.5]])

emission_probability=np.array([[0.5, 0.5],
  [0.4, 0.6],
  [0.7, 0.3]])

#构建模型
model=hmm.MultinomialHMM(n_components=n_states)
model.startprob_=start_probability
model.transmat_=transition_probability
model.emissionprob_=emission_probability

#预测
seen=np.array([[0,1,0]]).T
logprob,box=model.decode(seen,algorithm='viterbi')

pick_ball=[]
for i in range(len(seen)):
    pick_ball.append(observations[seen[i][0]])
hidden_box=[]
for i in box:
    hidden_box.append(states[i])
print "The ball picked:",pick_ball
print "The hidden box", ", ", hidden_box
#要注意的是score函数返回的是以自然对数为底的对数概率值，
maxprob= math.exp(logprob)
print '该路径的概率为：',maxprob

#概率计算：
logprob=model.score(seen)
prob=math.exp(logprob)
print '当前序列的概率为：',prob

实验结果为（是一致的）：

The ball picked: ['red', 'white', 'red']
The hidden box ,  ['box3', 'box3', 'box3']
该路径的概率为： 0.0147
当前序列的概率为： 0.130218

《完》

                所谓的不平凡就是平凡的N次幂。
                                        ----By Ada