sn06.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 隐式马尔可夫模型基础\n",
    "\n",
    "## 摘要\n",
    "\n",
    "我们如何将机器学习应用在基于时间序列的观测数据上？举个例子，比如我们可能想要从某人演讲的音频记录中得到他使用的词语序列；或者根据词性标注对一系列单词做注解等等。在本文中，我们将对马尔可夫模型（在时间序列上的状态的形式化推理）以及隐式马尔可夫模型（从观测数据的序列中找回原来的状态序列）的概念从数学方面做一个全面的介绍。本文的最后一节包含了本文中某些观点的来源并提供深入学习的材料。\n",
    "\n",
    "## 1. 马尔可夫模型\n",
    "\n",
    "对于给定的状态集$S=\\left\\{s_1,s_2,\\cdots,s_{\\lvert S\\rvert}\\right\\}$，我们可以沿时间序列观测到一些列数据$\\vec z\\in S^T$。举个例子，一个从天气系统中得到的状态集$S=\\left\\{sun,cloud,rain\\right\\}$，集合的大小为$\\lvert S\\rvert=3$，我们也在天气系统中得到一组观测数据$\\left\\{z_1=s_{sun},z_2=s_{cloud},z_3=s_{cloud},z_4=s_{rain},z_5=s_{cloud}\\right\\}$，观测数据集大小为$T=5$。\n",
    "\n",
    "天气例子的的观测状态向我们展示了一个在时间上的随机过程的输出。如果不做进一步假设，$t$时刻的状态$s_j$可能是一个关于任意多个变量的函数，其中包括了从时刻$1$到$t-1$的所有状态、以及一些我们甚至都没有建模的变量。不过，我们可以做出两个能够使时间序列更已处理的**马尔可夫假设（Markov Assumptions）**。\n",
    "\n",
    "**限制时域假设（Limited Horizon Assumption）**指的是到达$t$时刻状态的概率只与$t-1$时刻的状态有关。这条假设的直观理解就是$t$时刻的状态“足以概括”之前所有状态，并可以合理的预测未来的状态，形式化的表示为：\n",
    "$$\n",
    "P\\left(z_t\\mid z_{t-1},z_{t-2},\\cdots,z_1\\right)=P\\left(z_t\\mid z_{t-1}\\right)\n",
    "$$\n",
    "\n",
    "**平稳过程假设（Stationary Process Assumption）**指的是在当前状态下，下个状态的条件分布不随时间变化而改变，形式化的表示为：\n",
    "$$\n",
    "P\\left(z_t\\mid z_{t-1}\\right)=P\\left(z_2\\mid z_1\\right);\\ t\\in2\\cdots T\n",
    "$$\n",
    "\n",
    "按照习惯，我们也会假设初始状态和初始观测数据满足$z_0\\equiv s_0$，$s_0$表示在$0$时刻的状态的初始概率分布。这个习惯可以使我们将“观测第一个真实状态$z_1$的先验概率”的观点通过$P\\left(z_1\\mid z_0\\right)$表达出来。因为我们定义对于任意状态序列都有$z_0=s_0$，所以$P\\left(z_t\\mid z_{t-1},z_{t-2},\\cdots,z_1\\right)=P\\left(z_t\\mid z_{t-1},z_{t-2},\\cdots,z_1,z_0\\right)$。（其他HMM的表示方法有时通过向量$\\pi\\in\\mathbb R^{\\lvert S\\rvert}$表达其先验的观点。）\n",
    "\n",
    "我们定义一个状态转换矩阵$\\displaystyle A\\in\\mathbb R^{(\\lvert S\\rvert+1)\\times(\\lvert S\\rvert+1)}$来参数化这些转换，$A_{ij}$就表示在任意时刻$t$下从状态$i$转换到状态$j$的概率。（可以参考[线性代数中马尔科夫矩阵的例子](http://nbviewer.jupyter.org/github/zlotus/notes-linear-algebra/blob/master/chapter24.ipynb)）在天气系统的例子中，就可以定义如下的状态转换矩阵：\n",
    "$$\n",
    "A=\\begin{matrix}\n",
    "&s_0&s_{sun}&s_{cloud}&s_{rain}\\\\\n",
    "s_0&0&.33&.33&.33\\\\\n",
    "s_{cloud}&0&.2&.6&.2\\\\\n",
    "s_{rain}&0&.1&.2&.7\n",
    "\\end{matrix}\n",
    "$$\n",
    "\n",
    "容易观察到这些数字（虽然是编造的）表示的直观的天气的自相关性：如果当前是晴天，则天气倾向于保持晴朗；如果是多云则倾向于保持多云，等等。这种模式经常在马尔可夫模型中出现，体现在状态转移矩阵中就是对角线元素远大于其他元素。在本例中，初始状态$s_0$是转换到这三个天气状态上的均匀概率分布。\n",
    "\n",
    "### 1.1 马尔可夫模型的两个问题\n",
    "\n",
    "结合马尔可夫假设和状态转换参数$A$，我们可以回答在马尔科夫链中关于状态序列的两个基本问题：\n",
    "* 状态$\\vec z$对应的一个特定序列的概率是什么？\n",
    "* 对于观测到的序列$\\vec z$，如何估计模型中的参数$A$的才能最大化其似然性？\n",
    "\n",
    "#### 1.1.1 状态序列的概率\n",
    "\n",
    "我们可以利用概率论中的[链式法则](https://en.wikipedia.org/wiki/Chain_rule_(probability))状态$\\vec z$对应的一个特定序列的概率：\n",
    "$$\n",
    "\\begin{align}\n",
    "P\\left(\\vec z\\right)&=P\\left(z_t,z_{t-1},\\cdots,z_1;A\\right)\\\\\n",
    "&=P\\left(z_t,z_{t-1},\\cdots,z_1,z_0;A\\right)\\\\\n",
    "&=P\\left(z_t\\mid z_{t-1},z_{t-2},\\cdots,z_1;A\\right)P\\left(z_{t-1}\\mid z_{t-2},\\cdots,z_1;A\\right)\\cdots P(z_1\\mid z_0;A)\\\\\n",
    "&=P\\left(z_t\\mid z_{t-1};A\\right)P\\left(z_{t-1}\\mid z_{t-2};A\\right)\\cdots P\\left(z_2\\mid z_1;A\\right)P\\left(z_1\\mid z_0;A\\right)\\\\\n",
    "&=\\prod_{t=1}^TP\\left(z_{t}\\mid z_{t-1};A\\right)\\\\\n",
    "&=\\prod_{t=1}^TA_{z_{t-1}\\ z_t}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "式中的第二步，我们根据前面关于$z_0$的定义，将$z_0$引入联合分布。第三步是由对联合分布应用链式法则，或重复迭代贝叶斯发展推导出的。第四步使用了马尔可夫模型的限制时域假设。而最后一步是将这些乘积项用其对应在状态转换矩阵$A$中的元素表示出来。\n",
    "\n",
    "我们来计算一下天气系统例子中的序列：如果要计算序列$P\\left(z_1=s_{sun},z_2=s_{cloud},z_3=s_{rain},z_4=s_{rain},z_5=s_{cloud}\\right)$的概率，则可以分解为$P\\left(s_{sun}\\mid s_0\\right)P\\left(s_{cloud}\\mid s_{sun}\\right)P\\left(s_{rain}\\mid s_{cloud}\\right)P\\left(s_{rain}\\mid s_{rain}\\right)P\\left(s_{cloud}\\mid s_{rain}\\right)=.33\\times.1\\times.2\\times.7\\times.2$。\n",
    "\n",
    "#### 1.1.2 参数的最大似然估计\n",
    "\n",
    "从学习的角度，我们想要知道能够“最大化观测到的序列$\\vec z$的对数似然性”的参数$A$的取值。这相当于比较“从晴到阴”与“从晴到晴”等状态的似然性，以使得某个观测序列更有可能发生。我们定义马尔可夫模型的对数似然函数：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "l(A)&=\\log P\\left(\\vec z;A\\right)\\\\\n",
    "&=\\log\\prod_{t=1}^TA_{z_{t-1}\\ z_t}\\\\\n",
    "&=\\sum_{t=1}^T\\log A_{z_{t-1}\\ z_t}\\\\\n",
    "&=\\sum_{i=1}^{\\lvert S\\rvert}\\sum_{j=1}^{\\lvert S\\rvert}\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_{t}=s_j\\right\\}\\log A_{ij}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "我们在式中的最后一步使用了指示函数（其参数成立时函数值取$1$，否则取$0$）来从状态转移矩阵$A$中选出从时间序列中观测到的每一步状态转移对应的概率。在解这个优化问题时要注意，解出的参数$A$必须满足一个状态转移矩阵的要求——于每个状态$i$对应的分布，其概率之和必须为$1$；矩阵的所有元素均为非负数。我们可以使用拉格朗日乘数法来求解这个优化问题：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\max_A&\\ l(A)\\\\\n",
    "\\mathrm{s.t.}&\\ \\sum_{j=1}^{\\lvert S\\rvert}A_{ij}=1,\\ i=1,\\cdots,\\lvert S\\rvert\\\\\n",
    "&\\ A_{ij}\\geq0,\\ i,j=1,\\cdots,\\lvert S\\rvert\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "这个带有约束条件的优化问题可以使用拉格朗日乘数法直接求出其解析解。我们会将等式约束条件（即第一个约束条件）引入拉格朗日算子，而不等式约束条件则可以方向的忽略——因为最优解中的$A_{ij}$的元素都是非负的。于是，构造拉格朗日算子：\n",
    "\n",
    "$$\n",
    "\\mathcal L(A,\\alpha)=\\sum_{i=1}^{\\lvert S\\rvert}\\sum_{j=1}^{\\lvert S\\rvert}\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_{t}=s_j\\right\\}\\log A_{ij}+\\sum_{i=1}^{\\lvert S\\rvert}\\alpha_i\\left(1-\\sum_{j=1}^{\\lvert S\\rvert}A_{ij}\\right)\n",
    "$$\n",
    "\n",
    "对拉格朗日算子求关于$A_{ij}$的偏导，并将偏导置为$0$得到：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{\\partial\\mathcal L(A,\\alpha)}{\\partial A_{ij}}&=\\frac{\\partial}{\\partial A_{ij}}\\left(\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_{t}=s_j\\right\\}\\log A_{ij}\\right)+\\frac{\\partial}{\\partial A_{ij}}\\alpha_i\\left(1-\\sum_{j=1}^{\\lvert S\\rvert}A_{ij}\\right)\\\\\n",
    "&=\\frac{1}{A_{ij}}\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_{t}=s_j\\right\\}-\\alpha_{i}\\equiv0\\\\\n",
    "&\\Downarrow\\\\\n",
    "A_{ij}&=\\frac{1}{\\alpha_i}\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_{t}=s_j\\right\\}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "再求关于$\\alpha$的偏导，并将偏导置为$0$且将上式代入得到：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{\\partial\\mathcal L(A,\\alpha)}{\\partial\\alpha_i}&=1-\\sum_{j=1}^{\\lvert S\\rvert}A_{ij}\\\\\n",
    "&=1-\\sum_{j=1}^{\\lvert S\\rvert}\\frac{1}{\\alpha_i}\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_{t}=s_j\\right\\}\\equiv0\\\\\n",
    "\\Downarrow\\\\\n",
    "\\alpha_i&=\\sum_{j=1}^{\\lvert S\\rvert}\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_{t}=s_j\\right\\}\\\\\n",
    "&=\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\right\\}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "再将得到的$\\alpha_i$代回前面$A_{ij}$的表达式，就能得到参数$\\hat A_{ij}$的最大似然估计：\n",
    "\n",
    "$$\n",
    "\\hat A_{ij}=\\frac{\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_{t}=s_j\\right\\}}{\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\right\\}}\n",
    "$$\n",
    "\n",
    "这个式子的意义非常直观——系统从状态$i$转移到状态$j$的最大似然概率为“系统从状态$i$转移到状态$j$的次数”除以“系统处于状态$i$的次数”。换句话说，最大似然参数应为系统处于状态$i$时转移到状态$j$的概率。\n",
    "\n",
    "## 2. 隐式马尔可夫模型\n",
    "\n",
    "马尔可夫模型对基于时间序列的数据来说是一种强大的抽象，但是却不能表达一些很常见的场景。我们如何在“无法观测到状态本身，只能得到关于这些状态的某些概率性函数”的情况下，对状态序列进行讨论？这样的场景包括：我们只能观测到单词序列，而不能得到其词性的情况下如何进行词性标注？或是在只能得到音频序列的情况下，如何找出生成这些音频序列的单词序列？我们使用Jason Eisner推荐的例子——冰淇淋气候学——来展开讲解：\n",
    "\n",
    "> 假设我们是生活在2799年的气候学家，正在研究全球气候变暖的历史，但是找不到巴尔的摩的气象记录。不过却找到了Jason Eisner的日记——JE在日记中坚持不懈的记录了他每天吃掉的冰淇淋的数量。*那么，通过这本日记，我们能够得到关于天气的什么信息呢？*\n",
    "\n",
    "隐式马尔可夫模型就可以应用在这一类场景中——我们并不能观测到实际的状态序列（例子中每天的天气），而只能得到由这些状态导致的某种结果（每天吃掉的冰淇淋数量）。\n",
    "\n",
    "形式化的表达就是，HMM是这样一种马尔可夫模型——我们能够得到的是*可观测*输出序列$x=\\left\\{x_1,x_2,\\cdots,x_T\\right\\}$，它们来自输出集合$V=\\left\\{v_1,v_2,\\cdots,v_{\\lvert V\\rvert}\\right\\}$，即$x_t\\in V,t=1\\cdots T$。在上一节中，我们也假设存在状态序列$z=\\left\\{z_1,z_2,\\cdots,z_T\\right\\}$，它们来自状态集合$S=\\left\\{s_1,s_2,\\cdots,s_{\\lvert S\\rvert}\\right\\}$，即$z_t\\in S,t=1\\cdots T$，不过在这个场景中，这些状态是*不可观测*的。而状态$i,j$之间的转换同样表达在一个状态转换矩阵的元素$A_{ij}$中。\n",
    "\n",
    "我们同样会将输出观测值建模为一个关于隐藏状态的函数，不过，我们先要做一个**输出独立性假设**——定义$P\\left(x_t=v_k\\mid z_t=s_j\\right)=P\\left(x_t=v_k\\mid x_1,\\cdots,x_T,z_1,\\cdots,z_T\\right)=B_{jk}$。矩阵$B$表示系统在某时刻处于状态$s_j$的条件下生成可观测输出$v_k$的概率。回到天气的例子中，假设我们得到了一个四天内的冰淇淋消耗日志：$\\vec x=\\left\\{x_1=v_3,x_2=v_2,x_3=v_1,x_4=v_2\\right\\}$。其中输出集合中有这些观测值代码所表示的意义，即$V=\\left\\{v_1=1\\textit{ice cream},v_2=2\\textit{ice cream},v_3=3\\textit{ice cream}\\right\\}$。那么，HMM可以让我们解答什么问题呢？\n",
    "\n",
    "### 2.1 隐式马尔可夫模型的三个问题\n",
    "\n",
    "关于HMM，我们有三个基本问题想要知道：\n",
    "* 可观测序列的概率是什么（我们有多大可能得到一个$3,2,1,2$的冰淇淋消耗记录）？\n",
    "* 最可能生成这个可观测序列的状态序列是什么样的（这四天到底是什么样的天气）？\n",
    "* 如何从给定观测到的数据中学到HMM模型中的参数$A,B$？\n",
    "\n",
    "### 2.2 可观测序列的概率：前向算法\n",
    "\n",
    "在HMM中，我们假设得到的观测数据的生成过程为：假设在时间序列上存在一组状态$\\vec z$，而我们认为这组状态序列是由以参数$A$为状态转移矩阵的马尔可夫模型生成的。其中，对于每一个时刻$t$，我们都将$x_t$看做是一个关于$z_t$的函数。于是，要计算观测值序列的概率，就需要将在每一个可能的状态序列下生成观测值序列$\\vec x$的概率相加。\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "P\\left(\\vec x; A, B\\right)&=\\sum_{\\vec z}P\\left(\\vec x, \\vec z; A, B\\right)\\\\\n",
    "&=\\sum_{\\vec z}P\\left(\\vec x\\mid\\vec z; A, B\\right)P\\left(\\vec z; A, B\\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "这个式子在任意概率分布下恒成立，而且，如果我们使用HMM假设，就能够进一步化简该式：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "P\\left(\\vec x; A, B\\right)&=\\sum_{\\vec z}P\\left(\\vec x\\mid\\vec z; A, B\\right)P\\left(\\vec z; A, B\\right)\\\\\n",
    "&=\\sum_{\\vec z}\\left(\\prod_{t=1}^T P\\left(\\vec x_t\\mid \\vec z_t; B\\right)\\right)\\left(\\prod_{t=1}^T P\\left(\\vec z_t\\mid \\vec z_{t-1}; A\\right)\\right)\\\\\n",
    "&=\\sum_{t=1}^T\\left(\\prod_{t=1}^T B_{z_t\\ x_t}\\right)\\left(\\prod_{t=1}^T A_{z_{t-1}\\ z_t}\\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "好消息是，该式是一个由参数项组成的简单表达式。上面的推导使用的HMM假设有：输出独立性假设、马尔可夫假设、及平稳过程假设，式子的第二步就是由这些假设推导出的。坏消息是，这个求和是对于所有可能的$\\vec z$，而对于每一个时刻，$z_t$可以取到$\\lvert S\\rvert$种不同的值，于是这个求和就需要消耗$O\\left(\\lvert S\\rvert^T\\right)$次操作。\n",
    "\n",
    "幸运的是，一个称为**前向算法（Forward Procedure）**的动态规划算法可以让我们快速计算$P\\left(\\vec x; A, B\\right)$。首先我们定义$\\alpha_i(t)=P(x_1,x_2,\\cdots,x_t,z_t=s_i; A, B)$，这个量表示：系统在$t$时刻时处于$s_i$状态时，从初始时刻到$t$时刻中对应的所有观测值的总概率（计算时$t$时刻可以赋任意状态）。如果我们能够得到这个量，则观测值全集的概率$P\\left(\\vec x\\right)$可以表示为：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "P\\left(\\vec x; A, B\\right)&=P\\left(x_1,x_2,\\cdots,x_T; A, B\\right)\\\\\n",
    "&=\\sum_{i=1}^{\\lvert S\\rvert}P\\left(x_1,x_2,\\cdots,x_T,z_T=s_i; A, B\\right)\\\\\n",
    "&=\\sum_{i=1}^{\\lvert S\\rvert}\\alpha_i(T)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "该算法介绍了一种高效计算$\\alpha_i(t)$的方法，对于系统的每一个时刻只需要进行$O\\left(\\lvert S\\rvert\\right)$次操作。那么，对于计算可观测序列的总概率$P\\left(\\vec x; A, B\\right)$，这个算法的总体复杂度为$O\\left(\\lvert S\\rvert\\cdot T\\right)$。\n",
    "\n",
    "$$\n",
    "\\begin{array}\n",
    "{l}\n",
    "\\hline\n",
    "\\textbf{算法 1 }\\textrm{使用前向算法计算}\\alpha_i(t)\\\\\n",
    "\\hline\n",
    "\\textrm{1. 初始状态: }\\alpha_i(0)=A_{0i},\\ i=1,\\cdots,\\lvert S\\rvert\\\\\n",
    "\\textrm{2. 递推公式: }\\displaystyle\\alpha_j(t)=\\sum_{i=1}^{\\lvert S\\rvert}\\alpha_i(t-1)A_{ij}B_{jx_t},\\ j=1,\\cdots,\\lvert S\\rvert, t=1,\\cdots,T\\\\\n",
    "\\hline\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "还有一个类似的算法称为**后向算法（Backward Procedure）**，可以通过与前向算法相反的次序来计算相应的概率$\\beta_i(t)=P(x_T,x_{T-1},\\cdots,x_{t+1},z_t=s_i; A, B)$。\n",
    "\n",
    "### 2.3 寻找最大似然状态序列：Viterbi算法\n",
    "\n",
    "关于隐式马尔可夫模型，最常见的问题就是：对于给定的可观测序列$\\vec x\\in V^T$，生成这个可观测序列的最可能的状态序列$\\vec z\\in S^T$是什么？用公式表达就是：\n",
    "\n",
    "$$\n",
    "\\arg\\max_{\\vec z} P\\left(\\vec z\\mid\\vec x; A, B\\right)=\\arg\\max_{\\vec z}\\frac{P\\left(\\vec x,\\vec z; A, B\\right)}{\\sum_{\\vec z} P\\left(\\vec x, \\vec z; A, B\\right)}=\\arg\\max_{\\vec z} P\\left(\\vec x,\\vec z; A, B\\right)\n",
    "$$\n",
    "\n",
    "式子的第一步化简使用了贝叶斯定理，第二步化简则是通过观察得知——分母并不直接依靠$\\vec z$得出。简单的说就是我们将所有可能的状态序列$\\vec z$代入模型，然后选择联合概率最大的那个状态序列。不过这又将导致穷举所有可能的状态序列，而使得计算复杂度上升到$O\\left(\\lvert S\\rvert^T\\right)$。此时，我们会自然的想到上一节前向算法中介绍的动态规划解决方案，而且容易发现，如果将式中的$\\displaystyle\\arg\\max_{\\vec z}$换成$\\displaystyle\\sum_{\\vec z}$，则当前的问题就真的变成了前向算法中要解决的问题了。\n",
    "\n",
    "$$\n",
    "\\begin{array}\n",
    "{l}\n",
    "\\hline\n",
    "\\textbf{算法 2 }\\textrm{将EM算法应用在HMM上}\\\\\n",
    "\\hline\n",
    "\\textrm{重复直到收敛\\{}\\\\\n",
    "\\qquad\\textrm{（E步骤）对每个可能的状态标识}\\vec z\\in S^T\\textrm{令：}\\\\\n",
    "\\qquad\\qquad Q\\left(\\vec z\\right):=p\\left(\\vec z\\mid\\vec x; A, B\\right)\\\\\n",
    "\\qquad\\textrm{（M步骤）令：}\\\\\n",
    "\\qquad\\qquad\n",
    "\\begin{eqnarray}\n",
    "A,B&:=\\arg\\max_{A,B}&\\sum_{\\vec z}Q\\left(\\vec z\\right)\\log\\frac{P\\left(\\vec x,\\vec z; A,B\\right)}{Q\\left(\\vec z\\right)}\\\\\n",
    "&\\mathrm{s.t.}&\\sum_{j=1}^{\\lvert S\\rvert}A_{ij}=1,i=1,\\cdots,\\lvert S\\rvert;\\ A_{ij}\\geq0,\\ i,j=1,\\cdots,\\lvert S\\rvert\\\\\n",
    "& &\\sum_{j=1}^{\\lvert S\\rvert}B_{ik}=1,i=1,\\cdots,\\lvert S\\rvert;\\ B_{ik}\\geq0,\\ i=1,\\cdots,\\lvert S\\rvert,\\ k=1,\\cdots,\\lvert V\\rvert\n",
    "\\end{eqnarray}\\\\\n",
    "\\mathrm{\\}}\\\\\n",
    "\\hline\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "**Viterbi算法**与前向算法十分相似，只不过现在我们不再跟踪观测到某时刻$t$的状态序列的总概率（到时刻$t$为止的整个状态序列生成观测值序列的概率），而仅跟踪最大概率并记录其对应的状态序列。\n",
    "\n",
    "### 2.4 参数学习：EM算法在HMM上的应用\n",
    "\n",
    "关于HMM的最后一个问题就是，对于给定的观测值序列，在模型中对应的状态转换概率矩阵$A$、和可观测值的输出生成概率矩阵$B$最可能是什么？举个例子，在语音识别问题中，用语音识别数据集解出最大似然参数，可以让我们在不需要知道“某段音频信号对应的的最可能状态序列”的情况下有效的训练HMM。\n",
    "\n",
    "在本节，我们将推导出应用在隐式马尔可夫模型上的最大期望算法。这个推导将使用本课程（CS229）[讲义](chapter13.ipynb)中给出的EM算法的公司。上一节的算法给出了基础的EM算法的思路，可以看出M步骤中已经包含了关于合理的$A,B$的概率约束条件。类似在1.1.3节中介绍的（非隐式）马尔可夫模型的最大似然解法，而对于现在的优化问题，我们仍然可以使用拉格朗日乘数法解决。还需要留意的是，在E步骤和M步骤中，都需要穷举所有$\\lvert S\\rvert^T$个可能的状态序列标识$\\vec z$。我们将会使用前面提到的前向/后向算法求出E步骤和M步骤的充分统计量。\n",
    "\n",
    "首先，使用马尔可夫假设重新整理目标函数：\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "A,B&=&\\arg\\max_{A,B}\\sum_{\\vec z}Q\\left(\\vec z\\right)\\log\\frac{P\\left(\\vec x,\\vec z; A, B\\right)}{Q\\left(\\vec z\\right)}\\\\\n",
    "&=&\\arg\\max_{A,B}\\sum_{\\vec z}Q\\left(\\vec z\\right)\\log P\\left(\\vec x,\\vec z;A,B\\right)\\\\\n",
    "&=&\\arg\\max_{A,B}\\sum_{\\vec z}Q\\left(\\vec z\\right)\\log\\left(\\prod_{t=1}^TP\\left(\\vec x_t\\mid\\vec z_t;B\\right)\\right)\\left(\\prod_{t=1}^TP\\left(\\vec z_t\\mid\\vec z_{t-1};A\\right)\\right)\\\\\n",
    "&=&\\arg\\max_{A,B}\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T\\log B_{z_t\\ x_t}+\\log A_{z_{t-1}\\ z_t}\\\\\n",
    "&=&\\arg\\max_{A,B}\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{i=1}^{\\lvert S\\rvert}\\sum_{j=1}^{\\lvert S\\rvert}\\sum_{k=1}^{\\lvert V\\rvert}\\sum_{t=1}^{\\lvert T\\rvert}1\\left\\{z_t=s_j\\land x_t=v_k\\right\\}\\log B_{jk}+1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}\\log A_{ij}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "我们将第一步$\\log$中的除法转为减法并发现参数$A,B$并不影响该项的取值，所以在求最大化参数的过程中可以舍弃；而在第三步中应用了马尔可夫假设；第五步化为指示函数，用状态来索引矩阵$A,B$。与显式马尔可夫模型中最大似然参数的推导一样，我们可以舍弃不等式约束，因为求出的解会自然符合这些非负约束。构造拉格朗日算子：\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\mathcal L(A,B,\\sigma,\\epsilon)&=&\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{i=1}^{\\lvert S\\rvert}\\sum_{j=1}^{\\lvert S\\rvert}\\sum_{k=1}^{\\lvert V\\rvert}\\sum_{t=1}^{\\lvert T\\rvert}1\\left\\{z_t=s_j\\land x_t=v_k\\right\\}\\log B_{jk}+1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}\\log A_{ij}\\\\\n",
    "&&+\\sum_{j=1}^{\\lvert S\\rvert}\\epsilon_j\\left(1-\\sum_{k=1}^{\\lvert V\\rvert}B_{jk}\\right)+\\sum_{i=1}^{\\lvert S\\rvert}\\delta_i\\left(1-\\sum_{j=1}^{\\lvert S\\rvert}A_{ij}\\right)\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "求关于参数$A,B$的偏导并将其置为零：\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\frac{\\partial\\mathcal L(A,B,\\sigma,\\epsilon)}{\\partial A_{ij}}&=&\\sum_{\\vec z}Q\\left(\\vec z\\right)\\frac{1}{A_{ij}}\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}-\\sigma_i\\equiv0\\\\\n",
    "A_{ij}&=&\\frac{1}{\\sigma_i}\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}\\\\\n",
    "\\frac{\\partial\\mathcal L(A,B,\\sigma,\\epsilon)}{\\partial B_{jk}}&=&\\sum_{\\vec z}Q\\left(\\vec z\\right)\\frac{1}{B_{jk}}\\sum_{t=1}^T1\\left\\{z_t=s_j\\land x_t=v_k\\right\\}-\\epsilon_j\\equiv0\\\\\n",
    "B_{jk}&=&\\frac{1}{\\epsilon_j}\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_t=s_j\\land x_t=v_k\\right\\}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "再求关于拉格朗日乘数的偏导并将其置为零，然后再带入刚才求得的$A_{ij},B_{jk}$：\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\frac{\\partial\\mathcal L(A,B,\\sigma,\\epsilon)}{\\partial \\sigma_i}&=&1-\\sum_{j=1}^{\\lvert S\\rvert}A_{ij}\\\\\n",
    "&=&1-\\sum_{j=1}^{\\lvert S\\rvert}\\frac{1}{\\sigma_i}\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}\\equiv0\\\\\n",
    "\\sigma_i&=&\\sum_{j=1}^{\\lvert S\\rvert}\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}\\\\\n",
    "&=&\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\right\\}\\\\\n",
    "\\frac{\\partial\\mathcal L(A,B,\\sigma,\\epsilon)}{\\partial \\epsilon_j}&=&1-\\sum_{k=1}^{\\lvert V\\rvert}B_{jk}\\\\\n",
    "&=&1-\\sum_{k=1}^{\\lvert V\\rvert}\\frac{1}{\\epsilon_j}\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_t=s_j\\land x_t=v_k\\right\\}\\equiv0\\\\\n",
    "\\epsilon_j&=&\\sum_{k=1}^{\\lvert V\\rvert}\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_t=s_j\\land x_t=v_k\\right\\}\\\\\n",
    "&=&\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_t=s_j\\right\\}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "再将求得的拉格朗日乘数表达式带回$A_{ij},B_{jk}$的表达式，便可求得能够最大化目标函数的参数$\\hat A,\\hat B$：\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\hat A_{ij}&=&\\frac{\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}}{\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\right\\}}\\\\\n",
    "\\hat B_{jk}&=&\\frac{\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_t=s_j\\land x_t=v_k\\right\\}}{\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_t=s_j\\right\\}}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "不幸的是，式中关于$\\vec z$的求和都需要计算所有可能的状态序列$\\vec z\\in S^T$。回忆$Q\\left(\\vec z\\right)$项在E步骤中定义为$P\\left(\\vec z\\mid\\vec x; A,B\\right)$，它出现在状态序列的最后一步，其参数为$A,B$。现在来考虑如何用前向/后向算法中的$\\alpha_i(t),\\beta_j(t)$将$\\hat A_{ij}$的分子表示出来。\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}&=&\\sum_{t=1}^T\\sum_{\\vec z}1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}Q\\left(\\vec z\\right)\\\\\n",
    "&=&\\sum_{t=1}^T\\sum_{\\vec z}1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}P\\left(\\vec z\\mid\\vec x; A,B\\right)\\\\\n",
    "&=&\\frac{1}{P\\left(\\vec x; A,B\\right)}\\sum_{t=1}^T\\sum_{\\vec z}1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}P\\left(\\vec z, \\vec x; A,B\\right)\\\\\n",
    "&=&\\frac{1}{P\\left(\\vec x; A,B\\right)}\\sum_{t=1}^T\\alpha_i(t)A_{ij}B_{jx_t}\\beta_j(t+1)\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "在前两步中，我们交换了求和顺序，并代入了$Q$的定义；在第三步中使用了贝叶斯公式；第四步中代入了$\\alpha,\\beta,A,B$的定义。类似的，分母也可以表示为分子在$j$到$\\lvert S\\rvert$上的求和：\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\right\\}&=&\\sum_{j=1}^{\\lvert S\\rvert}\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}\\\\\n",
    "&=&\\frac{1}{P\\left(\\vec x;A,B\\right)}\\sum_{j=1}^{\\lvert S\\rvert}\\sum_{t=1}^T\\alpha_i(t)A_{ij}B_{jx_t}\\beta_j(t+1)\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "将上面求得的表达式结合起来，就可以描述出最大似然状态转换参数$\\hat A_{ij}$，而且不需要穷举所有可能的状态序列：\n",
    "\n",
    "$$\n",
    "\\hat A_{ij}=\\frac{\\sum_{t=1}^T\\alpha_i(t)A_{ij}B_{jx_t}\\beta_j(t+1)}{\\sum_{j=1}^{\\lvert S\\rvert}\\sum_{t=1}^T\\alpha_i(t)A_{ij}B_{jx_t}\\beta_j(t+1)}\n",
    "$$\n",
    "\n",
    "类似的，我们可以将$\\hat B_{jk}$的分子表示为：\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_t=s_j\\land x_t=v_k\\right\\}&=&\\frac{1}{P\\left(\\vec x;A,B\\right)}\\sum_{t=1}^T\\sum_{\\vec z}1\\left\\{z_t=s_j\\land x_t=v_k\\right\\}P\\left(\\vec z,\\vec x;A,B\\right)\\\\\n",
    "&=&\\frac{1}{P\\left(\\vec x;A,B\\right)}\\sum_{i=1}^{\\lvert S\\rvert}\\sum_{t=1}^T\\sum_{\\vec z}1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\land x_t=v_k\\right\\}P\\left(\\vec z,\\vec x;A,B\\right)\\\\\n",
    "&=&\\frac{1}{P\\left(\\vec x;A,B\\right)}\\sum_{i=1}^{\\lvert S\\rvert}\\sum_{t=1}^T1\\left\\{x_t=v_k\\right\\}\\alpha_i(t)A_{ij}B_{jx_t}\\beta_j(t+1)\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "$\\hat B_{jk}$的分母表示为：\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\sum_{\\vec z}Q\\left(\\vec z\\right)\\sum_{t=1}^T1\\left\\{z_t=s_j\\right\\}&=&\\frac{1}{P\\left(\\vec x;A,B\\right)}\\sum_{i=1}^{\\lvert S\\rvert}\\sum_{t=1}^T\\sum_{\\vec z}1\\left\\{z_{t-1}=s_i\\land z_t=s_j\\right\\}P\\left(\\vec z,\\vec x;A,B\\right)\\\\\n",
    "&=&\\frac{1}{P\\left(\\vec x;A,B\\right)}\\sum_{i=1}^{\\lvert S\\rvert}\\sum_{t=1}^T\\alpha_i(t)A_{ij}B_{jx_t}\\beta_j(t+1)\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "结合上面的表达式，就可以得到最大似然可观测值生成参数$\\hat B_{jk}$：\n",
    "\n",
    "$$\n",
    "\\hat B_{jk}=\\frac{\\sum_{i=1}^{\\lvert S\\rvert}\\sum_{t=1}^T1\\left\\{x_t=v_k\\right\\}\\alpha_i(t)A_{ij}B_{jx_t}\\beta_j(t+1)}{\\sum_{i=1}^{\\lvert S\\rvert}\\sum_{t=1}^T\\alpha_i(t)A_{ij}B_{jx_t}\\beta_j(t+1)}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{array}\n",
    "{l}\n",
    "\\hline\n",
    "\\textbf{算法 3 }\\textrm{HMM参数学习上的前向-后向算法}\\\\\n",
    "\\hline\n",
    "\\textrm{初始化：令}A,B\\textrm{为随机的有效概率矩阵，其中}A_{i0}=0,B_{0k}=0,\\ i=1,\\cdots,\\lvert S\\rvert,k=1,\\cdots,\\lvert V\\rvert\\textrm{。}\\\\\n",
    "\\textrm{重复直到收敛\\{}\\\\\n",
    "\\qquad\\textrm{（E步骤）运行前向及后向算法，计算}\\alpha_i,\\beta_i,\\ i=1,\\cdots,\\lvert S\\rvert\\textrm{令：}\\\\\n",
    "\\qquad\\qquad \\gamma_t(i,j):=\\alpha_i(t)A_{ij}B_{jx_t}\\beta_j(t+1)\\\\\n",
    "\\qquad\\textrm{（M步骤）重新计算参数的最大似然估计：}\\\\\n",
    "\\qquad\\qquad\n",
    "\\begin{eqnarray}\n",
    "A_{ij}&=&\\frac{\\sum_{t=1}^T\\gamma_t(ij)}{\\sum_{j=1}^{\\lvert S\\rvert}\\sum_{t=1}^T\\gamma_t(ij)}\\\\\n",
    "B_{jk}&=&\\frac{\\sum_{i=1}^{\\lvert S\\rvert}\\sum_{t=1}^T1\\left\\{x_t=v_k\\right\\}\\gamma_t(ij)}{\\sum_{i=1}^{\\lvert S\\rvert}\\sum_{t=1}^T\\gamma_t(ij)}\n",
    "\\end{eqnarray}\\\\\n",
    "\\mathrm{\\}}\\\\\n",
    "\\hline\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "这个算法是另一种形式的前向-后向算法，也被称作**Baum-Welch算法**。在E步骤中，我们并不会直接计算每一个可能的$\\vec z\\in S^T$对应的$Q\\left(\\vec z\\right)$，而是求出一个充分统计量$\\gamma_t(ij)=\\alpha_i(t)A_{ij}B_{jx_t}\\beta_j(t+1)$。这个量对应的是“在$t$时刻，对于所有给定的观测序列$\\vec x$，从状态$s_i$转换到状态$s_j$”的概率。而且，算法中推导出的关于$A_{ij},B_{jk}$的表达式的意义也很直观。其中$A_{ij}$是“从状态$s_i$转换到状态$s_j$的次数”除以“状态$s_i$出现的次数”；$B_{jk}$是“由状态$s_j$生成可观测值$v_k$的次数”除以“状态s_j出现的次数”。\n",
    "\n",
    "同其他EM算法处理的问题一样，HMM上的参数学习也是一个非凸问题，存在着很多局部最大值。EM算法可能会根据初始值的不同，收敛于不同的最大值处。所以，用不同的初始值多运行几次是一种比较常见的做法。另外，对由$A,B$表示的概率分布做平滑处理也很重要，这样可以避免某个“状态转换的过程”或“状态生成观测值的”是$0$概率。\n",
    "\n",
    "### 2.5 深入学习材料\n",
    "\n",
    "网上有很多隐式马尔可夫模型的学习资源，对于自然语言处理（NLP）参考Jurafsky & Martin的《Speech and Language Processing》第二版；或是Manning & Schütze的《Foundations of Statistical Natural Language Processing》；或是Eisner的《HMM in a spreadsheet》——一个仅靠spreadsheet建立的轻量级交互式HMM（主要用于为前向-后向算法提供一个直观的例子）。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}