Daniel@0: 
Daniel@0: clear all
Daniel@0: B0 = 1; Rtriple = 2; Damnio = 3;
Daniel@0: B1 = 4; Ramnio = 5; Dabort = 6; 
Daniel@0: B2 = 7; U = 8;
Daniel@0: 
Daniel@0: N = 8;
Daniel@0: dag = zeros(N,N);
Daniel@0: dag(B0, [Rtriple B1 Ramnio]) = 1;
Daniel@0: dag(Rtriple, [Damnio Dabort]) = 1;
Daniel@0: dag(Damnio, [B1 Ramnio]) = 1;
Daniel@0: dag(B1, B2) = 1;
Daniel@0: dag(Ramnio, [Dabort U]) = 1;
Daniel@0: dag(Dabort, B2) = 1;
Daniel@0: dag(B2, U) = 1;
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: ns = zeros(1,N);
Daniel@0: ns(B0) = 2;
Daniel@0: ns(B1) = 3;
Daniel@0: ns(B2) = 4;
Daniel@0: ns(Rtriple) = 2;
Daniel@0: ns(Ramnio) = 3;
Daniel@0: ns(Damnio) = 2;
Daniel@0: ns(Dabort) = 2;
Daniel@0: ns(U) = 1;
Daniel@0: 
Daniel@0: limid = mk_limid(dag, ns, 'chance', [B0 B1 B2], ...
Daniel@0: 		 'decision', [Damnio Dabort], 'utility', [U]);
Daniel@0: 
Daniel@0: % states of nature
Daniel@0: healthy = 1; downs = 2; miscarry = 3; aborted = 4;
Daniel@0: % test results
Daniel@0: pos = 1; neg = 2; unk = 3;
Daniel@0: % actions
Daniel@0: yes = 1; no = 2;
Daniel@0: 
Daniel@0: % Prior probability baby has downs syndrome
Daniel@0: tbl = zeros(2,1);
Daniel@0: p = 1/1000; % from www.downs-syndrome.org.uk figure
Daniel@0: p = 24/10000; % www-personal.umich.edu/~bobwolfe/560/review/Downs.pdf (for women agen 35-40)
Daniel@0: tbl(healthy) = 1-p;
Daniel@0: tbl(downs) = p;
Daniel@0: limid.CPD{B0} = tabular_CPD(limid, B0, tbl);
Daniel@0: 
Daniel@0: % Reliability of triple screen test
Daniel@0: % Unreliable sensor
Daniel@0: % B0 -> Rtriple
Daniel@0: tbl = zeros(2,2); % Rtriple = pos, neg
Daniel@0: p = 0.5;  % high false positive rate (guess)
Daniel@0: tbl(healthy, :) = [p 1-p];
Daniel@0: p = 0.6; % low detection rate (march of dimes figure)
Daniel@0: tbl(downs, :) = [p 1-p]; 
Daniel@0: limid.CPD{Rtriple} = tabular_CPD(limid, Rtriple, tbl);
Daniel@0: 
Daniel@0: limid.CPD{Damnio} = tabular_decision_node(limid, Damnio);
Daniel@0: 
Daniel@0: % Effect of amnio on baby  B0,Damnio -> B1
Daniel@0:  % 1/200 risk of miscarry 
Daniel@0: p = 1/200; % (march of dimes figure)
Daniel@0: tbl = zeros(2, 2, 3); % B1 = healthy, downs, miscarry
Daniel@0: tbl(healthy, no, :) =  [1     0     0];
Daniel@0: tbl(downs, no, :) =    [0     1     0];
Daniel@0: tbl(healthy, yes, :) = [1-p     0   p];
Daniel@0: tbl(downs, yes, :) =   [0     1-p   p];
Daniel@0: limid.CPD{B1} = tabular_CPD(limid, B1, tbl);
Daniel@0: 
Daniel@0: % Reliability of amnio  B0, Damnio -> Ramnio
Daniel@0: % Perfect sensor
Daniel@0: tbl = zeros(2,2,3); % Ramnio = pos, neg, unk
Daniel@0: tbl(:, no, :) =        repmat([0 0 1], 2 ,1);
Daniel@0: tbl(healthy, yes, :) = [0 1 0]; 
Daniel@0: tbl(downs, yes, :) =   [1 0 0]; 
Daniel@0: limid.CPD{Ramnio} = tabular_CPD(limid, Ramnio, tbl);
Daniel@0: 
Daniel@0: limid.CPD{Dabort} = tabular_decision_node(limid, Dabort);
Daniel@0: 
Daniel@0: % Effect of abortion on baby  B1, Dabort -> B2
Daniel@0: tbl = zeros(3, 2, 4); % B2 = healthy, downs, miscarry, aborted
Daniel@0: tbl(:, yes, :) =       repmat([0 0 0 1], 3, 1);
Daniel@0: tbl(healthy, no, :) =  [1 0 0 0];
Daniel@0: tbl(downs, no, :) =    [0 1 0 0];
Daniel@0: tbl(miscarry, no, :) = [0 0 1 0];
Daniel@0: limid.CPD{B2} = tabular_CPD(limid, B2, tbl);
Daniel@0: 
Daniel@0: % Utility U(Ramnio, B2)
Daniel@0: tbl = zeros(3, 4);
Daniel@0: tbl(:, healthy) = 5000;
Daniel@0: tbl(:, downs) = -50000;
Daniel@0: tbl(:, miscarry) = -1000;
Daniel@0: tbl(:, aborted) = -1000;
Daniel@0: 
Daniel@0: if 0
Daniel@0: %tbl(unk, miscarry) = 0; % this case is impossible
Daniel@0: tbl(pos, miscarry) = -1;
Daniel@0: tbl(neg, miscarry) = -1000;
Daniel@0: if 1
Daniel@0:   tbl(unk, aborted) = -100;
Daniel@0:   tbl(pos, aborted) = -1;
Daniel@0:   tbl(neg, aborted) = -500;
Daniel@0: else % pro-life utility fn
Daniel@0:   tbl(unk, aborted) = -500000;
Daniel@0:   tbl(pos, aborted) = -500000;
Daniel@0:   tbl(neg, aborted) = -500000;
Daniel@0: end 
Daniel@0: end
Daniel@0: 
Daniel@0: limid.CPD{U} = tabular_utility_node(limid, U,  tbl);
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: engine = jtree_limid_inf_engine(limid);
Daniel@0: [strategy, MEU] = solve_limid(engine);
Daniel@0: 
Daniel@0: % Rtriple U(Damnio=1=yes)  U(Damnio=2=no)
Daniel@0: % 1=pos    0               1
Daniel@0: % 2=neg    0               1
Daniel@0: dispcpt(strategy{Damnio})
Daniel@0: if isequal(strategy{Damnio}(1,:), strategy{Damnio}(2,:))
Daniel@0:   % Rtriple result irrelevant
Daniel@0:   doAmnio = argmax(strategy{Damnio}(1,:))
Daniel@0: else
Daniel@0:   doAmnio = 1;
Daniel@0: end
Daniel@0: 
Daniel@0: % Rtriple Ramnio U(Dabort=yes=1) U(Dabort=no=2)
Daniel@0: % 1=pos   1=pos  1               0
Daniel@0: % 2=neg   1=pos  1               0
Daniel@0: % 1=pos   2=neg  0               1
Daniel@0: % 2=neg   2=neg  0               1
Daniel@0: % 1=pos   3=unk  0               1
Daniel@0: % 2=neg   3=unk  0               1
Daniel@0: dispcpt(strategy{Dabort})
Daniel@0: