视频课程|R语言bnlearn包：贝叶斯网络的构造及参数学习的原理和实例
贝叶斯网络（BN）是一种基于有向无环图的概率模型，它描述了一组变量及其相互之间的条件依赖性。贝叶斯网络在信息不完备的情况下通过可以观察随机变量推断不可观察的随机变量，对于解决复杂的不确定性和关联性问题有很强的优势。
有向无环图（DAG）：DAG的结点V包括随机变量（类别和特征），有向连接E（A->B）表示结点A是结点B的parent，且B与A是有依赖关系的（不独立）。条件概率表集合：引入了一个条件性独立（conditionalindependence）概念，即图中任意结点v在给定v的parent结点的情况下，与图中其它结点都是独立的，也就是说P(v|par(v),x1,x2...,xn)=P(v|par(v))。这里par(v)表示v的parent结点集，x1,x2,...,xn表示图中其它结点。
1、确定随机变量间的拓扑关系，形成DAG。2、训练贝叶斯网络参数——估计出各节点的条件概率表。这一步也就是要完成条件概率表的构造。
1、基于评分搜索的方法：将结构学习视为组合优化问题，首先通过定义评分函数对BN结构空间中的不同元素与样本数据的拟合程度进行度量，然后利用搜索算法确定评分最高的网络结构，即与数据拟合最好的网络结构。
2、基于约束的方法：利用统计或信息论的方法定量的分析变量间的依赖关系以获取最优地表达这些关系的网络结构。首先对训练数据集进行统计测试，尤其是条件独立性测试，确定出变量之间的条件独立性，然后利用变量之间的条件独立性构造一个有向无环图，以尽可能多的涵盖这些条件独立性。
3、基于评分搜索和约束相混合的方法：首先采用条件独立性检验缩减搜索空间，然后将此搜索空间作为输入以限制基于评分的搜索。
将结构学习视为组合优化问题，首先通过定义评分函数对BN结构空间中的不同元素与样本数据的拟合程度进行度量，然后利用搜索算法确定评分最高的网络结构，即与数据拟合最好的网络结构。
1、爬山法（hc）：从一个初始的网络结构开始，通过加边、减边和转边三个搜索算子对当前网络结构进行修改，并对修改后的网络结构进行评分，最后找出评分最高的结构
2、禁忌搜索（tabu）:是一种亚启发式随机搜索算法，通过对局部领域逐步搜索得到全局最优解。为了避免陷入局部最优和重复迭代，禁忌搜索算法采用了一种禁忌表，利用该表对解的移动搜索过程进行记录和选择，当搜索的最优解已经在禁忌表中时，则会自动放弃该解，避免重复搜索。
MMHC算法结合了条件独立性测试、搜索和评分等方法,该算法分为2个阶段,在第1阶段利用MMPC(max-minparentsandchildren)算法构建贝叶斯网络结构的框架,在第2阶段执行评分搜索确定网络结构的边以及边的方向.
定义
贝叶斯网络（BNs）的定义是：
网络结构的主要作用是通过图形分离来表达模型中各变量之间的条件独立性关系，从而指定全局分布的因子化。
每个局部分布都有自己的参数集
;而
要比
小得多，因为许多参数是固定的，因为它们所属的变量是独立的。
R实现了以下学习算法。
基于约束的：PC,GS,IAMB,MMPC,Hilton-PC基于分数的：爬山算法、TabuSearch配对的：ARACNE,Chow-Liu混合：MMHC,RSMAX2
我们使用基于分数的学习算法，希尔算法。首先，我们将先为本教程生成简单的数据集。
接下来，我们将创建学习结构。
我们可以在一个图中看到结构。
>plot(hc_simd)
在这个图中，状态、元素、接受、类型和颜色被称为节点。节点之间的方向用弧线描述，弧线是一个包含从元素到元素方向数据的矩阵。
simd_fitted
基于上述训练数据，我们可以进行条件概率查询。
问题：受第三类错牙合畸形影响的患者（以下牙弓突出为特征），其骨骼不平衡在生命早期就产生，在青春期和骨骼成熟前会变得更加明显。在单个III类患者中早期预测治疗的成功或失败，使其更容易矫正，但仅从少量的形态决定因素中预测是很难做到的。原因是III类错颌畸形很少是单一颅面部件异常的结果，所以单个的临床和放射学测量值可能不如测量值本身的相互作用具有指示性。
任务:
我们将使用的数据集包含143名患者，在T1和T2年龄段有两组测量数据（以年为单位），用于以下变量。
所有的测量都是通过X射线扫描得出的，使用一套参考点建立的图，如以下。
>str(data)
首先，我们创建一个数据框架，其中包括所有变量的差异以及增长和治疗。
生长和治疗变量带有关于病人预后的冗余信息，这一点从TB和TG之间生长良好的病人比例的差异中可以看出。
为了避免在模型中包括这两个变量所导致的混杂，我们将治疗重新编码为一个二元变量，0表示NT，1表示TB或TG。同样地，我们对成长进行重新编码，0表示坏，1表示好。
由于我们将使用高斯BN进行分析，检查这些变量是否是正态分布；从下面的图来看，似乎并非所有的变量都是如此。
+hist(x,prob=TRUE)+lines(density(x),lwd=2)
这些变量是通过线性关系联系起来的吗？其中一些是，但不是全部。
>pairs(diff[,setdiff(names(diff)
最后，我们可以看看这些变量是否以任何方式聚在一起，因为聚在一起的变量更有可能在BN中发生联系。
>heatmap(rho)
我们可以在热图中看到两个集群：第一个集群包括dCoGo、dGoPg和dCoA，第二个集群包括Treatment、dANB和dCoA。第一个聚类在临床上很有意思，因为它包括治疗和两个都与唐氏A点有关的变量，这为治疗的主要效果提供了一些线索。
plot(ug)
在这里，我们使用保存在diff中的差异来为数据建模，而不是原始值；我们将使用GBN处理，因为所有变量都是数字。对差异进行建模会导致局部分布，其形式为回归模型
其中
对于其他回归因子，以此类推。我们可以将这种回归改写为
此外，包括增长变量意味着我们可以有以下形式的回归模型
从而允许不同的变化率，这取决于病人是否在畸形中表现出积极的发展，以及他是否正在接受治疗。
学习BN的第一步是学习其结构，即DAG
.我们可以使用数据（来自不同的数据框架）结合先验知识来做这件事；结合后者可以减少我们必须探索的模型空间，并生成更强大的BN。一个直接的方法是将那些编码我们知道不可能/真实的关系的弧列入黑名单;并将那些编码我们知道存在的关系的弧列入白名单。
黑名单只是一个矩阵（或一个数据框），其中有from和to两列，列出了我们不希望在BN中出现的弧。
白名单的结构与黑名单相同。
一个简单的学习
方法是在整个数据上找到具有最佳拟合度的网络结构。例如，使用hc()与默认分数（BIC）和整个diff数据框架。
至于绘图，关键函数是plot()。
plot(dag,,highlight)
然而，dag的质量关键取决于变量是否是正态分布，以及连接它们的关系是否是线性的；从探索性分析来看，并不清楚所有的变量都是如此。我们也不知道哪些弧线代表强关系，也就是说，它们能抵抗数据的扰动。我们可以用boot来解决这两个问题。
booth(diff,R=200)
boot.strength()的返回值包括，对于每一对节点，连接它们的弧的强度（例如，我们观察到dANB→dPPPM或dPPPM→dANB的频率）及其方向的强度（例如，当我们观察到dANB和dPPPM之间有弧时，我们观察到dANB→dPPPM的频率）。
因此，averaged.network()取所有强度至少为0.585的弧，并返回一个平均的共识网络，除非指定不同的阈值。
>avg.diff=averaged.network(str.diff)
纳入我们现在拥有的关于弧线强度的信息。
我们如何将平均的网络（avg.diff）与我们最初从所有数据中学习到的网络（dag）进行比较？最定性的方法是将两个网络并排绘制，节点位置相同，并突出显示一个网络中出现而另一个网络中没有的弧，或者出现的方向不同的弧。
我们可以看到，Treatment→dIMPa、dANB→dGoPg和dCoGo→dPPPM这些弧线只出现在平均网络中，而dPPPM→dANB只出现在我们从所有数据中学到的网络中。我们可以假设，前三个弧被数据的噪声加上小样本量和偏离常态的情况所隐藏。编程可以返回真阳性（出现在两个网络中的弧）和假阳性/阴性（只出现在两个网络中的一个的弧）的数量。
>compare
或弧=TRUE。
但是，考虑到网络是用BIC学习的，而BIC是等价的，那么所有的弧线方向是否都很确定？看一下dag和avg.diff的CPDAGs（并考虑到白名单和黑名单），我们看到没有无方向的弧。所有弧的方向都是唯一的。
最后，我们可以结合来进行原则性的比较，如果两个弧被唯一确定为不同，我们就说它们是不同的。
也可以看一下关于弧长分布的阈值：平均的网络是相当密集的（9个节点有17个弧），很难阅读。
因此，把阈值提高一点，多剔除几个弧就好了。看一下上面的图，由于弧长分布的差距，较高的阈值的两个自然选择是0.75（红色虚线）和0.85（蓝色虚线）。
我们通过在network()中设置阈值=0.85得到的更简单的网络如下所示；从定性的角度来看，它当然更容易推理。
在学习了结构之后，我们现在可以学习参数。由于我们正在处理连续变量，我们选择用GBN来建模。因此，如果我们使用最大似然估计来拟合网络的参数，我们就会发现每个局部分布都是一个典型的线性回归。
fit(avg,diff)
我们可以通过比较bn.fit()和lm()产生的模型，例如dANB，很容易确认这是事实。
>summary(lm(dANB~Growth+Treatment,data=diff))
我们会不会有共线性的问题？理论上是可能的，但在实践中，从数据中学习的网络结构大多不是问题。原因是，如果两个变量
和
是共线性的，在增加（比如说）Xi←Xj之后，那么Xj←Xk将不再显著提高BIC，因为Xj和Xk（在某种程度上）提供了关于Xi的相同信息。
如果参数估计因任何原因出现问题，我们可以用一组新的、来自不同方法的估计值来取代它们。
dANBdANB=penalized(dANB)dANB
有两种主要的方法来验证一个BN。
我们可以用通常的方法来衡量我们所选择的学习策略的预测准确性，即交叉验证。实现了：
对于:
首先，我们检查Growth，它编码了错牙合畸形的演变（0表示坏，1表示好）。我们检查它，把它转回离散变量并计算预测误差。
cv(diff)>for(iin1:10){+err[i]=(sum(tt)-sum(diag(tt)))/sum(tt)+}>
在这两种情况下，我们使用损失函数的变体，它使用从所有其他变量计算的后验期望值进行预测。基本的损失函数（cor,mse,pred）仅仅从它们的父代来预测每个节点的值，这在处理很少或没有父代的节点时是没有意义的。
确认BN是否有意义的另一种方法是把它当作的工作模型，看看它是否表达了关于关键事实，这些事实在学习过程中没有作为先验知识使用。否则，我们将只是拿回我们放在先验中的信息）。一些例子。
我们用于这个模型的数据是我们在分析开始时存储到正交的原始数据。然而，我们选择使用治疗变量而不是生长变量作为变量来表达受试者可能正在接受医疗的事实。原因是生长变量是一个衡量第二次测量时的预后的变量，它的值在第一次测量时是未知的；而治疗变量在两次测量时都是相同的。
>t2.variables
然后我们引入一个黑名单，其中。
相比之下，我们只对T1→T2的弧线进行白名单，因为第二次测量的年龄显然取决于第一次测量的年龄。
最后我们可以用bl和wl来学习BN的结构。
>dyn.dag
很明显，这个BN比前一个更复杂：它有更多的节点（16对9），更多的弧（27对19），因此有更多的参数（218对37）。
绘制这个新模型的最好方法是用plot()开始。
plot(dyn,render=FALSE)
然后，我们对变量进行分组，以方便区分const、t1.variables和t2.variables；我们选择从左到右而不是从上到下绘制网络。
与前一个模型一样，治疗作用于ANB：从治疗出去的唯一弧是治疗→ANB2和治疗→CoA2。同样，这两个子节点都与Down的A点有关。
我们想评估这个动态BN结构的稳定性，就像我们之前对静态BN所做的那样，我们可以再次做到这一点。
>boot(ortho)>plot(dyn)
avernet(dyn.str)
平均下来的avg和dag几乎是一样的：它们只相差两道弧。这表明结构学习产生了一个稳定的输出。
compare(dag,avg)tpfpfn2611
由于Treatment是一个离散变量，BN是一个CLGBN。这意味着以Treatment为父节点的连续节点的参数化与其他节点不同。
fit(dynavg)
我们可以对这个新模型提出另一组问题
相比之下，这是一个未经治疗的病人在相同初始条件下的模拟轨迹。
对CoA的模拟轨迹是比较现实的：它随着年龄的增长而减慢。这与ANB不同，它的发生是因为CoA2同时取决于T1和T2。(ANB2则两者都不依赖)。
>for(iinseq(nrow(interv)){+#进行联合预测，目前用predict()无法实现。+dist(dyn.fitted,nodes=c(),+intervals[i,]=weighted.mean(ANB2,weights)+intervals[i,]=weighted.mean(CoA2,weights)
在此对ChangZhang对本文所作的贡献表示诚挚感谢，她在河海大学完成了信息管理专业学位，专注数据挖掘、数据分析、自然语言处理。
若您是个人学习者或培训负责人，希望学习拓端的线上/线下课程，欲了解更多拓端学堂课程及信息，请点击文末“阅读原文”咨询。
最受欢迎的见解
1.matlab使用贝叶斯优化的深度学习
2.matlab贝叶斯隐马尔可夫hmm模型实现
3.R语言Gibbs抽样的贝叶斯简单线性回归仿真
4.R语言中的blockGibbs吉布斯采样贝叶斯多元线性回归
5.R语言中的Stan概率编程MCMC采样的贝叶斯模型
6.R语言贝叶斯Poisson泊松-正态分布模型分析职业足球比赛进球数
7.R语言使用贝叶斯层次模型进行空间数据分析
8.R语言随机搜索变量选择SSVS估计贝叶斯向量自回归（BVAR）模型