wolffd@0
|
1 % make a factor graph corresponding to an HMM with Gaussian outputs, where we absorb the
|
wolffd@0
|
2 % evidence up front
|
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 = 2;
|
wolffd@0
|
11 cts_obs = 1;
|
wolffd@0
|
12 param_tying = 1;
|
wolffd@0
|
13 bnet = mk_hmm_bnet(T, Q, O, cts_obs, param_tying);
|
wolffd@0
|
14 N = 2*T;
|
wolffd@0
|
15 onodes = bnet.observed;
|
wolffd@0
|
16 hnodes = mysetdiff(1:N, onodes);
|
wolffd@0
|
17
|
wolffd@0
|
18 data = sample_bnet(bnet);
|
wolffd@0
|
19
|
wolffd@0
|
20 init_factor = bnet.CPD{1};
|
wolffd@0
|
21 obs_factor = bnet.CPD{3};
|
wolffd@0
|
22 edge_factor = bnet.CPD{2}; % trans matrix
|
wolffd@0
|
23
|
wolffd@0
|
24 nfactors = T;
|
wolffd@0
|
25 nvars = T; % hidden only
|
wolffd@0
|
26 G = zeros(nvars, nfactors);
|
wolffd@0
|
27 G(1,1) = 1;
|
wolffd@0
|
28 for t=1:T-1
|
wolffd@0
|
29 G(t:t+1, t+1)=1;
|
wolffd@0
|
30 end
|
wolffd@0
|
31
|
wolffd@0
|
32 node_sizes = Q*ones(1,T);
|
wolffd@0
|
33
|
wolffd@0
|
34 % We tie params as follows:
|
wolffd@0
|
35 % the first hidden node use init_factor (number 1)
|
wolffd@0
|
36 % all hidden nodes on the backbone use edge_factor (number 2)
|
wolffd@0
|
37 % all observed nodes use the same factor, namely obs_factor
|
wolffd@0
|
38
|
wolffd@0
|
39 small_fg = mk_fgraph_given_ev(G, node_sizes, {init_factor, edge_factor}, {obs_factor}, data(onodes), ...
|
wolffd@0
|
40 'equiv_class', [1 2*ones(1,T-1)], 'ev_equiv_class', ones(1,T));
|
wolffd@0
|
41
|
wolffd@0
|
42 small_bnet = fgraph_to_bnet(small_fg);
|
wolffd@0
|
43
|
wolffd@0
|
44 % don't pre-process evidence
|
wolffd@0
|
45 % big_fg = bnet_to_fgraph(bnet); % can't handle Gaussian node
|
wolffd@0
|
46
|
wolffd@0
|
47
|
wolffd@0
|
48 engine = {};
|
wolffd@0
|
49 engine{1} = jtree_inf_engine(bnet);
|
wolffd@0
|
50 engine{2} = belprop_fg_inf_engine(small_fg, 'max_iter', 2*T);
|
wolffd@0
|
51 engine{3} = jtree_inf_engine(small_bnet);
|
wolffd@0
|
52 nengines = length(engine);
|
wolffd@0
|
53
|
wolffd@0
|
54
|
wolffd@0
|
55 % on BN, use the original evidence
|
wolffd@0
|
56 evidence = cell(1, 2*T);
|
wolffd@0
|
57 evidence(onodes) = data(onodes);
|
wolffd@0
|
58 tic; [engine{1}, ll(1)] = enter_evidence(engine{1}, evidence); toc
|
wolffd@0
|
59
|
wolffd@0
|
60
|
wolffd@0
|
61 % on small_fg, we have already included the evidence
|
wolffd@0
|
62 evidence = cell(1,T);
|
wolffd@0
|
63 tic; [engine{2}, ll(2)] = enter_evidence(engine{2}, evidence); toc
|
wolffd@0
|
64
|
wolffd@0
|
65
|
wolffd@0
|
66 % on small_bnet, we must add evidence to the dummy nodes
|
wolffd@0
|
67 V = small_fg.nvars;
|
wolffd@0
|
68 dummy = V+1:V+small_fg.nfactors;
|
wolffd@0
|
69 N = max(dummy);
|
wolffd@0
|
70 evidence = cell(1, N);
|
wolffd@0
|
71 evidence(dummy) = {1};
|
wolffd@0
|
72 tic; [engine{3}, ll(3)] = enter_evidence(engine{3}, evidence); toc
|
wolffd@0
|
73
|
wolffd@0
|
74
|
wolffd@0
|
75
|
wolffd@0
|
76 marg = zeros(T, nengines, Q); % marg(t,e,:)
|
wolffd@0
|
77 for t=1:T
|
wolffd@0
|
78 for e=1:nengines
|
wolffd@0
|
79 m = marginal_nodes(engine{e}, t);
|
wolffd@0
|
80 marg(t,e,:) = m.T;
|
wolffd@0
|
81 end
|
wolffd@0
|
82 end
|
wolffd@0
|
83 marg(:,:,1)
|