|
wolffd@0
|
1 % Make a QMR-like network
|
|
wolffd@0
|
2 % This is a bipartite graph, where the top layer contains hidden disease nodes,
|
|
wolffd@0
|
3 % and the bottom later contains observed finding nodes.
|
|
wolffd@0
|
4 % The diseases have Bernoulli CPDs, the findings noisy-or CPDs.
|
|
wolffd@0
|
5 % See quickscore_inf_engine for references.
|
|
wolffd@0
|
6
|
|
wolffd@0
|
7 pMax = 0.01;
|
|
wolffd@0
|
8 Nfindings = 10;
|
|
wolffd@0
|
9 Ndiseases = 5;
|
|
wolffd@0
|
10 %Nfindings = 20;
|
|
wolffd@0
|
11 %Ndiseases = 10;
|
|
wolffd@0
|
12
|
|
wolffd@0
|
13 N=Nfindings+Ndiseases;
|
|
wolffd@0
|
14 findings = Ndiseases+1:N;
|
|
wolffd@0
|
15 diseases = 1:Ndiseases;
|
|
wolffd@0
|
16
|
|
wolffd@0
|
17 G = zeros(Ndiseases, Nfindings);
|
|
wolffd@0
|
18 for i=1:Nfindings
|
|
wolffd@0
|
19 v= rand(1,Ndiseases);
|
|
wolffd@0
|
20 rents = find(v<0.8);
|
|
wolffd@0
|
21 if (length(rents)==0)
|
|
wolffd@0
|
22 rents=ceil(rand(1)*Ndiseases);
|
|
wolffd@0
|
23 end
|
|
wolffd@0
|
24 G(rents,i)=1;
|
|
wolffd@0
|
25 end
|
|
wolffd@0
|
26
|
|
wolffd@0
|
27 prior = pMax*rand(1,Ndiseases);
|
|
wolffd@0
|
28 leak = 0.5*rand(1,Nfindings); % in real QMR, leak approx exp(-0.02) = 0.98
|
|
wolffd@0
|
29 %leak = ones(1,Nfindings); % turns off leaks, which makes inference much harder
|
|
wolffd@0
|
30 inhibit = rand(Ndiseases, Nfindings);
|
|
wolffd@0
|
31 inhibit(not(G)) = 1;
|
|
wolffd@0
|
32
|
|
wolffd@0
|
33
|
|
wolffd@0
|
34 % first half of findings are +ve, second half -ve
|
|
wolffd@0
|
35 % The very first and last findings are hidden
|
|
wolffd@0
|
36 pos = 2:floor(Nfindings/2);
|
|
wolffd@0
|
37 neg = (pos(end)+1):(Nfindings-1);
|
|
wolffd@0
|
38
|
|
wolffd@0
|
39 % Make the bnet in the straightforward way
|
|
wolffd@0
|
40 tabular_leaves = 0;
|
|
wolffd@0
|
41 obs_nodes = myunion(pos, neg) + Ndiseases;
|
|
wolffd@0
|
42 big_bnet = mk_qmr_bnet(G, inhibit, leak, prior, tabular_leaves, obs_nodes);
|
|
wolffd@0
|
43 big_evidence = cell(1, N);
|
|
wolffd@0
|
44 big_evidence(findings(pos)) = num2cell(repmat(2, 1, length(pos)));
|
|
wolffd@0
|
45 big_evidence(findings(neg)) = num2cell(repmat(1, 1, length(neg)));
|
|
wolffd@0
|
46
|
|
wolffd@0
|
47 %clf;draw_layout(big_bnet.dag);
|
|
wolffd@0
|
48 %filename = '../public_html/Bayes/Figures/qmr.rnd.jpg';
|
|
wolffd@0
|
49 %% 3x3 inches
|
|
wolffd@0
|
50 %set(gcf,'units','inches');
|
|
wolffd@0
|
51 %set(gcf,'PaperPosition',[0 0 3 3])
|
|
wolffd@0
|
52 %print(gcf,'-djpeg','-r100',filename);
|
|
wolffd@0
|
53
|
|
wolffd@0
|
54
|
|
wolffd@0
|
55 % Marginalize out hidden leaves apriori
|
|
wolffd@0
|
56 positive_leaves_only = 1;
|
|
wolffd@0
|
57 [bnet, vals] = mk_minimal_qmr_bnet(G, inhibit, leak, prior, pos, neg, positive_leaves_only);
|
|
wolffd@0
|
58 obs_nodes = bnet.observed;
|
|
wolffd@0
|
59 evidence = cell(1, Ndiseases + length(obs_nodes));
|
|
wolffd@0
|
60 evidence(obs_nodes) = num2cell(vals);
|
|
wolffd@0
|
61
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 clear engine;
|
|
wolffd@0
|
64 engine{1} = quickscore_inf_engine(inhibit, leak, prior);
|
|
wolffd@0
|
65 engine{2} = jtree_inf_engine(big_bnet);
|
|
wolffd@0
|
66 engine{3} = jtree_inf_engine(bnet);
|
|
wolffd@0
|
67
|
|
wolffd@0
|
68 %fname = '/home/cs/murphyk/matlab/Misc/loopybel.txt';
|
|
wolffd@0
|
69 global BNT_HOME
|
|
wolffd@0
|
70 fname = sprintf('%s/loopybel.txt', BNT_HOME);
|
|
wolffd@0
|
71
|
|
wolffd@0
|
72
|
|
wolffd@0
|
73 max_iter = 6;
|
|
wolffd@0
|
74 engine{4} = pearl_inf_engine(bnet, 'protocol', 'parallel', 'max_iter', max_iter);
|
|
wolffd@0
|
75 %engine{5} = belprop_inf_engine(bnet, 'max_iter', max_iter, 'filename', fname);
|
|
wolffd@0
|
76 engine{5} = belprop_inf_engine(bnet, 'max_iter', max_iter);
|
|
wolffd@0
|
77
|
|
wolffd@0
|
78 E = length(engine);
|
|
wolffd@0
|
79 exact = 1:3;
|
|
wolffd@0
|
80 loopy = [4 5];
|
|
wolffd@0
|
81
|
|
wolffd@0
|
82 ll = zeros(1,E);
|
|
wolffd@0
|
83 tic; engine{1} = enter_evidence(engine{1}, pos, neg); toc
|
|
wolffd@0
|
84 tic; [engine{2}, ll(2)] = enter_evidence(engine{2}, big_evidence); toc
|
|
wolffd@0
|
85 tic; [engine{3}, ll(3)] = enter_evidence(engine{3}, evidence); toc
|
|
wolffd@0
|
86 tic; [engine{4}, ll(4), niter(4)] = enter_evidence(engine{4}, evidence); toc
|
|
wolffd@0
|
87 tic; [engine{5}, niter(5)] = enter_evidence(engine{5}, evidence); toc
|
|
wolffd@0
|
88
|
|
wolffd@0
|
89 ll
|
|
wolffd@0
|
90
|
|
wolffd@0
|
91 post = zeros(E, Ndiseases);
|
|
wolffd@0
|
92 for e=1:E
|
|
wolffd@0
|
93 for i=diseases(:)'
|
|
wolffd@0
|
94 m = marginal_nodes(engine{e}, i);
|
|
wolffd@0
|
95 post(e, i) = m.T(2);
|
|
wolffd@0
|
96 end
|
|
wolffd@0
|
97 end
|
|
wolffd@0
|
98
|
|
wolffd@0
|
99 for e=exact(:)'
|
|
wolffd@0
|
100 for i=diseases(:)'
|
|
wolffd@0
|
101 assert(approxeq(post(1, i), post(e, i)));
|
|
wolffd@0
|
102 end
|
|
wolffd@0
|
103 end
|
|
wolffd@0
|
104
|
|
wolffd@0
|
105 a = zeros(Ndiseases, 2);
|
|
wolffd@0
|
106 for ei=1:length(loopy)
|
|
wolffd@0
|
107 for i=diseases(:)'
|
|
wolffd@0
|
108 a(i,ei) = approxeq(post(1, i), post(loopy(ei), i));
|
|
wolffd@0
|
109 end
|
|
wolffd@0
|
110 end
|
|
wolffd@0
|
111 disp('is the loopy posterior correct?');
|
|
wolffd@0
|
112 disp(a)
|
