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

1.bayesian networks的一些基本概念 

贝叶斯网bayesian networks是一种有向无环图模型（DAG），可表示为G=(V,A)。其中V是节点的集合，节点表示随机变量；A是弧（或称为边）的集合，弧的箭头表示随机变量之间的概率相依性。有向无环图DAG定义了一个因子化的V中全体节点的联合概率分布，称为全局概率分布；相对的，与每个随机变量关联的，为局部概率分布。

因子化的形式由贝叶斯网的马尔科夫性质给出，对每个随机变量，其概率只依赖于其父代：

可以通过学习算法得到贝叶斯网的结构。学习算法首先是学习网络结构，然后在此基础上估计局部分布函数的参数。 尽管全局和局部分布的选择形式很多，最常用的还是下面的两种分布： 多项分布（离散数据） 多元正态分布（连续分布）

得到网络结构，可以利用概率的条件相依，对相关专业问题进行推断。

R中有多个包可以实现贝叶斯网的创建、学习和推断，详见下面表格。

前面我曾简单介绍过gRain包，本篇讨论最常用的bnlearn包。

2.网络创建和操作 bnlearn包自带数据集marks，88学生5门课的成绩。这个数据最早由Mardia et al(1979)研究过。

library(bnlearn)

## Warning: package 'bnlearn' was built under R version 3.0.1

data(marks)
str(marks)

## 'data.frame': 88 obs. of 5 variables:
## $ MECH: num 77 63 75 55 63 53 51 59 62 64 ...
## $ VECT: num 82 78 73 72 63 61 67 70 60 72 ...
## $ ALG : num 67 80 71 63 65 72 65 68 58 60 ...
## $ ANL : num 67 70 66 70 70 64 65 62 62 62 ...
## $ STAT: num 81 81 81 68 63 73 68 56 70 45 ...

创建一个空网络，节点对应于marks的变量。然后通过指派一个两列的矩阵来添加边。

生成一个无向图：
ug <- empty.graph(names(marks))
arcs(ug, ignore.cycles = TRUE) = matrix(c("MECH", "VECT", "MECH", "ALG", "VECT", 
    "MECH", "VECT", "ALG", "ALG", "MECH", "ALG", "VECT", "ALG", "ANL", "ALG", 
    "STAT", "ANL", "ALG", "ANL", "STAT", "STAT", "ALG", "STAT", "ANL"), ncol = 2, 
    byrow = TRUE, dimnames = list(c(), c("from", "to")))
ug

## 
## Random/Generated Bayesian network
## 
## 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 
## 
## generation algorithm: Empty

这个ug对象属于bn类，这个类用于在bnlearn包中管理网络结构。 这个对象包括三个方面的信息：
 (1)learning:结构的学习 (2)node 节点 (3)arc 边

生成一个有向图：

dg <- empty.graph(names(marks))
arcs(dg) = matrix(c("VECT", "MECH", "ALG", "MECH", "ALG", "VECT", "ANL", "ALG", 
    "STAT", "ALG", "STAT", "ANL"), ncol = 2, byrow = TRUE, dimnames = list(c(), 
    c("from", "to")))
dg

## 
## Random/Generated Bayesian network
## 
## model:
## [STAT][ANL|STAT][ALG|ANL:STAT][VECT|ALG][MECH|VECT:ALG] 
## 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 
## 
## generation algorithm: Empty

有时候也可以从邻接矩阵（adjacency matrix）来生成有向图dg。

mat <- matrix(c(0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 
    0, 0, 0, 0, 0), nrow = 5, dimnames = list(nodes(dg), nodes(dg)))
mat

## MECH VECT ALG ANL STAT
## MECH 0 0 0 0 0
## VECT 1 0 0 0 0
## ALG 1 1 0 0 0
## ANL 0 0 1 0 0
## STAT 0 0 1 1 0

dg2 <- empty.graph(nodes(dg))
amat(dg2) <- mat
all.equal(dg, dg2)

## [1] TRUE

手工修改一个已经存在的网络，利用加边（set.arc）、去边（drop.arc）、颠倒（rev.arc）这几个操作，也可以得到一个需要的网络。

dg3 <- empty.graph(nodes(dg))
dg3 <- set.arc(dg3, "VECT", "MECH")
dg3 <- set.arc(dg3, "ALG", "MECH")
dg3 <- set.arc(dg3, "ALG", "VECT")
dg3 <- set.arc(dg3, "ANL", "ALG")
dg3 <- set.arc(dg3, "STAT", "ALG")
dg3 <- set.arc(dg3, "STAT", "ANL")
all.equal(dg, dg3)

## [1] TRUE

all.equal(ug, moral(dg)) #moral图

## [1] TRUE

我们全面希望了解网络的结构，可以使用存储在每个节点的信息。 
（1）节点的拓扑顺序

node.ordering(dg)

## [1] "STAT" "ANL" "ALG" "VECT" "MECH"

(2)节点的邻居(nbr)和Markov毯(mb) 
Markov毯是指对一个节点A，它的所有父节点，子节点以及和A有相同子节点的其它节点 。

nbr(dg, "STAT")

## [1] "ALG" "ANL"

mb(dg, "STAT")

## [1] "ALG" "ANL"

"ANL" %in% mb(dg, "STAT")

## [1] TRUE

"STAT" %in% mb(dg, "ANL")

## [1] TRUE

(3)某个给定节点的子代（child）和父代(parents)和子代的其它父代（o.par）

child <- children(dg, "STAT")
pa <- parents(dg, "STAT")



2.绘制网络

 bnlearn包有两种对bn对象绘制网络结构的方法，一种是利用graph包和Rgraphviz包提供的接口给出的一些绘图函数，比如graphviz.plot。再一种方法是利用plot函数。下面对于dg，来绘制网络图。 使用graphviz.plot的好处在于它返回一个graph对象，便于对网络的进一步操作。

library(Rgraphviz)

graphviz.plot(dg, layout = "fdp")

plot(dg, radius = 200, arrow = 30)