wolffd@0: function marginal = marginal_nodes(engine, nodes)
wolffd@0: % MARGINAL_NODES Compute the marginal on the specified query nodes (likelihood_weighting)
wolffd@0: % marginal = marginal_nodes(engine, nodes)
wolffd@0: 
wolffd@0: bnet = bnet_from_engine(engine);
wolffd@0: ddom = myintersect(nodes, bnet.dnodes);
wolffd@0: cdom = myintersect(nodes, bnet.cnodes);
wolffd@0: nsamples = size(engine.samples, 1);
wolffd@0: ns = bnet.node_sizes;
wolffd@0: 
wolffd@0: %w = normalise(engine.weights);
wolffd@0: w = engine.weights;
wolffd@0: if mysubset(nodes, ddom)
wolffd@0:   T = 0*myones(ns(nodes));
wolffd@0:   P = prod(ns(nodes));
wolffd@0:   indices = ind2subv(ns(nodes), 1:P);
wolffd@0:   samples = reshape(cat(1, engine.samples{:,nodes}), nsamples, length(nodes));
wolffd@0:   for j = 1:P
wolffd@0:     rows = find_rows(samples, indices(j,:));
wolffd@0:     T(j) = sum(w(rows));
wolffd@0:   end
wolffd@0:   T = normalise(T);
wolffd@0:   marginal.T = T;
wolffd@0: elseif subset(nodes, cdom)
wolffd@0:   samples = reshape(cat(1, engine.samples{:,nodes}), nsamples*sum(ns(nodes)), length(nodes));
wolffd@0:   [marginal.mu, marginal.Sigma] =  wstats(samples', normalise(w));
wolffd@0: else
wolffd@0:   error('can''t handle mixed marginals yet');
wolffd@0: end
wolffd@0: 
wolffd@0: marginal.domain = nodes;
wolffd@0: 
wolffd@0: %%%%%%%%%
wolffd@0: 
wolffd@0: function rows = find_rows(M, v)
wolffd@0: % FINDROWS Find rows which are equal to a specified vector
wolffd@0: % rows = findrows(M, v)
wolffd@0: % Each row of M is a sample
wolffd@0: 
wolffd@0: temp = abs(M - repmat(v, size(M, 1), 1));
wolffd@0: rows = find(sum(temp,2) == 0);      
wolffd@0: 
wolffd@0: %%%%%%%%
wolffd@0: 
wolffd@0: function [mu, Sigma] = wstats(X, w)
wolffd@0: 
wolffd@0: % Computes the weighted mean and weighted covariance matrix for a given
wolffd@0: % set of observations X(:,i), and a set of normalised weights w(i).
wolffd@0: % Each column of X is a sample.
wolffd@0: 
wolffd@0: d = X - repmat(X * w', 1, size(X, 2));
wolffd@0: mu = sum(X .* repmat(w, size(X, 1), 1), 2);
wolffd@0: Sigma = d * diag(w) * d';