wolffd@0: function bnet  = mk_incinerator_bnet(ns)
wolffd@0: % MK_INCINERATOR_BNET The waste incinerator emissions example from Cowell et al p145
wolffd@0: % function bnet  = mk_incinerator_bnet(ns)
wolffd@0: % 
wolffd@0: % If ns is omitted, we use the scalars and binary nodes and the original params.
wolffd@0: % Otherwise, we use random params of the desired size.
wolffd@0: %
wolffd@0: % Lauritzen, "Propogation of Probabilities, Means and Variances in Mixed Graphical Association Models", 
wolffd@0: % JASA 87(420): 1098--1108
wolffd@0: % This example is reprinted on p145 of "Probabilistic Networks and Expert Systems",
wolffd@0: % Cowell, Dawid, Lauritzen and Spiegelhalter, 1999, Springer. 
wolffd@0: % For a picture, see http://www.cs.berkeley.edu/~murphyk/Bayes/usage.html#cg_model
wolffd@0: 
wolffd@0: % node numbers
wolffd@0: F = 1; W = 2; E = 3; B = 4; C = 5; D = 6; Min = 7; Mout = 8; L = 9;
wolffd@0: names = {'F', 'W', 'E', 'B', 'C', 'D', 'Min', 'Mout', 'L'};
wolffd@0: n = 9;
wolffd@0: dnodes = [F W B];
wolffd@0: cnodes = mysetdiff(1:n, dnodes);
wolffd@0: 
wolffd@0: % node sizes - all cts nodes are scalar, all discrete nodes are binary
wolffd@0: if nargin < 1
wolffd@0:   ns = ones(1, n);
wolffd@0:   ns(dnodes) = 2;
wolffd@0:   rnd = 0;
wolffd@0: else
wolffd@0:   rnd = 1;
wolffd@0: end
wolffd@0:   
wolffd@0: % topology (p 1099, fig 1)
wolffd@0: dag = zeros(n);
wolffd@0: dag(F,E)=1;
wolffd@0: dag(W,[E Min D]) = 1;
wolffd@0: dag(E,D)=1;
wolffd@0: dag(B,[C D])=1;
wolffd@0: dag(D,[L Mout])=1;
wolffd@0: dag(Min,Mout)=1;
wolffd@0: 
wolffd@0: % params (p 1102)
wolffd@0: bnet = mk_bnet(dag, ns, 'discrete', dnodes, 'names', names);
wolffd@0: 
wolffd@0: if rnd
wolffd@0:   for i=dnodes(:)'
wolffd@0:     bnet.CPD{i} = tabular_CPD(bnet, i);
wolffd@0:   end
wolffd@0:   for i=cnodes(:)'
wolffd@0:     bnet.CPD{i} = gaussian_CPD(bnet, i);
wolffd@0:   end
wolffd@0: else
wolffd@0:   bnet.CPD{B} = tabular_CPD(bnet, B, 'CPT', [0.85 0.15]); % 1=stable, 2=unstable
wolffd@0:   bnet.CPD{F} = tabular_CPD(bnet, F, 'CPT', [0.95 0.05]); % 1=intact, 2=defect
wolffd@0:   bnet.CPD{W} = tabular_CPD(bnet, W, 'CPT', [2/7 5/7]); % 1=industrial, 2=household
wolffd@0:   bnet.CPD{E} = gaussian_CPD(bnet, E, 'mean', [-3.9 -0.4 -3.2 -0.5], ...
wolffd@0: 			     'cov', [0.00002 0.0001 0.00002 0.0001]);
wolffd@0:   bnet.CPD{D} = gaussian_CPD(bnet, D, 'mean', [6.5 6.0 7.5 7.0], ...
wolffd@0: 			     'cov', [0.03 0.04 0.1 0.1], 'weights', [1 1 1 1]);
wolffd@0:   bnet.CPD{C} = gaussian_CPD(bnet, C, 'mean', [-2 -1], 'cov', [0.1 0.3]);
wolffd@0:   bnet.CPD{L} = gaussian_CPD(bnet, L, 'mean', 3, 'cov', 0.25, 'weights', -0.5);
wolffd@0:   bnet.CPD{Min} = gaussian_CPD(bnet, Min, 'mean', [0.5 -0.5], 'cov', [0.01 0.005]);
wolffd@0:   bnet.CPD{Mout} = gaussian_CPD(bnet, Mout, 'mean', 0, 'cov', 0.002, 'weights', [1 1]);
wolffd@0: end