function marginal = marginal_from_samples(samples, weights, nodes, ns, dnodes, cnodes)
% MARGINAL_FROM_SAMPLES Compute an empiricial distribution from a weighted set of samples
% marginal = marginal_from_samples(samples, weights, nodes, node_sizes, discrete_nodes, cts_nodes)
%
% samples(i,:) is the i'th example
% nodes is the set of nodes we want to compute the distribution over

ddom = intersect(nodes, dnodes);
cdom = intersect(nodes, cnodes);
if subset(nodes, ddom)
  marginal = mk_table(nodes, ns(nodes));
  for j = 1:prod(ns(nodes))
    rows = findrows(samples(:, nodes), ind2subv(ns(nodes), j));
    marginal.T(j) = sum(weights(rows));
  end
  marginal.T = normalise(marginal.T);
elseif subset(nodes, cdom)
  marginal = mk_gpot(nodes, ns(nodes), 'm');
  [marginal.mu, marginal.Sigma] =  wstats(samples(:, nodes)', weights);
else 
  error('can''t handle mixed marginals yet');
end


%%%%%%%%

function [mu, Sigma] = wstats(X, w)

% Computes the weighted mean and weighted covariance matrix for a given
% set of observations X(:,i), and a set of normalised weights w(i).
% size(X, 1) = dimension of observation.
% size(X, 2) = number of observations.

d = X - repmat(X * w', 1, size(X, 2));
mu = sum(X .* repmat(w, size(X, 1), 1), 2);
Sigma = d * diag(w) * d';
