wolffd@0: function m = marginal_family(engine, n, add_ev)
wolffd@0: % MARGINAL_FAMILY Compute the marginal on i's family (loopy)
wolffd@0: % m = marginal_family(engine, n, add_ev)
wolffd@0: 
wolffd@0: if nargin < 3, add_ev = 0; end
wolffd@0: 
wolffd@0: bnet = bnet_from_engine(engine);
wolffd@0: ns = bnet.node_sizes;
wolffd@0: ps = parents(bnet.dag, n);
wolffd@0: dom = [ps n];
wolffd@0: CPD = bnet.CPD{bnet.equiv_class(n)};
wolffd@0: 
wolffd@0: switch engine.msg_type
wolffd@0:   case 'd',
wolffd@0:    % The method is similar to the following HMM equation:
wolffd@0:    % xi(i,j,t) = normalise( alpha(i,t) * transmat(i,j) * obsmat(j,t+1) * beta(j,t+1) )
wolffd@0:    % where xi(i,j,t) = Pr(Q(t)=i, Q(t+1)=j | y(1:T))   
wolffd@0:    % beta == lambda, alpha == pi, alpha from each parent = pi msg
wolffd@0:    % In general, if A,B are parents of C,
wolffd@0:    % P(A,B,C) = P(C|A,B) pi_msg(A->C) pi_msg(B->C) lambda(C)
wolffd@0:    % where lambda(C) = P(ev below and including C|C) = prod incoming lamba_msg(children->C)
wolffd@0:    % and pi_msg(X->C) = P(X|ev above) etc
wolffd@0:    
wolffd@0:    T = dpot(dom, ns(dom), CPD_to_CPT(CPD));
wolffd@0:    for j=1:length(ps)
wolffd@0:      p = ps(j);
wolffd@0:      pi_msg = dpot(p, ns(p), engine.msg{n}.pi_from_parent{j});
wolffd@0:      T = multiply_by_pot(T, pi_msg);
wolffd@0:    end         
wolffd@0:    lambda = dpot(n, ns(n), engine.msg{n}.lambda);
wolffd@0:    T = multiply_by_pot(T, lambda);
wolffd@0:    T = normalize_pot(T);
wolffd@0:    m = pot_to_marginal(T);
wolffd@0:    if ~add_ev
wolffd@0:      m.T = shrink_obs_dims_in_table(m.T, dom, engine.evidence);
wolffd@0:    end
wolffd@0:  case 'g',
wolffd@0:   if engine.disconnected_nodes_bitv(n)
wolffd@0:     m.T = 1;
wolffd@0:     m.domain = dom;
wolffd@0:     if add_ev
wolffd@0:       m = add_ev_to_dmarginal(m, engine.evidence, ns)
wolffd@0:     end
wolffd@0:     return;
wolffd@0:   end
wolffd@0: 
wolffd@0:   [m, C, W] = gaussian_CPD_params_given_dps(CPD, dom, engine.evidence);
wolffd@0:   cdom = myintersect(dom, bnet.cnodes);
wolffd@0:   pot = linear_gaussian_to_cpot(m, C, W, dom, ns, cdom, engine.evidence); 
wolffd@0:   % linear_gaussian_to_cpot will set the effective size of observed nodes to 0,
wolffd@0:   % so we need to do this explicitely for the messages, too,
wolffd@0:   % so they are all the same size.
wolffd@0:   obs_bitv = ~isemptycell(engine.evidence);
wolffd@0:   ps = parents(engine.msg_dag, n);
wolffd@0:   for j=1:length(ps)
wolffd@0:     p = ps(j);
wolffd@0:     msg = engine.msg{n}.pi_from_parent{j};
wolffd@0:     if obs_bitv(p)
wolffd@0:       pi_msg = mpot(p, 0);
wolffd@0:     else
wolffd@0:       pi_msg = mpot(p, ns(p), 0, msg.mu, msg.Sigma);
wolffd@0:     end
wolffd@0:     pot = multiply_by_pot(pot, mpot_to_cpot(pi_msg));
wolffd@0:   end         
wolffd@0:   msg = engine.msg{n}.lambda;
wolffd@0:   if obs_bitv(n)
wolffd@0:     lambda = cpot(n, 0);
wolffd@0:   else
wolffd@0:     lambda = cpot(n, ns(n), 0, msg.info_state, msg.precision);
wolffd@0:   end
wolffd@0:   pot = multiply_by_pot(pot, lambda);
wolffd@0:   m = pot_to_marginal(pot);
wolffd@0:   if add_ev
wolffd@0:     m = add_evidence_to_gmarginal(m, engine.evidence, bnet.node_sizes, bnet.cnodes);
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: