机器学习之概率图模型

---

参考资料：

1. 《机器学习》——周志华
2. 《统计学习方法》——李航
3. https://blog.csdn.net/continueoo/article/details/77893587
4. https://blog.csdn.net/continueoo/article/details/77893587

0. 背景介绍

判别式模型,对条件分布进行建模；生成式模型，对联合分布进行建模。概率图模型是一类利用图来表达变量间相关关系的概率模型，既包含判别式模型，也包含生成式模型。

1. 隐马尔科夫模型（Hidden Markov Model, HMM）

1.1. HMM的数学定义

HMM是关于时序的概率模型，描述由一个隐藏的马尔科夫链随机生成不可观测的状态随机序列，再由各个状态生成一个可观测的随机序列的过程。

它是一种生成式模型，因为它需要先生成状态变量，然后在状态的条件下生成观测变量，根据条件概率的公式，等价于生成状态变量和观测变量的联合概率。所以是生成式模型。

HMM模型常用于序列标注问题建模，假设时间长度为n，
对于$i=1,2,\cdots ,n$时刻，HMM中的变量可以分为两组，用 X = { x 1 , x 2 , ⋯   , x n } , x i ∈ O = { o 1 , o 2 , ⋯   , o M } X=\{x_1,x_2,\cdots,x_n\},x_i\in O=\{o_1,o_2,\cdots,o_M\} X={x1​,x2​,⋯,xn​},xi​∈O={o1​,o2​,⋯,oM​}表示观测变量序列（每个观测变量有$M$个可能的取值）； Y = { y 1 , y 2 , ⋯   , y n } , y i ∈ S = { s 1 , s 2 , ⋯   , s N } Y=\{y_1,y_2,\cdots,y_n\},y_i\in S=\{s_1,s_2,\cdots,s_N\} Y={y1​,y2​,⋯,yn​},yi​∈S={s1​,s2​,⋯,sN​}（每个状态变量有$N$个可能的取值）表示状态变量序列，也叫隐变量。
并且有如下假设：在任一时刻，观测变量的取值仅仅依赖于该时刻的状态变量，即$x_t$仅由$y_t$决定；
$$\begin{gathered} P(x_t|y_T,x_T,y_{T-1},x_{T-1},\cdots ,y_1,x_1)=P(x_t|y_t)) \\ \text{\hspace{1em}(1-1)} \end{gathered}$$
同时，$y_t$仅依赖于$y_{t-1}$，即状态变量序列是一个马尔科夫链。由此可以得到所有变量的联合概率分布为（注意根据假设，条件概率的条件省去了无关变量）：
$$\begin{gathered} P(x_1,y_1,\cdots ,x_n,y_n)=P(y_1)P(x_1|y_1)\prod _{i=2}^{n}P(y_i|y_{i-1})P(x_i|y_i) \\ \text{\hspace{1em}(1-2)} \end{gathered}$$
除了前面提到的状态空间$Y$和观测空间$X$，要想确定一个HMM模型还需要如下三组参数：

• 状态转移概率矩阵（隐藏状态到隐藏状态之间的概率转换）：$A=[a_{ij}{]}_{N\times N}$，其中$a_{ij}=P(y_{t+1}=s_j|y_t=s_i),1\le i,j\le N$.为了便于后文公式推导中的符号说明，我们引入$a_{s_i{s}_j}=P(y_{t+1}=s_j|y_t=s_i),1\le i,j\le N$。有时候也用省略的写法$a_{ij}$
• 观测概率矩阵（也叫发射概率，隐藏状态到观测状态之间的概率转换）：$B=[b_{ij}{]}_{N\times M}$，其中$b_{ij}$表示$P(x_t=o_j|y_t=s_i),1\le i\le N,1\le j\le M$.为了便于后文公式推导中的符号说明，我们引入$b_{s_i{o}_j}$表示$P(x_t=o_j|y_t=s_i),1\le i\le N,1\le j\le M$。有时候也用省略的写法$b_{ij}$
• 初始状态概率向量：$\pi =({\pi }_1,{\pi }_2,\cdots ,{\pi }_N)$，其中${\pi }_i=P(y_1=s_i)$

记一个HMM为$\lambda =[A,B,\pi ]$，则可按照如下过程产生观测序列：

1. 设置t=1，根据初始状态概率$\pi$选择初始状态$y_1$（比如根据概率最大的原则选取）;
2. 根据状态$y_t$和输出观测概率矩阵$B$确定观测变量$x_t$；
3. 根据状态$y_t$和状态转移矩阵$A$确定$y_{t+1}$；
4. 若$t<n$，设置$t=t+1$，并跳转到2；否则停止。

1.2. HMM的实际应用问题描述

在实际应用中，我们常常关注HMM的下面三个问题：

1. 观测序列概率计算问题：给定模型$\lambda =[A,B,\pi ]$和观测序列 X = { x 1 , x 2 , ⋯   , x n } X=\{x_1,x_2,\cdots,x_n\} X={x1​,x2​,⋯,xn​}，如何有效计算出观测序列出现的概率$P(X|\lambda )$。换言之，如何评估模型与观测序列之间的匹配程度（因为计算出观测序列的时候，可能会有多个选择要选择最匹配的那个结构）。例如，根据时间序列的历史值预测未来值（已知n时刻及之前的x，对n+1时刻可能的x都拿来计算观测序列的概率，挑使得观测序列概率最大的状态作为n+1时刻的x）；这个过程称为前向算法
2. 状态序列预测问题算法：给定模型$\lambda =[A,B,\pi ]$和观测序列 X = { x 1 , x 2 , ⋯   , x n } X=\{x_1,x_2,\cdots,x_n\} X={x1​,x2​,⋯,xn​}，如何找到与此序观测序列最匹配的状态序列 Y = { y 1 , y 2 , ⋯   , y n } , y i ∈ { s 1 , s 2 , ⋯   , s N } Y=\{y_1,y_2,\cdots,y_n\},y_i\in\{s_1,s_2,\cdots,s_N\} Y={y1​,y2​,⋯,yn​},yi​∈{s1​,s2​,⋯,sN​}。换言之，如何根据观测序列推断出隐藏状态。例如，在语音识别ASR中，如何根据语音信号（观测值）确定文字序列（隐藏状态）；主流的算法是维特比算法
3. 参数学习问题：给定观测序列 X = { x 1 , x 2 , ⋯   , x n } X=\{x_1,x_2,\cdots,x_n\} X={x1​,x2​,⋯,xn​}，如何调整模型参数$\lambda =[A,B,\pi ]$使得该序列出现的概率$P(X|\lambda )$最大。换言之，如何训练模型使其能最好滴描述观测数据。例如，根据梯度下降的方式在序列样本上学习得到最优的模型参数；EM算法 前后向算法

1.2.1. 观测序列概率计算问题

给定模型$\lambda =[A,B,\pi ]$和观测序列 X = { x 1 , x 2 , ⋯   , x n } X=\{x_1,x_2,\cdots,x_n\} X={x1​,x2​,⋯,xn​}，如何有效计算出观测序列出现的概率$P(X|\lambda )$。

1.2.1.1. 直接暴力计算法

此方法是最简单粗暴的算法，但是复杂度过高，实际中一般不采用。其他算法一般都是在此算法上优化得到的，因此有必要了解一下直接计算法的原理。
直接计算法的思路是，通过列举所有可能的长度为T的状态序列 Y = { y 1 , y 2 , ⋯   , y T } Y=\{y_1,y_2,\cdots,y_T\} Y={y1​,y2​,⋯,yT​}，求各个状态序列Y与观测序列 X = { x 1 , x 2 , ⋯   , x T } X=\{x_1,x_2,\cdots,x_T\} X={x1​,x2​,⋯,xT​}的联合概率$P(X,Y|\lambda )$，然后利用全概率对所有可能的状态序列求和得到$P(X|\lambda )$。
对于某个具体的状态序列 Y = { y 1 , y 2 , ⋯   , y T } Y=\{y_1,y_2,\cdots,y_T\} Y={y1​,y2​,⋯,yT​}，其概率为：
$$P(Y|\lambda )={\pi }_{y_1}{a}_{y_1{y}_2}{a}_{y_2{y}_3}\cdots a_{y_{T-1}{y}_T}\text{\hspace{1em}(1-3)}$$
给定上述状态序列，观测序列 X = { x 1 , x 2 , ⋯   , x T } X=\{x_1,x_2,\cdots,x_T\} X={x1​,x2​,⋯,xT​}的概率为：
$$P(X|Y,\lambda )=b_{y_1,x_1}{b}_{y_2,x_2}\cdots b_{y_T{x}_T}\text{\hspace{1em}(1-4)}$$
则状态序列和观测序列的联合概率可以表示为：
$$\begin{array}{rl}P(X,Y|\lambda )& =P(X|Y,\lambda )P(Y|\lambda )\\ & ={\pi }_{y_1}{b}_{y_1{x}_1}{a}_{y_1{y}_2}{b}_{y_2{x}_2}{a}_{y_2{y}_3}\cdots b_{y_{T-1}{x}_{T-1}}{a}_{y_{T-1}{y}_T}{b}_{y_T{x}_T}\end{array}\text{\hspace{1em}(1-5)}$$
则观测序列的概率可以有全概率公式求得:
$$\begin{array}{rl}P(X|\lambda )& =\sum _{Y}P(X|Y,\lambda )P(Y|\lambda )\\ & =\sum _{y_1\in S,y_2\in S,\cdots ,y_T\in S}{\pi }_{y_1}{b}_{y_1{x}_1}{a}_{y_1{y}_2}{b}_{y_2{x}_2}{a}_{y_2{y}_3}\cdots b_{y_{T-1}{x}_{T-1}}{a}_{y_{T-1}{y}_T}{b}_{y_T{x}_T}\end{array}\text{\hspace{1em}(1-6)}$$
上式的计算复杂度高达$O(T{N}^T)$

1.2.1.2. 基于动态规划的前向算法

直接计算法复杂度过高的原因在于存在大量冗余计算。例如，当t时刻之前的状态序列固定时，穷举t时刻之后的状态序列的各种取值情况时，都会重复计算t时刻之前的各种取值情况。因此可以采用动态规划的思想，求得计算过程中的递推式，减少重复计算。
定义到t时刻时，观测序列 { x 1 , x 2 , ⋯   , x t } \{x_1,x_2,\cdots,x_t\} {x1​,x2​,⋯,xt​}且状态为$s_i$的概率为前向概率，记作：
$${\alpha }_t(s_i)=P(x_1,x_2,\cdots ,x_t,y_t=s_i|\lambda )\text{\hspace{1em}(1-7)}$$
进行如下递推计算过程：
初值t=1：
$${\alpha }_1(s_i)={\pi }_{s_i}{b}_{s_i{x}_i},i=1,2,\cdots ,N\text{\hspace{1em}(1-8)}$$
对于$t=1,2,\cdots ,T-1$:
$${\alpha }_{t+1}(s_i)=\left[\sum _{j=1}^N{\alpha }_t(s_j)a_{s_j{s}_i}\right]b_{s_i{x}_{t+1}},i=1,2,\cdots ,N\text{\hspace{1em}(1-9)}$$
终止值：
$$P(O|\lambda )=\sum _{i=1}^N{\alpha }_T(s_i)\text{\hspace{1em}(1-10)}$$

1.2.1.3. 《统计学习方法》例10.2python代码

项目地址：https://github.com/lankuohsing/machine-learning-in-python/tree/main/Probabilistic-Graphical-Model

# 状态 1 2 3之间的转移概率矩阵
A = [[0.5,0.2,0.3],
     [0.3,0.5,0.2],
     [0.2,0.3,0.5]]
# 初始状态概率
pi = [0.2,0.4,0.4]
# 每个袋子里，红白求的概率
# red white
B = [[0.5,0.5],
     [0.4,0.6],
     [0.7,0.3]]

# In[]
#前向算法
def hmm_forward(A,B,pi,O):
    T = len(O)# 观测序列长度
    N = len(A[0])# 状态个数
    #step1 初始化
    alpha = [[0]*T for _ in range(N)]# 每行代表不同的状态，每列代表不同的观测时刻
    for i in range(N):
        alpha[i][0] = pi[i]*B[i][O[0]]

    #step2 计算alpha(t)
    for t in range(1,T):
        for i in range(N):
            temp = 0
            for j in range(N):
                temp += alpha[j][t-1]*A[j][i]
            alpha[i][t] = temp*B[i][O[t]]

    #step3
    proba = 0
    for i in range(N):
        proba += alpha[i][-1]
    return proba,alpha

A = [[0.5,0.2,0.3],
     [0.3,0.5,0.2],
     [0.2,0.3,0.5]]
B = [[0.5,0.5],
     [0.4,0.6],
     [0.7,0.3]]
pi = [0.2,0.4,0.4]
O = [0,1,0]
proba,alpha=hmm_forward(A,B,pi,O)  #结果为 0.130218
print(proba)

1.2.2. 状态序列预测问题算法

给定模型$\lambda =[A,B,\pi ]$和观测序列 X = { x 1 , x 2 , ⋯   , x n } X=\{x_1,x_2,\cdots,x_n\} X={x1​,x2​,⋯,xn​}，如何找到与此序观测序列最匹配的状态序列 Y = { y 1 , y 2 , ⋯   , y n } , y i ∈ { s 1 , s 2 , ⋯   , s N } Y=\{y_1,y_2,\cdots,y_n\},y_i\in\{s_1,s_2,\cdots,s_N\} Y={y1​,y2​,⋯,yn​},yi​∈{s1​,s2​,⋯,sN​}。

1.2.2.1. 直接暴力计算法

穷举所有的状态序列的可能情况，分别计算它们的概率并求最大。此方法时间复杂度过高，基本没有实用价值

1.2.2.2. Viterbi算法

（注：下文的起点，都是默认指0时刻从起点出发的意思）
该算法是基于动态规划的最短路径算法，在其他最短路径问题中也具有广泛应用。它的思想是：记末尾时刻为$T$，如果最优路径在时刻t通过结点$i_t^{\ast }$，那么这条路径从结点$i_t^{\ast }$到终点$i_T^{\ast }$的部分路径，对于从$i_t^{\ast }$到$i_T^{\ast }$的所有子路径来说必须是最优的，同理，这条路径中从起点到结点$i_t^{\ast }$的部分路径也是从起点到结点$i_t^{\ast }$的所有子路径中最优的。根据这一原理，我们只需从时刻$t=1$开始，递推地计算从起点到在时刻t状态为i的结点的各条子路径中的最优路径，并记录下每个结点的上一时刻的最优结点，当计算到终点时就完成了整个最优路径的计算。

定义从起点到时刻t状态为i的所有子路径$(i_1,i_2,\cdots ,i_t)$中概率最大值为：
$${\delta }_t(i)=\underset{i_1,i_2,\cdots ,i_{t-1}}{\max }P(i_t=i,i_{t-1},\cdots ,i_1,x_t,\cdots ,x_1|\lambda ),i=1,2,\cdots ,N\text{\hspace{1em}(1-11)}$$
那么可以得到${\delta }_{t+1}(i)$的递推式：
$$\begin{array}{rl}{\delta }_{t+1}(i)& =\underset{i_1,i_2,\cdots ,i_t}{\max }P(i_t=i,i_t,\cdots ,i_1,x_{t+1},\cdots ,x_1|\lambda )\\ & =\underset{1\le j\le M}{\max }[{\delta }_t(j)a_{ji}]b_{i{x}_{t+1}},i=1,2,\cdots ,N;t=1,2,\cdots ,T-1\end{array}\text{\hspace{1em}(1-12)}$$
定义从起点到在时刻t状态为i的结点所有单个路径 { ( i 1 , i 2 , ⋯   , i t − 1 , t ) } \{(i_1,i_2,\cdots,i_{t-1},t)\} {(i1​,i2​,⋯,it−1​,t)}中概率最大的路径的第$t-1$时刻的结点为：
$${\psi }_t(i)=arg\underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]\text{\hspace{1em}(1-13)}$$

那么viterbi算法的流程就呼之欲出了：
(1) 赋予初值：
$$\begin{array}{rl}{\delta }_1(i)& ={\pi }_i{b}_{i{x}_1},i=1,2,\cdots ,N\\ {\psi }_i(i)& =0,i=1,2,\cdots ,N\end{array}\text{\hspace{1em}(1-14)}$$
(2) 递推式,$t=2,3,\cdots ,T$：
$$\begin{array}{rl}{\delta }_t(i)& =\underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}]b_{i{x}_{t+1}},i=1,2,\cdots ,N\\ {\psi }_t(i)& =arg\underset{1\le j\le N}{\max }[{\delta }_{t-1}(j)a_{ji}],i=1,2,\cdots ,N\end{array}\text{\hspace{1em}(1-15)}$$
(3) 终止，计算完整的最优路径和概率：
$$\begin{array}{rl}{P}^{\ast }& =\underset{1\le i\le N}{\max }{\delta }_T(i)\\ i_T^{\ast }& =arg\underset{1\le j\le N}{\max }[{\delta }_t(i)]\\ i_t^{\ast }& ={\psi }_{t+1}(i_{t+1}^{\ast }),t=T-1,T-2,\cdots ,1\end{array}\text{\hspace{1em}(1-16)}$$

1.3. HMM实际例子

1.4. 马尔科夫随机场（Markov Random Field, MRF）

MRF是典型的马尔科夫网，是一种著名的无向图模型，图中每个结点表示一个或一组变量，结点之间的边表示两个变量之间的依赖关系。定义势函数（亦称“因子”）为变量子集上的非负实函数，主要用于定义概率分布函数。
对于结点的一个子集，若其中任意两点结点之间都有边连接，则称该节点子集为一个“团”。若在一个团中加入任何一个节点都不再形成团，则称该团为“极大团”。对于n各变量 X = { x 1 , x 2 , ⋯   , x n } X=\{x_1,x_2,\cdots,x_n\} X={x1​,x2​,⋯,xn​}，记其所有团构成的集合为$C$，与团$Q\in C$对应的变量集合记为$X_Q$，则联合概率$P(X)$定义为
$$P(X)=\frac{1}{Z}\prod _{Q\in C}\psi (X_Q)\text{\hspace{1em}(1-2)}$$
其中$\psi (X_Q)$为与团$Q$对应的是函数，$Z={\sum }_X{\prod }_{Q\in C}\psi (X_Q)$ 为规范化因子，用以确保上式的值域为0~1.
为了减少需要计算的团的数目，有如下推论：若$Q$不是极大团，则它必被一个极大团$Q^{\ast }$，于是$P(X)$可以基于极大团来定义。假设所有极大团构成的集合为$C^{\ast }$，则有
$$P(X)=\frac{1}{Z^{\ast }}\prod _{Q\in C^{\ast }}\psi (X_Q)\text{\hspace{1em}(1-3)}$$
假设结点集合A中的结点到结点集B中的点都必须经过结点集C中的结点，则称结点集A和B被结点集C分离，C称为“分离集”。对于马尔科夫随机场，有“全局马尔科夫性”：给定两个变量子集的分离集，则这两个变量子集条件独立，也即有
$$P(x_A,x_B|x_C)=P(x_A|x_C)P(x_B|x_C)\text{\hspace{1em}(1-4)}$$
具体证明见《机器学习》（周志华）324.
由全局马尔科夫性可以得到两个推论：

1. 局部马尔科夫性：给定某变量的邻接变量，则该变量条件独立与其他变量
2. 成对马尔科夫性：给定所有其他变量，两个非邻接变量条件独立。

2. 条件随机场（Conditional Random Field， CRF）

2.1. CRF的数学定义

令 X = { x 1 , x 2 , ⋯   , x n } , x i ∈ { o 1 , o 2 , ⋯   , o M } X=\{x_1,x_2,\cdots,x_n\},x_i\in \{o_1,o_2,\cdots,o_M\} X={x1​,x2​,⋯,xn​},xi​∈{o1​,o2​,⋯,oM​}为观测序列， Y = { y 1 , y 2 , ⋯   , y n } , y i ∈ { s 1 , s 2 , ⋯   , s N } Y=\{y_1,y_2,\cdots,y_n\},y_i\in\{s_1,s_2,\cdots,s_N\} Y={y1​,y2​,⋯,yn​},yi​∈{s1​,s2​,⋯,sN​}为状态序列（标记或者叫标签），CRF的目的是构建条件概率模型$P(Y|X)$。注意，$Y$的各分量之间可以是相关的，例如nlp的词性标注任务中，观测数据为句子的单词序列，标记数据为相应的词性序列。
令$G=⟨V,E⟩$表示结点与标记变量$Y$中元素一一对应的无向图，$y_v$表示与结点$v$对应的标记变量，$n(v)$表示结点$v$的邻接节点，若图$G$中每个变量$y_v$都满足马尔科夫性，即
P ( y v ∣ X , Y V \ { v } ) = P ( y v ∣ X , Y n ( v ) ) (2-1) P\left(y_v|X,Y_{V\backslash\{v\}}\right)=P(y_v|X,Y_n(v))\tag{2-1} P(yv​∣X,YV\{v}​)=P(yv​∣X,Yn​(v))(2-1)
则$(Y,X)$构成一个条件随机场。

以链式条件随机场为例，CRF的条件概率被定义为
$$P(Y|X)=\frac{1}{Z}exp\left(\sum _j\sum _{i=1}^{n-1}{\lambda }_j{t}_j(y_{i+1},y_i,X,i)+\sum _k\sum _{i=1}^n{\mu }_k{s}_k(y_i,X,i)\right)\text{\hspace{1em}(2-2)}$$
其中$t_j(y_{i+1},y_i,X,i)$是定义在观测序列的两个相邻标记位置上的转移特征函数，用于刻画向量标记变量之间的相关关系以及观测序列对它们的影响；$s_k(y_i,X,i)$是定义在观测序列的标记位置$i$上的状态特征函数，用于刻画观测序列对标记变量的影响；${\lambda }_j$和${\mu }_k$为参数，$Z$为规范化因子。

2.2. 特征函数的例子

以词性标注为例，
$t_j(y_{i+1},y_i,X,i)=1,if{y}_{i+1}=[P],y_i=[V]and{x}_i="knock";0otherwise$

$s_k(y_i,X,i)=1,if{y}_i=[V]and{x}_i="knock";0otherwise$