隐马尔可夫模型参数训练

上一篇我们讲了隐马尔可夫模型(HMM)的推导计算。对于一个单高斯隐马尔可夫模型而言,其参数集合包含 \(A=[a_{ij}]\) 转移概率矩阵,高斯分布均值向量 \(\boldsymbol{\mu}_{i,m}\),高斯分布的协方差 \(\boldsymbol{\Sigma}_{i,m}\)。 这一篇我们来讲讲如何从数据中训练得到这些参数。

我们假定我们的完整数据是 \(\boldsymbol{y}=\{\boldsymbol{o},\boldsymbol{h}\}\)\(\boldsymbol{o}\) 是观测数据,\(\boldsymbol{h}\) 是无法观测的隐马尔可夫状态。设我们需要训练的参数是 \(\theta\),那么我们就是需要最大化 \(\log p(\boldsymbol{y};\theta )\)。显然因为隐藏量 \(\boldsymbol{h}\) 的存在,我们不能直接去求取 \(\theta\) 。但假如我们已经有了关于 \(\theta\) 的一个比较好的估计,我们就可以计算在该估计值和观测数据条件下的 \(\log p(\boldsymbol{y};\theta)\) 如下( \(\theta_0\) 是上一个估计)

\[Q(\theta | \theta_0)=E_{\boldsymbol{h}|\boldsymbol{o}}[\log p(\boldsymbol{y};\theta)|\boldsymbol{o};\theta_0] = E[\log p(\boldsymbol{o};\boldsymbol{h})|\boldsymbol{o};\theta_0]\]

我们的隐藏量 \(\boldsymbol{h}\) 是离散的状态量,我们有

\[Q(\theta|\theta_0)=\sum_{\boldsymbol{h}} P(\boldsymbol{h}|\boldsymbol{o};\theta_0) \log p(\boldsymbol{y};\theta)\]

不断地迭代操作 E 步骤(计算期望)/ M 步骤(最大化期望),便可以找到 \(\theta\) 的局部最优解。

Baum-Whlch 算法

Baum-Whlch 算法是用于解决 GMM-HMM 参数学习的算法。根据上面的叙述,我们先写出辅助函数(或者说条件期望值)。

根据上一篇,我们有

\[Q(\theta|\theta_0)=E[\log P(\boldsymbol{o}^T,\boldsymbol{q}^T|\theta)| \boldsymbol{o}^T,\theta_0]\]

\(\boldsymbol{q}^T\) 是隐藏的状态量,无法观测。

E 步骤

先写出条件期望表达式

\[Q(\theta|\theta_0)=\sum_{\boldsymbol{q}_1^T}P(\boldsymbol{q}_1^T|\boldsymbol{o}_1^T,\theta_0) \log P(\boldsymbol{o}_1^T,\boldsymbol{q}_1^T|\theta)\]

状态 \(i\) 下的观测服从单高斯分布 \(P_{\boldsymbol{o}_t}(i) \sim N(\boldsymbol{\mu}_{i,m},\boldsymbol{\Sigma}_{i,m})\) 然后我们写出状态 \(i\) 的对数高斯 PDF

\[N_t(i)=\log P_{\boldsymbol{o}_t}(i) = -\frac D2\log (2 \pi) -\frac12 \log|\Sigma_i|-\frac12(\boldsymbol{o}_t-\boldsymbol{\mu}_i)^T \boldsymbol{\Sigma}^{-1}_i(\boldsymbol{o}_t-\boldsymbol{\mu}_i)\]

在继续之前我们引入一个概念叫克罗内克函数 \(\delta_{ij}\) 如下

\[\delta_{ij}=\begin{cases} 0 \ i\ne j \\ 1\ i=j \end{cases}\]

\(P(\boldsymbol{q}_1^T)=\prod^{T-1}_{t=1}a_{q_tq_{t+1}}\)\(P(\boldsymbol{o}_1^T,\boldsymbol{q}_1^T)=P(\boldsymbol{o}_1^T|\boldsymbol{q}_1^T)P(\boldsymbol{q}_1^T)\) 接下来我们有

\[\log P(\boldsymbol{o_1^T},\boldsymbol{q_1^T}|\theta)=\sum^T_{t=1}N_t{q_t}+\sum^{T-1}_{t=1}\log a_{q_tq_{t+1}}\]

所以我们可以重写 \(Q(\theta|\theta_0)\) 如下

\[\begin{aligned} Q(\theta | \theta_0) &= \sum_{\boldsymbol{q}_1^T}P(\boldsymbol{o_1^T},\boldsymbol{q_1^T}|\theta_0)\sum_{t=1}^TN_t(q_t)+\sum_{\boldsymbol{q}_q^T}P(\boldsymbol{o_1^T},\boldsymbol{q_1^T}|\theta_0)\sum_{t=1}^{T-1}\log a_{q_tq_{t+1}} \\ &=Q_1(\theta|\theta_0)+Q_2(\theta|\theta_0) \end{aligned}\]

其中

\[\begin{aligned} Q_1(\theta|\theta_0)&=\sum^N_{i=1}\{ \sum_{\boldsymbol{q}_1^T}P(\boldsymbol{q}_1^T|\boldsymbol{o}_1^T,\theta_0)\sum_{t=1}^TN_t(q_t) \}\delta_{q_t,i} \\ Q_2(\theta|\theta_0)&=\sum^N_{i=1}\sum^N_{j=1}\{\sum_{\boldsymbol{q}_1^T}P(\boldsymbol{q_1^T}|\boldsymbol{o_1^T},\theta_0)\sum^{T-1}_{t=1}\log a_{q_tq_{t+1}} \}\delta_{q_t,i} \delta_{q_{t+1},j} \end{aligned}\]

\(\sum_{\boldsymbol{q}_1^T}P(\boldsymbol{q}_1^T|\boldsymbol{o}_1^T,\theta_0)\delta_{q_t,i}=P(\boldsymbol{q}_1^T=i|\boldsymbol{q}_1^T,\theta_0)\) 带入\(Q_1(\theta|\theta_0)\)我们有

\[Q_1(\theta|\theta_0)=\sum_{i=1}^N\sum_{t=1}^TP(\boldsymbol{q}_1^T=i|\boldsymbol{q}_1^T,\theta_0)N_t(i)\]

类似的我们有

\[Q_2(\theta|\theta_0)=\sum_{i=1}^N\sum_{j=1}^N\sum_{t=1}^{T-1}P(q_t=i,q_{t+1}=j|\boldsymbol{o}_1^T,\theta_0)\log a_{ij}\]

在给定观测序列 \(\boldsymbol{q}_1^T\) 和参数 \(\theta_0\) 的情况下,在时间 \(t\) 状态是 \(i\),在时间 \(t+1\) 状态是 \(j\) 的概率是

\[\begin{aligned} \xi_t(i,j)&=P(q_t=i,q_{t+1}=j|\boldsymbol{o}_1^T,\theta_0)=\frac{P(q_t=i,q_{t+1}=j,\boldsymbol{o}_1^T|\theta_0)}{P(\boldsymbol{o}_1^T|\theta_0)} \\ &=\frac{\alpha_t(i)\beta_{t+1}(j)a_{ij}\exp(N_{t+1}(j))}{P(\boldsymbol{o}_1^T|\theta_0)} \end{aligned}\]

这里文字解释一下(我这里推导卡了好久,MMP)推导过程。第一个等号用的贝叶斯条件概率 \(P(A|B)=\frac{P(AB)}{P(B)}\) 就可以。第二个等号主要是 \(P(q_t=i,q_{t+1}=j,\boldsymbol{o}_1^T|\theta_0)=\alpha_t(i)\beta_{t+1}(j)a_{ij}\exp(N_{t+1}(j))\) 解释一下:\(\alpha_t(i)\) 指的是时刻 \(t\) 的状态为 \(i\), 时刻 1 到 \(t\) 的观测概率,\(\beta_{t+1}(j)\) 则是在时刻 \(t+1\) 的状态为 \(j\),时刻 \(t+2\)\(T\) 的观测概率,\(a_{ij}\exp(N_{t+1}(j))\) 则是从状态 \(i\) 变化到 \(j\),时刻 \(t+1\) 的观测概率。不过要注意这里 \(\xi_T(i,j)\) 没有定义。

在给定观测序列 \(\boldsymbol{q}_1^T\) 和参数 \(\theta_0\) 的情况下,在时间 \(t\) 状态是 \(i\) 的概率是

\[\gamma_t(i)=\sum_{j=1}^N\xi_t(i,j)\]

而对于 \(t=T\) 的情况,我们可以用它的特定定义得到

\[\gamma_T(i)=P(q_T=i|\boldsymbol{o}_1^T,\theta_0)=\frac{P(q_T=i,\boldsymbol{o}_1^T|\theta_0)}{P(\boldsymbol{o}_1^T|\theta_0)}=\frac{\alpha_T(i)}{P(\boldsymbol{o}_1^T|\theta_0)}\]

M 步骤

最大化 \(Q_2\)

我们要最大化 \(Q_2\) ,同时需要满足条件 \(\sum^N_{j=1}a_{ij}=1\)。可以应用拉格朗日乘子法使得重估公式变为

\[\hat{a}_{ij}=\frac{\sum^{T-1}_{t=1}\xi_t(i,j)}{\sum^{T-1}_{t=1}\gamma_t(i)}\]

什么鬼!让我们还是先来复习一下什么叫拉格朗日乘子法吧。设函数 \(f(x,y)\)\(A\) 点处有极值 \(\kappa\),且在 \(A\) 点的邻域内连续,则在 \(A\) 点有 \(f(x,y)=\kappa\)。另有常值函数 \(g(x,y)=c\),则两个函数在 \(A\) 点的全微分如下

\[\mathrm{d}f=\frac{\partial f}{\partial x}\mathrm{d}x+\frac{\partial f}{\partial y}\mathrm{d}y\]

\[\mathrm{d}g=\frac{\partial g}{\partial x}\mathrm{d}x+\frac{\partial g}{\partial y}\mathrm{d}y\]

由于 \(\mathrm{d}x\)\(\mathrm{d}y\) 是任意无穷小量,故该线性方程组的系数成比例

\[\frac{\frac{\partial f}{\partial x}}{\frac{\partial g}{\partial x}}=\frac{\frac{\partial f}{\partial y}}{\frac{\partial g}{\partial y}}=-\lambda\]

\[\frac{\partial f}{\partial x}+\lambda\frac{\partial g}{\partial x}=0\]

\[\frac{\partial f}{\partial y}+\lambda\frac{\partial g}{\partial y}=0\]

将上二式分别乘以 \(\mathrm{d}x\)\(\mathrm{d}y\),再相加积分,得到拉格朗日函数如下,求原函数极值的问题就转化为求该函数极值的问题。

\[L(x,y,\lambda)=f(x,y)+\lambda (g(x,y)-c) \]

回到刚刚的求 \(Q_2\) 极值的问题,我们应用拉格朗日乘子法如下

\[\sum_{i=1}^N\sum_{j=1}^N\sum_{t=1}^{T-1}\xi_t(i,j)\log a_{ij}+\lambda(\sum_j a_{ij}-1)\]

\(a_{ij}\) 求偏导,然后令其为 0

\[\begin{aligned} \frac{\sum_{t=1}^{T-1}\xi_t(i,j)}{a_{ij}} + \lambda &= 0 \\ \sum_{t=1}^{T-1}\xi_t(i,j)+{a_{ij}} \lambda &= 0 \\ \sum_{j=1}^N( \sum_{t=1}^{T-1}\xi_t(i,j)+{a_{ij}} \lambda) &= 0 \\ \sum_{t=1}^{T-1} \gamma_t(i) + \lambda &= 0 \end{aligned}\]

再代入原来的式子,就有

\[\hat{a}_{ij} = \frac{\sum_{t=1}^{T-1} \xi_t(i,j)}{\sum_{t=1}^{T-1}\gamma_t(i)}\]

最大化 \(Q_1\)

去掉 \(Q_1\) 中与优化无关的式子,得到下面的目标函数

\[Q_1(\boldsymbol{\mu}_i,\boldsymbol{\Sigma}_i)=\sum_{i=1}^N\sum_{t=1}^T\gamma_t(i)(\boldsymbol{o}_t-\boldsymbol{\mu}_i)^T \boldsymbol{\Sigma_t^{-1}}(\boldsymbol{o}_t-\boldsymbol{\mu}_i)-\frac12\log |\boldsymbol{\Sigma}_i|\]

同样根据 \(\frac{\partial Q_1}{\partial \boldsymbol{\Sigma}_i}=0\) 来得到重估计公式如下(懒得去推导了。。。抄书抄书)

\[\hat{\boldsymbol{\Sigma}}_i=\frac{\sum_{t=1}^T \gamma_t(i)(\boldsymbol{o}_t-\hat{\boldsymbol{\mu}}_i)(\boldsymbol{o}_t-\hat{\boldsymbol{\mu}}_i)^T}{\sum_{t=1}^T \gamma_t(i)}\]

\[\hat{\boldsymbol{\mu}}_i=\frac{\sum_{t=1}^T\gamma_t(i)\boldsymbol{o}_t}{\sum_{t=1}^T \gamma_t(i)}\]