经典算法: 隐马尔可夫模型，HMM

1. 引言

隐马尔可夫模型，Hidden markov model, HMM,是一种关于时序的概率模型。

例子:
假设我手里有三个不同的骰子。第一个骰子是我们平常见的骰子（称这个骰子为D6），6个面，每个面（1，2，3，4，5，6）出现的概率是1/6。第二个骰子是个四面体（称这个骰子为D4），每个面（1，2，3，4）出现的概率是1/4。第三个骰子有八个面（称这个骰子为D8），每个面（1，2，3，4，5，6，7，8）出现的概率是1/8。

假设我们开始掷骰子，我们先从三个骰子里挑一个，挑到每一个骰子的概率都是1/3。然后我们掷骰子，得到一个数字，1，2，3，4，5，6，7，8中的一个。不停的重复上述过程，我们会得到一串数字，每个数字都是1，2，3，4，5，6，7，8中的一个。例如我们可能得到这么一串数字（掷骰子10次）：1 6 3 5 2 7 3 5 2 4

这串数字叫做可见状态链。但是在隐马尔可夫模型中，我们不仅仅有这么一串可见状态链，还有一串隐含状态链。在这个例子里，这串隐含状态链就是你用的骰子的序列。
比如，隐含状态链有可能是：D6 D8 D8 D6 D4 D8 D6 D6 D4 D8

一般来说，HMM中说到的马尔可夫链其实是指隐含状态链，因为隐含状态（骰子）之间存在转换概率（transition
probability）。在我们这个例子里，D6的下一个状态是D4，D6，D8的概率都是1/3。D4，D8的下一个状态是D4，D6，D8的转换概率也都一样是1/3。这样设定是为了最开始容易说清楚，但是我们其实是可以随意设定转换概率的。比如，我们可以这样定义，D6后面不能接D4，D6后面是D6的概率是0.9，是D8的概率是0.1。这样就是一个新的HMM。

同样的，尽管可见状态之间没有转换概率，但是隐含状态和可见状态之间有一个概率叫做输出概率（emission probability）。
就我们的例子来说，六面骰（D6）产生1的输出概率是1/6。产生2，3，4，5，6的概率也都是1/6。我们同样可以对输出概率进行其他定义。比如，我有一个被赌场动过手脚的六面骰子，掷出来是1的概率更大，是1/2，掷出来是2，3，4，5，6的概率是1/10。

隐马尔可夫模型是生成模型，可以看作是有向图中的一种序列模型。

隐马尔可夫模型，有3个基本问题：

概率计算问题：求已知观察序列的发生概率；

学习问题：求模型的参数；

预测问题：求状态序列；

下面将分别介绍这三种问题的算法:

2. 问题1：概率计算算法

2.1 直接计算法

计算量过大，应用少。

2.2 前向算法

思想：
基于上一观察结果求出的概率，来计算下一观察序列发生概率。
前向算法，是一种递推的思想。

2.3 后向算法

后向算法的递归公式推导：

3. 问题2：学习算法

监督学习：

Baum-Welch算法：

4. 问题3：预测算法

4.1 近似算法

4.2 维特比算法（Viterbi）

着重于递归。

思想：

在前一步$S_i$计算的概率是最大的前提下，后一步$S_{i+1}$要想得到最大概率，前面的状态肯定不变。
通过，计算从前一步Si到现在这步$S_{i+1}$的转移概率$a_{i(i+1)}$及在该状态下发生观察结果的概率$B_{i+1}$，来确定这一步的状态。

5. 实践

"""
隐马尔科夫模型
三类问题：1.概率计算 2.学习问题（参数估计） 3.预测问题（状态序列的预测）
"""
import numpy as np
from itertools import accumulate

class GenData:
    """
    根据隐马尔科夫模型生成相应的观测数据
    """

    def __init__(self, hmm, n_sample):
        self.hmm = hmm
        self.n_sample = n_sample

    def _locate(self, prob_arr):
        # 给定概率向量，返回状态
        seed = np.random.rand(1)
        for state, cdf in enumerate(accumulate(prob_arr)):
            if seed <= cdf:
                return state
        return

    def init_state(self):
        # 根据初始状态概率向量，生成初始状态
        return self._locate(self.hmm.S)

    def state_trans(self, current_state):
        # 转移状态
        return self._locate(self.hmm.A[current_state])

    def gen_obs(self, current_state):
        # 生成观测
        return self._locate(self.hmm.B[current_state])

    def gen_data(self):
        # 根据模型产生观测数据
        current_state = self.init_state()
        start_obs = self.gen_obs(current_state)
        state = [current_state]
        obs = [start_obs]
        n = 0
        while n < self.n_sample - 1:
            n += 1
            current_state = self.state_trans(current_state)
            state.append(current_state)
            obs.append(self.gen_obs(current_state))
        return state, obs

class HMM:
    def __init__(self, n_state, n_obs, S=None, A=None, B=None):
        self.n_state = n_state  # 状态的个数n
        self.n_obs = n_obs  # 观测的种类数m
        self.S = S  # 1*n, 初始状态概率向量
        self.A = A  # n*n, 状态转移概率矩阵
        self.B = B  # n*m, 观测生成概率矩阵

    def _alpha(self, obs, t):
        # 计算时刻t各个状态的前向概率
        b = self.B[:, obs[0]]
        alpha = np.array([self.S * b])  # n*1
        for i in range(1, t + 1):
            alpha = (alpha @ self.A) * np.array([self.B[:, obs[i]]])
        return alpha[0]

    def forward_prob(self, obs):
        # 前向算法计算最终生成观测序列的概率, 即各个状态下概率之和
        alpha = _alpha(self, obs, len(obs) - 1)
        return np.sum(alpha)

    def _beta(self, obs, t):
        # 计算时刻t各个状态的后向概率
        beta = np.ones(self.n_state)
        for i in reversed(range(t + 1, len(obs))):
            beta = np.sum(self.A * self.B[:, obs[i]] * beta, axis=1)
        return beta

    def backward_prob(self, obs):
        # 后向算法计算生成观测序列的概率
        beta = _beta(self, obs, 0)
        return np.sum(self.S * self.B[:, obs[0]] * beta)

    def fb_prob(self, obs, t=None):
        # 将前向和后向合并
        if t is None:
            t = 0
        res = _alpha(self, obs, t) * _beta(self, obs, t)
        return res.sum()

    def _gamma(self, obs, t):
        # 计算时刻t处于各个状态的概率
        alpha = _alpha(self, obs, t)
        beta = _beta(self, obs, t)
        prob = alpha * beta
        return prob / prob.sum()


    def point_prob(self, obs, t, i):
        # 计算时刻t处于状态i的概率
        prob = _gamma(self, obs, t)
        return prob[i]


    def _xi(self, obs, t):
        alpha = np.mat(_alpha(self, obs, t))
        beta_p = _beta(self, obs, t + 1)
        obs_prob = self.B[:, obs[t + 1]]
        obs_beta = np.mat(obs_prob * beta_p)
        alpha_obs_beta = np.asarray(alpha.T * obs_beta)
        xi = alpha_obs_beta * self.A
        return xi / xi.sum()

    def fit(self, obs_data, maxstep=100):
        # 利用Baum-Welch算法学习
        self.A = np.ones((self.n_state, self.n_state)) / self.n_state
        self.B = np.ones((self.n_state, self.n_obs)) / self.n_obs
        self.S = np.random.sample(self.n_state)  # 初始状态概率矩阵（向量），的初始化必须随机状态，否则容易陷入局部最优
        self.S = self.S / self.S.sum()
        step = 0
        while step < maxstep:
            xi = np.zeros_like(self.A)
            gamma = np.zeros_like(self.S)
            B = np.zeros_like(self.B)
            S = _gamma(self, obs_data, 0)
            for t in range(len(obs_data) - 1):
                tmp_gamma = _gamma(self, obs_data, t)
                gamma += tmp_gamma
                xi += _xi(self, obs_data, t)
                B[:, obs_data[t]] += tmp_gamma

            # 更新 A
            for i in range(self.n_state):
                self.A[i] = xi[i] / gamma[i]
            # 更新 B
            tmp_gamma_end = _gamma(self, obs_data, len(obs_data) - 1)
            gamma += tmp_gamma_end
            B[:, obs_data[-1]] += tmp_gamma_end
            for i in range(self.n_state):
                self.B[i] = B[i] / gamma[i]
            # 更新 S
            self.S = S
            step += 1
        return self

    def predict(self, obs):
        # 采用Viterbi算法预测状态序列
        N = len(obs)
        nodes_graph = np.zeros((self.n_state, N), dtype=int)  # 存储时刻t且状态为i时， 前一个时刻t-1的状态，用于构建最终的状态序列
        delta = self.S * self.B[:, obs[0]]  # 存储到t时刻，且此刻状态为i的最大概率
        nodes_graph[:, 0] = range(self.n_state)
        for t in range(1, N):
            new_delta = []
            for i in range(self.n_state):
                temp = [self.A[j, i] * d for j, d in enumerate(delta)]  # 当前状态为i时， 选取最优的前一时刻状态
                max_d = max(temp)
                new_delta.append(max_d * self.B[i, obs[t]])
                nodes_graph[i, t] = temp.index(max_d)
            delta = new_delta

        current_state = np.argmax(nodes_graph[:, -1])
        path = []
        t = N
        while t > 0:
            path.append(current_state)
            current_state = nodes_graph[current_state, t - 1]
            t -= 1
        return list(reversed(path))


if __name__ == '__main__':
    # 预测
    S = np.array([0.5, 0.5])
    A = np.array([[0.8, 1], [0.8, 0.8]])
    B = np.array([[0.2, 0.0, 0.8], [0, 0.8, 0.2]])
    hmm = HMM(2, 3, S, A, B)
    g = GenData(hmm, 200)
    state, obs = g.gen_data()
    print(obs)
    path = hmm.predict(obs)
    score = sum([int(i == j) for i, j in zip(state, path)])
    print(score / len(path))

参考：

1. 作者：Yang Eninala zhihu;
2. 隐马尔科夫模型(HMM)以及Python实现;