Daniel@0: % Do the example from Satnam Alag's PhD thesis, UCB ME dept 1996 p46
Daniel@0: 
Daniel@0: % Make the following polytree, where all arcs point down
Daniel@0: 
Daniel@0: % 1   2
Daniel@0: %  \ /
Daniel@0: %   3
Daniel@0: %  / \
Daniel@0: % 4   5
Daniel@0: 
Daniel@0: N = 5;
Daniel@0: dag = zeros(N,N);
Daniel@0: dag(1,3) = 1;
Daniel@0: dag(2,3) = 1;
Daniel@0: dag(3, [4 5]) = 1;
Daniel@0: 
Daniel@0: ns = [2 1 2 1 2];
Daniel@0: 
Daniel@0: bnet = mk_bnet(dag, ns, 'discrete', []);
Daniel@0: 
Daniel@0: bnet.CPD{1} = gaussian_CPD(bnet, 1, 'mean', [1 0]', 'cov', [4 1; 1 4]);
Daniel@0: bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', 1, 'cov', 1);
Daniel@0: B1 = [1 2; 1 0]; B2 = [2 1]';
Daniel@0: bnet.CPD{3} = gaussian_CPD(bnet, 3, 'mean', [0 0]', 'cov', [2 1; 1 1], ...
Daniel@0: 			   'weights', [B1 B2]);
Daniel@0: H1 = [1 1];
Daniel@0: bnet.CPD{4} = gaussian_CPD(bnet, 4, 'mean', 0, 'cov', 1, 'weights', H1);
Daniel@0: H2 = [1 0; 1 1];
Daniel@0: bnet.CPD{5} = gaussian_CPD(bnet, 5, 'mean', [0 0]', 'cov', eye(2), 'weights', H2);
Daniel@0: 
Daniel@0: engine = {};
Daniel@0: engine{end+1} = jtree_inf_engine(bnet);
Daniel@0: engine{end+1} = pearl_inf_engine(bnet, 'protocol', 'tree');
Daniel@0: engine{end+1} = pearl_inf_engine(bnet, 'protocol', 'parallel');
Daniel@0: E = length(engine);
Daniel@0: 
Daniel@0: if 1
Daniel@0: % no evidence
Daniel@0: evidence = cell(1,N);
Daniel@0: ll = zeros(1,E);
Daniel@0: for e=1:E
Daniel@0:   [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence);
Daniel@0:   add_ev = 1;
Daniel@0:   m = marginal_nodes(engine{e}, 3, add_ev);
Daniel@0:   assert(approxeq(m.mu, [3 2]'))
Daniel@0:   assert(approxeq(m.Sigma, [30 9; 9 6]))
Daniel@0: 
Daniel@0:   m = marginal_nodes(engine{e}, 4, add_ev);
Daniel@0:   assert(approxeq(m.mu, 5))
Daniel@0:   assert(approxeq(m.Sigma, 55))
Daniel@0: 
Daniel@0:   m = marginal_nodes(engine{e}, 5, add_ev);
Daniel@0:   assert(approxeq(m.mu, [3 5]'))
Daniel@0:   assert(approxeq(m.Sigma, [31 39; 39 55]))
Daniel@0: end
Daniel@0: end
Daniel@0: 
Daniel@0: if 1
Daniel@0: % evidence on leaf 5
Daniel@0: evidence = cell(1,N);
Daniel@0: evidence{5} = [5 5]';
Daniel@0: for e=1:E
Daniel@0:   [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence);
Daniel@0:   add_ev = 1;
Daniel@0:   m = marginal_nodes(engine{e}, 3, add_ev);
Daniel@0:   assert(approxeq(m.mu, [4.4022 1.0217]'))
Daniel@0:   assert(approxeq(m.Sigma, [0.7011 -0.4891; -0.4891 1.1087]))
Daniel@0: 
Daniel@0:   m = marginal_nodes(engine{e}, 4, add_ev);
Daniel@0:   assert(approxeq(m.mu, 5.4239))
Daniel@0:   assert(approxeq(m.Sigma, 1.8315))
Daniel@0: 
Daniel@0:   m = marginal_nodes(engine{e}, 1, add_ev);
Daniel@0:   assert(approxeq(m.mu, [0.3478 1.1413]'))
Daniel@0:   assert(approxeq(m.Sigma, [1.8261 -0.1957; -0.1957 1.0924]))
Daniel@0: 
Daniel@0:   m = marginal_nodes(engine{e}, 2, add_ev);
Daniel@0:   assert(approxeq(m.mu, 0.9239))
Daniel@0:   assert(approxeq(m.Sigma, 0.8315))
Daniel@0: 
Daniel@0:   m = marginal_nodes(engine{e}, 5, add_ev);
Daniel@0:   assert(approxeq(m.mu, evidence{5}))
Daniel@0:   assert(approxeq(m.Sigma, zeros(2)))
Daniel@0: end
Daniel@0: end
Daniel@0: 
Daniel@0: if 1
Daniel@0: % evidence on leaf 4 (non-info-state version is uninvertible)
Daniel@0: evidence = cell(1,N);
Daniel@0: evidence{4} = 10;
Daniel@0: for e=1:E
Daniel@0:   [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence);
Daniel@0:   add_ev = 1;
Daniel@0:   m = marginal_nodes(engine{e}, 3, add_ev);
Daniel@0:   assert(approxeq(m.mu, [6.5455 3.3636]'))
Daniel@0:   assert(approxeq(m.Sigma, [2.3455 -1.6364; -1.6364 1.9091]))
Daniel@0: 
Daniel@0:   m = marginal_nodes(engine{e}, 5, add_ev);
Daniel@0:   assert(approxeq(m.mu, [6.5455 9.9091]'))
Daniel@0:   assert(approxeq(m.Sigma, [3.3455 0.7091; 0.7091 1.9818]))
Daniel@0: 
Daniel@0:   m = marginal_nodes(engine{e}, 1, add_ev);
Daniel@0:   assert(approxeq(m.mu, [1.9091 0.9091]'))
Daniel@0:   assert(approxeq(m.Sigma, [2.1818 -0.8182; -0.8182 2.1818]))
Daniel@0: 
Daniel@0:   m = marginal_nodes(engine{e}, 2, add_ev);
Daniel@0:   assert(approxeq(m.mu, 1.2727))
Daniel@0:   assert(approxeq(m.Sigma, 0.8364))
Daniel@0: end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: if 1
Daniel@0: % evidence on leaves 4,5 and root 2
Daniel@0: evidence = cell(1,N);
Daniel@0: evidence{2} = 0;
Daniel@0: evidence{4} = 10;
Daniel@0: evidence{5} = [5 5]';
Daniel@0: for e=1:E
Daniel@0:   [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence);
Daniel@0:   add_ev = 1;
Daniel@0:   m = marginal_nodes(engine{e}, 3, add_ev);
Daniel@0:   assert(approxeq(m.mu, [4.9964 2.4444]'));
Daniel@0:   assert(approxeq(m.Sigma, [0.6738 -0.5556; -0.5556 0.8889]));
Daniel@0: 
Daniel@0:   m = marginal_nodes(engine{e}, 1, add_ev);
Daniel@0:   assert(approxeq(m.mu, [2.2043 1.2151]'));
Daniel@0:   assert(approxeq(m.Sigma, [1.2903 -0.4839; -0.4839 0.8065]));
Daniel@0: end
Daniel@0: end
Daniel@0: 
Daniel@0: if 1
Daniel@0:   [time, engine] = cmp_inference_static(bnet, engine, 'maximize', 0, 'check_ll', 0, ...
Daniel@0: 				     'singletons_only', 0, 'observed', [1 3 5]);
Daniel@0: end