隐马尔可夫模型学习笔记（之二，学习算法）

隐马尔可夫模型的学习，根据训练数据是包括观测序列和状态序列还是只有观测序列，可以分别由监督学习与非监督学习实现。由于监督学习需要使用训练数据，而人工标注训练数据往往代价很高，有时就会利用非监督学习的方法，即Baum-Welch算法（也就是EM算法)。在介绍学习算法之前，先介绍一些概率和期望值的计算。这些计算会成为Baum-Welch算法公式的基础。

一些概率和期望值的计算

利用前向概率和后向概率，可以得到关于单个状态和两个状态概率的计算公式。

1. 给定模型$\lambda$和观测$O$，在时刻$t$处于状态$q_i$的概率。记为
   $${\gamma }_t(i)=P(i_t=q_i|O,\lambda )$$
   先分解为分数形式
   $${\gamma }_t(i)=\frac{P(i_t=q_i,O|\lambda )}{P(O|\lambda )}\text{\hspace{1em}(1)}$$
   根据前向概率的定义可以做以下变换
   $${\alpha }_t(i)=P(o_1,o_2...o_t,i_t=q_t|\lambda )=P(i_t=q_t|\lambda )P(o_1,o_2...o_t|i_t=q_t,\lambda )$$
   后向概率的定义如下
   $${\beta }_t(i)=P(o_{t+1},o_{t+2}...,o_T|i_t=q_t,\lambda )$$
   将这两者相乘得到

$${\alpha }_t(i)\ast {\beta }_t(i)=P(i_t=q_t,O|\lambda )\text{\hspace{1em}(2)}$$

以上结果从两者的定义上也很好理解。
对变量$i$在范围$i=1,2,...N$上求和
$$\sum _{i=1}^{N}P(i_t=q_t,O|\lambda )=P(O|\lambda )\text{\hspace{1em}(3)}$$
将式$(2),(3)$代入$(1)$可以得到
$${\gamma }_t(i)=\frac{{\alpha }_t(i)\ast {\beta }_t(i)}{{\sum }_{j=1}^N{\alpha }_t(j)\ast {\beta }_t(j)}\text{\hspace{1em}(4)}$$
3. 给定模型$\lambda$和观测$O$，在时刻$t$处于状态$q_i$且在时刻$t+1$处于状态$q_j$的概率。记为
$${\xi }_t(i,j)=P(i_t=q_i,i_{t+1}=q_j|O,\lambda )$$
通过前向后向概率计算：
$${\xi }_t(i)=\frac{P(i_t=q_i,i_{t+1}=q_j,O|\lambda )}{P(O|\lambda )}=\frac{P(i_t=q_i,i_{t+1}=q_j,O|\lambda )}{{\sum }_{i=1}^N{\sum }_{j=1}^{N}P(i_t=q_i,i_{t+1}=q_j,O|\lambda )}$$
分子可以用前向后向概率表示
$$P(i_t=q_i,i_{t+1}=q_j,O|\lambda )={\alpha }_t(i)a_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)$$
则${\xi }_t(i)$可以表示为
$${\xi }_t(i)=\frac{{\alpha }_t(i)a_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)}{{\sum }_{i=1}^N{\sum }_{j=1}^N{\alpha }_t(i)a_{ij}{b}_j(o_{t+1}){\beta }_{t+1}(j)}$$
4. 将${\gamma }_t(i)$和${\xi }_t(i,j)$对各个时刻求和，可以得到一些有用的期望值。
(1) 观测$O$下，状态$i$出现的期望值
$$\sum _{t=1}^T{\gamma }_t(i)$$
将每一个时刻下，出现状态$i$的概率相加
(2) 观测$O$下，由状态$i$转移的期望值
$$\sum _{t=1}^{T-1}{\gamma }_t(i)$$
能够从状态$i$转移的时刻是$1,2...T-1$，比上一个求和公式少了时刻$T$
(3) 观测$O$下，由状态$i$转移到状态$j$的期望值
$$\sum _{t=1}^{T-1}{\xi }_t(i,j)$$

Baum-Welch模型

参数估计公式

推导的过程，尤其是拉格朗日对偶，我暂时还不十分理解，先直接给出训练方法，公式和代码。Baum-Welch算法（Baum-Welch algorithm)，它是EM算法在隐马尔可夫模型学习过程中的具体实现，由Baum和Welch提出。

(1)初始化
对$n=0$，选取$a_{ij}^0，b_j(k{)}^0，{\pi }_i^0$，得到模型${\lambda }^0=(a_{ij}^0，b_j(k{)}^0，{\pi }_i^0)$

(2)递推。对$n=1,2,...$
$$a_{ij}^{n+1}=\frac{{\sum }_{t=1}^{T-1}{\xi }_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)$$
公式右端按照观测$O=(o_1,o_2,...o_T)$和模型${\lambda }^n=(a_{ij}^n，b_j(k{)}^n，{\pi }_i^n)$代入计算
这两个训练的公式还是比较好理解。$a_{ij}^{n+1}$是转移概率，分母代表当前观测$O$下，由状态$i$转移的期望值，而分子代表观测$O$下，由状态$i$转移到状态$j$的期望值。两者相除即为$a_{ij}^{n+1}$。
(3)终止，得到模型${\lambda }^{n+1}=(a_{ij}^{n+1}，b_j(k{)}^{n+1}，{\pi }_i^{n+1})$

###Baum-Welch算法的Python实现

def baum_welch_train(self, observations, criterion=0.05):
    n_states = self.A.shape[0]
    n_samples = len(observations)
 
    done = False
    while not done:
        # alpha_t(i) = P(O_1 O_2 ... O_t, q_t = S_i | hmm)
        # Initialize alpha
        alpha = self._forward(observations)
 
        # beta_t(i) = P(O_t+1 O_t+2 ... O_T | q_t = S_i , hmm)
        # Initialize beta
        beta = self._backward(observations)
 
        xi = np.zeros((n_states,n_states,n_samples-1))
        for t in range(n_samples-1):
            denom = np.dot(np.dot(alpha[:,t].T, self.A) * self.B[:,observations[t+1]].T, beta[:,t+1])
            for i in range(n_states):
                numer = alpha[i,t] * self.A[i,:] * self.B[:,observations[t+1]].T * beta[:,t+1].T
                xi[i,:,t] = numer / denom
 
        # gamma_t(i) = P(q_t = S_i | O, hmm)
        gamma = np.sum(xi,axis=1)
        # Need final gamma element for new B
        prod =  (alpha[:,n_samples-1] * beta[:,n_samples-1]).reshape((-1,1))
        gamma = np.hstack((gamma,  prod / np.sum(prod))) #append one more to gamma!!!
 
        newpi = gamma[:,0]
        newA = np.sum(xi,2) / np.sum(gamma[:,:-1],axis=1).reshape((-1,1))
        newB = np.copy(self.B)
 
        num_levels = self.B.shape[1]
        sumgamma = np.sum(gamma,axis=1)
        for lev in range(num_levels):
            mask = observations == lev
            newB[:,lev] = np.sum(gamma[:,mask],axis=1) / sumgamma
 
        if np.max(abs(self.pi - newpi)) < criterion and \
                        np.max(abs(self.A - newA)) < criterion and \
                        np.max(abs(self.B - newB)) < criterion:
            done = 1
 
        self.A[:],self.B[:],self.pi[:] = newA,newB,newpi