wolffd@0
|
1 function [marginal, msg, loglik] = filter_evidence_old(engine, evidence)
|
wolffd@0
|
2 % [marginal, msg, loglik] = filter_evidence(engine, evidence) (pearl_dbn)
|
wolffd@0
|
3
|
wolffd@0
|
4 [ss T] = size(evidence);
|
wolffd@0
|
5 bnet = bnet_from_engine(engine);
|
wolffd@0
|
6 bnet2 = dbn_to_bnet(bnet, T);
|
wolffd@0
|
7 ns = bnet2.node_sizes;
|
wolffd@0
|
8 hnodes = mysetdiff(1:ss, engine.onodes);
|
wolffd@0
|
9 hnodes = hnodes(:)';
|
wolffd@0
|
10
|
wolffd@0
|
11 [engine.parent_index, engine.child_index] = mk_pearl_msg_indices(bnet2);
|
wolffd@0
|
12
|
wolffd@0
|
13 msg = init_msgs(bnet2.dag, ns, evidence);
|
wolffd@0
|
14 msg = init_ev_msgs(engine, evidence, msg);
|
wolffd@0
|
15
|
wolffd@0
|
16 verbose = 1;
|
wolffd@0
|
17 if verbose, fprintf('\nold filtering\n'); end
|
wolffd@0
|
18
|
wolffd@0
|
19 for t=1:T
|
wolffd@0
|
20 % update pi
|
wolffd@0
|
21 for i=hnodes
|
wolffd@0
|
22 n = i + (t-1)*ss;
|
wolffd@0
|
23 ps = parents(bnet2.dag, n);
|
wolffd@0
|
24 if t==1
|
wolffd@0
|
25 e = bnet.equiv_class(i,1);
|
wolffd@0
|
26 else
|
wolffd@0
|
27 e = bnet.equiv_class(i,2);
|
wolffd@0
|
28 end
|
wolffd@0
|
29 msg{n}.pi = compute_pi(bnet.CPD{e}, n, ps, msg);
|
wolffd@0
|
30 %if verbose, fprintf('%d computes pi\n', n); disp(msg{n}.pi); end
|
wolffd@0
|
31 msg{n}.pi = normalise(msg{n}.pi(:) .* msg{n}.lambda_from_self(:));
|
wolffd@0
|
32 if verbose, fprintf('%d recomputes pi\n', n); disp(msg{n}.pi); end
|
wolffd@0
|
33 end
|
wolffd@0
|
34 % send pi msg to children
|
wolffd@0
|
35 for i=hnodes
|
wolffd@0
|
36 n = i + (t-1)*ss;
|
wolffd@0
|
37 cs = children(bnet2.dag, n);
|
wolffd@0
|
38 for c=cs(:)'
|
wolffd@0
|
39 j = engine.parent_index{c}(n); % n is c's j'th parent
|
wolffd@0
|
40 pi_msg = normalise(compute_pi_msg(n, cs, msg, c, ns));
|
wolffd@0
|
41 msg{c}.pi_from_parent{j} = pi_msg;
|
wolffd@0
|
42 if verbose, fprintf('%d sends pi to %d\n', n,c); disp(pi_msg); end
|
wolffd@0
|
43 end
|
wolffd@0
|
44 end
|
wolffd@0
|
45 end
|
wolffd@0
|
46
|
wolffd@0
|
47
|
wolffd@0
|
48 marginal = cell(ss,T);
|
wolffd@0
|
49 lik = zeros(1,ss*T);
|
wolffd@0
|
50 for t=1:T
|
wolffd@0
|
51 for i=1:ss
|
wolffd@0
|
52 n = i + (t-1)*ss;
|
wolffd@0
|
53 %[bel, lik(n)] = normalise(msg{n}.pi .* msg{n}.lambda);
|
wolffd@0
|
54 [bel, lik(n)] = normalise(msg{n}.pi);
|
wolffd@0
|
55 marginal{i,t} = bel;
|
wolffd@0
|
56 end
|
wolffd@0
|
57 end
|
wolffd@0
|
58
|
wolffd@0
|
59 loglik = sum(log(lik));
|
wolffd@0
|
60
|
wolffd@0
|
61
|
wolffd@0
|
62
|
wolffd@0
|
63 %%%%%%%
|
wolffd@0
|
64
|
wolffd@0
|
65 function lambda = compute_lambda(n, cs, msg, ns)
|
wolffd@0
|
66 % Pearl p183 eq 4.50
|
wolffd@0
|
67 lambda = prod_lambda_msgs(n, cs, msg, ns);
|
wolffd@0
|
68
|
wolffd@0
|
69 %%%%%%%
|
wolffd@0
|
70
|
wolffd@0
|
71 function pi_msg = compute_pi_msg(n, cs, msg, c, ns)
|
wolffd@0
|
72 % Pearl p183 eq 4.53 and 4.51
|
wolffd@0
|
73 pi_msg = msg{n}.pi .* prod_lambda_msgs(n, cs, msg, ns, c);
|
wolffd@0
|
74
|
wolffd@0
|
75 %%%%%%%%%
|
wolffd@0
|
76
|
wolffd@0
|
77 function lam = prod_lambda_msgs(n, cs, msg, ns, except)
|
wolffd@0
|
78
|
wolffd@0
|
79 if nargin < 5, except = -1; end
|
wolffd@0
|
80
|
wolffd@0
|
81 %lam = msg{n}.lambda_from_self(:);
|
wolffd@0
|
82 lam = ones(ns(n), 1);
|
wolffd@0
|
83 for i=1:length(cs)
|
wolffd@0
|
84 c = cs(i);
|
wolffd@0
|
85 if c ~= except
|
wolffd@0
|
86 lam = lam .* msg{n}.lambda_from_child{i};
|
wolffd@0
|
87 end
|
wolffd@0
|
88 end
|
wolffd@0
|
89
|
wolffd@0
|
90
|
wolffd@0
|
91 %%%%%%%%%%%
|
wolffd@0
|
92
|
wolffd@0
|
93 function msg = init_msgs(dag, ns, evidence)
|
wolffd@0
|
94 % INIT_MSGS Initialize the lambda/pi message and state vectors (pearl_dbn)
|
wolffd@0
|
95 % msg = init_msgs(dag, ns, evidence)
|
wolffd@0
|
96 %
|
wolffd@0
|
97 % We assume all the hidden nodes are discrete.
|
wolffd@0
|
98
|
wolffd@0
|
99 N = length(dag);
|
wolffd@0
|
100 msg = cell(1,N);
|
wolffd@0
|
101 observed = ~isemptycell(evidence(:));
|
wolffd@0
|
102
|
wolffd@0
|
103 for n=1:N
|
wolffd@0
|
104 ps = parents(dag, n);
|
wolffd@0
|
105 msg{n}.pi_from_parent = cell(1, length(ps));
|
wolffd@0
|
106 for i=1:length(ps)
|
wolffd@0
|
107 p = ps(i);
|
wolffd@0
|
108 msg{n}.pi_from_parent{i} = ones(ns(p), 1);
|
wolffd@0
|
109 end
|
wolffd@0
|
110
|
wolffd@0
|
111 cs = children(dag, n);
|
wolffd@0
|
112 msg{n}.lambda_from_child = cell(1, length(cs));
|
wolffd@0
|
113 for i=1:length(cs)
|
wolffd@0
|
114 c = cs(i);
|
wolffd@0
|
115 msg{n}.lambda_from_child{i} = ones(ns(n), 1);
|
wolffd@0
|
116 end
|
wolffd@0
|
117
|
wolffd@0
|
118 msg{n}.lambda = ones(ns(n), 1);
|
wolffd@0
|
119 msg{n}.pi = ones(ns(n), 1);
|
wolffd@0
|
120
|
wolffd@0
|
121 msg{n}.lambda_from_self = ones(ns(n), 1);
|
wolffd@0
|
122 end
|
wolffd@0
|
123
|
wolffd@0
|
124
|
wolffd@0
|
125 %%%%%%%%%
|
wolffd@0
|
126
|
wolffd@0
|
127 function msg = init_ev_msgs(engine, evidence, msg)
|
wolffd@0
|
128 % Initialize the lambdas with any evidence
|
wolffd@0
|
129
|
wolffd@0
|
130 [ss T] = size(evidence);
|
wolffd@0
|
131 bnet = bnet_from_engine(engine);
|
wolffd@0
|
132 pot_type = 'd';
|
wolffd@0
|
133 t = 1;
|
wolffd@0
|
134 hnodes = mysetdiff(1:ss, engine.onodes);
|
wolffd@0
|
135 for i=hnodes(:)'
|
wolffd@0
|
136 c = engine.obschild(i);
|
wolffd@0
|
137 if c > 0
|
wolffd@0
|
138 fam = family(bnet.dag, c);
|
wolffd@0
|
139 e = bnet.equiv_class(c, 1);
|
wolffd@0
|
140 CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,1));
|
wolffd@0
|
141 temp = pot_to_marginal(CPDpot);
|
wolffd@0
|
142 n = i;
|
wolffd@0
|
143 msg{n}.lambda_from_self = temp.T;
|
wolffd@0
|
144 end
|
wolffd@0
|
145 end
|
wolffd@0
|
146 for t=2:T
|
wolffd@0
|
147 for i=hnodes(:)'
|
wolffd@0
|
148 c = engine.obschild(i);
|
wolffd@0
|
149 if c > 0
|
wolffd@0
|
150 fam = family(bnet.dag, c, 2);
|
wolffd@0
|
151 e = bnet.equiv_class(c, 2);
|
wolffd@0
|
152 CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,t-1:t));
|
wolffd@0
|
153 temp = pot_to_marginal(CPDpot);
|
wolffd@0
|
154 n = i + (t-1)*ss;
|
wolffd@0
|
155 msg{n}.lambda_from_self = temp.T;
|
wolffd@0
|
156 end
|
wolffd@0
|
157 end
|
wolffd@0
|
158 end
|