HMM代码实现

按照网上的代码，自己敲了一下，改了一点点，理解加深了一下。

还有训练HMM的EM算法没看懂，下次接着看；

参考连接：http://www.cnblogs.com/hanahimi/p/4011765.html

# -*- coding: utf-8 -*-

'''
function: 根据网上的代码，学习实现 HMM，前向计算概率，后向预测状态序列，学习算法参数
date: 2017.8.9
'''

import numpy as np

class HMM(object):
    """docstring Ann, Bnm, piln HMM"""
    def __init__(self, Ann, Bnm, piln):
        self.A = Ann
        self.B = Bnm
        self.pi = piln
        self.N = self.A.shape[0]  #状态的种类个数
        self.M = self.B.shape[1]  #观测序列的长度

    #打印hmm的信息
    def printHmm(self):
        print("=======================================")
        print('hmm content N = ', self.N, ' M = ', self.M)
        for i in range(self.N):
            if i == 0:
                print('hmm.A ', self.A[i,:],'hmm.B', self.B[i,:])
            else:
                print('      ', self.A[i,:], '    ', self.B[i,:])
        print('hmm.pi', self.pi)
        print('========================================')

    '''
    function: 维特比算法
    input: A,B,pi,O
    output: P(o|lambda) 最大的时候，状态的路径序列
    '''
    def viterbi(self, O):
        T = len(O) #观察序列的长度
        
        #初始化,这里从0~T-1
        #delta t行 n列，代表有t个时间点，每个时间点可能有n种状态
        delta = np.zeros((T, self.N), np.float) #二维数组记录计算的所有概率,包括了最有的点
        phi = np.zeros((T, self.N), np.int) #记录概率最大路径的前一个状态
        I = np.zeros(T, np.int) #这里如果不显示表明类型为 np.int,就是float?
        for i in range(self.N):
            delta[0,i] = self.pi[i] * self.B[i,O[0]] #t = 0时刻，各个状态的起始概率
            phi[0,i] = 0 #t=0时刻前缀状态都是0

        #递推
        for t in range(1,T): #从 1 开始
            for i in range(self.N):
                delta[t,i] = self.B[i,O[t]] * np.array([delta[t-1,j]*self.A[j,i] 
                    for j in range(self.N)]).max()
                phi[t,i] = np.array([delta[t-1, j]*self.A[j,i] for j in range(self.N)]).argmax()
        #结束
        prob = delta[T-1,:].max() #T-1时刻是最后时刻，哪个状态在最后时刻概率最大就是最优路径的起始点
        I[T-1] = phi[T-1,:].argmax() #最优路径的起点状态编号
        #状态序列获取
        for t in range(T-2, -1, -1): #从 T-1 到 -1(不包括-1)，间隔是-1，即递减
            I[t] = phi[t+1, I[t+1]]
        return I, prob

    '''
    function: 前向算法计算擦观察序列 O 出现的概率
    input: A,B,pi,O
    output: prob
    '''
    def forward(self, O):
        T = len(O)
        alpha = np.zeros((T, self.N), np.float) #暂存计算的所有概率，按照时间点向前推进
        #初始化
        for i in range(self.N):
            alpha[0,i] = self.pi[i] * self.B[i, O[0]]
        
        #迭代计算
        for t in range(T-1):
            for i in range(self.N): #这里B[i,O[t]]也可以放在for的外面乘
                alpha[t+1,i] = np.array([alpha[t, j]*self.A[j,i]*self.B[i,O[t+1]] for 
                    j in range(self.N)]).sum()

        #终止
        prob = np.array([alpha[T-1, j] for j in range(self.N)]).sum()
        return prob

    '''
    function: 后向算法，计算观测序列出现的概率

    '''
    def backword(self, O):
        T = len(O)
        beta = np.zeros((T, self.N), np.float) #暂存计算的概率

        #初始化
        for i in range(self.N):
            beta[T-1, i] = 1  #从后向前
        #迭代计算
        for t in range(T-2, -1, -1):
            for i in range(self.N):
                beta[t,i] = np.array([[A[j,i] * B[j,O[t+1]] * beta[t+1,j]] for 
                    j in range(self.N)]).sum()

        prob = np.array([self.pi[j] * self.B[i,O[1]] * beta[1,j] for j in range(self.N)]).sum()

        return prob



if __name__ == 'main':

    print('python my HMM')

    #HMM模型的参数
    A = [[0.8125,0.1875],[0.2,0.8]]
    B = [[0.875,0.125], [0.25,0.75]] #每一行的和是 1
    pi = [0.5,0.5]
    hmm = HMM(A,B,pi) #构建HMM

    print(hmm)

print('python my HMM')

#HMM模型的参数
A = np.mat([[0.8125,0.1875],[0.2,0.8]])
B = np.mat([[0.875,0.125], [0.25,0.75]]) #每一行的和是 1
pi = [0.5,0.5]
O = [[1,0,0,1,1,0,0,0,0],
    [1,1,0,1,0,0,1,1,0],
    [0,0,1,1,0,0,1,1,1]]

hmm = HMM(A,B,pi) #构建HMM

#计算前向概率，产生特定观测序列O的概率
prob = hmm.forward(O[0])
print('前向算法产生 O 序列的概率是: ' + str(prob))

#后向算法计算观测序列的概率
prob = hmm.backword(O[0])
print('后向算法概率是： ' + str(prob))
#计算隐含概率，维特比算法
path, prob2 = hmm.viterbi(O[0])
print('产生 O 序列最大概率路径是： ' + str(path))
print('概率是： ' + str(prob2))

hmm.printHmm()