Daniel@0: % decision theoretic version of asia network
Daniel@0: % Cowell et al, p177
Daniel@0: % We explicitely add the no-forgetting arcs.
Daniel@0: 
Daniel@0: Smoking = 1;
Daniel@0: VisitToAsia = 2;
Daniel@0: Bronchitis = 3;
Daniel@0: LungCancer = 4;
Daniel@0: TB = 5;
Daniel@0: Do_Xray = 6;
Daniel@0: TBorCancer = 7;
Daniel@0: Util_Xray = 8;
Daniel@0: Dys = 9;
Daniel@0: posXray = 10;
Daniel@0: Do_Hosp = 11;
Daniel@0: Util_Hosp = 12;
Daniel@0: 
Daniel@0: n = 12;
Daniel@0: dag = zeros(n);
Daniel@0: dag(Smoking, [Bronchitis LungCancer]) = 1;
Daniel@0: dag(VisitToAsia, [TB Do_Xray Do_Hosp]) = 1;
Daniel@0: dag(Bronchitis, Dys) = 1;
Daniel@0: dag(LungCancer, [Util_Hosp TBorCancer]) = 1;
Daniel@0: dag(TB, [Util_Hosp TBorCancer Util_Xray]) = 1;
Daniel@0: dag(Do_Xray, [posXray Util_Xray Do_Hosp]) = 1;
Daniel@0: dag(TBorCancer, [Dys posXray]) = 1;
Daniel@0: dag(Dys, Do_Hosp) = 1;
Daniel@0: dag(posXray, Do_Hosp) = 1;
Daniel@0: dag(Do_Hosp, Util_Hosp) = 1;
Daniel@0: 
Daniel@0: dnodes = [Do_Xray Do_Hosp];
Daniel@0: unodes = [Util_Xray Util_Hosp];
Daniel@0: cnodes = mysetdiff(1:n, [dnodes unodes]); % chance nodes
Daniel@0: ns = 2*ones(1,n);
Daniel@0: ns(unodes) = 1;
Daniel@0: limid = mk_limid(dag, ns, 'chance', cnodes, 'decision', dnodes, 'utility', unodes);
Daniel@0: 
Daniel@0: % 1 = yes, 2 = no
Daniel@0: limid.CPD{VisitToAsia} = tabular_CPD(limid, VisitToAsia, [0.01 0.99]);
Daniel@0: limid.CPD{Bronchitis} = tabular_CPD(limid, Bronchitis, [0.6 0.3  0.4 0.7]);
Daniel@0: limid.CPD{Dys} = tabular_CPD(limid, Dys, [0.9 0.7 0.8 0.1  0.1 0.3 0.2 0.9]);
Daniel@0: limid.CPD{TBorCancer} = tabular_CPD(limid, TBorCancer, [1 1 1 0  0 0 0 1]);
Daniel@0: 
Daniel@0: limid.CPD{LungCancer} = tabular_CPD(limid, LungCancer, [0.1 0.01  0.9 0.99]);
Daniel@0: limid.CPD{Smoking} = tabular_CPD(limid, Smoking, [0.5 0.5]);
Daniel@0: limid.CPD{TB} = tabular_CPD(limid, TB, [0.05 0.01  0.95 0.99]);
Daniel@0: limid.CPD{posXray} = tabular_CPD(limid, posXray, [0.98 0.5 0.05 0.5  0.02 0.5 0.95 0.5]);
Daniel@0: 
Daniel@0: limid.CPD{Util_Hosp} = tabular_utility_node(limid, Util_Hosp, [180 120 160 15  2 4 0 40]);
Daniel@0: limid.CPD{Util_Xray} = tabular_utility_node(limid, Util_Xray, [0 1 10 10]);
Daniel@0: 
Daniel@0: for i=dnodes(:)'
Daniel@0:   limid.CPD{i} = tabular_decision_node(limid, i);
Daniel@0: end
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);
Daniel@0: 
Daniel@0: exact = [1 2];
Daniel@0: %approx = 3;
Daniel@0: approx = [];
Daniel@0: 
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: 
Daniel@0: tol = 1e-2;
Daniel@0: for e=1:length(engines)
Daniel@0:   [strategy{e}, MEU(e), niter(e)] = solve_limid(engines{e});
Daniel@0: end
Daniel@0: 
Daniel@0: for e=exact(:)'
Daniel@0:   assert(approxeq(MEU(e), 47.49, tol))
Daniel@0:   assert(isequal(strategy{e}{Do_Xray}(:)', [1 0 0 1]))
Daniel@0:   
Daniel@0:   % Check the hosptialize strategy is correct (p180)
Daniel@0:   % We assume the patient has not been to Asia and therefore did not have an Xray.
Daniel@0:   % In this case it is optimal not to hospitalize regardless of whether the patient has
Daniel@0:   % dyspnoea or not (and of course regardless of the value of pos_xray).
Daniel@0:   asia = 2;
Daniel@0:   do_xray = 2;
Daniel@0:   for dys = 1:2
Daniel@0:     for pos_xray = 1:2
Daniel@0:       assert(argmax(squeeze(strategy{e}{Do_Hosp}(asia, do_xray, dys, pos_xray, :))) == 2)
Daniel@0:     end
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: for e=approx(:)'
Daniel@0:   approxeq(strategy{exact(1)}{Do_Xray}, strategy{e}{Do_Xray})
Daniel@0:   approxeq(strategy{exact(1)}{Do_Hosp}, strategy{e}{Do_Hosp})
Daniel@0: end
Daniel@0: 
Daniel@0:  
Daniel@0: