隐马尔科夫模型

概率图模型是一类用图来表达变量相关关系的概率模型

概率图模型一般可以分为两种：使用有向无环图DAG的称为贝叶斯网，使用无向图的称为马尔可夫网

隐马尔科夫模型

隐马尔科夫模型(Hidden Markov Model, HMM)是一种最简单的动态贝叶斯网

马尔可夫链的相关基础知识在马尔可夫链基础中已提到，此处不再赘述

令${Xn,n=1,2,…}$是转移概率为$P{ij}$，初始状态概率为$\pi_i=P(X_1=i)$的马尔可夫链，记$S$是一个有限的信号集合

若马尔可夫链在每一时刻都根据当前状态$X_n$发射一个信号$Y_n\in S$，满足若$X_n=j$则发射信号$s_i$的概率为

$P{Yn=s_i|X_n=j,X_1,Y_1,…X{n-1},Y_{n-1}}=P{Y_n=s_i|X_n=j}$

即发射信号的概率条件独立于过往时刻的状态和信号，只依赖于当前状态，这样的模型称为隐马尔科夫模型

隐马尔可夫模型中，只有信号序列$Y_1,Y_2,…$能被观测到，状态序列$X_1,X_2,…$是不可被观测的

记HMM的状态转移矩阵为$A$，观测值生成概率矩阵为$B$，初始状态分布为$\pi$，其中$A{ij}$表示马尔可夫链由状态$i$转移到状态$j$的概率，$B{ij}$表示处于状态$i$的马尔可夫链发射观测信号$s_j$的概率

HMM主要关注三个问题

• 观测概率问题：假设$\lambda=[A,B,\pi]$已知，给定观测序列$Y=(Y_1,Y_2,…,Y_T)$，如何计算该观测序列的出现概率$P(Y|\lambda)$，也即如何评估模型与观测序列的相似程度
• 预测问题：假设$\lambda=[A,B,\pi]$已知，给定观测序列$Y=(Y_1,Y_2,…,Y_T)$，如何找到使得该观测序列的产生概率最大的状态转移序列$X=(X_1,X_2,…,X_T)$
• 学习问题：给定观测序列$Y=(Y_1,Y_2,…,Y_T)$，何如确定$\lambda=[A,B,\pi]$使得该观测序列的出现概率最大，也即如何训练模型使其能更好的描述观测数据

前后向算法

对于观测概率计算问题，若使用直接计算方法，则有

$$\begin{aligned} P(Y|\lambda)&=\sum_{X}P(X|\lambda)P(Y|X,\lambda)\\ &=\sum_X(\pi_{X_1}A_{X_1X_2}A_{X_2X_3}...A_{X_{T-1}X_T})\cdot(B_{X_1,Y_1}B_{X_2,Y_2}...B_{X_T,Y_T}) \end{aligned}$$

假设马尔可夫链的状态空间大小为$N$，则该方法复杂度为$O(TN^T)$，显然计算量过大

前后向算法使用了动态规划思想，首先说前向算法，令$\alpha_t(i)$表示$T$时刻时观测序列为$Y$，且状态为$i$的概率

$$\alpha_t(i)=P(Y_1,Y_2,...,Y_t,X_t=i|\lambda)$$

则有递推公式

$$\alpha_{t+1}(i)=B_{i, Y_{t+1}}\sum_{j}^N\alpha_t(j)A_{ji},\quad i=1,2,...,N$$

其含义不难理解，就是枚举上一时刻的可能状态进行概率求和。特别的，初值为

$$\alpha_1(i)=\pi_iB_{i, Y_1}$$

最终答案为

$$P(Y|\lambda)=\sum_{i=1}^N\alpha_T(i)$$

显然该算法的复杂度只有$O(TN^2)$

类似的，可以定义后向算法，令$\beta_t(i)$表示已知$t$状态为$X_t=i$时，$t+1$到$T$时刻的观测序列为$Y$的概率

$$\beta_t(i)=P(Y_{t+1},Y_{t+2},...,Y_T|X_t=i,\lambda)$$

则初值和递推公式为

$$\begin{aligned} &\beta_T(i)=1 \\ &\beta_t(i)=\sum_{j=1}^NA_{ij}B_{j,Y_{t+1}}\beta_{t+1}(j), \quad i=1,2,...,N \end{aligned}$$

最终答案为

$$P(Y|\lambda)=\sum_{i=1}^N\pi_iB_{i,Y_1}\beta_1(i)$$

结合前后向算法，对于任意时刻$t=1,2,…,T$有

$$P(Y|\lambda)=\sum_{i=1}^N\alpha_t(i)\beta_t(i)$$

前后向算法的思想很常见，例如CTC损失也用了类似的算法

Viterbi算法

对于前述预测问题，朴素的想法是直接选择每个时刻最有可能出现的状态，即时刻$t$处于状态$i$的概率为

$$P(X_t=i|Y,\lambda)=\frac{\alpha_t(i)\beta(i)}{P(Y|\lambda)}=\frac{\alpha_t(i)\beta(i)}{\sum_{j=1}^N \alpha_t(j)\beta(j)}$$

则时刻$t$最有可能的状态是$Xt^*=\arg\max{1\leq i\leq N}(P(X_t=i|Y,\lambda))$，从而得到状态序列$(X_1^,X_2^,…,X_T^*)$

该算法有一个明显缺点，即各个状态的预测是独立的，得到的状态序列不一定是整体上概率最大的

解决预测问题的另一个方法是维特比(Viterbi)算法，其基本思想是动态规划

令$f(t,i)$表示时刻 $t$ 处于状态 $i$ 的所有路径中概率最大的路径对应的概率，即

$$f(t,i)=\max P(X_t=i,X_{t-1},...,X_1,Y_t,...,Y_1|\lambda), \quad i=1,2,...N$$

令$g(t,i)$表示时刻 $t$ 处于状态 $i$ 的所有路径中概率最大的路径$t-1$时刻的状态，则有递推公式

$$f(t,i)=\max_{1\leq j\leq N}[f(t-1,j)A_{ji}]B_{i, Y_t}\\ g(t,i)=\arg \max_{1\leq j\leq N}[f(t-1,j)A_{ji}]$$

边界值为

$$f(1,i)=\pi_iB_{i,Y_1}\\ g(1,i)=0$$

则最终答案$(X_1^,X_2^,…,X_T^*)$为

$$X_T^*=\arg \max_{i} (f(T,i)) \\ X_t^*=g(t+1,X_{t+1}^*), \quad t=T-1,T-2,...,1$$

该算法复杂度为$O(NT)$

Baum-Welch算法

给定观测序列$Y=(Y_1,Y_2,…,Y_T)$，并将状态序列$X=(X_1,X_2,…,X_T)$视为隐变量，则隐马尔可夫模型参数$\lambda=[A,B,\pi]$的学习可以由EM算法实现，称其为Baum-Welch算法

首先易知

$$P(X,Y|\lambda)=\pi_{X_1}A_{X_1,X_2}...A_{X_{T-1},X_T}B_{X_1,Y_1}...B_{X_T,Y_T}$$

定义EM算法中E-step的Q函数

$$\begin{aligned} Q(\lambda|\lambda^{(i)}) &=\sum_Y\sum_X P(X|Y,\lambda^{(i)})\ln P(Y,X|\lambda)\\ &=\sum_Y\sum_X \frac{P(X,Y|\lambda^{(i)})}{P(Y|\lambda^{(i)})}\ln P(Y,X|\lambda) \end{aligned}$$

注意到$P(Y|\lambda^{(i)})$是与隐变量$X$和参数$\lambda$无关的常量，对于M-step求极值来说可以消去

$$\begin{aligned} Q(\lambda|\lambda^{(i)}) =&\sum_Y\sum_X P(X,Y|\lambda^{(i)})\ln P(Y,X|\lambda)\\ =&\sum_Y\sum_X P(X,Y|\lambda^{(i)})\ln \pi_{X_1}+\sum_Y\sum_X P(X,Y|\lambda^{(i)})\sum_{t=1}^{T-1}\ln A_{X_t,X_{t+1}}\\&+\sum_Y\sum_X P(X,Y|\lambda^{(i)})\sum_{t=1}^T\ln B_{X_t,Y_t} \end{aligned}$$

此时三个参数$\pi,A,B$分离到三个项中，可以分别对各项求极值

---

对于参数$\pi$

$$\begin{aligned} \pi^{(i+1)} =&\arg\max_\pi \sum_Y\sum_X P(X,Y|\lambda^{(i)})\ln \pi_{X_1}\\ =&\arg\max_\pi \sum_Y\sum_{X_1=1}^N\ln \pi_{X_1}\sum_{X/X_1}P(X_1,X/X_1,Y,\lambda^{(i)})\\ =&\arg\max_\pi \sum_Y\sum_{i=1}^N P(X_1=i,Y,\lambda^{(i)})\ln \pi_{i}\\ \end{aligned}$$

其中$X/X_1$表示状态序列$X$除去$X_1$后的子状态序列，推导中第2步到第3步实际上就是求边缘概率

此时加入限制条件$\sum_{i=1}^N\pi_i=1$，则根据拉格朗日乘子法设

$$f=\sum_Y\sum_{i=1}^N P(X_1=i,Y,\lambda^{(i)})\ln \pi_{i}+\gamma(\sum_{i=1}^N \pi_i-1)$$

求极值

$$\frac{\partial f}{\partial \pi_i}=\frac{\sum_Y P(X_1=i,Y|\lambda^{(i)})}{\pi_i}+\gamma=0$$

得

$$\begin{aligned} \sum_Y P(X_1=i,Y|\lambda^{(i)})+\gamma\pi_i&=0\\ \sum_Y \sum_{i=1}^N P(X_1=i,Y|\lambda^{(i)})+\gamma\sum_{i=1}^N\pi_i&=0\\ \sum_Y P(Y|\lambda^{(i)})+\gamma=0 \end{aligned}$$

将$\gamma$带回得

$$\pi^{(i+1)}=\pi_i=\frac{\sum_Y P(X_1=i,Y|\lambda^{(i)})}{\sum_Y P(Y|\lambda^{(i)})}$$

---

对于参数$A$，首先有

$$\sum_X P(X,Y|\lambda^{(k)})\sum_{t=1}^{T-1}\ln A_{X_t,X_{t+1}}=\sum_{i=1}^N\sum_{j=1}^N\sum_{t=1}^{T-1}\ln A_{ij}P(Y,X_t=i,X_{t+1}=j|\lambda^{k})$$

加入限制条件$\sum{j=1}^NA{ij}=1$，由拉格朗日乘子法设

$$f=\sum_{i=1}^N\sum_{j=1}^N\sum_{t=1}^{T-1}\ln A_{ij}P(Y,X_t=i,X_{t+1}=j|\lambda^{k})+\gamma(\sum_{j=1}^NA_{ij}-1)$$

求极值

$$\frac{\partial f}{\partial A_{ij}}=\sum_{t=1}^{T-1}\frac{P(Y,X_t=i,X_{t+1}=j|\lambda^{k})}{A_{ij}}+\gamma=0$$

得

$$\begin{aligned} \sum_{t=1}^{T-1}P(Y,X_t=i,X_{t+1}=j|\lambda^{k})+A_{ij}\gamma&=0\\ \sum_{j=1}^N\sum_{t=1}^{T-1}P(Y,X_t=i,X_{t+1}=j|\lambda^{k})+\sum_{j=1}^NA_{ij}\gamma&=0\\ \sum_{t=1}^{T-1}P(Y,X_t=i|\lambda^{k})+\gamma&=0 \end{aligned}$$

带回得

$$A_{ij}^{(k+1)}=A_{ij}=\frac{\sum_Y\sum_{t=1}^{T-1}P(Y,X_t=i,X_{t+1}=j|\lambda^{k})}{\sum_Y\sum_{t=1}^{T-1}P(Y,X_t=i|\lambda^{k})}$$

---

对于$B_{ij}$有

$$\sum_X P(X,Y|\lambda^{(k)})\sum_{t=1}^T\ln B_{X_t,Y_t}=\sum_{i=1}^N\sum_{t=1}^TP(Y,X_t=i|\lambda^{(k)})\ln B_{iY_t}$$

注意到当且仅当$j=Yt$时该式对$B{ij}$求偏导才不为0，因此需要引入指示函数

$$I(j,Y_t)=\begin{cases}1,\quad j=Y_t\\0, \quad j\neq Y_t\end{cases}$$

同理由$\sum{i=1}^MB{ij}=1$根据拉格朗日乘子法得

$$B_{ij}^{(k+1)}=\frac{\sum_{t=1}^TP(Y,X_t=i|\lambda^{(k)})I(j,Y_t)}{\sum_{t=1}^T P(X_t=i,Y|\lambda^{(k)})}$$