17. 马尔科夫与隐马尔科夫模型

17.1 前言

本章，我们将讨论任意长度 $T$ 的序列观测数据 $X_1,…,X_T$ 的概率模型。这类模型可应用在计算机生物学，自然语言处理以及时序预测等领域。我们重点关注观测值在离散“时间”下产生的情况，尽管此处的”时间”也可以指一个序列中的不同位置。

17.2 马尔科夫模型

根据10.2.2节的内容，马尔科夫链背后的基本思想是：假设 $Xt$ 包含所有用来预测未来信息的知识(换句话说，我们可以将它作为一个充分统计量)。如果我们假设时间是离散的，那么观测序列的联合概率分布可以表示为: $$ p(X{1:T})=p(X1)p(X_2|X_1)p(X_3|X_2)…=p(X_1)\prod{t=2}^Tp(Xt|X{t-1})\tag{17.1} $$ 上式被称为马尔科夫链(Markov chain)或马尔科夫模型(Markov model)。

如果假设转移函数$p(Xt|X{t-1})$与时间无关，那么上述模型被称为是同质的(homogeneous), 静态的(stationary)或者时不变的(time-invariant)。这是一个参数共享(parameter tying)的例子，因为同样的参数被用于多个变量的生成。这个假设允许我们可以使用固定数量的参数去生成任意数量的变量，这种模型被称为随机过程(stochastic processes)。

如果假设我们观测到的是离散变量，也就是说$X_t \in {1,…,K}$，此时的马尔科夫链又被称为离散状态或者有限状态马尔科夫链。本章后面的内容将始终遵循该假设。

17.2.1 转移矩阵

如果 $Xt$ 是离散变量，即 $X_t\in{1,…,K}$，那么条件分布 $p(X_t|X{t-1})$ 的形式可以写成一个大小为 $K\times K$ 的转移矩阵(transition matrix) $\bf{A}$，其中 $A{ij}=p(X_t=j|X{t-1}=i)$ 为状态 $j$ 转移到状态 $i$ 的概率。矩阵的每一行满足 $\sumjA{ij}=1$，所以该矩阵又被称为随机矩阵(stochastic matrix)。

一个静态的(参数共享)，有限状态(观测值离散)的马尔科夫链等价于一个随机自动机(stochastic automaton)。通常情况下我们会使用一个有向图来可视化这样的自动机，其中节点代表状态，箭头代表合法的状态转移(对应于矩阵 $\bf{A}$ 中的非零元素)。该有向图被称为状态转移图(state transition diagram)。图中每一条弧线对应的权重代表转移的概率值。举例来说，图17.1(左)表示了一个 $2-\rm{state}$ 马尔科夫链: $$ \mathbf{A}= \begin{pmatrix} 1-\alpha & \alpha \ \beta & 1-\beta \end{pmatrix} \tag{17.2} $$ 图17.1(右)表示了一个 $3-\rm{state}$ 的马尔科夫链: $$ \mathbf{A}= \begin{pmatrix} A{11} & A{12} & 0 \ 0 & A{22} & A{23} \ 0 & 0 & 1 \end{pmatrix} \tag{17.3} $$ 上式被称为$\mathbf{\rm{left-to-right\ transition\ matrix}}$, 在语音识别中被普遍使用(参考17.6.2节)。

转移矩阵中的元素$A{ij}$表示在状态 $i$ 经过单步到达状态 $j$ 的概率。$n$ 步转移矩阵 ${\bf{A}}({\it{n}})$ 定义为: $$ A{ij}(n)\triangleq p(X{t+n}=j|X_t=i) \tag{17.4} $$ 上式表示状态 $i$ 经过 $n$ 步转移后到达状态 $j$ 的概率值。 显然 $\mathbf{A}(1)=\mathbf{A}$。根据$\mathbf{Chapman-Kolmogorov}$等式，有: $$ A{ij}(m+n)=\sum{k=1}^KA{ik}(m)A_{kj}(n) \tag{17.5} $$ 换句话说，状态 $i$ 经过 $m+n$ 步到达状态 $j$ 的概率等于从状态 $i$ 经过 $m$ 步到达状态 $k$, 再从状态 $k$ 经过 $n$ 步到达状态 $j$ ，并对所有的状态 $k$ 求和的概率。我们可以将上述公式写成矩阵乘法的形式: $$ \mathbf{A}(m+n)=\mathbf{A}(m)\mathbf{A}(n) \tag{17.6} $$ 所以 $$ \mathbf{A}(n)=\mathbf{A}\mathbf{A}(n-1)=\mathbf{A}\mathbf{A}\mathbf{A}(n-2)=…=\mathbf{A}^n \tag{17.7} $$ 所以我们可以通过使用转移矩阵的阶乘模拟马尔科夫链的多步转移。

17.2.2 应用：语言模型

马尔科夫的一个重要应用是构建统计语言模型（language models），即对序列单词的概率分布进行建模。我们定义状态空间为英语(或其他语言)中的所有单词。边际概率 $p(Xt=k)$ 被称为 $\mathbf{unigram}$ 统计量。如果我们使用一阶马尔科夫模型，那么 $p(X_t=k|X{t-1}=j)$ 被称为 $\mathbf{bigram}$ 模型。如果我们使用二阶马尔科夫模型，那么 $p(Xt=k|X{t-1}=j,X_{t-2}=i)$ 被称为 $\mathbf{trigram}$ 模型。以此类推，这类模型被称为 $\mathbf{n-gram}$模型。举例来说，图17-2展示了达尔文书籍《物种起源》中字母${a,…,z,-}$ ($-$ 表示空格) 的 $1-\rm{gram}$ 和 $2-\rm{grams}$ 的统计频次。

语言模型可以用来作以下几件事情：

  • 句子补全(Sentence completion)：语言模型可以在给定一个句子的历史单词的情况下预测下一个单词。这可以用来减轻拼写的负担，这对于残疾人士特别重要(可以参考 David Mackay 的 Dasher系统)，或者用在移动设备上。
  • 数据压缩(Data compression): 任何密度模型都可以用来定义一个编码策略，通过给那些高频字符串赋予短的编码词（译者注:参考 赫夫曼编码）。预测模型的精度越高，存储数据所需要的比特数越少。
  • 文本分类(Text classification): 任何密度模型可以用作一个类条件密度模型，进而转换成一个(生成式)分类器。值得注意的是，如果使用 $0-\rm{gram}$ 类条件密度模型(换句话说: 只使用 $\rm{unigram}$ 统计量)，对应的分类器等价于朴素贝叶斯分类器(见3.5节)。
  • 自动论文书写(Automatic essay writing): 我们可以从概率分布$p(x_{1:t})$中进行采样，从而人工生成文本。这是一种评估模型质量的方法。在表17.1中，我们给出了一个从 $4-\rm{gram}$ 模型中生成的一段文本，该模型基于400 million个单词的语料库训练得到。(Mikolov等人描述了一个更好的语言模型，该模型基于一个循环神经网络，可以生成在语义上更合理的文本。)

17.2.2.1 马尔科夫语言模型的最大似然解

我们现在讨论一种从训练集中估计转移矩阵的简单方法。长度为 $T$ 的任意序列的概率分布可以表示为: $$ p(x{1:T}|\mathbf{\theta}) = \pi(x_1)A(x_1,x_2)…A(x{T-1},xT) =\prod{j=1}^{K}(\pij)^{\mathbb{I}(x_i=j)}\prod{t=2}^{T}\prod{j=1}^{K}\prod{k=1}^{K}(A{jk})^{\mathbb{I}(x_t=k,x{t-1}=j)}\tag{*} $$ 因此对于一个序列集合 $\mathcal{D}=(\mathbf{x}1,…\mathbf{x}_N)$，其中 $\mathbf{x}_i=(x{i1},…,x{i,T_i})$ 为长度为 $T_i$，该集合的对数似然为: $$ {\rm{log}}\ p(\mathcal{D}|\mathbf{\theta})=\sum{i=1}^N {\rm{log}}\ p(\mathbf{x}i|\mathbf{\theta}) = \sum{j} Nj^1 {\rm{log}} \ \pi_j + \sum_j\sum_k N{jk} {\rm{log}}\ A{jk} \tag{17.10} $$ 其中我们定义： $$ N_j^1 \triangleq \sum{i=1}^N \mathbb{I}(x{i1}=j),N{jk}\triangleq \sum{i=1}^N\sum{t=1}^{Ti-1}\mathbb{I}(x{i,t}=j, x{i,t+1}=k) \tag{17.11} $$ 所以最大似然估计的结果为: $$ \hat{\pi}_j=\frac{N_j^1}{\sum_jN_j^1},\hat{A}{jk}=\frac{N{jk}}{\sum_kN{jk}} \tag{17.12} $$

上述结果可以直接推广至更高阶的马尔科夫模型。然而，当变量的状态值 $K$ 或者模型阶数 $n$ 很大时，零计数($\rm{zero-counts}$) 问题会变得非常严重。一个 $n-\rm{gram}$ 模型的参数量为 $O(K^n)$。如果在整个语料库中有近$K \sim 50,000$ 个单词，那么 $\rm{bi-gram}$ 模型将具有近 $2.5\ billion$ 个自由参数，对应于所有可能的单词配对。显然我们没有办法在训练集中涵盖所有的配对。然而，我们并不希望将一个在训练集中没有出现的单词预测为完全不可能——这将导致严重的过拟合现象。

一个简单的解决方法是$\rm{add-one}$平滑，在这种方法中，我们在归一化之前，将所有的经验计数加 $1$。章节3.3.4.1给出了这个方法的贝叶斯学派的证明。然而，$\rm{add-one}$ 平滑方法假设所有的 $\rm{n-grams}$ 具备同等的可能性，显然这不符合实际。一个更加复杂的贝叶斯方法将在17.2.2.2节介绍。

一个可选方案是使用更好的先验知识去收集更多的数据。举例来说，基于从网站上收集的 $\rm{one}\ \rm{trillion}$ 个单词，谷歌已经训练了一个 $\rm{n-gram}$ 模型($n=1:5$)。他们的数据解压后超过$100 \rm{GB}$并且已经开源。以下给出了一些数据库中 $\rm{4-grams}$ 的统计数据：

$\rm{serve\ as\ the\ incoming\ 92}$

$\rm{serve\ as\ the\ incubator\ 99}$

$\rm{serve\ as\ the\ independent\ 794}$

$\rm{serve\ as\ the\ index\ 223}$

$\rm{serve\ as\ the\ indication\ 72}$

$\rm{serve\ as\ the\ indicator\ 120}$

$\rm{serve\ as\ the\ indicators\ 45}$

$\rm{serve\ as\ the\ indispensable\ 111}$

$\rm{serve\ as\ the\ indispensible\ 40}$

$\rm{serve\ as\ the\ individual\ 234}$

$…$

尽管这种依靠 “brute force and ignorance”的方式可能是成功的，但却不尽如人意，因为显然这并不是人类本身学习的方式。一种更加优雅的贝叶斯方法需要更少的数据，将在17.2.2.2节介绍。

17.2.2.2 删除插值法的经验贝叶斯版本

一种常见的解决稀疏数据的启发式方法被称为删除插值法(deleted interpolation)。该方法将转移矩阵定义为$\rm{bigram}$频率$f{jk}=N{jk}/Nj$和$\rm{unigram}$频率$f_k=N_k/N$的凸组合： $$ A{jk}=(1-\lambda)f{jk}+\lambda f_k \tag{17.13} $$ 式中$\lambda$通常使用交叉验证的方法得到。然后另一种相关的技术被称为$\mathbf{backoff\ smoothing}$; 其思想是如果$f{jk}$的值过小，我们将”back off”一个更加可靠的估计量$f_k$。

17.3 隐马尔可夫模型

根据10.2.2节的介绍，一个隐马尔可夫模型(hidden Markov model, HMM)由一个离散时间，离散状态马尔科夫链，隐状态$zt\in{1,…,K}$和观测模型$p(\mathbf{x}_t|z_t)$组成。对应的联合概率分布具备如下形式： $$ p(\mathbf{z}{1:T},\mathbf{x}{1:T})=p(\mathbf{z}{1:T})p(\mathbf{x}{1:T}|\mathbf{z}{1:T})=\left[ p(z1)\prod{t=2}^Tp(zt|z{t-1}) \right] \left[ \prod_{t=1}^Tp(\mathbf{x}_t|z_t) \right] \tag{17.39} $$ HMM模型中的观测值可以是离散或是连续的。如果是离散变量，那么观测模型对应于一个观测矩阵: $$ p(\mathbf{x}_t=l|z_t=k,\mathbf{\theta})=B(k,l) \tag{17.40} $$ 如果观测值是连续的，那么观测模型通常是一个条件高斯分布: $$ p(\mathbf{x}_t|z_t=k,\mathbf{\theta})=\mathcal{N}(\mathbf{x}_t|\mathbf{\mu}_t,\mathbf{\Sigma}_k) \tag{17.41} $$ 图17.7展示了一个具备3种状态的模型，每种状态对应一个不同的高斯分布。最终的模型与一个高斯混合模型相似，不同的是前者簇的中心具备$\rm{Markovian\ dynamics}$。(的确，HMMs有时也被称为 $\bf{Markov\ switching\ models}$)。我们发现，我们倾向于在同一个位置得到多个观测值，然后再跳转到一个新的簇中心。

17.3.1 HMMs的应用

HMMs可以作为序列数据的黑箱密度模型。相较于马尔科夫模型，HMMs的优势在于，它可以间接地通过隐变量，来表示观测变量在较长范围内的相关性。特别地，在HMMs中，并不需要观测变量本身满足马尔科夫性质。这样的黑箱模型很适用于时间序列的预测。他们也可以用来在一个生成式分类器中定义类条件密度。

然而，更常见的情况是，我们赋予隐变量一些特定的含义，并且尝试通过观测值来估计隐变量，即：在线场景下计算 $p(zt|\mathbf{x}{1:t})$，或者离线场景下计算 $p(zt|\mathbf{x}{1:T})$(17.4.1节将讨论这两种场景的区别)。接下来我们将介绍一些HMMs基于该方法的应用案例：

  • 自动语音识别(Automatic speech recognition): 在这种场景下，$\mathbf{x}t$表示从语音信号中提取出来的特征，$z_t$表示语音对应的单词。转移模型(transition model) $p(z_t|z{t-1})$ 表示语言模型(language model)，观测模型 $p(\mathbf{x}_t|z_t)$表示声学模型。
  • 行为识别(Activity recognition): 此处 $\mathbf{x}_t$ 表示从一个视频帧中提取的特征，$z_t$ 表示视频中包含的行为的类型(比如: 奔跑，行走，坐下等等)。
  • 词性标注(Part of speech tagging): 此处 $\mathbf{x}_t$ 表示一个单词，$z_t$ 表示它的词性(名词，动词，形容词等等)。19.6.2.1节将介绍更多关于POS和相关任务的细节。
  • 基因识别(Gene finding): 此处$x_t$表示DNA核苷酸 (A,C,G,T)，$z_t$表示我们是否在一个基因编码区内。
  • 蛋白质序列比对(Protein sequence alignment): $x_t$ 代表氨基酸，$z_t$ 代表该氨基酸是否与该位置的潜在共有序列匹配。此模型称为配置HMM(profile HMM)，如图17.8所示。 HMM具有3种状态，分别称为匹配，插入和删除。如果$z_t$是匹配状态，则$x_t$等于共识的第$t$个值。如果$z_t$是插入状态，则$x_t$由与共有序列无关的均匀分布生成。如果$z_t$是删除状态，则$x_t =-$。这样，我们可以生成不同长度的共有序列的嘈杂副本。在图17.8(a)中，共识是 “$\rm{AGC}$”，我们在下面看到其各种版本。通过状态转换图的路径（如图17.8（b）所示）指定如何使序列与共识一致，例如，对于而言，最可能的路径是D，D，I，I，I，M。这意味着我们删除共有序列的A和G部分，插入3个A，然后匹配最终的C。我们可以通过计算此类过渡的数量以及每种类型的排放的数量来估算模型参数状态，如图17.8（c）所示。有关训练HMM的更多信息，请参见第17.5节；有关配置文件HMM的详细信息，请参见（Durbin等1998）。

值得注意的是，在上面的任务中，HMMs的判别式版本，条件随机场可能更加适合。第19章将介绍更多细节。

17.4 HMMs的推理

我们现在讨论在参数已知的情况下，如何预测HMM模型中隐状态序列。同样的算法也可以用于其他链式结构的图模型，比如链式条件随机场(见19.6.1节)。在第20章，我们将这些方法泛化到任意的图。在17.5.2节，我们将展示如何在参数估计的过程中利用这些推理输出。

17.4.1 时间模型中推理问题的类型

对于HMM(或者更一般的状态空间模型,SSM)模型，有几种不同的推理任务。为了说明这些区别，我们会考虑一个叫$\bf{occasionally\ dishonest\ casino}$的例子。在这个模型中，$x_t\in{1,2,…,6}$代表骰子朝上的数字，$z_t$代表正在使用的骰子的编号。在大部分情况下，会使用一个公平的骰子，即$z=1$，但有的时候，在较短时间内也会使用一个不平衡的骰子，即$z=2$。如果$z=1$，那么观测值将服从一个均匀的多伯努利分布，分布的支撑集为${1,…,6}$。如果$z=2$, 观测值的分布将倾向于 $6$ (见图17.9)。如果我们从这个模型中采样，我们会得到如下的观察值：

此处 “rolls” 表示观测到的符号，”die”表示隐状态($\rm{L}$表示不平衡，$\rm{F}$表示平衡的骰子)。因此，我们看到该模型生成了一个符号序列，但是分布的统计信息时不时地突然变化。在一个特定的应用场景中，我们仅仅观察到rolls，然后想推测使用了哪个骰子。但有几种不同的推测方法，总结如下：

  • 滤波(Filtering)表示随着数据流的输入，按在线或循环模式计算置信状态(belief state) $p(zt|\mathbf{x}{1:t})$。这种方式被称为”滤波”，因为它减轻噪声的影响，而不仅仅是简单地基于当前的估计值 $p(z_t|\mathbf{x}_t)$ 去估计隐变量的值。接下来我们将看到在序列模型中如何只是使用贝叶斯法则实现滤波。图17.10(a)给出了一个例子。

  • 平滑(Smoothing) 表示在离线状态下给定所有的线索，计算$p(zt|\mathbf{x}{1:T})$。图17.10(b)给出了一个例子。通过使用过去和未来的数据，我们预测的不确定性可以大幅度地降低。为了直观地理解这一点，考虑一个侦探试图找出谁犯了罪。 当他路过犯罪现场时，不确定性很高，直到找到关键线索为止。 然后他有一个“啊哈”时刻，他的不确定性降低了，事后看来，以前所有令人困惑的观察结果都很容易解释。

  • 固定延迟平滑(Fixed lag smoothing)是在线和离线估算之间的有趣折衷，它涉及计算 $p(z{t-l}|\mathbf{x}{1:t})$，其中$l\gt0$称为滞后值(tag)。 这样可以提供比滤波方法更好的性能，但会产生一些延迟。 通过改变滞后的大小，我们可以权衡准确性与延迟。

  • 预测(Prediction) 我们不希望像固定滞后平滑那样预测给定未来的过去，而是希望预测给定过去的未来，即计算 $p(z{t+h}|\mathbf{x}{1:t})$，其中$h> 0$称为预测范围(horizon) 。 例如，假设$h = 2$； 然后我们有 $$ p(z{t+2}|\mathbf{x}{1:t})=\sum{z{t+1}}\sum{z{t}}p(z{t+2}|z{t+1})p(z{t+1}|z_t)p(z_t|\mathbf{x}{1:t}) \tag{17.42} $$ 执行此计算非常简单：我们只需使用转移矩阵并将其应用于当前的置信状态。$p(z{t+h}|\mathbf{x}{1:t})$是关于未来隐藏状态的预测。 可以使用以下方法将其转换为对未来观测的预测: $$ p(\mathbf{x}{t+h}|\mathbf{x}{1:t})=\sum{z{t+h}}p(\mathbf{x}{t+h}|z{t+h})p(z{t+h}|\mathbf{x}{1:t}) \tag{17.43} $$ 上式为后验预测密度，可以用来作时间序列预测。图17.11给出了滤波，平滑和预测之间关系的示意图。

  • MAP估计 是指计算 ${\rm{argmax}}{\mathbf{z}{1:T}}p(\mathbf{z}{1:T}|\mathbf{x}{1:T})$，即概率最高的状态序列。在HMMs中，上述问题被称为$\rm{Viterbi\ decoding}$ (见17.4.4节)。图17.10说明了滤波，平滑以及MAP解码之间的区别。不难发现，使用离线平滑估计得到的结果确实比在线滤波方法得到的结果更加平滑。如果我们使用0.5作为预测的截断阈值，并且与真实的序列进行对比，我们发现滤波方法的误差为$71/300$，平滑方法的误差为$49/300$，MAP方法的误差为$60/300$。平滑方法优于Viterbi的结论并不意外，因为最小化$\rm{bit-error}$率的最优方式是基于后验边际分布使用阈值(见5.7.1.1节)。然而，在某些应用下，我们更倾向于使用Viterbi解码，我们将在17.4.4节进行讨论。

  • 后验采样：如果存在多种关于数据的合理解释，我们可以基于后验分布 $\mathbf{z}{1:T} {\sim} p(\mathbf{z}{1:T}|\mathbf{x}_{1:T})$进行采样。相较于基于使用平滑方法得到的边际分布采样的序列，这些序列包含更多的信息。

  • 证据的概率(Probability of the evidence): 通过对所有的隐变量序列求和 $p(\mathbf{x}{1:T})=\sum{\mathbf{z}{1:T}}p(\mathbf{z}{1:T},\mathbf{x}{1:T})$得到证据的概率$p(\mathbf{x}{1:T})$。对于$\rm{model-based}$聚类，$\rm{anomaly\ detection}$等任务，该值可以用于对序列进行分类 (举例来说，如果HMM被用来作为类条件密度函数)。

17.4.2 前向算法

我们现在描述在HMM中如何递归地计算滤波算法中的边际分布$p(zt|\mathbf{x}{1:t})$。

算法分两个步骤。首先是预测阶段，在这个阶段，我们计算 $\rm{one-step-ahead\ predictive\ density}$，该值可以用作 $t$ 时刻的最新的先验分布： $$ p(zt=j|\mathbf{x}{1:t-1})=\sum{i}p(z_t=j|z{t-1}=i)p(z{t-1}|\mathbf{x}{1:t-1}) \tag{17.44} $$ 接下来是更新步骤，在这个步骤中，我们使用贝叶斯法则将 $t$ 时刻的观测值吸收进来： $$ \begin{align} \alphat(j) & \triangleq p(z_t=j|\mathbf{x}{1:t})=p(zt=j|\mathbf{x}_t,\mathbf{x}{1:t-1}) \tag{17.45} \ & = \frac{1}{Zt}p(\mathbf{x}_t|z_t=j,\cancel{\mathbf{x{1:t-1}}})p(zt=j|\mathbf{x}{1:t-1}) \tag{17.46} \end{align} $$ 其中归一化常数为: $$ Zt \triangleq p(\mathbf{x}_t|\mathbf{x}{1:t-1})=\sumj p(z_t=j|\mathbf{x}{1:t-1})p(\mathbf{x}t|z_t=j) \tag{17.47} $$ 该过程被称为$\rm{predict-update\ cycle}$。分布 $p(z_t|\mathbf{x}{1:t})$ 被称为 $t$ 时刻的(滤波)置信状态，是一个长度为$K$的向量，通常表示为$\mathbf{\alpha}t$。在矩阵向量的符号表达中，我们可以将更新部分写成如下简单的形式： $$ \mathbf{\alpha}_t \propto \psi_t \odot (\mathbf{\Psi}^T\mathbf{\alpha}{t-1}) \tag{17.48} $$ 其中$\psit(j)=p(\mathbf{x}_t|z_t=j)$表示 $t$ 时刻的局部证据，$\Psi(i,j)=p(z_t=j|z{t-1}=i)$表示转移矩阵，$\mathbf{u} \ \odot\ \mathbf{v}$ 表示$\rm{Hadamard\ product}$，表示逐元素向量乘法。算法17.1给出了伪代码。

除了计算隐状态，我们可以使用该算法计算证据的对数概率分布： $$ {\rm{log}}\ p(\mathbf{x}{1:T}|\mathbf{\theta})=\sum{t=1}^T {\rm{log}}\ p(\mathbf{x}t|\mathbf{x}{1:t-1})=\sum_{t=1}^T {\rm{log}}\ Z_t \tag{17.49} $$ (之所以使用对数是为了避免数值溢出。)

17.4.3 前向-后向算法