2.2　概率图模型

概率图模型是通过以图为工具表达变量关系的概率模型，一般一个节点表示一个或一组随机变量，变量之间的概率关系通过节点间的边表示。概率图模型可以通过有向无环图和无向图表示变量间的依赖关系，前者称为有贝叶斯网或有向图模型，后者称为马尔可夫网或无向图模型。

2.2.1　隐马尔可夫模型

马尔可夫过程表示数学中具有马尔可夫性质的离散随机过程，在该过程中，下一时刻的状态只与当前状态有关，与上一时刻状态无关。最简单的马尔可夫过程就是一阶过程，每一个状态的转移只依赖其之前的那一个状态。

假定一天的天气状况有晴（高温状态）、阴（低温状态）、雨（凉爽状态）三种，这三件事按图2.1所示的方向转移：

这个过程就是一个简单的马尔可夫过程，但马尔可夫过程和确定性系统不同，状态之间的转移是有概率的，天气的状态是经常变化的，而且会在两个状态间任意切换，如图2.2所示。

图2.2中箭头表示天气从一个状态切换到另一个状态的概率，保持晴天的概率是0.5，由晴天转为阴天的概率为0.4，由晴天转为下雨的概率为0.1。

当前状态的转移依赖之前的n个状态，当n为1时即为马尔可夫假设；并由此可以得出马尔可夫链的如下定义：

马尔可夫链是一组具有马尔可夫性质的离散随机变量的集合S=｛St∶t＞0｝，随机变量的所有可能取值的集合被称为状态空间，St的值则是在时间t的状态。若St+1对于过去状态的条件概率分布仅是St的一个函数，即称S为马尔可夫链。

这里定义的马尔可夫链是离散时间马尔可夫链，具有连续指数集的情形称为马尔可夫过程。

隐马尔可夫模型是描述含有隐含未知参数的马尔可夫过程，从可观察的参数中确定隐含参数，然后利用这些参数来做进一步的统计与预测，是最简单的贝叶斯网，多用于动态时序数据建模，被建模的系统被认为是一个马尔可夫过程与未观测到的（隐藏的）状态统计。

隐马尔可夫模型包含两个状态集合，分别为可观测的状态集和隐藏的状态集，通过可观测状态预测隐藏状态，如图2.3所示。

由此可以通过观测状态温度、大气压和降水等观测状态O=｛OT，OA，OR｝来判断隐藏状态S=｛SF，SC，SR｝。在晴天状态下，观测到温度高的概率高，观测到气压低的概率低，观测到降水的概率低；如果观测到温度低的概率高，气压低的概率高，降水的概率较低的时候，则大概率是阴天的状态；如果观测到温度低的概率低，气压低的概率低，降水的概率高的时候，则大概率是下雨的状态。由此可知，天气可观测的状态序列和隐藏的状态序列是概率相关的，则可以将此类型的过程建模为一个隐藏的马尔可夫过程和一个与之概率相关的可观测状态集合，即隐马尔可夫模型。

隐马尔可夫模型可以通过混淆矩阵表示，矩阵的行表示隐藏状态，每一行的概率值之和为1，列表示可观测状态，图2.3的混淆矩阵为：

上述问题的状态转移矩阵为：

隐马尔可夫模型可以用5元组｛N，M，π，A，B｝表示，N表示隐藏状态的数量，M表示可观测状态的数量，π=｛πi｝表示初始隐藏状态的概率，A表示隐藏状态的转移矩阵，B表示混淆矩阵。

图2.3中隐藏状态N=3，可观测状态M=3。

在实际应用中，给定状态转移矩阵、混淆矩阵和初始状态概率可确定隐马尔可夫模型，可按照如下过程产生相应的观测序列：

（1）设置时间步为1，根据初始状态的概率选择初始状态值；

（2）根据当前状态值和观测概率选择观测变量；

（3）根据当前状态值和状态转移矩阵选择下一时间步的状态值；

（4）直至满足确定的时间步要求。

常见的和隐马尔可夫模型相关问题如下：

（1）给定隐马尔可夫模型的初始状态、状态转移矩阵和混淆矩阵，计算产生给定观测序列的概率。即根据已有的观测序列推测当前时刻最可能出现的观测值。

（2）给定隐马尔可夫模型的初始状态、状态转移矩阵、混淆矩阵和观测序列。计算与之对应的隐藏状态序列，即根据观测信号推断最大可能的隐藏状态序列。

（3）根据给定的观测序列，计算隐马尔可夫模型的参数，且使得出现给定序列的概率最大，常用于通过学习样本来得到模型参数，以便更好地描述观测数据。

2.2.2　贝叶斯网络

贝叶斯网络也称信念网络或有向无环图模型，是一种概率图模型，于1985年由Judea Pearl首先提出，是一种模拟人类推理过程中因果关系的不确定性处理模型。有向无环图可表示其网络拓扑结构，图中的节点表示随机变量，可以表示可观察到的变量、隐变量、未知参数等。有因果关系或非条件独立的变量和命题在有向无环图中通过箭头对表示“因”和“果”的两个节点进行连接，连接在一起的两个节点会产生一个条件概率值。例如，假定节点A直接影响到节点B，通过A→B来表示，可用从节点A指向节点B的箭头来构建A到B的有向弧（A，B），权值可用条件概率来表示，用圆圈表示随机变量，用箭头表示变量之间的条件依赖，如图2.4所示。

将某系统中涉及的全部随机变量，根据变量之间是否条件独立表示在一个有向图中，用来表示随机变量之间的条件依赖，即构成了贝叶斯网络。G=（I，E）表示有向无环图，式中I表示有向无环图中所有的节点的集合，E表示有向连接弧的集合，假设X=（Xj），j∈I表示图2.4中的某节点j表示的随机变量，则节点X的联合概率可以表示为

X即为相对于有向无环图G的贝叶斯网络。

概率的先验独立性和条件独立性如下：

先验独立性：P（X，Y）=P（X）×P（Y）

条件独立性：

采用多个条件概率的乘积来计算联合概率的方法叫作链式法则，每一个贝叶斯网络都对应了一套链式法则来计算联合概率，同时每一个链式法则在忽略弱依赖关系时，对于任意随机变量，联合概率可由各局部条件概率分布相乘求得：

简单贝叶斯网络如图2.5所示：

因为A导致B，A和B导致C，所以有下式：

贝叶斯网络共有三种结构形式，分别为head-to-head、tail-to-tail和head-to-tail。

贝叶斯网络如图2.6所示。

从图2.6可以比较直观地看出：

x1，x2，…，x7的联合分布为：

（1）；

（2）x1、x2独立；

（3）x4和x5在x6给定的条件下独立。

拓扑结构1：汇连结构（head-to-head）。汇连结构形式如图2.7所示：

则有成立，化简可得：

可得：

在x7未知的条件下，x1、x2是被阻断的，称为head-to-head条件独立。

拓扑结构2：分连结构（tail-to-tail）。分连结构形式如图2.8所示：

考虑C未知和C已知这两种情况：

在x6未知时，有，此时无法得到p（x4，x5）=p（x4）p（x5），即x6未知时，x4、x5不独立。

在x6已知时，有，将p（x4，x5，x6）=代入式中，可以得到：即x6已知时，x4、x5独立。因此，当x6已知时，x4、x5是被阻断的，称为tail-to-tail条件独立，即x4和x5在x6给定的条件下独立。

拓扑结构3：直连结构（head-to-tail）。直连结构形式如图2.9所示：

当x6未知时，有，无法推出p（x1，x4）=p（x1）p（x4），即x6未知时，x1、x4不独立。

当x6已知时，有，根据p（x1，，可得：

即在x6已知时，x1、x4被阻断（blocked），称为head-to-tail条件独立。

显然，head-to-tail其实就是一个链式网络，在x6给定的条件下，x4的分布和x1条件独立，也就是x4的分布状态只和x6有关，和其他变量条件独立。即当前状态只跟上一状态有关，这一随机过程叫作马尔可夫链。

例2.1：煤炭是工业的食粮，是人类的生产生活必不可缺的能量来源之一，煤炭的供应也关系到各国的工业乃至整个社会发展的稳定。造成煤炭价格波动的因素较为复杂，同时煤炭价格的波动也会影响其他相关行业的健康发展。当煤炭的价格上涨时，有80%的可能性会引起当地电价的上涨，并且有40%的可能性会带来失业率的上升；如果煤炭价格稳定，引起电价上涨和失业率提高的可能性只有10%。若煤炭的产量持续增加且不出现长期的高温天气，煤炭有70%的可能性保持价格稳定；如果产量下降或出现持续高温天气，煤炭有60%的可能性保持价格稳定。当产量下降的同时又伴随持续高温，则煤炭只有30%的可能性保持价格稳定。根据统计数据得知，煤炭产量持续增加的可能性为50%，同时有20%的可能性出现持续高温天气。

上述问题可以通过图2.10有向无环图表示，椭圆节点表示变量，有向边表示变量之间的依赖关系，起始节点表示父节点，终止节点为子节点，比如图2.10中Y节点是C节点的父节点，C节点为Y节点的子节点。每个节点的概率值以条件概率表的形式给出。条件概率表是在父节点约束下的子节点条件分布，变量以二值的形式取值，比如Y节点取值T时表示煤炭产量稳定，取值为F时表示产量不稳定，则该问题各变量的联合概率分布为P（Y，H，C，E，R）。根据条件概率公式和变量的独立性，可以得到如下的等量变换：

2.2.3　贝叶斯网络推理

贝叶斯网络推理算法一般分为精确推理算法和近似推理算法。其中精确推理算法希望能计算出目标变量的边际分布或条件分布的精确值，但此类算法的计算复杂度随着极大团规模的增长呈指数增长，不适用于规模较大的贝叶斯网络，因此当贝叶斯网络的规模较大时，一般采用近似推理，以便在较低时间复杂度下获得问题的近似解。常见的贝叶斯精确推理算法有：多树传播（Polytree Propagation）推理算法；团树传播（Clique Tree Propagation）推理算法，比如联结树（Junction Tree Propagation）推理算法；基于组合优化的求解算法，如符号推理（Symbolic ProbabilisticInference）和桶消元推理（Bucket Elemination Inference）算法。近似推理算法主要有：基于搜索的（Search-Based）方法、Monte Carlo算法等。

2.2.3.1　精确推理

贝叶斯网络精确推理的常用方法有变量消元法、枚举推理法等，本质上是通过动态规划算法来计算目标概率值，通过给定的证据变量计算查询变量的后验概率分布。变量消元法是一种基于变量的条件独立性的精确推理算法，通过计算边际概率消去非计算目标的概率值，如图2.11所示。

变量消元法就是对联合概率求和消去其他变量得到边缘分布，如图2.11首先把z消元，得到仅含x和y的表，然后再对y进行求和得到只含x的概率表。

以图2.12中的贝叶斯网为例，计算P（D2），则有：

为了降低推理计算的复杂度，可以分解联合分布的因子。与变量A有关的因子是P（A）和，与变量B有关的因子是和，与变量C有关的因子是和，因此，有：

消元变量A只涉及自身和变量B，消元变量B只涉及自身和变量C，消元变量D1只涉及自身和变量C，消元变量C只涉及自身和变量D2，因此可以通过局部运算有效降低推理的复杂度，可以如下描述消元法的算法。

（1）设贝叶斯网中概率分布的集合为Ω，有：

（2）将集合Ω中含变量A的因子P（A）和合并求和得到新因子：

并将新因子f1（B）加入集合Ω，同时消去因子P（A）和，达到消元变量A的目的，得到更新的集合Ω：

（3）将集合Ω中含变量B的因子f1（B）和合并求和得到新因子：

并将新因子f2（C）加入集合Ω，同时消去因子f1（B）和，达到消元变量B的目的，得到更新的集合Ω：

（4）将集合Ω中含变量D1的因子计算得到新因子：

并将新因子f（C）加入集合Ω，同时消去因子f1（B）和，达到消元变量D1的目的，得到更新的集合Ω：

（5）将集合Ω中含变量C的因子f2（C）、f3（C）和合并求和得到新因子：

达到消元变量C的目的。

（6）返回f4（D2），即问题的目标P（D2），如图2.12所示。

变量消元法计算的复杂度和消元的顺序有关系，而最优的消元顺序是NP-hard问题，所以实际计算中一般采取启发式的规则进行计算，常用的算法有最大势搜索和最小缺边搜索。

图2.13是包含8个节点的无向图，最大势搜索是对图中所有节点按照如下规则进行编号，在t步中选择具有最多已编号相邻节点的未编号节点，若有多个则选择其中一个并将编号定为8-t+1。当所有节点被编号后，按照从小到大的顺序排序节点，即可得到变量消元顺序。

如图2.13，任意选择节点H，并编号为8。下一可选节点为有1个已编号邻接节点H的E、F或G，任意选择节点F，并编号为7。下一节点则只能选择有2个已编号邻接节点H和F的E，编号为6。下一节点选择有2个已编号邻接节点E和H的G，编号为5。下一可选节点有1个已编号邻接节点的C、D或A，任意选择节点A，并编号为4。下一节点选择有2个已编号邻接节点A和E的D，编号为3。下一节点选择有2个已编号邻接节点D和E的C，编号为2。最后将节点B编号为1，得到最终的消元顺序为＜B，C，D，A，G，E，F，H＞。

变量消元算法的缺点主要表现在需要计算多个边缘分布，多次执行变量消元算法，中间会有很多结果可以复用，但是变量消元算法需计算多次，导致效率低下。

2.2.3.2　近似推理

精确推理需要大量的计算开销，在实际应用中面对较大规模的问题时近似推理更为常用。近似推理的核心思想是通过随机采样得到相应的样本集，进而通过样本的频率来逼近查询概率。蒙特卡洛方法也被称为统计模拟方法、统计试验法、随机抽样技术等，该方法于20世纪40年代由美国在第二次世界大战中研制原子弹的“曼哈顿计划”的成员S.M.乌拉姆和J.冯·诺伊曼首先提出，冯·诺伊曼用赌城摩纳哥的Monte Carlo来命名。广泛应用于人工智能、金融工程学、宏观经济学、生物医学、计算物理学，诸如粒子输运计算、量子热力学计算、空气动力学计算、核工程等领域。常见的采样方法有直接采样方法、接受—拒绝采样方法、重要性采样方法、MCMC采样方法、Metropolis-Hastings采样方法等。

（1）直接采样方法。

蒙特卡洛方法是通过随机模拟的方式求解不太容易求解的和或者积分问题，使用采样+平均核心思想估计无法计算的值。

直接采样的方法是根据概率分布进行采样。对一个已知概率密度函数与累积概率密度函数的概率分布，我们可以直接从累积分布函数进行采样。

在例2.1中，考虑如下推理问题：假设观察到失业率提高，计算该结果由于煤炭产量下降引起的概率是多少。针对此问题，首先可以根据煤炭产量和高温天气的先验概率生成一个样本，并记录样本的取值；其次以样本的取值为条件，按照煤炭价格上涨的条件概率生成样本，并记录样本的取值；最后根据独立向分别计算在煤炭价格为条件时生成电价和失业率的样本。产生大量样本后，统计取值＜*，r，*，y，*＞和＜*，r，*，*，*＞的样本个数，分别记为N（*，r，*，y，*）和N（*，r，*，*，*），根据概率的定义，当采集的样本量较大时，则N（*，r，*，y，*）/N（*，r，*，*，*）即可作为近似结果输出，当样本总量趋于无穷时，结果即为问题的概率。在实际应用中，由于很多分布无法写出概率密度函数与累积分布函数，因此直接采样方法的适用范围较窄。

例2.2：应用蒙特卡洛方法求圆周率，如图2.14所示。圆的半径为r=1，正方形的边长为1，根据圆面积的计算公式可知圆的面积为πr2=π，同时正方形内部的相切圆的面积为整个圆面积的1/4，也就是。正方形的面积为1。然后向正方形中随机打点，则点就会以一定的概率落在1/4圆中，而点落在1/4圆中的概率为：

圆的面积/正方形面积=，

然后即可推出圆周率的计算公式为：

（2）接受—拒绝采样方法。

接受—拒绝采样方法用于累积分布函数未知的问题。如图2.15所示，p（z）是希望采样的分布，很多实际问题中，p（z）因太复杂在程序中没法直接采样，需要借助其他的手段来采样，因此设置提议的程序可抽样分布q（z），比如高斯分布，且kq（z）＞p（z）。在kq（z）中按照直接采样的方法采样，然后判断这个样本的取值是否符合给定观测，不符合给定观测的样本予以拒绝，符合给定观测的样本予以接受，最终得到符合p（z）的N个样本。样本量越大，得到的查询概率越接近真实值，当样本量区域无穷大时得到的查询概率和真实值严格相等。接受—拒绝采样方法因为在采样过程中要丢弃大量不符合观测值的样本，所以计算效率不高。

（3）重要性采样方法。

接受—拒绝采样解决了累积分布函数未知的采样问题，但采样的效果取决于提议分布的选择，如果提议分布选择得不好，则会出现计算效率低下且难以得到符合观测值的样本。与直接采样与接受—拒绝采样将每个样本默认相同的权重不同，重要性采样给予每个样本不同的权重，使用加权平均的方法来计算期望。

通过从提议分布q（z）中采样样本z1，z2，…，zn，每个样本的权重是，通过加权平均的方式计算期望，如：

（4）MCMC抽样方法。

马尔可夫链是概率论和数理统计中具有马尔可夫性质且存在于离散的指数集和状态空间内的随机过程，也就是通过时间串联的一系列随机变量，可通过转移矩阵和转移图定义。马尔可夫链可被应用于蒙特卡洛方法中，形成马尔可夫链蒙特卡洛采样方法，即MCMC方法。

举例说明马尔可夫链。按照十年河东十年河西的说法，穷人和富人之间会在一定概率下相互转变。为便于计算，假定穷人人口和富人人口每年会发生一次转变，且转变概率是固定的。假定每年富人转变为穷人的概率是4%，穷人转变为富人的概率是8%。如果某一年，某区域的富人和穷人的人口分别是5000和35000，那么该区域下一年穷人和富人的分布如何呢？

首先构造转移矩阵：

则富人和穷人的分布为：

因此，每个人作为一个个体，其个人的穷富状态其实根据个人的奋斗情况在转变：比如t=1时，是穷人；t=2时，经过奋斗则会变成富人。马尔可夫链就是生成这样一段状态序列的随机过程，其中城市和农村互相流动的矩阵，叫作迁移矩阵。马尔可夫链的这个随机过程满足马尔可夫性质，也就是某一个状态的值，只跟前一个状态相关，即一阶齐次马尔可夫链。用公式表示就是：

在迁移矩阵中，每一个元素对应一个条件概率值，经过一定阶段的迭代之后，最终穷人和富人的人数相对固定了，富人是26667，穷人是13333。对于该区域而言，只要总人口是40000，无论初始的穷人和富人比例如何，最终都会收敛到富人是26667，穷人是13333。简单分析原因，穷人经过奋斗会变得富有，而富人如果堕落则会变得穷困。随着时间的推移，穷人逐渐减少，而富人逐渐增加，但最终穷人变富和富人变穷的数量是相等的，即两类人群的流入和流出是相等的，达到一种相对平稳的分布，这个分布就是马尔可夫链的平稳分布。

马尔可夫链收敛于平稳分布π，π是方程πP=π唯一非负解。把这个解用向量方式表示，即：

且，。

因此，对于一个非周期的马尔可夫链转移矩阵P，其任何两个状态是连通的，存在，并且与i独立，记为，则有：

π是方程πP=π唯一非负解：

且，π=［π（1），π（2），…，π（j），…］，

则π成为马尔可夫链的平稳分布。

MCMC方法的核心是对于任意给定的目标分布，找到以该目标分布为稳态分布的马尔可夫链，并基于马尔可夫链采样方法对其近似采样。

（5）Metropolis-Hastings采样方法。

若要基于给定的概率分布高效生成对应的样本，由于马尔可夫链可以收敛到平稳分布，最直观的需要是构造一个转移矩阵为P的马尔可夫链，使得该马尔可夫链的平稳分布是π（x），则从任何一个初始状态均可以得到相应的转移序列，收集马尔可夫链在收敛之后的样本。

所以，基于马尔可夫链采样的关键问题是构造转移矩阵，使得转移矩阵满足一定的条件，使平稳分布恰好是目标分布P（x）。

转移矩阵需要满足的这个条件就是细致平稳条件，细致平稳条件是指对于给定的马尔可夫链的转移矩阵P和分布π，以及马尔可夫链中的任意两个状态x和x*，如果有：

则分布π即为该马尔可夫链的平稳分布。

简单分析细致平稳条件，对于条件中等式两侧关于状态x求和，则有：

则有，可得到平稳分布的条件πP=π，即满足细致平稳条件的目标分布是基于转移矩阵P的马尔可夫链平稳分布。

Metropolis-Hastings采样方法提供了找到满足细致平稳条件的转移矩阵P的思路。

在一般情况下，对于给定的目标概率π，任意转移矩阵P不满足π（x），即此时的转移矩阵不满足细致平稳条件，Metropolis-Hastings方法试图通过引入α（x，x*）和α（x*，x）因子来满足细致平稳条件，即：

为方便表示，将记为p（x，x*），将记为p（x*，x）。为了使上式成立，按照对称性，只需：

于是具备了设置因子α（x，x*）和α（x*，x）的可能性，从而等式π（x）p（x，x*）α（x，x*）=π（x*）p（x*，x）α（x*，x）可以得到满足，对其进行如下改写：

则有：

于是就将原转移矩阵改造成了转移矩阵Q，且Q满足细致平稳条件，引入的因子被称为接受率，物理意义是以此接受率为概率接受转移。以接受率α（x，x*）为例，当从状态x转移到状态x*时，以α（x，x*）为概率接受这个转移，则新的马尔可夫链的转移概率由原来的p（x，x*）变为p（x，x*）α（x，x*）。如图2.16（a）所示。

状态2以0.4的概率向状态1转变，以0.5的概率向状态3转变，以0.1的概率保持在状态2。叠加上接受率α（2，1）=0.6和α（2，3）=0.8之后，如图2.16（b）所示。

叠加上接受率之后图示的转移概率发生了改变，状态2到状态1的转移概率原为0.4，因为叠加了接受率α（2，1）=0.6，则状态2到状态1的转移概率变为0.4×0.6=0.24，有0.4×0.4=0.16的概率仍保留在状态2；状态2到状态3的转移概率原为0.5，因为叠加了接受率α（2，3）=0.8，则状态2到状态3的转移概率变为0.5×0.8=0.4，有0.5×0.2=0.1的概率仍保留在状态2；而保留状态2的概率则由原来的概率0.1增大到0.36（0.1+0.16+0.1=0.36）。

在马尔可夫链蒙特卡洛方法中的Metropolis-Hastings方法中进一步明确了接受概率α的表达式：

以图2.16为例，假定接受率α（2，1）=0.1和α（2，3）=0.2，在实际采样过程中效率会很低，这是因为在等式π（x）p（x，x*）×0.1=π（x*）p（x*，x）×0.2中等号两端的接受率过小，此时如果对等号两端的接受率等比例放大，比如：

显然仍满足细致平稳条件，而采样的效率可以得到大幅提高。