|
wolffd@0
|
1 % Conditional Gaussian network
|
|
wolffd@0
|
2 % The waste incinerator emissions example from Lauritzen (1992),
|
|
wolffd@0
|
3 % "Propogation of Probabilities, Means and Variances in Mixed Graphical Association Models",
|
|
wolffd@0
|
4 % JASA 87(420): 1098--1108
|
|
wolffd@0
|
5 %
|
|
wolffd@0
|
6 % This example is reprinted on p145 of "Probabilistic Networks and Expert Systems",
|
|
wolffd@0
|
7 % Cowell, Dawid, Lauritzen and Spiegelhalter, 1999, Springer.
|
|
wolffd@0
|
8 %
|
|
wolffd@0
|
9 % For a picture, see http://www.cs.berkeley.edu/~murphyk/Bayes/usage.html#cg_model
|
|
wolffd@0
|
10
|
|
wolffd@0
|
11 ns = 2*ones(1,9);
|
|
wolffd@0
|
12 %bnet = mk_incinerator_bnet(ns);
|
|
wolffd@0
|
13 bnet = mk_incinerator_bnet;
|
|
wolffd@0
|
14
|
|
wolffd@0
|
15 engines = {};
|
|
wolffd@0
|
16 %engines{end+1} = stab_cond_gauss_inf_engine(bnet);
|
|
wolffd@0
|
17 engines{end+1} = jtree_inf_engine(bnet);
|
|
wolffd@0
|
18 engines{end+1} = cond_gauss_inf_engine(bnet);
|
|
wolffd@0
|
19 nengines = length(engines);
|
|
wolffd@0
|
20
|
|
wolffd@0
|
21 F = 1; W = 2; E = 3; B = 4; C = 5; D = 6; Min = 7; Mout = 8; L = 9;
|
|
wolffd@0
|
22 n = 9;
|
|
wolffd@0
|
23 dnodes = [B F W];
|
|
wolffd@0
|
24 cnodes = mysetdiff(1:n, dnodes);
|
|
wolffd@0
|
25
|
|
wolffd@0
|
26 evidence = cell(1,n); % no evidence
|
|
wolffd@0
|
27 ll = zeros(1, nengines);
|
|
wolffd@0
|
28 for e=1:nengines
|
|
wolffd@0
|
29 [engines{e}, ll(e)] = enter_evidence(engines{e}, evidence);
|
|
wolffd@0
|
30 end
|
|
wolffd@0
|
31 %assert(approxeq(ll(1), ll)))
|
|
wolffd@0
|
32 ll
|
|
wolffd@0
|
33
|
|
wolffd@0
|
34 % Compare to the results in table on p1107.
|
|
wolffd@0
|
35 % These results are printed to 3dp in Cowell p150
|
|
wolffd@0
|
36
|
|
wolffd@0
|
37 mu = zeros(1,n);
|
|
wolffd@0
|
38 sigma = zeros(1,n);
|
|
wolffd@0
|
39 dprob = zeros(1,n);
|
|
wolffd@0
|
40 addev = 1;
|
|
wolffd@0
|
41 tol = 1e-2;
|
|
wolffd@0
|
42 for e=1:nengines
|
|
wolffd@0
|
43 for i=cnodes(:)'
|
|
wolffd@0
|
44 m = marginal_nodes(engines{e}, i, addev);
|
|
wolffd@0
|
45 mu(i) = m.mu;
|
|
wolffd@0
|
46 sigma(i) = sqrt(m.Sigma);
|
|
wolffd@0
|
47 end
|
|
wolffd@0
|
48 for i=dnodes(:)'
|
|
wolffd@0
|
49 m = marginal_nodes(engines{e}, i, addev);
|
|
wolffd@0
|
50 dprob(i) = m.T(1);
|
|
wolffd@0
|
51 end
|
|
wolffd@0
|
52 assert(approxeq(mu([E D C L Min Mout]), [-3.25 3.04 -1.85 1.48 -0.214 2.83], tol))
|
|
wolffd@0
|
53 assert(approxeq(sigma([E D C L Min Mout]), [0.709 0.770 0.507 0.631 0.459 0.860], tol))
|
|
wolffd@0
|
54 assert(approxeq(dprob([B F W]), [0.85 0.95 0.29], tol))
|
|
wolffd@0
|
55 %m = marginal_nodes(engines{e}, bnet.names('E'), addev);
|
|
wolffd@0
|
56 %assert(approxeq(m.mu, -3.25, tol))
|
|
wolffd@0
|
57 %assert(approxeq(sqrt(m.Sigma), 0.709, tol))
|
|
wolffd@0
|
58 end
|
|
wolffd@0
|
59
|
|
wolffd@0
|
60 % Add evidence (p 1105, top right)
|
|
wolffd@0
|
61 evidence = cell(1,n);
|
|
wolffd@0
|
62 evidence{W} = 1; % industrial
|
|
wolffd@0
|
63 evidence{L} = 1.1;
|
|
wolffd@0
|
64 evidence{C} = -0.9;
|
|
wolffd@0
|
65
|
|
wolffd@0
|
66 ll = zeros(1, nengines);
|
|
wolffd@0
|
67 for e=1:nengines
|
|
wolffd@0
|
68 [engines{e}, ll(e)] = enter_evidence(engines{e}, evidence);
|
|
wolffd@0
|
69 end
|
|
wolffd@0
|
70 assert(all(approxeq(ll(1), ll)))
|
|
wolffd@0
|
71
|
|
wolffd@0
|
72 for e=1:nengines
|
|
wolffd@0
|
73 for i=cnodes(:)'
|
|
wolffd@0
|
74 m = marginal_nodes(engines{e}, i, addev);
|
|
wolffd@0
|
75 mu(i) = m.mu;
|
|
wolffd@0
|
76 sigma(i) = sqrt(m.Sigma);
|
|
wolffd@0
|
77 end
|
|
wolffd@0
|
78 for i=dnodes(:)'
|
|
wolffd@0
|
79 m = marginal_nodes(engines{e}, i, addev);
|
|
wolffd@0
|
80 dprob(i) = m.T(1);
|
|
wolffd@0
|
81 end
|
|
wolffd@0
|
82 assert(approxeq(mu([E D C L Min Mout]), [-3.90 3.61 -0.9 1.1 0.5 4.11], tol))
|
|
wolffd@0
|
83 assert(approxeq(sigma([E D C L Min Mout]), [0.076 0.326 0 0 0.1 0.344], tol))
|
|
wolffd@0
|
84 assert(approxeq(dprob([B F W]), [0.0122 0.9995 1], tol))
|
|
wolffd@0
|
85 end
|
|
wolffd@0
|
86
|
