Bayes Net toolbox 的使用实例

Example 1、离散贝叶斯网络模型

•  Cloudy, Sprinklet, Rain, WetGrass模型，每个节点的条件（先验）概率分布如图所示

%% ------------------ 模型设置 ------------------%%

% model structure

N = 4;

dag = zeros(N,N);

C = 1; S = 2; R = 3; W = 4;

dag(C,[R S]) = 1;

dag(R,W) = 1;

dag(S,W)=1;

% Make a bnet

discrete_nodes = 1:N;

node_sizes = 2*ones(1,N);

bnet = mk_bnet(dag, node_sizes, 'names', {'C','S','R','W'}, 'discrete', discrete_nodes);

% parameters setting, define CPDs

bnet.CPD{C} = tabular_CPD(bnet, C, [0.5 0.5]);

bnet.CPD{R} = tabular_CPD(bnet, R, [0.8 0.2 0.2 0.8]);

bnet.CPD{S} = tabular_CPD(bnet, S, [0.5 0.9 0.5 0.1]);

bnet.CPD{W} = tabular_CPD(bnet, W, [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);

% rearrange CPT

CPT = reshape([1 0.1 0.1 0.01 0 0.9 0.9 0.99], [2 2 2]);

%% ------------------- 估计模型参数 -------------------%%

% randomly generate samples

nsamples = 100;

samples = cell(N, nsamples);

for i=1:nsamples

  samples(:,i) = sample_bnet(bnet);

end

% Make a random bnet

bnet2 = mk_bnet(dag, node_sizes);

seed = 0;

rand('state', seed);  % uniform distribution

bnet2.CPD{C} = tabular_CPD(bnet2, C);% tabular_CPD(bnet, node) 产生随机CPT

bnet2.CPD{R} = tabular_CPD(bnet2, R);

bnet2.CPD{S} = tabular_CPD(bnet2, S);

bnet2.CPD{W} = tabular_CPD(bnet2, W);

% parameter estimation using Maximum Likelihood (ML) Estimation

% samples中包括所有节点的观测数据（不包含缺失数据）

bnet3 = learn_params(bnet2, samples);

% show estimated parameter

CPT3 = cell(1,N);

for i=1:N

  s=struct(bnet3.CPD{i});  % violate object privacy

  CPT3{i}=s.CPT;

end

disp('Display estimated CPT parameters...');

dispcpt(CPT3{4})

% real parameter shown

disp('Display real CPT parameters...');

dispcpt(CPT)

% 1 1 : 1.0000 0.0000

% 2 1 : 0.1000 0.9000

% 1 2 : 0.1000 0.9000

% 2 2 : 0.0100 0.9900

%% ------------------- 缺失部分数据时的模型估计 -------------------%%

% hide half the values

samples2 = samples;

hide = rand(N,nsamples) > 0.5;    % 将随机生成值大于0.5的节点全部隐藏

[I,J] = find(hide);

for k=1:length(I)

    samples2{I(k),J(k)} = [];   % hide means set to []

end

% parameter estimation using Maximum Likelihood (ML) Estimation

% samples中包括缺失数据

max_iter = 10;

engine2 = jtree_inf_engine(bnet2);

[bnet4, LLtrace] = learn_params_em(engine2, samples, max_iter);

% plot the loglikelihood at iteration i;

plot(LLtrace,'x-');

% show estimated parameters

CPT4 = cell(1,N);

for i=1:N

  s=struct(bnet4.CPD{i});  % violate object privacy

  CPT4{i}=s.CPT;

end

celldisp(CPT4);

Example 2、HMM模型

•  HMM模型可以用2-slice BN来建模，模型如图所示，假设有2个hidde state，2个离散的观测符号
• 示例文件：..\bnt\BNT\examples\dynamic\dhmm1.m

%% ------------------ 模型设置 ------------------%%

% intra-slice BN

intra = zeros(2); 

intra(1,2) = 1;

% inter-slice BN

inter = zeros(2);

inter(1,1) = 1;

Q = 2;         % num hidden states

O = 2;         % num observable symbols

ns = [Q O];    % number of states

dnodes = 1:2;

%onodes = [1:2]; % only possible with jtree, not hmm

onodes = [2];

% 声明等价类

eclass1 = [1 2];    % 'eclass1' - equiv class for slice 1 [1:n]

eclass2 = [3 2];    % 'eclass2' - equiv class for slice 2 [tie nodes with equivalent parents to slice 1]

eclass = [eclass1 eclass2];     % eclass(4) == eclass(2) = 2，因此只有3个节点需要指定CPD

% 新建一个指定结构的DBN

bnet = mk_dbn(intra, inter, ns, 'discrete', dnodes, 'observed', onodes, 'eclass1', eclass1, 'eclass2', eclass2); 

%% 各个节点的CPD（real parameters）

prior1 = normalise(rand(Q,1))

transmat1 = mk_stochastic(rand(Q,Q))

obsmat1 = mk_stochastic(rand(Q,O))

bnet.CPD{1} = tabular_CPD(bnet, 1, prior1);

bnet.CPD{2} = tabular_CPD(bnet, 2, obsmat1);

bnet.CPD{3} = tabular_CPD(bnet, 3, transmat1);

% 建立一个相同结构的bnet ，用于模型预测

bnet1 = mk_dbn(intra, inter, ns, 'discrete', dnodes, 'observed', onodes, 'eclass1', eclass1, 'eclass2', eclass2); 

for i=1:3

    % 随机生成CPD初始值

  bnet1.CPD{i} = tabular_CPD(bnet1, i);

end

% 推荐的两种推理算法

% 要求所有隐藏节点都是离散的，并且在使用hmm_2TBN_inf_engine前需要指明观测节点

% engine = smoother_engine(hmm_2TBN_inf_engine(bnet));

engine = smoother_engine(jtree_2TBN_inf_engine(bnet1));

%% ---------------- 随机生成训练样本 ------------------- %%

ss = 2;            %slice size(ss)

ncases = 10;    %number of examples

T=10;

cases = cell(1, ncases); 

for i=1:ncases

  ev = sample_dbn(bnet, T);

  cases{i} = cell(ss,T);

  cases{i}(onodes,:) = ev(onodes, :);

end

%% ------------------- EM算法进行模型训练 -------------------------%%

max_iter=10;      %iterations for EM 

[bnet2, LLtrace] = learn_params_dbn_em(engine, cases, 'max_iter', max_iter);

% show estimated parameters

CPT = cell(1,3);

for i=1:3

  s=struct(bnet2.CPD{i});  % violate object privacy

  CPT{i}=s.CPT;

end

celldisp(CPT);

%% ------------------- State 预测（Viterbi） -------------------------%%

% Create a sequence

T = 10;

ev = sample_dbn(bnet, T);

evidence = cell(2,T);

evidence(2,:) = ev(2,:);     % extract observed component

data = cell2num(ev(2,:));

% Viterbi算法得到的路径

obslik = multinomial_prob(data, obsmat1);

path = viterbi_path(prior1, transmat1, obslik);  

% Most Probable Explanation

engine = {};

engine{end+1} = smoother_engine(jtree_2TBN_inf_engine(bnet));

mpe = find_mpe(engine{1}, evidence);

assert(isequal(cell2num(mpe(1,:)), path))     % extract values of hidden nodes

Example 3、GMM-HMM模型

•  GMM-HMM模型可以用2TBN来建模，模型如图所示
• 假设有2个hidde state，2个离散的观测符号
• 示例文件：..\bnt\BNT\examples\dynamic\mhmm1.m

% continuous HMM learning

% Make an HMM with mixture of Gaussian observations

%    Q1 ---> Q2

%  /  |   /  |

% M1  |  M2  |

%  \  v   \  v

%    Y1     Y2

% where Pr(m=j|q=i) is a multinomial and Pr(y|m,q) is a Gaussian    

% 第1个time slice有3个节点

intra = zeros(3);

intra(1,[2 3]) = 1;

intra(2,3) = 1;

% 第1个和第2个time slice之间的节点连接关系

inter = zeros(3);

inter(1,1) = 1;

n = 3;

% GMM-HMM topology

Q = 2; % num hidden states

O = 2; % size of observed vector

M = 2; % num mixture components per state

% Node_Sizes

ns = [Q M O];

% discete nodes

dnodes = [1 2];

% observable nodes

onodes = [3];

% 等价类（按照topology排序）

eclass1 = [1 2 3];

eclass2 = [4 2 3];

% 建立指定结构的HMM

bnet = mk_dbn(intra, inter, ns, 'discrete', dnodes, 'eclass1', eclass1, 'eclass2', eclass2, ...

       'observed', onodes);

% GMM-HMM 的已知分布      

prior0 = normalise(rand(Q,1));

transmat0 = mk_stochastic(rand(Q,Q));

mixmat0 = mk_stochastic(rand(Q,M));

mu0 = rand(O,Q,M);

Sigma0 = repmat(eye(O), [1 1 Q M]); % 对角Cov

bnet.CPD{1} = tabular_CPD(bnet, 1, prior0);

bnet.CPD{2} = tabular_CPD(bnet, 2, mixmat0);    % M 节点存储的是GMM的weight

bnet.CPD{3} = gaussian_CPD(bnet, 3, 'mean', mu0, 'cov', Sigma0);

bnet.CPD{4} = tabular_CPD(bnet, 4, transmat0);

%% ----------------------- 建立相同结构的DBN，用于预测HMM ------------------------- %%

bnet1 = mk_dbn(intra, inter, ns, 'discrete', dnodes, 'eclass1', eclass1, 'eclass2', eclass2, ...

       'observed', onodes);

bnet1.CPD{1} = tabular_CPD(bnet1, 1);

bnet1.CPD{2} = tabular_CPD(bnet1, 2);

bnet1.CPD{3} = gaussian_CPD(bnet1, 3, 'cov_type', 'diag');

bnet1.CPD{4} = tabular_CPD(bnet1, 4);

% 得到 inference engine

engine = smoother_engine(jtree_2TBN_inf_engine(bnet1));

ss = 3;             % slice size(ss)

ncases = 50;    %number of examples

T = 10;

%% ----------------------- 模型预测 ------------------------------ %%

max_iter = 10;    %iterations for EM

cases = cell(1, ncases);

for i=1:ncases

  ev = sample_dbn(bnet, T);

  cases{i} = cell(ss,T);

  cases{i}(onodes,:) = ev(onodes, :);

end

[bnet2, LLtrace] = learn_params_dbn_em(engine, cases, 'max_iter', max_iter);

% show estimated parameters

prior0

s=struct(bnet2.CPD{1});

s.CPT

clc

transmat0

s=struct(bnet2.CPD{4});

s.CPT

clc

mixmat0

s=struct(bnet2.CPD{2});

s.CPT

clc

s=struct(bnet2.CPD{3});

mu0

s.mean

Sigma0

s.cov

%% ------------------- State 预测（Viterbi） -------------------------%%

% Create a sequence

test_case = cell(ss,T);

test_case = sample_dbn(bnet, T);

evidence = cell(ss,T);

evidence(onodes,:) = test_case(onodes,:);

data = cell2num(test_case(onodes,:));

% Viterbi算法得到的路径

B = mixgauss_prob(data,mu0,Sigma0,mixmat0);

path = viterbi_path(prior0, transmat0, B);  

% Most Probable Explanation

engine = {};

engine{end+1} = smoother_engine(jtree_2TBN_inf_engine(bnet));

mpe = find_mpe(engine{1}, evidence);

assert(isequal(cell2num(mpe(1,:)), path))     % extract values of hidden nodes