隐马尔可夫模型参数训练
上一篇我们讲了隐马尔可夫模型（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)}\]