个人学习笔记：HMM算法代码

本篇文章主要是个人在jupyter上实现的HMM算法中维特比算法实现的代码记录。

import numpy as np
import pandas as pd

class HiddenMM:
    """
    创建隐马尔可夫模型，定义三个主要方法分别对应解决三个问题：
    （1）概率计算问题：即已知模型参数 $$\lambda$$ =（A,B,pi）,
                        获得观测序列O后，估计该观测出现的概率P(O|$$\lambda$$)；
    （2）学习问题：即已知观测序列O=(O_1,O_2,O_3,...,O_T)，
                    对模型参数$$\lambda$$进行估计，
                    使得P(O|$$\lambda$$)最大，采用EM算法进行估计；
    （3）预测问题（decoding）：即已知模型参数$$\lambda$$，
                                对于给定的观测序列O=(O_1,O_2,O_3,...,O_T)，
                                求其最可能的隐藏状态序列P(I|O)；
    
    """
    def __init__(self,initial_prob = None,
                 transition_prob = None,
                 observation_prob = None):
        """
        创建离散变量的隐马尔可夫模型
        
        初始化参数
        ---------------
        initial_prob:<numpy.ndarray>,(n_states,) or default None
                    即一维矩阵，指示初始状态概率,即模型参数pi；
        transition_prob:<numpy.ndarray>,(n_states,n_states) or default None
                    即二维矩阵，指示状态转移概率矩阵，模型参数A；
                    其中元素(i,j)表示状态i转移到状态j的概率；
        observation_prob:<numpy.ndarray>,(n_states,n_observation) or default None
                    即二维矩阵，指示的是每个状态下取不同观测值的概率，模型参数B；
                    其中n_observation指示的是观测值唯一取值个数；
                    而在连续变量下，该矩阵将由不同分布进行刻画； 
                    在离散状态下，(i,j)表示状态i下观察变量取值为j的概率；
        属性
        ---------------
        initial_prob
        transition_prob
        observation_prob
        n_states
        n_observation
        """
        ###校验初始概率矩阵
        if not initial_prob is None:
            if initial_prob.ndim ==1:
                self.initial_prob = initial_prob
            elif (initial_prob.ndim == 2) and (initial_prob.shape[1] == 1):
                self.initial_prob = initial_prob.ravel()
            else:
                print( initial_prob.shape[1] == 1)
                raise ValueError("状态初始概率分布格式错误，应该为(n_states,)")
            self.n_states = self.initial_prob.size
            initial_judge = True
        else:
            self.initial_prob = None
            self.n_states = None
            initial_judge = False
        
        ###校验状态转移矩阵
        if not transition_prob is None:
            if transition_prob.shape[0] == transition_prob.shape[1]:
                self.transition_prob = transition_prob
                transition_judge = True
            else:
                raise ValueError("转移概率矩阵为非方阵！")
        else:
            self.transition_prob = None
            transition_judge = False
            
        ###校验观测概率分布矩阵
        if not observation_prob is None:
            if observation_prob.ndim == 2:
                self.observation_prob = observation_prob
                observation_judge = True
            else:
                raise ValueError("观测值概率分布矩阵非2维矩阵!")
        else:
            self.observation_prob = None
            observation_judge = False
        
        ###一致性校验
        if initial_judge & transition_judge & observation_judge:
            if self.n_states != self.transition_prob.shape[0]:
                raise ValueError("概率转移矩阵与初始化概率分布维度不一致！")
            if self.n_states != self.observation_prob.shape[0]:
                raise ValueError("观测值概率矩阵与初始化概率分布维度不一致！")                
                
    def fit_hmm(self,obseration_seq,
                initial_A = None,
                initial_B = None,
                initial_pi = None,
                max_iter = 1000,
                tol = 0.0001):
        """
        基于EM算法对给定的观测序列进行无监督拟合学习，估计模型参数，
        即Baum_Welch算法；
        对于有人工标注的数据，则直接使用样本数据进行参数计算即可；
        
        参数
        -------------
        observation_seq:<numpy.ndarray>,(n_sample,len_seq) or default None
                        为观测序列样本，其中n_sample为样本的数量，
                        len_seq指示序列长度，即参数T；
        initial_pi:<numpy.ndarray>,(n_states,) or default None
                EM算法的迭代初始化参数，初始状态分布概率；
        initial_A:<numpy.ndarray>,(n_states,n_states) or default None
                EM算法的迭代初始化参数，状态转移转移概率矩阵；
        initial_B:<numpy.ndarray>,(n_states,n_observation) or default None
                EM算法的迭代初始化参数，观测值概率分布矩阵；
        max_iter:int,default 100
                EM算法的最大迭代停止次数；
        tol:float,default 0.001
                EM算法的迭代停止条件，当参数精度达到停止条件时，
                或者达到最大迭代次数时，停止迭代。
        """
        pass
    
    
    def _Estep():
        """
        求完全对数似然函数在给定观测与前一轮参数下关于隐含状态条件概率分布的期望；
        """
        pass
    
    def _Mstep():
        """
        期望最大化
        """
        pass
    
    def _forward(self,obs_seq,mode:("loop","vector")="loop"):
        """
        前向算法，即前向递推求解给定序列出现的概率P(obs_seq|model_para);
        """
        len_seq = len(obs_seq)
        n_states = self.n_states
        alpha_forward = np.zeros((n_states,len_seq))
        pi = self.initial_prob
        A = self.transition_prob
        B = self.observation_prob
        
        if mode == "loop":
            ###################使用循环遍历进行计算###########################
            ###参数初始化
            obs_value = obs_seq[0]
            for s in range(n_states):
                alpha_forward[s,0] = pi[s]*B[s,obs_value]

            for obs_index in range(1,len_seq):
                obs_value = obs_seq[obs_index]
                for s in range(n_states):
                    alpha_forward[s,obs_index] = sum(
                    [
                        alpha_forward[pre_s,obs_index-1]*A[pre_s,s]
                        for pre_s in range(n_states)
                    ]
                    )*B[s,obs_value]
            ###################使用循环遍历进行计算—>结束###########################
        elif mode == "vector":
            ####################使用向量、矩阵乘法进行计算##################################
            alpha_forward[:,0] = pi*B[:,0]
            for obs_index in range(1,len_seq):
                obs_value = obs_seq[obs_index]
                alpha_forward[:,obs_index] = np.dot(alpha_forward[:,obs_index-1],A)*B[:,obs_value]
            ####################使用向量、矩阵乘法进行计算—>结束###########################
        else:
            raise ValueError("mode参数取值错误！")
            
        obs_condition_prob_gamma = alpha_forward.sum(axis = 0)[-1]
        return obs_condition_prob_gamma
        
    def _backward(self,obs_seq,mode:("loop","vector")= "loop"):
        """
        后向算法,即与前向算法的视角相反，前向算法对于当前结点的概率计算，是计算所有
        前一阶段的所有可能状态结点转移到当前状态结点的概率和，基于该概率和再结合该状态
        下取值指定观测值的概率，即获得了前向概率；
        而后向算法则向后看，即基于当前状态结点，向下一阶段进行状态转移时，可以基于状态转移
        矩阵计算获得到达下一阶段并取值为特定观测值的所有可能状态的对应概率，这些概率和即该
        阶段该结点的后向概率。
        """
        len_seq = len(obs_seq)
        n_states = self.n_states
        A = self.transition_prob
        B = self.observation_prob
        pi = self.initial_prob
        beta_backward = np.zeros((n_states,len_seq))
        
        if mode == "loop":
            ###############进行循环计算###########
            for s in range(n_states):
                beta_backward[s,len_seq-1] = 1
            for obs_index in range(len_seq-2,-1,-1):
                obs_value_next = obs_seq[obs_index+1]
                for s in range(n_states):
                    beta_backward[s,obs_index] = sum([
                        A[s,next_s]*B[next_s,obs_value_next]*beta_backward[next_s,obs_index+1]
                        for next_s in range(n_states)
                    ])
            ###############循环计算—>结束#########
        elif mode == "vector":
            ################矩阵、向量算法计算####################
            #逻辑相对复杂
            #（1）初始化
            beta_backward[:,len_seq-1] = 1
            #（2）因牵扯到自身的依赖计算，无法一次计算
            for obs_index in range(len_seq-2,-1,-1):
                obs_value_next = obs_seq[obs_index+1]
                temp_trans = (A*B[:,obs_value_next])
                temp_back = beta_backward[:,obs_index+1].reshape(-1,1)
                beta_backward[:,obs_index]= np.dot(temp_trans,temp_back).ravel()
            ################矩阵、向量算法计算####################
        else:
            raise ValueError("mode参数取值错误")
        
        #print(beta_backward)
        obs_0 = obs_seq[0]
        prob = pi*B[:,obs_0]*beta_backward[:,0]
        return prob.sum()
    
    def decode(self,obs_seq,mode:("loop","vector") = "loop"):
        """
        解码算法，即使用维特比算法对给定观测样本进行隐含状态求解。
        
        参数
        -------
        obs_seq:<numpy.ndarray>,(len_seq,)
                一个观测序列样本，其中：
                len_seq为序列的长度；
        返回值
        -------
        hidden_states:list of length 'len_seq'
                    该观测最有可能对应隐含状态序列
        max_prob:float
                该最大概率隐含序列的概率
        """
        eps = np.finfo(float).eps
        if obs_seq.ndim == 1:
            ##转换为(1,len_seq)矩阵数据
            obs_seq = obs_seq.reshape(1,-1)
        len_seq = obs_seq.shape[1]
        ###输入的观测序列的个数，应该为1
        num_seq = obs_seq.shape[0]
        if num_seq != 1:
            raise ValueError("只能处理单一观测序列！")
        
        A = self.transition_prob
        B = self.observation_prob
        pi = self.initial_prob
        
        ###进行算法参数初始化
        n_states = self.n_states
        viterbi_mat = np.zeros((n_states,len_seq))
        path_mat = np.zeros((n_states,len_seq),dtype = int)
        
        ###获取序列的第一个参数
        obs_0 = obs_seq[0,0]

        ###viterbi_mat实际为动态规划表
        if mode == "loop":
            ####################以循环的方式进行计算################################
            for start in range(n_states):
                path_mat[start,0] = 0
                ##初始状态，计算状态start下取值obs_value的概率
                ##为防止出现log(0)错误，引入eps
                viterbi_mat[start,0] = np.log(pi[start]+eps) + np.log(B[start,obs_0]+eps)
            ###对序列中的其余观测进行计算
            for obs_index in range(1,len_seq):
                obs_value = obs_seq[0,obs_index]
                for s in range(n_states):
                    ###在观测obs_value条件下前一阶段所有可能的状态转移到现在状态的概率
                    state_trans_prob = [
                        viterbi_mat[pre_s,obs_index -1]
                        +np.log(A[pre_s,s]+eps)
                        +np.log(B[s,obs_value]+eps)
                        for pre_s in range(n_states)
                    ]
                    ###获得最大概率
                    viterbi_mat[s,obs_index] = np.max(state_trans_prob)
                    ###该最大概率的前置节点索引,哪一个pre_s
                    path_mat[s,obs_index] = np.argmax(state_trans_prob)
                    #print(path_mat[s,obs_index])
            ####################以循环的方式进行计算—>结束##########################
        elif mode == "vector":
            #######################以矩阵运算的方式计算###########################
            #初始化，观测序列转换为1维序列
            obs_seq = obs_seq.ravel()
            viterbi_mat[:,0] = pi*B[:,obs_0]
            path_mat[:,0] = 0
            for obs_index in range(1,len_seq):
                obs_value = obs_seq[obs_index]
                ##向量方式计算基于前一阶段状态的概率转移，结果为(n_states,n_states)
                state_trans_mat = viterbi_mat[:,[obs_index-1]]*A
                ###从每列选出最大的概率值，即前一状态转移到该状态的最大概率及其编号
                max_pre_state = state_trans_mat.argmax(axis = 0)
                max_state_prob = state_trans_mat.max(axis = 0)
                ###计算当前阶段做出决策（选择状态）后的总的路径概率，
                ###遍历每一个状态，做出决策（即B[state,obs_value]）,
                ###计算该决策下的最佳结果（即概率最大）
                ###这样每个状态都对应了一个该阶段选择该状态的最大概率，
                ###并记录使该阶段在该状态/决策（特定状态）下达到最大概率的前一阶段状态
                viterbi_mat[:,obs_index] = max_state_prob*B[:,obs_value]
                path_mat[:,obs_index]= max_pre_state
            ####################以矩阵运算的方式计算—>结束##########################
        else:
            raise ValueError("mode参数输入错误！")
        
        print("动态规划矩阵：\n",viterbi_mat)
        print("回溯路径矩阵：\n",path_mat)
        ###获得最大的似然概率、最大概率
        max_seq_logprob = viterbi_mat[:,len_seq-1].max()
        ###进行路径回溯
        path_ = viterbi_mat[:,len_seq-1].argmax()
        max_seq_logprob_path = [path_]
        #print("动态规划表格：\n",viterbi_mat)
        #print("路径节点记录表格：\n",path_mat)
        for obs_index in range(len_seq-1,0,-1):
            path_ = path_mat[path_,obs_index]
            max_seq_logprob_path.append(path_)
        max_seq_logprob_path = max_seq_logprob_path[::-1]
        if mode == "loop":
           max_seq_logprob = np.power(np.e,max_seq_logprob)
        return max_seq_logprob,max_seq_logprob_path

Q = [0,1,2]
V = [0,1]
A = np.array([[0.5,0.2,0.3],
             [0.3,0.5,0.2],
              [0.2,0.3,0.5]])
B = np.array([[0.5,0.5],
             [0.4,0.6],
              [0.7,0.3]])
pi = np.array([[0.2],[0.4],[0.4]])
model = HiddenMM(initial_prob=pi,
                 transition_prob=A,
                 observation_prob=B)
model.decode(obs_seq=np.array([0,1,0]))

动态规划矩阵：
 [[-2.30258509 -3.57555077 -4.88488409]
 [-1.83258146 -2.9877641  -4.59720202]
 [-1.27296568 -3.17008566 -4.21990779]]
回溯路径矩阵：
 [[0 2 1]
 [0 2 1]
 [0 2 2]]





(0.01470000000000004, [2, 2, 2])

A1 = np.array([[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]])
B1 = np.array([[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]])
PI = np.array([0.2, 0.4, 0.4])
model = HiddenMM(initial_prob=PI,
                 transition_prob=A1,
                 observation_prob=B1)
model.decode(obs_seq=np.array([0,1,0,1]))

动态规划矩阵：
 [[-2.30258509 -3.57555077 -4.88488409 -6.27117845]
 [-1.83258146 -2.9877641  -4.59720202 -5.80117482]
 [-1.27296568 -3.17008566 -4.21990779 -6.11702777]]
回溯路径矩阵：
 [[0 2 1 0]
 [0 2 1 1]
 [0 2 2 2]]





(0.003024000000000012, [2, 1, 1, 1])

model._forward(obs_seq=np.array([0,1,0]))

0.130218

model._backward(obs_seq=np.array([0,1,0]))

0.130218

model._backward(obs_seq=np.array([0,1,0]))

0.130218

A1 = np.array([[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]])

x = np.array([[1],[2],[3]])

y[0,0]=100

y.argmax(axis = 0)

array([0, 1, 2], dtype=int64)

y[range(len(y.argmax(axis = 0))),y.argmax(axis = 0)]

array([100. ,   1. ,   1.5])

np.repeat(x,x.shape[0],axis=1)

array([[1, 1, 1],
       [2, 2, 2],
       [3, 3, 3]])

A1 = np.array([[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]])
B1 = np.array([[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]])
PI = np.array([0.2, 0.4, 0.4])

x  = np.array([1,2,3]).reshape(-1,1)

def test(A,B,pi,obs_seq):
        eps = np.finfo(float).eps
        T = len(obs_seq)
        N = len(pi)
        # initialize the backward trellis
        backward = np.zeros((len(pi), T))

        for s in range(N):
            backward[s, T - 1] = 1

        for t in reversed(range(T - 1)):
            ot1 = obs_seq[t + 1]
            for s in range(N):
                backward[s, t] = sum(
                    [
                        (A[s, s_] + eps)
                        * (B[s_, ot1] + eps)
                        * backward[s_, t + 1]
                        for s_ in range(N)
                    ]  # noqa: C812
                )
        return backward

test(A,B,pi,obs_seq=np.array([0,1,0]))

array([[0.2451, 0.54  , 1.    ],
       [0.2622, 0.49  , 1.    ],
       [0.2277, 0.57  , 1.    ]])

x = np.array([1,2,3,4,5,6])

x[range(4,-1,-1)]

array([5, 4, 3, 2, 1])