|
wolffd@0
|
1 function gibbs_test1()
|
|
wolffd@0
|
2
|
|
wolffd@0
|
3 disp('gibbs test 1')
|
|
wolffd@0
|
4
|
|
wolffd@0
|
5 rand('state', 0);
|
|
wolffd@0
|
6 randn('state', 0);
|
|
wolffd@0
|
7
|
|
wolffd@0
|
8 %[bnet onodes hnodes qnodes] = gibbs_ex_1;
|
|
wolffd@0
|
9 [bnet onodes hnodes qnodes] = gibbs_ex_2;
|
|
wolffd@0
|
10
|
|
wolffd@0
|
11 je = jtree_inf_engine(bnet);
|
|
wolffd@0
|
12 ge = gibbs_sampling_inf_engine (bnet, 'T', 50, 'burnin', 0, ...
|
|
wolffd@0
|
13 'order', [2 2 1 2 1]);
|
|
wolffd@0
|
14
|
|
wolffd@0
|
15 ev = sample_bnet(bnet);
|
|
wolffd@0
|
16
|
|
wolffd@0
|
17 evidence = cell(length(bnet.dag), 1);
|
|
wolffd@0
|
18 evidence(onodes) = ev(onodes);
|
|
wolffd@0
|
19 [je lj] = enter_evidence(je, evidence);
|
|
wolffd@0
|
20 [ge lg] = enter_evidence(ge, evidence);
|
|
wolffd@0
|
21
|
|
wolffd@0
|
22
|
|
wolffd@0
|
23 mj = marginal_nodes(je, qnodes);
|
|
wolffd@0
|
24
|
|
wolffd@0
|
25 [mg ge] = marginal_nodes (ge, qnodes);
|
|
wolffd@0
|
26 for t = 1:100
|
|
wolffd@0
|
27 [mg ge] = marginal_nodes (ge, qnodes, 'reset_counts', 0);
|
|
wolffd@0
|
28 diff = mj.T - mg.T;
|
|
wolffd@0
|
29 err(t) = norm (diff(:), 1);
|
|
wolffd@0
|
30 end
|
|
wolffd@0
|
31 clf
|
|
wolffd@0
|
32 plot(err);
|
|
wolffd@0
|
33 %title('error vs num. Gibbs samples')
|
|
wolffd@0
|
34
|
|
wolffd@0
|
35
|
|
wolffd@0
|
36 %%%%%%%
|
|
wolffd@0
|
37
|
|
wolffd@0
|
38 function [bnet, onodes, hnodes, qnodes] = gibbs_ex_1
|
|
wolffd@0
|
39 % bnet = gibbs_ex_1
|
|
wolffd@0
|
40 % a simple network to test the gibbs sampling engine
|
|
wolffd@0
|
41 % 1
|
|
wolffd@0
|
42 % / | \
|
|
wolffd@0
|
43 % 2 3 4
|
|
wolffd@0
|
44 % | | |
|
|
wolffd@0
|
45 % 5 6 7
|
|
wolffd@0
|
46 % \/ \/
|
|
wolffd@0
|
47 % 8 9
|
|
wolffd@0
|
48 % where all arcs point downwards
|
|
wolffd@0
|
49
|
|
wolffd@0
|
50 N = 9;
|
|
wolffd@0
|
51 dag = zeros(N,N);
|
|
wolffd@0
|
52 dag(1,2)=1; dag(1,3)=1; dag(1,4)=1;
|
|
wolffd@0
|
53 dag(2,5)=1; dag(3,6)=1; dag(4,7)=1;
|
|
wolffd@0
|
54 dag(5,8)=1; dag(6,8)=1; dag(6,9)=1; dag(7,9) = 1;
|
|
wolffd@0
|
55
|
|
wolffd@0
|
56 onodes = 8:9;
|
|
wolffd@0
|
57 hnodes = 1:7;
|
|
wolffd@0
|
58 qnodes = [1 2 6];
|
|
wolffd@0
|
59 ns = [2 3 4 3 5 2 4 3 2];
|
|
wolffd@0
|
60
|
|
wolffd@0
|
61 eclass = [1 2 3 2 4 5 6 7 8];
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 bnet = mk_bnet (dag, ns, 'equiv_class', eclass);
|
|
wolffd@0
|
64
|
|
wolffd@0
|
65 for i = 1:3
|
|
wolffd@0
|
66 bnet.CPD{i} = tabular_CPD(bnet, i);
|
|
wolffd@0
|
67 end
|
|
wolffd@0
|
68
|
|
wolffd@0
|
69 for i = 4:8
|
|
wolffd@0
|
70 bnet.CPD{i} = tabular_CPD(bnet, i+1);
|
|
wolffd@0
|
71 end
|
|
wolffd@0
|
72
|
|
wolffd@0
|
73
|
|
wolffd@0
|
74
|
|
wolffd@0
|
75 %%%%%%%
|
|
wolffd@0
|
76
|
|
wolffd@0
|
77 function [bnet, onodes, hnodes, qnodes] = gibbs_ex_2
|
|
wolffd@0
|
78 % bnet = gibbs_ex_2
|
|
wolffd@0
|
79 % a very simple network
|
|
wolffd@0
|
80 %
|
|
wolffd@0
|
81 % 1 2
|
|
wolffd@0
|
82 % \ /
|
|
wolffd@0
|
83 % 3
|
|
wolffd@0
|
84
|
|
wolffd@0
|
85 N = 3;
|
|
wolffd@0
|
86 dag = zeros(N,N);
|
|
wolffd@0
|
87 dag(1,3)=1; dag(2,3)=1;
|
|
wolffd@0
|
88
|
|
wolffd@0
|
89 onodes = 3;
|
|
wolffd@0
|
90 hnodes = 1:2;
|
|
wolffd@0
|
91 qnodes = 1:2;
|
|
wolffd@0
|
92 ns = [2 4 3];
|
|
wolffd@0
|
93
|
|
wolffd@0
|
94 eclass = [1 2 3];
|
|
wolffd@0
|
95
|
|
wolffd@0
|
96 bnet = mk_bnet (dag, ns, 'equiv_class', eclass);
|
|
wolffd@0
|
97
|
|
wolffd@0
|
98 for i = 1:3
|
|
wolffd@0
|
99 bnet.CPD{i} = tabular_CPD(bnet, i);
|
|
wolffd@0
|
100 end
|
|
wolffd@0
|
101
|
|
wolffd@0
|
102
|
|
wolffd@0
|
103
|
|
wolffd@0
|
104
|
|
wolffd@0
|
105
|
|
wolffd@0
|
106
|
|
wolffd@0
|
107
|
