EM & HMM

Algorithm

1、EM（Expectation Maximization）

1.1 EM算法简介

• EM算法是一种迭代算法，用于含有隐变量（hidden variable）的概率模型参数的极大似然估计（Maximum Likelihood Estimation)）。
• EM算法分为 E步（求期望）和 M步（求极大）。
• EM算法的基本思想是：
  
  1. 若参数$\theta$已知，则可根据训练数据推断出最优隐变量$Z$的值（E步）
  2. 若隐变量$Z$已知，则可方便地对参数$\theta$进行MLE；（M步）
  3. 说人话：EM算法首先固定其中的第一个参数，然后使用MLE计算第二个变量值。接着通过固定第二个变量，再使用MLE估测第一个变量值，依次迭代，直至收敛到局部最优解。
• EM算法的伟大之处：
  突破了MLE的限制，即当存在hidden variable时，也可以一定能找到最优解（收敛性），大概率是局部最优。

1.2 EM 算法的理论推导

1.2.1 三硬币模型

假设有三枚硬币A、B、C，每个硬币正面出现的概率$\pi、p、q$。进行如下的掷硬币实验：先掷硬币A，正面向上选B，反面选C；然后掷选择的硬币，正面记1，反面记0。独立的进行10次实验，结果如下：1，1，0，1，0，0，1，0，1，1。假设只能观察最终的结果(0 or 1)，而不能观测掷硬币的过程(不知道选的是B or C)，问如何估计三硬币的正面出现的概率？

• 首先针对某个输出值$y$，它在参数$\theta (\theta=(\pi, p, q))$下的概率分布为：
  $\qquad \begin{align*} P(y|\theta ) &= \sum_{z}P(y,z|\theta) \\ &= \sum_{z}P(z|\theta)P(y|z, \theta) \\ &= \prod _{j=1} ^{n} \pi p^y_j (1-p)^{1-y_j} + (1-\pi) q^y_j (1-q)^{1-y_j} \end{align*}$

• 容易得到，对于观测数据$Y=(y_1, y_2, \cdot\cdot\cdot, y_n)^T$的似然函数为：
  $\qquad \begin{align*} P(Y|\theta ) &= \sum_{z}P(Y,z|\theta) \\ &= \sum_{z}P(z|\theta)P(Y|z, \theta) \\ &= \prod _{j=1} ^{n} \pi p^y_j (1-p)^{1-y_j} + (1-\pi) q^y_j (1-q)^{1-y_j} \end{align*}$
  这里有一个trick： $\pi \in \{0,1\}$

• 因此本题的目标是求解参数$\theta$的极大似然估计，$\hat{\theta} = \underset{\theta }{argmax}logP(Y|\theta)$。

• 由于目标函数里有 hidden variable，无法用常规的 MLE 来求解，故这个问题需要用 EM 算法来求解。

  • E步：当第 t+1 步，隐变量的条件概率分布，也就是每一个$P(z|y_j, \theta)$，$\mu^{i+1}$表示在已知的模型参数$\theta^i$下观测数据$y_j$来自掷硬币B的概率，相应的来自掷C的概率就是$1-\mu^{i+1}$。
    $\qquad \mu ^{i+1} = \frac {\pi^i ({p^i})^{y_j}(1-p^i)^{1-y_j}} {\pi^i ({p^i})^{y_j}(1-p^i)^{1-y_j} + (1-\pi^i) ({q^i})^{y_j} (1-q^i)^{1-y_j}}$
  • M步：根据Q函数的定义，可得其表达式为：
    $\qquad % <![CDATA[ \begin{equation} \begin{aligned} Q(\theta,\theta^{(i)}) &=E_z[logP(Y,Z|\theta)|Y,\theta^{(i)}]\\ &=\Sigma_zlogP(Y,Z|\theta)P(Z|Y,\theta^{(i)}) \end{aligned} \end{equation} %]]>$
  • 通过推导，可以得出Q函数最后的形式如下：（推导过程见夹在《统计学习方法》中的A4纸）
    $\qquad Q(\theta, \theta^i)=\sum_{j=1}^N \mu_j^{t+1} \cdot log[\pi p^{y_j}(1-p)^{1-y_j}] + (1-\mu_j^{t+1}) \cdot log[(1-\pi) q^{y_j}(1-q)^{1-y_j}]$
  • 根据Q函数求解各个参数的解析解；（详见夹在《统计学习方法》中的A4纸）

• 小结
  

  1. EM算法将不完全的数据补全成完全数据，其中E步是将未观测的数据的所有补全可能都计算出对应的概率值，从而对所有可能补全计算出它们的期望值，作为下一步的未观测数据。
     
  2. 上述中求期望值的原因：
     既然$Z$是 hidden variable，而下一次的迭代又必须已知其值才能进行，故使用期望值更为合理。
     
  3. Q函数想要表达的逻辑和思想：
     
     • 在当前迭代中要找出一个新的$\theta$，使得以新的$\theta$为参数的某分布优于上一轮迭代的结果。
     • 但因为存在 hidden variable，且又找不到更好的办法知道 hidden variable 的分布，故只能选用期望。

1.2.2 EM算法的基本步骤

需要注意的点：
1. EM算法参数的初值可以任意选择，但EM算法对初值是敏感的；
2. 常用的解决方法是选取几个不同的初值进行迭代，然后对得到的各个估计值进行比较，从中选择最好的。

1.2.3 EM算法一般性推导

• EM算法的目标是极大化观测数据（不完全数据）$Y$关于参数$\theta$的对数似然函数，即极大化下式
  $L（\theta）=\log \sum_zP(Y \mid Z,\theta)P(Z \mid \theta)$
• 我们希望新的估计值$\theta$能使$L(\theta)$增加，即$L(\theta)>L(\hat{\theta})$，为此考虑二者的差
  
  其中第（3）步计算利用了Jensen不等式
• 对于上式，令
  
  EM算法通过不断通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。

1.3 EM 算法的 Python 实现

#!/usr/bin/env python
#coding:utf-8
"""
从工程上来看, 算法写的很一般。
但可以帮助自己理解如何用计算机实现一个复杂的机器学习算法。
代码非原创，来自github，地址如下：
https://github.com/jindc/em
"""
import sys
import math
class em():
    def run(self,data,loopcnt=900000):
        print data
        pi=0.4
        p=0.6
        q=0.7
        datanum=float(len(data))
        def cpi(j):
            tmp1 = pi*math.pow(p,data[j])*math.pow(1-p,1-data[j])
            tmp2 = (1-pi)*math.pow(q,data[j])*math.pow(1-q,1-data[j])
            return tmp1/(tmp1+tmp2)
        for i in range(loopcnt):
            pi = 1/datanum *\
                 sum([cpi(j)for j in range(int(datanum))])
            p = sum([cpi(j)*data[j] for j in range(int(datanum))]) /(sum([cpi(j)for j in range(int(datanum))]))
            q = sum([(1-cpi(j))*data[j] for j in range(int(datanum))] )/(sum([(1-cpi(j)) for j in range(int(datanum))]))
            if str(pi)[0:6]=='0.4064'or str(pi)[0:6]=='0.4063':
                if str(p)[0:6]=='0.5368' or str(p)[0:6]=='0.5367':
                    if str(q)[0:6]=='0.6432' or str(q)[0:6]=='0.6431':
                        break
            # print i,pi ,p, q
        self.pi =pi
        self.p = p
        self.q = q
if __name__=='__main__':
    data=[1,1,0,1,0,0,1,0,1,1]
    emi = em()
    emi.run(data)
    # print '---------'
    print emi.pi,emi.p,emi.q

2、HMM（Hidden Marcov Model）

2.1 HMM 的基本概念

• 隐马尔可夫模型的是处理序列问题的统计学模型，描述的过程为：由隐马尔科夫链随机生成不可观测的状态随机序列，然后各个状态分别生成一个观测，从而产生观测随机序列。
• 在这个过程中，不可观测的序列称为状态序列(state sequence), 由此产生的序列称为观测序列(observation sequence).
• 该过程可通过下图描述
  
  上图中，$X_1,X_2,…X_T$是隐含序列，而$O_1, O_2,..O_T$是观察序列。

• 隐马尔可夫模型由三个概率确定：

  1. 初始概率分布（向量），即初始的隐含状态的概率分布，记为$\pi$。
  2. 状态转移概率分布（矩阵），即隐含状态间的转移概率分布, 记为$A$。
  3. 观测概率分布（矩阵），即由隐含状态生成观测状态的概率分布, 记为$B$。
     
     • 状态转移概率分布$A$和初初始状态概率向量$\pi$，确定了隐藏的马尔科夫链，生成不可观测的状态序列。
     • 观测概率分布$B$确定了如何从状态生成观测，与状态序列综合确定了如何产生观测序列。
     • $\pi$和$A$决定状态序列，$B$决定观测序列。
     • 以上的三个概率分布可以说就是隐马尔可夫模型的参数，而根据这三个概率，根据这三个概率，能够确定一个隐马尔可夫模型，$\lambda = (A, B, \pi)$

• HMM的两个假设：

  1. 齐次马尔可夫性假设。任意时刻 t 的状态, 只依赖于其前一刻的状态, 与其他时刻的状态及观测无关, 也与时刻 t 无关。
     $\qquad P(i_t|i_{t-1},o_{t-1},i_{t-2},o_{t-2},...,i_{1},o_{1}) = P(i_t|i_{t-1})$
  2. 观测独立性假设。时刻的观测结果只与时刻的状态结果有关，和其他无关。
     $\qquad P(o_t|i_t,i_{t-1},o_{t-1},...,i_{1},o_{1}) = P(o_t|i_t)$

2.2 HMM的三个基本问题

1. 概率计算问题。给定模型$\lambda = (A, B, \pi)$和观测序列 $O$，计算在模型$\lambda$下观测序列出现的最大概率$P(O|\lambda)$。
2. 学习问题。给定观测序列$O$，估计模型的参数$\lambda$, 使得在该参数下观测序列出现的概率最大，即$P(O|\lambda)$最大。【易见：问题1是问题2的基础】
3. 预测问题。给定模型$\lambda = (A, B, \pi)$和观测序列 $O$，计算最有可能产生这个观测序列的隐含序列$X$, 即使得概率 $P(X|O, \lambda)$最大的隐含序列$X$。

2.3 概率计算问题

• 前向算法

  1. 前向概率：给定隐马尔可夫模型$\lambda$，定义为到时刻$t$为止观测序列为 $o_1,o_2,o_3…o_t$，且状态为$q_i$的概率，记为：
     $\qquad \alpha_t(i) = P(o_1,o_2,o_3…o_t,x_t = q_i|\lambda)$
  2. 算法流程：
     输入：隐马尔可夫模型$\lambda$，观测序列$O$。
     输出：观测序列概率$P(O|\lambda)$
     （1）初值：$\alpha_1(i)=\pi_i b_i(o_1),\quad i=1,2,\cdots,N$
     （2）递推：对于$t=1,2,\cdots,T-1$，
     $\qquad \alpha_{t+1}(i)=\big[\sum\limits_{j=1}^N\alpha_t(j)a_{ji}\big]b_i(o_{t+1})\quad i=1,2,\cdots,N$
     （3）终止
     $\qquad P(O\mid \lambda)=\sum\limits_{i=1}^N\alpha_T(i)$
  3. 算法理解：
     
     • 为什么是累加求和？
       
       1. 因为可以产生观测$O$的状态有很多种，所以求出它们的边缘概率才符合客观事实。
       2. $\alpha_T(i)$代表着某一种状态序列产生观测序列$O$的概率，而状态取值有$N$种，故应该求和。
     • 一张图搞定递推式子的本质理解
       

• 后向算法

  1. 后向概率：给定隐马尔可夫模型$\lambda$，定义在时刻$t$状态为$q_i$的条件下，从$t+1$到$T$的部分观测序列为$o_{t+1},o_{t+2},\dots,o_T$的概率，记为：
     $\qquad \beta_t(i) = P(o_{t+1}, o_{t+2},….o_T, x_t = q_i|\lambda)$
  2. 算法流程：
     输入：隐马尔可夫模型$\lambda$，观测序列$O$。
     输出：观测序列概率$P(O|\lambda)$
     （1）初值：$\beta_T(i)=1,\quad i=1,2,\cdots,N$
     （2）递推：对于$t=T-1,T-2,\cdots,1$
     $\qquad \beta_t(i)=\sum\limits_{j=1}^Na_{ij}b_j(o_{t+1})\beta_{t+1}(j),\quad i=1,2,\cdots,N$
     （3）终止
     $\qquad P(O\mid \lambda)=\sum_{i=1}^N\pi_i b_i(o_1)\beta_1(i)$
  3. 算法理解：
     
     • 一张图搞定递推式子的本质理解
       
     • 根据上述两幅图来写各自的递推公式会让逻辑更清晰，堪称 $easy$ $job$。
     • 能够产生观测序列$O$的状态有很多种，故一定是根据状态变量求累加（边缘概率）。

• 小结

  1. 利用前向概率和后向概率的定义，可以将观测序列概率$P(O|\lambda)$统一写成：
     $\qquad P(O|\lambda)=\sum\limits_{i=1}^N\sum\limits_{j=1}^N\alpha_t(i)a_{ij}\beta_{t+1}(j) · b_j(o_{t+1}),\quad i=1,2,\cdots,T-1$
  2. 公式理解：
     
     • 点乘前面的部分，保证了$t$时刻的状态为$q_j$；
     • $\alpha_t$和$\beta_{t+1}$保证了除$t$时刻外剩余观测序列的正确性；

• 一些概率与期望值的计算

  1. 给定模型$\lambda$和观测$O$，在时刻$t$处于状态$q_i$的概率：
     $\qquad \begin{equation} \begin{aligned} \gamma_t(i)=\frac{\alpha_t(i)\beta_t(i)}{\sum\limits_{j=1}^N\alpha_t(j)\beta_t(j)}\ \end{aligned} \end{equation}$
  2. 给定模型$\lambda$和观测$O$，在时刻$t$处于状态$q_i$且时刻$t+1$处于状态$q_j$的概率
     $\qquad \begin{equation} \begin{aligned} \zeta_t(i,j) &=\frac{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{\sum\limits_{i=1}^N\sum\limits_{j=1}^N \alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)} \end{aligned} \end{equation}$
  3. 将$\gamma_t(i)$和$\zeta_t(i,j)$对各个时刻$t$求和，可以得到一些有用的期望值：
     
     • 在观测$O$下状态$i$出现的期望值：
       $\qquad\qquad\qquad \sum\limits_{t=1}^{T}\gamma_t(i)$
     • 在观测$O$下由状态$i$转移的期望值：
       $\qquad\qquad\qquad \sum\limits_{t=1}^{T-1}\gamma_t(i)$
     • 在观测$O$下由状态$i$转移到状态$j$的期望值：
       $\qquad\qquad\qquad \sum\limits_{t=1}^{T-1}\zeta_t(i,j)$

2.4 学习问题

1. 学习问题是根据观测序列来推导模型参数，这一类的问题对应到概率论中的极大似估计问题。但是这里是有隐含变量的极大似然估计，因此无法使用 MLE 来求解，而要通过 EM 算法来求解。
2. Baum-Welch 是 EM 算法在 HMM 求解学习问题的具体应用，具体的求解公式如下：
   $\qquad\qquad\qquad a_{ij} = \frac{\sum\limits_{t=1}^T \xi_t(i,j)}{\sum\limits_{t=1}^T \gamma_t(i)}$
   $\qquad\qquad\qquad b_j(k) = \frac{\sum\limits_{t=1,o_t = v_k}^T \gamma_t(j)}{\sum\limits_{t=1}^T \gamma_t(j)}$
   $\qquad\qquad\qquad \pi_i = \gamma_1(i)$

2.5 预测问题

1. 维特比（Viterbi Algorithm）是用动态规划（dynamtic programming）求解 HMM 的预测问题，即求概率最大路径（最优路径），这条路径对应着一个状态序列。
2. 由于最优路径的子路径也是最优路径，故可以用迭代思想求解出最优路径。
3. Viterbi Algo 的算法流程：
   
4. 对 Viterbi Algo 的理解
   
   • 它是基于迭代思想找到最优路径；
   • 首先，$\zeta_t(i)$控制着“个体”之间的“最优解”，即$t$时刻$\zeta_t(i)$最大的$i$就是当前时刻最棒的“个体”；
   • 然后，基于$\psi_t(i)$从后往前推，找出最优路径，即在考虑前一个时刻的情况（$a_{ji}$）后，“最优个体”是否会发生变化。
   • 总结起来就是，$\zeta_t(i)$定最优终点，$\psi_t(i)$回溯出最优路径。