【学习笔记】统计学习方法——HMM

摘要： 隐马尔可夫模型（hidden Markov model，HMM）是可用于标注问题的统计学习模型，描述由隐藏的马尔可夫链随机生成观测序列的过程，属于生成模型。隐马尔可夫模型在语音识别，自然语言处理，生物信息，模式识别等领域有着广泛的应用。

1、马尔科夫模型与HMM

  要讲隐马尔科夫模型，需要先从马尔科夫模型讲起。已知N个有序随机变量，根据贝叶斯定理，他们的联合分布可以写成条件分布的连乘积：
$$P(x_1,x_2,\cdots ,x_N)=\prod _{n=1}^{N}P(x_n|x_{n-1},\cdots ,x_1)\text{\hspace{1em}(1)}$$
  马尔科夫模型是指，假设序列中的任何一个随机变量在给定它的前一个变量时的分布与更早的变量无关，即：
$$P(x_n|x_{n-1},\cdots ,x_1)=P(x_n|x_{n-1})\text{\hspace{1em}(2)}$$
  在此假设下，N个随机变量的联合分布可以简化为：
$$P(x_1,x_2,\cdots ,x_N)=\prod _{n=2}^{N}P(x_n|x_{n-1})\text{\hspace{1em}(3)}$$

  一阶马尔科夫性只能表达当前变量与前一个变量的关系，然而很多实际问题往往没有这么简单。为了表达当前变量与更早的变量之间的关系，可以引入高阶马尔科夫性。概括来说，M阶马尔科夫性是指当前随机变量在给定其之前的M个变量时与更早的变量无关，用公式表达就是：
$$P(x_n|x_{n-1},\cdots ,x_1)=P(x_n|x_{n-1},\cdots ,x_{n-M})\text{\hspace{1em}(4)}$$
  高阶马尔科夫性虽然达到了关联当前变量与更早的变量的目的，但有一个问题就是指数爆炸问题，即参数数量随着M的增大呈指数增长。假设每个随机变量有K种状态，对于一阶马尔科夫模型而言，要表达条件分布$P(x_n|x_{n-1})$，因为对于$x_{n-1}$的每个取值（共有K个取值）需要K个$x_n$的取值，根据条件概率的定义，这K个值的和为1，因此实际有个K-1个自由参数，因此共需要K(K-1)个参数。同理，对于M阶马尔科夫模型而言，要表达条件分布$P(x_n|x_{n-1},\cdots ,x_{n-M})$，则需要$K^M(K-1)$个参数（对于$x_n,\cdots ,x_{n-M}$的每一个取值的组合（共$K^M$个），都需要K-1个$x_n$的取值）。

  那么，有没有一种方法即能将当前变量与更早的变量关联起来，又不需要那么多参数呢？当然，这里有一种非常强大的手段，即引入隐变量，这是机器学习中用简单的基础部件表达丰富的模型的有力手段。这里如果假设隐变量构成一阶马尔科夫模型，而每个观测变量与一个隐变量关联，则可以得到一类模型的基础结构，即状态空间模型。如下图，$z_n$为隐藏变量，$x_n$为观测变量。

该类模型的关键是隐藏变量之间满足如下条件独立性，即在给定$z_n$时，$z_{n-1}$和$z_{n+1}$条件独立：

这类模型的联合分布可以表示为：
$$P(x_1,\cdots ,x_N,z_1,\cdots ,z_N)=P(z_1)$$
  可见，看似很复杂的模型被分解成了简单的$P(z_1)$、$P(z_n|z_{n-1})$和$P(x_n|z_n)$部分，这三者分别叫做初始概率模型、转移概率模型和观测概率模型，对状态空间模型建模实际就是对这三者进行建模。而且此时观测变量之间不再具有任何马尔科夫性，因为此时$x_n$的分布与其之前所有的观测变量都相关，无法从的$P(x_n|x_{n-1},\cdots ,x_1)$条件变量中拿掉任何一个变量，这样就将一个变量与其之前所有的变量关联起来了。这就是隐变量的能力。

  当$z_n$为离散变量时，该状态空间模型即为隐马尔科夫模型。现在可以解释一下“隐马尔科夫”这个名字的含义了。“马尔科夫”自然是表示随机变量之间构成一个马尔科夫链了。而“隐”是指我们要推测的变量是未知的、隐藏的。正是这些隐藏的变量构成了马尔科夫链，所以就叫“隐马尔科夫”了。

2、隐马尔可夫模型

  隐藏的马尔可夫链随机生成的状态的序列，称为状态序列；每个状态生成一个一个观测，而由此产生的观测的随机序列，称为观测序列。

  设Q是所有可能的状态的集合，V是所有可能的观测的集合：
Q = { q 1 , q 2 , ⋯   , q N } ， V = { v 1 , v 2 , ⋯   , v M } Q=\lbrace q_1,q_2,\cdots,q_N \rbrace，V=\lbrace v_1,v_2,\cdots,v_M \rbrace Q={q1​,q2​,⋯,qN​}，V={v1​,v2​,⋯,vM​}
其中，N是所有可能的状态的集合，V是所有可能的观测的集合。

  I是长度为T的状态序列，O是对应的观测序列：
$$I=(i_1,i_2,\cdots ,i_T)，O=(o_1,o_2,\cdots ,o_T)$$
  A是状态转移概率矩阵：
$$A=[a_{ij}{]}_{N\times N}$$
其中，
$$a_{ij}=P(i_{t+1}=q_j|i_t=q_i)，i=1,2,\cdots ,N;j=1,2,\cdots ,N$$
是在时刻$t$处于状态$q_i$的条件下在时刻$t+1$转移到状态$q_j$的概率。

  B是观测概率矩阵：
$$B=[b_j(k){]}_{N\times M}$$
其中，
$$b_j(k)=P(o_t=v_k|i_t=q_j)，k=1,2,\cdots ,M;j=1,2,\cdots ,N$$
是在时刻$t$处于状态$q_j$的条件下生成观测$v_k$的概率。

  $\pi$是初始状态概率向量：
$$\pi =({\pi }_i)$$
其中，
$${\pi }_i=P(i_1=q_i)，i=1,2,\cdots ,N$$
是时刻$t=1$处于状态$q_i$的概率。

3、隐马尔可夫模型的3个基本问题

1. 概率计算问题
2. 学习问题
3. 预测问题

为了解决上述的三个问题，由此产生了一些算法。

3.1、概率计算算法

• 直接计算法（计算复杂度达到了$T{N}^T$，复杂度太高，因此在实际中不使用）
• 前向算法
• 后向算法

1）前向算法：

定义：给定马尔可夫模型$\lambda$，定义到时刻$t$部分观测序列为$o_1,o_2,\cdots ,o_t$且状态为$q_i$的概率为前向概率，记作
$${\alpha }_t(i)=P(o_1,o_2,\cdots ,o_t,i_t=q_i|\lambda )$$
2）后向算法：

定义：给定马尔可夫模型$\lambda$，定义在时刻$t$状态为$q_i$的条件下，从$t+1$到$T$的 部分观测序列为$o_{t+1},o_{t+2},\cdots ,o_T$的概率为后向概率，记作
$${\beta }_t(i)=P(o_{t+1},o_{t+2},\cdots ,o_T|i_t=q_i,\lambda )$$

3.2、学习算法

• 监督学习算法
• 无监督学习算法——Baum-Welch算法（也就是EM算法）

1）监督学习算法：

使用极大似然估计法来估计模型参数
$$\begin{gathered} {\hat{a}}_{ij}=\frac{A_{ij}}{{\sum }_{j=1}^N{A}_{ij}}，i=1,2,\cdots ,N;j=1,2,\cdots ,N \\ {\hat{b}}_j(k)=\frac{B_{jk}}{{\sum }_{k=1}^M{B}_{jk}}，j=1,2,\cdots ,N;k=1,2,\cdots ,M \end{gathered}$$
2）无监督学习算法：
$$\begin{gathered} a_{ij}^{(n+1)}=\frac{{\sum }_{t=1}^{T-1}{\zeta }_t(i,j)}{{\sum }_{t=1}^{T-1}{\gamma }_t(i)} \\ {b}_j(k{)}^{(n+1)}=\frac{{\sum }_{t=1,o_t=v_k}^T{\gamma }_t(j)}{{\sum }_{t=1}^T{\gamma }_t(j)} \\ {\pi }_i^{(n+1)}={\gamma }_1(i) \end{gathered}$$

3.3、预测算法

• 近似算法
• 维特比算法

1）近似算法：
$$\begin{gathered} {\gamma }_t(i)=\frac{{\alpha }_t(i){\beta }_t(i)}{P(O|\lambda )}=\frac{{\alpha }_t(i){\beta }_t(i)}{{\sum }_{j=1}^N{\alpha }_t(j){\beta }_t(j)} \\ {i}_t^{\ast }=arg{\max }_{1\le i\le N}[{\gamma }_t(i)]，t=1,2,\cdots ,T \\ {I}^{\ast }=(i_1^{\ast },i_2^{\ast },\cdots ,i_T^{\ast }) \end{gathered}$$
2）维特比算法：

这个算法就要隆重的来学习一下了。维特比算法实际是用动态规划解隐马尔可夫模型预测问题，即用动态规划求解概率最大路径（最优路径）。（动态规划看另外一个博客）

1）如果概率最大的路径p(或者说最短路径)经过某个点，比如途中的$X_{22}$，那么这条路径上的起始点S到X22的这段子路径Q，一定是S到$X_{22}$之间的最短路径。否则，用S到$X_{22}$的最短路径R替代Q，便构成一条比P更短的路径，这显然是矛盾的。证明了满足最优性原理。

2）从S到E的路径必定经过第i个时刻的某个状态，假定第i个时刻有k个状态，那么如果记录了从S到第i个状态的所有k个节点的最短路径，最终的最短路径必经过其中一条，这样，在任意时刻，只要考虑非常有限的最短路即可。

3）结合以上两点，假定当我们从状态i进入状态i+1时，从S到状态i上各个节的最短路径已经找到，并且记录在这些节点上，那么在计算从起点S到第i+1状态的某个节点$X_{i+1}$的最短路径时，只要考虑从S到前一个状态i所有的k个节点的最短路径，以及从这个节点到$X_{i+1}$，j的距离即可。

4、代码实践

class HiddenMarkov:
    def forward(self, Q, V, A, B, O, PI):  # 使用前向算法
        N = len(Q)  #可能存在的状态数量
        M = len(O)  # 观测序列的大小
        alphas = np.zeros((N, M))  # alpha值
        T = M  # 有几个时刻，有几个观测序列，就有几个时刻
        for t in range(T):  # 遍历每一时刻，算出alpha值
            indexOfO = V.index(O[t])  # 找出序列对应的索引
            for i in range(N):
                if t == 0:  # 计算初值
                    alphas[i][t] = PI[t][i] * B[i][indexOfO]  # P176（10.15）
                    print(
                        'alpha1(%d)=p%db%db(o1)=%f' % (i, i, i, alphas[i][t]))
                else:
                    alphas[i][t] = np.dot(
                        [alpha[t - 1] for alpha in alphas],
                        [a[i] for a in A]) * B[i][indexOfO]  # 对应P176（10.16）
                    print('alpha%d(%d)=[sigma alpha%d(i)ai%d]b%d(o%d)=%f' %
                          (t, i, t - 1, i, i, t, alphas[i][t]))
                    # print(alphas)
        P = np.sum([alpha[M - 1] for alpha in alphas])  # P176(10.17)
        # alpha11 = pi[0][0] * B[0][0]    #代表a1(1)
        # alpha12 = pi[0][1] * B[1][0]    #代表a1(2)
        # alpha13 = pi[0][2] * B[2][0]    #代表a1(3)

    def backward(self, Q, V, A, B, O, PI):  # 后向算法
        N = len(Q)  # 可能存在的状态数量
        M = len(O)  # 观测序列的大小
        betas = np.ones((N, M))  # beta
        for i in range(N):
            print('beta%d(%d)=1' % (M, i))
        for t in range(M - 2, -1, -1):
            indexOfO = V.index(O[t + 1])  # 找出序列对应的索引
            for i in range(N):
                betas[i][t] = np.dot(
                    np.multiply(A[i], [b[indexOfO] for b in B]),
                    [beta[t + 1] for beta in betas])
                realT = t + 1
                realI = i + 1
                print(
                    'beta%d(%d)=[sigma a%djbj(o%d)]beta%d(j)=(' %
                    (realT, realI, realI, realT + 1, realT + 1),
                    end='')
                for j in range(N):
                    print(
                        "%.2f*%.2f*%.2f+" % (A[i][j], B[j][indexOfO],
                                             betas[j][t + 1]),
                        end='')
                print("0)=%.3f" % betas[i][t])
        # print(betas)
        indexOfO = V.index(O[0])
        P = np.dot(
            np.multiply(PI, [b[indexOfO] for b in B]),
            [beta[0] for beta in betas])
        print("P(O|lambda)=", end="")
        for i in range(N):
            print(
                "%.1f*%.1f*%.5f+" % (PI[0][i], B[i][indexOfO], betas[i][0]),
                end="")
        print("0=%f" % P)

    def viterbi(self, Q, V, A, B, O, PI):
        N = len(Q)  #可能存在的状态数量
        M = len(O)  # 观测序列的大小
        deltas = np.zeros((N, M))
        psis = np.zeros((N, M))
        I = np.zeros((1, M))
        for t in range(M):
            realT = t + 1
            indexOfO = V.index(O[t])  # 找出序列对应的索引
            for i in range(N):
                realI = i + 1
                if t == 0:
                    deltas[i][t] = PI[0][i] * B[i][indexOfO]
                    psis[i][t] = 0
                    print('delta1(%d)=pi%d * b%d(o1)=%.2f * %.2f=%.2f' %
                          (realI, realI, realI, PI[0][i], B[i][indexOfO],
                           deltas[i][t]))
                    print('psis1(%d)=0' % (realI))
                else:
                    deltas[i][t] = np.max(
                        np.multiply([delta[t - 1] for delta in deltas],
                                    [a[i] for a in A])) * B[i][indexOfO]
                    print(
                        'delta%d(%d)=max[delta%d(j)aj%d]b%d(o%d)=%.2f*%.2f=%.5f'
                        % (realT, realI, realT - 1, realI, realI, realT,
                           np.max(
                               np.multiply([delta[t - 1] for delta in deltas],
                                           [a[i] for a in A])), B[i][indexOfO],
                           deltas[i][t]))
                    psis[i][t] = np.argmax(
                        np.multiply(
                            [delta[t - 1] for delta in deltas],
                            [a[i]
                             for a in A])) + 1  #由于其返回的是索引，因此应+1才能和正常的下标值相符合。
                    print('psis%d(%d)=argmax[delta%d(j)aj%d]=%d' %
                          (realT, realI, realT - 1, realI, psis[i][t]))
        print(deltas)
        print(psis)
        I[0][M - 1] = np.argmax([delta[M - 1] for delta in deltas
                                 ]) + 1  #由于其返回的是索引，因此应+1才能和正常的下标值相符合。
        print('i%d=argmax[deltaT(i)]=%d' % (M, I[0][M - 1]))
        for t in range(M - 2, -1, -1):
            I[0][t] = psis[int(I[0][t + 1]) - 1][t + 1]
            print('i%d=psis%d(i%d)=%d' % (t + 1, t + 2, t + 2, I[0][t]))
        print("状态序列I：", I)

Reference:

• https://www.zhihu.com/people/yi-zhi-luo-shui-de-zhu
• 统计学习方法（第二版）
• 数学之美
• https://github.com/fengdu78/lihang-code