1、贝叶斯网络基础

首先复习一下贝叶斯公式

例题:分别有 AB 两个容器,在容器 A 里分别有 7 个红球和 3 个白球,在容器 B 里有 1 个红球和 9 个白球,现已知从这两个容器里任意抽出了一个球,且是红球,问这个红球是来自容器 A 的概率是多少?

则有:P() = 8/20P(A) = 1/2P(|A) = 7/10,其中P()表示整体上摸出红球的概率,P(A)表示选中A容器的概率,P(|A)表示从A容器条件下摸出红球的概率

按照公式,则有:P(A|) = P(|A)*P(A) / P() = (7/10)*(1/2) / (8/20) = 0.875

1.1、朴素贝叶斯理论

朴素贝叶斯的思想基础是这样的:对于给出的待分类项,求解在此项出现的条件下各个类别出现的概率,哪个最大,就认为此待分类项属于哪个类别。 

1.2、朴素贝叶斯理论的局限性

朴素贝叶斯理论假设各个特征属性条件独立,但是实际情况中这比较难满足,问题随之也来了,由于特征属性间存在依赖关系(如头像是否真实与好友密度之间),使得朴素贝叶斯分类不适用了。既然这样,需要寻找另外的解决方案。也就是贝叶斯网络

1.3、贝叶斯网络

一个贝叶斯网络定义包括一个有向无环图(Directed AcyclicGraph )和一个条件概率表集合。

DAG中每一个节点表示一个随机变量,可以是可直接观测变量或隐藏变量,而有向边表示随机变量间的条件依赖;条件概率表中的每一个元素对应DAG中唯一的节点,存储此节点对于其所有直接前驱节点的联合条件概率。

2、贝叶斯网络MATLAB编程实现

2.1、环境搭建

安装MATLAB,添加FULLBNT工具箱具体参考这里(http://blog.sina.com.cn/s/blog_6c7b434d01013ufz.html

注意MATLAB中字体需要设置为中文字体,才能在编辑器和命令窗口支持中文,否则汉字会显示为小方框

2.2、习题:是窃贼还是地震

福尔摩斯先生在他的办公室工作时接到了他邻居华生的电话。华生告诉他:他的家里可能进了窃贼,因为他家的警铃响了

被告知有窃贼闯入,福尔摩斯迅速开车回家。在路上,他听广播得知他家那里发生了地震。地震也有可能引起警报。这样,请问福尔摩斯先生应该回家抓贼还是迅速撤离该地区以躲避地震?

条件如下

题目分析:

简单讲,在路上的holmes需要判断是盗贼还是地震导致警铃?如果是前者,他需要回去抓贼,若是后者,则要逃离地震区。

所以图中虽然有5个节点,地震并不100%导致警铃,警铃也不100%导致华生的信号。

但是我们在得到信号,听到警铃的情况下,可以通过计算盗贼导致警铃的概率p1,和地震导致警铃的概率pp1来进行决策,也可以计算在地震发生条件下,盗贼导致警铃的概率p2

如果p2p1小,说明新添加的条件E才是导致A的主要原因。

2.3、具体实现:

%1、建立贝叶斯网络结构
N = 3; %三个节点,分别是B、E、A
dag = zeros(N,N);
B = 1; E = 2; A = 3; 
%节点之间的连接关系
dag(B,A) = 1;
dag(E,A) = 1; 
discrete_nodes = 1:N; %离散节点
node_sizes = 2*ones(1,N);%节点状态数
bnet =mk_bnet(dag,node_sizes,'names',{ 
    'BB','EE','AAA'
    },'discrete',discrete_nodes);

bnet.CPD{B} = tabular_CPD(bnet,B,[0.9 0.1]);%手动输入的条件概率
bnet.CPD{E} = tabular_CPD(bnet,E,[0.99 0.01]);
bnet.CPD{A} = tabular_CPD(bnet,A,[0.99 0.1 0.1 0.01 0.01 0.9 0.9 0.99]);
%2、画出建好的贝叶斯结构

draw_graph(dag);

%3、使用联合树引擎对贝叶斯网络进行推断
engine = jtree_inf_engine(bnet);
%4、求解边缘分布假设,
%我们要计算盗窃导致响铃的概率
evidence = cell(1,N);
evidence{A} = 2;

[engine, loglik] = enter_evidence(engine, evidence);
marg = marginal_nodes(engine, B);
p1 = marg.T(2);%算出p1=0.8412
%现在我们添加地震的证据观察它有什么不同
evidence{E} = 2;
[engine, loglik] = enter_evidence(engine, evidence);
marg = marginal_nodes(engine, B);
p2 = marg.T(2);%算出p2=0.1089
%结论是地震更能解释响铃这个主要事实
%联合概率分布
evidence = cell(1,N);
[engine, ll] = enter_evidence(engine, evidence);
m = marginal_nodes(engine, [B E A]);

1.     建好的贝叶斯网络

接下来添加RW节点

2.4、修改上述代码

%1、建立贝叶斯网络结构
N = 5; %三个节点,分别是B、E、A、R、W
%分别代表盗贼、地震、警铃、广播、华生致电holmes;
dag = zeros(N,N);
B = 1; E = 2; A = 3; 
R = 4; W=5;
%节点之间的连接关系
dag(B,A) = 1;
dag(E,A) = 1; 
dag(E,R) = 1;
dag(A,W) = 1;
discrete_nodes = 1:N; %离散节点
node_sizes = 2*ones(1,N);%节点状态数
bnet =mk_bnet(dag,node_sizes,'names',{ 
    ' BB','EE','AAA','RR','WWW'
    },'discrete',discrete_nodes);
bnet.CPD{B} = tabular_CPD(bnet,B,[0.9 0.1]);%手动输入的条件概率
bnet.CPD{E} = tabular_CPD(bnet,E,[0.99 0.01]);
bnet.CPD{A} = tabular_CPD(bnet,A,[0.99 0.1 0.1 0.01 0.01 0.9 0.9 0.99]);
bnet.CPD{R} = tabular_CPD(bnet,R,[0.999 0.01 0.001 0.99]);
%概率表输入顺序:本节点状态不变,条件变化……
bnet.CPD{W} = tabular_CPD(bnet,W,[0.99 0.35 0.01 0.65]);
%2、画出建好的贝叶斯结构
draw_graph(dag);
%3、使用联合树引擎对贝叶斯网络进行推断
engine = jtree_inf_engine(bnet);
%4、求解边缘分布假设,
%我们要计算盗窃导致响铃的概率
evidence = cell(1,N);
evidence{A} = 2;
[engine, loglik] = enter_evidence(engine, evidence);
marg = marginal_nodes(engine, B);
p1 = marg.T(2);%算出p1=0.8412
%现在我们添加地震的证据观察它有什么不同
evidence{E} = 2;
[engine, loglik] = enter_evidence(engine, evidence);
marg = marginal_nodes(engine, B);
p2 = marg.T(2);%算出p2=0.1089
%结论是地震更能解释响铃这个主要事实
%联合概率分布
evidence = cell(1,N);
[engine, ll] = enter_evidence(engine, evidence);
m = marginal_nodes(engine, [B E A]);

运行结果也可看出RW节点对决策无影响

3BNT参数、结构学习

Bayesian 网络学习指的是通过分析数据儿获得Bayesian网的过程,它包括以下两种情况

①参数学习:已知网络结构,确定网络参数。

②结构学习:既要确定网络结构,又要确定网络参数。

我们可以一个贝叶斯网络看作一个黑箱,那参数学习,结构学习就是创造一个白箱{可能是我造的名词}尽可能模拟这个黑箱!

在前面的基础之上,我们可以轻松画出贝叶斯网络bnet{!当然,其每个节点的CPT都是题目给出的}。经过bnet处理之后的样本数据会输出得到训练数据data

在①情况下,构建一个和bnet节点个数,箭头指向相同的bnet1,接着借助BNT中的工具,将databnet1作为参数,就可以学习得到新的贝叶斯网络bnet2

在②情况下,借助BNT中的工具,将data和其他数据{!具体参数见代码}作为参数,便可学习得到新的贝叶斯网络bnet3

接下来使用如下语句便可以方便的输出一个贝叶斯网络中各个节点的CPT{!第一节输出CPT的方法容易理解但有些繁琐}

%输出bnet的CPT
CPT = cell(1,N);
for i=1:N
  s=struct(bnet.CPD{i});
  CPT{i}=s.CPT;
end
dispcpt(CPT{1})
dispcpt(CPT{2})
dispcpt(CPT{3})
dispcpt(CPT{4})
dispcpt(CPT{5})

详见代码

%这里使用的例子依然是“福尔摩斯模型”。
%首先我先如上面实验那样建立好贝叶斯网bnet,并手动构造条件概率表CPT。
%然后使用BNT里的函数sample_bnet(bnet)来产生nsamples个数据样本,
%nsamples分别取值20,200,2000。
%然后,再重新建立一个不知道条件概率表的贝叶斯网bnet2(结构和bnet相同),
%并把得到的样本作为训练集代入learn_params()函数进行学习,
%把学习到的条件概率表CPT2与手动构造的CPT进行了比较。参数学习部分代码:
 
clc;
clear;
%1、构造样本数据data
nsamples = 2000;
samples = cell(N, nsamples);
for i = 1:nsamples
    samples(:,i) = sample_bnet(bnet);
end
data = cell2num(samples);
 
 
%2、手动构造和bnet结构相同的bnet1
bnet1 = mk_bnet(dag,node_sizes,'discrete',discrete_nodes);
seed = 0;
rand('state',seed);
bnet1.CPD{B} = tabular_CPD(bnet2,B);
bnet1.CPD{E} = tabular_CPD(bnet2,E);
bnet1.CPD{W} = tabular_CPD(bnet2,W);
bnet1.CPD{R} = tabular_CPD(bnet2,R);
bnet1.CPD{A} = tabular_CPD(bnet2,A);
%参数学习得到新的贝叶斯网络bnet2
bnet2 = learn_params(bnet1,data);
 
%3、结构学习
order=[1 2 3 4 5];
ns=[2 2 2 2 2];
max_fan_in=2;
%结构学习函数
dag2 = learn_struct_K2(data ,ns, order,'max_fan_in', max_fan_in);
figure
draw_graph(dag2);
bnet3=mk_bnet(dag2,ns,'discrete',order);

利用输出语句CPT的语句可以得到

原贝叶斯网络bnetCPT

1 : 0.9000 
2 : 0.1000 
1 : 0.9900 
2 : 0.0100 
1 1 : 0.9900 0.0100 
2 1 : 0.1000 0.9000 
1 2 : 0.1000 0.9000 
2 2 : 0.0100 0.9900 
1 : 0.9990 0.0010 
2 : 0.0100 0.9900 
1 : 0.9900 0.0100 
2 : 0.3500 0.6500 

新贝叶斯网络bnet2CPT

1 : 0.9030 
2 : 0.0970 
1 : 0.9890 
2 : 0.0110 
1 1 : 0.9888 0.0112 
2 1 : 0.0990 0.9010 

1 2 : 0.1500 0.8500 
2 2 : 0.0000 1.0000 
1 : 0.9995 0.0005 
2 : 0.0000 1.0000 
1 : 0.9877 0.0123 
2 : 0.3302 0.6698 

和(http://my.oschina.net/SnifferApache/blog/343756#OSC_h2_7)所示的条件概率表相比,相差无几

新贝叶斯网络bnet3CPT

Referenceto non-existent field 'CPT'.

shit,只能用第一节中的旧方法了……

engine2 = jtree_inf_engine(bnet3);
%新证据
evidence2 = cell(1,N);
[engine2, loglike] = enter_evidence(engine2, evidence2);
 
%计算节点B的CPT
marg2 = marginal_nodes(engine2,B);
marg2.T   
%计算节点E的CPT
marg2 = marginal_nodes(engine2,E);
marg2.T 
%计算节点W的CPT
marg2 = marginal_nodes(engine2,[W A]);
marg2.T 
%计算节点R的CPT
marg2 = marginal_nodes(engine2,[R E]);
marg2.T 
%计算节点A的CPT
marg2 = marginal_nodes(engine2,[A B E]);
marg2.T

输出如下

Error usingjtree_inf_engine/enter_evidence (line 55)

must define CPD 1

事实证明结构学习得到的bnet3不包含CPT字段,导致无法计算CPT,此处挖坑待高人解答吧。

{!补充:BNT中参数学习函数:

最大似然性估计learn_params() ;

贝叶斯方法bayes_update_params()

结构学习函数:

•K2算法learn_struct_K2() ;

贪婪搜索GS(GreedySearch) 算法learn_struct_gs()

爬山HC ( HillClimbing) 算法learn_struct_hc()

BNT提供的多种推理引擎:

联合树推理引擎jtree_inf_engine();

全局联合树推理引擎global_joint_inf_engine();

.信念传播推理引擎belprop_inf_engine();

变量消元推理引擎var_elim_inf_engine()

}

4、两道人工智能练习题

4.1、试采用matlab的贝叶斯工具箱BNT完成以下过劳死问题的建模。要求代码中使用BNT工具箱中的贝叶斯结构、参数学习命令。给出网络学习代码及运行结果。

仅供参考

clear;
clc
N=5;
%构造网络结构
dag=zeros(N,N);
C=1;U=2;W=3;B=4;D=5;
dag(C,U)=1;
dag(U,[W,B])=1;
dag([W,B],D)=1;
discrete_nodes=1:N;
node_sizes=2*ones(1,N);
bnet=mk_bnet(dag,node_sizes,'names',{'country','university','workhard','badbody','die'},'discrete',discrete_nodes);
%构造概率表
bnet.CPD{C}= tabular_CPD(bnet,C,[0.5,0.5]);
bnet.CPD{U}= tabular_CPD(bnet,U,[0.99,0.05,0.01,0.95]);
bnet.CPD{W}= tabular_CPD(bnet,W,[0.95,0.10,0.05,0.90]);
bnet.CPD{B}= tabular_CPD(bnet,B,[0.99,0.70,0.01,0.30]);
bnet.CPD{D}= tabular_CPD(bnet,D,[1.00,0.95,0.70,0.665,0.00,0.05,0.30,0.335]);
figure 
draw_graph(dag);
 
%输出bnet的CPT
CPT = cell(1,N);
for i=1:N
  s=struct(bnet.CPD{i});
  CPT{i}=s.CPT;
end
dispcpt(CPT{1})
dispcpt(CPT{2})
dispcpt(CPT{3})
dispcpt(CPT{4})
dispcpt(CPT{5})
%输出结束
 
engine=jtree_inf_engine(bnet);
evidence=cell(1,N);
engine=enter_evidence(engine,evidence);
%计算节点C的CPT
marg=marginal_nodes(engine,C);
marg.T
%计算节点D的CPT
evidence{W}=2;
evidence{B}=2;
[engine,loglik]=enter_evidence(engine,evidence);
marg1=marginal_nodes(engine,[W B D]);
marg1.T(2)
 
%构造样本数据
nsamples=2000;
samples=cell(N,nsamples);
for i=1:nsamples
    samples(:,i)=sample_bnet(bnet);
end
data=cell2num(samples);
bnet1=mk_bnet(dag,node_sizes,'discrete',discrete_nodes);
%手动构造条件概率表CPT
seed=0;
rand('state',seed);
bnet1.CPD{C}=tabular_CPD(bnet1,C);
bnet1.CPD{U}=tabular_CPD(bnet1,U);
bnet1.CPD{W}=tabular_CPD(bnet1,W);
bnet1.CPD{B}=tabular_CPD(bnet1,B);
bnet1.CPD{D}=tabular_CPD(bnet1,D);
%进行参数学习
bnet2=learn_params(bnet1,data);
%验证学习结果
engine2=jtree_inf_engine(bnet2);
evidence2=cell(1,N);
engine2=enter_evidence(engine2,evidence2);
%计算节点C的CPT
marg=marginal_nodes(engine2,C);
marg.T
%计算节点D的CPT
evidence2{W}=2;
evidence2{B}=2;
[engine2,ll]=enter_evidence(engine2,evidence2);
m=marginal_nodes(engine2,[W B D]);
m.T(2) 
 
%结构学习
order=[1 2 3 4 5];
ns=[2 2 2 2 2];
max_fan_in=2;
%结构学习函数
dag2 = learn_struct_K2(data ,ns, order,'max_fan_in', max_fan_in);
figure
draw_graph(dag2);
bnet3=mk_bnet(dag2,ns,'discrete',order);

4.2、试采用核函数方法对下列数据进行非线性分类。给出matlab具体代码,及采用训练样本进行测试得到的准确率结果。     

       x=[0 1 0 1 2 -1];

       y=[0 0 1 1 2 -1];

       z=[-1 1 1 -1 1 1;

      其中,(x,y)代表二维的数据点,z 表示相应点的类型属性。

仅供参考

clc;
clear;
%手工录入+1类和-1类的数据
sp = [1,0; 0,1; 2,2; -1,-1] %positive
sn = [0,0; 1,1;]%negative
data = [sp;sn]
datalabel=[1 1 1 1 0 0]; %分类
%groups为logical类型
groups = ismember(datalabel,1);
%交叉标记data为train和test两部分
[train,test]=crossvalind('holdOut',groups);
%核函数 mlp,预测准确率为0.3
svmStruct=svmtrain(data(train,:),groups(train),'kernel_function','mlp','showplot',true);
%预测准确率为0
%svmStruct2=svmtrain(data,groups,'kernel_function','mlp','showplot',true);
%加入参数[0.5 -0.5]之后,预测准确率为0.6
svmStruct3 = svmtrain(data,groups,'kernel_function','mlp','mlp_params',[0.5 -0.5],'showplot',true);
 
title(sprintf('核函数'))
cp=classperf(groups);%用来评价分类器
classes = svmclassify(svmStruct3,data(test,:),'showplot',true);
classperf(cp,classes,test);
cp.CorrectRate
5、参考资料

[1]贝叶斯网络-使用MATLAB工具集,http://www.cnblogs.com/xjx-user/archive/2013/04/14/3020273.html

6、贝叶斯网络简述

* Bayes Net最初的目的是给一个专家系统增加概率。

6.1网络构造

在实际中,变量之间可能存在依赖关系。贝叶斯网络提供一种因果关系的图形模型,训练后的贝叶斯网络可以用于分类。

Bayes Net由两部分组成,有向无环图(DAG,directed acyclic graph)和条件概率表(CPT,conditionalprobability distribution tabular)

Dag中每个节点表示一个随机变量,变量可以是离散的也可以连续,每条边表示一个概率依赖。如A---B有一条边,则ABparentsBA的后代。给定parents,每个变量条件独立于图中它的非后代。

对每个变量,Bayes Net提供一个条件概率表(CPT)。

6.2训练网络

Bayes Net训练包括两部分,网络拓扑和网络参数

网络拓扑可以人工构造,或数据导出。

网络变量可以使可观测的或隐藏的。(可观测意思是可以直接给出变量值)

如果网络拓扑一直,变量可观测,则训练网络参数就是计算CPT表目。与朴素贝叶斯分类的概率计算类似。

## How to use the Bayes Network Toolbox

定义一个Bayes网络,要确定两点,网络结构和参数。

教程地址:http://www.cs.ubc.ca/~murphyk/Software/BNT/usage.html#inference

1.构造Bayes Net

Bayes Net是一个有向非循环图(directed acyclicgraph-dag),我们用邻接矩阵来表示一个dag

* 对节点标记顺序,父节点序号在子节点之前。

* 如果节点的取值是离散的,那么节点的大小就是总共可以取值的数目

       discrete_nodes=1:N; %标注那些节点是离散的,这里就是1:N的都是离散的

       node_sizes=2*ones(1,N); %node_sizes(i) 是节点 i 可以取得值的数目,在这个例子中,每个节点都只有两个取值

* 构造BN可以用如下命令:

        bnet=mk_bnet(dag,node_sizes)

* 如果可以确定哪个变量(节点)是可以观察到的,使用下面的代码。如果没有,可以空缺

        onodes=[];

       bnet=mk_bnet(dag,node_sizes,'discrete',discrete_nodes,onodes);

* 可选参数

        'names' 一个cell向量,对于每个节点的名字

       bnet=mk_bnet(dag,node_sizes,'names',{'A','B','C'},'discrete',1:3);

2. 网络拓扑结构学习

1. 中,我们假设已知网络的拓扑结构,因此可以直接构造网络。但如果我们不知道网络的拓扑,那么就需要对拓扑结构进行学习。

* 在网络拓扑结构学习中,有两个主要的路线

1. 基于约束条件的学习:将网络全连接,然后去掉条件独立的边

2. 搜索打分学习:在可能的DAGs空间中搜索,对每个进行打分

* 但是要注意,可能的DAG非常之多,是随节点数目指数级增长的。

* 如果节点数比较小(小于等于6个)可以用暴力搜索,因为6个节点的DAG共有3781503个(7个节点的话就是10^9数量级个)。

       dags=mk_all_dags(N);

       score=score_dags(data,ns,dags);

        % data(i,m)是节点i在第mcase中的取值 ns(i)是节点i的大小

K2 algorithm

如果我们知道所有节点的顺序,那么就可以独立的给每个节点选择它们的父节点

K2是一种贪婪搜索算法。一开始每个节点都没有父节点,然后逐渐添加父节点,每个新添加的都是使得当前分数最大增加的。

       dag=learn_struct_K2(data,node_sizes,order)

        其中 order 是节点的顺序信息,父节点要在子节点前

3.Bayes Net的参数

参数用CPDobjects(CPD=Conditional Probability Distribution)表示。CPD表示给定父节点时,该节点的概率分布。当所有节点的取值都是离散的时候,CPD可以用一个table来表示。即为 CPTs (Tabular CPDs)

* Tabular CPDs(扁平化的CPD)成为CPTs

CPT是一个多维向量,比如节点A有两个父节点C,D,那么ACPT就是4个维度,第三个维度为A的取值可能,第四个维度是每种取值的概率。

* 工具箱中默认定义:false==1,true==2,在对参数赋值中应该保持一致性。

如果已知CPD,构造CPT的函数是

       bnet.CPD{A}=tabular_CPD(bnet,A,[a1,a2,a3,a4...]);

4.CPD

CPD的定义是:P(X(i)|X(Pa(i)))X(Pa(i))是节点i的父节点,CPDnode i的父节点发生的情况下i发生的概率。

学习CPD

已有12个参数的BN,网络拓扑结构dag

如果我们不指定CPT,那么算法会构造随机参数,每一列CPT会从均匀分布中提取。

       

 bnet=mk_bnet(dag,node_sizes);
        seed=0;%为了保证实验的可重复性
        rand('state',seed);
        for i=1:12
          bnet.CPD{i}=tabular_CPD(bnet,i);
        end
        bnet2=learn_params(bnet,data);%参数学习完毕
        for i=1:12
            s=struct(bnet2.CPD{i});%展示参数
            CPT3{i}=s.CPT;
        end
        dispcpt(CPT3{1})%展示第一个node的参数

5.Inference推理

BNT 有许多不同的Inference"engines"比如联合树引擎(junction tree engine)。可以通过如下方式使用:

        engine =jtree_inf_engine(bnet);

举个例子,在之后的讨论中也将沿用这个例子

    脑卒中数据挖掘中,颈部血管筛查的结果,总结出7个节点。按照顺序(父节点在子节点之前)列如下:内膜增厚(TMT),斑块数目,斑块形态,斑块表面,溃疡,回声,是否脑卒中(dfStroke)。

    我们已经构造出网络拓扑,bnet,并且学习了参数。下面使用它进行推断

假设我们知道斑块数量,斑块形形态以及是否溃疡三个证据,我们想知道被检查者患有脑卒中的概率(也是计算边缘分布),将使用到BayesNet中的Inference

代码如下:

engine=jtree_inf_engine(bnet2);%使用联合树搜索引擎
evidence=cell(1,7);%由于有7个可观测节点(包含是否脑卒中节点)
evidence{2}=2;
evidence{3}=1;
evidence{4}=2;%分别给出观测结果.为赋值的evidence是隐藏的,即evidence{i}=[]
[engine,loglik]=enter_evidence(engine,evidence);%The first return argument contains the modified engine, which incorporates theevidence. The second return argument contains the log-likelihood of theevidence
marg=marginal_nodes(engine,7)%计算第7个节点(dfStroke)的边缘分布

 

Logo

权威|前沿|技术|干货|国内首个API全生命周期开发者社区

更多推荐