代码拉取完成，页面将自动刷新
在前一章中，我们已经了解到，随着图模型的树宽增加，精确的推断变得不可行。追求近似推理的动机来自真实世界的网络，在那里精确的推理是难以处理的。
在本章中，我们将学习计算近似推理的方法。我们将重温前一章中提到的消息传递算法，并了解当网络不是树形结构时，它们如何也可以用于计算近似推理。我们将从优化的角度探索推理，并学习采样方法。我们还将查看一些代码示例，以了解算法如何实现近似推理。
近似推理方法有两大类。在第一部分中，我们将探索作为优化的推理，在第二部分中，我们将探索基于粒子的推理方法。
在近似推理中，我们寻求找到或构建目标分布的近似。假设我们有真实的目标分布
，我们寻求找到一组容易运行推理的分布 Q，然后在这些容易的分布中搜索“接近”目标分布的那个。所使用的方法寻求优化用于测量 Q 和 p 之间的距离或相似性的函数
将推理问题归结为一个优化问题，可以让我们借鉴约束优化领域中研究得很好的方法。最常用的方法之一是生成一组描述目标函数最优值的方程。在图模型的上下文中，这种方法采用一组等式的形式，其中每个变量都是根据其他变量定义的。事实证明，将问题转化为一种约束优化的方程相当于图中的消息传递。
在这一节中，我们将探索信念传播算法。在前一章中，我们在树形结构的网络上运行了信念传播，而在这里，我们希望将讨论扩展到具有其他类型的图结构的网络。这个算法被称为 Loopy 信念传播 ( LBP )之所以这样命名，是因为聚类图可以包含循环，有时被称为 Loopy。
信念传播是用于执行推理的消息传递算法之一。尽管它最初被设计为在树形结构的网络上运行，但人们发现它也可以用于具有循环或循环的一般图。
在图结构上运行算法与在树结构上运行算法有一些不同。
|树形网络
|
图结构网络
| | --- | --- | | 一个节点有首先从它所有的叶子节点接收消息，然后只向它的父节点发送一条消息，整个过程向下重复。 | 虽然树形结构网络中的节点可以在向其父节点发送消息之前等待其子节点的消息，但是在有环路的网络中，节点在发送消息之前等待消息并没有简单的过程。 | | 每个节点和它的父节点之间只需要传递两条消息，就可以实现信念的收敛。 | 消息传递一直持续到达到收敛。 | | 消息传递导致集群和分离集中信念的收敛。 | 大多数情况下会发生收敛，但这并不能保证。有些网络可能会出现振荡，有些可能永远不会收敛。 | | 估计边缘收敛到真实边缘。 | 即使在收敛时，边缘也可能是不正确的，或者集群信念不一定等于真正的边缘。 |
运行 LBP 的第一个任务是当给定一组因子时，创建一个聚类图。每个因子必须分配给一个聚类，但是一个聚类可以分配多个因子。从一组因子到聚类图的转换不是唯一的。但是，对于执行信念传播的集群图，它必须满足以下属性:
让我们看一个集群图创建的例子。
我们有如下公式所示的一组因子(以及每个因子范围内的变量):
下图显示了一种可能的集群图表示。分离集中的变量在连接每个聚类的边上标注。
注意前面提到的两个属性在集群图中是如何满足的。
运行交叉点:集群 5 、 2 和 3 在其范围内包含变量 D 。运行交叉口要求，如果 5 和 3 在其范围内包含 D ，它们之间的簇也应该包含 D (簇 2 )。
家族保存:聚类因子赋值如下表所示(注意聚类的范围包括其赋值因子范围内的所有变量):
|簇
|
工厂
|
| --- | --- |
| 1:甲、乙、丙 | 
|
| 2:C、D | 
|
| 3 :D、F、H | 
|
| 4:英、法 | 
|
| 5:乙、戊、丁 | 
|
在传递任何消息之前，每个集群的初始信念只是分配给它的所有因素的产物。
消息传递以类似于我们在上一章中看到的方式进行。首先，所有消息都被赋值 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 现有的信仰为例
将传入消息和相乘
Summing out variable D and calculating the message by using the following formula:
发送到 C1
请注意，当消息从集群 2 返回到 1 时，它避免发送集群 1 首先发送的信息，否则从集群 1 发送的信念只是来回呼应(重复计数)，并且在这个过程中变得更强。
通过的消息以前述相同的方式继续，并在收敛时停止。当跨时间步长的分离集信念变化在预定的容差值内时，可以假设聚类图已经收敛。此外，每对集群信念与分离集信念一致(对于位于这对集群之间的边缘)。
我们可以将在 LBP 算法中的步骤总结如下:
将因子分配给集群
构建初始势(将一个簇中的所有因子相乘)
初始化每个集群的信念
重复消息传递步骤:选择一个集群对，并在它们之间传递消息
在 n 次重复后，测试收敛是否已经发生或停止
Summarize beliefs at each cluster by multiplying all the messages received with their initial beliefs as follows:
这里，第 I 个集群中的信念是初始信念
的产物，也是来自所有 k 相邻集群的所有消息
的产物。
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。
在您尝试图像分割之前，让我们尝试使用 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函数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。
从状态的转变被定义在转变矩阵中，并且转变矩阵的特征向量对应于该组状态上的唯一平稳分布。
为了使马尔可夫链收敛到平稳分布，它需要遵守以下一组属性:
建立马尔可夫链给了我们什么？使用蚱蜢藏匿点分布类比，我们可以将我们的初始信念想象为无信息先验分布(比如说，由于蚱蜢有五个藏匿点，在任何藏匿点找到蚱蜢的概率相等，也就是 20%)，在马尔可夫链达到收敛(蚱蜢已经徘徊足够长的时间)之后，我们可以从真实的后验分布中抽取样本。
此时一个合乎逻辑的问题就是“走得够长”到底有多长。不幸的是，没有简单的测试来确定我们是否达到了平稳分布(当我们达到平稳分布时，我们可以说链已经“混合”了)。然而，我们可以使用一些统计测试来找出一个链是否没有混合，当这些测试返回否定答案时，我们可以相信这个链确实混合了。
使用马尔可夫链对进行采样，我们首先让链老化，这是我们认为链还没有混合的时期，一旦我们确信混合已经发生，我们就可以收集样本并计算统计数据。
综上所述，使用马尔可夫链采样易于实现，具有良好的理论性质，但在实践中，该算法有许多可调性。有时候收敛很慢，判断一个链条是否混过，有点艺术。
吉布斯链是 MCMC 稳定的通用采样器，可用于从图模型中提取样本。当每个变量的条件分布已知或易于采样时，此方法适用。我们将在下表中演示取样程序的一个步骤，该步骤描述了从一个样品( x )移动到一个新样品( x' )的过程。
我们从吉布斯分布开始，即因子的乘积，或者来自贝叶斯网络，或者来自马尔可夫网络。状态空间是所有随机变量的完整赋值集。下面是我们打算采样的示例网络，它包含三个连接成 V 型结构的随机变量:
以下是吉布斯取样程序中涉及的步骤:
变量，而其他变量保持其先前的值在下表中，我们遵循采样过程的每个步骤来获得一个吉布斯样本。我们假设三个随机变量
是二进制值，当前赋值是
。灰色单元格表示该值相对于前一行是固定的，白色单元格正在被采样。
循环
|
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 之类的方法，这些方法允许我们使用样本来获得期望的后验分布。
此处可能存在不合适展示的内容，页面不予展示。您可通过相关编辑功能自查并修改。
如您确认内容无涉及 不当用语 / 纯广告导流 / 暴力 / 低俗色情 / 侵权 / 盗版 / 虚假 / 无价值内容或违法国家有关法律法规的内容，可点击提交进行申诉，我们将尽快为您处理。