|
wolffd@0
|
1 function m = marginal_family(engine, n, add_ev)
|
|
wolffd@0
|
2 % MARGINAL_FAMILY Compute the marginal on i's family (loopy)
|
|
wolffd@0
|
3 % m = marginal_family(engine, n, add_ev)
|
|
wolffd@0
|
4
|
|
wolffd@0
|
5 if nargin < 3, add_ev = 0; end
|
|
wolffd@0
|
6
|
|
wolffd@0
|
7 bnet = bnet_from_engine(engine);
|
|
wolffd@0
|
8 ns = bnet.node_sizes;
|
|
wolffd@0
|
9 ps = parents(bnet.dag, n);
|
|
wolffd@0
|
10 dom = [ps n];
|
|
wolffd@0
|
11 CPD = bnet.CPD{bnet.equiv_class(n)};
|
|
wolffd@0
|
12
|
|
wolffd@0
|
13 switch engine.msg_type
|
|
wolffd@0
|
14 case 'd',
|
|
wolffd@0
|
15 % The method is similar to the following HMM equation:
|
|
wolffd@0
|
16 % xi(i,j,t) = normalise( alpha(i,t) * transmat(i,j) * obsmat(j,t+1) * beta(j,t+1) )
|
|
wolffd@0
|
17 % where xi(i,j,t) = Pr(Q(t)=i, Q(t+1)=j | y(1:T))
|
|
wolffd@0
|
18 % beta == lambda, alpha == pi, alpha from each parent = pi msg
|
|
wolffd@0
|
19 % In general, if A,B are parents of C,
|
|
wolffd@0
|
20 % P(A,B,C) = P(C|A,B) pi_msg(A->C) pi_msg(B->C) lambda(C)
|
|
wolffd@0
|
21 % where lambda(C) = P(ev below and including C|C) = prod incoming lamba_msg(children->C)
|
|
wolffd@0
|
22 % and pi_msg(X->C) = P(X|ev above) etc
|
|
wolffd@0
|
23
|
|
wolffd@0
|
24 T = dpot(dom, ns(dom), CPD_to_CPT(CPD));
|
|
wolffd@0
|
25 for j=1:length(ps)
|
|
wolffd@0
|
26 p = ps(j);
|
|
wolffd@0
|
27 pi_msg = dpot(p, ns(p), engine.msg{n}.pi_from_parent{j});
|
|
wolffd@0
|
28 T = multiply_by_pot(T, pi_msg);
|
|
wolffd@0
|
29 end
|
|
wolffd@0
|
30 lambda = dpot(n, ns(n), engine.msg{n}.lambda);
|
|
wolffd@0
|
31 T = multiply_by_pot(T, lambda);
|
|
wolffd@0
|
32 T = normalize_pot(T);
|
|
wolffd@0
|
33 m = pot_to_marginal(T);
|
|
wolffd@0
|
34 if ~add_ev
|
|
wolffd@0
|
35 m.T = shrink_obs_dims_in_table(m.T, dom, engine.evidence);
|
|
wolffd@0
|
36 end
|
|
wolffd@0
|
37 case 'g',
|
|
wolffd@0
|
38 if engine.disconnected_nodes_bitv(n)
|
|
wolffd@0
|
39 m.T = 1;
|
|
wolffd@0
|
40 m.domain = dom;
|
|
wolffd@0
|
41 if add_ev
|
|
wolffd@0
|
42 m = add_ev_to_dmarginal(m, engine.evidence, ns)
|
|
wolffd@0
|
43 end
|
|
wolffd@0
|
44 return;
|
|
wolffd@0
|
45 end
|
|
wolffd@0
|
46
|
|
wolffd@0
|
47 [m, C, W] = gaussian_CPD_params_given_dps(CPD, dom, engine.evidence);
|
|
wolffd@0
|
48 cdom = myintersect(dom, bnet.cnodes);
|
|
wolffd@0
|
49 pot = linear_gaussian_to_cpot(m, C, W, dom, ns, cdom, engine.evidence);
|
|
wolffd@0
|
50 % linear_gaussian_to_cpot will set the effective size of observed nodes to 0,
|
|
wolffd@0
|
51 % so we need to do this explicitely for the messages, too,
|
|
wolffd@0
|
52 % so they are all the same size.
|
|
wolffd@0
|
53 obs_bitv = ~isemptycell(engine.evidence);
|
|
wolffd@0
|
54 ps = parents(engine.msg_dag, n);
|
|
wolffd@0
|
55 for j=1:length(ps)
|
|
wolffd@0
|
56 p = ps(j);
|
|
wolffd@0
|
57 msg = engine.msg{n}.pi_from_parent{j};
|
|
wolffd@0
|
58 if obs_bitv(p)
|
|
wolffd@0
|
59 pi_msg = mpot(p, 0);
|
|
wolffd@0
|
60 else
|
|
wolffd@0
|
61 pi_msg = mpot(p, ns(p), 0, msg.mu, msg.Sigma);
|
|
wolffd@0
|
62 end
|
|
wolffd@0
|
63 pot = multiply_by_pot(pot, mpot_to_cpot(pi_msg));
|
|
wolffd@0
|
64 end
|
|
wolffd@0
|
65 msg = engine.msg{n}.lambda;
|
|
wolffd@0
|
66 if obs_bitv(n)
|
|
wolffd@0
|
67 lambda = cpot(n, 0);
|
|
wolffd@0
|
68 else
|
|
wolffd@0
|
69 lambda = cpot(n, ns(n), 0, msg.info_state, msg.precision);
|
|
wolffd@0
|
70 end
|
|
wolffd@0
|
71 pot = multiply_by_pot(pot, lambda);
|
|
wolffd@0
|
72 m = pot_to_marginal(pot);
|
|
wolffd@0
|
73 if add_ev
|
|
wolffd@0
|
74 m = add_evidence_to_gmarginal(m, engine.evidence, bnet.node_sizes, bnet.cnodes);
|
|
wolffd@0
|
75 end
|
|
wolffd@0
|
76 end
|
|
wolffd@0
|
77
|
|
wolffd@0
|
78
|
|
wolffd@0
|
79
|
|
wolffd@0
|
80
|
