annotate toolboxes/FullBNT-1.0.7/bnt/examples/static/Belprop/belprop_polytree_gauss.m @ 0:cc4b1211e677 tip

initial commit to HG from Changeset: 646 (e263d8a21543) added further path and more save "camirversion.m"
author Daniel Wolff
date Fri, 19 Aug 2016 13:07:06 +0200
parents
children
rev   line source
Daniel@0 1 % Do the example from Satnam Alag's PhD thesis, UCB ME dept 1996 p46
Daniel@0 2
Daniel@0 3 % Make the following polytree, where all arcs point down
Daniel@0 4
Daniel@0 5 % 1 2
Daniel@0 6 % \ /
Daniel@0 7 % 3
Daniel@0 8 % / \
Daniel@0 9 % 4 5
Daniel@0 10
Daniel@0 11 N = 5;
Daniel@0 12 dag = zeros(N,N);
Daniel@0 13 dag(1,3) = 1;
Daniel@0 14 dag(2,3) = 1;
Daniel@0 15 dag(3, [4 5]) = 1;
Daniel@0 16
Daniel@0 17 ns = [2 1 2 1 2];
Daniel@0 18
Daniel@0 19 bnet = mk_bnet(dag, ns, 'discrete', []);
Daniel@0 20
Daniel@0 21 bnet.CPD{1} = gaussian_CPD(bnet, 1, 'mean', [1 0]', 'cov', [4 1; 1 4]);
Daniel@0 22 bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', 1, 'cov', 1);
Daniel@0 23 B1 = [1 2; 1 0]; B2 = [2 1]';
Daniel@0 24 bnet.CPD{3} = gaussian_CPD(bnet, 3, 'mean', [0 0]', 'cov', [2 1; 1 1], ...
Daniel@0 25 'weights', [B1 B2]);
Daniel@0 26 H1 = [1 1];
Daniel@0 27 bnet.CPD{4} = gaussian_CPD(bnet, 4, 'mean', 0, 'cov', 1, 'weights', H1);
Daniel@0 28 H2 = [1 0; 1 1];
Daniel@0 29 bnet.CPD{5} = gaussian_CPD(bnet, 5, 'mean', [0 0]', 'cov', eye(2), 'weights', H2);
Daniel@0 30
Daniel@0 31 engine = {};
Daniel@0 32 engine{end+1} = jtree_inf_engine(bnet);
Daniel@0 33 engine{end+1} = pearl_inf_engine(bnet, 'protocol', 'tree');
Daniel@0 34 engine{end+1} = pearl_inf_engine(bnet, 'protocol', 'parallel');
Daniel@0 35 E = length(engine);
Daniel@0 36
Daniel@0 37 if 1
Daniel@0 38 % no evidence
Daniel@0 39 evidence = cell(1,N);
Daniel@0 40 ll = zeros(1,E);
Daniel@0 41 for e=1:E
Daniel@0 42 [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence);
Daniel@0 43 add_ev = 1;
Daniel@0 44 m = marginal_nodes(engine{e}, 3, add_ev);
Daniel@0 45 assert(approxeq(m.mu, [3 2]'))
Daniel@0 46 assert(approxeq(m.Sigma, [30 9; 9 6]))
Daniel@0 47
Daniel@0 48 m = marginal_nodes(engine{e}, 4, add_ev);
Daniel@0 49 assert(approxeq(m.mu, 5))
Daniel@0 50 assert(approxeq(m.Sigma, 55))
Daniel@0 51
Daniel@0 52 m = marginal_nodes(engine{e}, 5, add_ev);
Daniel@0 53 assert(approxeq(m.mu, [3 5]'))
Daniel@0 54 assert(approxeq(m.Sigma, [31 39; 39 55]))
Daniel@0 55 end
Daniel@0 56 end
Daniel@0 57
Daniel@0 58 if 1
Daniel@0 59 % evidence on leaf 5
Daniel@0 60 evidence = cell(1,N);
Daniel@0 61 evidence{5} = [5 5]';
Daniel@0 62 for e=1:E
Daniel@0 63 [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence);
Daniel@0 64 add_ev = 1;
Daniel@0 65 m = marginal_nodes(engine{e}, 3, add_ev);
Daniel@0 66 assert(approxeq(m.mu, [4.4022 1.0217]'))
Daniel@0 67 assert(approxeq(m.Sigma, [0.7011 -0.4891; -0.4891 1.1087]))
Daniel@0 68
Daniel@0 69 m = marginal_nodes(engine{e}, 4, add_ev);
Daniel@0 70 assert(approxeq(m.mu, 5.4239))
Daniel@0 71 assert(approxeq(m.Sigma, 1.8315))
Daniel@0 72
Daniel@0 73 m = marginal_nodes(engine{e}, 1, add_ev);
Daniel@0 74 assert(approxeq(m.mu, [0.3478 1.1413]'))
Daniel@0 75 assert(approxeq(m.Sigma, [1.8261 -0.1957; -0.1957 1.0924]))
Daniel@0 76
Daniel@0 77 m = marginal_nodes(engine{e}, 2, add_ev);
Daniel@0 78 assert(approxeq(m.mu, 0.9239))
Daniel@0 79 assert(approxeq(m.Sigma, 0.8315))
Daniel@0 80
Daniel@0 81 m = marginal_nodes(engine{e}, 5, add_ev);
Daniel@0 82 assert(approxeq(m.mu, evidence{5}))
Daniel@0 83 assert(approxeq(m.Sigma, zeros(2)))
Daniel@0 84 end
Daniel@0 85 end
Daniel@0 86
Daniel@0 87 if 1
Daniel@0 88 % evidence on leaf 4 (non-info-state version is uninvertible)
Daniel@0 89 evidence = cell(1,N);
Daniel@0 90 evidence{4} = 10;
Daniel@0 91 for e=1:E
Daniel@0 92 [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence);
Daniel@0 93 add_ev = 1;
Daniel@0 94 m = marginal_nodes(engine{e}, 3, add_ev);
Daniel@0 95 assert(approxeq(m.mu, [6.5455 3.3636]'))
Daniel@0 96 assert(approxeq(m.Sigma, [2.3455 -1.6364; -1.6364 1.9091]))
Daniel@0 97
Daniel@0 98 m = marginal_nodes(engine{e}, 5, add_ev);
Daniel@0 99 assert(approxeq(m.mu, [6.5455 9.9091]'))
Daniel@0 100 assert(approxeq(m.Sigma, [3.3455 0.7091; 0.7091 1.9818]))
Daniel@0 101
Daniel@0 102 m = marginal_nodes(engine{e}, 1, add_ev);
Daniel@0 103 assert(approxeq(m.mu, [1.9091 0.9091]'))
Daniel@0 104 assert(approxeq(m.Sigma, [2.1818 -0.8182; -0.8182 2.1818]))
Daniel@0 105
Daniel@0 106 m = marginal_nodes(engine{e}, 2, add_ev);
Daniel@0 107 assert(approxeq(m.mu, 1.2727))
Daniel@0 108 assert(approxeq(m.Sigma, 0.8364))
Daniel@0 109 end
Daniel@0 110 end
Daniel@0 111
Daniel@0 112
Daniel@0 113 if 1
Daniel@0 114 % evidence on leaves 4,5 and root 2
Daniel@0 115 evidence = cell(1,N);
Daniel@0 116 evidence{2} = 0;
Daniel@0 117 evidence{4} = 10;
Daniel@0 118 evidence{5} = [5 5]';
Daniel@0 119 for e=1:E
Daniel@0 120 [engine{e}, ll(e)] = enter_evidence(engine{e}, evidence);
Daniel@0 121 add_ev = 1;
Daniel@0 122 m = marginal_nodes(engine{e}, 3, add_ev);
Daniel@0 123 assert(approxeq(m.mu, [4.9964 2.4444]'));
Daniel@0 124 assert(approxeq(m.Sigma, [0.6738 -0.5556; -0.5556 0.8889]));
Daniel@0 125
Daniel@0 126 m = marginal_nodes(engine{e}, 1, add_ev);
Daniel@0 127 assert(approxeq(m.mu, [2.2043 1.2151]'));
Daniel@0 128 assert(approxeq(m.Sigma, [1.2903 -0.4839; -0.4839 0.8065]));
Daniel@0 129 end
Daniel@0 130 end
Daniel@0 131
Daniel@0 132 if 1
Daniel@0 133 [time, engine] = cmp_inference_static(bnet, engine, 'maximize', 0, 'check_ll', 0, ...
Daniel@0 134 'singletons_only', 0, 'observed', [1 3 5]);
Daniel@0 135 end