|
wolffd@0
|
1 % make an unrolled HMM, convert to factor graph, and check that
|
|
wolffd@0
|
2 % loopy propagation on the fgraph gives the exact answers.
|
|
wolffd@0
|
3
|
|
wolffd@0
|
4 seed = 1;
|
|
wolffd@0
|
5 rand('state', seed);
|
|
wolffd@0
|
6 randn('state', seed);
|
|
wolffd@0
|
7
|
|
wolffd@0
|
8 T = 3;
|
|
wolffd@0
|
9 Q = 3;
|
|
wolffd@0
|
10 O = 3;
|
|
wolffd@0
|
11 cts_obs = 0;
|
|
wolffd@0
|
12 param_tying = 1;
|
|
wolffd@0
|
13 bnet = mk_hmm_bnet(T, Q, O, cts_obs, param_tying);
|
|
wolffd@0
|
14
|
|
wolffd@0
|
15 data = sample_bnet(bnet);
|
|
wolffd@0
|
16
|
|
wolffd@0
|
17 fgraph = bnet_to_fgraph(bnet);
|
|
wolffd@0
|
18 big_bnet = fgraph_to_bnet(fgraph);
|
|
wolffd@0
|
19 % converting factor graph back does not recover the structure of the original bnet
|
|
wolffd@0
|
20
|
|
wolffd@0
|
21 max_iter = 2*T;
|
|
wolffd@0
|
22
|
|
wolffd@0
|
23 engine = {};
|
|
wolffd@0
|
24 engine{1} = jtree_inf_engine(bnet);
|
|
wolffd@0
|
25 engine{2} = belprop_inf_engine(bnet, 'max_iter', max_iter);
|
|
wolffd@0
|
26 engine{3} = belprop_fg_inf_engine(fgraph, 'max_iter', max_iter);
|
|
wolffd@0
|
27 engine{4} = jtree_inf_engine(big_bnet);
|
|
wolffd@0
|
28 nengines = length(engine);
|
|
wolffd@0
|
29
|
|
wolffd@0
|
30 big_engine = 4;
|
|
wolffd@0
|
31 fgraph_engine = 3;
|
|
wolffd@0
|
32
|
|
wolffd@0
|
33
|
|
wolffd@0
|
34 N = 2*T;
|
|
wolffd@0
|
35 evidence = cell(1,N);
|
|
wolffd@0
|
36 onodes = bnet.observed;
|
|
wolffd@0
|
37 evidence(onodes) = data(onodes);
|
|
wolffd@0
|
38 hnodes = mysetdiff(1:N, onodes);
|
|
wolffd@0
|
39
|
|
wolffd@0
|
40 bigN = length(big_bnet.dag);
|
|
wolffd@0
|
41 big_evidence = cell(1, bigN);
|
|
wolffd@0
|
42 big_evidence(onodes) = data(onodes);
|
|
wolffd@0
|
43 big_evidence(N+1:end) = {1}; % factors are observed to be 1
|
|
wolffd@0
|
44
|
|
wolffd@0
|
45 ll = zeros(1, nengines);
|
|
wolffd@0
|
46 for i=1:nengines
|
|
wolffd@0
|
47 if i==big_engine
|
|
wolffd@0
|
48 tic; [engine{i}, ll(i)] = enter_evidence(engine{i}, big_evidence); toc
|
|
wolffd@0
|
49 else
|
|
wolffd@0
|
50 tic; [engine{i}, ll(i)] = enter_evidence(engine{i}, evidence); toc
|
|
wolffd@0
|
51 end
|
|
wolffd@0
|
52 end
|
|
wolffd@0
|
53
|
|
wolffd@0
|
54 % compare all engines to engine{1}
|
|
wolffd@0
|
55
|
|
wolffd@0
|
56 % the log likelihood values may be bogus...
|
|
wolffd@0
|
57 for i=2:nengines
|
|
wolffd@0
|
58 %assert(approxeq(ll(1), ll(i)));
|
|
wolffd@0
|
59 end
|
|
wolffd@0
|
60
|
|
wolffd@0
|
61
|
|
wolffd@0
|
62 marg = zeros(T, nengines, Q); % marg(t,e,:)
|
|
wolffd@0
|
63 for t=1:T
|
|
wolffd@0
|
64 for e=1:nengines
|
|
wolffd@0
|
65 m = marginal_nodes(engine{e}, t);
|
|
wolffd@0
|
66 marg(t,e,:) = m.T;
|
|
wolffd@0
|
67 end
|
|
wolffd@0
|
68 end
|
|
wolffd@0
|
69 marg
|
|
wolffd@0
|
70
|
|
wolffd@0
|
71
|
|
wolffd@0
|
72 m = cell(nengines, T);
|
|
wolffd@0
|
73 for i=1:T
|
|
wolffd@0
|
74 for e=1:nengines
|
|
wolffd@0
|
75 m{e,i} = marginal_nodes(engine{e}, hnodes(i));
|
|
wolffd@0
|
76 end
|
|
wolffd@0
|
77 for e=2:nengines
|
|
wolffd@0
|
78 assert(approxeq(m{e,i}.T, m{1,i}.T));
|
|
wolffd@0
|
79 end
|
|
wolffd@0
|
80 end
|
|
wolffd@0
|
81
|
|
wolffd@0
|
82 mpe = {};
|
|
wolffd@0
|
83 ll = zeros(1, nengines);
|
|
wolffd@0
|
84 for e=1:nengines
|
|
wolffd@0
|
85 if e==big_engine
|
|
wolffd@0
|
86 mpe{e} = find_mpe(engine{e}, big_evidence);
|
|
wolffd@0
|
87 mpe{e} = mpe{e}(1:N); % chop off dummy nodes
|
|
wolffd@0
|
88 else
|
|
wolffd@0
|
89 mpe{e} = find_mpe(engine{e}, evidence);
|
|
wolffd@0
|
90 end
|
|
wolffd@0
|
91 end
|
|
wolffd@0
|
92
|
|
wolffd@0
|
93 % fgraph can't compute loglikelihood for software reasons
|
|
wolffd@0
|
94 % jtree on the big_bnet gives the wrong ll
|
|
wolffd@0
|
95 for e=2:nengines
|
|
wolffd@0
|
96 %assert(approxeq(ll(1), ll(e)));
|
|
wolffd@0
|
97 assert(approxeq(cell2num(mpe{1}), cell2num(mpe{e})))
|
|
wolffd@0
|
98 end
|
