使用-Python-构建概率图模型-全-

使用 Python 构建概率图模型（全）

原文：Building probabilistic graphical models with Python

协议：CC BY-NC-SA 4.0

零、前言

在这本书里，我们从图模型的基础，它们的类型，为什么使用它们，以及它们解决什么样的问题的探索之旅开始。然后，我们在图模型的上下文中探索子问题，例如它们的表示、构建它们、学习它们的结构和参数，并使用它们来回答我们的推理查询。

这本书试图给出足够的理论信息，然后使用代码样本来了解一些算法是如何实现的。代码示例还提供了一个方便的模板来构建图模型和回答我们的概率查询。在文献中描述的多种图模型中，这本书主要关注离散贝叶斯网络，偶尔有来自马尔可夫网络的例子。

这本书涵盖了什么

第 1 章、概率，涵盖了理解图模型所需的概率概念。

第 2 章、定向图模型提供了关于贝叶斯网络的信息，以及它们与独立性、条件独立性和 D-分离相关的属性。本章使用代码片段加载贝叶斯网络并了解其独立性属性。

第三章、无向图模型，涵盖了马尔可夫网络的性质，它们与贝叶斯网络的区别，以及它们的独立性。

第 4 章、结构学习，涵盖了使用数据集推断贝叶斯网络结构的多种方法。我们还学习了结构学习的计算复杂性，并在本章中使用代码片段来学习采样数据集中给出的结构。

第 5 章、参数学习，介绍了使用来自 PyMC 的代码样本进行参数学习的最大似然和贝叶斯方法。

第 6 章、使用图模型的精确推理，解释了用于精确推理的变量消除算法，并探索了使用相同算法回答我们的推理查询的代码片段。

第 7 章、近似推理方法探讨了对于太大而无法进行精确推理的网络的近似推理。我们还将浏览在马尔可夫网络上使用循环信念传播运行近似推理的代码样本。

附录、参考文献，包括所有有助于轻松理解书中章节的链接和网址。

这本书你需要什么

要运行本书中的代码示例，您需要一台安装了 IPython 的笔记本电脑或台式机。我们在本书中使用了几个软件包，大部分都可以使用pip或easy_install等 Python 安装程序进行安装。在某些情况下，软件需要从源代码编译，可能需要 C++编译器。

这本书是给谁的

本书面向熟悉 Python 并希望使用代码示例探索图模型细微差别的开发人员。

这本书对于已经从理论上介绍了图模型，并且希望实现图模型的实现，并对不同(图模型)库处理现实世界模型的能力有所了解的学生来说也是理想的。

熟悉分类和回归模型的机器学习实践者，以及希望探索和实验图模型可以解决的问题类型的人，也会发现这本书是一个无价的资源。

这本书将图模型视为一种工具，可以用来解决机器学习领域的问题。此外，它并没有试图解释图模型的数学基础，也没有详细说明所使用的每种算法的步骤。

惯例

在这本书里，你会发现许多区分不同种类信息的文本风格。以下是这些风格的一些例子，以及对它们的含义的解释。

文本中的码字、数据库表名、文件夹名、文件名、文件扩展名、路径名、伪 URL、用户输入和 Twitter 句柄如下所示:“我们可以通过创建一个TfidfVectorizer对象来完成同样的工作。”

代码块设置如下:

clf = MultinomialNB(alpha=.01)
print "CrossValidation Score: ", np.mean(cross_validation.cross_val_score(clf,vectors, newsgroups.target, scoring='f1'))
CrossValidation Score:  0.954618416381

注

警告或重要提示会出现在这样的框中。

型式

提示和技巧是这样出现的。

读者反馈

我们随时欢迎读者的反馈。让我们知道你对这本书的看法——你喜欢或可能不喜欢什么。读者反馈对我们开发您真正能从中获得最大收益的标题非常重要。

要给我们发送一般反馈，只需向<[feedback@packtpub.com](mailto:feedback@packtpub.com)>发送电子邮件，并通过您的消息主题提及书名。

如果你对某个主题有专业知识，并且对写作或投稿感兴趣，请参阅我们在www.packtpub.com/authors上的作者指南。

客户支持

现在，您已经自豪地拥有了一本书，我们有许多东西可以帮助您从购买中获得最大收益。

下载示例代码

您可以从您在http://www.packtpub.com的账户中下载您购买的所有 Packt 书籍的示例代码文件。如果您在其他地方购买了这本书，您可以访问http://www.packtpub.com/support并注册，以便将文件直接通过电子邮件发送给您。

勘误表

尽管我们尽了最大努力来确保我们内容的准确性，但错误还是会发生。如果你在我们的某本书里发现了错误——可能是文本或代码中的错误——如果你能向我们报告，我们将不胜感激。通过这样做，你可以让其他读者免受挫折，并帮助我们改进这本书的后续版本。如果您发现任何勘误表，请访问http://www.packtpub.com/submit-errata，选择您的书籍，点击勘误表提交表链接，并输入您的勘误表的详细信息。一旦您的勘误表得到验证，您的提交将被接受，勘误表将上传到我们的网站上，或添加到该标题的勘误表部分下的任何现有勘误表列表中。通过从http://www.packtpub.com/support中选择您的标题，可以查看任何现有的勘误表。

盗版

互联网上版权材料的盗版是所有媒体的一个持续问题。在 Packt，我们非常重视版权和许可证的保护。如果您在互联网上遇到任何形式的我们作品的非法拷贝，请立即向我们提供位置地址或网站名称，以便我们寻求补救。

请通过<[copyright@packtpub.com](mailto:copyright@packtpub.com)>联系我们，获取疑似盗版资料的链接。

我们感谢您在保护我们作者方面的帮助，以及我们为您带来有价值内容的能力。

问题

如果您对本书的任何方面有问题，可以在<[questions@packtpub.com](mailto:questions@packtpub.com)>联系我们，我们将尽最大努力解决。

一、概率

在我们踏上穿越图模型的旅程之前，我们必须为自己配备一些有助于理解的工具。我们将首先介绍概率及其概念，如随机变量和分布类型。

然后，我们将试图理解概率可以帮助我们回答的问题类型以及概率的多种解释。最后，我们将快速了解一下贝叶斯规则，它有助于我们理解概率之间的关系，还将了解附带的条件概率和链式规则的概念。

概率论

我们经常会遇到这样的情况:我们不得不锻炼自己对事件发生的主观信念；例如，天气或交通等本质上随机的事件。概率也可以理解为主观信念的程度。

当我们谈论天气(例如，今天晚上)时，可以理解为天气可以有多种结果，如下雨、晴天或多云。所有可能结果的空间被称为事件(也称为样本空间)。例如，掷骰子的结果是一组从 1 到 6 的数字。在处理可测量的结果时，例如掷骰子或今天的天气(可能是雨天、晴天或阴天)，我们可以为每个结果分配一个概率值，以概括我们对这些结果的信任程度。用来表达我们信念的记数法的一个例子是 P(雨天)=0.3 ，可以解读为下雨的概率是 0.3%或者 30%。

Kolmogorov 提出的概率公理陈述如下:

• The probability of an event is a non-negative real number (that is, the probability that it will rain today may be small, but nevertheless will be greater than or equal to 0). This is explained in mathematical terms as follows:

  

• The probability of the occurrence of some event in the sample space is 1 (that is, if the weather events in our sample space are rainy, sunny, and cloudy, then one of these events has to occur), as shown in the following formula:

  

• The sum of the probabilities of mutually exclusive events gives their union, as given in the following formula:

  

当我们讨论骰子或抛硬币的公平(或不公平)时，我们讨论的是概率的另一个关键方面，即模型参数。公平硬币的概念转化为这样一个事实，即控制参数的值为 0.5，有利于头部，这也转化为这样一个事实，即我们假设所有结果的可能性相等。在这本书的后面，我们将研究需要多少参数来完全指定一个概率分布。然而，我们正在超越自己。首先让我们了解一下概率分布。

一个概率分布由与每个可测量结果相关的概率组成。在离散结果(如掷骰子或掷硬币)的情况下，分布由概率质量函数指定，而在连续结果(如一个班级中学生的身高)的情况下，分布由概率密度函数指定。

让我们用一个的例子来看看离散分布。抛硬币有两种结果:正面和反面，一枚公平的硬币会给所有结果分配相同的概率。这意味着概率分布很简单——头部为 0.5，尾部为 0.5。像这样的分布(例如，正面 0.3，反面 0.7)将是对应于有偏硬币的分布。下图显示了掷出两个骰子时数值总和的离散概率分布:

给所有结果分配相等概率的分布称为均匀分布。这是我们将探索的许多发行版之一。

让我们来看看与连续结果相关的常见分布之一，即高斯或正态分布，它呈钟形，因此被称为钟形曲线(尽管还有其他形状类似于钟形的分布)。以下是一些真实世界的例子:

• 一个班级学生的身高是对数正态分布的(如果我们取学生身高的对数并作图，得到的分布是正态分布的)
• 物理实验中的测量误差

高斯分布有两个参数:平均值()和方差()。参数均值和方差分别决定了中间点和远离均值的分布离散度。

下图显示了均值和方差值不同的多个高斯分布。可以看出，方差越大，分布越广，而平均值会移动 x 轴上的峰值，如下图所示:

概率推理的目标

既然我们已经理解了概率的概念，我们必须问问自己这是如何使用的。我们问的这类问题分为以下几类:

• 第一个问题是参数估计，比如，一个硬币有偏差还是公平？如果有偏差，参数值是多少？

• The second question is that given the parameters, what is the probability of the data? For example, what is the probability of five heads in a row if we flip a coin where the bias (or parameter) is known.

  上述问题取决于数据(或缺乏数据)。如果我们有一组硬币翻转的观测值，我们可以估计控制参数(即参数估计)。如果我们有一个参数的估计，我们想估计硬币翻转产生的数据的概率(第二个问题)。然后，有时我们会来回改进模型。

• 模型是否适合这个问题，是我们可以询问的第三个问题。有没有一个单一的参数控制投币实验的结果？当我们希望对一个复杂的现象(如交通或天气预测)建模时，模型中肯定存在几个参数，其中数百甚至数千个参数并不罕见。在这种情况下，我们想问的问题是，哪种模型更适合数据？我们将在后面的章节中看到一些关于模型拟合不同方面的例子。

条件概率

让我们用一个具体的例子，我们有一群申请工作的候选人。一个事件( x )可能是所有获得录用的候选人的集合，而另一个事件( y )可能是所有经验丰富的候选人的集合。我们可能想推理一下联合事件的集合()，这是一组有经验的候选人得到了一个报价(联合事件的概率也写成)。提出的问题是，如果我们知道一个事件已经发生，是否会改变另一个事件发生的概率。在这种情况下，如果我们确定某个候选人得到了一份工作，这能告诉我们他们的经历吗？

条件概率正式定义为，可以理解为 x 的概率，假设 y 发生。分母是联合分布所有可能结果的总和，取值 x 相加，即。

链式法则

链式法则允许我们利用一组随机变量的条件概率来计算它们的联合分布。换句话说，联合分布是单个条件概率的乘积。从开始，如果是事件，。

我们将在图模型中详细回到这一点，其中链式规则通过将一个大问题(计算联合分布)拆分成更小的问题(条件概率)来帮助我们分解它。

贝叶斯法则

贝叶斯法则是概率论的基础之一，这里就不多赘述了。它来自条件概率的定义，如下式所示:

从公式中，我们可以推断出关于贝叶斯规则的以下内容——我们对我们正在推理的问题抱有先验信念。这就是所谓的前一项。当我们开始看到数据时，我们的信念发生变化，这就产生了我们最终的信念(称为后验)，如下式所示:

让我们用一个例子来看看贝叶斯规则背后的直觉。艾米和卡尔正站在火车站等火车。在过去的一年里，艾米每天都乘同一列火车，这是卡尔第一天来车站。他们之前对火车准点的看法是什么？

在过去的一年里，艾米每天都在赶火车，她总是看到火车在预定的出发时间两分钟内到达。因此，她坚信火车最多会晚点两分钟。由于这是卡尔的第一天，他不知道火车会准时。然而，卡尔在过去的一年里一直在环游世界，去过火车不准时的地方。因此，他不太相信火车会晚点 30 分钟。

第一天，火车晚点了 5 分钟。这个观察对艾米和卡尔的影响是不同的。由于艾米有很强的优先性，她的信念稍微修改了一下，接受了火车可能晚点 5 分钟的事实。卡尔的信念现在改变了，认为这里的火车相当准时。

换句话说，后验信念受到多种方式的影响:当具有强先验的人看到一些观察结果时，他们的后验信念与前验信念相比变化不大。另一方面，当先验较弱的人看到大量的观测值(强似然)时，他们的后验信念会发生很大的变化，并且很大程度上受观测值(似然)而不是前验信念的影响。

让我们看一个贝叶斯规则的数值例子。 D 是运动员使用兴奋剂 ( PEDs ) 的赛事。 T 是药检返回阳性的事件。在整个讨论中，我们使用质数( ' )符号来表示事件没有发生；例如， D' 代表运动员没有使用 PEDs 的项目。

P(D|T) 是在药检返回阳性的情况下，运动员使用 PEDs 的概率。 P(T|D) 是在运动员使用 PEDs 的情况下，药检返回阳性的概率。

进行药物测试的实验室声称，它可以在 90%的时间内检测到 PEDs。我们还了解到假阳性率(测试呈阳性但未使用 PEDs 的运动员)为 15%，10%的运动员使用 PEDs。如果药检结果呈阳性，运动员使用 PEDs 的概率有多大？

根据贝叶斯规则的基本形式，我们可以写出下面的公式:

现在，我们有以下数据:

• P(T|D) :等于 0.90

• P(T|D') :这等于 0.15(假设运动员没有使用 PEDs，测试返回阳性)

• P(D): This is equal to 0.1

  当我们代入这些值时，我们得到的最终值为 0.4，如下式所示:

  

这个结果似乎有点违反直觉，因为尽管 PEDs 检测呈阳性，运动员使用 PEDs 的几率只有 40%。这是因为 PEDs 本身的使用率很低(只有 10%的运动员使用 PEDs)，而且假阳性率相对较高(0.15%)。

概率的解释

在前面的例子中，我们注意到了我们是如何拥有先验信念的，并且观测数据的引入可以改变我们的信念。然而，这个观点是对概率的多种解释之一。

第一种(我们已经讨论过了)是贝叶斯解释，认为概率是一种信念程度，信念程度在考虑证据前后变化。

第二种观点被称为频率论解释，其中概率衡量结果的比例，并假设先前的信念是一个没有数据支持的不正确的概念。

为了用一个例子来说明这一点，让我们回到抛硬币实验，我们希望了解硬币的偏差。我们做了两个实验，分别翻转硬币 10 次和 10000 次。在第一个实验中，我们得到 7 个头，在第二个实验中，我们得到 7000 个头。

从频率主义者的观点来看，在这两个实验中，获得头部的概率是 0.7 (7/10 或 7000/10000)。然而，我们可以很容易地说服自己，我们对第二个实验的结果比第一个实验更有信心。这是因为第一个实验的结果有贝叶斯的观点，如果我们有先验的信念，第二个实验的观察会压倒先验，这在第一个实验中是不太可能的。

对于以下部分的讨论，让我们考虑一个公司面试求职者的例子。在邀请候选人参加面试之前，会根据候选人的经验以及候选人在毕业结果中获得的平均绩点对候选人进行筛选。如果候选人通过筛选，他将被要求参加面试。一旦候选人被面试，公司可能会向候选人发出工作邀请(这是基于候选人在面试中的表现)。候选人也在评估攻读研究生学位，候选人是否被录取他选择的研究生学位课程取决于他在学士学位中的成绩。下图直观地展示了我们对影响工作选择(和研究生学位录取)标准的因素之间的关系的理解:

随机变量

随机变量的经典概念是其值会因偶然性而变化(维基百科)。大多数程序员都遇到过编程语言标准库中的随机数。从程序员的角度来看，与普通变量不同，随机变量每次读取其值时都会返回一个新值，其中变量的值可能是新调用随机数生成器的结果。

我们早些时候已经看到了事件的概念，以及我们如何从一组可测量的事件中考虑单个事件发生的概率。然而，考虑结果的属性可能是合适的。

在求职者的求职例子中，求职者的属性之一是他的经验。该属性可以采用多个值，例如高度相关或不相关。讨论不同结果中的属性及其值的形式机制被称为随机变量[柯勒等人，2.1.3.1]。

随机变量可以采用分类值(如硬币投掷结果的 {Heads，Tails} )或真实值(如一个班级学生的身高)。

边际分布

我们已经看到求职示例(在上图中描述)有五个随机变量。分别是成绩经历面试offer录取。这些随机变量有一组对应的事件。

现在，让我们考虑随机变量的子集 X ，其中 X 只包含经验随机变量。该子集包含高度相关的事件和不相关的事件。**

**如果我们在子集 X 中登记所有事件的概率，这将被称为边际分布，其示例可以在下面的公式中找到:

像所有有效分布一样，概率总和应该是 1。

随机变量集(描述边际分布的)可以只包含一个变量(如前面的例子)，也可以包含几个变量，如{经验、等级、面试} 。

联合配送

我们已经看到边际分布是描述随机变量子集的分布。接下来，我们将讨论描述集合中所有随机变量的分布。这叫做联合分配。让我们看看求职示例中涉及学位分数和经验随机变量的联合分布:

| **度评分** | **相关经验** |   | |   | 高度相关 | 不相关 |   | | 贫穷的；贫困的 | Zero point one | Zero point one | Zero point two | | 平均的 | Zero point one | Zero point four | Zero point five | | 优秀的 | Zero point two | Zero point one | Zero point three | |   | Zero point four | Zero point six | one |

深灰色单元格中的值是联合分布，浅灰色单元格中的值是边缘分布(有时称为边缘分布，因为它们写在页边距上)。可以观察到，单个边际分布的总和为 1，就像正态概率分布一样。

一旦描述了联合分布，就可以通过对单个行或列求和来找到边际分布。在上表中，如果我们将各列相加，第一列给出了高度相关的概率，第二列给出了不相关的概率。可以看出，对行应用类似的策略给了我们学位得分的概率。

独立

独立事件的概念可以通过看一个例子来理解。假设我们有两个骰子，当我们一起掷骰子时，一个骰子的分数为 2，另一个骰子的分数为 3。不难看出，这两个事件是相互独立的，因为每个骰子的结果不会影响另一个或与之相互作用。

我们可以用多种方式定义独立的概念。假设我们有两个事件 a 和 b ，这两个事件结合的概率只是它们概率的乘积，如下式所示:

如果我们把的概率写成(即在事件 a 已经发生的情况下， a 的概率和 b 的概率的乘积)，如果事件 a 和 b 是独立的，它们就分解为。

指定独立性的另一种方式是说 a 给 b 的概率只是 a 的概率，也就是说 b 的出现不影响 a 的概率，如下式所示:

可见独立性是对称的，即。虽然这种独立性的定义是在事件的背景下，但同样的概念可以推广到随机变量的独立性。

条件独立

给定两个事件，确定它们是否独立并不总是明显的。考虑一个求职者，他申请了两家公司的工作，脸书和谷歌。它可能导致两个事件，第一个是候选人接到谷歌的面试电话，另一个事件是他接到脸书的面试电话。知道第一个事件的结果是否能告诉我们第二个事件的概率？是的，确实如此，因为我们可以推断，如果一个候选人足够聪明，能够接到谷歌的电话，他就是一个有前途的候选人，而且接到脸书电话的概率相当高。

到目前为止已经确定的是，这两个事件不是独立的。假设我们得知公司根据候选人的成绩决定发送面试邀请，我们得知候选人的成绩是 A，由此我们推断候选人相当聪明。我们可以推断，由于候选人相当聪明，知道候选人接到谷歌的面试电话不会告诉我们更多关于他感知的智力的信息，也不会改变来自脸书的面试电话的概率。这可以被正式注释为给定谷歌面试呼叫和 A 级的脸书面试呼叫的概率等于给定 A 级的脸书面试呼叫的概率，如下式所示:

换句话说，我们是说来自脸书的邀请是有条件独立于来自谷歌的邀请的，给定的候选人具有 a 级

查询类型

了解了联合概率分布和条件概率分布之后，让我们将注意力转向我们可以对这些分布提出的查询类型。

概率查询

这是最常见的查询类型，由以下两部分组成:

• 证据:这是已经观察到的随机变量的子集 E
• 查询:这个是随机变量的子集 Y

我们希望计算概率的值，它是后验概率或超过 Y 的边际概率。再次使用求职者的例子，我们可以计算面试电话的边际分布，条件是学位分数= A 级。

地图查询

最大后验概率(MAP )是对变量的某些子集的最高概率联合赋值。在概率查询的情况下，重要的是概率的值。在地图的情况下，计算联合分配的精确概率值是次要的，相比于寻找所有随机变量的联合分配的任务。

如果概率值相等，可以返回多个联合赋值。我们将从一个例子中看到，在联合赋值的情况下，每个边际值的最高概率可能不是最高的联合赋值(下面的例子来自柯勒等人)。

考虑两个非独立随机变量 X 和 Y ，其中 Y 依赖于 X 。下表显示了 X 上的概率分布:

|

X0

|

X1

|
| --- | --- |
| Zero point four | Zero point six |

我们可以看到随机变量 X 的 MAP 赋值是 X1 ，因为它有更高的值。下表显示了 X 和 Y 的边际分布:

|

P(Y|X)

|

Y0

|

Y1

|
| --- | --- | --- |
| X0 | Zero point one | Zero point nine |
| X1 | Zero point five | Zero point five |

X 和 Y 上的联合分布如下表所示:

|

分配

|

价值

|
| --- | --- |
| X0，Y0 | Zero point zero four |
| X0，Y1 | Zero point three six |
| X1，Y0 | Zero point three |
| X1，Y1 | Zero point three |

在上表所示的联合分布中，随机变量(X，Y)的 MAP 赋值为(X0，Y1)，而 X (X1)的 MAP 赋值不是联合赋值的 MAP 的一部分。综上所述，MAP 赋值不能简单地取每个随机变量在边际分布中的最大概率值。

另一种不同类型的 MAP 查询是边际 MAP 查询，在边际 MAP 查询中，我们只有构成查询的变量子集，而不是联合分布。在前面的示例中，边际 MAP 查询将是 MAP (Y)，这是 MAP 分配给随机变量 Y 的最大值，可以通过查看联合分布并对 x 的值求和来读取。从下表中，我们可以读取最大值并确定 MAP (Y)是 Y1:

|

分配

|

价值

|
| --- | --- |
| Y0 | Zero point three four |
| Y1 | Zero point six six |

注

边际查询的数据由萨古尔·斯里哈里通过查询联合概率分布获得。你可以在找到它。

总结

在这一章中，我们研究了基本概率、随机变量和贝叶斯定理的概念。我们还通过一个求职者的例子学习了链式法则、联合分布和边际分布，我们将在后面的章节中回到这个例子。掌握了这些主题后，我们现在可以在接下来的章节中继续探索贝叶斯和马尔可夫网络，我们将正式描述这些网络来回答我们在本章中讨论的一些概率查询。虽然这一章完全是理论性的，但从下一章开始，我们将实现 Python 代码来寻找问题的答案。**

二、有向图模型

在这一章中，我们将学习定向图模型，也称为贝叶斯网络。我们从什么(我们试图解决的问题)、如何(图表示)、为什么(因式分解以及 CPD 和图因式分解的等价性)开始，然后继续使用 Libpgm Python 库来玩一个小的贝叶斯网。

图术语

在我们进入贝叶斯网络之前，让我们学习一些图术语。图 G 由一组节点(也称为顶点)和另一组边组成。连接一对节点的边可以有两种类型:有向的(由表示)和无向的(由表示)。图也可以表示为邻接矩阵，在无向图的情况下，如果位置 G(i，j) 包含 1 ，则表示 i 和 j 顶点之间的边。在有向图的情况下， 1 或 -1 的值表示边的方向。

在许多情况下，我们对所有边都是有向或无向的图感兴趣，导致它们分别被称为有向图或无向图。

有向图中节点的父节点是具有终止于的输出边的一组节点。

节点的子节点是具有离开的输入边的节点集

节点的度是它参与的边的数量。

团是一组节点，其中每对节点由一条边连接。最大团是指如果包含任何其他节点，则失去团属性的团。

如果存在从一个节点经过其他节点后返回自身的路径，则称之为循环或循环。

一个有向无环图 ( DAG ) 是一个没有循环的图。

一个部分有向无环图 ( PDAG ) 是一个既可以包含有向边又可以包含无向边的图。

森林是一组树。

Python 题外话

我们将很快开始探索使用 Python 的 GMs，这是回顾您的 Python 安装的好时机。本书示例的推荐基础 Python 安装是 IPython，所有平台都可以使用。有关特定平台的文档，请参考 IPython 网站。

我们还将使用多个 Python 库来探索图模型的各个领域。除非另有说明，安装 Python 库的常用方法是使用pip install <packagename>或easy_install <packagename>。

要使用本章中的代码，请安装 Lipgm(https://pypi.python.org/pypi/libpgm)和 scipy(http://scipy.org/)。

独立性和独立参数

图模型解决的关键问题之一是定义联合分布。我们来看一下求职面试的例子，有一定经验和教育程度的求职者在找工作。这位候选人也在申请进入高等教育项目。

我们正试图充分说明工作机会上的联合分配，这(根据我们的直觉)取决于工作面试的结果、候选人的经验和他的成绩(我们假设候选人被研究生院录取与工作机会无关)。三个随机变量 {Offer，Experience，gradies }取两个值(比如对于工作 Offer 取 yes 和 no，对于工作经历取高度相关和不相关)，面试取三个值，联合分布会用一个 24 行的表来表示(即 2×2×2×3)。

每行包含该行中随机变量赋值的概率。虽然表的不同实例可能有不同的概率分配，但我们需要 24 个参数(每行一个)来编码表中的信息。然而，为了计算的目的，我们将只需要 23 个独立的参数。为什么我们要移除一个？由于概率之和等于 1，最后一个参数可以通过从已经找到的 23 个参数之和中减去 1 来计算。

|   |

经验

|

等级

|

采访

|

提供

|

可能性

|
| --- | --- | --- | --- | --- | --- |
| Zero | Zero | Zero | Zero | Zero | 0.7200 |
| one | Zero | Zero | Zero | one | 0.0800 |
| Two | Zero | Zero | one | Zero | 0.0720 |
| three | Zero | Zero | one | one | 0.1080 |
| 。 |   |   |   |   |   |
| 。 |   |   |   |   |   |
| 行省略 |
| Twenty-one | one | one | one | one | 0.1200 |
| Twenty-two | one | one | Two | Zero | 0.0070 |
| Twenty-three | one | one | Two | one | 0.6930 |

前面的联合分布是printjointdistribution.ipynb IPython 笔记本的输出，它打印随机变量的经验、成绩、面试和报价的所有排列，以及它们的概率。

观察上表应该清楚，由于以下原因，很难获得完全指定的联合分布:

• 从计算的角度来看，它太大了，无法存储和操作
• 我们将需要联合分布的每个赋值的大量数据来正确地得出概率
• 大的联合分布中的个体概率变得非常小，对人类的理解不再有意义

如何才能避免必须指定联合分配？我们可以通过使用独立参数的概念来做到这一点，我们将在下面的例子中探讨这一点。

等级和录取上的联合分布 P 如下(在本书中，上标 0 和 1 表示低分和高分):

|

等级

|

准许进入

|

概率

|
| --- | --- | --- |
| S ⁰ | A ⁰ | Zero point six six five |
| S ⁰ | A ¹ | Zero point zero three five |
| S ¹ | A ⁰ | Zero point zero six |
| S ¹ | A ¹ | Zero point two four |

当我们从因果关系的角度对研究生录取进行推理时，很明显录取取决于成绩，也可以用条件概率表示如下:

上式要求的参数个数为三个，一个，和各两个。由于这是一个简单的分布，所以条件分布和联合分布所需的参数数量是相同的，但是让我们观察整个网络，看看条件参数化是否有所不同:

我们如何计算前图中贝叶斯网的参数数量？让我们一次一个参数地浏览每个条件概率表。体验和等级取两个值，因此每个需要一个独立的参数。面试表有 12 个(3×4)参数。然而，每行加起来是 1，因此，每行需要两个独立的参数。整个表需要 8 个(2×4)独立参数。类似地， Offer 表有 6 个条目，但每行只需要 1 个独立参数，这样就有 3 个(1×3)独立参数。因此，参数总数为 1(经验)+ 1(等级)+ 12(面试)+ 3(报价)总计为 17，这比 24 个参数要少得多，以充分指定联合分布。因此，贝叶斯网络中的独立性假设有助于我们避免指定联合分布。

贝叶斯网络

贝叶斯网络是一种可以表示为有向无环图的结构，它所包含的数据可以从以下两个角度来看:

• 它允许使用贝叶斯网络的链式规则对联合分布进行紧凑和模块化的表示
• 它允许观察顶点之间的条件独立性假设

我们将探讨到目前为止已经看到的求职面试例子中的两个想法(顺便说一下，这是一个贝叶斯网络)。

贝叶斯网络的模块化结构是一组局部概率模型，表示每个变量对其父变量的依赖性质(柯勒等人 3.2.1.1)。经验和等级各有一个概率分布，条件概率分布 ( CPD )分别存在于面试和报价中。给定分配给其父代的所有组合，CPD 指定随机变量的分布。因此，给定贝叶斯网络的模表示是每个随机变量的 CPD 集合。

条件独立的观点来自不同随机变量之间的边缘(我们的直觉得出)，在那里，我们假设工作面试的要求必须取决于候选人的经验以及他在学位课程中获得的分数，而获得工作机会的概率仅取决于工作面试的结果。

链式法则

链式法则允许我们将联合分布定义为因素的乘积。在求职面试的例子中，利用概率的链式法则，我们可以写出以下公式:

在上式中， E 、 G 、 I 、 O 分别代表体验、等级、面试、 Offer 。然而，我们可以使用由图编码的条件独立性假设来重写联合分布，如下所示:

这是贝叶斯网络的链式规则的一个例子。更一般地说，我们可以这样写:

这里，是图 G 中的节点，是图 G 中节点的父节点

推理模式

在本节中，我们将研究贝叶斯网络中使用的不同类型的推理。我们将使用 Libpgm 库创建一个贝叶斯网络。Libpgm 从具有特定格式的 JSON 格式的文件中读取网络信息，如节点、边和与每个节点相关联的 CPD 概率。这个 JSON 文件被读入NodeData和GraphSkeleton对象，以创建一个离散贝叶斯网络(顾名思义，这是一个贝叶斯网络，其中 CPD 取离散值)。TableCPDFactorization对象是包装离散贝叶斯网络的对象，允许我们查询网络中的 CPD。本例中的 JSON 文件job_interview.txt应该与 IPython 笔记本放在同一个文件夹中，以便自动加载。

下面的讨论使用整数 0、1 和 2 作为每个随机变量的离散结果，其中 0 是最差的结果。例如面试= 0 表示面试结果最差，面试= 2 表示结果最好。

因果推理

我们将要探索的第一种推理叫做因果推理。最初，我们观察一个事件的先验概率不受任何证据的制约(对于这个例子，我们将集中于提供的随机变量)。然后我们引入对父变量之一的观察。与我们的逻辑推理一致，我们注意到，如果观察到一个事件的父母之一(相当于原因)，那么我们对孩子的随机变量有更强的信念( Offer )。

我们首先定义一个函数来读取 JSON 数据文件，并创建一个可以用来运行概率查询的对象。以下代码来自Bayes net-Causal Reasoning.ipynb IPython 笔记本:

from libpgm.graphskeleton import GraphSkeleton
from libpgm.nodedata import NodeData
from libpgm.discretebayesiannetwork import DiscreteBayesianNetwork
from libpgm.tablecpdfactorization import TableCPDFactorization

def getTableCPD():
    nd = NodeData()
    skel = GraphSkeleton()
    jsonpath="job_interview.txt"
    nd.load(jsonpath)
    skel.load(jsonpath)
    # load bayesian network
    bn = DiscreteBayesianNetwork(skel, nd)
    tablecpd=TableCPDFactorization(bn)
    return tablecpd

我们现在可以使用specificquery函数在我们定义的网络上运行推理查询。获得T4 优惠的先验概率是多少？请注意，概率查询采用两个字典参数:第一个是查询，第二个是证据集，由空字典指定，如以下代码所示:

tcpd=getTableCPD()
tcpd.specificquery(dict(Offer='1'),dict())

以下是前面代码的输出:

0.432816

大约是 43%，如果我们现在引入候选人成绩差的证据，如何改变获得录用的概率？我们将评估的值，如下代码所示:

tcpd=getTableCPD()
tcpd.specificquery(dict(Offer='1'),dict(Grades='0'))

以下是前面代码的输出:

0.35148

不出所料，这降低了被录取的概率，因为我们推断成绩差的学生不太可能被录取。添加候选人经验也低的进一步证据，我们评估，如下代码所示:

tcpd=getTableCPD()
tcpd.specificquery(dict(Offer='1'),dict(Grades='0',Experience='0'))

以下是前面代码的输出:

0.2078

正如预期的那样，在额外的证据上，它下降得更低，从 35%下降到 20%。

我们所看到的是，观察到的父随机变量的引入强化了我们的信念，这就导致了因果推理这个名字。

在下图中，我们可以看到因果推理和证据推理所采取的不同路径:

证据推理

证据推理是当我们观察一个孩子变量的值时，我们希望推理它如何加强我们对它的父母的信念。我们将评估高Experience 的先验概率，如下代码所示:

tcpd=getTableCPD()
tcpd.specificquery(dict(Experience='1'),dict())

前面代码的输出如下:

0.4

我们现在引入候选人面试良好的证据，并评估P(Experience = 1 | Interview = 2)的值，如下代码所示:

tcpd=getTableCPD()
print tcpd.specificquery(dict(Experience='1'),dict(Interview='2'))

前面代码的输出如下:

0.864197530864

我们看到，如果候选人在面试中得分较高，那么候选人经验丰富的概率就会增加，这遵循的推理是候选人必须有良好的经验或学历，或者两者都有。在证据推理中，我们从结果推理到原因。

因果间推理

顾名思义，因果间推理是一种推理类型，其中单个效应的多个原因相互作用。我们首先确定具有高的相关经验的先验概率；因此，我们将评估 P(Experience=1) ，如下代码所示:

tcpd=getTableCPD()
tcpd.specificquery(dict(Experience='1'),dict())

以下是前面代码的输出:

0.4

通过介绍面试进行得非常顺利的证据，我们认为候选人一定很有经验。我们现在将评估的值，如以下代码所示:

tcpd=getTableCPD()
tcpd.specificquery(dict(Experience='1'),dict(Interview='2'))

以下是前面代码的输出:

0.864197530864

贝叶斯网络证实了我们认为的是真的(候选人是有经验的)，高经验的概率从 0.4 上升到 0.86。现在，如果我们引入证据，证明候选人没有取得好成绩，但仍然设法在面试中获得好成绩，我们可能会得出结论，候选人一定非常有经验，他的成绩根本不重要。我们将评估的值，如下代码所示:

tcpd=getTableCPD()
tcpd.specificquery(dict(Experience='1'),dict(Interview='2',Grades='0'))

前面代码的输出如下:

0.909090909091

这证实了我们的预感，即使高经验的概率只上升了一点点，但它加强了我们对候选人高经验的信念。这个例子展示了求职面试节点的两个家长之间的相互作用，分别是经验和学位分数，并向我们展示了如果我们知道一个结果背后的原因之一，就会降低另一个原因的重要性。换句话说，我们在观察候选人的经历时解释了成绩差的原因。这种现象俗称解释走。

下图显示了因果推理中涉及的节点之间的交互路径:

贝叶斯网络通常是将独立事件画在上面，影响从上到下流动(类似于求职面试的例子)。回忆因果推理从上到下流动，证据推理从下到上流动，因果间推理横向流动，这可能是有用的。

D-分离

了解了箭头的方向表示在贝叶斯网络中一个节点可以影响另一个节点之后，让我们看看在贝叶斯网络中如何确切地影响流量。我们可以看到，成绩最终会影响工作机会，但是在一个非常大的贝叶斯网络的情况下，说明叶节点受到贝叶斯网络顶部所有节点的影响是没有帮助的。是否存在影响不流动的情况？我们将在下表中看到解释影响力流动的简单规则:

|

没有观察到变量

|

观察到了 y

|
| --- | --- |
| | |
| | |
| | |
| | |

上表描述了三个节点 X 、 Y 和 Z 之间的开放和封闭活动轨迹。在第一列中，没有观察到变量，而在第二列中，观察到了 Y 。

我们将首先考虑没有观察到随机变量的情况。考虑上表第一列中的节点链。请注意，前三行中的规则允许影响从第一个节点流向最后一个节点。

影响可以沿着边缘的路径流动，即使这些链延伸了更长的序列。必须指出，影响的流动不受连接它们的环节的方向性的限制。

但是，有一种情况我们要注意，那就是所谓的 V 型结构，——之所以这么叫，可能是因为边缘的方向指向内部。

在这种情况下，影响不能从 X 流向 Z，因为它被 y 阻挡。在更长的链中，影响将流动，除非它被 V 形结构阻挡。

在这种情况下，由于在 X 处的 V 型结构，影响可以从和流出，但不能穿过节点 X

我们现在可以陈述活动轨迹(影响力)的概念。如果没有观察到证据，踪迹如果不包含 V 型结构，则是活动的。如果两个变量之间存在多个轨迹，如果没有任何轨迹是活动的，则它们是有条件独立的。

现在让我们看看第二种情况，我们确实观察到了证据变量。如果我们将案例与之前的链条进行比较，更容易理解，在这里我们观察到随机变量 Y 。

上表中显示的笑脸轨迹表示活动轨迹，其他轨迹表示受阻轨迹。可以观察到，证据的引入只是否定了先前活跃的踪迹，如果先前阻塞的 V 型结构存在，它就会打开。

我们现在可以声明，如果证据集中存在中间节点或其任何 V 型结构的后代(例如， Y 或其后代在)时，给定证据的踪迹将是活动的。换句话说，观察 Y 或者它的任何一个子体，都会打开阻塞的 V 型结构，使其成为一个活跃的踪迹。此外，从下表中可以看出，公开的线索会因证据的引入而受阻，反之亦然。

D-分离的例子

在本节中，我们将使用求职者示例来了解 D-分离。在执行因果推理的过程中，我们将查询工作机会，并将在工作机会的父项中引入观察到的变量，以验证我们在上一节中看到的活动轨迹的概念。以下代码来自D-separation.ipynb IPython 笔记本。

我们首先在没有其他观察变量的情况下查询工作机会，如以下代码所示:

getTableCPD().specificquery(dict(Offer='1'),dict())

前面代码的输出如下:

0.432816

从主动跟踪规则中我们知道，观察Experience应该会改变报价的概率，如下代码所示:

getTableCPD().specificquery(dict(Offer='1'),dict(Experience='1'))

前面代码的输出如下:

0.6438

根据输出，它发生变化。现在，让我们添加Interview观察变量，如下代码所示:

getTableCPD().specificquery(dict(Offer='1'),dict(Interview='1'))

前面代码的输出如下:

0.6

我们得到的Offer概率略有不同。从 D 分离规则中我们知道，观察Interview应该会阻挡从Experience到Offer的活动轨迹，如下代码所示:

getTableCPD().specificquery(dict(Offer='1'),dict(Interview='1',Experience='1'))

前面代码的输出如下:

0.6

观察Offer的概率从0.6没有变化，尽管添加了Experience变量被观察。我们可以添加Interview对象的父变量的其他值，如下面的代码所示:

query=dict(Offer='1')
results=[getTableCPD().specificquery(query,e) for e in [dict(Interview='1',Experience='0'),
dict(Interview='1',Experience='1'),
dict(Interview='1',Grades='1'),
dict(Interview='1',Grades='0')]]
print results

前面代码的输出如下:

[0.6, 0.6, 0.6, 0.6]

前面的代码显示，一旦观察到Interview变量，Experience和Offer之间的活动轨迹就会被阻断。因此Experience和Offer在给出Interview的时候是有条件独立的，也就是说观察面试父母的价值观，Experience和Grades并不会有助于改变报价的概率。

阻断和解除阻断一个 V 型结构

让我们来看看网络中唯一的 V 型结构，体验 → 面试☆等级，看看观察到的证据对活跃的踪迹有什么影响。

getTableCPD().specificquery(dict(Grades='1'),dict(Experience='0'))
getTableCPD().specificquery(dict(Grades='1'),dict())

前面代码的结果如下:

0.3
0.3

根据的 D-分离规则，面试节点在Experience和Grades之间是一个 V 型结构，它阻挡了它们之间的主动踪迹。前面的代码显示了观测变量Experience的引入对等级的概率没有影响。

getTableCPD().specificquery(dict(Grades='1'),dict(Interview='1'))

以下是前面代码的输出:

0.413016270338

以下代码应激活Experience和Grades之间的轨迹:

getTableCPD().specificquery(dict(Grades='1'),dict(Interview='1',Experience='0'))
getTableCPD().specificquery(dict(Grades='1'),dict(Interview='1',Experience='1'))

前面代码的输出如下:

0.588235294118
0.176470588235

前面的代码现在显示了在Experience和Grades之间存在活动轨迹，其中改变观察到的Experience值会改变Grades的概率。

因式分解和 I-映射

到目前为止，我们已经理解了图 G 是分布 P 的表示，我们可以用下面的方式正式定义图 G 和分布 P 之间的关系。

如果 G 是随机变量上的图，我们可以说分布 P 在 G 上分解，如果。这里，是的父节点。换句话说，联合分布可以定义为每个随机变量的乘积，当它的父变量被给定时。

因式分解和独立性之间的相互作用是一个有用的现象，它允许我们声明，如果分布在一个图上因式分解，并且给定两个节点是 D 分离的，则分布满足这些独立性()。

或者，我们可以声明图 G 是分布 P 的独立图(I-图)，如果 P 在 G 上进行因子分解，那么我们可以从图中读取独立，而不管参数如何。一个 I-map 可能不会编码分布中的所有独立性。但是，如果该图满足分布中的所有依赖关系，则称之为完美图 ( P 图)。工作面试的图就是一个 I 图的例子。

朴素贝叶斯模型

我们可以这样总结:从以下两个角度可以看到一个图:

• 因式分解:这是图允许分布被表示的地方
• I-map :这是图编码的独立性在分布中存在的地方

朴素贝叶斯模型是一个简单的独立性假设。我们使用朴素贝叶斯模型来执行二进制分类这里，我们被给予一组实例，其中每个实例由一组特征和一个类组成。分类中的任务是预测的正确类别，而其余的特征...都给了。

例如，给我们一组来自两个新闻组的新闻组帖子。给定一个特定的帖子，我们想要预测这个特定的帖子来自哪个新闻组。每个帖子都是一个由一袋单词组成的实例(我们假设单词的顺序无关紧要，只考虑单词的有无)，因此特征表示单词的有无。

这里，我们将把朴素贝叶斯模型看作一个分类器。

朴素贝叶斯和求职者示例之间的区别在于，朴素贝叶斯之所以这么叫，是因为它做出了天真的条件独立性假设，并且模型将先验概率和单个条件概率的乘积分解，如下式所示:

虽然左边的术语是需要大量独立参数的联合分布(2 ⁿ⁺¹ - 1，如果每个特征都是二进制值)，但是右边的朴素贝叶斯表示只需要 2n+1 个参数，从而将参数的数量从指数(在典型的联合分布中)减少到线性(在朴素贝叶斯中)。

在新闻组示例的上下文中，我们有一组从alt.atheism和sci.med新闻组中提取的单词，如{无神论者、医学、宗教、解剖学} 。在这个模型中，你可以说每个单词出现的概率是只依赖于类别(即新闻组)，而不依赖于帖子中的其他单词。很明显，这是一个过于简化的假设，但是已经证明它在功能数量大而实例数量少的领域中有相当好的性能，比如文本分类，我们将在 Python 程序中看到。

一旦我们看到特征之间的强相关性，分层贝叶斯网络就可以被认为是朴素贝叶斯模型的进化版本。

朴素贝叶斯的例子

在朴素贝叶斯示例中，我们将使用 Scikit-learn 的朴素贝叶斯实现——一个机器学习库来对新闻组帖子进行分类。我们从 Scikit-learn 提供的数据集中选择了两个新闻组(alt.atheism和sci.med)，我们将使用朴素贝叶斯来预测特定帖子来自哪个新闻组。以下代码来自Naive Bayes.ipynb文件:

from sklearn.datasets import fetch_20newsgroups
import numpy as np
from sklearn.naive_bayes import MultinomialNB
from sklearn import metrics,cross_validation
from sklearn.feature_extraction.text import TfidfVectorizer

cats = ['alt.atheism', 'sci.med']
newsgroups= fetch_20newsgroups(subset='all',remove=('headers', 'footers', 'quotes'), categories=cats)

我们首先使用 Scikit-learn 提供的实用函数加载新闻组数据(这将从互联网下载数据集，可能需要一些时间)。newsgroup对象是一张地图，新闻组帖子保存在data键下，目标变量在newsgroups.target，如下代码所示:

newsgroups.target

前面代码的输出如下:

array([1, 0, 0, ..., 0, 0, 0], dtype=int64)

由于特征是单词，我们使用术语频率-逆文档频率 ( Tfidf )将它们转换为另一种表示。Tfidf 的目的是不强调所有帖子中出现的单词(如“The”、“by”和“for”)，而是强调特定类别独有的单词(如宗教和神创论，来自alt.atheism新闻组)。我们可以通过创建一个TfidfVectorizer对象，然后将所有新闻组数据转换为矢量表示来完成同样的工作，如下面的代码所示:

vectorizer = TfidfVectorizer()
vectors = vectorizer.fit_transform(newsgroups.data)

向量现在包含我们可以用作朴素贝叶斯分类器的输入数据的特征。一个形状查询显示它包含 1789 个实例，每个实例包含大约 24000 个特征，如下面的代码所示。然而，这些特征中的许多可以是 0，这表明这些单词没有出现在那个特定的帖子中:

vectors.shape

前面代码的输出如下:

(1789, 24202)

Scikit-learn 提供了几个版本的朴素贝叶斯分类器，我们将使用的那个叫做MultinomialNB。由于使用分类器通常包括将数据集分为训练集、测试集和验证集，然后在训练集上进行训练，并在验证集上测试有效性，因此我们可以使用 Scikit-learn 提供的实用程序为我们做同样的事情。cross_validation.cross_val_score函数自动将数据分割成多个集合，并返回f1分数(衡量分类器准确性的指标)，如以下代码所示:

clf = MultinomialNB(alpha=.01)
print "CrossValidation Score: ", np.mean(cross_validation.cross_val_score(clf,vectors, newsgroups.target, scoring='f1'))

前面代码的输出如下:

CrossValidation Score:  0.954618416381

我们可以看到，尽管假设当给定类时，所有特征都是有条件独立的，但分类器保持了 95%的良好f1分数。

总结

在这一章中，我们学习了条件独立性如何允许联合分布被表示为贝叶斯网络。然后，我们参观了推理的类型，了解了影响如何通过贝叶斯网络流动，并使用 Libpgm 探索了相同的概念。最后，我们使用一个简单的贝叶斯网络(朴素贝叶斯)来解决一个真实世界的文本分类问题。

在下一章中，我们将学习无向图模型或马尔可夫网络。

三、无向图模型

我们现在来看看模型，其中变量之间的相互作用使得方向性不能归因于它们。我们将看到如何根据它们的分布和 D-分离来表示这些无向模型及其性质。

成对马尔可夫网络

我们将用下面的例子来推动本章的讨论。我们有四个同事( Amar 、 Bob 、 Charlie 和 Deepak )在办公室聚会上相遇。聚会上有传言说一名高级经理可能会离开组织。这四位同事偶遇，一对一对话，可能是关于谣言。我们知道查理和鲍勃·邦德真的很好，但查理和迪帕克就不一样了。简而言之，每一对都会影响另一对，对彼此有一定程度的喜欢或不喜欢。我们希望该模型能够捕捉每个人是否知道或不知道谣言的实例。在下图的图表中可以看到它们相互作用的模型。我们将使用马尔可夫网络来模拟这四个同事是否知道这个谣言。

在建模问题时，您可能会想起贝叶斯网络，我们可以尝试使用它来建模聚会示例。看起来变量之间的相互作用是对称的；人们不能把一个方向归因于两个人之间的亲密关系。我们可以使用马尔可夫网络对这种类型的交互进行建模，其中节点代表随机变量，边代表两个节点之间的亲缘关系、交互或喜欢/不喜欢。

两个变量之间的相互作用可以用一个称为因子的函数来表示，该函数将量化其范围内变量之间的好恶程度。因子有一组输入变量，称为因子的范围。因子中的值是实数集中的，其中零表示亲和力低，反之亦然。

我们可以用因子和来表示前面的网络，它们的值可以分别被认为是对 Amar 和 Bob、Bob 和 Charlie、Charlie 和 Deepak、Deepak 和 Amar 的相容因子(或局部快乐)。下表中的值描述了 Amar 和 Bob bond 都很好地分享了信息，因此他们两个知道或不知道谣言的可能性很高，一个人知道而不与另一个人分享的可能性很低:

|

|
| --- |
| | B₀T3】 | 20 |
| | B₁T3】 | 2 |
| | B₀T3】 | 5 |
| | B₁T3】 | 10 |

所以，我们可以看到图被分解成小的因素，但是我们如何定义整体的联合分布呢？我们将通过使用定义如下的因子乘积来实现这一点:

然而，你可能会说，这不是一个有效的概率分布，因为它不等于 1，这是正确的。这是通过将所有项除以来纠正的，T1 也被称为分区函数，定义如下:

这里， Z 简单来说就是所有因子的值之和，如下式所示:

既然我们已经了解了贝叶斯网络，一个自然的问题可能会出现——因素是某种边际概率还是某种条件概率或者它们的某种组合？不幸的是，这些边际概率或条件概率与因子概率之间没有直接的对应关系。我们所能陈述的是，因素是构成马尔可夫网络的其他因素的集合。因此，在马尔可夫网络编码的概率分布和网络构成因素之间没有直接的映射，这与我们在贝叶斯网络中遇到的情况略有不同。

我们可以将成对马尔可夫网络定义为由表示随机变量的节点组成的无向图，这些节点由两个节点()之间的边连接。这表示两个节点之间的因子(或势)，即。

吉布斯分布

在成对马尔可夫网络中，我们的因子在其范围内只有两个变量。如果我们扩展网络，使每个节点连接到网络中的所有其他节点，我们将获得如下图所示的边:

成对马尔可夫网络中的是否代表随机变量上的任何概率分布否，因为该网络中的参数不足以描述 n 随机变量上的联合概率分布(如果每个随机变量取 d 值，则需要参数)。细心的读者可能已经注意到对于前面的成对马尔可夫网络(有四个变量)，参数的数量(与联合分布相比)确实相等。让我们看看下表，了解当每个随机变量取二进制值时所需的参数数量:

|

随机变量的数量

|

边数(因素)

|

成对马尔可夫网络中的参数

|

联合分布中的参数

|
| --- | --- | --- | --- |
| 4 (2 x 2 网格) | four | Sixteen | =16 |
| 9 (3 x 3 网格) | Twelve | Forty-eight | =512 |

我们可以看到，在一般情况下，成对马尔可夫网络中的参数数量不足以表示联合分布。解决这个问题的一个方法是允许一个因子在其范围内有无限数量的变量，而不是两两两的情况。这被称为吉布斯分布，定义如下:

每个因子在其范围内可能有几个变量，吉布斯分布就是一组这样的因子。正如我们已经看到的成对情况，非标准化分布是其因素的乘积，如下式所示:

为了使它成为一个加起来为 1 的概率分布(因为因子的取值不受约束)，我们必须通过除以配分函数 Z 进行归一化。我们得到，这里是。

那么，这个吉布斯分布到底给了我们什么呢？它只是给了我们一个有效的概率分布(加起来是 1)，它看起来类似于一个概率分布，可以映射到一个图，就像贝叶斯网络一样。

诱导马尔可夫网络

如果我们有涉及多个变量的因素，马尔可夫网络会是什么样子？在成对案例中，因素意味着随机变量 A 和 B 之间的边缘。假设因子在所有对中具有优势确实是合乎逻辑的。这叫做诱导马尔可夫网络，图中包含两个因子和的边，如下图所示:

因式分解

在贝叶斯网络的表示部分，我们遇到了两种表示:一种是图，另一种是概率分布。我们问自己一些问题，比如它们是否等价，当我们从一种观点切换到另一种观点时，我们是失去了信息还是获得了信息？我们现在将在马尔可夫网络的背景下研究这些问题。

分布 D 和图 G 的等价性是什么？D 什么时候对 G 进行因子分解？换句话说，什么时候可以用 G 来表示 D？理解因式分解的一种方法是将其视为分解问题。我们有一个问题(例如，一个巨大的联合分布)，我们想把它分解成更小的部分(比如贝叶斯网络中的条件概率分布)。

如果我们有一组因子(它是它的单个因子的乘积)，我们可以说分布 D 在 G 上分解，G 是这组因子的诱导图。

然而，与贝叶斯网络不同的是，如果给我们一个图 G，就没有一组唯一的因子可以从图中读取。让我们看一下下图:

我们知道因子范围内的所有变量(对)都是由边连接的。例如在因子中，之间有一条边。所有其他因素也是如此。

现在考虑两组因素:和。我们可以看到，在这两个集合中，每个单个因子都满足其范围内的所有变量都由一条边连接的条件。如果我们用四个顶点画一个图，然后按照因子集或的描述连接边，我们会发现实体化的图是一样的。

所以不同的组因子可以归纳出同一个图，不能从图中读出因式分解，或者拆分成单个因子。这不应该被理解为与贝叶斯网络相比不够好，而是一个图可以有多个因子分解，所有这些都是正确的。

影响流

影响如何沿着马尔可夫网络流动？我们已经看到，在给定的图 G 中不存在单一唯一的因子分解，影响的流不依赖于因子的形式。

两个变量可以相互影响，只要它们通过一组边连接起来。在上图中， Amar 、 Bob 、 Charlie 通过一组边连接在一起，可以相互影响，因式分解，不管是还是，其实并不重要。

主动跟踪和分离

与贝叶斯网络中活动轨迹的详细规则不同，在马尔可夫网络中，规则很简单:影响可以在任意两个通过边连接的随机变量之间流动。如果观察到某个节点，则该节点会被阻止。我们可以在下图的示例中看到这一点。

因此，我们可以陈述马尔可夫网络中分离的概念(注意，我们不称之为 D-分离，因为这不是一个有向图)。给定证据集 Z，如果在 X 和 Y 之间的 G 中没有活动踪迹，则随机变量 X 和 Y 在图 G 中是分开的。也就是说，在连接 X 和 Y 的踪迹之一中应该没有证据或观察到的节点。

结构化预测

在机器学习分类器的典型应用中，分类器预测目标变量的类别，例如垃圾邮件/非垃圾邮件，以对电子邮件进行分类。通常，每个实例(如单个电子邮件)都独立于下一封或上一封电子邮件。

然而，有几类应用的目标变量与其邻居相关。以图像分割为例，其中一个典型的问题是，在牧场上的牛的图片中，我们希望将每个像素(或超级像素，如图像处理文献所述，它是一组连续的像素)分类为牛或草。每个超级像素都有几个邻居，如果超级像素 Sp 在牛的胃部区域，被其他已经归类为牛的超级像素包围，很明显超级像素 Sp 也应该归类为牛。超级像素 Sp 的分类任务通过使用局部结构而不是通过独立分类，忽略附近的超级像素而变得容易得多。

同样，在自然语言处理领域，名为词性标注 ( 词性标注)的任务涉及用名词、代词或动词等词类标注句子中的每个单词。在这里，如果我们看整个句子，而不是孤立地看每个单词，将一个单词标记为名词或动词会容易得多。

在这种情况下，当我们希望执行特定任务的预测，并且存在我们可以利用的局部结构时，我们可以使用诸如隐马尔可夫模型 ( HMMs )和条件随机场 ( CRFs )的模型。局部结构可以是位于单词序列中的附近的单词或位于像素网格中的附近的像素。

相关特征问题

如果我们试图预测超级像素的类别，特征(例如每个超级像素的颜色直方图)彼此高度相关。在像朴素贝叶斯这样的生成模型中，同一个特征被计数多次，这就给了一个高度倾斜的概率给一个类。校正独立性假设需要在要素之间添加边，因为相关的要素不是条件独立的。很难弄清楚它们是如何关联的，模型之间会紧密连接，这意味着很多顶点之间会有很多边。

这个问题的一个解决方案是建模而不是()，也就是建模条件分布而不是联合分布(其中 Y 是目标变量的集合， X 是观测变量的集合)。我们不关心特性如何变化(或者不一起变化)；我们只关心目标变量。通用报告格式不同于 MRF，它试图模拟条件分布而不是联合分布()。

通用报告格式表示

一个 CRF 可以用一组因子来表示，看起来类似于吉布斯分布。不同的是归一化常数 Z 只对 y 的所有值求和

如果联合分布有两个变量 A 和 B，如下式所示:

然后通过除以对条件分布进行归一化，如下式所示:

这里，是 A 值的联合分布，求和出来。

朴素贝叶斯模型是生成模型的一个例子(模拟联合分布的模型)，它的条件等价物是逻辑回归模型(模拟条件分布的模型)。

通用报告格式示例

通过这个例子可以更好地理解通用报告格式的使用:

给定一组与文字中手写字符相对应的图像，我们需要预测每个图像中的字符。我们可以想象，我们构建的解决方案将从最简单的开始，然后努力改进它。我们可以从一个最简单的模型开始，它有单因素，也就是说，每个角色的预测值将只取决于图像特征，如下图所示。已经观察到图像 I1 到 I4 ，我们想识别类变量 c1 到 c4 (可以取 26 个值，从 A 到 Z 每个字符一个)。

通过使用强分类器，如逻辑回归或随机森林，我们可以使用单例特征获得相当好的预测分数。

但是，我们可以看到这种模式是可以改进的。考虑上图中的字符 l2 。我们知道字母 q 后面最常见的字符是 u ，其他字符极其少见。同样，跟在 qui 后面的字符很可能是 t 或 d (quid 或 quit)。我们可以通过添加成对特征(在两个字符类之间)、三元组特征(在三个字符类之间)甚至更多来编码知识。观察下图，该图描绘了添加的特征的边缘:

捕捉类 C 和图像 I 之间的势的单因素可以表示为，成对和三元因素可以表示为和。英语字符的成对特征适用于所有组合，如{A - A}、{A - B}、{A - C}一直到{Z - Z}，同样适用于三元组因子{ A-A }、{ A-A-B }，等等。假设我们在英语中有 26 个字符，我们将分别需要 26 个 ₂ 和 26 个 ₃ 因子来表示成对和三元组因子的所有组合。我们必须从数据中了解这些因素的参数。例如，为了学习{X-Y-Z}的参数，我们可以计算字符{X-Y-Z}在一个大型单词语料库中的出现次数。我们可能会预料到许多这样的因素极少发生(例如组合 xmq )。

添加范围内有三个以上变量的较大因子可能很有诱惑力，但考虑到指数级的增长，在因子数量呈指数级的网络上存储、操作和运行推理变得很困难。

因式分解-独立探戈

在前一章中，我们了解到图和分布的两个概念都编码在一个图模型中。我们现在转向因子分解和独立性的等价性，我们想知道在马尔可夫网的背景下，它们在两种观点中是否都受到尊重。

以下是我们希望解决的问题:

• 第一个问题是，如果图中的两个节点是有条件独立的，分布是否尊重这种独立性？
• 第二个问题是，将一个分布分解成一个图是有效的分解吗？

这个定理与贝叶斯网络有相似之处。如果一个分布 P 在图 G 上因式分解，并且假设两个随机变量 X 和 Y 被分开(在图 G 中)给定，那么分布 P 满足独立性陈述，并且 X 在给定 z 的情况下条件独立于 Y。

换句话说，由图 H 通过分离属性定义的独立性被表示该图的分布 P 所保留。类似于贝叶斯网络，图 H 是分布 p 的 I 图

这种关系的逆是因式分解的独立性——如果分布 P 满足图 G 所描述的独立性，则 P 在图 G 上进行因式分解。与贝叶斯网络的情况不同，这仅适用于正分布，也就是说，随机变量的所有赋值的概率都必须大于零。这种关系源自哈默斯利-克利福德定理。

总而言之，正如在贝叶斯网络的情况下，我们有两个马尔可夫网络的等价视图。因式分解使图 G 能够表示(正)分布，而 I-map 使图 G 描述的独立性在(正)分布中得到尊重。

总结

在这一章中，我们学习了马尔可夫随机场，它是图模型的一种替代表示。我们了解了它们的属性，它们与贝叶斯网络的相似之处和不同之处。我们还研究了一些使用磁流变液的实际应用。我们还了解了条件随机场以及它们如何用于结构化预测。

在学习了表示 PGM 的概念和方法之后，在本书的下一部分，我们将学习如何将表示转换为可以被软件工具使用来解决现实问题的工件。

四、结构学习

在使用 PGM 解决问题时，我们可能采取的常见路径是从数据集开始，然后在该数据集上运行一些类型的推理查询。我们可能有一些网络结构的领域知识，以及这个网络的参数的知识。通常，我们可能没有网络结构或参数先验的先验知识，我们所拥有的只是一个数据集。

给定一个数据集，我们手头有以下任务:

• 学习网络的结构。这可以通过单独从数据中学习来实现，也可以通过提供一些领域知识(边缘之间的联系)来实现。
• 了解网络的参数。同样，领域知识可以帮助我们获得参数先验，或者它可以完全从数据中学习。
• 使用推理引擎运行条件概率查询或映射查询。

隐含着对一种工具的需求，这种工具可以帮助我们完成前面所有的任务，结构和参数学习，以及推理。

在本书的下一部分，我们将按照同样的顺序探索结构、参数学习，最后是推理。

在本章中，我们将学习如何学习 PGM 的结构。

当我们希望使用 PGM 时，一个任务是确定图模型的结构。有时，我们有一个领域专家，他可以安排一个合适的层次模型。不幸的是，情况并非总是如此，在没有任何领域知识的情况下，我们可能希望了解给定数据的网络结构。即使有一些可用的领域专业知识，我们也希望验证或改进网络结构。

结构学习的另一个用途是当目标是发现结构时(而不一定只是为了运行推理查询而发现结构)。例如，因果贝叶斯网络已经被用于建模蛋白质信号网络，其中发现网络结构被用于理解药物相互作用和患病细胞中的功能障碍信号。

结构学习景观

结构学习之旅是这样进行的:我们讨论的算法可以分为两个领域，一个领域使用约束，另一个领域关注基于分数的方法。

我们还将讨论已发现的结构类型，如树木、森林和图表。我们将探讨每组中算法的优缺点，并尝试一些简单的例子来练习一些算法。

基于约束的结构学习

在这种方法中，我们从一组表示数据中随机变量的顶点开始，然后测试数据中的条件相关性(和独立性)。这种方法的目标是从贝叶斯网络结构中读取数据的条件依赖和独立性。约束本质上是对随机变量之间条件独立性的测试。

该算法在逻辑上可以分为三个部分。

第一部分

对于每个变量 X _i ，算法试图找到见证变量的子集(比如说， X ₁ 到 X _n )，其中 X _i 独立于其他变量。然而，检查随机变量的所有子集将需要搜索指数。为了补救这一点，该算法通过执行以下任务仅使用多项式次数的独立性测试:

它限制每个顶点的父顶点的数量。这取决于我们对网络规模的直觉。更大的网络可能需要更大的见证变量子集。

用于确定条件独立性的程序也值得讨论。这个问题通常在使用假设检验的统计文献中得到解决。假设我们在变量 A 和 B 上有一个分布，我们希望知道联合分布是否是个体概率的乘积，。变量不独立的分布可能是。

对于离散值随机变量，接受或拒绝零假设(零假设表明两个变量是独立的)的合适检验是卡方检验。当两个变量独立时，卡方检验很可能返回值 0，当两个变量不独立时，卡方检验将返回高值。皮尔逊卡方检验统计量通过使用以下公式计算:

在该公式中，参数如下:

• :这是皮尔逊累积检验统计量
• :这是观察到的频率
• :这是零假设给出的预期或理论频率
• :这是自由度

通过卡方检验，我们可以在网络中建立一组独立性。然而，当运行多个假设测试时，由于数据稀疏和统计噪声，统计误差会累积并导致不正确的独立性假设。在有大量实例的小型网络上，这些错误不会影响结构的正确性，但对于大型网络来说，这是不能说的。数据不足的大型网络可能会产生虚假的独立性，这将导致构建的结构不正确。

第二部分

众所周知，分布 P 中独立(I-map)的图不是唯一的，这意味着有多个图 G可以表示该分布，这些图可以称为 I-等价图。该算法的目标是找到属于 G的 I 等价类的任何一个网络。

在第二部分，我们将创建一个无向图或图骨架，如果 X 和 Y 节点在图中相邻，我们将在它们之间添加边。通过邻接，我们意味着最终的图将具有形式或的有向边。我们可以使用形式的独立查询来推断它们是否相邻。如果我们发现任何见证节点 U ，它会阻挡 X 和 Y 之间的活动踪迹，从而得出 X 和 Y 是有条件独立的结论。另一方面，如果我们没有找到任何见证节点 U ，我们得出结论，它们必须是条件依赖的，因此我们添加了一条连接 X 和 Y 的边。

假设我们有一组独立，我们使用前面提到的测试在不独立的节点之间添加边。

第三部分

现在我们有了一个骨架，我们希望把无向边转换成有向边。一个定理(柯勒 3.8)假设所有 I 等价于 G*的有向图都具有相同的不道德集。假设一个节点有两个父节点；如果父母之间没有边缘联系(可能类似于未婚)，那他们就是不道德的。该算法继续测试和验证骨架中潜在的不道德性。潜在的不道德，目前被表示为，必须改变为以下之一:

• 
• 
• 
• 

对于这四个目的地中的每一个，我们都有规则来确定给定一个无向三元组节点的正确目的地。有关这些规则的描述，请参考柯勒·3.4.3.2 的文章。

我们现在来看看constraint-based.ipynb IPython 笔记本，在这里我们试图通过使用基于约束的方法来学习贝叶斯网络的结构。我们将首先加载网络，如以下代码所示:

from libpgm.nodedata importNodeData
from libpgm.graphskeleton importGraphSkeleton
from libpgm.discretebayesiannetwork importDiscreteBayesianNetwork
from libpgm.pgmlearner importPGMLearner

nd=NodeData()
skel=GraphSkeleton()
fpath="job_interview.txt"
nd.load(fpath)
skel.load(fpath)
skel.toporder()

bn=DiscreteBayesianNetwork(skel,nd)

我们正在加载一个具有现有结构和参数(在job_interview.txt文件中定义)的网络，这可能看起来很奇怪。对于本例，我们将使用合成数据，即从现有网络中提取的样本。这有助于我们将结果与开始时的已知网络进行比较。首先，我们将从job_interview网络中随机抽取两个样本，我们在之前的章节中已经看到过，如下代码所示:

bn.randomsample(2)

前面代码的输出如下:

[{u'Admission': u'admitted',
u'Experience': u'high',
u'Grades': u'poor',
u'Interview': u'good',
u'Offer': u'no'},
 {u'Admission': u'admitted',
u'Experience': u'low',
u'Grades': u'poor',
u'Interview': u'poor',
u'Offer': u'yes'}]

我们可以看到，随机样本是每个随机变量的一个特定赋值，它是从联合分布中抽取的。它也可以被认为是网络中所有节点的随机分配。

让我们讨论一下算法是如何进行的。我们首先查询所有节点对的条件独立性。这是通过运行卡方检验实现的。零假设指出节点 X 和 Y 是条件独立的，给定 z

下面的discrete_condind方法返回卡方值和 p 值，这是独立是由于偶然的概率。我们为 p 值选择一个阈值(比如 0.05)。如果卡方检验统计返回大于阈值的 p 值，则意味着 X 和 Y 之间的独立概率太高，不可能是偶然发生的。所以，我们可以得出结论，X 和 Y 确实是独立的。

learner = PGMLearner()
data = bn.randomsample(200)

X,Y='Grades','Offer'
c,p,w=learner.discrete_condind(data,X,Y,[])
print "independence between X and Y: ",c," p-value",p," witness node: ",w

前面代码的输出如下:

independence between X and Y:  8184619.56996  p-value 0.0  witness node:  []

我们将从下图中描述的求职面试网络运行独立性查询:

在前面代码片段的查询中，我们可以看到 p 值小于 0.05，Grades和Offer变量可以认为是不独立的。由于 D 分离规则规定，给定工作面试网络，Grades和Offer之间有一条活跃的轨迹，如果观察到Interview变量，这条轨迹就会被阻断。如果我们观察位于Grades和Offer之间的Interview变量会发生什么？：

X,Y='Grades','Offer'
c,p,w=learner.discrete_condind(data,X,Y,['Interview'])
print "Independence between X and Y: ",c," p-value",p," witness node ",w

前面代码的输出如下:

现在，Grades和Offer是条件独立的(因为 p 值远大于 0.05)。在给定其他见证变量的情况下，算法的第一阶段基本上试图确定网络中所有节点对的条件独立性。然后，我们剩下一组节点之间的无向依赖关系。

Independence between X and Y:  2.79444519518  p-value 0.993172910586  witness node  ['Interview']

算法的第二和第三阶段基本上包含在discrete_constraint_estimatestruct方法中，其中依赖集被转换成无向图。然后，方向性被解决。这在以下代码中显示:

result=learner.discrete_constraint_estimatestruct(data)
print result.E

前面代码的输出如下:

[[u'Grades', u'Admission'], [u'Experience', u'Interview'], [u'Grades', u'Interview'], [u'Interview', u'Offer']]

result.E变量返回网络中学习到的有向边，每条边都是包含该边的起点和终点节点的列表。我们可以看到，求职面试网络中的四条原始边已经找到了。然而，这是一个小网络。这种算法会扩展到更大的网络吗？

我们现在来看看alarm_network.ipynb IPython 笔记本，其中我们尝试在更大的数据集上使用基于约束的方法。逻辑警报减少机制网络是一个贝叶斯网络，旨在提供一个警报消息系统来监控患者。这个网络有 37 个顶点和 46 条边，比我们目前使用的 5 个顶点和 4 条边的求职面试网络要大得多。

数据集可以在http://www.cs.ru.nl/~peterl/BN/alarm.csv找到，通常称为报警网络。关于数据集的更多信息(如列描述)可在http://www.bnlearn.com/documentation/man/alarm.html找到。

让我们使用pandas库加载alarm.csv文件，如下代码所示:

import pandas as pd
import numpy as np
df=pd.read_csv("alarm.csv")

alarm.csv文件中有记录，我们应该将其转换为libpgm可以使用的格式。每个实例应该是一个字典，其中键是列名，值是列值。下面的函数执行格式转换，并返回一个字典列表，这些字典是从原始数据集中采样而不替换的。这在以下代码中显示:

from random import randint,sample

def rand_index(dframe,n_samples=100):
rindex =  np.array(sample(xrange(len(dframe)) ,n_samples if n_samples<=len(dframe) else len(dframe)))
    return [{i:j.values()[0] for i,j in dframe.iloc[[k]].to_dict().items()} for k in rindex ]

#Lets examine a single sample:
rand_index(df,n_samples=1)

前面的代码给出了以下输出:

 [{'Anaphylaxis': 'b',
 'ArtCOb': 'b',
 'BP': 'c',
 'CO': 'a',
 'CVP': 'a',
 ….<rows elided>
}]

让我们加载数据，创建一个学习者对象的实例，并用少量样本(100 个)估计结构，如下面的代码所示:

from libpgm.nodedata import NodeData
from libpgm.graphskeleton import GraphSkeleton
from libpgm.discretebayesiannetwork import DiscreteBayesianNetwork
from libpgm.pgmlearner import PGMLearner

data=rand_index(df,n_samples=100)
learner = PGMLearner()
result=learner.discrete_constraint_estimatestruct(data)
print result.E

前面代码的输出如下:

[['VentLung', 'MinVol'], ['HR', 'HRBP'], ['LVFailure', 'History'], ['BP', 'TPR'], ['SaOb', 'Shunt'], ['Intubation', 'Shunt'], ['PCWP', 'LVEDVolume'], ['Hypovolemia', 'LVEDVolume']]

为了比较libpgm学习的结构与正确报警网络的性能，让我们加载parent-child.txt文件中描述的正确报警网络。该文件的每一行都包含父顶点，后跟子顶点(在某些情况下，某些节点是没有子节点的叶节点)。这在以下代码中显示:

file = open('parent-child.txt', 'r')

def edges(line):
st=line.strip('\n').strip(' ').split(' ')
    #print st
    return [[st[0],i] for i in st[1:] ]

all_edges=[l for line in file for l in edges(line)]
#a set containing the correct edges
ground_truth=set([tuple(i) for i in all_edges])
print all_edges[:5]

前面代码的输出如下:

[['HISTORY', 'LVFAILURE'], ['CVP', 'LVEDVOLUME'], ['PCWP', 'LVEDVOLUME'], ['LVEDVOLUME', 'HYPOVOLEMIA'], ['LVEDVOLUME', 'LVFAILURE']]

让我们定义一个诊断函数，将找到的网络与正确的网络进行比较。我们希望查询两者中的边数，以及具有正确方向的边数和正确连接节点但方向错误的边数。代码如下:

def printdiag(result):
    found=set([tuple([j.upper() for j in i]) for i in result.E])
    correct=ground_truth.intersection(found)
    undirected_common_edges=[(i,j) for i,j in found for k,l in ground_truth if i.find(k)!=-1 and j.find(l)!=-1]
    print "Number of edges in learnt network ",len(found)
    print "Total number of edges in true network ",len(ground_truth)
    print "Number of edges with correct directionality ",len(correct)
    print "Number of edges with incorrect directionality ",len(undirected_common_edges)

让我们定义一个函数，当数据集选取特定数量的样本时，该函数打印结果统计，如下所示:

def learn_structure(n_samples):
    data=rand_index(df,n_samples)
    learner = PGMLearner()

result1=learner.discrete_constraint_estimatestruct(data)
printdiag(result1)

learn_structure(1000)

前面代码的输出如下:

Number of edges in learnt  36
Total number of edges in true network  46
Number of edges with correct directionality  6
Number of edges with incorrect directionality  7

让我们检查一下learn_structure方法在提供 1000 个样本时的性能，如下面的代码所示。它找到的边不到一半，只有几条边正确地连接了节点。增加样本量有帮助吗？(这可能需要几分钟来运行。):

learn_structure(10000)

前面代码的输出如下:

Number of edges in learnt network  30
Total number of edges in true network  46
Number of edges with correct directionality  2
Number of edges with incorrect directionality  3

我们可以看到，样本的数量仍然不足以提高正确识别的边的数量。算法没有学习到正确的结构并不奇怪，因为复杂性呈指数增长。精确方法的可行节点数据说约为 30 个(卡西欧·德·坎波斯，ICML，2009)。算法的复杂度为O(n<sup>d+2</sup> )运行时间，其中n为顶点数，d为见证集的上限。该算法的各种版本试图约束见证集以提高其性能(这在第 2 章、定向图模型(http://arxiv.org/pdf/1111.6925.pdf)中有所描述)。

基于约束的方法概述

基于约束的方法是一种可以学习有向和无向图模型的结构的方法。这些方法对数据中的噪声相当敏感，会产生错误的独立性假设，导致网络错过寻找 I-等价结构的机会。它们不能扩展到具有大量节点的网络，这可以通过将边界设置为见证集的大小来改善。

基于分数的学习

基于分数的方法是分配一个分数，该分数指示一个图与数据的匹配程度，然后在所有可能的网络结构的空间中搜索，以找到一个最大化分数的图。基于分数的方法加强了稀疏性(更少的边)。这就产生了一个优化问题，我们需要搜索指数级的网络结构。

可能性得分

给定特定的图结构 g，似然分数是使数据的似然性最大化的分数。最大似然估计 g 的参数。通常，使用似然分数的对数。可能性分数分解为以下公式:

等式的右侧由两项组成:

第一项是图中节点与其父节点之间的互信息之和(由表示)。互信息是一个信息论标准，可以理解为两个随机变量X₁T9】和X₂T13】的联合分布和边际分布的乘积之间的平均距离。**

第二项是每个随机变量的熵之和。可以观察到它独立于它的双亲(因此独立于图结构)。 M 为样本数。

可能性得分背后的直觉是，如果一个变量与其父变量相关，那么网络的结构会更好。换句话说，得分较高的图会将变量及其父变量放在一起。

让我们比较一下包含两个随机变量的简单网络的得分，如下图所示:

为了比较两个网络之间的分数，我们可以减去它们的分数(使用前面的等式)。

我们可以看到，当我们减去分数后，(第二个)熵项抵消，只剩下第一个项，这是 X 和 Y 之间的互信息。这显示在以下公式中:

互信息得分总是大于或等于零。当两个随机变量 X 和 Y 独立时，互信息应该为零。然而，在现实世界中，我们可能会在数据集中遇到统计噪声，给定足够的样本，随机变量之间总是会有一些依赖关系，从而导致互信息的非零值。

为什么这对我们很重要？这很重要，因为互信息项奖励网络增加一条边，很少奖励网络移除一条边，因此密集连接的网络比稀疏连接的网络获得更高的分数。

这种现象被称为过拟合，正如我们从机器学习文献中所知，这是需要避免的(因为它使 PGM 中的参数学习变得复杂)。避免过拟合的方法包括限制父对象和参数的数量。其他评分函数(我们将在后面的章节中看到)规定惩罚过度拟合，并且可以在足够连接的网络和它与数据的拟合之间建立折衷。

贝叶斯信息准则得分

试图提高可能性分数的评分函数试图以各种方式惩罚复杂性。贝叶斯信息标准(【BIC】)和阿卡克信息标准(【AIC】)是衡量统计模型与数据拟合质量的信息论得分。BIC 评分采用的方法是从可能性(第一项)中减去惩罚项(以下公式中的第二项)，如下式所示:

这里， M 是训练实例的数量，Dim【G】是分布中独立参数的数量(在联合分布中，独立参数的数量是联合分布中的总行数减一)。第二项试图将似然项添加边的倾向限制为参数和训练实例数量的函数。

第二个术语只是平衡模型复杂性和它对数据的适应性之间的折衷的一种方法。如果把换成 1 ，就叫 AIC。BIC 评分的否定被称为最小描述长度 ( MDL )标准。

已经表明，BIC 分数具有一种称为一致性的属性。随着数据实例的数量逐渐增加，图 G 的正确结构(或任何其他 I 等价结构)获得最大得分。

贝叶斯得分

另一类评分函数使用贝叶斯原理对可能的结构空间和可能的参数进行平均，以找到最适合数据和参数的图结构。像BDE和BDEU(贝叶斯狄利克雷的变体)这样的评分函数存在于贝叶斯保护伞下，但是它们使用不同的结构和参数优先级。

与所有贝叶斯方法一样，当数据实例数量较低时，先验(结构和参数)将产生较大的影响。但是，随着实例数量的增加，数据量会降低先前的强度。因此，对于少量的数据实例，贝叶斯得分低于数据。

贝叶斯方法的渐近行为产生了类似于 BIC 分数的分数，并且它们很少过度拟合数据。

在下面的笔记本(bnfinder.ipynb)中，我们现在将使用BNFinder库，使用基于分数的方法来学习贝叶斯网络的结构。

BNFinder工具实现了多种基于分数的方法，如BDE、MIT和MDL。可以使用pip install BNFinder命令安装 BNFinder。

在使用数据集找到关系之前，BNFinder工具允许我们将诸如潜在父母和已知边等先验信息添加到模型中。

让我们加载必要的导入并创建我们将使用的评分函数(BDe)。我们将使用默认参数值，如以下代码所示:

fromBNfinder.BDEimportBDE
fromBNfinder.dataimportdataset

score=eval("BDE")(data_factor=1.0,chi_alpha=.9999,sloops=False)
score

前面代码的输出如下:

<BNfinder.BDE.BDE instance at 0x00000000124EB688>

我们创建了一个数据集对象，并从一个文件中加载它。job_interview_samples.txt文件包含从libpgm'srandomsample()方法生成的样本。生成样本的代码可以在 IPython 笔记本作业_interview_samples.ipynb中找到。

BNFinder工具使用的数据集格式是普通数据集的转置。如果普通数据集有三列，则每个实例行包含三个值。但是BNFinder数据集格式将有三行，一行将包含该列的所有值。

让我们创建以下方法，读入数据集，执行结构学习，并将输出保存为.bif格式，并写出 CPDs。其他进行贝叶斯推理的工具使用贝叶斯交换格式 ( BIF )。

def learn_structure(sample_data,dataset_name):
    d = dataset(dataset_name).fromNewFile(open(sample_data))
    score2,g,subpars = d.learn(score=score,data_factor=1.0)
    d.write_bif(g,dataset_name+".bif")
    d.write_cpd(g,file(dataset_name+"_cpd.txt","w"))
    return score2,g,subpars

我们现在将尝试从求职面试网络的样本数据中学习结构，如以下代码所示:

s,g,sp=learn_structure("job_interview_samples.txt","job_interview")
g

图表信息如下:

 Admission(Admission) => Grades(-),
 Experience(Experience) => Grades(+), Interview(+),
 Grades(Grades) => Admission(-), Experience(+), Interview(+),
 Interview(Interview) => Experience(+), Grades(+), Offer(+),
 Offer(Offer) => Interview(+),

图的输出是边的列表，父在左边，子在右边。加号/减号表示节点之间的正/负相关性。

虽然代码很快终止，但它找不到正确的边，除了面试-提供和成绩-录取之间的边。BNFinder工具确实提供了通过使用两种方法将我们的直觉/领域知识添加到结构发现过程中的机会:我们可以为给定顶点定义已知的父节点列表，或者指定一个调节器(来自BNFinder工具在生物网络重建中的根的技术术语)，这将约束网络仅允许列出的节点作为根节点或所有顶点的潜在父节点。

我们可以将调节器添加到样本数据文件的序言部分，该部分保存在job_interview_samples_preamble1.txt中。如果您的系统外壳中有 head 命令，以下命令应该可以工作。

preamble段第一行提示Experience和Grades在网络的根，没有父母，如下代码所示:

!head -1 job_interview_samples_preamble1.txt
#regulators Experience Grades

s,g,sp=learn_structure("job_interview_samples_preamble1.txt","job_interview")
g

图表信息如下:

 Admission(Admission) =>
 Experience(Experience) => Interview(+), Offer(+),
 Grades(Grades) => Admission(-), Interview(+), Offer(+),
 Interview(Interview) =>
 Offer(Offer) =>

这次我们看到其他种类的错误，Interview和Offer之间的边缘消失了。我们将以下行添加到preamble部分(从父宏开始)。指定Interview顶点的父节点为Experience和Grades，如下图所示:

!head -2 job_interview_samples_preamble2.txt
#regulators Experience Grades
#parents Interview Experience Grades

s,g,sp=learn_structure("job_interview_samples_preamble2.txt","job_interview")
g

图表信息如下:

 Admission(Admission) =>
 Experience(Experience) => Interview(+), Offer(+),
 Grades(Grades) => Admission(-), Interview(+), Offer(+),
 Interview(Interview) =>
 Offer(Offer) =>

获得的边缺少面试-录用的边。

我们可以添加三个约束，指定Experience和Offer的父以及监管者。代码如下:

!head -3 job_interview_samples_preamble3.txt
#regulators Experience Grades
#parents Interview Experience Grades
#parents Offer Interview

s,g,sp=learn_structure("job_interview_samples_preamble3.txt","job_interview")
#save the file to open in cytoscape
net_str=g.to_SIF()
f=open("job_interview_sif.txt","w")
f.write(net_str)
f.close()
g

图表信息如下:

 Admission(Admission) =>
 Experience(Experience) => Interview(+),
 Grades(Grades) => Admission(-), Interview(+),
 Interview(Interview) => Offer(+),
 Offer(Offer) =>

在preamble部分指定了前面的元素之后，我们就到达了正确的网络。

我们可以使用细胞景观工具查看创建的网络，这也允许我们导出找到的节点和边。文件在细胞景观中的查看和保存是离线完成的，在这里展示。细胞角读取。sif文件创建(在前面的片段中)并可以导出一个.png图像，如下代码所示:

Image("job_net.png")

我们也可以将BNFinder工具基于分数的方法应用于更大的数据集。它比 libpgm 中基于约束的方法更有效，但它需要手握(通过添加已知边)来防止无环性并更快完成。

在下一个片段中，我们将尝试学习真实世界网络的结构。这个数据集由不同扰动下蛋白质信号网络状态的测量组成。虽然在应用结构学习算法之前，一些信令关系是已知的，但是所学习的结构被用来阐明传统报道的关系，以及以贝叶斯网络的形式推断新的网络因果关系。

以下示例摘自BNFinder项目文档:

s,g,sp=learn_structure("sachs.inp","sachs")
net_str=g.to_SIF()
f=open("sachs_cpd.sif","w")
f.write(net_str)
f.close()

以下是细胞角为网络生成的图像。该网络从真实数据集中的 17 条边中提取 11 条正确的边。

Image(filename='sachs_network.png')

我们可以得出结论，从现实世界网络中导出的数据集学习完美的网络结构是不现实的，尤其是当随机变量的数量大于一个小数字时。结构学习算法通常通过首先结合领域知识(节点之间的已知边)，然后使用数据或者经验地验证这些边的合理性，或者学习领域专家未知的新边来使用。

基于分数的学习总结

从计算的角度来看，基于分数的方法使用的优化角度是 NP 难的，相比之下，基于约束的方法具有低多项式时间复杂度。避免 NP 难度复杂性的其他方法是启发式搜索算法(如贪婪爬山、禁忌搜索和模拟退火)来执行图结构的局部搜索。这些方法使用相似的评分函数来评估结构对数据的适合度，并且它们的效率来源于要搜索的图结构的空间的减少。

由于我们经常使用学习的结构来运行推理查询，我们必须记住贝叶斯网络中的推理是 NP 难的。即使是近似推理也是 NP 难的。因此，我们努力学习图结构，这将使我们能够快速推断。因此，PGM 社区倾向于学习特定类别的结构(如树和森林)，对于这些结构，结构学习和推理都更容易管理。

树木因其理想的特性而受到青睐。即使在高维空间中，也可以使用优化方法来学习它们的结构，并且参数的数量可以保持较低。

总结

在这一章中，我们首先看了从数据中学习网络结构。我们研究了两种方法，一种基于约束，另一种基于评分函数。此外，我们运行了一些 Python 示例来查看这些算法的运行情况。我们了解到，为了找到最佳结构，我们必须在指数数量的结构上进行搜索，因此树和森林是优于一般图的结构。

在下一章中，我们将学习参数学习以及如何为我们已经知道的图结构定义参数。

五、参数学习

在本章中，我们将学习估计 PGM 参数的方法。我们从玩具的例子开始，比如估算抛硬币实验中的偏差。

到目前为止，我们在 PGM 上的旅程可以比作一个销售人员(让我们称他为杰克)试图向一家大型公司销售软件包的任务。他可能会尝试识别所涉及的不同人员，例如软件的最终用户、经理和采购部门等。这类似于在图模型中寻找随机变量(包括潜在或隐藏的随机变量)。杰克将试图与可能感兴趣的人建立联系，或者通过识别，例如，谁影响谁，谁是决策者，以及什么是组织层级。这类似于在图模型中建立随机变量之间的联系，这是我们在前一章中探讨过的。

考虑到有些人比其他人施加更大程度的影响，对杰克来说，人与人之间的影响程度也很重要。他发现购买决定是鲍勃做出的，爱丽丝和查理与鲍勃有个人联系，但鲍勃对爱丽丝的意见信任度远远高于查理(换句话说，爱丽丝对鲍勃的影响程度大于查理)。因此，如果杰克向爱丽丝提供他正在销售的软件的演示，他就有更高的销售机会，而不是查理。

确定连接之间影响程度的任务称为参数学习，我们将在本章后面讨论。有向(贝叶斯)和无向(马尔可夫)网络都需要参数学习。

我们将探索的参数学习领域包括两个广泛的领域。第一种叫做最大似然估计 ( 极大似然估计)。尽管名字听起来很复杂，但通过计算事件发生的次数(取决于分布的类型)，最大似然估计通常很容易计算。第二种方法使用贝叶斯统计。方法的选择是有争议的(在网上可以找到许多关于Bayesian和Frequentist的参考资料)，每种方法都有其优缺点。

我们从一个玩具例子开始。我们用贝叶斯和最大似然方法多次投掷图钉来学习它落在头上的概率。下图显示了图钉:

图钉的投掷可以使用伯努利分布来建模，伯努利分布是一种概率分布，取概率为的 1(正面)和概率为的 0(反面)的值，其中是图钉的偏差。

如果我们把图钉扔 100 次并记录结果，我们会得到类似于的 100 个事件的序列{ H，H，T，H，H……。} 。这一结果序列被称为独立且同分布的 ( 内径)样本，它们来自相同的分布并且彼此独立。

图钉翻转的贝叶斯网络如下图所示。(类似于第二章、有向图模型中的朴素贝叶斯网络)。

我们预计，随着足够多的翻转，我们会接近偏差的真实估计。让我们考虑以下结果的情况: {H，T，H，T} 。第一个正面的概率是θ，，第二个反面的概率是。既然 tosses 是独立的，。

因此，达到包含两个头和两个尾的序列 {H，T，H，T} 的概率为。我们可以进一步推广，说明包含头部( h )和尾部( t )的序列的概率等于。

似然函数

在概率上，我们从一个已知的参数开始，预测数据。可能性的概念是从数据开始并预测参数。最大似然是最大化数据可能性的参数值，或者换句话说，最大化似然函数值。

伯努利分布情况下的最大似然值如下:

我们看到序列的顺序没有影响，我们所需要的只是头部和尾部的数量。

最大似然估计法有一个缺点:参数估计不能表明我们对参数的信心。不管翻转的次数是多少，θ的值都可以是 0.7(10 次翻转 7 个头，10000 次翻转 7000 个头)。很明显，当结果是从 10000 次翻转而不是 10 次翻转中获得时，我们更有信心；然而，MLE 并没有给我们这样的画面。

图钉(伯努利)实验得出的头尾数称为充分统计。这是因为要得到最大似然估计，我们只需要数(头尾)而不需要序列。不同类型的分布需要不同的充分统计。例如，对于多项式分布(如掷出六面骰子)，足够的统计数据是结果为以下值之一的次数: {1，2，3，4，5，6} 。

使用最大似然估计的参数学习示例

在Chap5_thumbtack_mle.ipynb IPython 笔记本中，我们将使用最大似然估计来检查图钉的参数。

我们首先导入bernoulli配送。向后操作，我们设置图钉的偏置(通过为参数选择一个值)，然后多次翻转图钉(通过对随机值进行采样)，如下代码所示:

from scipy.stats import bernoulli
import numpy as np
import matplotlib.pyplot as plt
import matplotlib

param_theta=0.3
num_flips=100

为了生成图钉的翻转，我们只需在bernoulli对象上调用rvs函数(随机变量)。为了生成10翻转，我们运行以下代码:

bernoulli.rvs(0.3,size=10)

前面代码的输出如下:

array([0, 1, 0, 0, 1, 1, 1, 1, 1, 0])

估计参数所需的足够统计数据仅仅是正面和反面的计数。让我们编写一个函数，生成一些翻转，计算头部的数量，并通过将头部的数量除以翻转的总数，给出参数的估计值。

下面的函数返回一个元组。第一个是累积估计，第二个是最终估计。

def estimate_param(param, num_flips):
    res=bernoulli.rvs(param,size=num_flips)

    r_est=[np.sum(res[:i])/float(i) for i in range(1,len(res)+1)]

    final_estimate=ones/float(num_flips)

    return r_est,final_estimate

print estimate_param(param_theta,num_flips=10)

前面代码的输出如下:

([0.0, 0.0, 0.0, 0.0, 0.2, 0.16, 0.14, 0.125, 0.11, 0.1], 0.1)

由于翻转过程的随机性，执行该操作时看到的实际值可能会有很大不同。下面的代码确定了越来越多样本的偏差。生成的图显示，随着样本数量的增加，估计值更接近真实值(0.3)，如水平线所示。每个单独的子情节表明，随着样本数量的增加，估计值会提高。

x=[10,100,1000]

colors = matplotlib.rcParams['axes.color_cycle'] 

f, axarr = plt.subplots(len(x), 1)
f.tight_layout()
plt.subplots_adjust(hspace = 1)

for i, samples in enumerate(x):
    est,res=estimate_param(param_theta,samples)
    ax1 = plt.subplot(len(x),1,i+1)
    ax1.plot(range(samples),est , label=samples, alpha=0.4, lw=3)
    msg="number of samples ",samples

    ax1.set_title(msg)
    ax1.set_xlabel("$x$")

    ax1.set_ylim([0,0.8])
    ax1.axhline(y=res, xmin=0, xmax=1,linewidth=1, color='r')
    ax1.axhline(y=param_theta, xmin=0, xmax=1,linewidth=1, color='b')

plt.show()

前面代码的输出如下:

贝叶斯网络的最大似然估计

现在我们已经看到了如何使用最大似然估计进行估计，我们将尝试在贝叶斯网络中应用相同的方法。我们已经看到贝叶斯网络的结构如何允许我们将一个大的联合分布分解成一个更小的 CPD 集，并且在参数估计过程中同样的可分解性对我们有所帮助。

让我们以一个有两个二元随机变量的小贝叶斯网络为例， X 和 Y ，它们是这样连接的。该网络有两个由(针对 X0 和 X1)和(针对给定 X 的 Y 的 CPD 中的四个条目)参数化的 CPD。我们得到的数据集是一组表单的实例。

我们的目标(在 MLE 中)是找到最大化数据可能性的(一组)参数。

原来似然函数(即给定 D 的参数的似然)分解成两项的集合，每个随机变量 X 和 Y 各一项。每个项都被称为局部似然函数，它衡量给定父项如何预测它。

对于随机变量 X ，类似于我们在前面章节中看到的，参数值如下:

其中 M 为充分统计，例如为我们看过的次数。

同样的，四个参数的值也是通过取相似的计数得到的。的具体数值如下:

这里是我们在数据集中看到的次数，是我们看到的次数，不考虑 Y 的值。

对于更大的网络，通过使用局部似然函数独立估计每个 CPD 的参数来获得最大似然参数估计。然后可以将这些结合起来，形成整个网络的最大似然估计解决方案。

对于每个随机变量 X 加上父代 U ，MLE 参数为，其中。

使用最大似然估计的贝叶斯参数学习示例

在job_interview_libpgm.ipynb IPython 笔记本中，我们将使用最大似然估计的 libpgm 实现来学习工作面试网络中 CPDs 的参数。

下面是从libpgm加载 CPDs 的代码，如前所示:

from libpgm.graphskeleton import GraphSkeleton
from libpgm.nodedata import NodeData
from libpgm.discretebayesiannetwork import DiscreteBayesianNetwork
from libpgm.tablecpdfactor import TableCPDFactor
from libpgm.pgmlearner import PGMLearner
import pandas as pd
nd = NodeData()
skel = GraphSkeleton()
jsonpath="job_interview.txt"
nd.load(jsonpath)
skel.load(jsonpath)
skel.toporder()

我们可以创建贝叶斯网络，并使用以下代码获得随机样本:

bn = DiscreteBayesianNetwork(skel, nd)
samples=bn.randomsample(2000)

在下面的代码中，我们实例化了PGMLearner类。discrete_mle_estimateparams方法已经知道网络的结构。如前一节所述，每个持续专业发展的估计只需要来自父母的信息，这种分解使得学习每个持续专业发展的参数成为可能:

learner = PGMLearner()
result = learner.discrete_mle_estimateparams(skel, samples)

以下是来自Interview节点样本的 CPD 结果。由值0和1组成的最左边的列是来自家长分配的值([0,0])表示经验= 0 和等级= 0)，标题为0、1和2的列表示Interview节点取值0、1和2的概率。

pd.DataFrame(result.Vdata['Interview']['cprob']).transpose()

前面代码的输出如下:

 0	   1	         2
['0', '0']	 0.809582	 0.165848	 0.024570
['0', '1']	 0.321678	 0.396853	 0.281469
['1', '0']	 0.323204	 0.591160	 0.085635
['1', '1']	 0.115079	 0.182540	 0.702381

以下结果是来自原始网络的 CPD 参数。结果中的一些值是彼此相当接近的近似值，而其他值的误差在 20%的范围内。

pd.DataFrame(bn.Vdata['Interview']['cprob']).transpose()

前面代码的输出如下:

 0	  1	  2
['0', '0']	 0.8	 0.18	 0.02
['0', '1']	 0.3	 0.60	 0.10
['1', '0']	 0.3	 0.40	 0.30
['1', '1']	 0.1	 0.20	 0.70

数据碎片化

对于具有大量随机变量的网络，每个随机变量具有大量的值，我们看到每个唯一的赋值给父变量和子变量的数据实例的数量变得非常少。在某些情况下，我们根本没有实例。因此，从几个可用的数据实例中估计的参数通常会对数据进行过度填充，或者在没有可用实例的情况下返回值零。

因此，随着父集合的维数增加，可靠地估计参数变得越来越困难。这成为我们从数据中学习贝叶斯网络能力的一个限制因素。

数据碎片化对参数估计的影响

在下面的data_segmentation.ipynb IPython Notebook 中，我们将看到数据分割对使用最大似然的参数估计的影响。我们有一个在small_network.txt中定义的小网络，它有两个随机变量，，由一条弧线连接。父母 X 取 5 个值，子女 Y 取 2 个值。

我们首先从文件加载网络，创建一个DiscreteBayesianNetwork文件，如下代码所示:

from libpgm.graphskeleton import GraphSkeleton
from libpgm.nodedata import NodeData
from libpgm.discretebayesiannetwork import DiscreteBayesianNetwork
from libpgm.pgmlearner import PGMLearner

nd = NodeData()
skel = GraphSkeleton()
jsonpath="small_network.txt"
nd.load(jsonpath)
skel.load(jsonpath)
skel.toporder()

bn = DiscreteBayesianNetwork(skel, nd)

在下面的代码中，我们将创建一个函数，该函数使用从中采样的数据来学习网络的参数。抽取50样本后，我们打印赋值X = 3和Y = 0的估计参数值。我们运行同一个函数几次来比较我们得到的结果。由于采样是随机的，因此当您运行以下命令时，可能会得到不同的结果:

def learn_param(num_samp=50):
    data = bn.randomsample(num_samp)
    # instantiate learner 
    learner = PGMLearner()

    # estimate parameters from data and skeleton
    result = learner.discrete_mle_estimateparams(skel, data)
    numer=len([1 for m in data if m["X"]=='3' and m["Y"]=='0'])
    denom=len([1 for m in data if m["X"]=='3'])

    print "numerator:",numer," denominator:",denom," result=",numer/float(denom)

[learn_param() for _ in range(5)]

前面代码的输出如下:

numerator: 1  denominator: 6  result= 0.166666666667
numerator: 2  denominator: 6  result= 0.333333333333
numerator: 2  denominator: 10  result= 0.2
numerator: 1  denominator: 4  result= 0.25
numerator: 2  denominator: 12  result= 0.166666666667

我们看到结果变化很大，因为采样数据集中感兴趣样本(X == 3和Y == 0)的数量很少。我们在文件中设定的实际值是0.2。

nd.Vdata["Y"]["cprob"]["['3']"][0]

以下是前面代码的输出:

0.2

只有当我们增加样本数量时，我们才能得到接近的值，如下面的代码所示:

[learn_param(5000) for _ in range(3)]

以下是前面代码的输出:

numerator: 160  denominator: 720  result= 0.222222222222
numerator: 156  denominator: 763  result= 0.204456094364
numerator: 156  denominator: 757  result= 0.20607661823

虽然这是一个只有一个父节点的小网络，但是父节点承担了越来越多的离散值，我们得到的最大似然估计很差。通过对网络建模，可以避免数据碎片，如下所示:

• 限制给定节点的父节点数量
• 尽可能限制父节点的离散值的数量

贝叶斯参数估计

在最大似然估计的例子中，数据单独用于估计被观测的参数。然而，在许多情况下，我们对观察到的参数有一个相当好的概念。如果问我们一个硬币的公平性，我们往往相当确定参数的值是 0.5，也就是正面和反面的可能性相等。贝叶斯统计允许我们将这种先验直觉考虑在内，并找到由先验和数据共同决定的后验。即使我们认为硬币是公平的，但如果我们在 10 万次翻转中得到 3 万个头，我们将确信参数接近 0.3，而不是之前猜测的 0.5。

我们通过颠倒我们的假设开始分析，即每个翻转都是独立的，是一个固定的量。我们假设是一个随机变量，每个连续的翻转告诉我们更多关于的值。我们假设翻转是有条件独立的。

tosses 和的联合分布如下:

假设两次投掷之间有条件独立，接下来是右手项:

这里 H 和 T 分别是正面和反面的数量。

我们可以用这个联合分布来指定上的后验，如公式所示:

分子中的前两个项是先验和似然，而分母是归一化常数，使右侧成为适当的密度函数。

当然，贝叶斯方法允许我们在对概率分布有强烈直觉的情况下使用信息先验(不像均匀先验)。

参数学习的贝叶斯方法示例

在thumbtack_Bayesian.ipynb IPython 笔记本中，我们将使用贝叶斯方法确定图钉的偏差，从参数获得。

同样，我们将使用优秀的 PyMC 库。请使用easy_install PyMC安装 PyMC。

目前我们将避免讨论 PyMC 的细节。PyMC 的优秀文档可以在 PyMC 文档中找到，也可以在《黑客概率编程和贝叶斯方法》一书中找到。

为了估计θ，我们将首先从伯努利分布生成 30 个样本，其中控制参数取值0.3(对应于生成值 1 的大约 30%，其余为 0)。

我们试图找到后验分布的参数，并了解到后验既依赖于数据，也依赖于先验。我们有很多选择。在这种情况下，我们将使用统一优先。由于我们对图钉的偏差知之甚少，我们相信它可以在 0 到 1 之间。

PyMC 模型包含两个变量。第一个是统一先验，它代表了我们的信念，即参数的值可以在 0 到 1 之间的任何地方。第二个是伯努利变量，我们向它提供数据。

这两个变量以父子关系联系在一起；一致先验是伯努利变量的指定母体。

我们使用以下类型的 PyMC 变量:

• 随机:这是uni_prior变量，可以根据参数θ取不同的值
• 确定性:这个是bern变量，它的值是由其父代决定的

最后，模型中的所有变量都被包装在一个model对象中。

对于所有值不是true的变量，PyMC 的模拟将在模拟过程中略微影响变量的值。(在我们的例子中，uni_prior变量的值)将开始逼近它的后值，如下面的代码所示:

from pymc import *
from scipy.stats import bernoulli
import matplotlib.pyplot as plt
import pymc.Matplot as plott

def create_model(data):
    #create a uniform prior, the lower and upper limits of which are 0 and 1
    uni_prior = Uniform('uni_prior', lower=0,upper=1.0 )
    bern = Bernoulli('bern',p=uni_prior, value=data,observed=True)
    model=Model([uni_prior,bern])
    return model

从先验分布中抽取的样本显示先验在 x 轴上均匀分布在0和1之间，如下代码所示:

uni_prior = Uniform('uni_prior', lower=0,upper=1.0 )

samples = [uni_prior.random() for i in range(20000)]
plt.hist(samples, bins=100, normed=True )
plt.title("Prior distribution for $\ theta$")
plt.show()

前面代码的输出如下:

在下面的片段中，我们使用采样方法马尔可夫链蒙特卡罗 ( MCMC )。MCMC 是一种用于从后验分布中抽取样本的方法。我们绘制5000样本和样本的直方图(用 MCMC 的说法叫做轨迹)。在贝叶斯统计中，参数θ被表示为随机变量，而不是单个值。

sample_size=30

def get_traces(sample_size):
    data=bernoulli.rvs(0.3,size=sample_size)
    model=create_model(data)
    model.seed()
    mc1 = MCMC(model)
    mc1.sample(iter=5000,burn=1000)
    return mc1,mc1.trace('uni_prior')[:]

mc1,traces=get_traces(sample_size)
plott.histogram(traces,"uni_prior")

前面代码的输出如下:

我们绘制了参数θ的后验分布。注意分布有相当多的方差，分布的峰值(用黑色垂直线表示)不对应0.3，这是参数θ的真值。分布峰的峰值在贝叶斯术语中称为点估计，如果我们想将参数θ表示为单个值，这类似于最佳估计。

在下面的代码片段中，我们绘制了越来越多样本的后验分布。对于初始模型，我们仅使用30样本。增加样本数量是否有助于改善参数的估计？

num_samples=[20,50,100,500,5000]
for i in num_samples:
    m,traces=get_traces(i)
    plott.histogram(traces,"num samples = "+str(i),datarange=(0,0.6))

前面代码的输出如下:

我们看到越来越多的样本使分布变得更清晰，这表明它对参数θ的估计越来越有信心。

我们使用 PyMC 的绘图功能来创建层次模型的图像。我们还需要安装一些 python 库，比如pydot、pygraphviz(也要安装base graphviz)。

我们使用这个简单的单行代码创建点格式图表示:

gdot=pymc.graph.dag(mc1)

查看文件也相当简单。我们只需将dot对象写入一个.png文件并查看它，如下代码所示:

from IPython.display import Image
gdot.write_png('thumbtack_graph.png')
Image(filename='thumbtack_graph.png')

前面代码的输出如下:

贝叶斯网络的贝叶斯估计

类似于单个参数估计的贝叶斯情况，使用贝叶斯框架需要我们指定所有数据实例和未知参数的联合分布。

对于我们试图估计的参数，如果我们决定使参数先验独立(这可能不适用于所有情况)，那么计算后验变得更容易(这类似于最大似然估计中的似然分解)。如果我们有一个包含两个节点()的网络，那么我们可以独立于上的后验来计算的后验，同样的可分解性可以推广到更大的网络。

贝叶斯估计的例子

在learn_cpd_Bayesian.ipynbIPython 笔记本中，我们将使用 PyMC 从求职面试网络中估计一个 CPD 的参数，如下代码所示:

from libpgm.graphskeleton import GraphSkeleton
from libpgm.nodedata import NodeData
from libpgm.discretebayesiannetwork import DiscreteBayesianNetwork
import pandas as pd

from pymc import *
import matplotlib.pyplot as plt
from pymc.Matplot import plot

nd = NodeData()
skel = GraphSkeleton()
jsonpath="job_interview.txt"
nd.load(jsonpath)
skel.load(jsonpath)
skel.toporder()
bn = DiscreteBayesianNetwork(skel, nd)

我们要学习的 CPD 是Interview CPD，如下代码所示:

pd.DataFrame(bn.Vdata['Interview']['cprob']).transpose()

前面代码的输出如下:

	 		  0	  1	  2
['0', '0']	 0.8	 0.18	 0.02
['0', '1']	 0.3	 0.60	 0.10
['1', '0']	 0.3	 0.40	 0.30
['1', '1']	 0.1	 0.20	 0.70

CPD 有 12 个唯一的概率。前面输出中的每一行对应于对其父项(Experience和Grades)的一个特定赋值，并且可以用多项式分布来表示。在这种特殊情况下，Interview取三个值，因此代表每行的多项式可以比作一个三面骰子。

由于父变量有四个唯一的赋值，我们需要学习四个Multinomial变量的参数。

最后，我们在这个模型上运行 MCMC 采样并观察后验，如下所示:

def create_model(data,num_draws):
    partial_dirich = pymc.Dirichlet(name="partial_dirich",theta=[1.0, 1.0, 1.0])
    full_dirich=pymc.CompletedDirichlet(name="full_dirich",D=partial_dirich)
    multi = pymc.Multinomial(
            name="multi",
            value=data,
            n=num_draws,
            p=full_dirich,
            observed=True)
    model=Model([partial_dirich,full_dirich,multi])
    return model

def run_mcmc(model,**kwargs):
    mcmc = pymc.MCMC(model)
    mcmc.sample(**kwargs)
    return mcmc

前面的片段包含一个 PyMC 模型，由树中的以下三个变量组成:

• Dirichlet:先验由这个分布表示，它允许我们有一个初始先验，包含一组三个概率值，加起来是 1。在 PyMC 中，Dirichlet随机变量仅存储 k - 1 概率(对于具有 k 结果的多项式),因为可以计算 kth 值，因为概率总和为 1。
• CompletedDirichlet:为了简单起见，我们将其定义为Dirichlet变量的子代，这样更容易表示所有 k 值(本例中 k = 3)。
• Multinomial:这是最后一个变量，是CompletedDirichlet对象的子对象。我们将从求职面试网络中抽取数据并将其传递给Multinomial变量。

注意，在图钉估计问题中，我们只学习了一个伯努利分布的估计。在这种情况下，我们必须学习四个单独的多项式分布的参数估计。

因此，我们对求职面试网络进行了四次抽样调查，每次都给父母分配一份工作。我们先从第一个家长作业开始: { Experience = 0，Grades = 0} ，帮助我们学习 CPD 第一行的参数。

PyMC 可以消费的样本需要特定的格式，由一系列实验的结果组成。

我们先来理解实验这个术语。在一个实验的单次运行中，从网络中进行采样 n 次(类似于 n 次掷出三面骰子)。这将为我们提供长度为 0、1 和 2 的序列 n 。我们修剪这个序列以获得 0、1 和 2 的出现频率(在get_counts方法中执行)，如以下代码所示:

def get_counts(vals):
    b = {0:0,1:0,2:0}
    for item in vals:
        b[item] = b.get(item, 0) + 1
    return np.array(b.values())
get_counts([0,1,2,0,1])

前面代码的输出如下:

array([2, 2, 1])

产生的数组告诉我们0和1出现了两次，2出现了一次。

我们进行了几个这样的实验(在这个特定的例子中是 10 个)，整理频率，并将其作为数据传递给我们的 PyMC 模型中的Multinomial变量。

然后，我们在 PyMC 模型上运行 MCMC 采样，就像在随机游走期间一样，PyMC 模型中full_dirich变量的值开始逼近真正的多项式分布，如以下代码所示:

def get_relevant_samples(experiments,evidence,num_samples=10):
    ''' for n experiments, sample and add the frequencies obtained.'''
    res=[]
    for i in xrange(experiments):
        vals=[float(i['Interview']) for i in bn.randomsample(num_samples,evidence)]
        res.append(get_counts(vals))
    return res

def plot_traces(traces):
    colors = ["#348ABD", "#A60628","#884732"]

    plt.plot(np.arange(len(traces)), traces[:,:, 0], c=colors[0])
    plt.plot(np.arange(len(traces)), traces[:,:, 1], c=colors[1])
    plt.plot(np.arange(len(traces)), traces[:,:, 2], c=colors[2])
    plt.title("traces of posterior ")
    plt.show()

def estimate_parameters(evidence,**kwargs):
    ''' run a few experiments, get the data, and estimate the parameter '''
    experiments=10
    samples=get_relevant_samples(experiments,evidence,num_samples=10)
    model=create_model(samples,experiments)    
    mc=run_mcmc(model,**kwargs)   

    traces=mc.trace('full_dirich')[:]
    return [np.mean(traces[:, :,0]),np.mean(traces[:,:, 1]),np.mean(traces[:,:, 2])],traces

我们首先确定父分配的参数， {Experience = 0，Grades = 0} ，如下代码所示，并绘制轨迹。我们看到三条独立的情节线。这些是多项式分布中三个概率的值，加起来应该是 1。

means,traces=estimate_parameters(dict(Grades='0',Experience='0'),iter=5000)
plot_traces(traces)

前面代码的输出如下:

在前面的剧情中，我们取前 5000 条痕迹。我们看到后部最初四处移动，然后稍微下沉。如果我们在采样前给 PyMC 添加一个burn参数，算法将丢弃由burn值指定的第一个 n 样本，并收集下一个轨迹。由于 MCMC 算法需要时间从先前的开始进行随机游走，指定一个burn值允许算法丢弃表面上不是来自真实后验分布的轨迹，如下面的代码所示:

means,traces=estimate_parameters(dict(Grades='0',Experience='0'),iter=20000,burn=1000)
plot_traces(traces)

前面代码的输出如下:

我们看到个体概率已经收敛。我们首先确定了父作业{经验= 0，等级= 0}的参数。为了学习父变量所有联合赋值的值，我们运行以下代码片段:

all_parent_assignments=[['0','0'],['0','1'],['1','0'],['1','1']]
res=[]
for i,j in all_parent_assignments:
    evidence=dict(Grades=i,Experience=j)
    means,traces=estimate_parameters(evidence,iter=50000,burn=10000)
    res.append([[i,j]]+means)

pd.DataFrame(res,columns=['parent_assignment',0,1,2])

前面代码的输出如下:

 par_assign	   0	          1	          2
0	 [0, 0]	 0.835853	 0.144636	 0.019511
1	 [0, 1]	 0.303323	 0.464297	 0.232380
2	 [1, 0]	 0.291015	 0.563842	 0.145143
3	 [1, 1]	 0.126254	 0.223568	 0.650178

我们将前面的后验均值与后面的工作面试网络中指定的真实概率进行比较。一些学到的价值观非常接近真实的价值观，而另一些则与基本事实相差高达 16%。

pd.DataFrame(bn.Vdata['Interview']['cprob']).transpose()

前面代码的输出如下:

par_assign    0	  1	  2
['0', '0']	 0.8	 0.18	 0.02
['0', '1']	 0.3	 0.60	 0.10
['1', '0']	 0.3	 0.40	 0.30
['1', '1']	 0.1	 0.20	 0.70

总结

在这一章中，我们学习了参数估计问题。我们首先使用两种方法估计图钉的参数，最大似然法和贝叶斯法。然后，我们扩展这些基础，使用最大似然估计来确定贝叶斯网络的参数。

已经完成了结构和参数学习的基础知识，在接下来的章节中，我们将使用网络来回答我们的查询，使用各种类型的推理方法。

六、使用图模型的精确推理

到目前为止，我们已经了解了构建图模型的方法。在本章中，我们将使用各种推理机的完全指定的图模型来获得问题的答案。

在我们开始推理之前，我们将了解推理问题的复杂性。在特定的环境中使用不同类型的推理是合适的，我们将学习哪种类型适用于哪里，然后用这两种方法进行实验。

推理的复杂性

图模型可用于回答概率查询和地图查询。使用该模型最简单的方法是生成联合分布，并将所有变量相加，除了我们感兴趣的变量。然而，我们需要确定和指定指数爆炸发生的联合分布。

在最坏的情况下，我们需要确定 NP-hard 中的确切推断。精确这个词的意思是以一定的精度指定概率值(比如，小数点后五位数)。假设我们降低精度要求(例如，小数点后最多两位数)。现在，(近似)推理任务变得容易了吗？不幸的是没有——即使是近似推理也是 NP 难的，也就是说，获得值远比随机猜测(50%或 0.5 的概率)好，后者需要指数时间。

看起来推理是一项无望的任务，但这只是最坏的情况。在一般情况下，我们可以使用精确推理来解决某些类别的现实世界问题(例如具有少量离散随机变量的贝叶斯网络)。当然，对于更大的问题，我们不得不求助于近似推理。

由于这本书是关于 Python 中的图模型的，让我们离题到我们可以用来运行推理的工具的选择。

现实世界的问题

由于推理是一个 NP 难的任务，所以推理机是用尽可能接近裸机的语言编写的；通常用 C 或 C++。由于这是一本关于 Python 中 PGM 的书，我们有几个选择:

• 使用推理算法的 Python 实现。这些完整和成熟的软件包并不常见。
• 使用有 Python 接口的推理引擎，比如斯坦(mc-stan.org)。这个选择很好地平衡了运行 Python 代码和快速推理实现之间的关系。
• 使用没有 Python 接口的推理引擎，这对于大多数推理引擎来说都是正确的。一份相当全面的名单可以在http://en.wikipedia.org/wiki/Bayesian_network#Software找到。这里 Python 的使用仅限于创建一个文件，该文件以推理引擎可以使用的格式描述模型。

在关于推理的章节中，我们将坚持列表中的前两个选择。我们将在运行玩具大小的问题时使用(推理算法的)本机 Python 实现来窥视推理算法的内部，然后使用带有 Python 接口的外部推理引擎来尝试一个更真实的问题。

使用变量消除算法

在部分，我们将学习变量消去算法。在VarElim_asia.ipynb IPython 笔记本中，我们将使用亚洲网络来了解该算法的细节。

亚洲网络(http://www.bnlearn.com/bnrepository/#asia)是一个用于患者诊断的玩具贝叶斯网络。该网络试图根据去亚洲旅行、吸烟史和 x 光结果等因素来引出肺癌或肺结核的可能性。有关亚洲网络的更多详情，请访问http://www.norsys.com/tutorials/netica/secA/tut_A1.htm。

以下为亚洲网示意图(由http://www.bnlearn.com/bnrepository/asia/asia.png提供):

这个代码片段的目标是比较两种推理方法。第一种方法是蛮力，我们首先列出整个联合分布，然后询问我们的概率。第二种方法使用变量消除，我们将研究它如何改进蛮力方法。

假设我们希望推断患者患有支气管炎的概率；从今以后，这将被写成。

我们先来了解一下我们是如何利用积法则得出的值的，如果我们知道，就可以计算。此外，可以用类似的方式计算，然后概率必须归一化。

我们将从存储网络定义的文件中加载网络，如下所示:

from libpgm.graphskeleton import GraphSkeleton
from libpgm.nodedata import NodeData
from libpgm.discretebayesiannetwork import DiscreteBayesianNetwork
from libpgm.tablecpdfactor import TableCPDFactor
import itertools
import pandas as pd
from libpgm.tablecpdfactorization import TableCPDFactorization

def loadbn(jsonpath):
    nd = NodeData()
    skel = GraphSkeleton()
    nd.load(jsonpath)
    skel.load(jsonpath)
    skel.toporder()

    bn = DiscreteBayesianNetwork(skel, nd)
    return bn

bn=loadbn("asia1.txt")

我们第一次尝试计算边际概率将使用所有变量的联合分布，然后将我们不需要的变量相加。为了列出联合分布，我们将所有因素相乘。根据因子乘积计算的联合分布列在以下代码中，其完整规范有 255 个条目:

#a method that prints the distribution as a table.
def printdist(jd,bn,normalize=False):
    x=[bn.Vdata[i]["vals"] for i in jd.scope]
    zipover=[i/sum(jd.vals) for i in jd.vals] if normalize else jd.vals
    #creates the cartesian product
    k=[a + [b] for a,b in zip([list(i) for i in itertools.product(*x[::-1])],zipover)]
    df=pd.DataFrame.from_records(k,columns=[i for i in reversed(jd.scope)]+['probability'])
    return df

#instantiate TableCPDs for all the nodes and multiply them. The result is in the factor that calls the methods.
jc=TableCPDFactor("asia",bn)    
[jc.multiplyfactor(TableCPDFactor(i,bn)) for i in bn.V if i != "asia"]
df=printdist(jc,bn,normalize=True)
print "values in joint distribution ",len(jc.vals)
#print the first few values in the table
df.head()
values in joint distribution  256

前面代码的输出如下:

Rows	 dysp	 xray	either	tub	lung	bronc	smoke	asia	probability
0	 yes	 yes	 yes	 yes	 yes	 yes	 yes	 yes	 0.000013
1	 yes	 yes	 yes	 yes	 yes	 yes	 yes	 no	 0.000262
2	 yes	 yes	 yes	 yes	 yes	 yes	 no	 yes	 0.000001
3	 yes	 yes	 yes	 yes	 yes	 yes	 no	 no	 0.000013
4	 yes	 yes	 yes	 yes	 yes	 no	 yes	 yes	 0.000007

上表列出了联合分布中 255 行中的前五行。

在整个讨论中，我们还会稍微增加符号。由于所有变量都是二元的，表示证据变量smoke为真，除非另有说明，即与相同

边缘化不相关的因素

从联合分布中，我们可以通过执行以下操作获得期望的边际概率。我们只对bronc栏的数据感兴趣。因此，我们首先折叠其他列(以及其中的数据)，这是一个称为因子边缘化的操作。

我们将使用Pandas数据分析库来操作类似于表格的 CPD。将联合分布加载到Pandas数据框后，我们将只对我们想要的列进行分组，并将probability列中的值相加，如下代码所示:

t2=df.groupby(['bronc','smoke'],as_index=False)
t3=t2['probability'].sum()
t3

前面代码的输出如下:

Rows	bronc	smoke	probability
0	 no	 no	 0.35
1	 no	 yes	 0.20
2	 yes	 no	 0.15
3	 yes	 yes	 0.30

这个操作是怎么发生的？让我们以上表中的第 0 行(bronc == no，smoke == no)为例，看看该值是如何获得的。0.35的值是将其他列的所有值相加得到的。我们将选择表格中所有bronc和smoke值都为no的行，如下代码所示:

df.loc[(df["smoke"] == 'no') & (df["bronc"] == 'no'), :].head()

前面代码的输出如下:

 dysp	xray	either	tub	lung	bronc	smoke	asia	probability
6	 yes	 yes	 yes	 yes	 yes	 no	 no	 yes	 0.000001
7	 yes	 yes	 yes	 yes	 yes	 no	 no	 no	 0.000024
14	 yes	 yes	 yes	 yes	 no	 no	 no	 yes	 0.000119
15	 yes	 yes	 yes	 yes	 no	 no	 no	 no	 0.002353
22	 yes	 yes	 yes	 no	 yes	 no	 no	 yes	 0.000023

我们将从上表的中选择所有行，只选择probability列。然后，我们将总结这些值，如下面的代码所示:

df.loc[(df["bronc"] == 'no') & (df["smoke"] == 'no'), "probability"].sum()

前面代码的输出如下:

0.34999999999999998

我们可以看到因子边缘化操作如何给我们赋值的概率0.35(bronc = = no，smoke ==no )。

过滤证据的因子约简

让我们稍微离题一下来观察查询观察到变量如的情况。然后，我们只需要前一个表的第一行和第三行(带有列bronc、smoke和probability)。这些行对应

该操作称为因子减少，因为其他行没有smoke的具体赋值(即值不是yes)，所以被删除，如下代码所示:

t4=t3.loc[ (t3["smoke"] == 'yes'), :]
t4

前面代码的输出如下:

Rows	bronc	smoke	probability
1	 no	 yes	 0.2
3	 yes	 yes	 0.3

这几乎是我们需要的概率集，只是概率加起来不等于 1。因此，我们将通过将每个概率除以所有概率的总和来进行归一化，如以下代码所示:

psum=t4['probability'].sum()
t4['probability']=t4['probability']/psum
t4

前面代码的输出如下:

 bronc	smoke	probability
1	 no	 yes	 0.4
3	 yes	 yes	 0.6

让我们回到最初的查询——。我们已经看到上表表示了smoke存在时bronc的条件概率分布。现在，我们简单地将smoke的值相加，得到的值，如下代码所示:

t5=t3.groupby(['bronc'],as_index=False)
t5['probability'].sum()

前面代码的输出如下:

 bronc	probability
0	 no	 0.55
1	 yes	 0.45

上表给出了推理查询的期望结果。

蛮力方法的缺点

我们之前讨论的是得出概率查询值的蛮力方法。这种方法效率低下，原因如下:

• 它不使用模型中存在的任何条件独立性。假设我们想知道的概率，以及节点 X 和 Y 是否有条件独立于 A 和 B，那么计算联合分布就没有意义了，联合分布包括 X 和 Y，然后求和出值。
• 我们已经知道，获取联合分布中的概率值并不是一件容易的事情，存储和操作大的联合分布也不是一件容易的事情。对于报警示例中的仅 8 个二进制变量，联合分布有 255 行(或 2 个 ^{8 个}，对于更大的网络，学习、存储和操作联合分布变得不可能。即使有了联合分布，将其他值相加也有一个时间复杂度。

使用变量消去法

我们现在来看看变量消去算法，这是一种计算概率的有效方法。效率源于避免重复操作以及避免计算条件独立的概率。即使存在任何证据/观察变量，该算法也可用于计算边际概率。

我们将使用因子，类似于我们之前看到的表格。因子积(是两组元素之间的笛卡儿积)因子边缘化****因子约简的运算是的演算，通过我们可以操纵因子。

每个因素都有一个范围。如果我们把一个因子看作一个表，那么作用域中的变量就是表的列名。asia.txt定义的八个 CPD 为初始因子。例如因子范围内有变量tub和asia，对应 CPD。类似地，因子对应于 CPD。有多个父节点的 CPD，如，转换为因子，该因子的范围内有所有的父节点和子节点。

我们在下面的列表中列出了asia网络中的因素，以及它们所源自的相应的 CPD。包含所有变量的联合分布将是列表中各因素的乘积。

• 
• 
• 
• 
• 
• 
• 
• 

我们将使用变量消去算法来推断的条件概率，它对应于一个在其范围内具有相同变量的因子:

在扫描前面的因素列表时，我们会看到没有一个因素在其范围内有。因此，为了创建，我们必须使用因子乘积运算，其中结果因子范围内的变量应该出现在被相乘的因子之一中。

以为例，这组因素在其组合范围内有三个变量:。因此，因子的乘积将给出一个因子，其中概率值是通过将因子中与赋值匹配的行相乘得到的。

顾名思义，在变量消去法中，我们一次消去一个变量。每个步骤包括以下操作:

• 倍增因子
• 边缘化一个变量(它在所有倍增因素的范围内)
• 产生新的因素

我们将通过以下方式得出因子。下面表的每一行都详细说明了算法的一个步骤，即因子积、淘汰变量、新因子创建，追求。

|

第一步

|

要素产品

|

消除变量

|

新因子已创建

|
| --- | --- | --- | --- |
| one | | asia | |
| Two | | smoke | |
| three | | tub | |
| four | | lung | |
| five | | either | |
| six | | dysp | |

在图中的八个节点中，每一步都有六个被淘汰，这就给我们留下了我们想要的因子。

我们将使用libpgm库来完成算法的每一步。libpgm库为我们提供了multiplyfactor和sumout方法来在每一步创建一个新的因子。我们已经介绍了以下代码的第一步:

asia=TableCPDFactor("asia",bn)
phi_1=TableCPDFactor("tub",bn)

phi_1.multiplyfactor(asia)
printdist(phi_1,bn)

前面代码的输出如下:

Rows	asia	tub	probability
0	 yes	 yes	 0.0005
1	 yes	 no	 0.0095
2	 no	 yes	 0.0099
3	 no	 no	 0.9801

由于起始因子只是 CPD，所以我们在创建TableCPDFactor("tub",bn)时，是同时涉及tub和asia的因子，这是从的 CPD 计算出来的。我们现在将使用以下代码消除asia:

phi_1.sumout("asia")
printdist(phi_1,bn)

前面代码的输出如下:

 tub	probability
0	 yes	 0.0104
1	 no	 0.9896

第二步，我们将因子相乘，消去smoke产生，如下代码所示:

phi_2=TableCPDFactor("smoke",bn)
[phi_2.multiplyfactor(TableCPDFactor(i,bn)) for i in ["lung","bronc"]]
phi_2.sumout("smoke")
printdist(phi_2,bn)

前面代码的输出如下:

 bronc	lung	probability
0	 yes	 yes	 0.0315
1	 yes	 no	 0.4185
2	 no	 yes	 0.0235
3	 no	 no	 0.5265

在第三步中，乘以因子并消除tub产生，如下代码所示:

phi_3=TableCPDFactor("either",bn)
phi_3.multiplyfactor(phi_1)
phi_3.sumout("tub")
printdist(phi_3,bn)

前面代码的输出如下:

Rows	lung	either	probability
0	 yes	 yes	 1.0000
1	 yes	 no	 0.0000
2	 no	 yes	 0.0104
3	 no	 no	 0.9896

在第四步中，乘以因子以消除lung并产生因子，如下代码所示我们没有打印 CPD，因为我们已经看到了它的样子。

phi_4=phi_3
phi_4.multiplyfactor(phi_2)
phi_4.sumout("lung")
print "variables in scope ",phi_4.scope

前面代码的输出如下:

variables in scope  ['either', 'bronc']

在第五步，将因子相乘至消除either，产生因子，如下代码所示:

phi_5=TableCPDFactor("xray",bn)
phi_5.multiplyfactor(phi_4)
phi_5.multiplyfactor(TableCPDFactor("dysp",bn))
phi_5.sumout("either")
print "variables in scope ",phi_5.scope

前面代码的输出如下:

variables in scope  ['xray', 'bronc', 'dysp']

在最后一步中，从因子中，我们将消除dysp以产生因子，如下代码所示:

phi_6=phi_5
phi_6.sumout("dysp")
printdist(phi_6,bn)

前面代码的输出如下:

Rows	bronc	xray	probability
0	 yes	 yes	 0.055843
1	 yes	 no	 0.394157
2	 no	 yes	 0.054447
3	 no	 no	 0.495553

前面的因子在其范围内有xray和bronc变量，这是我们需要的，由于我们只需要xray=yes的一个具体赋值，我们可以通过给定的证据来减少因子，如下代码所示:

phi_6.reducefactor("xray",'yes')
printdist(phi_6,bn)

前面代码的输出如下:

Rows	bronc	probability
0	 yes	 0.055843
1	 no	 0.054447

由于这是不是一个有效的概率分布，我们必须通过除以和来归一化概率，如下代码所示:

summ = sum(phi_6.vals)
phi_6.vals=[i/float(summ) for i in phi_6.vals]
printdist(phi_6,bn)

前面代码的输出如下:

 bronc	probability
0	 yes	 0.506326
1	 no	 0.493674

前面的代码片段详细描述了获得所需 CPD 的步骤。

使用libpgm库时，所有算法步骤都包含在condprobve方法中；因此，我们只需加载网络并使用该方法，如以下代码所示:

bn = loadbn("asia1.txt")
evidence = {"xray":'yes'}
query = {"bronc":'yes'}
fn = TableCPDFactorization(bn)
result = fn.condprobve(query, evidence)
printdist(result,bn)

前面代码的输出如下(注意，我们得到的值与逐步过程中得到的值相同):

Rows	bronc	probability
0	 yes	 0.506326
1	 no	 0.493674

变量消除算法可以总结如下:

• 每个节点上的 CPD 的起始因素
• 从因素中消除非查询变量 Z
• 将其余因素相乘
• 对所有非查询变量重复相同的步骤，直到只剩下查询变量(和证据/观察变量，如果有的话)

变量消除算法适用于贝叶斯网络和马尔可夫网络。一旦 Z 中的非查询变量集从所有因素的范围中消除，则算法完成。虽然 Z 中的变量可以以任何顺序消除(它产生相同的最终 CPD)，但在下一节中，我们将了解优化的消除顺序可以帮助算法快速终止。

变量消去的复杂性

在前一节的开始，我们声称对于推理查询，使用变量消除算法是对查询联合分布的改进。为了更好地理解为什么变量推理是一种改进，我们需要了解算法的算法复杂度。

我们将从 m 因子开始，在每个消除步骤中，我们将生成一个因子(通过消除一个非查询变量)。如果我们有 n 个变量，我们最多有 n 轮淘汰。生成的因子总数，即 m* ，将少于 m + n (我们从因子开始，除了淘汰的 n 变量)。

让 N 代表最大因子(其范围内变量数最多的因子)的大小。

算法中的每一步都包括推导出一个 Factor 乘积，然后求和出一个变量，这叫做和积运算。

因此，复杂度与和积运算的次数成正比。乘积运算小于N×m **，和运算为N×N*。因此，就 N 和 m* 而言，复杂度是线性的，最大因子的大小和生成因子的总数也是线性的。

虽然出现了线性项，但事实上，计算最大因子 N 需要指数时间。如果一个因子有四个变量，并且都是二进制值，那么它的复杂性就是。一般情况下，是计算一个因子的计算成本，如果 v 是一个变量在其范围内的最大值个数(称为其基数) k 是变量个数。

为什么消去顺序很重要？

尽管变量消除算法没有指定消除变量(非查询、非证据)的顺序，但消除顺序在复杂性中起着一定作用。我们来看看下图中的马尔可夫网络:

在前面的网络中，假设我们选择求和或消除变量 A 。我们首先需要有一个因子的乘积，然后将变量求和。创建的因子有一个范围，在 n 中是指数的。

相反，如果我们选择先将变量相加，因子乘积会产生一个范围为的因子，它只有三个变量。我们在前面的部分已经了解到复杂性包含术语 N (最大因子的大小)。假设所有变量都是二进制的，消去顺序的不同导致第一种情况下的和第二种情况下的的复杂性。

由于变量消除算法的复杂度很大程度上取决于生成的最大因子的大小(其范围是指数级的)，因此由消除排序来生成小的中间因子以改善变量消除算法的运行时间。本示例摘自 Coursera PGM 课程，可在https://www.coursera.org/course/pgm.访问

图透视

当我们忙于执行因子操作时，图结构也随着每次因子相乘和边缘化而改变。我们知道，因素和图只是同一信息的不同表示；那么，新因素的消除和创造是如何影响图的呢？

由于有向图和无向图在变量消除算法中的工作方式相同，我们可以通过假设图是无向图(即使对于贝叶斯网络)来继续分析。

我们将从上一节看一个亚洲网络中图变化的例子。当我们通过乘法来消除一个变量，也就是作用域中节点的父节点时，结果因子在一个叫做道德化的过程中增加了子节点之间的联系。

例如，在运行亚洲网络的变量消除算法的第二步中，我们将乘以因子并消除smoke以产生一个新的因子:。这个新因子表示两个节点之间增加了一个新的链接，如下图所示(左侧和右侧表示第二步完成前后的变量):

为什么要创建一个新的因子需要我们连接之前没有连接的节点？因为一个因子编码一些独立，同样的独立也必须存在于图中。因此，随着因子消除和(新的)因子创建的过程在变量消除算法中继续，我们向图中添加新的边来编码相同的独立性。由于道德化而产生的马尔可夫网络被称为诱导马尔可夫网络。

对于在变量消除算法中生成的每个因子，因子范围内的变量通过称为填充边的边连接(即，如果不存在边，则添加边)。对应于每个因子的全连通子图是该因子中变量分布的最小 I 图。你可以回想一下，如果满足以下条件，图 G 是分布 P 的最小 I 图:

• g 是 P 的一个 I 图
• 如果和不是 P 的 I 图

换句话说，一个最小的 I 图是一组独立的图，从 G 中移除任何边都会导致它不再是 I 图。

添加的边取决于变量消除的顺序(这也决定了创建的因素)。

从图结构中学习诱导宽度

在我们继续讨论诱导宽度之前，让我们离题提醒自己一些用来描述图结构的术语。团是一个极大的、完全连通的子图。让我们看看下面这个有四个节点的马尔可夫网络:

节点 B 、 C 和 D 由于都是相互连接的，所以形成了一个小团体。团是最大连接的，因为它不能添加更多的节点。将 A 添加到小团体将使完全连接的属性失效。

诱导宽度为什么重要？

诱导宽度是最大团的节点数减 1。最小诱导宽度是在所有 VE 排序上获得的最小诱导宽度，这将是最佳性能的下限。

事实证明，在价值工程算法运行过程中产生的每个新因素都是诱导马尔可夫网络中的一个小团体。因此，诱导图的小团体给了我们一个 ve 算法运行时间的快速近似。即使我们确实找到了最佳的 VE 排序(这本身就是一个 NP 难问题)，推理仍然需要指数时间，即使有最佳排序。因此，如果我们使用最优排序，并发现诱导图中的团在其范围内有许多(这是基于您的硬件的相对数量)变量，那么可能是时候抛弃精确推理方法，转而使用近似方法了。

寻找价值工程订单

贪婪算法是找到最佳 ve 排序的相当有效的机制。可以使用几个成本函数，例如首先选择最小的因子(邻居数量最少的节点)。

树算法

我们现在来看另一类基于消息传递的精确推理算法。

消息传递 是一种通用机制，消息传递算法存在多种变体。我们将看一小段团树消息传递算法(有时也称为连接树算法)。消息传递算法的其他版本也被用于近似推理。

我们通过澄清使用的一些术语开始讨论。

聚类图是一种网络排列，其中变量组被放置在聚类中。它类似于一个因子，其中每个聚类在其范围内有一组变量。

消息传递算法就是在集群之间传递消息。举个例子，想象一下聚会上的流言蜚语，雪莉和克莱尔正在谈话。如果 Shelly 认识 B、C 和 D，并且她正在和认识 D、E 和 F 的 Clair 聊天(注意，他们唯一共同认识的人是 D)，他们可以分享关于他们共同朋友 D 的信息(或传递消息)。

在消息传递算法中，两个集群通过一个分离集 ( 分离集)连接，该分离集包含两个集群共有的变量。使用前面的例子，两个集群和通过分离集连接，分离集包含两个集群共有的唯一变量。

在下一节中，我们将了解连接树算法的实现细节。我们将首先了解算法的四个阶段，然后使用代码片段从实现的角度了解它。

连接树算法的四个阶段

在本节中，我们将讨论连接树算法的四个阶段。

在第一个阶段，贝叶斯网络被转换成称为连接树的二级结构(文献中这种结构的替代名称是连接树、聚类树或团树)。从贝叶斯网络到连接树的转换按照以下步骤进行:

• 我们将通过把所有有向边变为无向边来构造一个道德图。进入所述节点的具有 V 型结构的所有节点的父节点都与一条边相连。我们已经看到了这个过程的一个例子(在 ve 算法中)叫做道德化，这是一个可能的参考来连接(显然未婚)有一个孩子(节点)的父母。
• 然后，我们将有选择地为道德图添加边，以创建一个三角图。三角图是一个无向图，其中节点之间的最大循环长度是 3。
• 从三角化的图中，我们将识别节点的子集(称为团)。
• 从团作为簇开始，我们将排列簇以形成一个称为连接树的无向树，它满足运行交集属性。该属性声明，如果一个节点出现在两个群中，它也应该出现在连接这两个群的路径上的所有节点中。

在第二阶段，初始化每个簇的电势。电位类似于 CPD 或表格。它们有一个值列表，对应于其范围内变量的每个赋值。团簇和分离集都包含一组势。势这一术语与概率相对使用，因为在马尔可夫网络中，与概率不同，势的值不必加起来等于 1。

这个阶段包括相邻集群之间的消息传递或信任传播。每条消息都包含一个信念，即集群拥有一个特定的变量。

每个消息都可以异步传递，但是它必须等待来自其他集群的信息，然后整理这些信息并将其传递给下一个集群。考虑一个树形结构的集群图会很有用，其中消息传递分为两个阶段:向上传递阶段和向下传递阶段。只有在节点从叶节点接收到消息后，它才会将消息发送给其父节点(在“向上传递”中)，只有在节点从其父节点接收到消息后，它才会将消息发送给其子节点(在“向下传递”中)。

当每个集群分离集具有一致的信念时，消息传递阶段就完成了。回想一下，连接到 sepset 的集群有公共变量。例如，聚类 C 和 sepset S 在其范围内有和变量。然后，从聚类或分离集获得的对的势具有相同的值，这就是为什么说聚类图具有一致的信念或集团被校准。

一旦整个聚类图具有一致的信念，第四个阶段是边缘化，在这里我们可以查询图中任何变量的边际分布。

我们现在将继续研究连接树算法的实现。

使用连接树算法进行推理

在JunctionTreeAlgorithm.ipynb IPython 笔记本中，我们将使用贝叶斯信念网络 ( BBN )库，使用连接树算法运行精确推理。在 Github(https://github.com/eBay/bayesian-belief-networks)上有库，在 Github 页面上提到了安装库的文档。

BBN 具有加载存储在贝叶斯交换格式 ( bif 中的网络的功能，该格式由贝叶斯社区开发，以促进不同推理工具之间更容易的数据共享。

我们将再次使用我们在本章前面看到的asia网络。

强制导入后，我们用bif_parser模块解析.bif格式文件，返回一个贝叶斯网络对象，如下代码所示:

import bif_parser
import prettytable
import pydot
from IPython.core.display import Image 
from bayesian.bbn import *

name = 'asia'

module_name = bif_parser.parse(name)
module = __import__(module_name)
bg = module.create_bbn()

我们可以使用 BBN 提供的graphviz功能查看贝叶斯网络(graphviz是一个图可视化工具)，如下面的代码所示:

def show_graphgiz_image(graphviz_data):
    graph = pydot.graph_from_dot_data(graphviz_data)
    graph.write_png('temp.png')
    return 'temp.png'

sf=bg.get_graphviz_source()
Image(filename=show_graphgiz_image(sf))

上图为我们展示了.bif文件中编码的网络结构。就是我们在本章前面看到的同一个asia网络。

阶段 1.1–道德化

在下面的片段中，我们将查看道德化阶段。请注意，V 型结构，例如，让他们的父母说教或加入了一个新的环节。

gu=make_undirected_copy(bg)
m1=make_moralized_copy(gu,bg)
s2=m1.get_graphviz_source()
Image(filename=show_graphgiz_image(s2)) 

阶段 1.2–三角测量

在三角化阶段，无向图是三角化的，如果有一个周期长度大于 4，则是节点间边的相加。注意，在上图中，形成了一个有四个节点的循环。在下面的片段中，在和之间添加了一个链接来对其进行三角测量，这将循环的最大长度减少到 3:

cliques, elimination_ordering = triangulate(m1, priority_func)
s2=m1.get_graphviz_source()
Image(filename=show_graphgiz_image(s2))

前面代码的输出如下:

阶段 1.3–构建连接树

接下来，我们将创建团和分离集，这完全在build_join_tree方法中完成。该方法根据前面的图创建小团体，并创建每对小团体之间的交叉点 sepsets。

请注意，生成的图中节点的命名约定并不十分美观。Clique_EITHERLUNGTUB团有变数either、lung、tub都砸在一起了。

jt=bg.build_join_tree()
sf=jt.get_graphviz_source()

Image(filename=show_graphgiz_image(sf)) 

我们可以从变量名中看出，分离集包含了它们所连接的小集团的变量的交点。例如，派系和中的变量就是分离集中的变量。

阶段 2–初始化电位

下一步是创建初始集群(通常是小集团)和初始化电势，如下代码所示:

assignments = jt.assign_clusters(bg)
jt.initialize_potentials(assignments,bg)

第 3 阶段——信息传递

在建立了运行消息传递所需的结构后，我们进入下一阶段。在这个阶段，每对拉帮结派之间会发送两条消息:一条是正向传递，另一条是反向传递。

消息传递的细节包含在连接树对象中的propagate()方法中。下图显示了网络中的消息序列:

Image(filename="../book/chapteimg/9004OS_06_05.png")

那么，当一条消息被传递时，到底会发生什么呢？一个消息传递有三个参与者:源集群、中间的分离集和目标集群。每一个都有一组势(非常类似于 CPD 以及相关的概率，除了它不需要是有效的概率分布)。

例如，从到的消息将传递分离集，其中变量在每个集群或分离集的范围内。

消息将修改和的电位，即分离电位和目的簇电位。

对分离集 R 的第一个赋值是源中的电位，分离集中没有的变量被边缘化，如下式所示:

对集群 Y 的第二次分配如下:

让我们使用propagate()方法运行消息传递位，如下面的代码所示。这将运行所有变量之间的消息传递，并在信念收敛时停止。

jt.propagate()

一旦所有的信息传递完成，我们就剩下一棵树，它的所有集群都有一致的信念。

在查询特定变量(例如bronc)时，我们只需要找到一个范围内有bronc的聚类(或一个分离集)，并边缘化其他变量。这里打印的集群范围内有bronc、dysp和either变量。我们打印了所有与这个集群相关的潜力。我们可以观察到每一行都列出了对bronc、dysp和either的具体分配。

bronc_clust=[i for i in jt.clique_nodes for v in i.variable_names if v =='bronc']
bronc_clust[0].potential_tt

前面代码的输出如下:

{(('bronc', 'no'), ('dysp', 'no'), ('either', 'no')): 0.4689219599,
 (('bronc', 'no'), ('dysp', 'no'), ('either', 'yes')): 0.008692680,
 (('bronc', 'no'), ('dysp', 'yes'), ('either', 'no')): 0.05210244,
 (('bronc', 'no'), ('dysp', 'yes'), ('either', 'yes')): 0.020282,
 (('bronc', 'yes'), ('dysp', 'no'), ('either', 'no')): 0.08282951999,
 (('bronc', 'yes'), ('dysp', 'no'), ('either', 'yes')): 0.003585240,
 (('bronc', 'yes'), ('dysp', 'yes'), ('either', 'no')): 0.3313180799,
 (('bronc', 'yes'), ('dysp', 'yes'), ('either', 'yes')): 0.032267160}

让我们通过使用以下代码边缘化集群中的dysp和either变量来尝试找到bronc的边际:

pot=bronc_clust[0].potential_tt

#a function to return the sum for a specific assignment, such as 'bronc,yes'
sum_assignments=lambda imap,tup:sum([v for k,v in imap.iteritems() for i in k if i == tup])

#get the sum for bronc=yes and bronc=no
yes,no=[sum_assignments(pot,('bronc',i)) for i in ['yes','no']]

print 'bronc: yes ', yes/float(yes+no)," no ", no/float(yes+no)

bronc: yes  0.45  no  0.55

既然我们声称连接树是一致的，那么包含bronc的其他聚类会返回相同的边际值吗？让我们使用第二个范围为bronc的集群，并忽略它的潜力，如下面的代码所示:

pot2=bronc_clust[1].potential_tt
yes,no=[sum_assignments(pot2,('bronc',i)) for i in ['yes','no']]

print 'bronc: yes ', yes/float(yes+no)," no ", no/float(yes+no)
bronc: yes  0.45  no  0.55

我们可以看到，我们可以从两个不同的聚类中获得相同的bronc的边际值，这表明所有变量的信念在所有聚类和分离集中是一致的。

注

关于连接树算法的详尽介绍，请参考黄和达维切在的程序指南。

在总结中，本章中给出的两种算法在原理上是相似的，它们使用相同的因子乘积和因子边缘化运算。团树中的小团体类似于变量消除中使用的因子。

连接树算法有一些优点，例如能够在校准树的单次计算中回答几个边缘查询。由于消息传递是异步的，所以算法可以并行化，并且可以更快地返回结果。

该算法的缺点是在空间方面更昂贵。在变量消去法中，中间因子不被存储，但是它们被存储在连接树中。

总结

我们首先探讨了推理问题，在这里我们研究了推理的类型。然后，我们了解到推理是 NP 难的，并了解到，对于大型网络，精确的推理是不可行的。

然后，我们探索了一些精确的推理算法，例如变量消除和消息传递，以及解释其基本原理的代码示例。我们研究了变量消除的复杂性和降低其复杂性的方法。

在下一章中，我们将研究运行近似推理算法的方法。

七、近似推理方法

在前一章中，我们已经了解到，随着图模型的树宽增加，精确的推断变得不可行。追求近似推理的动机来自真实世界的网络，在那里精确的推理是难以处理的。

在本章中，我们将学习计算近似推理的方法。我们将重温前一章中提到的消息传递算法，并了解当网络不是树形结构时，它们如何也可以用于计算近似推理。我们将从优化的角度探索推理，并学习采样方法。我们还将查看一些代码示例，以了解算法如何实现近似推理。

近似推理方法有两大类。在第一部分中，我们将探索作为优化的推理，在第二部分中，我们将探索基于粒子的推理方法。

优化视角

在近似推理中，我们寻求找到或构建目标分布的近似。假设我们有真实的目标分布，我们寻求找到一组容易运行推理的分布 Q，然后在这些容易的分布中搜索“接近”目标分布的那个。所使用的方法寻求优化用于测量 Q 和 p 之间的距离或相似性的函数

将推理问题归结为一个优化问题，可以让我们借鉴约束优化领域中研究得很好的方法。最常用的方法之一是生成一组描述目标函数最优值的方程。在图模型的上下文中，这种方法采用一组等式的形式，其中每个变量都是根据其他变量定义的。事实证明，将问题转化为一种约束优化的方程相当于图中的消息传递。

在这一节中，我们将探索信念传播算法。在前一章中，我们在树形结构的网络上运行了信念传播，而在这里，我们希望将讨论扩展到具有其他类型的图结构的网络。这个算法被称为 Loopy 信念传播 ( LBP )之所以这样命名，是因为聚类图可以包含循环，有时被称为 Loopy。

一般图中的信念传播

信念传播是用于执行推理的消息传递算法之一。尽管它最初被设计为在树形结构的网络上运行，但人们发现它也可以用于具有循环或循环的一般图。

在图结构上运行算法与在树结构上运行算法有一些不同。

|

树形网络

|

图结构网络

|
| --- | --- |
| 一个节点有首先从它所有的叶子节点接收消息，然后只向它的父节点发送一条消息，整个过程向下重复。 | 虽然树形结构网络中的节点可以在向其父节点发送消息之前等待其子节点的消息，但是在有环路的网络中，节点在发送消息之前等待消息并没有简单的过程。 |
| 每个节点和它的父节点之间只需要传递两条消息，就可以实现信念的收敛。 | 消息传递一直持续到达到收敛。 |
| 消息传递导致集群和分离集中信念的收敛。 | 大多数情况下会发生收敛，但这并不能保证。有些网络可能会出现振荡，有些可能永远不会收敛。 |
| 估计边缘收敛到真实边缘。 | 即使在收敛时，边缘也可能是不正确的，或者集群信念不一定等于真正的边缘。 |

创建运行 LBP 的集群图

运行 LBP 的第一个任务是当给定一组因子时，创建一个聚类图。每个因子必须分配给一个聚类，但是一个聚类可以分配多个因子。从一组因子到聚类图的转换不是唯一的。但是，对于执行信念传播的集群图，它必须满足以下属性:

• 运行交集属性:如果一个变量 X 存在于两个集群中，那么它也应该是存在于(所有集群)连接两个集群的唯一路径中。
• 家族保存:对于分配给一个聚类的每个因子，该聚类的范围必须包括该因子范围内的所有变量。

让我们看一个集群图创建的例子。

我们有如下公式所示的一组因子(以及每个因子范围内的变量):

下图显示了一种可能的集群图表示。分离集中的变量在连接每个聚类的边上标注。

注意前面提到的两个属性在集群图中是如何满足的。

• 运行交叉点:集群 5 、 2 和 3 在其范围内包含变量 D 。运行交叉口要求，如果 5 和 3 在其范围内包含 D ，它们之间的簇也应该包含 D (簇 2 )。

• 家族保存:聚类因子赋值如下表所示(注意聚类的范围包括其赋值因子范围内的所有变量):

  |

  簇

  |

  工厂

  |
  | --- | --- |
  | 1:甲、乙、丙 | |
  | 2:C、D | |
  | 3 😄、F、H | |
  | 4:英、法 | |
  | 5:乙、戊、丁 | |

LBP 中的消息传递

在传递任何消息之前，每个集群的初始信念只是分配给它的所有因素的产物。

消息传递以类似于我们在上一章中看到的方式进行。首先，所有消息都被赋值 1。我们将把表示为从源集群 X 传递到目的集群 Y 的消息。我们将详细检查在 C1 和 C2 集群之间传递的和消息以及在该上下文中发送的其他消息。

• :这些是关于从 C5 发送到 C1 的变量 B 的信念。
• :这些是关于从 C1 发往 C2 的变量 C 的信念。这是通过执行以下步骤来完成的:

  • 以 C1 现有的信仰为例

  • 将传入消息相乘

  • Summing out variables A and B and calculating the message as follows:

    

  • 发送到 C2

• 我们将快进到几个时间段后，在此期间发送了以下消息:
  • :这些是关于从 C5 发送到 C2 的变量 D 的信念
  • :这些是关于从 C3 发往 C2 的变量 D 的信念
• 我们现在准备从 C2 向 C1 发送一条返回消息。回想一下，早些时候，集群 1 将从其网络部分获知的信念(关于公共变量 C)发送给集群 2 ，后者现在以其关于变量 C 的信念做出响应。它以以下方式创建消息:

  • 以 C2 现有的信仰为例

  • 将传入消息和相乘

  • Summing out variable D and calculating the message by using the following formula:

    

  • 发送到 C1

请注意，当消息从集群 2 返回到 1 时，它避免发送集群 1 首先发送的信息，否则从集群 1 发送的信念只是来回呼应(重复计数)，并且在这个过程中变得更强。

通过的消息以前述相同的方式继续，并在收敛时停止。当跨时间步长的分离集信念变化在预定的容差值内时，可以假设聚类图已经收敛。此外，每对集群信念与分离集信念一致(对于位于这对集群之间的边缘)。

LBP 算法中的步骤

我们可以将在 LBP 算法中的步骤总结如下:

• 将因子分配给集群

• 构建初始势(将一个簇中的所有因子相乘)

• 初始化每个集群的信念

• 重复消息传递步骤:选择一个集群对，并在它们之间传递消息

• 在 n 次重复后，测试收敛是否已经发生或停止

• Summarize beliefs at each cluster by multiplying all the messages received with their initial beliefs as follows:

  

这里，第 I 个集群中的信念是初始信念的产物，也是来自所有 k 相邻集群的所有消息的产物。

提高 LBP 的收敛性

LBP 算法是消息传递算法通用框架中的算法之一。LBP 的变体试图通过调整以下内容来优化性能:

• 消息传递的持续时间。
• 消息的定时(同步/异步):与同步消息传递相比，异步消息传递具有更快的收敛速度。必须理解，收敛是一种局部性质。一些集群/分离集比其他集群/分离集收敛得更快。如果我们对给定集合或聚类的边缘感兴趣，我们不需要等待整个聚类图收敛。
• 发送消息的路径:像温赖特的 TRWBP(树加权的)和科尔莫戈罗夫的 TRW-S(顺序树加权的 BP)这样的 BP 变体试图在一个聚类图中找到一组最小生成树，对于每个生成树，只在这些树上传递消息。这些在某些条件下有收敛保证，这是对原始 LBP 的改进，原始 LBP 不提供任何收敛保证。
• 消息的平滑或衰减:这防止了振荡并增加了收敛的机会。

应用 LBP 分割图像

我们来看看LBP_image_segmentation.ipynb IPython 笔记本中应用 LBP 进行图像分割的情况。

图像分割是将一幅图像分割成多个片段(称为超像素)的过程。图像分割算法的结果将是为图像中的每个像素分配一个类别标签。图像分割在几个领域中是有用的，例如物体检测(行人检测)、医学成像(定位肿瘤)以及其他几个领域。

在本笔记本中，我们将使用循环信念传播的 OpenGM 库(https://github.com/opengm/opengm)实现来分割图像。

OpenGM 是一个带有 Python 包装器的 C++库，有几种推理算法的实现。安装说明请参考 OpenGM 网站。对于 Python 包装器，OpenGM 还需要安装 HDF5(http://www.hdfgroup.org/HDF5/)库。

我们将首先加载图像，并使用以下代码将其转换为灰度图像:

import opengm
from matplotlib import pyplot as plt
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import Image

fname = 'cow_image.jpg'
image = Image.open(fname).convert("L")
arr = np.asarray(image).astype(float)/255
plt.imshow(arr, cmap = cm.Greys_r)
plt.show()

下面是我们试图分割的原始图像。一个成功的分割算法应该能够区分分配给奶牛的像素和分配给背景的像素。

图像存储在numpy阵列中，其尺寸为 183 x 275 像素。每个像素包含一个介于 0 和 1 之间的值，其中 0 代表黑色，1 代表白色。

我们将离题一点，以了解 OpenGM 使用的基于能源的模型。

了解基于能源的模型

我们在图像中有 183×275 个像素，每个像素可以取两个标签中的一个。给像素指定一个标签叫做标记。我们希望从该图像的标签总数中选择正确的分割。给定图像的大小，可能的标签数量为。

给定观察到的图像特征，标记的概率是。我们希望使用以下公式找到最佳标注，即 MAP 估计值:

特定标记的概率定义如下:

这里，是配置的能量。

在这个讨论中，我们将使用以下几种势:

称为一元势的第一个势是一个在其范围内只有一个像素的因子。每个像素都有自己的一元电位。

第二个势是成对势，这是一个在其范围内具有相邻像素的因素。我们将在后面看到的图表将阐明这些潜力。

假设我们使用的是一元和成对的势，当给定一个特定的标记时，整个图像中的势之和代表图像的能量，如下式所示:

这里，右手边的第一项和第二项表示一元势和成对势。表示范围内分别有 1 个和 2 个像素的因子。

因此，给定给像素的特定标签分配，最佳标签或配置是具有最小能量的标签或配置。图像分割的目标是找到这种最佳配置。

关于图像分割用 MRFs 的详细介绍，请参考http://www . INF . u-szeged . Hu/~ kato/teaching/EMM/multi-layer-MRF . pdf。

在 3×3 的网格上可视化一元和成对的因子

在您尝试图像分割之前，让我们尝试使用 OpenGM 的模型可视化来理解马尔可夫网络的结构。

在下面的片段中，我们将创建一个马尔可夫网络，该网络有九个节点，排列在一个 3×3 的网格中。

下图是马尔可夫网络，其中从 0 到 8 的随机变量用白色未填充的圆圈表示，一元和成对的势用黑色方块表示:

grid2d2Order方法创建以下因素集:

• 第一组因素包括与每个变量相关联的一元电位，这些电位用从 0 到 8 的黑色方块表示，并附在每个节点上。一元因子的作用域只包含当前节点。

• 第二组因素包括连接相邻节点的成对电位。这些也用从 9 到 20 的黑色方块表示。成对因子的范围包含两个相邻的节点。比如图中的因子 13 可以用来表示，表示 2 和 5 是其范围内的节点。这个因素模拟了网格中两个相邻节点之间的交互。

  numLabels=2
  shape=(3,3)
  unaries=numpy.random.rand(shape[0],shape[1],numLabels)
  potts=opengm.PottsFunction([numLabels,numLabels],0.0,0.5)
  gm=opengm.grid2d2Order(unaries=unaries,regularizer=potts)
  opengm.visualizeGm( gm,plotFunctions=False,
                      plotNonShared=True,relNodeSize=0.9)
  

为了计算整个网络的能量，我们需要评估每个因素的能量。在 OpenGM 中，这是通过将一个函数与每个因子相关联来实现的。我们在前面的代码片段中使用了potts函数，如果输入变量(在因子的范围内)一致，则取值 0，如果不一致，则取值 0.5。

网络有21因子，边上的节点(0)有三个与之相关的因子，它的一元势(0)和它的成对势(9，10)。类似地，除了一元电势(4)之外，中间的节点(4)将具有四个成对电势(11，15，16，17)。

print "number of factors:",gm.numberOfFactors
print "number of factors of node 0: ",gm.numberOfFactorsOfVariable(0)
print "number of factors of node 4: ",gm.numberOfFactorsOfVariable(4)

前面代码的输出如下:

number of factors: 21
number of factors of node 0:  3
number of factors of node 4:  5

创建图像分割模型

我们将以与前面讨论的 3×3 网络相同的方式创建一个具有一元和成对势的马尔可夫网络。我们现在将有一个网格，而不是 3×3 的网格，该网格的大小=长度×宽度以像素为单位。

shape=img.shape
dimx,dimy=shape[0],shape[1]
numVar=dimx*dimy
numLabels=2
beta=0.1

numberOfStates=numpy.ones(numVar,dtype=opengm.index_type)*numLabels
gm=opengm.graphicalModel(numberOfStates)

#add Unary factors, assign Potts function to the factors
for y in range(dimy):
   for x in range(dimx):
      f=numpy.ones(2,dtype=numpy.float32)
      f[0]=img[x,y]
      f[1]=1.0-img[x,y]
      fid=gm.addFunction(f)
      gm.addFactor(fid,(x*dimy+y,))

#Adding binary function and factors

#create the pairwise function (Potts function)
f=numpy.ones(pow(numLabels,2),dtype=numpy.float32).reshape(numLabels,numLabels)*beta
for l in range(numLabels):
   f[l,l]=0
fid=gm.addFunction(f)

#create pairwise factors for the whole grid, and
#assign the Potts function created above, to each new factor.
for y in range(dimy):
   for x in range(dimx):
      #add a factor between each pair of neighboring nodes.
      if(x+1<dimx):
         #add a factor between a node and its neighbor on the right
         gm.addFactor(fid,numpy.array([x*dimy+y,(x+1)*dimy+y],dtype=opengm.index_type))
      if(y+1<dimy):
         #add a factor between a node and its neighbor above.
         gm.addFactor(fid,[x*dimy+y,x*dimy+(y+1)])

前面的代码片段执行了以下操作:

• 创建一个在其标签空间中具有图像高度x、图像宽度x和标签数量的图模型
• 创建一元因子，每个像素一个
• 给每个因子分配一个potts函数
• 创建成对因子，在 x 轴上相邻的每对像素一个因子，在 y 轴上相似
  • 给每个成对因子分配一个potts函数

我们几乎已经准备好开始在马尔可夫网络上运行推理了。OpenGM 允许我们创建一个可以观察推理过程的访问者类。在间隔调用的visit方法中，我们可以观察到正在进行的能量最小化。我们还可以观察到，随着马尔可夫网络的整体能量降低，每个像素的标签分配如何改善。

然后我们创建一个BeliefPropagation类的实例，并对其进行推理。LBP 算法将在像素之间传递消息，并在每个节点计算信念。调用回调方法时，我们使用以下代码查看当前配置的能量:

imgplot=[]

class PyCallback(object):
    def appendLabelVector(self,labelVector):
        #save the labels at each iteration, to examine later.
        labelVector=labelVector.reshape(self.shape)
        imgplot.append([labelVector])
    def __init__(self,shape,numLabels):
        self.shape=shape
        self.numLabels=numLabels
        matplotlib.interactive(True)
    def checkEnergy(self,inference):
        gm=inference.gm()
        #the arg method returns the (class) labeling at each pixel.
        labelVector=inference.arg()
        #evaluate the energy of the graph given the current labeling.
        print "energy  ",gm.evaluate(labelVector)
        self.appendLabelVector(labelVector)
    def begin(self,inference):
        print "beginning of inference"
        self.checkEnergy(inference)
    def end(self,inference):
        print "end of inference"
    def visit(self,inference):
        self.checkEnergy(inference)

inf=opengm.inference.BeliefPropagation(gm,parameter=opengm.InfParam(damping=0.05))
#parameter=opengm.InfParam(damping=0.1)
callback=PyCallback(shape,numLabels)
visitor=inf.pythonVisitor(callback,visitNth=1)

inf.infer(visitor)

前面代码的输出如下:

beginning of inference
energy   21002.9691286
energy   16100.9689678
energy   16082.888577
energy   16069.1768137
energy   16051.5650505
energy   16032.2630915
<rows elided>
..

energy   15824.9003709
end of inference

我们在每个中间步骤保存了标签，现在我们将查看前六个类别标签向量，以观察算法如何进行图像分割。我们可以看到，奶牛的形状很容易通过使用以下代码分配给每个像素的标签来识别:

fig = plt.figure(figsize=(16, 12))
for (counter, im) in enumerate(imgplot[0:6]):
    a=fig.add_subplot(3,2,counter+1)
    plt.imshow(im[0],cmap=cm.gray, interpolation="nearest")

plt.draw()

前面代码的输出如下:

最后，我们将绘制每个像素的最终标签分配，这是类标签对像素的 最大后验 ( MAP )分配，以及原始图像。这是对应于最小能量的配置。

fig = plt.figure(figsize=(16, 12))
a=fig.add_subplot(1,2,1)
plt.imshow(imgplot[-1][0],cmap=cm.gray, interpolation="nearest")
a=fig.add_subplot(1,2,2)
plt.imshow(img,cmap=cm.gray, interpolation="nearest")
plt.draw()

下面的图像是从前面的代码片段生成的:

枸杞多糖的应用

在 20 世纪 90 年代，高性能代码的新方法被发现，这些方法接近香农极限(对于特定的噪声水平，可靠通信可能达到的信道的理论最大承载能力)。原来，这种被称为 turbo 码(http://en.wikipedia.org/wiki/Turbo_code)的新代码正在实现 LBP，以实现这种改进的性能，这导致了人们对 LBP 及其变体在近似推理领域的兴趣激增。因此，LBP 现在被用于许多使用高性能代码的领域，例如移动电话标准、数字视频广播、移动电视和深空通信。

总之，消息传递是一种广泛使用的近似推理方法。对于连接不紧密的网络，推理相当有效。虽然像 LBP 这样的算法理论上不能保证收敛，但其经验性能往往相当不错。

基于采样的方法

我们将现在继续检查执行近似推理的另一种方法。

在基于抽样的方法中，我们使用从分布中抽取的样本来估计总体分布的统计量。抽取的样本是独立且同分布的。

在关于参数估计的一章中，我们从后验分布中抽取样本，使用最大似然和贝叶斯方法等方法来估计 CPD 的概率。

在本节中，我们将学习从贝叶斯网络中采样，这与采样分布略有不同。

正向采样

使用贝叶斯网络的拓扑排序的采样方法称为正向采样。贝叶斯网络中的拓扑排序是节点以形式的排序，这样对于每个边，我们都有。如果我们按照拓扑顺序遍历图，所有的边都指向前方。换句话说，我们先对父母进行采样，然后下降到叶子。

让我们以工作面试网络为例来看看抽样程序:

请注意，在上图中，采样的实体以粗体显示。

下面的表显示了一个取样程序的例子。表格的每一行列出:

• 我们从中取样的节点

• 父代的采样值(如果有)

• 采样的值为当前节点

  |

  数字

  |

  接触势差(接触电位差的缩写)

  |

  来自父母样本的证据

  |

  样本值

  |
  | --- | --- | --- | --- |
  | 一 | 经验 | - | |
  | 二 | 等级 | - | |
  | 三 | 采访 | 经验= 0,成绩= 0, | |
  | 四 | 提供 | 面试= 1 | |

因此，我们使用一次向前采样得到一个样本 。

要运行推理查询，如，我们可以使用样本使用最大似然估计边际概率。这可以通过用满足赋值的样本分数除以配分函数来实现。

接受-拒绝抽样方法

接受-拒绝抽样方法改进了前向抽样方法，加入了证据变量。

当我们引入证据变量时，我们这样修改抽样程序:我们首先抽取一个样本，测试它是否满足对证据变量的赋值，如果不满足则拒绝它(如果满足则接受它)。

我们知道，为了在一定的误差范围内推断边际概率，我们必须抽取大量样本。我们可以看到这种方法是浪费的，因为它抽取样本来拒绝(可能是大量的)样本。因此，积累大量满足证据的样本是一个耗时的过程。

让我们举一个癌症检测的例子。假设我们希望创建一个乳腺癌患者的样本。乳腺癌的患病率约为 125/10 万患者(http://seer.cancer.gov/statfacts/html/breast.html)。因此，如果我们抽取 100000 个样本，我们必须拒绝 99.875%的样本，因为它们不符合观察到的证据。如果我们有额外的观察变量，如 20 至 30 岁的年龄组、亚洲血统、男性等，找到满足这些任务的样本的机会变得微乎其微。

因此，采样进行推断是一种在低维度上提供可接受结果的方法。然而，随着维度数量的增加，或者如果抽样作业的概率较低，我们将需要转向更有效的推断方法。

马尔可夫链蒙特卡罗抽样过程

到目前为止，在我们的采样方法中，我们一直认为每个样本都是独立的。我们可以把镖靶想象成我们希望取样的表面，每一个随机取样都是一个被投掷在镖靶上的镖留下的洞。

我们知道，到目前为止，我们所学的一些采样方法效率很低，因此，设计有效的采样方法对于从难以采样的分布中推断任何统计数据至关重要。

我们可以使用一个叫做马尔可夫链蒙特卡罗 ( MCMC )的迭代采样过程。再用飞镖靶类比一下，想象一只蚱蜢在飞镖靶上跳跃。蚱蜢跳到的每一个新的地方都可以被认为是一个新的样本。使用这种迭代过程抽取的样本不是独立的，并且相互关联；然而，我们可以鼓励蚱蜢前往我们感兴趣的样本所在的区域。

在了解 MCMC 之前，我们先跑题了解一下马氏链的一些性质。

马尔可夫属性

马尔可夫属性表示，如果我们有一系列状态到，那么向状态的转换完全取决于状态。这就是所谓的一阶马尔可夫特性。

马尔可夫链

为了理解马尔可夫链，让我们以我们友好的蚱蜢为例。

蚱蜢有五堆大小不一的食物，经常从一堆食物到另一堆食物。每个转变(比如从 stash 1 到 stash 2)都有与之相关的概率，也有自我转变(蚱蜢可能从 stash 1 开始，四处游荡，然后回到 stash 1)。

下图描述了蚱蜢的过渡。每个边都有一个与之相关的非负转移概率。

你可能会问这样一个问题:在 1 号仓库找到蚱蜢的概率有多大？

达到稳定状态

事实证明如果我们观察蚱蜢足够长的时间，我们就会得到所谓的平稳分布。如果我们只观察蚱蜢一小段时间，我们可以看到它有 20%的时间在 stash 1，30%的时间在 stash 2，以此类推。然而，如果你观察蚱蜢足够长的时间，不管蚱蜢从哪个藏身地开始，它都会稳定下来，形成一个平稳分布，在时间 t 的概率几乎与时间 t + 1 相同(平稳分布也是谷歌 PageRank 算法的核心，在这个算法中，蚱蜢被一个网上冲浪者取代)。

蚱蜢是马尔可夫链的一个例子；它定义了所有状态的转换模型x。对于所有状态，它们的转移概率之和加起来是 1，等于 1。

从状态的转变被定义在转变矩阵中，并且转变矩阵的特征向量对应于该组状态上的唯一平稳分布。

为了使马尔可夫链收敛到平稳分布，它需要遵守以下一组属性:

• 它是有规律的，即从任意两对状态，比如说状态 x 到 y 行进的概率，正好在 k 步应该大于 0
• 从状态过渡出来的概率 x 应该大于 0，否则应该没有你不能过渡出来的节点
• 所有国家都应该有自我过渡。

使用马尔可夫链进行采样

建立马尔可夫链给了我们什么？使用蚱蜢藏匿点分布类比，我们可以将我们的初始信念想象为无信息先验分布(比如说，由于蚱蜢有五个藏匿点，在任何藏匿点找到蚱蜢的概率相等，也就是 20%)，在马尔可夫链达到收敛(蚱蜢已经徘徊足够长的时间)之后，我们可以从真实的后验分布中抽取样本。

此时一个合乎逻辑的问题就是“走得够长”到底有多长。不幸的是，没有简单的测试来确定我们是否达到了平稳分布(当我们达到平稳分布时，我们可以说链已经“混合”了)。然而，我们可以使用一些统计测试来找出一个链是否没有混合，当这些测试返回否定答案时，我们可以相信这个链确实混合了。

使用马尔可夫链对进行采样，我们首先让链老化，这是我们认为链还没有混合的时期，一旦我们确信混合已经发生，我们就可以收集样本并计算统计数据。

综上所述，使用马尔可夫链采样易于实现，具有良好的理论性质，但在实践中，该算法有许多可调性。有时候收敛很慢，判断一个链条是否混过，有点艺术。

吉布斯取样

吉布斯链是 MCMC 稳定的通用采样器，可用于从图模型中提取样本。当每个变量的条件分布已知或易于采样时，此方法适用。我们将在下表中演示取样程序的一个步骤，该步骤描述了从一个样品( x )移动到一个新样品( x' )的过程。

我们从吉布斯分布开始，即因子的乘积，或者来自贝叶斯网络，或者来自马尔可夫网络。状态空间是所有随机变量的完整赋值集。下面是我们打算采样的示例网络，它包含三个连接成 V 型结构的随机变量:

吉布斯取样程序中的步骤

以下是吉布斯取样程序中涉及的步骤:

• 从所有变量的当前赋值开始
• 对于 j = 1 的变量数，采样每个变量，而其他变量保持其先前的值

在下表中，我们遵循采样过程的每个步骤来获得一个吉布斯样本。我们假设三个随机变量是二进制值，当前赋值是。灰色单元格表示该值相对于前一行是固定的，白色单元格正在被采样。

|

循环

|

X1

|

X2

|

X3

|   |
| --- | --- | --- | --- | --- |
| 起始值(x) | Zero | Zero | Zero |   |
| 样品 | one | Zero | Zero | 抛硬币选择的值，我们得到 1 |
| 样品 | one | Zero | Zero | 抛硬币选择的值，我们得到 0 |
| 样品 | one | Zero | one | 抛硬币选择的值，我们得到 1 |
| x ' | one | Zero | one | 这是新的吉布斯样本 x’ |

最后一行是新状态x’，上表显示了 Gibbs 链生成下一个状态x’的单个步骤。

总之，吉布斯链是图模型最简单的马尔可夫链。吉布斯链不一定总是正则的，也就是说，它可能不会达到平稳分布，但在某些条件下(如正因子)，它保证是混合的。它是 PGM 最简单的马尔可夫链之一，但混合起来很慢，尤其是当概率达到峰值时。它只适用于我们可以从一个因素的产品中取样的情况。

吉布斯抽样的一个例子

让我们看一下Comparing Gibbs and Random Sampling IPython 笔记本中比较随机抽样和吉布斯抽样的代码片段。

在下面的代码片段中，我们将首先了解离散分布的平稳分布意味着什么，以及我们可以使用什么方法来更快地到达那里。

我们将使用熟悉的工作面试例子来主持讨论。

求职面试网络有五个变量，联合分布有 48 行(2×2×3×2×2，每个变量取值的个数)。我们对变量子集的边际分布感兴趣，我们也有一些观察到的证据。

我们对边际概率感兴趣，在边际概率中，我们观察了入场和体验随机变量的值。

在下面的代码片段中，我们使用精确推理(变量消除)来确定条件概率，然后打印相同的 CPD:

tcpd,bn,skel=getTableCPD()
query={'Offer':'0','Grades':'0','Interview':'0'}
evidence={'Admission':'0','Experience':'0'}
fac=tcpd.condprobve(query,evidence)
df=printdist(fac,bn)
df

前面代码的输出如下:

 Offer  Interview Grades Probability
0	 0	 0	 0	 0.641455
1	 0	 0	 1	 0.029455
2	 0	 1	 0	 0.064145
3	 0	 1	 1	 0.026182
4	 0	 2	 0	 0.000178
5	 0	 2	 1	 0.000109
6	 1	 0	 0	 0.071273
7	 1	 0	 1	 0.003273
8	 1	 1	 0	 0.096218
9	 1	 1	 1	 0.039273
10	 1	 2	 0	 0.017640
11	 1	 2	 1	 0.010800

要得到想要的分布，也就是我们首先要抽取样本，剔除不满足证据的。

libpgm库允许我们使用随机抽样和吉布斯抽样抽取样本。在这两种情况下，我们都可以通过证据进行条件。

在下面的代码中，我们使用吉布斯抽样和随机抽样抽取了 5000 个样本，并比较了从样本中学习到的边际概率:

def estimate_distrib(skel,samples):
    learner=PGMLearner()
    #learn the parameters of the network from the samples, given skeleton
    #returns a new bayes net.
    bayesnet=learner.discrete_mle_estimateparams(skel,samples)
    tablecpd=TableCPDFactorization(bayesnet)
    #run a conditional probability query for
    #P(Offer,Grades,Interview| Admission=0,Experience=0)
    fac=tablecpd.condprobve(query,evidence)
    #create a dataframe listing the marginals 
    df2=printdist(fac,bayesnet)
    return df2

#learn the marginals from gibbs samples
def gibbs_marginals(num_samples=5000):
    tcpd,bn,skel=getTableCPD()
    samples=tcpd.gibbssample(evidence,num_samples)
    df2=estimate_distrib(skel,samples)
    return df2['probability']

#learn the marginals from random samples
def random_sample_marginals(num_samples=5000):
    tcpd,bn,skel=getTableCPD()
    samples=bn.randomsample(num_samples,evidence)
    df2=estimate_distrib(skel,samples)
    return df2['probability']

df['prob from gibbs']=gibbs_marginals()
df['prob from random samples']=random_sample_marginals()
df

前面代码的输出如下:

Offer  Interview  Grades	probability	P: Gibbs P: random samples
0	 0	 0	 0	 0.641455	 0.645444	 0.078557
1	 0	 0	 1	 0.029455	 0.025156	 0.113443
2	 0	 1	 0	 0.064145	 0.065997	 0.058145
3	 0	 1	 1	 0.026182	 0.026203	 0.008655
4	 0	 2	 0	 0.000178	 0.000000	 0.013869
5	 0	 2	 1	 0.000109	 0.000000	 0.028531
6	 1	 0	 0	 0.071273	 0.072956	 0.048443
7	 1	 0	 1	 0.003273	 0.002844	 0.069957
8	 1	 1	 0	 0.096218	 0.096203	 0.504855
9	 1	 1	 1	 0.039273	 0.038197	 0.075145
10	 1	 2	 0	 0.017640	 0.016400	 0.000131
11	 1	 2	 1	 0.010800	 0.010600	 0.000269

最后三列列出了真实概率(从精确推断中获得)、来自吉布斯样本的概率和来自随机样本的概率。我们可以看到，吉布斯样本的概率相当接近真实的边缘值，随机样本与真实概率相差很大。

很明显吉布斯抽样是一个比随机抽样高效得多的抽样过程。然而，对于更大的维度，吉布斯抽样也将难以获得接近真实边际值的边际值。

总结

在这一章中，我们学习了通过精确推理解决棘手问题的方法。第一种方法是推理作为优化，其中之一是基于消息传递。

我们研究的第二种方法是基于粒子的推理(也称为采样)。我们学习了一旦维度增加，采样是如何挣扎的，我们还学习了诸如 MCMC 之类的方法，这些方法允许我们使用样本来获得期望的后验分布。

八、附录 a：参考文献

第一章

• 关于运动员使用药物的贝叶斯规则示例，请访问。
• 萨古尔·斯里哈里。查询联合概率分布(http://www . cedar . buffalo . edu/~ srihari/CSE 574/chap 8/Ch8-PGM-Directed/8 . 1 . 2-查询概率分布. pdf )。

第二章

• 如需了解 Libpgm，请访问https://pypi.python.org/pypi/libpgm

第三章

• Coursera PGM 课程可在https://www.coursera.org/course/pgm访问

第四章

• 周扬。概率图模型的结构学习:综合调查(http://arxiv.org/pdf/1111.6925.pdf)。
• 卡西欧·德·坎波斯、智曾加特、强吉。利用约束进行贝叶斯网络的结构学习(http://www . machine Learning . org/archive/icml 2009/papers/246 . pdf)。
• BNFinder 项目位置在http://bioputer.mimuw.edu.pl/software/bnf/。
• K .萨克斯、 O .佩雷斯、 D .佩尔、 D .劳芬伯格和 G .诺兰。“多参数单细胞数据衍生的因果蛋白-信号网络”(http://www . C2 B2 . Columbia . edu/danapeerlab/html/pub/science 2005 . pdf)。

第五章

• http://pymc-devs.github.io/pymc/的 PyMC 文件。
• 黑客的概率编程&贝叶斯方法。
• 参数估计章节、概率图模型。达芙妮·柯勒和尼尔·弗里德曼，麻省理工学院出版社。

第六章

• 信念网络中的推理:程序指南。的【黄】和的(交叉点树算法)。
• DAGs，I-图，因式分解，d-分离，最小 I-图，贝叶斯网络。尼尔·弗里德曼(http://classes.soe.ucsc.edu/cmps290c/Spring04/paps/nir2.pdf)。
• 贝叶斯信念网络项目页面位于https://github.com/eBay/bayesian-belief-networks。

第七章

• 有关 OpenGM 的信息，请访问https://github.com/opengm/opengm。
• 近似推理的循环信念传播:一项实证研究。凯文·墨菲、亚伊尔·韦斯和迈克尔·乔丹(http://www.cs.ubc.ca/~murphyk/Papers/loopy_uai99.pdf)。

其他参考文献

• 机器学习，概率视角。凯文·墨菲。
• 概率图模型。达芙妮·柯勒和尼尔·弗里德曼。
• https://www.coursera.org/course/pgm 的 Coursera PGM 课程。
• 贝叶斯人工智能。凯文·b·科尔布和安·e·尼克尔森在http://www.csse.monash.edu.au/bai/book/about.php。