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