贝叶斯网的R实现（ Bayesian networks in R）bnlearn(2)

3.结构学习

上面我们采用一个预先设定的结构建立了一个关于marks的贝叶斯网。这种方式在某些情况下（比如存在先验的专家知识）是合适的。但是对大多数的贝叶斯网络，我们需要从数据中学习网络。

3.1贝叶斯网的结构简介 
贝叶斯网关于节点（随机变量）的条件依赖或条件独立可以从图的角度讨论节点之间的连通与分割。

如果两个节点A,B直接相连，它们之间存在直接依赖关系。若两个节点不是直接相连的，要么A和B之间没有任何联系，即条件独立；要么可以通过其它节点在A，B之间建立起来联系，即条件依赖。

下面考虑两个节点A和B通过第三个节点C间接连接的情况，这有3种基本情形：

按图中顺序，自上到下分别成为： 分连（diverging connection） 顺连（serial connection） 汇连（converging connection）

对这三种连接形式，当给定节点C时，可以由贝叶斯网的Markov 性质解释条件依赖的情况：
 分连：P(A,B,C)=P(B|C)P(A|C)P（C） 
顺连：P(A,B,C)=P(B|C)P(C|A)P(A) 
汇连：P(A,B,C)=P(C|A,B)P(B)P(A) 对节点A,B分连和顺连是条件独立的，汇连是条件依赖的（C依赖于AB的联合分布）。 汇连这种情况也称为一个v-结构。

设X,Y,Z是3个不相交的节点集合，若在X,Y之间有一个节点V满足下列之一： （1）V在Z中且没有汇连 （2）V有一条汇连的边，且V或它的子代都不在Z中， 称Z d-分隔A,B（也称有向分隔）。 此时，若Z已知，则X和Y在给定Z时条件独立。

3.2结构学习算法 
结构学习算法可分为三类：基于约束的算法、基于得分的算法以及混合算法。

3.2.1基于约束的算法 基于约束的算法的基本想法是，节点A,B在节点C给定的条件下，是否满足条件独立。若是，则A和B被C d-分隔，A和B之间没有边；否则，A,B之间存在边。

基于约束的算法来源于Pearl关于因果图模型的工作：IC（inductive causation)算法。
 这个算法的基本步骤是这样的： 
（1）对成对的节点A和B寻找节点集合V的子集S：满足给定S是A,B独立且A,B不在S中。如果找不到这样的S，那么在AB之间加一条无向边。这一步可以看做是在一张完全图上利用条件独立的假设检验进行修剪的后向过程。
（2）对有共同邻居C的不相邻节点A和B，如果C也不在（1）的节点子集S中，那么把A-C和C-B由无向边改为有向边AB，CB（就是寻找v-结构）
 (3)迭代下面两步，对无向边设置方向： （a）如果A和B相连接，且有一个从A到B的有向路径，就把A-B改为A→B (b)如果A和B没连接，但A→C且C-B，则C-B改为C←B

当然，这个算法对于节点较多的真实问题，计算复杂度是很高的。所以应用的是一些基于IC算法发展的算法，在bnlearn包里使用的是这样一些算法（括号内是函数名）：
 Grow-Shrink (gs)
 Incremental Association (iamb) 
Fast Incremental Association (fast.iamb)
 Interleaved Incremental Association (inter.iamb)
 Max-Min Parents and Children (mmpc) 这些算法的细节不在赘述，需要指出的是，在bnlearn包中这些算法可以利用snow包实现并行计算。

在算法中为检验条件独立性，使用了一些统计或信息论中的方法。 对离散数据，基于CPT（条件概率表），可以计算：
 mutual information（mi） 
Pearson's X2（x2） 
Akaike Information Criterion（aict）
 对连续数据，则是基于偏相关系：
 linear correlation（cor/mc-cor） 
Fisher's Z（zf/mc-zf）
 mutual information (mi-g)

marks.bn <- gs(marks, undirected = FALSE)
marks.bn

## 
## Bayesian network learned via Constraint-based methods
## 
## model:
## [undirected graph]
## nodes: 5 
## arcs: 6 
## undirected arcs: 6 
## directed arcs: 0 
## average markov blanket size: 2.40 
## average neighbourhood size: 2.40 
## average branching factor: 0.00 
## 
## learning algorithm: Grow-Shrink 
## conditional independence test: Pearson's Linear Correlation 
## alpha threshold: 0.05 
## tests used in the learning procedure: 32 
## optimized: TRUE

marks.bn <- gs(marks, debug = TRUE) #debug=TRUE可以看到调试的过程
(结果略)


graphviz.plot(marks.bn, layout = "fdp")

利用gs（）得到的网络是个无向图（利用其它的基于约束的算法会得到类似的结果）。那么使用白名单（whitelist）和黑名单（blacklist）可以对gs（）得到的网络进行修正。黑名单是不该出现在图中的弧的集合而白名单则相反。 比如,稍作修正：

blacklisting <- data.frame(c("STAT", "STAT", "ALG"), c("ANL", "ALG", "MECH"))
marks.bn2 <- gs(marks, blacklist = blacklisting)
highlight.opts <- list(nodes = c("STAT", "ANL"), arcs = c("ANL", "STAT"), col = "blue", 
    fill = "grey", lwd = 2.5)
graphviz.plot(marks.bn2, layout = "fdp", highlight = highlight.opts)

3.2.2基于评分的算法 
运用优化中的贪婪算法，对每个“备选”的网络指定一个表示拟合优度的得分，然后取评分最大的网络。

bnlearn包里边使用评分法的函数只有一个：hc（），使用的是爬山法（Hill-Climbing）。
 这个方法的基本思路是： 
（1）选择一个网络结构G（通常是空的） 
（2）计算G的得分，定义score=score(G)
 (3)取maxscore=score(G) 
（4）随着maxscore的增加，重复下面的步骤： （i）对边的增加，删减或颠倒操作不会导致一个有环的网络。 （ii）计算修正的网络，更新maxscore 
（5）最后返回一个有向图

常用的评分方法有： 
likelihood (lik)/log-likelihood (loglik) Akaike (aic) 和 Bayesian (bic) Information Criterion Bayesian Dirichlet equivalent score (bde) 
K2 score (k2), 另一种 Dirichlet 后验密度

下面看hc( )的学习效果：

marks.bn1 <- hc(marks)
marks.bn1

## 
## Bayesian network learned via Score-based methods
## 
## model:
## [MECH][VECT|MECH][ALG|MECH:VECT][ANL|ALG][STAT|ALG:ANL] 
## nodes: 5 
## arcs: 6 
## undirected arcs: 0 
## directed arcs: 6 
## average markov blanket size: 2.40 
## average neighbourhood size: 2.40 
## average branching factor: 1.20 
## 
## learning algorithm: Hill-Climbing 
## score: 
## Bayesian Information Criterion (Gaussian) 
## penalization coefficient: 2.239 
## tests used in the learning procedure: 34 
## optimized: TRUE

graphviz.plot(marks.bn1, layout = "fdp")

score(marks.bn1, data = marks, type = "bic-g") #得分BIC

## [1] -1720

3.2.3混合算法 bnlearn包中的混合算法包括Max-Min Hill Climbing（mmhc）以及2-phase Restricted Maximization（rsmax2）