my HMM

HMM 模型

本节对隐马尔可夫模型进行了详细的公式推导，并实现了 HMM 模型中涉及到的三个经典问题。

公式推导

HMM 模型涉及到大量公式，本文以《统计学习方法》为基础对其中的公式进行了推导，为方便起见，以手写的形式进行了记录，字丑莫怪。

HMM 代码实现

本部分对上述公式进行了实现，涉及到了 HMM 模型中的前向算法、后向算法、baumwelch算法和viterbi算法，即分别对应第一部分中的概率计算、学习和预测。

import numpy as np


class HMM:
    def __init__(self, pi, A, B):
        self.pi = pi  # 初始隐状态概率矩阵
        self.A = A  # 转移概率矩阵
        self.B = B  # 发射概率矩阵
        self.alpha = None
        self.beta = None

    def forward(self, O):
        T = len(O)  # 观测序列长度
        N = len(self.A)  # 状态数
        alpha = np.zeros([T, N])  # 初始化alpha矩阵

        # alpha_1(i)=p_i*b_i(o_1)
        alpha[0, :] = self.pi * self.B[:, O[0]]
        for t in range(1, T):
            # alpha_t+1(i)=\sum_j{alpha_t(j)*a_ji}*b_i(o_t+1)
            alpha[t, :] = np.sum(alpha[t - 1, :] * self.A.T, axis=1) * self.B[:, O[t]]
        # print(alpha)
        self.alpha = alpha  # shape: N * T

        return np.sum(alpha[-1, :])

    def backward(self, O):
        T = len(O)  # 观测序列长度
        N = len(self.A)  # 状态数
        beta = np.zeros([T, N])  # 初始化beta矩阵

        # beta_T=1
        beta[-1, :] = 1
        for t in range(T - 2, -1, -1):
            # beta_t(i)=\sum_j{beta_t+1(j)*a_ij*b_j(o_t+1)}
            beta[t, :] = np.sum(beta[t + 1, :] * self.A * self.B[:, O[t + 1]], axis=1)
        # print(beta)
        self.beta = beta  # shape: N * T

        # \sum_i{pi_i*b_i(o_1)*beta_1(i}
        return np.sum(self.pi * self.B[:, O[0]] * beta[0, :])

    def baumwelch(self, O, criterion=0.001):
        T = len(O)
        N = len(self.A)
        while True:
            self.forward(O)
            self.backward(O)

            # 求gamma, shape: T * N
            molecular1 = self.alpha * self.beta
            denominator1 = np.sum(molecular1, axis=1).reshape(-1, 1)
            gamma = molecular1 / denominator1

            # 求xi, shape: T * N * N
            xi = np.zeros([T - 1, N, N])
            for t in range(T - 1):
                molecular2 = self.alpha[t, :].reshape(1, -1) * self.A.T * self.B[:, O[t + 1]].reshape(-1,
                                                                                                      1) * self.beta[
                                                                                                           t + 1,
                                                                                                           :].reshape(
                    -1, 1)
                denominator2 = np.sum(molecular2)
                xi[t, :, :] = molecular2.T / denominator2

            # 递推
            newpi = gamma[0, :]
            newA = np.sum(xi, axis=0) / np.sum(gamma[:-1, :], axis=0).reshape(-1, 1)
            newB = np.zeros(self.B.shape)
            for k in range(self.B.shape[1]):
                mask = (O == k)
                newB[:, k] = np.sum(gamma[mask, :], axis=0) / np.sum(gamma, axis=0)

            # 终止
            if np.max(abs(self.pi - newpi)) < criterion and \
                    np.max(abs(self.A - newA)) < criterion and \
                    np.max(abs(self.B - newB)) < criterion:
                break

            self.A, self.B, self.pi = newA, newB, newpi

    def viterbi(self, O):
        # 初始化T1表、T2表
        T, N = len(O), len(self.A)
        T1_table = np.zeros([N, T])
        T2_table = np.zeros([N, T])

        # 时刻1
        T1_table[:, 0] = self.pi * self.B[:, O[0]]
        T2_table[:, 0] = np.nan

        for i in range(1, T):
            # 时刻t
            curr_score = T1_table[:, i - 1].reshape(1, -1) * self.A.T * self.B[:, O[i]].reshape(-1, 1)

            # 存入T1 T2中
            T1_table[:, i] = np.max(curr_score, axis=-1)
            T2_table[:, i] = np.argmax(curr_score, axis=-1)

        # 回溯
        best_tag_id = int(np.argmax(T1_table[:, -1]))
        best_tags = [best_tag_id]
        for i in range(T - 1, 0, -1):
            best_tag_id = int(T2_table[best_tag_id, i])
            best_tags.append(best_tag_id)
        return list(reversed(best_tags))

参考资料

《统计学习方法》
https://cloud.tencent.com/developer/article/1526748
https://www.bilibili.com/video/BV1aE411o7qd?p=82