隐马尔科夫模型的简单实现

模型 HMM=(A, B, \pi)

初始化三个参数：隐状态初始状态概率向量\pi，状态转移概率矩阵A，隐状态生成观测状态概率矩阵B。

实现三个方法：

1、2. 计算给定观测状态序列向量的概率（前向和后向两种计算方法），

3. 给定观测状态序列，求出与该序列最匹配的隐状态序列及其概率（Viterbi算法）。

程序如下：

# -*- coding: utf-8 -*-
# @Author : ZhaoKe
# @Time : 2022-09-18 14:43
from typing import List

import numpy as np

class HMM():
    def __init__(self):
        # 对应三个隐状态初始化概率矩阵
        self.init_prob = np.array([0.4, 0.3, 0.3])

        # 隐状态会生成观测状态的概率矩阵
        self.gene_prob = np.array([
            [0.8, 0.2],
            [0.6, 0.4],
            [0.3, 0.7]
        ])

        # 状态转移概率矩阵
        self.tran_prob = np.array([
            [0.4, 0.5, 0.1],
            [0.3, 0.4, 0.3],
            [0.1, 0.5, 0.4]
        ])

    # 求解生成特定的序列的概率
    # 例如，输入[0, 1, 0]
    def forward_prob_generate(self, series: List):
        # 行数等于隐状态个数，列数等于序列长度
        # print(len(self.init_prob), len(series))
        res = np.zeros((len(self.init_prob), len(series)))
        res[:, 0] = np.multiply(self.init_prob, self.gene_prob[:, series[0]])
        for i in range(1, len(series)):
            res[:, i] = np.multiply(np.matmul(res[:, i-1], self.tran_prob), self.gene_prob[:, series[i]])
        print(res)

    # 后向生成特定序列的概率
    def backward_prob_generate(self, series: List):
        res = np.zeros((len(self.init_prob), len(series)))
        n = len(series)
        res[:, n-1] = 1
        for i in range(n-2, -1, -1):
            res[:, i] = np.matmul(self.tran_prob, np.multiply(self.gene_prob[:, series[i+1]], res[:, i+1]))
        print(res)
        print(np.sum(self.init_prob * self.gene_prob[:, 1] * res[:, 0]))

    # 维特比算法求解特定观测序列的最优隐状态向量及其概率
    def viterbi(self, series: List):
        # 动态规划
        dp = np.zeros((3, 4))
        dp[:, 0] = 1
        dp[:, 1] = self.init_prob * self.gene_prob[:, series[0]]
        per = list(map(int, np.zeros(len(series)+1)))
        per[1] = np.argmax(dp[:, 1])
        # print(dp)
        # print(per)
        for i in range(2, len(series)+1):
            for j in range(0, 3):
                tmp = dp[:, i-1] * self.tran_prob[:, j]
                dp[j, i] = np.max(tmp) * self.gene_prob[j][series[i-1]]
            print(dp[:, i])
            per[i] = np.argmax(dp[:, i])
        print("动态规划矩阵：\n", dp)
        print("最优路径：", per)
        print("概率计算：", dp[per[-1], 3])


if __name__ == '__main__':
    hmm = HMM()
    # hmm.forward_prob_generate([0, 1, 0])
    # hmm.backward_prob_generate([1, 0, 1])
    # hmm.viterbi([0,0,0])
    # hmm.viterbi([0,0,1])
    # hmm.viterbi([0,1,0])
    # hmm.viterbi([0,1,1])
    hmm.viterbi([1,0,0])
    # hmm.viterbi([1,0,1])
    # hmm.viterbi([1,1,0])
    # hmm.viterbi([1,1,1])

和例题对照，结果正确