Mercurial > hg > camir-aes2014

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

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