Daniel@0: % pigs model from Lauritzen and Nilsson, 2001
Daniel@0: 
Daniel@0: seed = 0;
Daniel@0: rand('state', seed);
Daniel@0: randn('state', seed);
Daniel@0: 
Daniel@0: % we number nodes down and to the right
Daniel@0: h = [1 5 9 13];
Daniel@0: t = [2 6 10];
Daniel@0: d = [3 7 11];
Daniel@0: u = [4 8 12 14];
Daniel@0: 
Daniel@0: N = 14;
Daniel@0: dag = zeros(N);
Daniel@0: 
Daniel@0: % causal arcs
Daniel@0: for i=1:3
Daniel@0:   dag(h(i), [t(i) h(i+1)]) = 1;
Daniel@0:   dag(d(i), [u(i) h(i+1)]) = 1;
Daniel@0: end
Daniel@0: dag(h(4), u(4)) = 1;
Daniel@0: 
Daniel@0: % information arcs
Daniel@0: fig = 2;
Daniel@0: switch fig
Daniel@0:  case 0,
Daniel@0:   % no info arcs
Daniel@0:  case 1,
Daniel@0:    % no-forgetting policy (figure 1)
Daniel@0:    for i=1:3
Daniel@0:      dag(t(i), d(i:3)) = 1;
Daniel@0:    end
Daniel@0:  case 2,
Daniel@0:   % reactive policy (figure 2)
Daniel@0:   for i=1:3
Daniel@0:     dag(t(i), d(i)) = 1;
Daniel@0:   end
Daniel@0:  case 7,
Daniel@0:   % omniscient policy (figure 7: di has access to hidden state h(i-1))
Daniel@0:   dag(t(1), d(1)) = 1;
Daniel@0:   for i=2:3
Daniel@0:     %dag([h(i-1) t(i-1) d(i-1)], d(i)) = 1;
Daniel@0:     dag([h(i-1) d(i-1)], d(i)) = 1; % t(i-1) is redundant given h(i-1)
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: ns = 2*ones(1,N);
Daniel@0: ns(u) = 1;
Daniel@0: 
Daniel@0: % parameter tying
Daniel@0: params = ones(1,N);
Daniel@0: uparam = 1;
Daniel@0: final_uparam = 2;
Daniel@0: tparam = 3;
Daniel@0: h1_param = 4;
Daniel@0: hparam = 5;
Daniel@0: dparams = 6:8;
Daniel@0: 
Daniel@0: params(u(1:3)) = uparam;
Daniel@0: params(u(4)) = final_uparam;
Daniel@0: params(t) = tparam;
Daniel@0: params(h(1)) = h1_param;
Daniel@0: params(h(2:end)) = hparam;
Daniel@0: params(d) = dparams;
Daniel@0: 
Daniel@0: limid = mk_limid(dag, ns, 'chance', [h t], 'decision', d, 'utility', u, 'equiv_class', params);
Daniel@0: 
Daniel@0: % h = 1 means healthy, h = 2 means diseased
Daniel@0: % d = 1 means don't treat, d = 2 means treat
Daniel@0: % t = 1 means test shows healthy, t = 2 means test shows diseased
Daniel@0: 
Daniel@0: if 0
Daniel@0:   % use random params
Daniel@0:   limid.CPD{final_uparam} = tabular_utility_node(limid, u(4));
Daniel@0:   limid.CPD{uparam} = tabular_utility_node(limid, u(1));
Daniel@0:   limid.CPD{tparam} = tabular_CPD(limid, t(1));
Daniel@0:   limid.CPD{h1_param} = tabular_CPD(limid, h(1));
Daniel@0:   limid.CPD{hparam} = tabular_CPD(limid, h(2));
Daniel@0: else
Daniel@0:   limid.CPD{final_uparam} = tabular_utility_node(limid, u(4), [1000 300]);
Daniel@0:   limid.CPD{uparam} = tabular_utility_node(limid, u(1), [0 -100]); % costs have negative utility!
Daniel@0:   
Daniel@0:   % h  P(t=1) P(t=2)
Daniel@0:   % 1  0.9   0.1
Daniel@0:   % 2  0.2   0.8
Daniel@0:   limid.CPD{tparam} = tabular_CPD(limid, t(1), [0.9 0.2 0.1 0.8]);
Daniel@0:   
Daniel@0:   % P(h1)
Daniel@0:   limid.CPD{h1_param} = tabular_CPD(limid, h(1), [0.9 0.1]);
Daniel@0:   
Daniel@0:   % hi di P(hj=1) P(hj=2),  j = i+1, i=1:3
Daniel@0:   % 1  1  0.8     0.2
Daniel@0:   % 2  1  0.1     0.9
Daniel@0:   % 1  2  0.9     0.1
Daniel@0:   % 2  2  0.5     0.5
Daniel@0:   limid.CPD{hparam} = tabular_CPD(limid, h(2), [0.8 0.1 0.9 0.5 0.2 0.9 0.1 0.5]);
Daniel@0: end
Daniel@0: 
Daniel@0: % Decision nodes get assigned uniform policies by default
Daniel@0: for i=1:3
Daniel@0:   limid.CPD{dparams(i)} = tabular_decision_node(limid, d(i));
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: fname = '/home/cs/murphyk/matlab/Misc/loopybel.txt';
Daniel@0: 
Daniel@0: engines = {};
Daniel@0: engines{end+1} = global_joint_inf_engine(limid);
Daniel@0: engines{end+1} = jtree_limid_inf_engine(limid);
Daniel@0: %engines{end+1} = belprop_inf_engine(limid, 'max_iter', 1*N, 'filename', fname, 'tol', 1e-3);
Daniel@0: 
Daniel@0: exact = [1 2];
Daniel@0: %approx = 3;
Daniel@0: approx = [];
Daniel@0: 
Daniel@0: max_iter = 1;
Daniel@0: order = d(end:-1:1);
Daniel@0: %order = d(1:end);
Daniel@0: 
Daniel@0: NE = length(engines);
Daniel@0: MEU = zeros(1, NE);
Daniel@0: niter = zeros(1, NE);
Daniel@0: strategy = cell(1, NE);
Daniel@0: for e=1:NE
Daniel@0:   [strategy{e}, MEU(e), niter(e)] = solve_limid(engines{e}, 'max_iter', max_iter, 'order',  order);
Daniel@0: end
Daniel@0: MEU
Daniel@0: 
Daniel@0: % check results match those in the paper (p. 22)
Daniel@0: direct_policy = eye(2); % treat iff test is positive
Daniel@0: never_policy = [1 0; 1 0]; % never treat
Daniel@0: tol = 1e-0; % results in paper are reported to 0dp
Daniel@0: for e=exact(:)'
Daniel@0:   switch fig
Daniel@0:    case 2, % reactive policy
Daniel@0:     assert(approxeq(MEU(e), 727, tol));
Daniel@0:     assert(approxeq(strategy{e}{d(1)}(:), never_policy(:)))
Daniel@0:     assert(approxeq(strategy{e}{d(2)}(:), direct_policy(:)))
Daniel@0:     assert(approxeq(strategy{e}{d(3)}(:), direct_policy(:)))
Daniel@0:    case 1, assert(approxeq(MEU(e), 729, tol));
Daniel@0:    case 7, assert(approxeq(MEU(e), 732, tol));
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: for e=approx(:)'
Daniel@0:   for i=1:3
Daniel@0:     approxeq(strategy{exact(1)}{d(i)}, strategy{e}{d(i)})
Daniel@0:     dispcpt(strategy{e}{d(i)})
Daniel@0:   end
Daniel@0: end
Daniel@0: