% ********************************************************************* %
% ************************ 贝叶斯网络实例 ***************************** %
% ********************************************************************* %
%
clear
clc
% *********************************************************************** %
% ……………………………………… 已知网络 ……………………………………………… %
% *********************************************************************** %
N = 5; % 4个节点
dag = zeros(N,N); % dag为网络结构边的关系
E = 1; B = 2; R = 3; A = 4; C = 5; % 节点编号
dag(E,[R A]) = 1; % 节点加边
dag(B,A) = 1; % 节点加边
dag(A,C)=1; % 节点加边
false = 1; true = 2;
ns = 2*ones(1,N); % 节点可能的取值个数（这里为2个，取值为：0或1）
bnet = mk_bnet(dag, ns); % 生成bnet，包括结构和参数（概率）
bnet.CPD{E} = tabular_CPD(bnet, E, [0.1 0.9]);% 给各个节点添加先验概率
bnet.CPD{B} = tabular_CPD(bnet, B, [0.01 0.99]);
bnet.CPD{R} = tabular_CPD(bnet, R, [0.65 0.01 0.35 0.99]);
bnet.CPD{A} = tabular_CPD(bnet, A, [0.95 0.8 0.3 0.001 0.05 0.2 0.7 0.999]);
bnet.CPD{C} = tabular_CPD(bnet, C, [0.7 0.05 0.3 0.95]);
% *********************************************************************** %
% ……………………………………… 生成数据 ……………………………………………… %
% *********************************************************************** %
seed = 0; %seed表示采用v4版本的随机数产生器，state是v5版本的随机数产生器
rand('state', seed); % 使用state指定状态，这样产生的随机结果相同，
randn('state', seed); % 也就是说：你我生成的随机数相同，不论时间和地点
ncases = 1000;
data = zeros(N, ncases);
for m=1:ncases
data(:,m) = cell2num(sample_bnet(bnet)); % 根据以上生成的bnet生成随机样本
end
% *********************************************************************** %
% ………………………………………1.结构学习……………………………………………… %
% *********************************************************************** %
order = [2,1,4,5,3]; % 定义节点序，节点序属于先验知识
max_fan_in = 2; % 定义最大父节点个数，不超过总节点个数N
% 结构学习（采用K2算法进行结构训练）
dag2 = learn_struct_K2(data, ns, order, 'max_fan_in', max_fan_in, 'verbose', 'yes');
% [dags,~,~] = learn_struct_mcmc(data, ns);
% Sum_dag = zeros(N,N);
% for i=1:size(dags,2)
% Sum_dag = Sum_dag+dags{1,i};
% end
% dag2 = triu(Sum_dag);
% for i=1:N
% for j=1:N
% if dag2(i,j)>300
% dag2(i,j) = 1;
% else
% dag2(i,j) = 0;
% end
% end
% end
% *********************************************************************** %
% ………………………………………2.参数学习……………………………………………… %
% *********************************************************************** %
bnet2 = mk_bnet(dag2, ns);
seed = 0;
rand('state', seed);
bnet2.CPD{E} = tabular_CPD(bnet2, E, 'clamped', 1, 'CPT', [0.1 0.9], 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
bnet2.CPD{B} = tabular_CPD(bnet2, B, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
bnet2.CPD{R} = tabular_CPD(bnet2, R, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
bnet2.CPD{A} = tabular_CPD(bnet2, A, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
bnet2.CPD{C} = tabular_CPD(bnet2, C, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
% 从完全观测数据找最大似然估计(maximumlikelihood estimator,MLE)
bnet3 = learn_params(bnet2, data);
% *********************************************************************** %
% ……………………………………………3.推理……………………………………………… %
% *********************************************************************** %
engine = jtree_inf_engine(bnet3); % 联合树推理引擎
ev = cell(1,N); % 添加证据
ev{C} = 1; % 节点C出现
ev{R} = 1; % 节点R出现
engine = enter_evidence(engine, ev); % 向网络添加指定的证据
mE = marginal_nodes(engine, E); % 计算节点E的边缘概率
mB = marginal_nodes(engine, B); % 计算节点B的边缘概率
% *********************************************************************** %
% …………………………………………… 结果显示 ………………………………………… %
% *********************************************************************** %
% （1）显示结构学习效果
h1 = view(biograph( dag )); % 画出已知网络结构图
h2 = view(biograph( dag2 )); % 画出学习网络结构图
% （2）显示参数学习效果
CPT3 = cell(1,N);
for i=1:N
s=struct(bnet.CPD{i});
CPT{i}=s.CPT;
s1=struct(bnet3.CPD{i});
CPT3{i}=s1.CPT;
end
AA = [CPT{1,4}(:,:,1) CPT{1,4}(:,:,2)]; % 取了其中一个节点比较（节点A）
BB = [CPT3{1,4}(:,:,1) CPT3{1,4}(:,:,2)];
data = [[AA(1,:),AA(2,:)]' [BB(1,:),BB(2,:)]'];
colnames = {'原始概率' '学习概率'};
figure;
uitable('position',[0 80 200 200],'data',data,'ColumnName',colnames)
% （3）显示推理结果
msgbox(['P(E|C,R)=',num2str(mE.T(1)),'；P(B|C,R)=',num2str(mB.T(1)),' www.ttin.top'])
