贝叶斯网络结构学习（基于BDAGL工具箱的MATLAB实现）

题目：贝叶斯网络结构学习（基于BDAGL工具箱的MATLAB实现）

        有关贝叶斯网络结构学习的一基本概念可以参考：贝叶斯网络结构学习方法简介

        有关函数输入输出参数的解释可以参考：贝叶斯网络结构学习若干问题解释

        注意到BDAGL工具箱比前面介绍的BNT工具箱在时间上还要早一些，但由于琢磨BDAGL一段时间后并没有找到任何头绪，后来只能另寻出路，所以才找了BNT工具箱。基本琢磨明白了BNT后，还是有些不甘心，所以回过头来继续琢磨BDAGL。在感觉没有希望想要放弃的时候，12月11号突然感觉找到了头绪。其实也并没有什么特别的事情，只是又看了看BDAGL的主页：

        http://www.cs.ubc.ca/~murphyk/Software/BDAGL/index.html，

包括子网页介绍的simpleDemo：

        http://www.cs.ubc.ca/~murphyk/Software/BDAGL/simpleDemo.htm

仅此而已，其实主页的例子已经说的很清楚了，那我们就开始吧。

        首先，按照主页的QuickStart流程走一遍：

        Download and unzip BDAGL.zip

        In Matlab, change the directory to the newly created BDAGL folder

        Run the script mkPath

        Type runDemos.

这个过程没有大问题就好，第3步运行mkPath可能会有问题，感觉安装一个Visual Studio也许就可以了吧，需要compiles Mex files，个人理解就是把C/C++写的文件编译成matlab可执行的文件。我后来直接将工具箱添加到了Matlab目录，方法参见《贝叶斯网络结构学习之K2算法（基于FullBNT-1.0.4的MATLAB实现）》的附录2。

        【注】若仅想使用该工具箱，此时可以调用文末的函数BDAGL_learn_struct_basic.m，此函数仅需输入数据集，输出即为学得的贝叶斯网络结构。

        接下来，打开目录\BDAGL\demos2内的cancerDemo.m，看这个例子不是关键，关键是通过这个例子学习一下工具箱中的函数如何使用。试着Run一下这个文件，应该是没问题的，若有问题可能是开始的步骤的原因。先不管运行结果，接下来解释一下这个文件所做的事情，后面自然就明白运行出的结果的含义了。

        第1步，得到aflp和aflml两个结果：

这两个结果的意思不用管，只知道接下来一直需要用它们就是了。

        第2步，关注四个结构学习结果：epDP、epMCMC、epEnumer、optimalDAG：

以上四个结果除了optimalDAG之外，其它均不可以直接使用，因为其它三个的结果是存储的是0~1之间的小数，运行结果的figure 1中的左上图(epDP)、右上图(epMCMC)、右下图(epEnumer)分别以不同的颜色显示了这三个结果。注意这三个结果是均是5*5的矩阵，因此图中四幅子图均为25个小格子，不同颜色对应0~1之间的值如子图右边的bar所示。左下图为真实网络结构，因为只有0和1两个值，所以图中只有两种颜色。

若想使用这三个结构学习结果，需要设置一个域值，类似例子如运行结果的Figure 2所示：

图中，左上图为真实网络结构，右上图为epDP域值为0的结构（即epDP矩阵中大于0的置为1，小于0的置为0），左下图为epDP域值为0.1的结构（即epDP矩阵中大于0.1的置为1，小于0.1的置为0），右下图为optimalDAG的结构。值得注意的是，右上图根本不是一个有向无环图，即根本不是一个符合要求的贝叶斯网；左下图的A和B、A和C之间有两条路径，颜色一深一淡，深的对应epDP的数值较大，淡的对应epDP的数值较小，必须分别舍掉其一方能符合有向无环图的要求。使用epMCMC、epEnumer类似，找一个合适的域值。

        简单介绍一下myDrawGraph函数：

        参数’labels’设置各节点的名称，’layout’设置各节点的坐标，在graphicalLayout矩阵中存储的是ABCDE五个节点的坐标：

A的坐标是(0.4,0.75)、B的坐标是(0.2,0.5)、C的坐标是(0.6,0.5)、D的坐标是(0.4,0.25)、E的坐标是(0.8,0.25)，坐标是将显示区域的横轴和纵轴都归一化到了0~1之间。myDrawGraph函数并不太好用，因为还要人为的设置各节点的坐标，当然可以使用默认设置，例如运行：

myDrawGraph(epDP, 'labels', labels,  'thresh', 0.1);

        输出结果如下：

        这一点在BDAGL的主页Graph layout部分也进行了说明：

并且作者推荐使用pajek，可以使用\BDAGL\demos2\adj2pajek2.m将学得的网络结构DAG转化为pajek需要的格式。

        那么什么是pajek呢？pajek是哪些单词的简称呢？打开百度翻译，将翻译语种设置为由斯洛文尼亚语到中文：

输入“pajek”：

        Pajek是大型复杂网络分析工具，是用于研究目前所存在的各种复杂非线性网络的有力工具。Pajek在Windows环境下运行，用于带上千乃至数百万个结点大型网络的分析和可视化操作。在斯洛文尼亚语中Pajek是蜘蛛的意思。（摘自百度百科词条“pajek”）

        具体参见链接：http://mrvar.fdv.uni-lj.si/pajek/

        使用说明中译本：《蜘蛛：社会网络分析技术》（林枫 译，世界图书出版公司）

        到此，似乎该说的都说完了，但应该还是缺点什么。对于普通人来说，我们需要的是给函数输入数据集，输出就是学得的贝叶斯网络结构。下面，我给出一个基于BDAGL工具箱的贝叶斯网络结构学习函数，可以达到这个要求： 

function [ DAG ] = BDAGL_learn_struct_basic( data_input )%调用工具箱BDAGL中的computeOptimalDag函数学习贝叶斯网络结构%该函数主要是对BDAGL进行二次封装%输入仅有data，输出即为学得贝叶斯网络结构DAG    data = data_input;    %节点个数    nNodes = size(data,1);    %初始化DAG    %DAG = zeros(nNodes,nNodes);    %得到每个节点在数据集中可能取值的最小值    data_low_limit = min(data,[],2);    %将数据集进行转化，使每个节点的可能取值从1开始    %这是函数mkAllFamilyLogMargLik的要求，否则会报错    for nn=1:nNodes        data(nn,:)=data(nn,:)+ (1-data_low_limit(nn));    end    %得到mkAllFamilyLogMargLik的输入参数nodeArity    %nodeArity: Number of values each node can take on (applicable only for    %multinomial case). **data for node i must be in range 1..nodeArity(i)**    nodeArity = max(data,[],2);    %调用mkAllFamilyLogMargLik得到aflml    %aflml(All Families Log Marginal Likelihood matrix)    aflml = mkAllFamilyLogMargLik(data,'nodeArity', nodeArity);    %计算最优网络结果： MLE, max p(D|G)    DAG = computeOptimalDag(aflml);end

        特别注意的是函数mkAllFamilyLogMargLik的输入参数nodeArity的设置，函数注释中明确要求“data for node i must be in range 1..nodeArity(i)”，试了几次设置不当的情况：

也不知究竟是何错误，竟然导致Matlab软件直接挂掉。因此，在BDAGL_learn_struct_basic内部，对data进行了处理，将其元素值转换到了1..nodeArity(i)范围内。

        函数BDAGL_learn_struct_basic仅简单使用了epDP、epMCMC、epEnumer、optimalDAG四种结果中的optimalDAG，若需使用其它结果，可自行改造。改造时需要注意的是首先要将epDP、epMCMC、epEnumer变成0/1二值，而且要保证结果是有向无环图（即先设置域值二值化，然后结合原值与二值化的结果进行二次处理保证其为有向无环图）。可以使用以下函数检测结果是否为有向无环图：

function NotDAG = IsDAG(DAG)%检测矩阵DAG是否表示一个有向无环图(Directed Acyclic Graph,DAG)%输出为1表示输入矩阵DAG不是一个有向无环图    tmp_DAG=DAG;    NotDAG = 0;    nNodes = size(DAG,1);    order=zeros(1,nNodes);        selected=[];        for i=1:nNodes        b=(sum(tmp_DAG)==0);        idx=find(b);        idx=setdiff(idx,selected);        if(isempty(idx))            NotDAG =1;            disp('The graph is not a directed acyclic graph, please check again.');            break;        else                        idx=idx(1);            order(1,idx)=i;        end        tmp_DAG(idx,:)=0;        selected=[selected,idx];    endend

若输出为0则表示是有向无环图，若输出为1则表示不是有向无环图。