Python-马尔科夫模型实用指南-全-
Python 马尔科夫模型实用指南（全）
原文：
annas-archive.org/md5/5e11750eef3992654a98af283d98a1e8译者：飞龙
前言
使用隐马尔可夫模型（HMMs）是一种用于建模具有未观察状态的马尔可夫过程的技术。它们是动态贝叶斯网络（DBNs）的一种特殊情况，但在许多问题中表现良好。HMMs 广泛应用于语音识别领域，因为它们能够提供一种非常自然的方式来建模语音数据。本书从 HMM 的理论基础开始，介绍了概率论的基本知识，随后讨论了 HMM 的不同应用。
本书适用读者
具备基本的概率论、线性代数和微积分知识将使阅读本书更加轻松。对于代码示例，期望读者具有基本的 Python 编程基础。
本书内容概述
第一章，马尔可夫过程简介，首先讨论了基本的概率论，然后介绍了马尔可夫链。该章节还讨论了基于连续或离散状态及时间间隔分类的不同类型的马尔可夫链。
第二章，隐马尔可夫模型，基于马尔可夫过程和动态贝叶斯网络（DBNs）的概念，引入了 HMM 的基本概念。
第三章，状态推断—预测状态，介绍了可以用来预测定义好的 HMM 状态的算法。章节介绍了前向算法、后向算法、前向后向算法和维特比算法。
第四章，基于最大似然的参数推断，讨论了最大似然学习的基础。然后，章节介绍了在 HMM 案例中应用最大似然学习，并介绍了维特比学习算法和鲍姆-韦尔奇算法。
第五章，基于贝叶斯方法的参数推断，首先介绍了贝叶斯学习的基本概念。然后，章节将这些概念应用于隐马尔可夫模型（HMMs）的案例，并讨论了用于贝叶斯学习的不同近似方法。
第六章，时间序列预测，讨论了 HMM 在时间序列数据中的应用。章节通过股票价格波动的例子，尝试使用 HMM 对其进行建模。
第七章，自然语言处理，讨论了 HMM 在语音识别领域的应用。章节重点讨论了两个主要应用领域：词性标注和语音识别。
第八章，用于图像处理的二维 HMM，介绍了二维 HMM 的概念，并讨论了它们在图像处理领域的应用。
第九章，马尔可夫决策过程，介绍了强化学习的基本概念，随后讲解了马尔可夫决策过程，并介绍了贝尔曼方程来求解它们。
为了充分利用本书的内容
你需要在你的计算机上安装 Python 3.4，以确保顺利完成本书的各章节内容。
下载示例代码文件
你可以从你的账户中下载本书的示例代码文件，网址为 www.packt.com。如果你从其他地方购买了本书，可以访问 www.packt.com/support，并注册后将文件直接通过电子邮件发送给你。
你可以按照以下步骤下载代码文件：
-
登录或注册 www.packt.com。
-
选择“支持”选项卡。
-
点击“代码下载与勘误”。
-
在搜索框中输入书名，按照屏幕上的指示操作。
一旦文件下载完成，请确保使用以下最新版本解压或提取文件夹：
-
WinRAR/7-Zip for Windows
-
Zipeg/iZip/UnRarX for Mac
-
7-Zip/PeaZip for Linux
本书的代码包也托管在 GitHub 上，地址为 github.com/PacktPublishing/Hands-On-Markov-Models-with-Python。如果代码有更新，它将会在现有的 GitHub 仓库中更新。
我们还提供了其他来自我们丰富书籍和视频目录的代码包，网址为 github.com/PacktPublishing/。快去看看吧！
下载彩色图像
我们还提供了一份 PDF 文件，其中包含本书中使用的截图/图表的彩色图像。你可以在此处下载： www.packtpub.com/sites/default/files/downloads/9781788625449_ColorImages.pdf。
使用的约定
本书中使用了多种文本约定。
CodeInText：表示文本中的代码词汇、数据库表名、文件夹名、文件名、文件扩展名、路径名、虚拟 URL、用户输入和 Twitter 用户名。例如：“现在，我们来定义MFCTagger类。”
代码块的格式如下：
from hmmlearn.hmm import GaussianHMM
import numpy as np
import matplotlib.pyplot as plt
任何命令行输入或输出均按以下格式书写：
pip install matplotlib datetime
粗体：表示新术语、重要词汇或你在屏幕上看到的词汇。例如，菜单或对话框中的词汇会像这样出现在文本中。举个例子：“结果的可能状态也被称为随机变量的定义域。”
警告或重要说明将以此格式出现。
提示和技巧如下所示。
联系我们
我们始终欢迎读者的反馈。
一般反馈：如果你对本书的任何内容有疑问，请在邮件主题中注明书名，并通过 customercare@packtpub.com 向我们发送邮件。
勘误：虽然我们已尽全力确保内容的准确性，但错误仍然可能发生。如果您在本书中发现任何错误，我们将感激不尽，如果您能向我们报告。请访问 www.packt.com/submit-errata，选择您的书籍，点击“勘误提交表单”链接，并填写详细信息。
盗版：如果您在互联网上遇到任何我们作品的非法复制品，无论以何种形式，我们将不胜感激，如果您能提供相关地址或网站名称。请通过copyright@packt.com与我们联系，并附上该材料的链接。
如果您有兴趣成为作者：如果您在某个领域拥有专业知识，并且有意写书或为书籍贡献内容，请访问 authors.packtpub.com。
评论
请留下评论。在您阅读并使用本书后，不妨在您购买该书的网站上留下评论。潜在的读者可以看到并参考您的公正意见来做出购买决策，我们在 Packt 也可以了解您对我们产品的看法，而我们的作者也能看到您对其书籍的反馈。谢谢！
欲了解更多关于 Packt 的信息，请访问 packt.com。
第一章：马尔可夫过程简介
在本章中，我们将展开理解隐马尔可夫模型（HMM）所需的基本概念。我们将涵盖以下主题：
-
随机过程
-
马尔可夫过程
-
马尔可夫链或离散时间马尔可夫过程
-
连续时间马尔可夫链
随机变量
就像我们在统计学中常做的那样，让我们从一个简单的掷骰子例子开始。如果我们考虑掷一个公平的骰子，骰子的结果可以是从 1 到 6 的任何一个数，且是随机的。为了表示这种情况（本例中是掷骰子的结果），在数学中我们使用随机变量的概念。在我们的日常生活中，我们会遇到很多这样的变量。另一个例子是点餐。在这种情况下，结果可以是菜单上的任何一道菜。一般来说，随机变量是一个其可能的值是随机现象结果的变量。结果的可能状态也被称为随机变量的定义域，而结果则是基于在随机变量定义域上定义的概率分布。
回到掷骰子的例子，随机变量的结果O的定义域是domain(O) = (1, 2, 3, 4, 5, 6)，其概率分布为均匀分布，P(o) = 1/6 ∀ o ∈ domain(O)。类似地，在餐厅的例子中，随机变量选择菜肴的定义域将是菜单上的每一项菜品，而概率分布则取决于你的食物偏好。在之前的两个例子中，随机变量的定义域是离散的；这种随机变量称为离散随机变量。但也有可能定义域是一个连续的空间。例如，考虑表示明天 Google 股票价格的随机变量。该随机变量的定义域将是所有正实数，并且大部分概率质量分布在今天价格的±5%范围内。这种随机变量称为连续随机变量。
随机过程
在上一节中，我们讨论了能够在数学上表示单一随机现象结果的随机变量。但如果我们想要表示这些随机事件在一段时间内或者实验持续过程中发生的情况，该怎么办呢？例如，假设我们想要表示一天内每小时的股价，或者想要表示一个从某高度在真空中自由下落的球在每秒钟的高度。对于这种情况，我们需要一组随机变量，每个随机变量表示在给定时间点的结果。这些表示随机现象在一段时间内的随机变量集合也被称为随机过程。值得注意的是，所有这些随机变量的定义域是相同的。因此，我们也可以把这个过程看作是状态的变化。
在这里，我们讨论的是不同时间点的随机变量，但并不一定每次都以时间为基础。它也可以是其他某种事件。但由于在大多数情况下，通常是时间，因此讨论随机过程时更容易使用时间来表示任何此类事件。如果它随着其他事件变化而非时间变化，创建模型时相同的概念也会适用。
现在让我们更详细地讨论之前的两个例子。首先是从真空中将球从一定高度掉落的例子，如果我们知道重力的准确数值和球体掉落的高度，那么我们可以通过牛顿运动定律，在每秒的时间间隔内精确地确定球的位置。
这类随机过程中，我们可以通过初始条件（在此案例中是掉落球，初速度为零）和系统参数（在此案例中是重力的数值），确定每个随机变量的状态，这类过程被称为确定性随机过程（通常称为确定性过程）。
现在让我们看第二个例子：表示股票价格随时间的变化。在这种情况下，即使我们知道当前价格以及下一小时的价格的精确概率分布，我们也无法确定性地计算出价格。这些随机过程，即使在给定初始条件和所有系统参数的情况下，我们仍然无法确定过程的状态，被称为随机随机过程（通常称为过程）。理解或感受随机过程的一种非常好的方式是将其视为与确定性过程的对立面。
马尔可夫过程
如果随机变量在下一时刻的状态仅依赖于当前时刻随机变量的结果，则该随机过程被称为马尔可夫过程。用简单的数学术语表示，对于一个随机过程，S = {R1, R[2], . . ., R[n]} = {R}[t=1, . . ., n]，要成为一个马尔可夫过程，它必须满足以下条件：
根据前面的条件，马尔可夫过程中的任何变量在给定时刻的概率分布是一个条件分布，只依赖于上一时刻的随机变量。系统的这一特性，即系统的未来状态仅依赖于当前状态，也被称为马尔可夫性质。满足马尔可夫性质的系统也被称为无记忆系统，因为它们不需要记住之前的状态来计算下一个状态的分布，换句话说，下一个状态仅依赖于系统的当前状态。
一个常用的例子来解释马尔可夫过程是一个醉汉在街上走路。我们假设，由于这个人醉了，他可以选择向后走一步、向前走一步，或者停留在当前位置，这些动作的概率分布是 [0.4, 0.4, 0.2]。现在，假设我们知道这个人在某一时刻的位置，那么他在下一个时刻的位置仅仅取决于他当前的位置以及系统的参数（他的步长和可能动作的概率分布）。因此，这是一个马尔可夫过程的例子。
在前面的例子中，我们假设醉汉在固定的时间间隔内采取行动（向前/向后走或者停留在原地），并且他的步长始终相同。考虑到这些因素，我们的例子中的马尔可夫过程具有离散的状态空间。此外，由于这个人是在固定的时间间隔后迈出步伐，我们可以把它看作是一个离散时间的过程。但马尔可夫过程不一定需要离散状态空间或离散时间间隔。考虑到离散与连续时间以及离散与连续状态空间，我们可以将马尔可夫过程分为四大类：
-
离散时间与离散状态空间
-
离散时间与连续状态空间
-
连续时间与离散状态空间
-
连续时间与连续状态空间
我们将在接下来的章节中更详细地讨论这些马尔可夫过程的类别。
安装 Python 和包
在继续之前，我们需要设置 Python 和所有运行代码示例所需的包。在本书中的所有代码示例都将使用 Python 3.4。书中的所有示例代码也可以在 GitHub 上找到：github.com/PacktPublishing/HandsOnMarkovModelswithPython。我们强烈建议使用 Miniconda 来设置运行示例所需的环境。Miniconda 可以从 conda.io/miniconda.html 下载。
在 Windows 上的安装
Miniconda 可以通过双击下载的 .exe 文件并按照安装说明进行安装来在 Windows 系统上安装。安装完成后，我们需要创建一个 conda 环境，并在环境中安装所有所需的包。要创建一个名为 hmm 的 Python 3.4 环境，请运行以下命令：
conda create -n hmm python=3.4
在创建环境之后，我们需要激活它并安装所需的包。这可以通过以下命令完成：
activate hmm
conda install numpy scipy
在 Linux 上的安装
在 Linux 上，下载 Miniconda 文件后，我们需要赋予它执行权限，然后进行安装。这可以通过以下命令完成：
chmod +x Miniconda.sh
./Miniconda.sh
执行文件后，我们可以简单地按照安装说明进行操作。一旦安装完成，我们需要创建一个新的环境并安装所需的软件包。我们可以使用以下命令创建一个名为 hmm 的 Python 3.4 环境：
conda create -n hmm python=3.4
环境创建完成后，我们需要激活该环境，并使用以下命令在其中安装软件包：
source activate hmm
conda install numpy scipy
马尔可夫链或离散时间马尔可夫过程
马尔可夫链是一种马尔可夫过程，其中时间是离散的。然而，研究人员对于哪些类别的马尔可夫过程应该被称为马尔可夫链存在较大分歧。但通常来说，它是用来指代离散状态空间的马尔可夫过程。因此，马尔可夫链是一个定义在离散状态空间上的随机过程，满足马尔可夫性质。更正式地说，我们可以说，离散时间马尔可夫链是一个随机变量序列 X[1]、X[2]、X[3]、...，这些随机变量满足马尔可夫性质，即从当前状态转移到下一个状态的概率仅依赖于当前状态，而不依赖于之前的任何状态。就概率分布而言，我们可以说，假设系统在时间点 n 时处于某状态，那么在下一个时间点 n + 1 的条件分布是条件独立于时间点 {1, 2, ..., n-1} 的系统状态，前提是知道时间点 n 的随机变量状态。这可以写成如下形式：
马尔可夫链通常使用有向图表示。有向图中的节点表示随机变量的不同可能状态，边表示系统在下一个时间点从一个状态转移到另一个状态的概率。为了更好地理解这种表示方法，假设我们要预测天气。我们考虑随机变量 Weather={晴天，雨天，雪天} 的三种可能状态，可能的马尔可夫链可以表示如下，如图 1.1所示：
图 1.1：一个简单的马尔可夫链，表示随机变量，表示随机变量 Weather={晴天，雨天，雪天}，并展示了该随机变量在下一个时间点切换到其他状态的概率。
理解马尔可夫链的一个关键点是，我们正在模拟一系列随机变量随时间变化的结果。这一点有时会让人感到困惑，因为模型是通过一个单一的图来表示的，而图中并没有提到任何关于时间的内容。因此，“状态转移”这个名称并不特别合适，因为任何随机变量的状态并没有发生变化；相反，我们是试图根据当前随机变量的观测状态来确定下一个随机变量的状态。回到我们的例子中，可以看到图中的节点表示随机变量天气的不同可能状态，而节点之间的边则显示在当前随机变量状态下，下一个随机变量转移到不同可能状态的概率。自循环则显示模型保持在当前状态的概率。在前面的马尔可夫链中，假设我们知道当前随机变量的观测状态是晴天，那么下一个时间点随机变量仍然是晴天的概率是0.8。它也可能是雨天，其概率为0.19，或者是雪天，其概率为0.01。这里需要注意的是，从任何状态出发的所有外向边的概率值之和应该等于 1，因为这是一个穷尽事件。
现在，让我们尝试编写这个简单的马尔可夫链。我们将首先定义一个简单的MarkovChain类，并在本章中逐步为这个类添加方法：
import numpy as np
class MarkovChain(object):
def __init__(self, transition_prob):
"""
Initialize the MarkovChain instance.
Parameters
----------
transition_prob: dict
A dict object representing the transition probabilities in
Markov Chain. Should be of the form: {'state1': {'state1':
0.1, 'state2': 0.4}, 'state2': {...}}
"""
self.transition_prob = transition_prob
self.states = list(transition_prob.keys())
def next_state(self, current_state):
"""
Returns the state of the random variable at the next time
instance.
Parameters
----------
current_state: str
The current state of the system.
"""
return np.random.choice(
self.states, p=[self.transition_prob[current_state][next_state]
for next_state in self.states])
def generate_states(self, current_state, no=10):
"""
Generates the next states of the system.
Parameters
----------
current_state: str
The state of the current random variable.
no: int
The number of future states to generate.
"""
future_states = []
for i in range(no):
next_state = self.next_state(current_state)
future_states.append(next_state)
current_state = next_state
return future_states
现在，我们可以尝试用这个MarkovChain类来运行我们的例子：
>>> transition_prob = {'Sunny': {'Sunny': 0.8, 'Rainy': 0.19,
'Snowy': 0.01},
'Rainy': {'Sunny': 0.2, 'Rainy': 0.7,
'Snowy': 0.1},
'Snowy': {'Sunny': 0.1, 'Rainy': 0.2,
'Snowy': 0.7}}
>>> weather_chain = MarkovChain(transition_prob=transition_prob)
>>> weather_chain.next_state(current_state='Sunny')
'Sunny'
>>> weather_chain.next_state(current_state='Snowy')
'Snowy'
>>> weather_chain.generate_states(current_state='Snowy', no=10)
['Snowy', 'Snowy', 'Snowy', 'Rainy', 'Snowy', 'Snowy', 'Rainy',
'Rainy', 'Snowy', 'Snowy']
在前面的代码示例中，你可能会发现你的输出与这里显示的不同。这是因为马尔可夫链本质上是概率性的，它基于概率分布选择下一个状态，这可能导致不同运行结果的不同输出。
到目前为止，我们讨论的假设是变量的概率空间在不同时间点之间是不变的。这种情况称为时间齐次马尔可夫链，但也可以存在时间不齐次马尔可夫链，这种情况在许多应用中也有用，但超出了本书的范围。
马尔可夫链的参数化
在上一节的马尔可夫链代码中，我们使用字典对马尔可夫链进行参数化，字典中包含了所有可能的状态转移的概率值。另一种表示状态转移的方式是使用转移矩阵。顾名思义，转移矩阵使用表格形式表示转移概率。图 1.1中的例子对应的转移矩阵如下表所示。
下表展示了图 1.1中马尔可夫链的转移矩阵。概率值表示系统从行中的状态转移到列中状态的概率：
| 状态 | 晴天 | 雨天 | 雪天 |
|---|---|---|---|
| 晴天 | 0.8 | 0.19 | 0.01 |
| 雨天 | 0.2 | 0.7 | 0.1 |
| 雪天 | 0.1+ | 0.2 | 0.7 |
转移矩阵以更紧凑的方式表示与字典相同的信息。因此，转移矩阵是表示马尔可夫链的标准方法。让我们修改我们的MarkovChain类，使其能够接受转移矩阵：
import numpy as np
class MarkovChain(object):
def __init__(self, transition_matrix, states):
"""
Initialize the MarkovChain instance.
Parameters
----------
transition_matrix: 2-D array
A 2-D array representing the probabilities of change of
state in the Markov Chain.
states: 1-D array
An array representing the states of the Markov Chain. It
needs to be in the same order as transition_matrix.
"""
self.transition_matrix = np.atleast_2d(transition_matrix)
self.states = states
self.index_dict = {self.states[index]: index for index in
range(len(self.states))}
self.state_dict = {index: self.states[index] for index in
range(len(self.states))}
def next_state(self, current_state):
"""
Returns the state of the random variable at the next time
instance.
Parameters
----------
current_state: str
The current state of the system.
"""
return np.random.choice(
self.states,
p=self.transition_matrix[self.index_dict[current_state], :])
def generate_states(self, current_state, no=10):
"""
Generates the next states of the system.
Parameters
----------
current_state: str
The state of the current random variable.
no: int
The number of future states to generate.
"""
future_states = []
for i in range(no):
next_state = self.next_state(current_state)
future_states.append(next_state)
current_state = next_state
return future_states
运行这段代码应该会得到与我们在前一节中得到的类似结果。使用转移矩阵可能看起来不是一个好主意，因为它要求我们创建额外的变量来存储索引。但在状态数目达到数百个时，使用转移矩阵比使用简单的字典实现要高效得多。在转移矩阵的情况下，我们可以简单地使用 NumPy 索引在next_state方法中获取概率值，而在前一个实现中，我们是循环遍历所有状态名称的：
>>> transition_matrix = [[0.8, 0.19, 0.01],
[0.2, 0.7, 0.1],
[0.1, 0.2, 0.7]]
>>> weather_chain = MarkovChain(transition_matrix=transition_matrix,
states=['Sunny', 'Rainy', 'Snowy'])
>>> weather_chain.next_state(current_state='Sunny')
'Sunny'
>>> weather_chain.next_state(current_state='Snowy')
'Sunny'
>>> weather_chain.generate_states(current_state='Snowy', no=10)
['Snowy', 'Rainy', 'Rainy', 'Rainy', 'Rainy', 'Rainy',
'Rainy', 'Rainy', 'Sunny', 'Sunny']
马尔可夫链的性质
在本节中，我们将讨论马尔可夫链的不同性质，即可约性、周期性、瞬时性和重现性、遍历性、稳态分析和极限分布。我们还将尝试一些简单的MarkovChain类的示例，以展示这些性质。
可约性
如果我们可以从马尔可夫链中的任何状态到达其他任何状态，则该马尔可夫链被称为不可约的。就状态而言，状态 j 被称为从另一个状态 i 可访问，如果从状态 i 启动的系统有非零的概率到达状态 j。更正式地说，状态 j 被称为从状态 i 可访问，如果存在一个整数 n[ij] ≥ 0，使得满足以下条件：
这里的n[ij] 基本上是从状态 i 到状态 j 所需的步数，它对于不同的 i 和 j 值对可能不同。此外，对于给定的状态 i，如果所有的 n[ij] = 0，这意味着马尔可夫链的所有状态都可以从它直接访问。可访问性关系是自反的和传递的，但不一定是对称的。我们可以通过一个简单的例子来理解这一性质：
图 1.2: 一个不可约马尔可夫链的示例
在之前的例子中，可以清楚地看到所有状态都可以从其他状态访问，因此是不可约的。
请注意，在图 1.2 和 图 1.3 中的例子，我们没有表示概率值为 0 的边。这有助于保持模型不那么复杂，更容易阅读。
在以下示例中，我们可以看到状态D无法从A、B或C访问。此外，状态C无法从A或B访问。但所有状态都可以从状态D访问，并且状态A和B可以从状态C访问：
图 1.3：一个可约马尔可夫链的示例
我们还可以为MarkovChain类添加几个方法，用于检查我们链中的哪些状态是可达的，并且检查我们的链是否是不可约的：
from itertools import combinations
def is_accessible(self, i_state, f_state):
"""
Check if state f_state is accessible from i_state.
Parameters
----------
i_state: str
The state from which the accessibility needs to be checked.
f_state: str
The state to which accessibility needs to be checked.
"""
reachable_states = [i_state]
for state in reachable_states:
if state == self.index_dict[f_state]:
return True
else:
reachable_states.append(np.nonzero(
self.transition_matrix[self.index_dict[i_state], :])[0])
return False
def is_irreducible(self):
"""
Check if the Markov Chain is irreducible.
"""
for (i, j) in combinations(self.states, self.states):
if not self.is_accessible(i, j):
return False
return True
让我们尝试使用图 1.2和图 1.3中的示例来进行实验：
>>> transition_irreducible = [[0.5, 0.5, 0, 0],
[0.25, 0, 0.5, 0.25],
[0.25, 0.5, 0, 0.25],
[0, 0, 0.5, 0.5]]
>>> transition_reducible = [[0.5, 0.5, 0, 0],
[0, 1, 0, 0],
[0.25, 0.5, 0, 0],
[0, 0, 0.25, 0.75]]
>>> markov_irreducible = MarkovChain(transition_matrix=transition_irreducible,
states=['A', 'B', 'C', 'D'])
>>> markov_reducible = MarkovChain(transition_matrix=transition_reducible,
states=['A', 'B', 'C', 'D'])
>>> markov_irreducible.is_accessible(i_state='A', f_state='D')
True
>>> markov_irreducible.is_accessible(i_state='B', f_state='D')
True
>>> markov_irreducible.is_irreducible()
True
>>> markov_reducible.is_accessible(i_state='A', f_state='D')
False
>>> markov_reducible.is_accessible(i_state='D', f_state='A')
True
>>> markov_reducible.is_accessible(i_state='C', f_state='D')
False
>>> markov_reducible.is_irreducible()
False
周期性
如果从状态i返回到状态i的任何可能路径的步数都是k的倍数，则状态i的周期为k。形式化地定义如下：
这里，gcd表示最大公约数（GCD）。基本上，k是从状态i返回到自身的所有可能路径的长度/步数的最大公约数。如果没有从状态i返回到自身的路径，则其周期未定义。我们还需要注意，k与返回起始状态所需的步数无关。例如，假设对于任意给定的状态，返回该状态所需的步数是(4, 6, 8, 12, 16)。在这种情况下，k=2，但返回所需的最小步数是4，而2甚至没有出现在可能的步数列表中。
对于马尔可夫链中的任何给定状态，如果k=1，则该状态称为非周期性。如果马尔可夫链的所有状态都是非周期性的，则该马尔可夫链称为非周期性的。需要注意的是，在不可约马尔可夫链的情况下，单个非周期性状态就足以说明所有状态都是非周期性的。让我们举一个简单的例子，检查不同状态的周期性：
图 1.4：马尔可夫链也是周期性的
在前面的例子中，我们可以很容易地看到，对于状态A，返回的可能路径是A -> B -> C -> A或A -> B -> C -> D -> E -> C -> A。这两条路径的长度分别是 3 和 6，因此状态A的周期是 3。同样，B、C、D和E在马尔可夫链中也都有 3 的周期，因此该马尔可夫链也是周期性的：
图 1.5：具有非周期性状态的马尔可夫链示例。
在这个例子中，我们增加了一些额外的边，因此A的可能路径长度现在是3, 5, 7, ...；B的可能路径长度是2, 3, 4, 5, ...。由于这些路径长度的最大公约数（GCD）为 1，因此A和B现在都是非周期性的。类似地，我们可以计算其他节点的周期，每个节点的周期也是 1，因此该马尔可夫链也是非周期性的。
现在让我们为MarkovChain类添加几个新方法，用于计算不同状态的周期并检查我们的模型是否是非周期性的：
def get_period(self, state):
"""
Returns the period of the state in the Markov Chain.
Parameters
----------
state: str
The state for which the period needs to be computed.
"""
return gcd([len(i) for i in all_possible_paths])
def is_aperiodic(self):
"""
Checks if the Markov Chain is aperiodic.
"""
periods = [self.get_period(state) for state in self.states]
for period in periods:
if period != 1:
return False
return True
现在让我们在我们的例子中尝试这些方法。在这个例子中，我们将随机分配不同过渡的概率值：
>>> transition_periodic = [[0, 1, 0, 0, 0],
[0, 0, 1, 0, 0],
[0.5, 0, 0, 0.5, 0],
[0, 0, 0, 0, 1],
[0, 0, 1, 0, 0]]
>>> transition_aperiodic = [[0, 1, 0, 0, 0],
[0, 0, 1, 0, 0],
[0.5, 0.25, 0, 0.25, 0],
[0, 0, 0, 0, 1],
[0, 0, 0.5, 0.5, 0]]
>>> markov_periodic = MarkovChain(transition_matrix=transition_periodic,
states=['A', 'B', 'C', 'D', 'E'])
>>> markov_aperiodic = MarkovChain(transition_matrix=transition_aperiodic,
states=['A', 'B', 'C', 'D', 'E'])
>>> markov_periodic.get_period('A')
3
>>> markov_periodic.get_period('C')
3
>>> markov_aperiodic.is_periodic()
False
>>> markov_aperiodic.get_period('A')
1
>>> markov_aperiodic.get_period('B')
1
>>> markov_aperiodic.is_periodic()
True
瞬态与重现性
假设我们从状态i开始，如果存在一个非零概率我们永远不会返回到状态i，则称该状态为瞬态。为了更正式地定义这一点，我们可以将随机变量T[i]看作是首次返回到状态i的时间：
现在我们定义另一个术语，
，作为系统在n步之后返回到状态i的概率：
现在我们可以定义，当以下条件满足时，任何给定状态i是瞬态的：
在前面的公式中，我们基本上是在检查返回到状态i的概率总和，步长小于
时是否小于1。如果总和小于1，则意味着T[i]返回到
的概率大于0，这就意味着状态i是瞬态的。如果状态i不是瞬态的，那么它被称为常返：
图 1.6：
在上面的例子中，我们可以看到状态A和B是瞬态的，因为A没有任何输入边。B确实有输入边，但它来自另一个瞬态状态，因此它也是瞬态的。因此，一旦系统离开状态A或B，就无法再返回。
检查一个给定状态是否是瞬态状态其实非常简单。我们只需检查是否有来自其他状态的输入边。如果没有，那么该状态就是瞬态的。让我们为MarkovChain类编写一个简单的方法来检查这一点：
def is_transient(self, state):
"""
Checks if a state is transient or not.
Parameters
----------
state: str
The state for which the transient property needs to be checked.
"""
if all(self.transition_matrix[~self.index_dict[state], self.index_dict[state]] == 0):
return True
else:
return False
现在我们可以在图 1.6中的示例中使用此方法来检查哪些节点是瞬态的：
>>> transient_matrix = [[0, 0.5, 0.5, 0],
[0, 0, 0.25, 0.75],
[0, 0, 0, 1],
[0, 0, 0.5, 0.5]]
>>> transient_markov = MarkovChain(transition_matrix=transient_matrix,
states=['A', 'B', 'C', 'D'])
>>> transient_markov.is_transient('A')
True
>>> transient_markov.is_transient('B')
True
>>> transient_markov.is_transient('C')
False
在接下来的小节中，我们将讨论随机变量T[i]的统计性质。
平均返回时间
初始状态i的首次返回时间也被称为击中时间。在上一节中，它通过随机变量T[i]表示。状态i的平均返回时间定义为其预期返回时间：
如果平均返回时间M[i]是有限的，则该状态被称为正常常返，否则称为无常返。
预期访问次数
显然，从名称上看，任何状态i的预期访问次数是系统预计将处于该状态的次数。此外，只有当状态i的预期访问次数为无限时，该状态才是常返的：
吸收状态
如果一旦系统达到某个状态后无法离开该状态，则称状态i为吸收状态。要成为吸收状态，停留在同一状态的概率应该是1，而所有其他状态的概率应该是0：
在马尔可夫链中，如果所有状态都是吸收状态，则我们称其为吸收马尔可夫链：
图 1.7：一个示例，展示了吸收状态 C，因为从状态 C 转移到 C 的概率是1
再次，我们可以为我们的MarkovChain类添加一个非常简单的方法来检查吸收状态：
def is_absorbing(self, state):
"""
Checks if the given state is absorbing.
Parameters
----------
state: str
The state for which we need to check whether it's absorbing
or not.
"""
state_index = self.index_dict[state]
if self.transition_matrix[state_index, state_index]
我们可以通过创建一个马尔可夫链并使用is_absorbing方法再次检查示例中的状态是否是吸收状态：
>>> absorbing_matrix = [[0, 1, 0],
[0.5, 0, 0.5],
[0, 0, 1]]
>>> absorbing_chain = MarkovChain(transition_matrix=absorbing_matrix,
states=['A', 'B', 'C'])
>>> absorbing_chain.is_absorbing('A')
False
>>> absorbing_chain.is_absorbing('C')
True
遗传性
状态i被称为具有遗传性的，如果它是重现的，周期为1，并且具有有限的平均重现时间。如果一个马尔可夫链的所有状态都是遗传的，那么它就是一个遗传马尔可夫链。一般而言，如果存在一个数字N，使得系统中的任何状态都可以通过大于或等于N步数的方式从任何其他状态到达，那么这个马尔可夫链就是遗传的。因此，在完全连接的转移矩阵中，所有转移具有非零概率时，这个条件通过N=1得以满足。
稳态分析与极限分布
在马尔可夫链中，如果向量π满足∀ j ∈ s，并且满足以下条件，则称其为定常分布：
定常分布是马尔可夫链最重要的性质之一，我们将在本章后续部分详细讨论它。
连续时间马尔可夫链
连续时间马尔可夫链与离散时间马尔可夫链非常相似，不同之处在于，在连续时间的情况下，我们显式地使用一个正值随机变量来建模状态之间的转换时间。此外，我们还会考虑系统在所有可能的时间值下的状态，而不仅仅是转换时间。
指数分布
随机变量x被称为具有分布率λ的指数分布，如果其概率密度函数定义如下：
这里，分布率λ需要大于0。我们还可以按如下方式计算X的期望：
我们看到 X 的期望与学习率成反比。这意味着，具有较高学习率的指数分布会有较低的期望值。指数分布通常用于建模涉及某个事件发生时间的问题。一个简单的例子是建模闹钟响铃之前的时间，或者餐厅服务员上桌前的时间。而且，正如我们所知 
，学习率越高，我们期望事件发生的时间越短，因此得名 学习率。
我们还可以计算指数分布的第二矩和方差：
使用第一矩和第二矩，我们可以计算该分布的方差：
图 1.x：指数分布的概率分布
现在我们将继续讨论一些与我们示例相关的指数分布的属性：
-
无记忆性：图 1.x 显示了指数分布的图像。在图中，我们可以清晰地看到，在任意给定点（此例为 a）之后的图形与原始分布完全相同。我们也可以说，条件为 (X > a) 的指数分布仍然是指数分布。如果从我们的示例来理解这一属性，这意味着如果我们有一个闹钟，在任意时间 t，我们检查它是否还没有响铃，我们仍然可以确定 t 之后的时间的分布，这将是相同的指数分布。指数分布的这一属性被称为 无记忆性，因为在任何给定的时间点，如果你知道系统的当前状态（例如闹钟没有响），你就可以确定未来时间的概率分布。指数分布的这一属性与马尔科夫链非常相似，正如你可能在前面章节中回忆到的。
-
最小值的概率：假设我们有 n 个独立的指数分布，随机变量分别为 X[0]，...，X[n]，学习率分别为 λ[0]，...，λ[n]。对于这些分布，我们可以证明，min(X[0], . . ., X[n]) 的分布也是一个指数分布，学习率为
![]()
：
：
我们将在后面的连续时间马尔科夫链示例中使用指数分布的这两个属性。
泊松过程
泊松过程是一个连续过程，可以有多种解释，从而导致不同的定义。在本节中，我们将从正式定义开始，逐步过渡到更简单、更直观的定义。一个连续时间随机过程N(t): t > 0，如果满足以下条件，则称其为泊松过程，速率为λ > 0：
-
N(0) = 0
-
它具有平稳和独立增量
-
N(t)的分布为泊松分布，均值为λt：
首先，我们需要定义什么是平稳增量和独立增量。对于一个连续时间随机过程X(t): t ≥ 0，增量被定义为系统在两个时间点之间状态的差异；也就是说，给定两个时间点s和t，且s < t，从时间s到时间t的增量是X(t) - X(s)。顾名思义，一个过程如果其增量的分布仅依赖于时间差，则称为具有平稳增量的过程。
换句话说，如果一个过程的增量分布依赖于时间差，而与具体时间点无关，则该过程称为具有平稳增量；即，如果X(t[1]) - X(s[1])的分布等于X(t[2]) - X(s[2])，且满足t[1] > s[1], t[2] > s[2]并且t[1] - s[1] = t[2] - s[2]，则该过程具有平稳增量。如果任意两个不重叠时间区间的增量是独立的，则该过程称为具有独立增量；也就是说，如果t[1] > s[1] > t[2] > s[2]，那么增量X(t[2]) - X(s[2])和X(t1) - X(s1)是独立的。
现在我们回到定义泊松过程。泊松过程本质上是一个计数过程，用来统计在时间t之前发生的事件数量。这个时间t之前发生的事件数量由N(t)给出，同样，在时间区间t到t + s之间发生的事件数量由N(t + s) - N(t)给出。N(t + s) - N(t)的值服从泊松分布，均值为λ[s]。我们可以看到，泊松过程在固定时间间隔内具有平稳增量，但如同 
，N(t)的值也将趋近于无穷大；也就是说，
。另一个值得注意的点是，随着λ值的增加，发生的事件数量也会增加，这就是为什么λ也被称为过程的速率。
这引出了我们对泊松过程的第二个简化定义。一个连续时间随机过程N(t): t ≥ 0，如果满足以下条件，则称其为泊松过程，学习速率为λ > 0：
-
N(0) = 0
-
它是一个计数过程；也就是说，N(T)给出了在时间t之前发生的事件数量。
-
事件之间的时间间隔是独立且同分布的，服从指数分布，学习速率为λ。
连续时间马尔科夫链示例
现在，既然我们对指数分布和泊松过程有了基本的理解，我们可以进入示例，构建一个连续时间马尔可夫链。在这个例子中，我们将展示如何利用指数分布的特性来构建通用的连续时间马尔可夫链。假设有一个酒店接待处，n个接待员并行工作。同时，假设客人按照泊松过程到达，速率为λ，每位客人的服务时间使用具有学习速率µ的指数随机变量表示。此外，如果所有接待员都忙碌，当新客人到达时，他/她将无法获得服务并离开。现在假设一个新客人到达并发现所有接待员都忙碌，接下来我们尝试计算在下一个时间区间内的忙碌接待员的期望数量。
假设T[k]表示下一个时间点上有k个忙碌的接待员。我们也可以使用T[k]表示如果当前时刻有k个接待员忙碌，下一个到达的客人找到的忙碌接待员的期望数量。这两种T[k]的表示是等效的，因为指数分布具有无记忆性。
首先，T[0]显然是0，因为如果当前没有忙碌的接待员，下一个到达的客人也一定会发现没有忙碌的接待员。现在考虑T[1]，如果当前有i个忙碌的接待员，下一个到达的客人会发现有 1 个忙碌的接待员，前提是下一次到达的时间小于忙碌接待员剩余的服务时间。从无记忆性属性出发，我们知道下一次到达时间服从指数分布，学习速率为λ，而剩余的服务时间也服从指数分布，学习速率为µ。因此，下一个客人发现一个接待员忙碌的概率是
，因此以下结论成立：
一般来说，我们考虑 k 个接待员忙碌的情况。我们可以通过条件概率获得 T[k] 的表达式，具体条件取决于哪个事件先发生。当我们有 k 个接待员忙碌时，我们可以认为有 k+1 个独立的指数分布：k 个指数分布，代表每个接待员的剩余服务时间，学习速率为 µ，以及 1 个指数分布，代表下一位到达的客人，学习速率为 λ。在我们的案例中，我们希望根据服务完成还是新客人到达哪个事件先发生来进行条件化。服务完成的时间将是 k 个指数分布中的最小值。这个完成时间同样遵循指数分布，学习速率为 kµ。现在，服务完成发生在下一位客人到达之前的概率为 
。类似地，下一件发生的事件是客人到达的概率为 
。
现在，基于这一点，我们可以说，如果下一个事件是服务完成，那么预期的忙碌接待员数量将是 T[k-1]。否则，如果客人先到，将会有 k 个忙碌的接待员。因此我们得出以下结论：
我们现在只需要解这个递归方程。 T[2] 将由以下方程给出：
如果我们继续这个模式，我们将得到 T[3]，其形式如下：
我们可以看到 T[1] 和 T[2] 的值中有规律可循，因此我们可以写出一个通项公式，如下所示：
让我们总结一下之前例子的观察结果：
-
在任何给定的时间点，如果有 i 个忙碌的接待员，且 i < n，则有 i + 1 个独立的指数分布，其中 i 个的速率为 µ，而 1 个的速率为 λ。直到过程发生跳跃的时间是指数分布，其速率为 iµ + λ。如果所有接待员都忙碌，那么只有与服务时间相关的 n 个指数分布能够触发跳跃，且过程发生跳跃的时间是指数分布，速率为 nµ。
-
当过程从状态 i 跳转时，若 i < n，则以概率
![]()
跳转到状态 i + 1，以概率 ![]()
跳转到状态 i - 1。 -
当过程从状态 i 发生跳跃时，我们可以启动一整套新的分布，分别对应我们跳转到的状态。这是因为，即使某些旧的指数分布尚未触发，也等同于重置或替换这些分布。
跳转到状态 i + 1，以概率 
跳转到状态 i - 1。每次我们跳转到状态 i 时，无论时间何时，停留在状态 i 中的时间分布以及离开状态 i 后跳转到下一个状态的概率都是相同的。换句话说，该过程是时间齐次的。
前面对连续时间随机过程的描述对应于一个连续时间马尔可夫链。在接下来的部分，我们将尝试以更正式的方式来定义它。
连续时间马尔可夫链
在前一部分中，我们展示了一个连续时间马尔可夫链的例子，以说明它是如何工作的。现在我们将继续正式定义它。在一个具有离散状态空间 S 的连续时间马尔可夫链中，对于每个状态 i ∈ S，我们有一组 n[i] 个独立的指数分布，其速率为 q[i], j[1], ..., q[i], j[n[i]]，其中 j[1], ..., j[n[i]] 是该过程离开状态 i 时可能跳转到的状态集合。当该过程进入状态 i 时，它在状态 i 中停留的时间服从速率为 v[i] = q[i]j[1] + ... + q[i]j[n[i]] 的指数分布，而当它离开状态 i 时，它将以概率 
跳转到状态 j[l]，其中 l = 1, ..., n[i]。
我们还可以将马尔可夫性从离散时间情况扩展到连续时间情况。
对于一个状态空间为 S 的连续时间随机过程 (X(t) : t ≥ 0)，我们说它具有马尔可夫性质，如果满足以下条件：
这里，0 ≤ t[1] ≤ t2 ≤ ... ≤ t[n-1] ≤ s ≤ t 是任何非递减的 n + 1 个时间点的序列，i[1, i]2, ..., i[n-1], i, j ∈ S 是状态空间中任意的 n + 1 个状态，对于任意整数 n ≥ 1。
类似地，我们可以将时间齐次性扩展到连续时间马尔可夫链的情况。我们说一个连续时间马尔可夫链是时间齐次的，如果对于任意 s ≤ t 和任意状态 i, j ∈ S，满足以下条件：
与离散时间马尔可夫链的情况一样，连续时间马尔可夫链不需要是时间齐次的，但非齐次的马尔可夫链超出了本书的讨论范围。关于非齐次马尔可夫链的更多细节，您可以参考黄成吉关于这一主题的论文：lib.dr.iastate.edu/cgi/viewcontent.cgi?article=8613&context=rtd。
现在让我们定义连续时间马尔可夫链的转移概率。就像连续时间马尔可夫链中的速率 q[ij] 是离散时间马尔可夫链中转移概率 p[ij] 的对应物一样，时间齐次的连续时间马尔可夫链也有与 n 步转移概率 pij 对应的概念，其定义如下：
摘要
在本章中，我们对马尔可夫链进行了详细介绍。我们讨论了不同类型的马尔可夫链，主要是具有离散状态空间的链，涵盖离散时间或连续时间的情况。我们还介绍了时间齐次和非时间齐次马尔可夫链的概念。我们详细讨论了马尔可夫链的不同属性，并提供了相关的示例和代码。
马尔可夫链及其属性是隐马尔可夫模型（HMM）构建的基础概念。在下一章，我们将更加详细地讨论 HMM。
第二章：隐马尔可夫模型
在上一章中，我们讨论了马尔可夫链，它对于建模跨时间的观察序列非常有帮助。在本章中，我们将研究隐马尔可夫模型（HMM），它也用于建模序列数据，但比马尔可夫链更加灵活。
在本章中，我们将涵盖以下主题：
-
马尔可夫模型
-
HMM
-
HMM 的评估
-
HMM 的扩展
马尔可夫模型
马尔可夫模型是一种随机模型，其中随机变量在下一个时间点的状态仅依赖于当前时间点的随机变量结果。最简单的马尔可夫模型是马尔可夫链，我们在第一章中讨论过，马尔可夫过程简介。
假设我们有一组顺序观察值（x[1],. . ., x[n])，遵循马尔可夫性质，那么我们可以表示N个观察值的联合概率分布如下：
一阶马尔可夫链的图示表示，其中当前观察值的分布依赖于前一个观察值的值
上述马尔可夫链的表示与我们之前看到的表示不同。在这种表示中，观察值呈现为节点，边表示两个观察值之间的条件概率，即Pr(x[n]|x[n-1])。这就是概率图模型的一般表示方式，其中节点表示随机变量，边表示这两个变量之间的条件概率分布。这种图形化表示帮助我们理解随机变量之间的因果关系。
状态空间模型
我们可以看到，简单的马尔可夫链是非常有限的，不能很好地处理我们预期多个连续观察值能提供预测下一个观察值所需的重要信息的情况。幸运的是，马尔可夫链可以进行调整，以支持这些情况：
二阶马尔可夫链，其中当前观察值的分布依赖于最后两个观察值的值
让我们考虑一个马尔可夫链，其中下一个状态的概率不仅依赖于当前状态，还依赖于前一个状态。这种类型的马尔可夫链被称为二阶马尔可夫链，其联合概率分布可以表示如下：
利用 d-分离性质，我们可以看到，给定X[n-1]和X[n-2]的条件下，X[n]的条件分布与所有观察值X[1], . . ., X[n-3]是独立的。
类似地，我们可以将其扩展到一个M^(th)阶马尔可夫链，其中特定观察值的条件分布依赖于前M个观察值。然而，现在我们需要为增加的灵活性付出代价，即需要更多的参数。
假设观察值是具有K状态的离散变量。那么在一阶马尔可夫链中，条件分布Pr(X[n]|X[n-1])将由每个K状态的X[n-1]的K-1个参数来指定，给出总共K(K-1)个参数。如果我们将其扩展到一个M^(th)阶马尔可夫链，其中联合分布由条件分布Pr(x[n]|x[n-M], . . ., x[n-1])构建，那么该模型中的参数数量将为K^(M-1)(K-1)。由于这个数量随着M的增加呈指数增长，它通常会使得对于较大的M值，这种方法变得不切实际。但如果我们想构建一个不受马尔可夫假设阶数限制的序列数据模型，并且它可以通过有限数量的参数来表示，该怎么办呢？
这种模型可以通过引入潜在（或隐藏）变量来创建。潜在变量使我们能够构建由简单组件组成的丰富模型类。假设对于每个观察值，x[n]，我们有一个潜在变量，z[n]（它可能与观察变量的维度相同，也可能不同）。如果这些潜在变量形成一个一阶马尔可夫链，那么这种类型的模型可以称为状态空间模型，其表示如下图所示：
表示分布的状态空间模型，其中每个观察值都基于一个潜在变量进行条件化
利用 d-separation 性质，我们可以看到，任何两个观察变量之间始终存在一条路径，x[n]和x[m]通过潜在变量连接，而且这条路径永远不会被阻塞。因此，给定所有前期观察值的Pr(x[n+1]|x[1], . . ., x[n])分布不具有任何条件独立性特性，因此，x[n+1]的预测依赖于所有前期观察值。
由于潜在变量形成了一个马尔可夫链，它们满足以下条件分布：
因此，这个模型的联合分布可以表述如下：
隐马尔可夫模型（HMM）
隐马尔可夫模型（HMM）是状态空间模型的一个特例，其中潜在变量是离散的多项式变量。从图形表示来看，我们也可以将 HMM 视为一个双重随机过程，包含一个我们无法直接观察的隐藏随机马尔可夫过程（潜在变量），以及另一个生成观察序列的随机过程，这个过程基于第一个过程。
在继续讨论参数化之前，我们先通过一个掷硬币的例子来了解它的工作原理。假设我们有两枚不公平的硬币，M[1] 和 M[2]，其中 M[1] 具有较高的正面概率（70%），而 M[2] 具有较高的反面概率（80%）。某人按顺序掷这两枚硬币，但我们并不知道是哪一枚。我们只能观察到结果，它可能是正面（H）或反面（T）：
我们可以将选中的不公平硬币视为潜在变量，硬币投掷的结果视为观察数据。为了预测下一次观察序列的结果，我们至少需要以下信息：首先选择了哪枚硬币、给定之前选中的硬币下一枚硬币的选择以及根据硬币得到 H 或 T 的概率。假设两枚硬币在首次选择时的优先级相等，并且每次硬币的选择都是根据之前选择的硬币以同等的概率进行选择，我们可以创建如下的状态图：
掷硬币 HMM 的状态图
在之前的图示中，z[1] 和 z[2] 代表潜在变量硬币的状态（其中 z[1] 代表硬币 M[1] 被选中，z[2] 代表硬币 M[2] 被选中）。弧线表示从一个潜在变量状态转移到另一个状态的转移概率，直线则表示给定潜在变量（选中的硬币）时观察到的变量（投掷结果）的概率。
HMM 的参数化
在上一节中，我们看到了一个 HMM 的例子来了解模型的工作原理。现在我们正式开始对 HMM 进行参数化。
由于 HMM 的潜在变量是离散的多项式变量，我们可以使用 1-of-K 编码方案来表示它，其中 z[n] 变量由一个 K 维的二进制变量向量表示，z[nk] ∈ {0,1}，当 z[n] 变量处于 k 状态时，z[nk] = 1 且 z[nj] = 0（对于 j ≠ k）。
考虑到这一点，我们可以创建一个包含转移概率矩阵 A 的矩阵，其中 A[ij] = Pr(Z[nj] = 1| z[n-1], i = 1)。由于 A[ij] 代表从状态 i 转移到状态 j 的概率，它具备 
的性质，并且可以通过 K(K-1) 参数表示。因此，我们可以如下表示条件概率分布：
转移矩阵通常通过状态转移图来表示，正如我们在第一章中看到的，马尔可夫过程简介。但我们可以将相同的表示方法展开到时间维度，得到一个网格或格子图，如下图所示。在接下来的章节中，我们将使用这种 HMM 的表示方法来学习参数并从模型中进行推断：
带有三个状态的 HMM 隐变量的格子图
由于初始隐节点z[1]没有父节点，它具有边际分布Pr(z[1])，可以通过概率向量π表示，π[k] = Pr(z[1k] = 1)，如图所示。这样，Pr(z[1]|π)*的概率可以表示如下：
参数化 HMM 所需的第三个也是最后一个参数是给定隐变量的观测变量的条件概率，即发射概率。它由条件分布Pr(x[n]| z[n], Φ)表示，且由一些参数Φ控制。如果观测变量x[n]是离散的，则发射概率可以采用条件概率表的形式（多项式 HMM）。类似地，如果观测变量x[n]是连续的，则该分布可能是高斯分布（高斯 HMM），其中
表示控制分布的参数集，即均值和方差。
因此，隐变量和观测变量的联合概率分布可以表示如下：
这里，X = {x[1], ..., x[N]}，Z = {z[1], ..., z[N]}，以及θ = {A, π, Φ}表示控制模型的参数集。
当所有控制隐变量的条件分布共享相同的转移矩阵A，且所有发射概率共享相同的参数Φ时，HMM 模型被称为齐次模型：
现在，让我们尝试编写一个简单的多项式隐马尔可夫模型（HMM）。我们将从定义一个简单的MultinomialHMM类开始，并在接下来的过程中不断添加方法：
import numpy as np
class MultinomialHMM:
def __init__(self, num_states, observation_states, prior_probabilities,
transition_matrix, emission_probabilities):
"""
Initialize Hidden Markov Model
Parameters
-----------
num_states: int
Number of states of latent variable
observation_states: 1-D array
An array representing the set of all observations
prior_probabilities: 1-D array
An array representing the prior probabilities of all the states
of latent variable
transition_matrix: 2-D array
A matrix representing the transition probabilities of change of
state of latent variable
emission_probabilities: 2-D array
A matrix representing the probability of a given observation
given the state of the latent variable
"""
# As latent variables form a Markov chain, we can use
# use the previous defined MarkovChain class to create it
self.latent_variable_markov_chain = MarkovChain(
transition_matrix=transition_matrix,
states=['z{index}'.format(index=index) for index in range(num_states)],
)
self.observation_states = observation_states
self.prior_probabilities = np.atleast_1d(prior_probabilities)
self.transition_matrix = np.atleast_2d(transition_matrix)
self.emission_probabilities = np.atleast_2d(emission_probabilities)
使用MultinomialHMM类，我们可以像之前讨论的那样定义 HMM 硬币模型，如下所示：
coin_hmm = MultinomialHMM(num_states=2,
observation_states=['H', 'T'],
prior_probabilities=[0.5, 0.5],
transition_matrix=[[0.5, 0.5], [0.5, 0.5]],
emission_probabilities=[[0.8, 0.2], [0.3, 0.7]])
生成观测序列
对于一个由{A, π, Φ}参数化的 HMM，我们可以通过以下步骤生成一系列观测值，{x[1], ..., x[N]}：
-
设置n = 1
-
根据先验分布π选择初始隐状态z[1]：
-
根据给定的z[1]值，通过采样由Φ控制的发射概率分布来选择观测值，x[1]：
-
根据状态转移概率矩阵 A，转移到潜变量 z[n+1] 的下一个状态
-
设置 n = n + 1 并重复第 3 步，如果 n ≤ N，否则终止
我们可以向之前定义的 MultinomialHMM 类中添加一个生成样本的方法，如下所示：
def observation_from_state(self, state):
"""
Generate observation for a given state in accordance with
the emission probabilities
Parameters
----------
state: int
Index of the current state
"""
state_index = self.latent_variable_markov_chain.index_dict[state]
return np.random.choice(self.observation_states,
p=self.emission_probabilities[state_index, :])
def generate_samples(self, no=10):
"""
Generate samples from the hidden Markov model
Parameters
----------
no: int
Number of samples to be drawn
Returns
-------
observations: 1-D array
An array of sequence of observations
state_sequence: 1-D array
An array of sequence of states
"""
observations = []
state_sequence = []
initial_state = np.random.choice(self.latent_variable_markov_chain.states,
p=self.prior_probabilities)
state_sequence.append(initial_state)
observations.append(self.observation_from_state(initial_state))
current_state = initial_state
for i in range(2, no):
next_state = self.latent_variable_markov_chain.next_state(current_state)
state_sequence.append(next_state)
observations.append(self.observation_from_state(next_state))
current_state = next_state
return observations, state_sequence
我们可以在 HMM 投币示例上使用 generate_samples 方法来生成观测序列：
>>> coin_hmm.generate_samples()
(['T', 'H', 'H', 'T', 'T', 'H', 'H', 'H', 'H'], ['z1', 'z0', 'z0', 'z1', 'z1', 'z0', 'z1', 'z1', 'z1'])
安装 Python 包
HMM 也可以拥有用于发射概率的高斯分布。就像 MultinomialHMM 一样，我们也可以从 GaussianHMM 中进行采样。在下一个代码示例中，我们使用 hmmlearn 库提供的 GaussianHMM 类，看看如何从这种类型的模型中生成样本：
source activate hmm
conda install scikit-learn
pip install hmmlearn
一旦安装了 Python 包，我们可以使用以下代码从 Gaussian HMM 生成样本：
from hmmlearn.hmm import GaussianHMM
import numpy as np
import matplotlib.pyplot as plt
startprob = np.array([0.6, 0.3, 0.1, 0.0])
# The transition matrix, note that there are no transitions possible
# between component 1 and 3
transmat = np.array([[0.7, 0.2, 0.0, 0.1],
[0.3, 0.5, 0.2, 0.0],
[0.0, 0.3, 0.5, 0.2],
[0.2, 0.0, 0.2, 0.6]])
# The means of each component
means = np.array([[0.0, 0.0],
[0.0, 11.0],
[9.0, 10.0],
[11.0, -1.0]])
# The covariance of each component
covars = .5 * np.tile(np.identity(2), (4, 1, 1))
# Build an HMM instance and set parameters
model = hmm.GaussianHMM(n_components=4, covariance_type="full")
# Instead of fitting it from the data, we directly set the estimated
# parameters, the means and covariance of the components
model.startprob_ = startprob
model.transmat_ = transmat
model.means_ = means
model.covars_ = covars
X, state_sequence = model.sample(n_samples=100)
plt.plot(X[:, 0], X[:, 1], ".-", label="observations", ms=6,
mfc="orange", alpha=0.7)
for i, m in enumerate(means):
plt.text(m[0], m[1], 'Component %i' % (i + 1),
size=12, horizontalalignment='center',
bbox=dict(alpha=.7, facecolor='w'))
plt.legend(loc='best')
plt.show()
从具有四状态潜变量 z 和高斯发射模型 p(x|z) 的 HMM 中进行采样
HMM 的评估
在前一节中，我们讨论了生成给定 HMM 的观测序列。但是，实际上，大多数时候我们并不关心生成观测序列，主要是因为我们一开始并不知道 HMM 的参数，因此无法生成观测。
对于给定的 HMM 表示，在大多数应用中，我们总是试图解决以下三个问题：
-
模型评估：给定模型的参数和观测序列，估计该序列的概率
-
预测最优序列：给定模型的参数和观测序列，估计产生这些观测的最可能的状态序列
-
参数学习：给定一系列观测，估计生成这些观测的 HMM 模型的参数
在本节中，我们将讨论第一个问题，即模型评估，接下来的章节中，我们将详细讨论其他两个问题。正如我们将在后续章节中看到的，评估模型是解决其他两个问题的基础。因此，高效地解决这个问题是迈向参数学习和推理的第一步。
让我们正式描述模型评估问题。给定一个由 θ = {A, π, Φ} 参数化的 HMM 和一个观测序列，X = {x[1], ..., x[N]}，我们需要计算 Pr(X|θ) 的概率。通过前一节的讨论，我们可以说，我们可以通过边缘化联合概率分布 Pr(X, Z|θ) 来计算 Pr(X|θ)，其中 Z = {z[1], ..., z[N]}，对 Z 进行边际化：
在前一节中，我们看到以下结论：
因此，Pr(X, Z|θ) 概率可以表示为如下：
对于具有 K 状态和观测长度 N 的模型，有 K^T 种可能的状态序列。求和中的每一项需要 2N 次运算。因此，评估的复杂度为 2N X K^T 的数量级。例如，如果我们考虑一个具有五个状态的模型，K = 5，并且观测序列的长度为 N = 100，则所需的运算次数为 10⁷² 级别，这使得即使对于一个非常小的 HMM，评估方法也变得不可处理。
HMM 的扩展
在前面的章节中，我们讨论了 HMM，从中进行采样并根据其参数评估给定序列的概率。在本节中，我们将讨论它的一些变种。
因式分解 HMM
让我们考虑建模一系列图像中多个物体的问题。如果图像中有 M 个物体，且每个物体有 K 种不同的位置和方向，则该图像背后的系统可能有 K^M 种状态。一个 HMM 将需要 K^M 个不同的状态来建模该系统。这种表示系统的方式不仅效率低，而且难以解释。我们更希望我们的 HMM 能够通过使用 M 个不同的 K 维变量来捕捉状态空间。
因式分解 HMM 就是这样的一种表示方法。在这个模型中，存在多个独立的马尔可夫链潜在变量，并且在任何给定时间，观测变量的分布都依赖于该时间点所有对应潜在变量的状态。该系统的图形模型可以表示如下：
包含两个马尔可夫链的因式分解 HMM
考虑因式分解 HMM 的动机可以通过注意到，在某一时间步表示 10 位信息时，标准 HMM 需要 K = 2¹⁰ = 1024 个潜在状态，而因式分解 HMM 可以使用 10 个二进制的潜在链。然而，这会带来额外的训练复杂度，正如我们将在后续章节中看到的那样。
树形结构 HMM
在前一节中，我们讨论了因式分解 HMM，其中潜在变量形成独立的马尔可夫链。通过引入潜在变量之间的耦合，可以放宽关于潜在变量的独立性假设。一种耦合潜在变量的方法是将它们排序，使得 
依赖于 
对于所有 1 ≤ l ≤ m。此外，如果所有的输出变量和潜在变量都依赖于某个随机输入变量 x[n]，则我们得到一个树形结构的 HMM，表示如下：
树形结构 HMM 的图示表示
该模型的架构可以解释为一个概率决策树。我们首先考虑时间切片n=1，并尝试理解这个模型是如何生成数据的。根据x[1]的值，顶节点，
，可以取K个值（假设隐藏变量有K个状态）。这将x空间划分为K个决策组。下一个节点，
，进一步划分为K个子区域，以此类推。输出，y[1]，是由输入x[1]和每个节点的 K 种决策生成的。在下一个时间切片中，同样的事情会发生，只不过树中的每个决策都依赖于前一个时间切片中该节点所做的决策。因此，该模型可以解释为一个具有马尔可夫动态的概率决策树。
总结
本章中，我们详细介绍了马尔可夫模型和 HMM。我们讨论了如何对 HMM 进行参数化，从中生成样本，以及它们的代码。我们讨论了如何估计观测的概率，这将构成推理的基础，下一章我们将讲解这一部分内容。我们还讨论了 HMM 的各种扩展。
在下一章中，我们将深入探讨 HMM 中的推理。
第三章：状态推理 - 预测状态
在前几章中，我们介绍了马尔可夫链和隐马尔可夫模型（HMM），并展示了如何用它们来建模问题。在本章中，我们将看到如何使用这些模型进行预测或向模型提问（称为推理）。用于计算这些值的算法称为推理算法。在本章中，我们将专门研究如何计算状态变量的概率分布。
本章将涵盖以下主题：
-
隐马尔可夫模型中的状态推理
-
动态规划
-
前向-后向算法
-
维特比算法
隐马尔可夫模型中的状态推理
让我们从一个简单的例子开始，展示我们可以对我们的 HMM 模型提出什么有趣的问题。我们以机器人定位为例。这个例子有很多变体，但我们假设机器人在一个二维网格中移动，如图 3.1所示。机器人上有四个传感器，这些传感器能检测机器人周围是否有墙壁，并且每个传感器都有一个特定的方向。
我们希望在以下网格中建模机器人的运动，并结合传感器的观测结果：
图 3.1：机器人在不同时间点的位置概率分布
在图 3.1中，我们可以看到不同时间点的观察如何改变机器人在网格中位置的概率。最初，我们从对网格中所有位置的均匀概率开始。现在，在时间t=1时，假设我们机器人的传感器显示出顶部和底部有墙壁。得到这个观察结果后，我们对机器人位置的感知会发生变化。现在，我们会更倾向于认为机器人处于 2、4、10 或 12 块区域，如图 3.1中的第二块所示。如果我们在时间实例t=2，机器人传感器显示只有顶部有墙壁，那么它最有可能处于 3 块区域，如图 3.1中的第三块所示。这是因为我们知道它上次最有可能的位置，并且将这个信息与当前的传感器读数结合起来，3 块区域就是机器人可能的位置。现在，如果我们在时间t=3，机器人传感器显示左右两侧都有墙壁，那么它最有可能位于 7 块区域。这一过程使我们能够仅凭传感器读数随着时间推移定位机器人在网格中的位置。
由于我们在模型中考虑了机器人在一段时间内的状态转移（网格中的位置）以及每个时刻的结果（传感器输出），因此 HMM 在这种情况下似乎是完美的模型。在这个示例中，我们假设我们知道机器人的位置转移概率。我们还假设我们知道网格的结构，因此我们也能知道发射概率。我们可以利用发射概率来建模传感器输出的不确定性，因为它们在某些时刻可能会给出噪声结果：
图 3.2：机器人定位的一个 HMM 示例
现在让我们考虑一下我们可能想要向模型提问的问题。我们可能有兴趣了解在某个时刻给定直到该时刻的所有观测值后，机器人的位置。我们可能还想知道在某个时刻，给定直到该时刻的所有位置，传感器输出的概率。我们也可能有兴趣计算我们观察到的变量和机器人位置的联合分布。所有这些值都可以使用 前向算法、后向算法或前向-后向算法 轻松计算。
现在，我们不再询问分布，而是可能对机器人走过的最可能路径感兴趣。为了计算最可能的路径，我们需要对机器人每个时刻的状态进行 MAP 推理。使用维特比算法可以高效地完成这项工作。
在接下来的章节中，我们将正式介绍这些算法，并看看如何实现它们。所有这些算法都依赖于一个非常重要的编程范式，称为动态规划。动态规划使我们能够在可处理的时间内运行这些 HMM 推理算法。我们将在下一节详细讨论动态规划。
动态规划
动态规划是一种编程范式，在这种范式中，我们将一个复杂问题划分为更小的子问题。我们解决这些子问题并存储结果。每当我们需要重新计算相同的子问题时，我们只需使用存储的结果，从而节省了计算时间，代价是使用了存储空间。这种缓存子问题结果的技术称为备忘录法。因此，使用动态规划通过备忘录法可以加速我们的计算，并且在某些情况下，它能将计算复杂度从指数级降低到线性级，正如我们将在以下示例中看到的那样。
使用动态规划优化的最简单示例之一是计算斐波那契数列的 n^(th) 项。斐波那契数列中的任何项都是前两项之和，可以正式定义为如下：
fib(0)=0
fib(1)=1
fib(n)=fib(n-1)+fib(n-2)
在这里，fib(n)表示斐波那契数列中的n^(th)项。根据定义，我们可以轻松地计算出斐波那契数列：0, 1, 1, 2, 3, 5, 8, 13。
现在，假设我们想写一个函数，它将返回斐波那契数列中的n^(th)项。写这个函数的一个简单方法是使用递归，如下所示：
def fibonacci(n):
"""
Returns the n-th number in the Fibonacci sequence.
Parameters
----------
n: int
The n-th number in the Fibonacci sequence.
"""
if n <= 1:
return n
else:
return fibonacci(n-1) + fibonacci(n-2)
在上述代码中，我们有一个简单的if ... else条件语句，其中如果n小于或等于 1，我们返回n的值；否则，我们使用递归来计算数列中前两个数字的和。现在让我们尝试确定对于一个小的n值，比如 5，fibonacci函数的调用次数。我们的函数调用可能类似于以下方式：
fibonacci(5) = fibonacci(4) + fibonacci(3)
fibonacci(5) = (fibonacci(3) + fibonacci(2)) + (fibonacci(2) + fibonacci(1))
fibonacci(5) = ((fibonacci(2) + fibonacci(1)) + (fibonacci(1) + fibonacci(0))) + ((fibonacci(1) + fibonacci(0)) + fibonacci(1))
fibonacci(5) = (((fibonacci(1) + fibonacci(0)) + fibonacci(1)) + (fibonacci(1) + fibonacci(0))) + ((fibonacci(1) + fibonacci(0)) + fibonacci(1))
对于如此小的n值，我们仍然可以看到在相同参数的函数调用中存在重复。我们可以看到fibonacci(1)被调用了五次，而fibonacci(0)被调用了三次。如果我们再往上看一层，我们会看到fibonacci(2)也被多次调用。在这种情况下，计算仍然是可处理的，但对于较大的n值，这个函数的运行时间将呈指数增长；其时间复杂度为O(2^n)。为了让大家了解运行时间的增长速度，使用该算法计算斐波那契数列的第 1,000 项，我们所需要的时间将比我们可观察到的宇宙的年龄还要长，前提是计算机上有宇宙中所有的电子。
由于我们已经证明，对于任何适中的n值，计算斐波那契数列的n^(th)项是不可能的，因此我们将研究另一种基于动态规划的算法，代码如下所示：
cache = {0: 0, 1: 1} # Initialize the first two values.
def fibonacci(n):
"""
Returns the n-th number in the Fibonacci sequence.
Parameters
----------
n: int
The n-th number in the Fibonacci sequence.
"""
try:
return cache[n]
except KeyError:
fib = fibonacci(n-1) + fibonacci(n-2)
cache[n] = fib
return fib
在这种情况下，我们将每次调用函数的结果存储在字典中，这样我们就可以在O(1)的时间内访问它。由于这个缓存，我们只需要计算斐波那契数列的每一项一次。对于每次调用，我们首先检查是否已经计算过该值。如果已经计算过，我们直接从字典中访问该值，否则我们进行计算。这个算法的运行时复杂度为O(n)，因为我们每次只计算一次数列中的每一项。因此，我们可以看到，使用动态规划可以在运行时复杂度和内存复杂度之间进行权衡，从而将运行时复杂度从指数级降低到线性级。
如果我们考虑一下我们编写斐波那契数列的方式，我们首先尝试计算n^(th)项，然后计算缺失的值。这可以看作是一种自顶向下的方法。另一种方法是自底向上的方法，我们从计算第 0 项开始，然后逐步计算第一项、第二项，依此类推。动态规划的概念在这两种情况下是相同的，但在如何编写代码上有一个小的差异，如下所示：
cache = [0, 1] # Initialize with the first two terms of Fibonacci series.
def fibonacci(n):
"""
Returns the n-th number in the Fibonacci sequence.
Parameters
----------
n: int
The n-th number in the Fibonacci sequence.
"""
for i in range(2, n):
cache.append(cache[i-1] + cache[i-2])
return cache[-1]
如你所见，前面示例中的代码要简洁得多，我们不需要检查是否已经计算过这些值。这对于需要利用所有先前值来计算下一个值的问题非常有效。然而，如果我们不需要所有的先前值，我们最终会在自下而上的方法中计算一些不必要的值。
在接下来的章节中，我们将看到 HMM 中的推理算法使我们能够使用动态规划将问题拆解成子问题，从而使得计算变得可处理。
前向算法
现在我们正式定义前向算法中的问题。在前向算法的情况下，我们试图计算机器人在任何时间点的联合分布，利用到那个时间点为止传感器的输出，如下图所示：
前向算法：P(Z[k], X[1:k])
图 3.3：HMM 显示两个时间片，k-1和k
为了计算这个概率分布，我们将尝试将联合分布项拆分成更小的已知项。正如我们将看到的，我们可以为该分布写出一个递归公式。我们首先通过在分布P(Z[k], X[1:k])中引入一个新变量Z[k-1]，如下所示：
概率的边缘化法则是：
概率的乘法法则是：
在这里，我们基本上是使用概率的边缘化法则来引入Z[k-1]，然后对其状态进行求和。在这种情况下，我们假设Z[k-1]有m个状态，因此现在我们可以使用概率的乘法法则来将此项拆分如下：
在上一章中，我们看到了 HMM 的 d-分离属性如何使变量彼此独立。在这里，我们将应用其中的一些条件来简化前面方程中的项。正如我们所知道的，给定隐藏状态时，观察值与之前时间实例中的所有项相互独立，我们也知道：
。将其应用到我们前面方程的第一项，我们可以将其写为如下：
同样，我们知道当前的隐藏状态依赖于上一个隐藏状态，并且与之前的隐藏状态相互独立。因此，在这种情况下，公式是！。利用这一属性，我们可以将方程中的第二项写为如下：
现在，如果我们将最后一项与我们试图计算的项进行比较，我们应该能够看到它们之间的相似性。我们定义一个新的函数，α，如下所示：
现在我们可以将原始方程重写为如下：
所以，我们现在有了一个很好的递归方程，并且我们对方程中的所有项都很熟悉。第一个项，P(X[k]|Z[k])，是 HMM 的发射概率。方程的第二个项，P(Z[k]|Z[k-1])，是转移概率，也是已知的。现在我们可以集中精力求解这个递归方程。在求解任何递归方程时，我们需要知道至少一个值，以便开始计算连续的项。在这个例子中，我们知道α(1)的值，给出如下：
这里，P(Z[1])是我们已知的机器人位置的初始概率，而P(X[1]|Z[1])是已知的隐马尔可夫模型（HMM）的发射概率。通过这两个值，我们可以计算出α(1)的值。一旦得到α(1)的值，我们可以利用递归公式计算所有的α值。
现在让我们来讨论一下这个算法的计算复杂度，看看推理是否可行。正如我们在计算每个α的公式中看到的那样，我们对所有状态Z[k-1]进行求和，而我们假设它有m个状态；在每一步中，我们进行m次乘法操作来计算P(X[k]|Z[k])P(Z[k]|Z[k-1])α(k-1)。因此，为了计算下一个α，我们需要进行m²次计算。如果我们想要计算α(n)，我们将需要nm²次计算，这样我们就得到了算法的计算复杂度为O(nm²)。
现在让我们尝试编写代码来解决我们的机器人定位问题，利用前向算法计算联合分布。为了简化起见，假设我们只有一个传感器，它检查机器人左侧是否有墙壁。为了定义模型，我们主要需要两个量：转移矩阵和发射矩阵。在我们的例子中，我们假设机器人在任意时刻可以保持原位置，或者以相等的概率0.33向任意方向移动一个单位。因此，如果机器人在任意时刻处于位置 1，它可以出现在位置 1、2 或 6，且概率相等。
这样，我们可以将从状态 1 的转移概率写成：
以类似的方式，我们可以写出从每个位置的转移概率P(Z[t+1]|Z[t])，并得到以下转移矩阵：
import numpy as np
transition_matrix = np.array([[0.33, 0.33, 0, 0, 0, 0.33, 0, 0, 0, 0, 0, 0, 0],
[0.33, 0.33, 0.33, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[ 0, 0.25, 0.25, 0.25, 0, 0, 0.25, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0.33, 0.33, 0.33, 0, 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0.33, 0.33, 0, 0, 0.33, 0, 0, 0, 0, 0],
[0.33, 0, 0, 0, 0, 0.33, 0, 0, 0.33, 0, 0, 0, 0],
[ 0, 0, 0.33, 0, 0, 0, 0.33, 0, 0, 0, 0.33, 0, 0],
[ 0, 0, 0, 0, 0.33, 0, 0, 0.33, 0, 0, 0, 0, 0.33],
[ 0, 0, 0, 0, 0, 0.33, 0, 0, 0.33, 0.33, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0, 0, 0.33, 0.33, 0.33, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.33, 0.33, 0.33, 0],
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.33, 0.33, 0.33],
[ 0, 0, 0, 0, 0, 0, 0, 0.33, 0, 0, 0, 0.33, 0.33]])
在这个问题中，关于发射概率P(X[t]|Z[t])，我们应该为左侧有墙壁的状态分配发射概率 1。因此，发射概率将如下所示：
emission = np.array([1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0])
然而，由于我们不知道机器人在t=0时的位置，我们假设它在所有可能的状态上均匀分布。因此，初始概率P(Z[0])可以写成如下形式：
init_prob = np.array([0.077, 0.077, 0.077, 0.077, 0.077, 0.077, 0.077,
0.077, 0.077, 0.077, 0.077, 0.077, 0.077])
在掌握了这些值后，我们应该能够运行前向算法，但在此之前，我们需要按如下方式编码算法：
def forward(obs, transition, emission, init):
"""
Runs forward algorithm on the HMM.
Parameters
----------
obs: 1D list, array-like
The list of observed states.
transition: 2D array-like
The transition probability of the HMM.
size = {n_states x n_states}
emission: 1D array-like
The emission probability of the HMM.
size = {n_states}
init: 1D array-like
The initial probability of HMM.
size = {n_states}
Returns
-------
float: Probability value for the obs to occur.
"""
n_states = transition.shape[0]
fwd = [{}]
for i in range(n_states):
fwd[0][y] = init[i] * emission[obs[0]]
for t in range(1, len(obs)):
fwd.append({})
for i in range(n_states):
fwd[t][i] = sum((fwd[t-1][y0] * transition[y0][i] * emission[obs[t]]) for y0 in
range(n_states))
prob = sum((fwd[len(obs) - 1][s]) for s in range(n_states))
return prob
我们可以通过对一些观察结果进行计算来尝试计算该概率，如下所示：
>>> obs = [0, 0, 0, 0] # Staying in the same location
>>> forward(obs, transition_matrix, emission, init_prob)
0.97381776799999997
>>> obs = [0, 10, 8, 6] # Should be 0 because can't jump from state 0 to 10.
>>> forward(obs, transition_matrix, emission, init_prob)
0.0
计算给定观察结果的隐藏状态的条件分布
使用前向算法，我们已经能够计算出P(Z[x], X)的值，因此我们可能会认为，我们可以轻松地使用以下乘法规则来计算条件分布P(Z[k], X)：
然而，计算分布P(X)是计算上不可行的，正如我们将看到的那样。我们可以将P(Z[k], X)表示为：
因此，我们可以计算P(X)如下：
如果我们查看使用上述方程计算P(X)的计算复杂度，它是O(m^n)，对于任何足够大的m和n值而言，这都是不可行的。因此，仅使用前向算法计算隐藏状态的条件分布是不可能的。在接下来的章节中，我们将介绍反向算法，并向您展示如何计算这些条件分布。
反向算法
让我们正式定义反向算法的问题描述。在这种情况下，我们试图计算给定当前状态下观察变量的概率：
反向算法：P(X[k+1:n]|Z[k])
图 3.4：HMM 展示了两个时间片，k和k+1
类似于我们在前向算法中所做的那样，我们将再次尝试将这个概率项转换为一个递归方程，以已知分布为基础，从而递归地计算不同时间点的概率。首先，我们引入一个新项Z[k+1]，在P(X[k+1:n]|Z[k])中，使用边际化规则：
在这里，我们通过对所有可能的状态求和来对Z[k+1]变量进行边际化，我们假设它的取值为m。现在，我们可以使用概率的乘法规则将上述方程分解为：
现在，从我们模型的 d-分离性质中，我们知道！。此外，我们从 HMM 的定义中知道！。利用这些独立条件，我们可以将我们的方程写成：
现在，我们方程中的项看起来很熟悉。第二项是发射概率，最后一项是我们 HMM 的转移概率。现在，为了将其表示为递归，让我们定义一个新的函数β，如以下所示：
我们可以在前一个方程中使用β来表示它作为递归：
现在，由于我们已经有了递归方程，我们可以开始计算不同的 β[k] 值。但在计算这些值时，我们至少需要知道递归中的一个项，所以我们先计算 β(1) 的值：
def backward(obs, transition, emission, init):
"""
Runs backward algorithm on the HMM.
Parameters
----------
obs: 1D list, array-like
The list of observed states.
transition: 2D array-like
The transition probability of the HMM.
size = {n_states x n_states}
emission: 1D array-like
The emission probabilitiy of the HMM.
size = {n_states}
init: 1D array-like
The initial probability of HMM.
size = {n_states}
Returns
-------
float: Probability value for the obs to occur.
"""
n_states = transition.shape[0]
bkw = [{} for t in range(len(obs))]
T = len(obs)
for y in range(n_states):
bkw[T-1][y] = 1
for t in reversed(range(T-1)):
for y in range(n_states):
bkw[t][y] = sum((bkw[t+1][y1] * transition[y][y1] * emission[obs[t+1]]) for y1 in
range(n_states))
prob = sum((init[y] * emission[obs[0]] * bkw[0][y]) for y in range(n_states))
return prob
我们还可以在相同的观察值上运行这个算法，看看结果是否正确：
>>> obs = [0, 0, 0, 0] # Staying in the same location
>>> backward(obs, transition_matrix, emission, init_prob)
0.97381776799999997
>>> obs = [0, 10, 8, 6] # Should be 0 because can't jump from state 0 to 10.
>>> backward(obs, transition_matrix, emission, init_prob)
0.0
前向-后向算法（平滑）
进入前向-后向算法，我们现在正试图计算给定观察值下隐藏状态的条件分布。
以我们的机器人定位为例，我们现在试图找到给定传感器读数时，机器人在某个时刻位置的概率分布：
前向-后向算法：P(Z[k]|X)
图 3.5：HMM 显示三个时间片，k-1，k 和 k+1。
现在，由于模型中所有观察到的变量都已给定，我们可以说 P(Z[k]|X) 的值将与 Z[k] 和 X 的联合分布成正比：
现在，我们知道可以将 X={X[1:k]， X[k+1:n]}。将其代入前面的方程中，我们得到：
我们可以应用前面的链式法则将其写为：
从我们的模型结构来看，我们知道 
，并利用这一独立性性质，我们可以将前面的方程写为：
现在，如果我们看一下前面的项，第一项是 P(X[k+1:n]|Z[k])，这是我们在后向算法中计算的结果。第二项是 P(Z[k], X[1:k])，这是我们在前向算法中计算的结果。因此，在计算 P(Z[k|]X) 时，我们可以使用前向和后向算法计算这两个项。但由于 P(Z[k|]X) 与这两个项的乘积成正比，我们需要对分布进行归一化。
维特比算法
到目前为止，我们一直在尝试计算模型中的不同条件和联合概率。但前向-后向算法无法做的一件事是，给定观察值后，找出模型中隐藏变量的最可能状态。形式上，我们可以将这个问题写为：我们知道观察变量、转移概率和网络的发射概率，并且我们希望计算 *Z^**，其定义为：
其中，
Z={Z[1]， Z[2]， …， Z[n]}
另外，
X={X[1]， X[2]， …， X[n]}
概率分布操作的性质：
当我们对概率分布进行操作（边缘化、最大化等）时，可以通过分布的独立项将操作推进去。我们可以在边缘化和 argmax 的例子中看到这些：
图 3.6：显示三个时间切片 k-2、k-1 和 k 的 HMM
由于我们看到 P(Z|X)∝ P(Z, X)，并且因为我们正在尝试计算 argmax，所以无论我们在哪两个项中计算，结果都不重要。因此我们可以说：
现在，我们将再次尝试将我们的方程形式化为递归，以便更容易进行计算。那么，接下来我们引入一个新的项，µ(k)，定义为：
再次，我们将尝试将此项拆分为已知项。使用链式法则，我们可以将其写为：
现在，我们开始使用该性质推入 argmax 参数（请参见信息框了解详细信息）。这将给出我们：
这些术语看起来很熟悉，P(X[k]|Z[k]) 是发射概率，P(Z[k]|Z[k-1]) 是转移概率，且 
是 µ(k-1)。现在我们有一个递归方程可以使用：
既然我们有了递归公式，如果我们知道第一个项的值，就可以计算任何 k 的值。那么，首先让我们来看一下递归的第一个项，即 µ(1)：
这里，第一个项是 P(Z[1])，即我们的初始概率，是已知的。第二项是 P(X[1]|Z[1])，即我们模型的发射概率：
import numpy as np
def viterbi(obs, transition, emission, init=None):
"""
Return the MAP estimate of state trajectory of Hidden Markov Model.
Parameters
----------
y : array (T,)
Observation state sequence. int dtype.
transition : array (K, K)
State transition matrix. See HiddenMarkovModel.state_transition for
details.
emission : array (K,)
Emission matrix. See HiddenMarkovModel.emission for details.
init: optional, (K,)
Initial state probabilities: Pi[i] is the probability x[0] == i. If
None, uniform initial distribution is assumed (Pi[:] == 1/K).
Returns
-------
x : array (T,)
Maximum a posteriori probability estimate of hidden state trajectory,
conditioned on observation sequence y under the model parameters.
T1: array (K, T)
the probability of the most likely path so far
T2: array (K, T)
the x_j-1 of the most likely path so far
"""
# Cardinality of the state space
K = transition.shape[0]
emission = np.repeat(emission[np.newaxis, :], K, axis=0)
# Initialize the priors with default (uniform dist) if not given by caller
init = init if init is not None else np.full(K, 1 / K)
T = len(obs)
T1 = np.empty((K, T), 'd')
T2 = np.empty((K, T), 'B')
# Initilaize the tracking tables from first observation
T1[:, 0] = init * emission[:, obs[0]]
T2[:, 0] = 0
# Iterate throught the observations updating the tracking tables
for i in range(1, T):
T1[:, i] = np.max(T1[:, i - 1] * transition.T * emission[np.newaxis, :, obs[i]].T, 1)
T2[:, i] = np.argmax(T1[:, i - 1] * transition.T, 1)
# Build the output, optimal model trajectory
x = np.empty(T, 'B')
x[-1] = np.argmax(T1[:, T - 1])
for i in reversed(range(1, T)):
x[i - 1] = T2[x[i], i]
return x, T1, T2
我们可以用相同的观测值来试试：
>>> x, T1, T2 = viterbi([0, 0, 0, 0], transition_matrix, emission, init_prob)
>>> print(x)
array([0, 0, 0, 0], dtype=uint8)
>>> print(T1)
array([[ 0.077, 0.02541, 0.0083853, 0.00276715],
[ 0.077, 0.02541, 0.0083853, 0.00276715],
[ 0.077, 0.02541, 0.0083853, 0.00276715],
[ 0.077, 0.02541, 0.0083853, 0.00276715],
[ 0.077, 0.02541, 0.0083853, 0.00276715],
[ 0.077, 0.02541, 0.0083853, 0.00276715],
[ 0.077, 0.02541, 0.0083853, 0.00276715],
[ 0.077, 0.02541, 0.0083853, 0.00276715],
[ 0.077, 0.02541, 0.0083853, 0.00276715],
[ 0.077, 0.02541, 0.0083853, 0.00276715],
[ 0.077, 0.02541, 0.0083853, 0.00276715],
[ 0.077, 0.02541, 0.0083853, 0.00276715],
[ 0.077, 0.02541, 0.0083853, 0.00276715]])
>>> print(T2)
array([[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 1, 1, 1],
[ 0, 3, 3, 3],
[ 0, 3, 3, 3],
[ 0, 0, 0, 0],
[ 0, 6, 6, 6],
[ 0, 4, 4, 4],
[ 0, 5, 5, 5],
[ 0, 8, 8, 8],
[ 0, 6, 6, 6],
[ 0, 10, 10, 10],
[ 0, 7, 7, 7]], dtype=uint8)
>>> x, T1, T2 = viterbi([0, 10, 8, 6], transition_matrix, emission, init_prob)
>>> print(x)
array([0, 0, 0, 0], dtype=uint8)
>>> print(T1)
array([[ 0.077, 0., 0., 0\. ],
[ 0.077, 0., 0., 0\. ],
[ 0.077, 0., 0., 0\. ],
[ 0.077, 0., 0., 0\. ],
[ 0.077, 0., 0., 0\. ],
[ 0.077, 0., 0., 0\. ],
[ 0.077, 0., 0., 0\. ],
[ 0.077, 0., 0., 0\. ],
[ 0.077, 0., 0., 0\. ],
[ 0.077, 0., 0., 0\. ],
[ 0.077, 0., 0., 0\. ],
[ 0.077, 0., 0., 0\. ],
[ 0.077, 0., 0., 0\. ]])
>>> print(T2)
array([[ 0, 0, 0, 0],
[ 0, 0, 0, 0],
[ 0, 1, 0, 0],
[ 0, 3, 0, 0],
[ 0, 3, 0, 0],
[ 0, 0, 0, 0],
[ 0, 6, 0, 0],
[ 0, 4, 0, 0],
[ 0, 5, 0, 0],
[ 0, 8, 0, 0],
[ 0, 6, 0, 0],
[ 0, 10, 0, 0],
[ 0, 7, 0, 0]], dtype=uint8)
总结
在本章中，我们介绍了用于推理我们隐马尔可夫模型（HMM）的算法。我们研究了前向-后向算法，用于根据观测结果预测我们的隐藏状态。我们还讨论了维特比算法，用于计算模型中最可能的状态。
在所有这些算法中，我们假设已经知道模型的转移概率和发射概率。但是在现实世界中的问题中，我们需要从数据中计算这些值。在下一章中，我们将介绍用于计算转移和发射概率的最大似然法算法。
第四章：使用最大似然估计进行参数学习
在上一章中，我们讨论了隐马尔可夫模型（HMM）中的状态推断。我们尝试使用前一个状态转移的信息预测 HMM 的下一个状态。但在每种情况下，我们假设我们已经知道模型的转移和发射概率。然而，在实际问题中，我们通常需要从观察数据中学习这些参数。
在本章中，我们将尝试通过从观察中收集到的数据来估计我们的 HMM 模型的参数。我们将覆盖以下主题：
-
最大似然学习，带有示例
-
HMM 中的最大似然学习
-
期望最大化算法
-
Baum-Welch 算法
最大似然学习
在深入了解 HMM 中的最大似然估计（MLE）之前，我们先来理解一下 MLE 在一般情况下的基本概念。顾名思义，MLE 试图选择最大化观察数据似然性的模型参数。对于给定参数的任何模型，似然性定义为获取观察数据的概率，可以写成如下形式：
这里，D={D[1], D[2], D[3], …, D[n]}是观察数据，θ是控制我们模型的参数集。在大多数情况下，为了简便起见，我们假设数据点是独立同分布（IID）的。在这个假设下，我们可以将似然函数的定义简化为如下形式：
在这里，我们使用了独立随机变量的乘法规则，将联合分布分解为单个数据点的乘积。
回到最大似然估计，MLE 试图找到使得P(D|θ)值最大化的θ值。所以，基本上我们现在面临的是一个优化问题：
在接下来的几个小节中，我们将尝试将最大似然估计（MLE）应用于一些简单的示例，以便更好地理解它。
硬币投掷中的 MLE
假设我们想通过投掷硬币获得的观察数据来学习一个硬币模型。由于硬币只有两种结果：正面或反面，因此可以用一个参数来建模。假设我们将该参数定义为θ，即投掷硬币时获得正面的概率。获得反面的概率将自动为1-θ，因为正面和反面是互斥事件。
我们已经准备好了模型，接下来我们来计算该模型的似然函数。假设我们得到了一些硬币投掷的观察数据，D={H,H,T,H,T,T}。对于给定的数据，我们可以将似然函数写成如下形式：
现在，我们希望找到能够最大化P(D|θ)的θ值。为此，我们需要对我们的似然函数进行求导，将其等于0，然后求解θ：
因此，我们的 MLE 估计器学到抛硬币得到正面的概率是0.5。根据我们的观察，由于正反面出现的次数相等，我们也预期这个概率应该是相同的。
现在让我们尝试编写代码，来学习模型的参数θ。但是，如我们所知，在计算机上寻找最优值时可能会遇到数值问题，是否有可能避免这些问题，直接计算θ[MLE]呢？如果我们仔细观察我们的似然方程，我们会发现我们可以为这个模型写出一个通用的似然公式。如果我们假设我们的数据中有n个正面和m个反面，我们可以将似然写成如下形式：
现在，我们实际上可以使用这个似然函数找到θ[MLE]的封闭形式解，避免依赖任何数值方法来计算最优值：
我们可以看到，我们已经能够找到一个关于 MLE 解θ的封闭形式解。现在，编写代码就是简单地计算前面的公式，如下所示：
import numpy as np
def coin_mle(data):
"""
Returns the learned probability of getting a heads using MLE.
Parameters
----------
data: list, array-like
The list of observations. 1 for heads and 0 for tails.
Returns
-------
theta: The learned probability of getting a heads.
"""
data = np.array(data)
n_heads = np.sum(data)
return n_heads / data.size
现在，让我们尝试对不同的数据点应用我们的函数：
>>> coin_mle([1, 1, 1, 0, 0])
0.59999999999999998
>>> coin_mle([1, 1, 1, 0, 0, 0])
0.5
>>> coin_mle([1, 1, 1, 0, 0, 0, 0])
0.42857142857142855
输出结果如我们预期，但 MLE 方法的一个缺点是它对数据中的随机性非常敏感，在某些情况下，这可能导致其学习到错误的参数。尤其是在数据集较小时，这种情况尤为明显。例如，假设我们抛一枚公平的硬币三次，每次都得到正面。在这种情况下，MLE 方法会学习到θ的值为 1，这显然是不正确的，因为我们抛的是一枚公平的硬币。输出结果如下：
>>> coin_mle([1, 1, 1])
1.0
在第五章《使用贝叶斯方法进行参数推断》中，我们将尝试通过从参数的先验分布开始来解决这个 MLE 问题，并且随着数据量的增加，它会调整其先验分布。
正态分布的最大似然估计（MLE）
在上一节中，我们讨论了一个包含单个参数的模型。在本节中，我们将把相同的概念应用到一个稍微复杂的模型中。我们将尝试从给定的观测数据中学习正态分布的参数（也称为高斯分布）。如我们所知，正态分布由其均值和标准差来参数化，分布的公式如下：
这里，µ是正态分布的均值，σ是标准差。
正如我们之前讨论的，使用 MLE 来估计参数时，我们需要一些观测数据，在本例中，我们假设这些数据来自正态分布（或者可以使用正态分布来逼近）。假设我们有一些观测数据：X = {x[1], x[2], ..., x[N]}。我们想要估计我们的模型的参数μ（均值）和σ²（方差）。
我们将遵循与上一节相同的步骤。我们将首先定义正态分布的似然函数。似然函数是给定参数的情况下，数据被观测到的概率。所以，给定观测数据，我们可以将似然函数写成如下形式：
在尝试处理小数值的乘积时，我们通常会遇到一个问题，那就是数字可能变得太小，计算机无法处理。为了避免遇到这个问题，我们转而使用对数似然，而不是简单的似然。由于对数是一个递增函数，最大对数似然函数的值与似然函数的最大值对应的参数是一样的。对数似然可以定义如下：
我们可以通过对每个变量分别求偏导数，令其等于0，并解方程来找到最大化对数似然函数的μ[MLE]和σ[MLE]的值。为了得到均值，我们需要对对数似然函数关于μ求偏导数，同时保持σ为常数，并设其为0，得到以下结果：
类似地，标准差σ²的 MLE 可以通过对对数似然函数关于σ²求偏导数，同时保持μ为常数，令其等于0，然后解出σ²：
正如我们所见，我们再次能够得出一个闭式解，因此在编程时无需依赖数值方法。让我们尝试编写代码，检查我们的 MLE 方法是否学习到了正确的参数：
import numpy as np
def gaussian_mle(data):
"""
Returns the learned parameters of the Normal Distribution using MLE.
Parameters
----------
data: list, array-like
The list of observed variables.
Returns
-------
\mu: The learned mean of the Normal Distribution.
\sigma: The learned standard deviation of the Normal Distribution.
"""
data = np.array(data)
mu = np.mean(data)
variance = np.sqrt(np.mean((data - mu)**2))
return mu, variance
我们已经准备好学习函数，现在可以从已知分布生成一些数据，并检查我们的函数是否能够从生成的数据中学习到相同的参数：
>>> from numpy.random import normal
>>> data = normal(loc=1, scale=2, size=10)
>>> data
array([ 1.8120102, 2.14363679, 1.49010868, -1.95531206, 1.62449155,
1.49345327, 1.48957918, -0.67536313, 4.31506202, 4.24883442])
>>> mu, sigma = gaussian_mle(data)
>>> mu
1.5986500906187573
>>> sigma
1.805051208889392
在这个例子中，我们可以看到学到的值并不是非常准确。这是因为 MLE 对于观测数据点过于敏感，正如我们在上一节中讨论的那样。让我们尝试使用更多的观测数据运行这个相同的例子：
>>> data = normal(loc=1, scale=2, size=1000)
>>> data[:10]
array([ 4.5855015, 1.55162883, -1.61385859, 0.52543984, 0.90247428,
3.40717092, 1.4078157, 0.01560836, -1.19409859, -0.01641439])
>>> mu, sigma = gaussian_mle(data)
>>> mu
1.0437186891666821
>>> sigma
1.967211026428509
在这种情况下，使用更多的数据，我们可以看到学到的值与原始值更接近。
隐马尔可夫模型的最大似然估计（MLE）
在基本了解 MLE 的基础上，我们现在可以将这些概念应用于 HMM 的情况。在接下来的几个小节中，我们将看到 HMM 学习的两种可能场景，即监督学习和无监督学习。
监督学习
在监督学习的情况下，我们使用通过对我们试图建模的过程进行采样生成的数据。如果我们试图使用简单的离散分布来参数化我们的 HMM 模型，我们可以直接应用最大似然估计（MLE）通过统计从任何给定状态到其他状态的转移次数来计算转移和发射分布。类似地，我们可以通过统计不同隐状态下的输出状态来计算发射分布。因此，转移和发射概率可以通过以下方式计算：
在这里，T(i,j) 是从状态 i 到状态 j 的转移概率。而 E(i,s) 是从状态 i 获取状态 s 的发射概率。
让我们举一个非常简单的例子来使这个概念更加清晰。我们希望建模天气以及是否会在一段时间内下雨。同时，我们假设天气可以有三种可能的状态：
-
晴天 (S)
-
多云 (C)
-
有风 (W)
并且 雨天 变量有两个可能的状态：下雨 (R) 或 不下雨 (NR)。一个 HMM 模型大致是这样的：
假设我们有一些观察数据，如 D={(S,NR), (S,NR), (C,NR), (C,R), (C,R), (W,NR), (S,NR), (W,R), (C,NR)}。在这里，每个数据点的第一个元素表示当天观察到的天气，而第二个元素表示当天是否下雨。现在，使用我们之前推导出的公式，我们可以轻松地计算转移和发射概率。我们将从计算 S 到 S 的转移概率开始：
同样，我们可以计算所有其他状态组合的转移概率：
因此，我们就得到了天气的所有可能状态下的完整转移概率。我们可以以表格形式表示它，使其看起来更清晰：
| 晴天(S) | 多云(C) | 有风(W) | |
|---|---|---|---|
| 晴天(S) | 0.33 | 0.33 | 0.33 |
| 多云(C) | 0 | 0.66 | 0.33 |
| 有风 (W) | 0.5 | 0.5 | 0 |
表 1：天气模型的转移概率
现在，计算发射概率时，我们可以再次按照之前推导的公式进行：
同样，我们可以计算分布中所有其他值：
因此，我们的发射概率可以写成表格形式如下：
| 晴天(S) | 多云(C) | 有风(W) | |
|---|---|---|---|
| 雨天 (R) | 0 | 0.5 | 0.5 |
| 无雨 (NR) | 1 | 0.5 | 0.5 |
表 2：天气模型的发射概率
在前面的示例中，我们看到如何使用最大似然估计（MLE）和一些简单计算来计算 HMM 的参数。但由于在此案例中，我们假设了转移和发射概率是简单的离散条件分布，因此计算更为简便。在更复杂的情况下，我们需要估计比在前一节中正态分布情况更多的参数。
代码
现在让我们尝试编写前述算法的代码：
def weather_fit(data):
"""
Learn the transition and emission probabilities from the given data
for the weather model.
Parameters
----------
data: 2-D list (array-like)
Each data point should be a tuple of size 2 with the first element
representing the state of *Weather* and the second element representing
whether it rained or not.
Sunny = 0, Cloudy = 1, Windy = 2
Rain = 0, No Rain = 1
Returns
-------
transition probability: 2-D array
The conditional distribution representing the transition probability
of the model.
emission probability: 2-D array
The conditional distribution representing the emission probability
of the model.
"""
data = np.array(data)
transition_counts = np.zeros((3, 3))
emission_counts = np.zeros((3, 2))
for index, datapoint in enumerate(data):
if index != len(data)-1:
transition_counts[data[index][0], data[index+1][0]] += 1
emission_counts[data[index][0], data[index][1]] += 1
transition_prob = transition_counts / np.sum(transition_counts, axis=0)
emission_prob = (emission_counts.T / np.sum(emission_counts.T, axis=0)).T
return transition_prob, emission_prob
让我们生成一些数据，并尝试使用之前的函数学习参数：
>>> import numpy as np
>>> weather_data = np.random.randint(low=0, high=3, size=1000)
>>> rain_data = np.random.randint(low=0, high=2, size=1000)
>>> data = list(zip(weather_data, rain_data))
>>> transition_prob, emission_prob = weather_fit(data)
>>> transition_prob
array([[ 0.3125, 0.38235294, 0.27272727],
[ 0.28125, 0.38235294, 0.36363636],
[ 0.40625, 0.23529412, 0.36363636]])
>>> emission_prob
array([[ 0.3125, 0.38235294, 0.27272727],
[ 0.28125, 0.38235294, 0.36363636],
[ 0.40625, 0.23529412, 0.36363636]])
无监督学习
在前一节中，我们看到当所有变量（包括隐藏变量）都被观察到时，如何使用监督学习。但在实际问题中，这种情况通常并不成立。对于这种情况，我们使用无监督学习来估计模型的参数。
用于此目的的两种主要学习算法如下：
-
维特比学习算法
-
Baum-Welch 算法
我们将在接下来的几个子节中讨论这些内容。
维特比学习算法
维特比学习算法（与用于状态估计的维特比算法不同）接收一组训练观测值 O^r，其中 1≤r≤R，并通过迭代计算维特比对齐来估计单个 HMM 的参数。当用于初始化一个新的 HMM 时，维特比分割被均匀分割替代（即，每个训练观测值被划分为 N 等分）用于第一次迭代。
除了在新模型的第一次迭代中，每个训练序列 O 都使用状态对齐过程进行分割，该过程是通过最大化以下公式得到的：
对于 1<i<N，其中：
初始条件由以下给出：
对于 1<j<N。在离散情况下，输出概率 
定义为：
其中 S 是流的总数，vs 是给定的输出，输入 O[st]，以及 P[js][v] 是状态 j 给定输出 v 的概率。
如果 A[ij] 表示在执行上述最大化时，从状态 i 到状态 j 的总转移次数，则可以从相对频率中估计转移概率：
最大化 ∅N 的状态序列意味着训练数据观测与状态的对齐。在每个状态内，还会进一步将观测值对齐到混合成分。通常，对于每个状态和每个输出流，可以使用两种机制：
-
使用聚类将每个观测值 O[st] 分配给一个 M[s] 集群
-
将每个观测值 O[st] 与具有最高概率的混合成分关联
在任何情况下，最终结果是每个观察值都与一个唯一的混合成分相关联。这个关联可以通过指示函数 
表示，当 
与状态为j的流s的混合成分m相关联时，值为1，否则为零。
均值和方差随后通过计算简单的均值来估计：
混合权重基于分配给每个成分的观察数量：
Baum-Welch 算法（期望最大化）
期望最大化（EM）算法（在应用于 HMM 时称为Baum-Welch）是一种迭代方法，用于求解统计模型中依赖于未观察到的潜在变量的参数的最大似然估计或最大后验估计（MAP）。EM 迭代在执行期望（E）步骤和最大化（M）步骤之间交替进行，E步骤创建一个期望的对数似然函数，该函数使用当前的参数估计进行评估，M步骤则计算最大化E步骤找到的期望对数似然的参数。这些参数估计随后用于确定在下一个E步骤中潜在变量的分布。
EM 算法从参数的初始值（θ(old)*）开始。在**E**步骤中，我们取这些参数并找到潜在变量的后验分布*P(Z|X,θ(old))。然后，我们使用这个后验分布来评估完整数据似然函数对数的期望，作为参数θ的函数，得到函数Q(θ,θ^(old))，其定义如下：
让我们引入一些术语，这些术语可以帮助我们将来使用。γ(Z[n])表示潜在变量的边际后验分布：
ξ(z[n-1], z[n])表示两个连续潜在变量的边际后验分布：
因此，对于每个n值，我们可以将 γ(Z[n]) 存储为一个K维的非负数向量，且这些数的和为1，同样我们可以使用一个K×K的非负数矩阵来存储ξ(z[n-1], z[n])，其和也为1。
正如我们在前几章中讨论的，潜在变量z[n]可以表示为K维的二元变量，当z[n]处于状态k时，z[nk] = 1。我们还可以用它表示z[nk] = 1的条件概率，同样，ξ(z[n-1], j, z[nk])可以表示z[n-1]、j = 1和z[nk] = 1的条件概率。由于二元随机变量的期望就是其值为1的概率，我们可以得出以下结论：
正如我们在上一章讨论的那样，HMM 的联合概率分布可以表示如下：
因此，我们可以将数据似然函数写成如下形式：
在E步骤中，我们尝试有效地评估量γ(z[n])和ξ(z[n-1], z[n])。为了高效地计算这两个项，我们可以使用前向后向算法或维特比算法，正如上一章所讨论的那样。在M步骤中，我们尝试对Q(θ, θ^(old)) 关于参数θ={A, π, Φ}进行最大化，其中我们将γ(z[n])和ξ(z[n-1], z[n])视为常数。
通过这样做，我们得到了参数的 MLE 值，如下所示：
如果我们假设发射分布是正态分布，令 
，那么对Q(θ, θ^(old)) 关于Φ[k]的最大化将得到如下结果：
EM 算法必须通过选择π和A的初始值来初始化，这些值当然应该是非负的，并且总和为1。
代码
参数估计的算法看起来相当复杂，但hmmlearn，一个用于处理 HMM 的 Python 包，提供了很好的实现。hmmlearn也托管在 PyPI 上，因此可以直接通过pip:pip install hmmlearn安装。对于代码示例，我们将以通过学习股票价格的高斯 HMM 进行股票价格预测为例。这个示例来自hmmlearn的示例页面。
在这个示例中，我们还需要安装matplotlib和datetime包，这些也可以通过pip安装：
pip install matplotlib datetime
进入代码部分，我们应从导入所有必要的包开始：
from __future__ import print_function
import datetime
import numpy as np
from matplotlib import cm, pyplot as plt
from matplotlib.dates import YearLocator, MonthLocator
try:
from matplotlib.finance import quotes_historical_yahoo_ochl
except ImportError:
# For Matplotlib prior to 1.5.
from matplotlib.finance import (
quotes_historical_yahoo as quotes_historical_yahoo_ochl
)
from hmmlearn.hmm import GaussianHMM
print(__doc__)
接下来，我们将从 Yahoo! Finance 获取我们的股票价格数据：
quotes = quotes_historical_yahoo_ochl(
"INTC", datetime.date(1995, 1, 1), datetime.date(2012, 1, 6))
# Unpack quotes
dates = np.array([q[0] for q in quotes], dtype=int)
close_v = np.array([q[2] for q in quotes])
volume = np.array([q[5] for q in quotes])[1:]
# Take diff of close value. Note that this makes
# ``len(diff) = len(close_t) - 1``, therefore, other quantities also
# need to be shifted by 1.
diff = np.diff(close_v)
dates = dates[1:]
close_v = close_v[1:]
# Pack diff and volume for training.
X = np.column_stack([diff, volume])
接下来，我们定义一个高斯 HMM 模型，并为我们的数据学习参数：
# Make an HMM instance and execute fit
model = GaussianHMM(n_components=4, covariance_type="diag", n_iter=1000).fit(X)
# Predict the optimal sequence of internal hidden state
hidden_states = model.predict(X)
现在我们可以打印出我们学到的参数：
print("Transition matrix")
print(model.transmat_)
print()
print("Means and vars of each hidden state")
for i in range(model.n_components):
print("{0}th hidden state".format(i))
print("mean = ", model.means_[i])
print("var = ", np.diag(model.covars_[i]))
print()
输出如下：
Transition matrix
[[ 9.79220773e-01 2.57382344e-15 2.72061945e-03 1.80586073e-02]
[ 1.12216188e-12 7.73561269e-01 1.85019044e-01 4.14196869e-02]
[ 3.25313504e-03 1.12692615e-01 8.83368021e-01 6.86228435e-04]
[ 1.18741799e-01 4.20310643e-01 1.18670597e-18 4.60947557e-01]]
Means and vars of each hidden state
0th hidden state
mean = [ 2.33331888e-02 4.97389989e+07]
var = [ 6.97748259e-01 2.49466578e+14]
1st hidden state
mean = [ 2.12401671e-02 8.81882861e+07]
var = [ 1.18665023e-01 5.64418451e+14]
2nd hidden state
mean = [ 7.69658065e-03 5.43135922e+07]
var = [ 5.02315562e-02 1.54569357e+14]
3rd hidden state
mean = [ -3.53210673e-01 1.53080943e+08]
var = [ 2.55544137e+00 5.88210257e+15]
我们还可以绘制出隐藏状态随时间变化的图：
fig, axs = plt.subplots(model.n_components, sharex=True, sharey=True)
colours = cm.rainbow(np.linspace(0, 1, model.n_components))
for i, (ax, colour) in enumerate(zip(axs, colours)):
# Use fancy indexing to plot data in each state.
mask = hidden_states == i
ax.plot_date(dates[mask], close_v[mask], ".-", c=colour)
ax.set_title("{0}th hidden state".format(i))
# Format the ticks.
ax.xaxis.set_major_locator(YearLocator())
ax.xaxis.set_minor_locator(MonthLocator())
ax.grid(True)
plt.show()
前述代码的输出如下：
图 1：隐藏状态随时间变化的图
总结
在本章中，我们介绍了用于估计给定 HMM 模型参数的算法。我们从 MLE 的基本概念入手，然后将这些概念应用于 HMM。对于 HMM 的训练，我们讨论了两种不同的场景：有监督训练，当我们有隐藏状态的观测值时；以及无监督训练，当我们只有输出观测值时。
我们还讨论了使用 MLE 进行估计时遇到的问题。在下一章中，我们将介绍使用贝叶斯方法进行参数估计的算法，该方法试图解决这些问题。
第五章：使用贝叶斯方法进行参数推断
在上一章中，我们讨论了使用最大似然法推断参数。在本章中，我们将通过贝叶斯方法探讨相同的问题。主要内容如下：
-
贝叶斯学习简介
-
隐马尔可夫模型中的贝叶斯学习
-
用于估计分布的近似算法
贝叶斯学习
在最大似然学习方法中，我们试图找到使得似然函数最大化的最优参数。但现实中的数据通常是非常嘈杂的，并且在大多数情况下，它并不能代表真实的底层分布。在这种情况下，最大似然方法会失败。例如，考虑投掷一枚公平的硬币几次。可能我们的所有投掷结果都是正面或反面。如果我们对这些数据使用最大似然方法，它会将正面或反面分别赋予概率 1，意味着我们永远不会得到另一面。或者，假设我们投掷硬币 10 次，得到了三次正面和七次反面。在这种情况下，最大似然方法会将正面的概率定为 0.3，反面的概率定为 0.7，这并不是公平硬币的真实分布。这个问题也通常被称为过拟合。
贝叶斯学习采用略有不同的方法来学习这些参数。我们首先为模型的参数分配一个先验分布。先验使我们对模型的假设变得显式。在投掷硬币的例子中，我们可以从一个将正面和反面赋予相等概率的先验开始。然后，我们应用贝叶斯定理根据数据计算参数的后验分布。这使我们能够根据数据将我们的信念（先验）向数据的指向转移，从而进行更为温和的参数估计。通过这种方式，贝叶斯学习能够解决最大似然方法的一个主要缺点。
更一般来说，在贝叶斯学习的情况下，我们尝试学习模型参数的分布，而不是学习一个能够最大化似然的单一参数。为了学习这种参数分布，我们使用贝叶斯定理，公式如下：
这里，P(θ) 是我们对模型参数的先验分布，P(D|θ) 是给定参数下数据的似然，P(D) 是观察到数据的概率。P(D) 也可以通过先验和似然表示如下：
现在让我们分别讨论这些术语，看看我们如何计算它们。先验，P(θ)，是一个参数的概率分布，代表我们对参数值的信念。例如，在抛硬币的情况下，我们可以将初始信念设为θ在 0 到 1 之间，并且是均匀分布的。似然项，P(D|θ)，是我们在第四章中尝试最大化的相同项，使用最大似然进行参数推断。它表示给定模型参数的情况下，观测数据的可能性。下一个术语，P(D)，是观测到数据的概率，它充当归一化项。由于它要求我们对所有可能的θ值求和，因此计算上非常困难，对于任何足够多的参数，计算会迅速变得不可处理。在本章的接下来的部分中，我们将看到可以用来近似这些值的不同算法。我们试图计算的术语，P(D|θ)，被称为后验。它表示给定我们的观测数据后，模型参数的最终概率分布。基本上，我们的先验通过使用似然项来更新，得到最终的分布。
贝叶斯学习解决的另一个问题是模型选择。由于贝叶斯学习为不同的可能模型提供了一个分布，而不是单一的模型，我们有几种方式可以从这些模型中进行预测。第一种方法是选择具有最大概率的特定模型，这也通常被称为最大后验（MAP）估计。另一种可能的方法是基于后验分布计算所有模型的预测期望。这使我们能够对预测进行正则化，因为我们是在对所有可能的模型进行期望计算。
选择先验
在进行贝叶斯学习时，一个常见的问题是如何选择合适的先验。正如大卫·麦凯所说，没有假设就没有推理，我们需要对先验做出一个猜测。我们的先验应该代表我们认为模型中最可能的参数。使用我们自己的先验的一个巨大好处是，我们使得关于模型的假设变得明确。一旦我们开始应用贝叶斯定理，使用我们的先验和观测数据，后验分布就会从我们的先验向一个更能代表数据的分布偏移。
从理论上讲，这听起来不错，因为我们可以选择非常复杂的先验分布来捕捉我们对模型的理解，但在应用贝叶斯定理时，我们需要将先验与似然相乘，而对于复杂的分布，这很快会变得计算上不可处理。因此，在实践中，我们通常选择与似然分布共轭的先验分布。共轭先验使我们能够得到贝叶斯定理的闭式解。正因为如此，高斯分布被用作先验和似然，因为将一个高斯分布与另一个高斯分布相乘的结果仍然是一个高斯分布。而且，计算两个高斯分布的乘积并不昂贵。
不可处理性
除了选择困难的先验分布外，贝叶斯学习中的另一个不可处理性的来源是贝叶斯定理中的分母项：
正如我们在前面的公式中看到的，为了计算P(D)，我们需要对所有可能的θ值求和，而θ是我们模型中所有参数的集合。如果我们模型中的参数很多，那么计算这个项在计算上是不可处理的，因为项的大小会随着参数数量的增加而呈指数增长。很多工作已经致力于近似这个值，正如我们在本章的下一节中将看到的那样。
HMM 中的贝叶斯学习
正如我们在前一节中看到的，在贝叶斯学习的情况下，我们假设所有变量都是随机变量，给它们分配一个先验，然后根据这个先验计算后验。因此，在 HMM 的情况下，我们可以对我们的转移概率、发射概率或观测状态的数量分配先验。
因此，我们需要解决的第一个问题是选择先验。从理论上讲，先验可以是模型参数上的任何分布，但在实践中，我们通常尝试使用与似然分布共轭的先验，这样我们就可以得到方程的闭式解。例如，在 HMM 输出为离散的情况下，常见的先验选择是 Dirichlet 分布。这主要有两个原因，第一个原因是 Dirichlet 分布是多项式分布的共轭分布，允许我们轻松地将它们相乘。
共轭分布：如果通过将先验与似然相乘得到的后验分布属于与先验分布相同的分布族，则称一族先验分布是与一族似然分布共轭的。
例如，由于给定π参数向量时，初始状态的似然是多项式分布：
如果π的先验概率是 Dirichlet 分布：
其中 u = [u[1], u[2], ..., u[K]] 是超参数向量，Z 是归一化常数。我们现在可以根据似然函数和先验来计算后验，其公式如下：
我们可以看到后验分布也是一个狄利克雷分布。因此，我们可以说狄利克雷分布是多项式分布的共轭先验。以类似的方式，我们可以为我们的转移矩阵和发射矩阵设定狄利克雷先验。
选择狄利克雷先验的第二个原因是它具有一个理想的性质，即其超参数可以解释为假设的观测次数。在前面的例子中，如果 u[i] = 2 且 u[j] = 1 （对于 j ≠ i），那么 π 的 MAP 估计将与最大似然估计相同，假设训练数据中有一个额外的数据点，其初始状态为 i。这种共轭性质使得我们可以通过对 Baum-Welch 算法做一个小的变动来进行 MAP 估计。这也为隐马尔可夫模型中看似临时但非常常见的正则化方法提供了理论依据，该方法只是将一个小的正数加到参数向量的所有元素上。
在前面的几段中，我们专门讨论了输出为离散的情况。但相同的概念也可以扩展到连续输出的情况。在连续分布的情况下也存在共轭分布。最常用的分布之一是高斯分布，因为它在进行不同操作后仍然保持在同一分布族中。
近似所需的积分
如前所述，贝叶斯方法将所有未知量视为随机变量。我们为这些变量分配先验分布，然后在观察到数据后估计这些变量的后验分布。在隐马尔可夫模型（HMM）的情况下，未知量包括 HMM 的结构，即状态的数量、网络的参数和隐藏状态。与最大似然估计或 MAP 估计不同，在这些方法中我们为这些参数找到点估计，而现在我们有这些参数的分布。这使得我们可以在模型结构之间进行比较，但为了做到这一点，我们需要对模型的参数和隐藏状态进行积分。这通常被称为贝叶斯积分。
由于这些积分在计算上难以处理，我们采用近似方法来计算这些值。在接下来的几个小节中，我们将概述其中的一些方法。这些方法的详细分析超出了本书的范围。
采样方法
采样方法是估计难以处理的分布的最常见方式之一。其基本思想是以某种方式从分布空间中采样，使得我们能从高概率区域获得更多的样本。然后基于这些样本，我们估计分布。
拉普拉斯近似
拉普拉斯近似使用中央极限定理，该定理从表现良好的先验和数据中得出，后验参数将在大量训练样本的极限中收敛到围绕参数 MAP 估计的高斯分布。为了使用拉普拉斯近似估计证据，首先在常规优化过程中找到 MAP 参数，然后在 MAP 估计处计算对数似然的海森矩阵。证据通过在 MAP 估计的 θ 处评估 P(θ,D)/P(θ|D) 比率来近似，分母使用高斯近似。拉普拉斯近似有几个缺点：
-
从参数计算海森矩阵通常非常昂贵。
-
高斯近似对于参数为正且总和为 1 的模型并不适用，特别是当相对于数据集的大小有许多参数时。
出于这些原因，拉普拉斯近似通常不用于隐马尔可夫模型（HMM）。
Stolke 和 Omohundro 的方法
在著名的论文 HMM 通过贝叶斯模型合并进行归纳 中，Stolke 和 Omohundro 提出了一种新的技术，用于近似 HMM 的贝叶斯积分。考虑一种情况，所有的 HMM 状态都是观察到的，且先验是狄利克雷分布。在这种情况下，当使用贝叶斯学习法学习参数时，后验也将是狄利克雷分布，然后证据积分可以表示为狄利克雷积分的乘积，这些积分可以轻松计算。因此，从某种意义上讲，我们可以说，证据积分难以处理的原因是状态和参数是隐藏的。
Stolke 和 Omohundro 提出了一种方法，旨在使用类似维特比算法的方法找到最可能的隐藏状态序列，并将此序列作为观察状态。使用这些观察值，我们可以轻松进行证据积分。该方法建议在这两个步骤之间迭代，逐步搜索模型结构，根据近似证据的比较合并或拆分状态。在他们的论文中，Stolke 和 Omohundro 表明，通过将对隐藏变量的积分替代为对参数的积分，这种方法能够获得良好的结果。
变分方法
变分方法是另一种非常常见的近似分布的方法。一般思路是从选择一个更简单的分布族开始，然后尝试找到这个分布的超参数，使得该分布尽可能接近我们的原始分布。用于确定两个分布接近度的度量有多种；最常用的度量是 Kullback-Leibler 散度。这种方法本质上将推理问题转化为一个优化问题，在这个问题中我们尝试最小化我们的散度度量。
在 HMMs 的情况下，我们通常假设隐藏状态与模型的参数是独立的。这使得我们可以同时近似隐藏状态和参数的分布。更具体地说，通过应用 Jensen 不等式两次，可以对证据进行下界估计：
变分贝叶斯方法通过迭代地最大化 
作为两个自由分布，Q(S) 和 Q(θ) 的泛函。在前述方程中，我们可以看到，这种最大化等价于最小化 Q(S)Q(θ) 和隐藏状态及 P(S,θ|D,M) 参数的联合后验之间的 KL 散度。David MacKay 首次提出了一种变分贝叶斯方法来进行 HMMs 中的学习。他假设先验是一个狄利克雷分布，并且假设参数与隐藏状态是独立的，他展示了最优的 Q(θ) 是一个狄利克雷分布。此外，他还展示了可以通过将前向-后向算法应用于具有伪参数的 HMM 来获得最优的 Q(S)，这些伪参数由 
给出，且可以对狄利克雷分布进行评估。因此，整个变分贝叶斯方法可以作为 Baum-Welch 算法的简单修改来实现。本质上，我们可以说，变分贝叶斯方法是 MAP 方法和 Stolke 与 Omohundro 方法的特殊情况的结合。这是非常有前景的，尤其是考虑到它在其他模型中已成功应用于非平凡的模型结构学习；它在 HMMs 及其扩展中的潜力尚未得到充分探索。
代码
目前，Python 中没有支持使用贝叶斯学习进行学习的包，而且编写完整的代码以适应本书的内容是非常困难的。尽管使用贝叶斯学习有很多优点，但在很多情况下，它通常在计算上是不可行的。基于这些原因，我们跳过了隐藏马尔可夫模型（HMMs）中的贝叶斯学习代码。
摘要
在本章中，我们讨论了在隐马尔可夫模型（HMMs）中应用贝叶斯学习来学习参数。贝叶斯学习相比最大似然估计器有一些优势，但除非我们有闭式解，否则它在计算上是相当昂贵的。闭式解只有在使用共轭先验时才有可能。接下来的章节中，我们将讨论隐马尔可夫模型在各种问题中的详细应用。
第六章：时间序列预测
在前几章中，我们详细讨论了 隐马尔可夫模型（HMMs）及与推断相关的各种算法。从本章开始，我们将讨论 HMM 的应用。
HMM（隐马尔可夫模型）能够预测和分析基于时间的现象。因此，它们可以用于语音识别、自然语言处理和金融市场预测等领域。在本章中，我们将探讨 HMM 在金融市场分析中的应用，主要是股票价格预测。
使用 HMM 进行股票价格预测
由于许多大公司对股票市场的浓厚兴趣，股票市场预测一直是过去活跃的研究领域之一。从历史上看，各种机器学习算法已被应用，且取得了不同程度的成功。然而，由于股票市场的非平稳性、季节性和不可预测性，股票预测仍然受到严重限制。仅凭前期的股票数据预测未来的股票走势更具挑战性，因为它忽略了多个外部因素。
如前所述，HMM 能够基于顺序观察到的数据建模隐藏状态转换。股票预测问题也可以看作是遵循相同的模式。股票价格取决于多种因素，这些因素通常对投资者来说是不可见的（隐藏变量）。这些潜在因素之间的转换会根据公司政策和决策、财务状况及管理决策发生变化，并且这些变化会影响股票价格（观察数据）。因此，HMM 非常适合股票价格预测问题。
在本章中，我们将尝试预测 Alphabet Inc.（GOOGL）、Facebook（FB）和 Apple Inc.（AAPL）的股票价格。
收集股票价格数据
我们将使用 pystock 数据 (data.pystock.com) 来获取历史股票价格数据。每天，在美国股市开盘前的 9:30 EST/EDT，pystock 爬虫会收集股票价格和财务报告，并将数据推送到仓库，包括前一天的开盘价、收盘价、最高价和最低价等信息。这些数据是按天收集的，意味着我们不会有小时或分钟级别的数据。
我们将尝试下载给定年份的 pystock 数据。由于数据集较大，我们将创建一个 Python 脚本来下载给定年份的数据，并且可以同时运行该程序来并行下载三个不同年份的所有数据：
"""
Usage: get_data.py --year=<year>
"""
import requests
import os
from docopt import docopt
# docopt helps parsing the command line argument in
# a simple manner (http://docopt.org/)
args = docopt(doc=__doc__, argv=None,
help=True, version=None,
options_first=False)
year = args['--year']
# Create directory if not present
year_directory_name = 'data/{year}'.format(year=year)
if not os.path.exists(year_directory_name):
os.makedirs(year_directory_name)
# Fetching file list for the corresponding year
year_data_files = requests.get(
'http://data.pystock.com/{year}/index.txt'.format(year=year)
).text.strip().split('\n')
for data_file_name in year_data_files:
file_location = '{year_directory_name}/{data_file_name}'.format(
year_directory_name=year_directory_name,
data_file_name=data_file_name)
with open(file_location, 'wb+') as data_file:
print('>>> Downloading \t {file_location}'.format(file_location=file_location))
data_file_content = requests.get(
'http://data.pystock.com/{year}/{data_file_name}'.format(year=year, data_file_name=data_file_name)
).content
print('<<< Download Completed \t {file_location}'.format(file_location=file_location))
data_file.write(data_file_content)
同时运行以下脚本来获取三个不同年份的数据：
python get_data.py --year 2015
python get_data.py --year 2016
python get_data.py --year 2017
一旦数据下载完成，我们将通过结合所有年份的相关数据来尝试获取每只前述股票的全部数据：
"""
Usage: parse_data.py --company=<company>
"""
import os
import tarfile
import pandas as pd
from pandas import errors as pd_errors
from functools import reduce
from docopt import docopt
args = docopt(doc=__doc__, argv=None,
help=True, version=None,
options_first=False)
years = [2015, 2016, 2017]
company = args['--company']
# Getting the data files list
data_files_list = []
for year in years:
year_directory = 'data/{year}'.format(year=year)
for file in os.listdir(year_directory):
data_files_list.append('{year_directory}/{file}'.format(year_directory=year_directory, file=file))
def parse_data(file_name, company_symbol):
"""
Returns data for the corresponding company
:param file_name: name of the tar file
:param company_symbol: company symbol
:type file_name: str
:type company_symbol: str
:return: dataframe for the corresponding company data
:rtype: pd.DataFrame
"""
tar = tarfile.open(file_name)
try:
price_report = pd.read_csv(tar.extractfile('prices.csv'))
company_price_data = price_report[price_report['symbol'] == company_symbol]
return company_price_data
except (KeyError, pd_errors.EmptyDataError):
return pd.DataFrame()
# Getting the complete data for a given company
company_data = reduce(lambda df, file_name: df.append(parse_data(file_name, company)),
data_files_list,
pd.DataFrame())
company_data = company_data.sort_values(by=['date'])
# Create folder for company data if does not exists
if not os.path.exists('data/company_data'):
os.makedirs('data/company_data')
# Write data to a CSV file
company_data.to_csv('data/company_data/{company}.csv'.format(company=company),
columns=['date', 'open', 'high', 'low', 'close', 'volume', 'adj_close'],
index=False)
运行以下脚本，以创建一个包含 GOOGL、FB 和 AAPL 股票所有历史数据的 .csv 文件：
python parse_data.py --company GOOGL
python parse_data.py --company FB
python parse_data.py --company AAPL
股票价格预测的特征
一旦我们拥有每个股票价格的数据，我们就可以开始预测股票的价格了。如前所述，我们每一天的特征非常有限，具体是该天的开盘价、收盘价、股票最高价和最低价。因此，我们将利用这些数据来计算股票价格。通常，我们希望在给定当天开盘价以及前几天的某些 d 天数据的情况下，计算该天的收盘价。我们的预测器会有 d 天的延迟。
让我们创建一个名为StockPredictor的预测器，它将包含所有预测特定公司在特定日期股票价格的逻辑。
我们不直接使用股票的开盘价、收盘价、最低价和最高价，而是尝试提取它们的分数变化，这些变化将用于训练我们的 HMM。随着我们深入研究，选择这些特征的原因会变得更加清晰。我们可以定义三个参数如下：
因此，对于股票价格预测 HMM，我们可以将单个观测值表示为一个向量，这些参数为 X[t] = < frac[change], frac[high], frac[low] >：
import pandas as pd
class StockPredictor(object):
def __init__(self, company, n_latency_days=10):
self._init_logger()
self.company = company
self.n_latency_days = n_latency_days
self.data = pd.read_csv(
'data/company_data/{company}.csv'.format(company=self.company))
def _init_logger(self):
self._logger = logging.getLogger(__name__)
handler = logging.StreamHandler()
formatter = logging.Formatter(
'%(asctime)s %(name)-12s %(levelname)-8s %(message)s')
handler.setFormatter(formatter)
self._logger.addHandler(handler)
self._logger.setLevel(logging.DEBUG)
@staticmethod
def _extract_features(data):
open_price = np.array(data['open'])
close_price = np.array(data['close'])
high_price = np.array(data['high'])
low_price = np.array(data['low'])
# Compute the fraction change in close, high and low prices
# which would be used a feature
frac_change = (close_price - open_price) / open_price
frac_high = (high_price - open_price) / open_price
frac_low = (open_price - low_price) / open_price
return np.column_stack((frac_change, frac_high, frac_low))
# Predictor for GOOGL stocks
stock_predictor = StockPredictor(company='GOOGL')
使用 HMM 预测价格
一旦我们从数据中提取了特征，就可以开始预测股票的价格了。我们希望预测某一天股票的收盘价，已知该天的开盘价和前几天的股票价格。
第一步是训练一个 HMM，从我们之前计算的给定观测序列中计算参数。由于观测值是连续随机变量的向量，我们必须假设发射概率分布是连续的。为了简化，我们假设它是一个多项式高斯分布，具有参数（μ 和 Σ）。因此，我们需要确定以下参数：过渡矩阵 A、先验概率 π，以及表示多项式高斯分布的 μ 和 Σ。
目前，假设我们有四个隐藏状态。在接下来的章节中，我们将探讨如何找到最优的隐藏状态数量。我们将使用hmmlearn包提供的GaussianHMM类作为我们的 HMM，并尝试使用它提供的fit()方法进行参数估计：
from hmmlearn.hmm import GaussianHMM
class StockPredictor(object):
def __init__(self, company, n_latency_days=10, n_hidden_states=4):
self._init_logger()
self.company = company
self.n_latency_days = n_latency_days
self.hmm = GaussianHMM(n_components=n_hidden_states)
self.data = pd.read_csv(
'data/company_data/{company}.csv'.format(company=self.company))
def fit(self):
self._logger.info('>>> Extracting Features')
feature_vector = StockPredictor._extract_features(self.data)
self._logger.info('Features extraction Completed <<<')
self.hmm.fit(feature_vector)
在机器学习中，我们将整个数据集划分为两类。第一类是训练数据集，用于训练模型。第二类是测试数据集，用于提供对最终模型在训练数据集上拟合的无偏评估。将训练数据集与测试数据集分开，避免了我们将数据过拟合到模型上。因此，在我们的案例中，我们也会将数据集分为两类，train_data用于训练模型，test_data用于评估模型。为此，我们将使用sklearn.model_selection模块提供的train_test_split方法：
train_test_split可以将数组或矩阵随机拆分成训练集和测试集。由于我们用的是顺序数据来训练 HMM，因此我们不希望随机拆分数据。为了防止随机拆分测试数据和训练数据，可以将shuffle=False作为参数传递。
from sklearn.model_selection import train_test_split
class StockPredictor(object):
def __init__(self, company, test_size=0.33,
n_latency_days=10, n_hidden_states=4):
self._init_logger()
self.company = company
self.n_latency_days = n_latency_days
self.hmm = GaussianHMM(n_components=n_hidden_states)
self._split_train_test_data(test_size)
def _split_train_test_data(self, test_size):
data = pd.read_csv(
'data/company_data/{company}.csv'.format(company=self.company))
_train_data, test_data = train_test_split(
data, test_size=test_size, shuffle=False)
self._train_data = _train_data
self._test_data = test_data
def fit(self):
self._logger.info('>>> Extracting Features')
feature_vector = StockPredictor._extract_features(self._train_data)
self._logger.info('Features extraction Completed <<<')
self.hmm.fit(feature_vector)
一旦模型训练完成，我们需要预测股票的收盘价。正如我们之前提到的，我们想预测给定开盘价的一天的收盘价。这意味着，如果我们能够预测给定日期的frac[change]，就能按如下方式计算收盘价：
因此，我们的问题归结为计算给定X[1],...,X[t]和 HMM 参数下的一天的X[t+1] = < frac[change], frac[high], frac[low] >观测向量的值！，即找到一个X[t+1]的值，使得后验概率P(X[t+1]|X[1],...,X[t],θ)最大化：
如你所见，一旦我们从最大化方程中去除所有与X[t+1]无关的参数，我们就剩下了一个问题：找到X[t+1]的值，该值能够优化P(X[1],...,X[t+1]|θ)的概率。我们在第四章《最大似然法中的参数学习》中遇到过这个问题，当时我们在评估给定模型参数下的序列概率。这个问题可以通过前向后向算法有效计算。
如果我们假设frac[change]是一个连续变量，问题的优化将是计算上非常困难的。因此，我们可以将这些分数变化划分为一些离散值，范围在两个有限变量之间（如下表所示），并尝试找到一组分数变化，< frac[change], frac[high], frac[low] >，使得概率P(X[1],...,X[t+1]|θ)最大化：
| 观测值 | 最小值 | 最大值 | 点数 |
|---|---|---|---|
| frac[change] | -0.1 | 0.1 | 20 |
| frac[high] | 0 | 0.1 | 10 |
| frac[low] | 0 | 0.1 | 10 |
因此，使用前面列出的离散值集合，我们需要进行（20 x 10 x 10 =) 2,000 次操作：
def _compute_all_possible_outcomes(self, n_steps_frac_change,
n_steps_frac_high, n_steps_frac_low):
frac_change_range = np.linspace(-0.1, 0.1, n_steps_frac_change)
frac_high_range = np.linspace(0, 0.1, n_steps_frac_high)
frac_low_range = np.linspace(0, 0.1, n_steps_frac_low)
self._possible_outcomes = np.array(list(itertools.product(
frac_change_range, frac_high_range, frac_low_range)))
现在，我们可以实现该方法来预测收盘价，如下所示：
def _get_most_probable_outcome(self, day_index):
previous_data_start_index = max(0, day_index - self.n_latency_days)
previous_data_end_index = max(0, day_index - 1)
previous_data = self._test_data.iloc[previous_data_end_index: previous_data_start_index]
previous_data_features = StockPredictor._extract_features(
previous_data)
outcome_score = []
for possible_outcome in self._possible_outcomes:
total_data = np.row_stack(
(previous_data_features, possible_outcome))
outcome_score.append(self.hmm.score(total_data))
most_probable_outcome = self._possible_outcomes[np.argmax(
outcome_score)]
return most_probable_outcome
def predict_close_price(self, day_index):
open_price = self._test_data.iloc[day_index]['open']
predicted_frac_change, _, _ = self._get_most_probable_outcome(
day_index)
return open_price * (1 + predicted_frac_change)
让我们尝试预测未来几天的收盘价格，并绘制两条曲线：
"""
Usage: analyse_data.py --company=<company>
"""
import warnings
import logging
import itertools
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from hmmlearn.hmm import GaussianHMM
from sklearn.model_selection import train_test_split
from tqdm import tqdm
from docopt import docopt
args = docopt(doc=__doc__, argv=None, help=True,
version=None, options_first=False)
# Supress warning in hmmlearn
warnings.filterwarnings("ignore")
# Change plot style to ggplot (for better and more aesthetic visualisation)
plt.style.use('ggplot')
class StockPredictor(object):
def __init__(self, company, test_size=0.33,
n_hidden_states=4, n_latency_days=10,
n_steps_frac_change=50, n_steps_frac_high=10,
n_steps_frac_low=10):
self._init_logger()
self.company = company
self.n_latency_days = n_latency_days
self.hmm = GaussianHMM(n_components=n_hidden_states)
self._split_train_test_data(test_size)
self._compute_all_possible_outcomes(
n_steps_frac_change, n_steps_frac_high, n_steps_frac_low)
def _init_logger(self):
self._logger = logging.getLogger(__name__)
handler = logging.StreamHandler()
formatter = logging.Formatter(
'%(asctime)s %(name)-12s %(levelname)-8s %(message)s')
handler.setFormatter(formatter)
self._logger.addHandler(handler)
self._logger.setLevel(logging.DEBUG)
def _split_train_test_data(self, test_size):
data = pd.read_csv(
'data/company_data/{company}.csv'.format(company=self.company))
_train_data, test_data = train_test_split(
data, test_size=test_size, shuffle=False)
self._train_data = _train_data
self._test_data = test_data
@staticmethod
def _extract_features(data):
open_price = np.array(data['open'])
close_price = np.array(data['close'])
high_price = np.array(data['high'])
low_price = np.array(data['low'])
# Compute the fraction change in close, high and low prices
# which would be used a feature
frac_change = (close_price - open_price) / open_price
frac_high = (high_price - open_price) / open_price
frac_low = (open_price - low_price) / open_price
return np.column_stack((frac_change, frac_high, frac_low))
def fit(self):
self._logger.info('>>> Extracting Features')
feature_vector = StockPredictor._extract_features(self._train_data)
self._logger.info('Features extraction Completed <<<')
self.hmm.fit(feature_vector)
def _compute_all_possible_outcomes(self, n_steps_frac_change,
n_steps_frac_high, n_steps_frac_low):
frac_change_range = np.linspace(-0.1, 0.1, n_steps_frac_change)
frac_high_range = np.linspace(0, 0.1, n_steps_frac_high)
frac_low_range = np.linspace(0, 0.1, n_steps_frac_low)
self._possible_outcomes = np.array(list(itertools.product(
frac_change_range, frac_high_range, frac_low_range)))
def _get_most_probable_outcome(self, day_index):
previous_data_start_index = max(0, day_index - self.n_latency_days)
previous_data_end_index = max(0, day_index - 1)
previous_data = self._test_data.iloc[previous_data_end_index: previous_data_start_index]
previous_data_features = StockPredictor._extract_features(
previous_data)
outcome_score = []
for possible_outcome in self._possible_outcomes:
total_data = np.row_stack(
(previous_data_features, possible_outcome))
outcome_score.append(self.hmm.score(total_data))
most_probable_outcome = self._possible_outcomes[np.argmax(
outcome_score)]
return most_probable_outcome
def predict_close_price(self, day_index):
open_price = self._test_data.iloc[day_index]['open']
predicted_frac_change, _, _ = self._get_most_probable_outcome(
day_index)
return open_price * (1 + predicted_frac_change)
def predict_close_prices_for_days(self, days, with_plot=False):
predicted_close_prices = []
for day_index in tqdm(range(days)):
predicted_close_prices.append(self.predict_close_price(day_index))
if with_plot:
test_data = self._test_data[0: days]
days = np.array(test_data['date'], dtype="datetime64[ms]")
actual_close_prices = test_data['close']
fig = plt.figure()
axes = fig.add_subplot(111)
axes.plot(days, actual_close_prices, 'bo-', label="actual")
axes.plot(days, predicted_close_prices, 'r+-', label="predicted")
axes.set_title('{company}'.format(company=self.company))
fig.autofmt_xdate()
plt.legend()
plt.show()
return predicted_close_prices
stock_predictor = StockPredictor(company=args['--company'])
stock_predictor.fit()
stock_predictor.predict_close_prices_for_days(500, with_plot=True)
输出结果如下：
概要
在本章中，我们使用隐马尔可夫模型（HMM）预测了股票价格。我们应用了参数估计和模型评估方法来确定股票的收盘价格。使用 HMM 进行股市分析只是 HMM 应用于时间序列数据分析的另一个例子。
在下一章中，我们将探讨 HMM 在自然语言处理领域中的一个有趣应用。
第七章：自然语言处理
自动语音识别有许多潜在的应用，例如音频转录、听写、音频搜索和虚拟助手。我相信现在每个人至少与一个虚拟助手互动过，无论是苹果的 Siri、亚马逊的 Alexa，还是谷歌的 Assistant。在所有这些语音识别系统的核心，都是一组基于语言中不同单词或声音的统计模型。由于语音具有时间结构，隐马尔可夫模型（HMMs）是最自然的建模框架。
隐马尔可夫模型（HMMs）几乎是所有语音识别系统的核心，建模中的核心概念长时间没有太大变化。但随着时间的推移，已经开发出许多复杂的技术来构建更好的系统。在接下来的部分，我们将尝试涵盖推动这些系统发展的主要概念。
词性标注
我们将要研究的第一个问题被称为词性标注（POS 标注）。根据维基百科，词性标注，也叫语法标注或词类别消歧义，是根据一个单词的定义及其上下文（即它与短语、句子或段落中相邻和相关单词的关系）标记该单词为对应某一特定词性的过程。这个过程的一个简化版本，通常在学校里教授，就是将单词分类为名词、动词、形容词等。
词性标注并不像听起来那么简单，因为同一个单词在不同的上下文中可以有不同的词性。一个简单的例子是单词dogs。单词dogs通常被认为是名词，但在以下句子中，它充当了动词：
The sailor dogs the hatch。
正确的语法标注将反映出dogs在这里作为动词使用，而不是作为更常见的复数名词。语法上下文是确定这一点的一种方式；语义分析也可以用来推断sailor和hatch使dogs暗示为：
-
在航海术语中
-
作用于物体hatch的动作
在英语教学中，通常只教九大词性：名词、动词、冠词、形容词、介词、代词、副词、连词和感叹词。但我们可以将单词划分为更多类别和子类别，以便进行更精细的标注。例如，名词可以细分为复数、所有格和单数。同样，动词可以根据时态、体态等进行子分类。通常，基于计算机的词性标注系统能够区分 50 到 150 个英语词性。对科伊内希腊语的随机方法标注工作已经使用了超过 1,000 个词性，并发现那里有很多单词和英语一样存在歧义。
代码
对于代码示例，我们将使用pomegranate库来构建一个用于词性标注的 HMM。可以通过在命令行运行以下命令来安装 Pomegranate：
pip install pomegranate
在这个例子中，我们不深入讨论统计 POS 标注器。我们使用的数据是 Brown 语料库的副本。Brown 语料库包含 500 个英语文本样本，总计约 100 万个单词，编纂自 1961 年在美国出版的作品。
获取数据
我们先通过定义一些函数来读取corpus文件中的数据：
# Imports
import random
from itertools import chain
from collections import Counter, defaultdict
import numpy as np
import matplotlib.pyplot as plt
from pomegranate import State, HiddenMarkovModel, DiscreteDistribution
Sentence = namedtuple("Sentence", "words tags")
def read_data(filename):
"""
Function to read tagged sentence data.
Parameters
----------
filename: str
The path to the file from where to read the data.
"""
with open(filename, 'r') as f:
sentence_lines = [l.split("\n") for l in f.read().split("\n\n")]
return OrderedDict(((s[0], Sentence(*zip(*[l.strip().split("\t") for l in s[1:]])))
for s in sentence_lines if s[0])
def read_tags(filename):
"""
Function to read a list of word tag classes.
Parameters
----------
filename: str
The path to the file from where to read the tags.
"""
with open(filename, 'r') as f:
tags = f.read().split("\n")
return frozenset(tags)
现在让我们定义几个类，Subset和Dataset，以便更方便地处理数据：
class Subset(namedtuple("BaseSet", "sentences keys vocab X tagset Y N stream")):
"""
Class to handle a subset of the whole data. This is required when we split the
data into training and test sets.
"""
def __new__(cls, sentences, keys):
word_sequences = tuple([sentences[k].words for k in keys])
tag_sequences = tuple([sentences[k].tags for k in keys])
wordset = frozenset(chain(*word_sequences))
tagset = frozenset(chain(*tag_sequences))
N = sum(1 for _ in chain(*(sentences[k].words for k in keys)))
stream = tuple(zip(chain(*word_sequences), chain(*tag_sequences)))
return super().__new__(cls, {k: sentences[k] for k in keys}, keys, wordset,
word_sequences, tagset, tag_sequences, N, stream.__iter__)
def __len__(self):
return len(self.sentences)
def __iter__(self):
return iter(self.sentences.items())
class Dataset(namedtuple("_Dataset", "sentences keys vocab X tagset Y" +
"training_set testing_set N stream")):
"""
Class to represent the data in structured form for easy processing.
"""
def __new__(cls, tagfile, datafile, train_test_split=0.8, seed=112890):
tagset = read_tags(tagfile)
sentences = read_data(datafile)
keys = tuple(sentences.keys())
wordset = frozenset(chain(*[s.words for s in sentences.values()]))
word_sequences = tuple([sentences[k].words for k in keys])
tag_sequences = tuple([sentences[k].tags for k in keys])
N = sum(1 for _ in chain(*(s.words for s in sentences.values())))
# split data into train/test sets
_keys = list(keys)
if seed is not None:
random.seed(seed)
random.shuffle(_keys)
split = int(train_test_split * len(_keys))
training_data = Subset(sentences, _keys[:split])
testing_data = Subset(sentences, _keys[split:])
stream = tuple(zip(chain(*word_sequences), chain(*tag_sequences)))
return super().__new__(cls, dict(sentences), keys, wordset, word_sequences, tagset,
tag_sequences, training_data, testing_data, N, stream.__iter__)
def __len__(self):
return len(self.sentences)
def __iter__(self):
return iter(self.sentences.items())
现在，让我们尝试初始化Dataset类并查看它是如何工作的：
>>> data = Dataset("tags-universal.txt", "brown-universal.txt", train_test_split=0.8)
>>> print("There are {} sentences in the corpus.".format(len(data)))
There are 57340 sentences in the corpus.
>>> print("There are {} sentences in the training set.".format(len(data.training_set)))
There are 45872 sentences in the training set.
>>> print("There are {} sentences in the testing set.".format(len(data.testing_set)))
There are 11468 sentences in the testing set.
探索数据
现在，让我们探索数据，更好地理解我们的类如何存储信息。我们随机选择了b100-38532键：
>>> key = 'b100-38532'
>>> print("Sentence: {}".format(key))
Sentence: b100–38532
>>> print("words: {!s}".format(data.sentences[key].words))\
words: ('Perhaps', 'it', 'was', 'right', ';', ';')
>>> print("tags: {!s}".format(data.sentences[key].tags))
tags: ('ADV', 'PRON', 'VERB', 'ADJ', '.', '.')
我们还可以检查corpus中的唯一元素：
>>> print("There are a total of {} samples of {} unique words in the corpus.".format(
data.N, len(data.vocab)))
There are a total of 1161192 samples of 56057 unique words in the corpus.
>>> print("There are {} samples of {} unique words in the training set.".format(
data.training_set.N, len(data.training_set.vocab)))
There are 928458 samples of 50536 unique words in the training set.
>>> print("There are {} samples of {} unique words in the testing set.".format(
data.testing_set.N, len(data.testing_set.vocab)))
There are 232734 samples of 25112 unique words in the testing set.
>>> print("There are {} words in the test set that are missing in the training set.".format(
len(data.testing_set.vocab - data.training_set.vocab)))
There are 5521 words in the test set that are missing in the training set.
我们还可以使用Dataset类的X和Y属性来访问单词及其对应的标签：
>>> for i in range(2):
... print("Sentence {}:".format(i + 1), data.X[i])
... print("Labels {}:".format(i + 1), data.Y[i], "\n")
Sentence 1: ('Mr.', 'Podger', ‘had', 'thanked', 'him', 'gravely', ',', 'and', 'now', 'he', 'made', 'use', 'of', 'the', 'advice', '.')
Labels 1: ('NOUN', 'NOUN', 'VERB', 'VERB', 'PRON', 'ADV', '.', 'CONJ', 'ADV', 'PRON', 'VERB', 'NOUN', 'ADP', 'DET', 'NOUN', '.')
Sentence 2: ('But', 'there', 'seemed', 'to', 'be', 'some', 'difference', 'of', 'opinion', 'as', 'to', 'how', 'far', 'the', 'board', 'should', 'go', ',', 'and', 'whose', 'advice', 'it', 'should', 'follow', '.')
Labels 2: ('CONJ', 'PRT', 'VERB', 'PRT', 'VERB', 'DET', 'NOUN', 'ADP', 'NOUN', 'ADP', 'ADP', 'ADV', 'ADV', 'DET', 'NOUN', 'VERB','VERB', '.', 'CONJ', 'DET', 'NOUN', 'PRON', 'VERB', 'VERB', '.')
我们还可以使用stream方法来遍历单词和标签的配对：
>>> for i, pair in enumerate(data.stream()):
... print(pair)
... if i > 3:
... break
('Podger', 'NOUN')
('had', 'VERB')
('thanked', 'VERB')
('him', 'PRON')
查找最频繁的标签
现在，为了比较我们 HMM 模型的性能，让我们构建一个最频繁类标注器（MFC Tagger）。我们从定义一个函数开始，来统计标签和单词的配对：
def pair_counts(tags, words):
d = defaultdict(lambda: defaultdict(int))
for tag, word in zip(tags, words):
d[tag][word] += 1
return d
tags = [tag for i, (word, tag) in enumerate(data.training_set.stream())]
words = [word for i, (word, tag) in enumerate(data.training_set.stream())]
现在，让我们定义MFCTagger类：
FakeState = namedtuple('FakeState', 'name')
class MFCTagger:
missing = FakeState(name = '<MISSING>')
def __init__(self, table):
self.table = defaultdict(lambda: MFCTagger.missing)
self.table.update({word: FakeState(name=tag) for word, tag in table.items()})
def viterbi(self, seq):
"""This method simplifies predictions by matching the Pomegranate viterbi() interface"""
return 0., list(enumerate(["<start>"] + [self.table[w] for w in seq] + ["<end>"]))
tags = [tag for i, (word, tag) in enumerate(data.training_set.stream())]
words = [word for i, (word, tag) in enumerate(data.training_set.stream())]
word_counts = pair_counts(words, tags)
mfc_table = dict((word, max(tags.keys(), key=lambda key: tags[key])) for word, tags in word_counts.items())
mfc_model = MFCTagger(mfc_table)
这里有一些辅助函数，用于从模型中进行预测：
def replace_unknown(sequence):
return [w if w in data.training_set.vocab else 'nan' for w in sequence]
def simplify_decoding(X, model):
_, state_path = model.viterbi(replace_unknown(X))
return [state[1].name for state in state_path[1:-1]]
>>> for key in data.testing_set.keys[:2]:
... print("Sentence Key: {}\n".format(key))
... print("Predicted labels:\n-----------------")
... print(simplify_decoding(data.sentences[key].words, mfc_model))
... print()
... print("Actual labels:\n--------------")
... print(data.sentences[key].tags)
... print("\n")
评估模型准确性
为了检查我们模型的表现，让我们评估模型的准确率：
def accuracy(X, Y, model):
correct = total_predictions = 0
for observations, actual_tags in zip(X, Y):
# The model.viterbi call in simplify_decoding will return None if the HMM
# raises an error (for example, if a test sentence contains a word that
# is out of vocabulary for the training set). Any exception counts the
# full sentence as an error (which makes this a conservative estimate).
try:
most_likely_tags = simplify_decoding(observations, model)
correct += sum(p == t for p, t in zip(most_likely_tags, actual_tags))
except:
pass
total_predictions += len(observations)
return correct / total_predictions
>>> mfc_training_acc = accuracy(data.training_set.X, data.training_set.Y, mfc_model)
>>> print("Training accuracy mfc_model: {:.2f}%".format(100 * mfc_training_acc))
Training accuracy mfc_model: 95.72%
>>> mfc_testing_acc = accuracy(data.testing_set.X, data.testing_set.Y, mfc_model)
>>> print("Testing accuracy mfc_model: {:.2f}%".format(100 * mfc_testing_acc))
Testing accuracy mfc_model: 93.01%
基于 HMM 的标注器
现在，我们将尝试使用 HMM 构建一个 POS 标注器，并希望它能提高我们的预测性能。我们将首先定义一些辅助函数：
def unigram_counts(sequences):
return Counter(sequences)
tags = [tag for i, (word, tag) in enumerate(data.training_set.stream())]
tag_unigrams = unigram_counts(tags)
def bigram_counts(sequences):
d = Counter(sequences)
return d
tags = [tag for i, (word, tag) in enumerate(data.stream())]
o = [(tags[i],tags[i+1]) for i in range(0,len(tags)-2,2)]
tag_bigrams = bigram_counts(o)
def starting_counts(sequences):
d = Counter(sequences)
return d
tags = [tag for i, (word, tag) in enumerate(data.stream())]
starts_tag = [i[0] for i in data.Y]
tag_starts = starting_counts(starts_tag)
def ending_counts(sequences):
d = Counter(sequences)
return d
end_tag = [i[len(i)-1] for i in data.Y]
tag_ends = ending_counts(end_tag)
现在让我们构建模型：
basic_model = HiddenMarkovModel(name="base-hmm-tagger")
tags = [tag for i, (word, tag) in enumerate(data.stream())]
words = [word for i, (word, tag) in enumerate(data.stream())]
tags_count=unigram_counts(tags)
tag_words_count=pair_counts(tags,words)
starting_tag_list=[i[0] for i in data.Y]
ending_tag_list=[i[-1] for i in data.Y]
starting_tag_count=starting_counts(starting_tag_list)#the number of times a tag occurred at the start
ending_tag_count=ending_counts(ending_tag_list) #the number of times a tag occurred at the end
to_pass_states = []
for tag, words_dict in tag_words_count.items():
total = float(sum(words_dict.values()))
distribution = {word: count/total for word, count in words_dict.items()}
tag_emissions = DiscreteDistribution(distribution)
tag_state = State(tag_emissions, name=tag)
to_pass_states.append(tag_state)
basic_model.add_states()
start_prob={}
for tag in tags:
start_prob[tag]=starting_tag_count[tag]/tags_count[tag]
for tag_state in to_pass_states :
basic_model.add_transition(basic_model.start,tag_state,start_prob[tag_state.name])
end_prob={}
for tag in tags:
end_prob[tag]=ending_tag_count[tag]/tags_count[tag]
for tag_state in to_pass_states :
basic_model.add_transition(tag_state,basic_model.end,end_prob[tag_state.name])
transition_prob_pair={}
for key in tag_bigrams.keys():
transition_prob_pair[key]=tag_bigrams.get(key)/tags_count[key[0]]
for tag_state in to_pass_states :
for next_tag_state in to_pass_states :
basic_model.add_transition(tag_state,next_tag_state,transition_prob_pair[(tag_state.name,next_tag_state.name)])
basic_model.bake()
>>> for key in data.testing_set.keys[:2]:
... print("Sentence Key: {}\n".format(key))
... print("Predicted labels:\n-----------------")
... print(simplify_decoding(data.sentences[key].words, basic_model))
... print()
... print("Actual labels:\n--------------")
... print(data.sentences[key].tags)
... print("\n")
>>> hmm_training_acc = accuracy(data.training_set.X, data.training_set.Y, basic_model)
>>> print("Training accuracy basic hmm model: {:.2f}%".format(100 * hmm_training_acc))
Training accuracy basic hmm model: 97.49%
>>> hmm_testing_acc = accuracy(data.testing_set.X, data.testing_set.Y, basic_model)
>>> print("Testing accuracy basic hmm model: {:.2f}%".format(100 * hmm_testing_acc))
Testing accuracy basic hmm model: 96.09%
我们可以看到，基于 HMM 的模型已经能够提高我们模型的准确率。
语音识别
在 1950 年代，贝尔实验室是语音识别的先锋。早期设计的系统仅限于单一发言者，并且词汇量非常有限。经过约 70 年的努力，当前的语音识别系统能够处理来自多个发言者的语音，并且能够识别数千个词汇，涵盖多种语言。由于每种技术都有足够的工作，已经足够成书，本书并不深入讨论所有使用的技术。
但是，语音识别系统的通用工作流程是首先通过麦克风将物理声音转换为电信号来捕捉音频。麦克风生成的电信号是模拟的，需要通过模拟-数字转换器将其转换为数字形式，以便存储和处理。一旦我们得到了数字化的语音，就可以应用算法来理解这些语音。
如前所述，大多数先进的语音识别系统仍然使用隐马尔可夫模型（HMM）作为核心。其基本假设是，语音信号在短时间内（几毫秒）是一个平稳的过程。因此，语音识别系统的第一步是将信号拆分为大约 10 毫秒的小片段。然后，每个片段的功率谱被映射到一个实数向量，这个向量被称为倒谱系数。这个向量的维度通常较小，尽管更精确的系统通常会使用超过 32 维的向量。HMM 模型的最终输出是一系列这些向量。一旦我们获得这些向量，这些向量组会与一个或多个音素匹配，音素是语音的基本单位。但为了有效地将这些向量组与音素匹配，我们需要训练系统，因为不同发音者以及同一发音者不同发音之间，音素的发音差异巨大。一旦我们得到音素的序列，系统会尝试猜测最有可能产生该音素序列的词。
如我们所想，这整个检测过程计算开销非常大。为了解决这一复杂性问题，现代语音识别系统使用神经网络进行特征转换和降维，然后再使用 HMM 进行识别。另一个常用的减少计算量的技术是使用语音活动检测器，它可以检测信号中包含语音的区域。利用这些信息，我们可以设计识别器只对包含语音的信号部分进行计算。
幸运的是，Python 拥有一个非常发达的生态系统，可以与语音识别技术进行配合。在下一部分中，我们将探讨用于语音识别的不同 Python 包。
语音识别的 Python 包
Python 包托管服务 PyPI 列出了许多语音识别系统。以下是一些最常用的：
-
SpeechRecognition (
pypi.org/project/SpeechRecognition/) -
apiai (
pypi.org/project/apiai/) -
assemblyai (
pypi.org/project/assemblyai/) -
pocketsphinx (
pypi.org/project/pocketsphinx/) -
google-cloud-speech (
pypi.org/project/google-cloud-speech/) -
watson-developer-cloud (
pypi.org/project/watson-developer-cloud/)
一些 Python 包（如 apiai）不仅提供语音识别，还实现了自然语言处理算法，用户可以通过这些算法从语音中识别出说话者的意图。其他包则专注于语音识别，能够将音频转换为文本。
在本章中，我们将使用 SpeechRecognition 包。我们选择 SpeechRecognition 有两个原因：
-
它提供了一个非常简单的 API，可以直接访问和处理音频信号。对于其他包，通常需要编写小脚本才能访问文件。
-
它是多个流行语音 API 的封装，因此非常灵活，且可以在不对代码做太多修改的情况下使用多个服务。
因此，要开始使用 SpeechRecognition，我们需要安装该包。由于它托管在 PyPI 上，因此可以通过 pip 直接安装：
pip install SpeechRecognition
语音识别基础
该包中最重要的类是 Recognizer 类，因为它处理大多数的识别任务。我们可以在初始化类时指定不同的设置和功能，以便从音频源中识别语音。
Recognizer 类可以非常简单地初始化，无需传递任何参数：
>>> import speech_recognition as sr
>>> r = sr.Recognizer()
Recognizer 类的每个实例有七种不同的方法可以用于将语音转换为文本。这些方法分别使用了特定的语音识别服务。七种方法如下：
-
recognize_bing：使用微软的 Bing 语音服务。 -
recognize_google：使用 Google 的 Web Speech API。 -
recognize_google_cloud：使用 Google 的 Cloud Speech。使用此方法需要安装google-cloud-speech，可以通过运行pip install google-cloud-speech来轻松安装。 -
recognize_houndify：使用 SoundHound 的 Houndify 服务。 -
recognize_ibm：使用 IBM 的语音转文本服务。 -
recognize_sphinx：使用 CMU 的 Sphinx。此方法依赖于PocketSphinx，可以通过运行pip install pocketsphinx来安装。 -
recognize_wit：使用 Wit.ai 服务。
使用这些方法时要记住的一个重要事项是，由于大多数语音识别服务是通过 Web API 由公司提供的，因此我们需要互联网连接才能访问这些服务。此外，一些服务只有在用户在线注册后才允许使用。在这七种方法中，只有 recognize_sphinx 可以离线工作。
在所有这些 Web API 中，只有 Google 的 Web Speech API 在没有注册或 API 密钥的情况下即可使用。因此，为了简化起见，我们将在本章剩余部分使用它。
如果服务器不可用、没有互联网连接或 API 配额限制已达，recognize 方法会抛出 RequestError。
接下来我们需要做任何识别的，是一些音频数据。SpeechRecognition提供了直接的功能，可以使用音频文件或附加的麦克风音频。在接下来的部分中，我们将探讨这两种方法。
从音频文件中进行语音识别
要开始使用音频文件，首先需要下载一个。对于下面的示例，我们将使用harvard.wav文件，它可以从raw.githubusercontent.com/realpython/python-speech-recognition/master/audio_files/harvard.wav下载。
确保将音频文件保存在 Python 解释器运行的同一目录中。否则，在接下来的代码中，文件的路径需要进行修改。
对于音频文件的处理，SpeechRecognition提供了AudioFile类，用于读取和操作音频文件。我们可以使用AudioFile的record方法来处理音频文件的内容，然后再与Recognizer类一起使用：
>>> harvard = sr.AudioFile('harvard.wav')
>>> with harvard as source:
... audio = r.record(source)
上下文管理器打开音频并将其内容记录到audio中，audio是AudioFile的一个实例。我们可以通过调用type方法检查它：
>>> type(audio)
<class 'speech_recognition.AudioData'>
现在，一旦我们拥有了一个AudioFile实例，就可以调用任何识别方法，并将其作为参数传递。识别方法会调用特定的网络 API，将音频文件中的语音转换成文本，并返回如下文本：
>>> r.recognize_google(audio)
'the stale smell of old beer lingers it takes heat
to bring out the odor a cold dip restores health and
zest a salt pickle taste fine with ham tacos al
Pastore are my favorite a zestful food is the hot
cross bun'
在这个案例中，我们转录了整个音频文件，但如果我们只想翻译音频文件的特定部分，该怎么办？可以通过将额外的参数传递给record方法来实现：
>>> with harvard as source:
... audio_part = r.record(source, duration=4)
>>> r.recognize_google(audio_part)
'the stale smell of old beer lingers'
record方法会在音频文件中保留一个指针，指向记录已经完成的位置。因此，如果我们在同一个文件上再录制四秒钟，它将从原始音频文件的四秒标记到八秒标记进行录音。
在前面的示例中，我们转录了音频文件的一部分，但起始点是文件的开始。如果我们希望从不同的时间点开始该怎么办？可以通过将另一个参数offset传递给record方法来实现：
>>> with harvard as source:
... audio_offset = r.record(source, offset=4, duration=3)
>>> recognizer.recognize_google(audio_offset)
'it takes heat to bring out the odor'
如果你听一下harvard.wav文件，你会发现录音是在没有外部噪音的完美环境下完成的，但在现实生活中的音频通常不是这样。我们来试着转录另一个音频信号jackhammer.wav，它可以从raw.githubusercontent.com/realpython/python-speech-recognition/master/audio_files/jackhammer.wav下载。如果你听这个音频文件，你会注意到它有很多背景噪音。让我们试着转录这个文件，看看识别器的表现如何：
>>> jackhammer = sr.AudioFile('jackhammer.wav')
>>> with jackhammer as source:
... audio_jack = r.record(source)
>>> r.recognize_google(audio_jack)
'the snail smell of old gear vendors'
如我们所见，转录结果差得很远。在这种情况下，我们可以使用Recognizer类中提供的adjust_for_ambient_noise方法来调整我们的音频信号与噪声的匹配。
adjust_for_ambient_noise默认使用前一秒的数据进行校准，但我们可以通过传递duration参数来改变这一点：
>>> with jackhammer as source:
... r.adjust_for_ambient_noise(source, duration=1)
... audio = r.record(source)
>>> r.recognize_google(audio)
'still smell of old beer vendors'
如果我们不想失去太多信息，我们可以减小duration参数的值，但这可能会导致较差的校准。如我们所见，转录结果仍然不完美，但比我们没有使用adjust_for_ambient_noise时要好得多。实际上，我们可以通过使用信号处理技术来清理音频中的噪声，从而获得更好的结果，但这超出了本书的范围。
在这种情况下，我们还可以做的一件事是查看识别器最可能的所有转录。可以通过在调用recognize方法时使用show_all参数来实现：
>>> r.recognize_google(audio, show_all=True)
{'alternative': [
{'transcript': 'the snail smell like old Beer Mongers'},
{'transcript': 'the still smell of old beer vendors'},
{'transcript': 'the snail smell like old beer vendors'},
{'transcript': 'the stale smell of old beer vendors'},
{'transcript': 'the snail smell like old beermongers'},
{'transcript': 'destihl smell of old beer vendors'},
{'transcript': 'the still smell like old beer vendors'},
{'transcript': 'bastille smell of old beer vendors'},
{'transcript': 'the still smell like old beermongers'},
{'transcript': 'the still smell of old beer venders'},
{'transcript': 'the still smelling old beer vendors'},
{'transcript': 'musty smell of old beer vendors'},
{'transcript': 'the still smell of old beer vendor'}
], 'final': True}
使用这个方法，我们可以选择最适合我们具体问题的转录结果。
使用麦克风进行语音识别
在上一节中，我们使用了识别器方法从音频文件中转录语音。在本节中，我们将使用从麦克风录制的语音进行转录。
但在我们开始之前，我们需要安装一个额外的软件包，叫做PyAudio。它也可以在 PyPI 上找到，并且可以直接使用pip: pip install PyAudio进行安装。
在上一节中，处理音频文件时我们使用了AudioFile类，但在处理麦克风时，我们需要使用Microphone类。大部分的识别 API 仍然保持不变。让我们通过一个简单的例子来理解它是如何工作的：
>>> import speech_recognition as sr
>>> r = sr.Recognizer()
>>> mic = sr.Microphone()
>>> with mic as source:
... r.adjust_for_ambient_noise(source)
... audio = r.listen(source)
与上一节初始化AudioFile类类似，这一次我们需要初始化Microphone类。此外，我们需要调用listen方法来录制音频，而不是record。当执行前面的代码块时，Python 解释器会等待一段时间来录制音频。试着对着麦克风说话。一旦解释器提示符返回，我们就可以调用recognize方法来转录录制的音频：
>>> r.recognize_google(audio)
'hello world' # This would depend on what you said in the microphone
总结
在这一章中，我们探讨了 HMM 的两个主要应用：词性标注和语音识别。我们使用最频繁标签算法编写了词性标注器，并使用pomegranate包基于 HMM 构建了一个词性标注器。我们比较了这两种方法的性能，发现基于 HMM 的方法优于最频繁标签法。然后，我们使用SpeechRecognition包通过 Google 的 Web 语音 API 将音频转录为文本。我们探讨了如何使用该包处理音频文件和来自麦克风的实时音频。
在下一章中，我们将探索 HMM 的更多应用，特别是在图像识别领域。
第八章：用于图像处理的二维隐马尔可夫模型（2D HMM）
本章我们将介绍隐马尔可夫模型（HMM）在图像分割中的应用。在图像分割中，我们通常将给定的图像拆分成多个相同大小的块，然后对每个块进行估计。然而，这些算法通常忽略了来自邻近块的上下文信息。为了解决这个问题，引入了二维隐马尔可夫模型（2D HMM），它通过一个潜在的二维马尔可夫网格考虑特征向量之间的依赖关系。在本章中，我们将讨论二维隐马尔可夫模型是如何工作的，并推导出相应的参数估计算法。本章将讨论以下主题：
-
假二维隐马尔可夫模型
-
二维隐马尔可夫模型（2D HMM）简介
-
二维隐马尔可夫模型中的参数学习
-
应用
一维隐马尔可夫模型回顾
让我们回顾一下前几章中讨论的如何使用一维隐马尔可夫模型（1D HMM）。我们看到，隐马尔可夫模型其实是马尔可夫链上的一个过程。在任何时刻，隐马尔可夫模型处于某个可能的状态，模型将转移到的下一个状态取决于当前状态以及模型的转移概率。
假设对于隐马尔可夫模型（HMM）有 M = {1, 2, ..., M} 个可能的状态，且从某个状态 i 转移到状态 j 的转移概率由 a[i,j] 给出。对于这样的模型，如果在时刻 t-1 模型处于状态 i，那么在时刻 t 它将以概率 a[i,j] 进入状态 j。这个概率被称为转移概率。另外，我们已经定义了模型中的观察变量，它仅依赖于隐藏变量的当前状态。我们可以将时刻 t 的观察变量定义为 u[t]，假设状态 i 对变量 u[t] 的发射分布为 bi。
我们还需要定义初始概率 π[i]，即在时刻 t = 1 时处于状态 i 的概率。通过给定这些值，我们可以确定观察到任何给定序列的似然，
，如下所示：
在大多数情况下，我们假设状态为高斯混合模型；在这种情况下，前面的方程可以进一步推广，如下所示：
在这里，u[i] 是均值，∑[i] 是协方差矩阵，k 是观察变量 u[t] 的维度。
现在我们已经定义了模型，可以继续进行估计方法。估计通常通过我们在第四章中看到的 Baum-Welch 算法来进行，最大似然参数推断，该算法执行最大似然估计。令 Li 为在时刻 t 处于状态 i 的条件分布，给定观察结果，Hi,j 为从状态 i 在时刻 t + 1 转移到状态 j 的条件概率，仍然是给定观察结果。然后，我们可以重新估计均值、协方差和转移概率，如下所示：
为了计算 Li 和 Hi,j 的值，我们使用前向-后向算法。前向算法给出了观察到前 t 个结果的概率，αi，以及在时刻 t 时处于状态 i 的概率。
这个概率可以通过以下方程组进行评估：
我们还定义了后向概率，βi，作为给定模型在时刻 t 处于状态 i 时，观察到 {u[r]}r=t + 1,...,T 的条件概率。后向概率 βi 可以通过以下方式计算：
使用这些值后，我们可以计算出 Li 和 Hi,j 可以通过以下方式求解：
我们可以通过假设每个观察结果来自最可能的隐藏状态来近似该算法。这使得我们可以简化 Baum-Welch 算法；这种方法通常也被称为 维特比训练算法。给定观察到的状态，并假设状态序列为 
，其具有最大条件概率，可以通过以下方式表示：
这里，
可以通过以下方式计算：
使用这些值，我们可以计算模型参数，如下所示：
这里，I 是指示函数，当函数的参数为真时返回 1，否则返回 0。
在本节中，我们简要回顾了当状态使用高斯混合模型参数化时的一维 HMM（隐马尔可夫模型）的基本概念，现在我们可以继续讨论二维 HMM。
二维 HMM
关于二维 HMM，已经做了很多工作，但最近的、受欢迎的工作是由 Jia Li、Amir Najmi 和 Robert Gray 在他们的论文《通过二维隐马尔可夫模型进行图像分类》中完成的。本节内容基于他们的工作编写。我们将从介绍他们提出的通用算法开始，随后在接下来的小节中，我们将看到该算法的具体运作方式。
算法
图像分类的算法如下：
-
训练：
-
将训练图像划分为大小相等且不重叠的块，并为每个块提取特征向量
-
选择 2D HMM 的状态数
-
根据特征向量和训练标签估计模型参数
-
-
测试：
-
与训练类似，为测试图像生成特征向量
-
根据训练模型，给定特征向量，搜索具有最大后验概率的类别集合
-
2D HMM 模型的假设
本节我们将快速浏览 2D HMM 模型的假设以及这些假设如何简化我们的方程。更详细的推导请参阅原始论文。
我们从将图像划分为较小的块开始，从中评估特征向量，并使用这些特征向量对图像进行分类。在 2D HMM 的情况下，我们假设特征向量是由一个马尔可夫模型生成的，并且每次块的变化都会发生状态变化。我们还基于块的前后关系定义块之间的关系。若位置为 (i', j') 的块出现在位置 (i, j) 之前，则满足 i' 或 i' = i 且 j' < j。假设有 M = {1, 2, ... M} 个状态，任何给定位置 (i, j) 的块的状态用 S[i,j] 表示，特征向量用 u[i,j] 表示，类别用 c[i,j] 表示。另一个需要注意的点是，块的顺序只是为了说明模型的假设，算法在进行分类时并不考虑块的顺序。
模型做出的第一个假设如下：
该假设声明，知道当前状态足以估计转移概率，这意味着 u 是条件独立的。同时，根据该假设，状态转移是一个二维的马尔可夫过程，系统进入任何特定状态的概率取决于前一时间和观察实例中模型在水平和垂直方向上的状态。我们还假设状态与类别之间是 1 对 1 的映射，因此一旦知道状态，就可以直接计算类别。
第二个假设是特征向量对于每个状态来说是一个高斯混合分布。我们知道，任何 M 成分的高斯混合可以分解为 M 个单一高斯分布的子状态；因此，对于状态 s 和特征向量 u 的块，分布的概率由以下公式给出：
这里，∑[s] 是协方差矩阵，μ[s] 是高斯分布的均值向量。
我们现在可以使用马尔可夫假设来简化状态的概率评估。图像中所有块的状态概率用 P{s[i,j] : (i,j) ∈ N} 表示，其中 N = {(i,j) : 0 ≤ i < w, 0 ≤ j < z}。但在我们使用马尔可夫假设有效地展开概率之前，我们需要证明，基于前两个假设，图像满足二维马尔可夫属性的旋转形式。
我们定义了 
的旋转关系，记为 
，它指定了 
，如果 j' < j，或 j' = j 且 i' < i。我们需要证明以下内容：
因此，我们使用前面定义的 
并引入以下新符号：
从前面的方程中，我们还可以看到这一点：
现在，为了简化符号，假设我们引入 
和 
。从这些符号中，我们可以进行以下推导：
展开条件概率，我们得到如下结果：
使用马尔可夫假设，我们得到如下结果：
最后，使用马尔可夫假设以及假设给定状态下，块的特征向量是条件独立的，我们得到以下结果：
其中 m = s[i-1,j]，n = s[i,j-1]，和 l = s[i,j]。
如果我们将 
替换为 
，并将 
替换为 
在推导过程中，所有的方程仍然成立。这给我们带来了以下结果：
由于前面的方程包含了原始的马尔可夫假设及其旋转版本，我们可以如下展示这两个假设的等价性：
现在我们可以简化 P{s[i,j] : (i,j) ∈ N} 的展开式，如下所示：
这里，w 和 z 分别是图像的行数和列数，T[i] 是对角线 i 上块的状态序列。接下来，我们需要证明 P(T[i]|T[i-1],...,T[0]) = P(T[i]|T[i-1])。假设 T[i] = {s[i,0], s[i-1,1], ..., s[0,i]}，这意味着 T[i-1] = {s[i-1,0], s[i-2,1], ..., s[0,i-1]}，因此我们得到以下结果：
因此，我们可以得出结论：
使用此方法，我们得到以下简化的方程：
使用 EM 进行参数估计
由于我们已经准备好模型，现在需要估计模型的参数。我们需要估计均值向量μ[m]；协方差矩阵∑[m]；以及转移概率a[m,n,l]，其中m,n,l = 1,..., M，M是状态的总数。我们将使用期望最大化（EM）算法。
正如我们在前面的章节中看到的，EM 是一个迭代算法，可以在缺失数据的情况下学习最大似然估计；也就是说，当我们的数据中有未观察到的变量时。假设我们的未观察到的变量x位于样本空间x中，而观察到的变量y位于样本空间y中。如果我们假设一个分布族f(x|Φ)，其参数Φ ∈ Ω，那么y的分布如下所示：
EM 算法将尝试找到一个值！，该值会最大化给定观察到的y时，g(y|Φ)的值。EM 迭代Φ^((p)) → Φ^((p+1))定义如下：
-
E 步：计算Q(Φ|**Φ^((p)))，其中Q(Φ'|Φ)是log f(x|Φ')的期望值。
-
M 步：选择Φ^((p+1))为Φ ∈ Ω中的一个值，以最大化Q(Φ|**Φ^((p)))。
现在让我们定义二维 HMM 所需的术语：
现在，让我们定义一个二维 HMM 中x的分布，如下所示：
从这个，我们可以计算出以下log f(x|Φ')：
我们现在知道，给定y，x只能有有限个与y的值一致的状态。因此，x的分布如下所示：
这里，α 是归一化常数，I 是指示函数。现在，对于 M 步，我们需要将Φ^((p+1))的值设为Φ'，以最大化以下公式：
由于前面的项有两部分相加，我们可以分别处理每一项，因为我们正在尝试最大化总值。考虑第一项：
该项在
中是凹函数；因此，使用拉格朗日乘数并对其求导，我们得到以下结果：
回到最大化
的第二项：
再次为了简化前面的方程，我们定义
，前面的项变成了如下形式：
在这种情况下，我们的高斯分布的最大似然估计（ML 估计）如下所示：
为了总结整个推导过程，我们可以将 EM 算法写成以下两个步骤：
- 给定模型估计Φ^((p)), 参数更新如下所示：
这里，项
可以按以下方式计算：
- 转移概率的更新如下所示：
这里，
的计算方法如下：
通过迭代应用前面两个方程，我们的算法将收敛到模型参数的最大似然估计。
总结
本章我们首先回顾了上一章介绍的 1D HMMs（隐马尔可夫模型）。随后，我们介绍了 2D HMMs（二维隐马尔可夫模型）的概念，并推导了为了简化计算，我们对 2D HMMs 所做的各种假设，这些假设可以应用于图像处理任务。接着，我们介绍了一种通用的 EM 算法，用于学习 2D-HMMs 中的参数。
在下一章中，我们将探讨 HMMs 在强化学习领域中的另一种应用，并介绍 MDP（马尔可夫决策过程）。
第九章：马尔可夫决策过程
在本章中，我们将讨论隐马尔可夫模型（HMMs）的另一种应用，即马尔可夫决策过程（MDP）。在 MDP 的情况下，我们为模型引入了奖励，任何由过程经历的状态序列都会产生特定的奖励。我们还将引入折扣的概念，这将使我们能够控制代理的短视或远视程度。代理的目标是最大化它能够获得的总奖励。
在本章中，我们将讨论以下主题：
-
强化学习
-
马尔可夫奖励过程
-
马尔可夫决策过程
-
代码示例
强化学习
强化学习是机器学习中的一种不同范式，其中代理通过做出决策/行动并观察结果来学习在定义的环境中如何最优地行为。因此，在强化学习中，代理并不是从某个给定的数据集中学习，而是通过与环境互动，代理尝试通过观察其行为的效果来学习。环境的定义方式使得代理在其行动使其接近目标时会获得奖励。
人类就是以这种方式学习的。例如，考虑一个站在壁炉前的孩子，孩子是代理，周围的空间是环境。现在，如果孩子把手靠近火炉，他会感受到温暖，这种感觉是愉快的，某种程度上，孩子（或代理）因将手靠近火炉而获得奖励。但如果孩子将手移得太靠近火炉，手就会被烧伤，因此会收到负奖励。通过这些奖励，孩子能够找出保持手与火炉之间的最优距离。强化学习正是模仿这种系统，训练代理学习如何在给定的环境中优化其目标。
更正式地说，要训练一个代理，我们需要有一个环境，这个环境代表了代理应该能够采取行动的世界。对于这些行动，环境应该返回包含奖励信息的观察结果，告诉代理这次行动的好坏。观察结果还应包含关于代理在环境中下一状态的信息。基于这些观察结果，代理尝试找出到达目标的最优方式。图 1 显示了代理与环境之间的互动：
强化学习与其他算法根本不同的地方在于没有监督者。相反，只有一个奖励信号，反馈代理它所采取的行动。这里还需要提到一个重要的事情，那就是环境可以被构建成奖励是延迟的，这可能会导致代理在达到目标之前四处游荡。也有可能代理在达到目标之前要经历大量的负反馈。
在监督学习中，我们给定的数据集基本上会告诉我们的学习算法在不同情况下的正确答案。然后，我们的学习算法会查看所有这些不同的情况以及这些情况下的解决方案，并试图基于这些情况进行归纳。因此，我们也期望给定的数据集是独立同分布（IID）的。但在强化学习中，数据不是 IID，生成的数据取决于代理所走的路径，因此取决于代理采取的行动。因此，强化学习是一个主动学习过程，代理所采取的行动影响环境，进而影响环境生成的数据。
我们可以通过一个非常简单的例子来更好地理解强化学习代理和环境的行为。考虑一个代理试图学习玩超级马里奥兄弟：
-
代理将从环境中接收初始状态。以超级马里奥为例，这将是当前的游戏画面。
-
在接收到状态信息后，代理将尝试采取行动。假设代理采取的行动是向右移动。
-
当环境接收到这个行动时，它将基于此返回下一个状态。下一个状态也会是一个画面，但如果前一个状态中马里奥旁边有敌人，那么画面可能会显示马里奥死亡。否则，画面只会显示马里奥向右移动了一步。环境还会根据行动返回奖励。如果马里奥右边有敌人，奖励可能是-5（因为该行动导致马里奥死亡），或者如果马里奥朝着完成关卡的方向移动，奖励可能是+1。
奖励假设
强化学习的整个概念是基于一种叫做奖励假设的东西。根据奖励假设，所有目标都可以通过最大化预期的累积奖励来定义。
在任意给定的时间点 t，总的累积奖励可以表示为：
但实际上，离当前状态较近的奖励比远离的奖励更可能发生。为了解决这个问题，我们引入了另一个术语，称为折扣率，
。折扣率的值始终在 0 和 1 之间。较大的折扣率
意味着较小的折扣，这使得智能体更关心长期奖励；而较小的折扣率
则使得智能体更关心短期奖励。通过引入折扣率项，我们现在可以定义我们的累积奖励为：
环境和智能体的状态
正如我们之前所看到的，智能体与环境在时间间隔t = 0, 1, 2, 3, .., 进行交互，在每个时间点，智能体根据环境的状态S[t]采取动作A[t]并获得奖励R[t+1]。这种随时间变化的状态、动作和奖励的序列被称为智能体的历史，表示为：
理想情况下，我们希望智能体所采取的动作能够基于其所有历史，但通常不可行，因为智能体的历史可能非常庞大。因此，我们定义状态的方式是将其作为智能体历史的总结：
基本上，我们将状态定义为智能体历史的一个函数。需要注意的是，环境状态是环境用来根据动作决定其下一个状态的状态，并且给出奖励。此外，环境状态是环境的私有信息。
另一方面，智能体状态是智能体用来决定下一步动作的状态。智能体的状态是其内部表示，可以是前述历史的任何函数。
我们使用马尔可夫状态来表示智能体的状态，这基本上意味着智能体的当前状态能够总结所有历史，智能体的下一步动作只会依赖于当前的状态。因此，
在本章中，我们仅考虑智能体能够直接观察环境状态的情况。这导致环境的观察既是智能体的当前状态，也是环境的状态。这种特殊情况通常被称为MDP。
智能体的组件
在本节中，我们将正式定义智能体可能具有的不同类型的组件。
-
策略：它是给定状态下对动作的条件概率分布。根据这个条件分布，智能体在任何给定状态下选择其动作。策略可以是确定性的：a = π(s)，也可以是随机的：
![]()
。 -
价值函数：价值函数试图预测在给定状态下采取某个动作所期望的奖励。它的表达式为：
![]()
。其中 E 是期望值，
![]()
是折扣因子。 -
模型：模型是智能体对环境的表征。它使用一个转移函数 P 来定义，预测环境的下一个状态：
![]()
。以及一个奖励函数，用于预测在给定状态下采取某个动作所关联的奖励：
![]()
。
。
。
是折扣因子。
。
。基于这些组件，智能体可以被划分为以下五类：
-
基于价值的智能体：拥有隐式的策略，智能体根据价值函数做出决策
-
基于策略的智能体：拥有显式的策略，智能体搜索最优的动作-价值函数
-
演员-评论员智能体：结合了基于价值和基于策略的智能体
-
基于模型的智能体：尝试基于环境构建一个模型
-
无模型智能体：不尝试学习环境，而是尝试学习策略和价值函数
马尔可夫奖励过程
在上一节中，我们介绍了 MDP。在这一节中，我们将正式定义问题陈述，并看到解决该问题的算法。
MDP（马尔可夫决策过程）用于定义强化学习中的环境，几乎所有的强化学习问题都可以使用 MDP 来定义。
为了理解 MDP，我们需要使用马尔可夫奖励过程（MRP）的概念。MRP 是一个随机过程，通过给每个状态添加奖励率，扩展了马尔可夫链。我们还可以定义一个附加变量来跟踪随时间积累的奖励。形式上，MRP 定义为 
，其中 S 是有限状态空间，P 是状态转移概率函数，R 是奖励函数，且 
是折扣率：
其中 
表示期望值。这里的 R[s] 表示状态 s 的期望奖励。
在 MRP 的情况下，我们可以定义从状态 s 开始时的期望回报为：
其中 G[t] 是我们在前一节中定义的累计增益：
现在，为了最大化累计奖励，智能体将尽力从它进入的每个状态中获取期望的奖励总和。为了做到这一点，我们需要找到最优的价值函数。我们将在下一节看到解决这个问题的算法。
贝尔曼方程
使用贝尔曼方程，我们将值函数分解为两个独立的部分，一个表示即时奖励，另一个表示未来奖励。根据我们之前的定义，即时奖励表示为R[t+1]，未来奖励表示为 
，其中：
现在让我们展开G[t]并将G[t+1]代入其中：
现在，由于我们知道：
使用这个恒等式，我们得到：
这为我们提供了 MRP 的贝尔曼方程：
MDP
现在我们对 MRP 有了基本的理解，可以继续讨论 MDP。MDP 是包含决策的 MRP。环境中的所有状态也是马尔可夫的，因此下一个状态只依赖于当前状态。形式上，MDP 可以使用 
来表示，其中S是状态空间，A是动作集，P是状态转移概率函数，R是奖励函数，γ是折扣率。状态转移概率函数P和奖励函数R被正式定义为：
我们也可以正式定义一个策略π为：
由于 MDP 中的状态被认为是马尔可夫的，MDP 策略只依赖于当前状态，这意味着策略是静态的，也就是说， 
。这意味着每当智能体进入相同的状态时，它将根据之前决定的相同策略做出决策。决策函数可以是随机的，这样智能体就不会一直做出相同的决策，从而能够探索环境。
现在，由于我们希望在 MDP 的情况下使用贝尔曼方程，我们将从 MDP 恢复一个 MRP。给定一个 MDP，
和一个策略 
，状态序列S[1], S[2], ...是策略π上的马尔可夫过程(S,P)。状态和奖励序列S1, R1, S2, R2, ...也是由 
给出的 MRP，其中：
我们也可以类似地将奖励函数表示为：
并且，由于我们知道 MDP 的状态值函数Vπ(s)是从状态S开始，然后遵循策略π的期望回报，值函数表示为：
同样，动作值函数可以表示为：
拥有这些值后，我们可以再次推导出 MDP 情况下的贝尔曼期望方程。我们再次通过将状态值函数分解为即时奖励和未来奖励来开始：
与 MRP 的情况类似，动作值函数也可以分解为：
由于我们从每个状态S有多个动作，并且策略定义了动作的概率分布，因此我们需要对其进行平均，以得到贝尔曼期望方程：
我们还可以对所有可能的动作值进行平均，以了解在给定状态S下有多好：
代码示例
在以下的代码示例中，我们实现了一个简单的 MDP：
import numpy as np
import random
class MDP(object):
"""
Defines a Markov Decision Process containing:
- States, s
- Actions, a
- Rewards, r(s,a)
- Transition Matrix, t(s,a,_s)
Includes a set of abstract methods for extended class will
need to implement.
"""
def __init__(self, states=None, actions=None, rewards=None, transitions=None,
discount=.99, tau=.01, epsilon=.01):
"""
Parameters:
-----------
states: 1-D array
The states of the environment
actions: 1-D array
The possible actions by the agent.
rewards: 2-D array
The rewards corresponding to each action at each state of the environment.
transitions: 2-D array
The transition probabilities between the states of the environment.
discount: float
The discount rate for the reward.
"""
self.s = np.array(states)
self.a = np.array(actions)
self.r = np.array(rewards)
self.t = np.array(transitions)
self.discount = discount
self.tau = tau
self.epsilon = epsilon
# Value iteration will update this
self.values = None
self.policy = None
def getTransitionStatesAndProbs(self, state, action):
"""
Returns the list of transition probabilities
"""
return self.t[state][action][:]
def getReward(self, state):
"""
Gets reward for transition from state->action->nextState.
"""
return self.r[state]
def takeAction(self, state, action):
"""
Take an action in an MDP, return the next state
Chooses according to probability distribution of state transitions,
contingent on actions.
"""
return np.random.choice(self.s, p=self.getTransitionStatesAndProbs(state, action))
def valueIteration(self):
"""
Performs value iteration to populate the values of all states in
the MDP.
"""
# Initialize V_0 to zero
self.values = np.zeros(len(self.s))
self.policy = np.zeros([len(self.s), len(self.a)])
policy_switch = 0
# Loop until convergence
while True:
# To be used for convergence check
oldValues = np.copy(self.values)
for i in range(len(self.s)-1):
self.values[i] = self.r[i] + np.max(self.discount * \
np.dot(self.t[i][:][:], self.values))
# Check Convergence
if np.max(np.abs(self.values - oldValues)) <= self.epsilon:
break
def extractPolicy(self):
"""
Extract policy from values after value iteration runs.
"""
self.policy = np.zeros([len(self.s),len(self.a)])
for i in range(len(self.s)-1):
state_policy = np.zeros(len(self.a))
state_policy = self.r[i] + self.discount* \
np.dot(self.t[i][:][:], self.values)
# Softmax the policy
state_policy -= np.max(state_policy)
state_policy = np.exp(state_policy / float(self.tau))
state_policy /= state_policy.sum()
self.policy[i] = state_policy
def simulate(self, state):
"""
Runs the solver for the MDP, conducts value iteration, extracts policy,
then runs simulation of problem.
NOTE: Be sure to run value iteration (solve values for states) and to
extract some policy (fill in policy vector) before running simulation
"""
# Run simulation using policy until terminal condition met
while not self.isTerminal(state):
# Determine which policy to use (non-deterministic)
policy = self.policy[np.where(self.s == state)[0][0]]
p_policy = self.policy[np.where(self.s == state)[0][0]] / \
self.policy[np.where(self.s == state)[0][0]].sum()
# Get the parameters to perform one move
stateIndex = np.where(self.s == state)[0][0]
policyChoice = np.random.choice(policy, p=p_policy)
actionIndex = np.random.choice(np.array(np.where(self.policy[state][:] == policyChoice)).ravel())
# Take an action, move to next state
nextState = self.takeAction(stateIndex, actionIndex)
print "In state: {}, taking action: {}, moving to state: {}".format(
state, self.a[actionIndex], nextState)
# End game if terminal state reached
state = int(nextState)
if self.isTerminal(state):
# print "Terminal state: {} has been reached. Simulation over.".format(state)
return state
使用这个 MDP，我们现在可以编写一个简单的投注游戏代码：
class BettingGame(MDP):
"""
Defines the Betting Game:
Problem: A gambler has the chance to make bets on the outcome of
a fair coin flip. If the coin is heads, the gambler wins as many
dollars back as was staked on that particular flip - otherwise
the money is lost. The game is won if the gambler obtains $100,
and is lost if the gambler runs out of money (has 0$). This gambler
did some research on MDPs and has decided to enlist them to assist
in determination of how much money should be bet on each turn. Your
task is to build that MDP!
Params:
pHead: Probability of coin flip landing on heads
- Use .5 for fair coin, otherwise choose a bias [0,1]
"""
def __init__(self, pHeads=.5, discount=.99, epsilon=.1, tau=.0001):
MDP.__init__(self,discount=discount,tau=tau,epsilon=epsilon)
self.pHeads = pHeads
self.setBettingGame(pHeads)
self.valueIteration()
self.extractPolicy()
# Edge case fix: Policy for $1
self.policy[1][:] = 0
self.policy[1][1] = 1.0
def isTerminal(self, state):
"""
Checks if MDP is in terminal state.
"""
return True if state is 100 or state is 0 else False
def setBettingGame(self, pHeads=.5):
"""
Initializes the MDP to the starting conditions for
the betting game.
Params:
pHeads = Probability that coin lands on head
- .5 for fair coin, otherwise choose bias
"""
# This is how much we're starting with
self.pHeads = pHeads
# Initialize all possible states
self.s = np.arange(102)
# Initialize possible actions
self.a = np.arange(101)
# Initialize rewards
self.r = np.zeros(101)
self.r[0] = -5
self.r[100] = 10
# Initialize transition matrix
temp = np.zeros([len(self.s),len(self.a),len(self.s)])
# List comprehension using tHelper to determine probabilities for each index
self.t = [self.tHelper(i[0], i[1], i[2], self.pHeads) for i,x in np.ndenumerate(temp)]
self.t = np.reshape(self.t, np.shape(temp))
for x in range(len(self.a)):
# Remember to add -1 to value it, and policy extract
# Send the end game states to the death state!
self.t[100][x] = np.zeros(len(self.s))
self.t[100][x][101] = 1.0
self.t[0][x] = np.zeros(len(self.s))
self.t[0][x][101] = 1.0
def tHelper(self, x, y, z, pHeads):
"""
Helper function to be used in a list comprehension to quickly
generate the transition matrix. Encodes the necessary conditions
to compute the necessary probabilities.
Params:
x,y,z indices
pHeads = probability coin lands on heads
"""
# If you bet no money, you will always have original amount
if x + y is z and y is 0:
return 1.0
# If you bet more money than you have, no chance of any outcome
elif y > x and x is not z:
return 0
# If you bet more money than you have, returns same state with 1.0 prob.
elif y > x and x is z:
return 1.0
# Chance you lose
elif x - y is z:
return 1.0 - pHeads
# Chance you win
elif x + y is z:
return pHeads
# Edge Case: Chance you win, and winnings go over 100
elif x + y > z and z is 100:
return pHeads
else:
return 0
return 0
总结
在本章中，我们首先简要介绍了强化学习。我们讨论了智能体、奖励和强化学习中的学习目标。在下一节中，我们介绍了 MRP，它是 MDP 的一个主要概念。在理解了 MRP 之后，我们接下来介绍了 MDP 的概念，并附上了一个代码示例。
浙公网安备 33010602011771号