Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/examples/limids/pigs1.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 % pigs model from Lauritzen and Nilsson, 2001
wolffd@0 2
wolffd@0 3 seed = 0;
wolffd@0 4 rand('state', seed);
wolffd@0 5 randn('state', seed);
wolffd@0 6
wolffd@0 7 % we number nodes down and to the right
wolffd@0 8 h = [1 5 9 13];
wolffd@0 9 t = [2 6 10];
wolffd@0 10 d = [3 7 11];
wolffd@0 11 u = [4 8 12 14];
wolffd@0 12
wolffd@0 13 N = 14;
wolffd@0 14 dag = zeros(N);
wolffd@0 15
wolffd@0 16 % causal arcs
wolffd@0 17 for i=1:3
wolffd@0 18 dag(h(i), [t(i) h(i+1)]) = 1;
wolffd@0 19 dag(d(i), [u(i) h(i+1)]) = 1;
wolffd@0 20 end
wolffd@0 21 dag(h(4), u(4)) = 1;
wolffd@0 22
wolffd@0 23 % information arcs
wolffd@0 24 fig = 2;
wolffd@0 25 switch fig
wolffd@0 26 case 0,
wolffd@0 27 % no info arcs
wolffd@0 28 case 1,
wolffd@0 29 % no-forgetting policy (figure 1)
wolffd@0 30 for i=1:3
wolffd@0 31 dag(t(i), d(i:3)) = 1;
wolffd@0 32 end
wolffd@0 33 case 2,
wolffd@0 34 % reactive policy (figure 2)
wolffd@0 35 for i=1:3
wolffd@0 36 dag(t(i), d(i)) = 1;
wolffd@0 37 end
wolffd@0 38 case 7,
wolffd@0 39 % omniscient policy (figure 7: di has access to hidden state h(i-1))
wolffd@0 40 dag(t(1), d(1)) = 1;
wolffd@0 41 for i=2:3
wolffd@0 42 %dag([h(i-1) t(i-1) d(i-1)], d(i)) = 1;
wolffd@0 43 dag([h(i-1) d(i-1)], d(i)) = 1; % t(i-1) is redundant given h(i-1)
wolffd@0 44 end
wolffd@0 45 end
wolffd@0 46
wolffd@0 47
wolffd@0 48 ns = 2*ones(1,N);
wolffd@0 49 ns(u) = 1;
wolffd@0 50
wolffd@0 51 % parameter tying
wolffd@0 52 params = ones(1,N);
wolffd@0 53 uparam = 1;
wolffd@0 54 final_uparam = 2;
wolffd@0 55 tparam = 3;
wolffd@0 56 h1_param = 4;
wolffd@0 57 hparam = 5;
wolffd@0 58 dparams = 6:8;
wolffd@0 59
wolffd@0 60 params(u(1:3)) = uparam;
wolffd@0 61 params(u(4)) = final_uparam;
wolffd@0 62 params(t) = tparam;
wolffd@0 63 params(h(1)) = h1_param;
wolffd@0 64 params(h(2:end)) = hparam;
wolffd@0 65 params(d) = dparams;
wolffd@0 66
wolffd@0 67 limid = mk_limid(dag, ns, 'chance', [h t], 'decision', d, 'utility', u, 'equiv_class', params);
wolffd@0 68
wolffd@0 69 % h = 1 means healthy, h = 2 means diseased
wolffd@0 70 % d = 1 means don't treat, d = 2 means treat
wolffd@0 71 % t = 1 means test shows healthy, t = 2 means test shows diseased
wolffd@0 72
wolffd@0 73 if 0
wolffd@0 74 % use random params
wolffd@0 75 limid.CPD{final_uparam} = tabular_utility_node(limid, u(4));
wolffd@0 76 limid.CPD{uparam} = tabular_utility_node(limid, u(1));
wolffd@0 77 limid.CPD{tparam} = tabular_CPD(limid, t(1));
wolffd@0 78 limid.CPD{h1_param} = tabular_CPD(limid, h(1));
wolffd@0 79 limid.CPD{hparam} = tabular_CPD(limid, h(2));
wolffd@0 80 else
wolffd@0 81 limid.CPD{final_uparam} = tabular_utility_node(limid, u(4), [1000 300]);
wolffd@0 82 limid.CPD{uparam} = tabular_utility_node(limid, u(1), [0 -100]); % costs have negative utility!
wolffd@0 83
wolffd@0 84 % h P(t=1) P(t=2)
wolffd@0 85 % 1 0.9 0.1
wolffd@0 86 % 2 0.2 0.8
wolffd@0 87 limid.CPD{tparam} = tabular_CPD(limid, t(1), [0.9 0.2 0.1 0.8]);
wolffd@0 88
wolffd@0 89 % P(h1)
wolffd@0 90 limid.CPD{h1_param} = tabular_CPD(limid, h(1), [0.9 0.1]);
wolffd@0 91
wolffd@0 92 % hi di P(hj=1) P(hj=2), j = i+1, i=1:3
wolffd@0 93 % 1 1 0.8 0.2
wolffd@0 94 % 2 1 0.1 0.9
wolffd@0 95 % 1 2 0.9 0.1
wolffd@0 96 % 2 2 0.5 0.5
wolffd@0 97 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 98 end
wolffd@0 99
wolffd@0 100 % Decision nodes get assigned uniform policies by default
wolffd@0 101 for i=1:3
wolffd@0 102 limid.CPD{dparams(i)} = tabular_decision_node(limid, d(i));
wolffd@0 103 end
wolffd@0 104
wolffd@0 105
wolffd@0 106 fname = '/home/cs/murphyk/matlab/Misc/loopybel.txt';
wolffd@0 107
wolffd@0 108 engines = {};
wolffd@0 109 engines{end+1} = global_joint_inf_engine(limid);
wolffd@0 110 engines{end+1} = jtree_limid_inf_engine(limid);
wolffd@0 111 %engines{end+1} = belprop_inf_engine(limid, 'max_iter', 1*N, 'filename', fname, 'tol', 1e-3);
wolffd@0 112
wolffd@0 113 exact = [1 2];
wolffd@0 114 %approx = 3;
wolffd@0 115 approx = [];
wolffd@0 116
wolffd@0 117 max_iter = 1;
wolffd@0 118 order = d(end:-1:1);
wolffd@0 119 %order = d(1:end);
wolffd@0 120
wolffd@0 121 NE = length(engines);
wolffd@0 122 MEU = zeros(1, NE);
wolffd@0 123 niter = zeros(1, NE);
wolffd@0 124 strategy = cell(1, NE);
wolffd@0 125 for e=1:NE
wolffd@0 126 [strategy{e}, MEU(e), niter(e)] = solve_limid(engines{e}, 'max_iter', max_iter, 'order', order);
wolffd@0 127 end
wolffd@0 128 MEU
wolffd@0 129
wolffd@0 130 % check results match those in the paper (p. 22)
wolffd@0 131 direct_policy = eye(2); % treat iff test is positive
wolffd@0 132 never_policy = [1 0; 1 0]; % never treat
wolffd@0 133 tol = 1e-0; % results in paper are reported to 0dp
wolffd@0 134 for e=exact(:)'
wolffd@0 135 switch fig
wolffd@0 136 case 2, % reactive policy
wolffd@0 137 assert(approxeq(MEU(e), 727, tol));
wolffd@0 138 assert(approxeq(strategy{e}{d(1)}(:), never_policy(:)))
wolffd@0 139 assert(approxeq(strategy{e}{d(2)}(:), direct_policy(:)))
wolffd@0 140 assert(approxeq(strategy{e}{d(3)}(:), direct_policy(:)))
wolffd@0 141 case 1, assert(approxeq(MEU(e), 729, tol));
wolffd@0 142 case 7, assert(approxeq(MEU(e), 732, tol));
wolffd@0 143 end
wolffd@0 144 end
wolffd@0 145
wolffd@0 146
wolffd@0 147 for e=approx(:)'
wolffd@0 148 for i=1:3
wolffd@0 149 approxeq(strategy{exact(1)}{d(i)}, strategy{e}{d(i)})
wolffd@0 150 dispcpt(strategy{e}{d(i)})
wolffd@0 151 end
wolffd@0 152 end
wolffd@0 153