cg1.m

Bayesian网络工具箱.

function cg1

% The waste incinerator emissions example from Lauritzen (1992),
% "Propogation of Probabilities, Means and Variances in Mixed Graphical Association Models", 
% JASA 87(420): 1098--1108
%
% This example is reprinted in chap 7 of "Probabilistic Networks and Expert Systems",
% Cowell, Dawid, Lauritzen and Spiegelhalter, 1999, Springer.


% node numbers
F = 1; W = 2; E = 3; B = 4; C = 5; D = 6; Min = 7; Mout = 8; L = 9;
n = 9;

% node sizes - all cts nodes are scalar, all discrete nodes are binary
ns = ones(1, n);
dnodes = [F W B];
cnodes = mysetdiff(1:n, dnodes);
ns(dnodes) = 2;

% topology (p 1099, fig 1)
dag = zeros(n);
dag(F,E)=1;
dag(W,[E Min D]) = 1;
dag(E,D)=1;
dag(B,[C D])=1;
dag(D,[L Mout])=1;
dag(Min,Mout)=1;

% params (p 1102)
bnet = mk_bnet(dag, ns, 'discrete', dnodes);

bnet.CPD{B} = tabular_CPD(bnet, B, 'CPT', [0.85 0.15]); % 1=stable, 2=unstable
bnet.CPD{F} = tabular_CPD(bnet, F, 'CPT', [0.95 0.05]); % 1=intact, 2=defect
bnet.CPD{W} = tabular_CPD(bnet, W, 'CPT', [2/7 5/7]); % 1=industrial, 2=household
bnet.CPD{E} = gaussian_CPD(bnet, E, 'mean', [-3.9 -0.4 -3.2 -0.5], ...
			   'cov', [0.00002 0.0001 0.00002 0.0001]);
bnet.CPD{D} = gaussian_CPD(bnet, D, 'mean', [6.5 6.0 7.5 7.0], ...
			   'cov', [0.03 0.04 0.1 0.1], 'weights', [1 1 1 1]);
bnet.CPD{C} = gaussian_CPD(bnet, C, 'mean', [-2 -1], 'cov', [0.1 0.3]);
bnet.CPD{L} = gaussian_CPD(bnet, L, 'mean', 3, 'cov', 0.25, 'weights', -0.5);
bnet.CPD{Min} = gaussian_CPD(bnet, Min, 'mean', [0.5 -0.5], 'cov', [0.01 0.005]);
bnet.CPD{Mout} = gaussian_CPD(bnet, Mout, 'mean', 0, 'cov', 0.002, 'weights', [1 1]);


engine1 = jtree_inf_engine(bnet);
engine2 = cond_gauss_inf_engine(bnet);

% Initial values
evidence = cell(1,n);
[engine1, ll1] = enter_evidence(engine1, evidence);
[engine2, ll2] = enter_evidence(engine2, evidence);

% Compare to the results in table on p1107.
% These results are printed to 3dp in Cowell p150

global USEC
if USEC
  tol = 1; % hack because strong_elim_order.c is broken!
else
  tol = 1e-2; % match to 2dp
end

[dprob, mu, sigma] = extract_posteriors(engine1, evidence);
assert(approxeq(mu([E D C L Min Mout]), [-3.25 3.04 -1.85 1.48 -0.214 2.83], tol))
assert(approxeq(sigma([E D C L Min Mout]), [0.709 0.770 0.507 0.631 0.459 0.860], tol))
assert(approxeq(dprob([B F W]), [0.85 0.95 0.29], tol))

[dprob, mu, sigma] = extract_posteriors(engine2, evidence);
assert(approxeq(mu([E D C L Min Mout]), [-3.25 3.04 -1.85 1.48 -0.214 2.83], tol))
assert(approxeq(sigma([E D C L Min Mout]), [0.709 0.770 0.507 0.631 0.459 0.860], tol))
assert(approxeq(dprob([B F W]), [0.85 0.95 0.29], tol))



% Add evidence (p 1105, top right)
evidence = cell(1,n);
evidence{W} = 1; % industrial
evidence{L} = 1.1;
evidence{C} = -0.9;

[engine1, ll1] = enter_evidence(engine1, evidence);
[engine2, ll2] = enter_evidence(engine2, evidence);

[dprob, mu, sigma] = extract_posteriors(engine1, evidence);
assert(approxeq(mu([E D C L Min Mout]), [-3.90 3.61 -0.9 1.1 0.5 4.11], tol))
assert(approxeq(sigma([E D C L Min Mout]), [0.076 0.326 0 0 0.1 0.344], tol))
assert(approxeq(dprob([B F W]), [0.0122 0.9995 1], tol))

[dprob, mu, sigma] = extract_posteriors(engine2, evidence);
assert(approxeq(mu([E D C L Min Mout]), [-3.90 3.61 -0.9 1.1 0.5 4.11], tol))
assert(approxeq(sigma([E D C L Min Mout]), [0.076 0.326 0 0 0.1 0.344], tol))
assert(approxeq(dprob([B F W]), [0.0122 0.9995 1], tol))

assert(approxeq(ll1, ll2, tol))

marg = marginal_nodes(engine1, [D Mout]);

%%%%%%%%%%%%%

function [T, mu, sigma] = extract_posteriors(engine, evidence)

F = 1; W = 2; E = 3; B = 4; C = 5; D = 6; Min = 7; Mout = 8; L = 9;
n = 9;
dnodes = [B F W];
cnodes = mysetdiff(1:n, dnodes);

observed = ~isemptycell(evidence);
mu = zeros(1,n);
sigma = zeros(1,n);
for i=cnodes(:)'
  m = marginal_nodes(engine, i);
  if observed(i)
    mu(i) = evidence{i};
    sigma(i) = 0;
  else
    mu(i) = m.mu;
    sigma(i) = sqrt(m.Sigma);
  end
end

T = zeros(1,n);
for i=dnodes(:)'
  m = marginal_nodes(engine, i);
  if observed(i)
    if evidence{i}==1
      T(i) = 1;
    else
      T(i) = 0;
    end
  else
    T(i) = m.T(1);
  end
end