|
wolffd@0
|
1 % Test jtree_compiled on a toy QMR network.
|
|
wolffd@0
|
2
|
|
wolffd@0
|
3 rand('state', 0);
|
|
wolffd@0
|
4 randn('state', 0);
|
|
wolffd@0
|
5 pMax = 0.01;
|
|
wolffd@0
|
6 Nfindings = 10;
|
|
wolffd@0
|
7 Ndiseases = 5;
|
|
wolffd@0
|
8
|
|
wolffd@0
|
9 N=Nfindings+Ndiseases;
|
|
wolffd@0
|
10 findings = Ndiseases+1:N;
|
|
wolffd@0
|
11 diseases = 1:Ndiseases;
|
|
wolffd@0
|
12
|
|
wolffd@0
|
13 G = zeros(Ndiseases, Nfindings);
|
|
wolffd@0
|
14 for i=1:Nfindings
|
|
wolffd@0
|
15 v= rand(1,Ndiseases);
|
|
wolffd@0
|
16 rents = find(v<0.8);
|
|
wolffd@0
|
17 if (length(rents)==0)
|
|
wolffd@0
|
18 rents=ceil(rand(1)*Ndiseases);
|
|
wolffd@0
|
19 end
|
|
wolffd@0
|
20 G(rents,i)=1;
|
|
wolffd@0
|
21 end
|
|
wolffd@0
|
22
|
|
wolffd@0
|
23 prior = pMax*rand(1,Ndiseases);
|
|
wolffd@0
|
24 leak = 0.5*rand(1,Nfindings); % in real QMR, leak approx exp(-0.02) = 0.98
|
|
wolffd@0
|
25 %leak = ones(1,Nfindings); % turns off leaks, which makes inference much harder
|
|
wolffd@0
|
26 inhibit = rand(Ndiseases, Nfindings);
|
|
wolffd@0
|
27 inhibit(not(G)) = 1;
|
|
wolffd@0
|
28
|
|
wolffd@0
|
29 % first half of findings are +ve, second half -ve
|
|
wolffd@0
|
30 % The very first and last findings are hidden
|
|
wolffd@0
|
31 pos = 2:floor(Nfindings/2);
|
|
wolffd@0
|
32 neg = (pos(end)+1):(Nfindings-1);
|
|
wolffd@0
|
33
|
|
wolffd@0
|
34 big = 1;
|
|
wolffd@0
|
35
|
|
wolffd@0
|
36 if big
|
|
wolffd@0
|
37 % Make the bnet in the straightforward way
|
|
wolffd@0
|
38 tabular_leaves = 1;
|
|
wolffd@0
|
39 obs_nodes = myunion(pos, neg) + Ndiseases;
|
|
wolffd@0
|
40 bnet = mk_qmr_bnet(G, inhibit, leak, prior, tabular_leaves, obs_nodes);
|
|
wolffd@0
|
41 evidence = cell(1, N);
|
|
wolffd@0
|
42 evidence(findings(pos)) = num2cell(repmat(2, 1, length(pos)));
|
|
wolffd@0
|
43 evidence(findings(neg)) = num2cell(repmat(1, 1, length(neg)));
|
|
wolffd@0
|
44 else
|
|
wolffd@0
|
45 % Marginalize out hidden leaves apriori
|
|
wolffd@0
|
46 positive_leaves_only = 1;
|
|
wolffd@0
|
47 [bnet, vals] = mk_minimal_qmr_bnet(G, inhibit, leak, prior, pos, neg, positive_leaves_only);
|
|
wolffd@0
|
48 obs_nodes = bnet.observed;
|
|
wolffd@0
|
49 evidence = cell(1, Ndiseases + length(obs_nodes));
|
|
wolffd@0
|
50 evidence(obs_nodes) = num2cell(vals);
|
|
wolffd@0
|
51 end
|
|
wolffd@0
|
52
|
|
wolffd@0
|
53 engine = {};
|
|
wolffd@0
|
54 engine{end+1} = jtree_inf_engine(bnet);
|
|
wolffd@0
|
55
|
|
wolffd@0
|
56 E = length(engine);
|
|
wolffd@0
|
57 exact = 1:E;
|
|
wolffd@0
|
58 ll = zeros(1,E);
|
|
wolffd@0
|
59 for e=1:E
|
|
wolffd@0
|
60 tic; [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence); toc
|
|
wolffd@0
|
61 end
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 assert(all(approxeq(ll(exact), ll(exact(1)))))
|
|
wolffd@0
|
64
|
|
wolffd@0
|
65 post = zeros(E, Ndiseases);
|
|
wolffd@0
|
66 for e=1:E
|
|
wolffd@0
|
67 for i=diseases(:)'
|
|
wolffd@0
|
68 m = marginal_nodes(engine{e}, i);
|
|
wolffd@0
|
69 post(e, i) = m.T(2);
|
|
wolffd@0
|
70 end
|
|
wolffd@0
|
71 end
|
|
wolffd@0
|
72 for e=exact(:)'
|
|
wolffd@0
|
73 for i=diseases(:)'
|
|
wolffd@0
|
74 assert(approxeq(post(1, i), post(e, i)));
|
|
wolffd@0
|
75 end
|
|
wolffd@0
|
76 end
|
|
wolffd@0
|
77
|
