Daniel@0: % make a factor graph corresponding to an  HMM with Gaussian outputs, where we absorb the
Daniel@0: % evidence up front 
Daniel@0: 
Daniel@0: seed = 1;
Daniel@0: rand('state', seed);
Daniel@0: randn('state', seed);
Daniel@0: 
Daniel@0: T = 3;
Daniel@0: Q = 3;
Daniel@0: O = 2;
Daniel@0: cts_obs = 1;
Daniel@0: param_tying = 1;
Daniel@0: bnet = mk_hmm_bnet(T, Q, O, cts_obs, param_tying);
Daniel@0: N = 2*T;
Daniel@0: onodes = bnet.observed;
Daniel@0: hnodes = mysetdiff(1:N, onodes);
Daniel@0: 
Daniel@0: data = sample_bnet(bnet);
Daniel@0: 
Daniel@0: init_factor = bnet.CPD{1};
Daniel@0: obs_factor = bnet.CPD{3};
Daniel@0: edge_factor = bnet.CPD{2}; % trans matrix
Daniel@0: 
Daniel@0: nfactors = T;
Daniel@0: nvars = T; % hidden only
Daniel@0: G = zeros(nvars, nfactors);
Daniel@0: G(1,1) = 1;
Daniel@0: for t=1:T-1
Daniel@0:   G(t:t+1, t+1)=1;
Daniel@0: end
Daniel@0: 
Daniel@0: node_sizes = Q*ones(1,T);
Daniel@0: 
Daniel@0: % We tie params as follows:
Daniel@0: % the first hidden node use init_factor (number 1)
Daniel@0: % all hidden nodes on the backbone use edge_factor (number 2)
Daniel@0: % all observed nodes use the same factor, namely obs_factor
Daniel@0: 
Daniel@0: small_fg = mk_fgraph_given_ev(G, node_sizes, {init_factor, edge_factor}, {obs_factor}, data(onodes), ...
Daniel@0: 			 'equiv_class', [1 2*ones(1,T-1)], 'ev_equiv_class', ones(1,T));
Daniel@0: 
Daniel@0: small_bnet = fgraph_to_bnet(small_fg);
Daniel@0: 
Daniel@0: % don't pre-process evidence
Daniel@0: % big_fg = bnet_to_fgraph(bnet); % can't handle Gaussian node
Daniel@0: 
Daniel@0: 
Daniel@0: engine = {};
Daniel@0: engine{1} = jtree_inf_engine(bnet);
Daniel@0: engine{2} = belprop_fg_inf_engine(small_fg, 'max_iter', 2*T);
Daniel@0: engine{3} = jtree_inf_engine(small_bnet);
Daniel@0: nengines = length(engine);
Daniel@0: 
Daniel@0: 
Daniel@0: % on BN, use the original evidence
Daniel@0: evidence = cell(1, 2*T);
Daniel@0: evidence(onodes) = data(onodes);
Daniel@0: tic; [engine{1}, ll(1)] = enter_evidence(engine{1}, evidence); toc
Daniel@0: 
Daniel@0: 
Daniel@0: % on small_fg, we have already included the evidence
Daniel@0: evidence = cell(1,T);
Daniel@0: tic; [engine{2}, ll(2)] = enter_evidence(engine{2}, evidence); toc
Daniel@0: 
Daniel@0: 
Daniel@0: % on small_bnet, we must add evidence to the dummy nodes 
Daniel@0: V = small_fg.nvars;
Daniel@0: dummy = V+1:V+small_fg.nfactors;
Daniel@0: N = max(dummy);
Daniel@0: evidence = cell(1, N);
Daniel@0: evidence(dummy) = {1};
Daniel@0: tic; [engine{3}, ll(3)] = enter_evidence(engine{3}, evidence); toc
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: marg = zeros(T, nengines, Q); % marg(t,e,:)
Daniel@0: for t=1:T
Daniel@0:   for e=1:nengines
Daniel@0:     m = marginal_nodes(engine{e}, t);
Daniel@0:     marg(t,e,:) = m.T;
Daniel@0:   end
Daniel@0: end
Daniel@0: marg(:,:,1)