matthiasm@8: function mcmc_post = mcmc_sample_to_hist(sampled_graphs, dags)
matthiasm@8: % MCMC_SAMPLE_TO_HIST Convert a set of sampled dags into a histogram over dags
matthiasm@8: % hist = mcmc_sample_to_hist(sampled_graphs, dags)
matthiasm@8: %
matthiasm@8: % sampled_graphs{m} is the m'th sampled dag
matthiasm@8: % dags{i} is the i'th dag in the hypothesis space
matthiasm@8: % hist(i) = Pr(model i | data)
matthiasm@8: 
matthiasm@8: ndags = length(dags);
matthiasm@8: nsamples = length(sampled_graphs);
matthiasm@8: nnodes = length(dags{1});
matthiasm@8: % sampled_bitv(m, :) is the m'th sampled graph represented as a vector of n^2 bits, computed
matthiasm@8: % by stacking the columns of the adjacency matrix vertically.
matthiasm@8: sampled_bitvs = zeros(nsamples, nnodes*nnodes);
matthiasm@8: for m=1:nsamples
matthiasm@8:   sampled_bitvs(m, :) = sampled_graphs{m}(:)';
matthiasm@8: end
matthiasm@8:   
matthiasm@8: [ugraphs, I, J] = unique(sampled_bitvs, 'rows');  % each row of ugraphs is a unique bit vector
matthiasm@8: sampled_indices  = subv2ind(2*ones(1,nnodes*nnodes), ugraphs+1);
matthiasm@8: counts = hist(J, 1:size(ugraphs,1)); % counts(i) = number of times graphs(i,:) occurs in the sample
matthiasm@8: 
matthiasm@8: mcmc_post = zeros(1, ndags);
matthiasm@8: for i=1:ndags
matthiasm@8:   bitv = dags{i}(:)';
matthiasm@8:   % Find the samples that corresponds to this graph by converting the graphs to bitvectors and
matthiasm@8:   % then to integers.
matthiasm@8:   ndx = subv2ind(2*ones(1,nnodes*nnodes), bitv+1);
matthiasm@8:   locn = find(ndx == sampled_indices);
matthiasm@8:   if ~isempty(locn)
matthiasm@8:     mcmc_post(i) = counts(locn);
matthiasm@8:   end
matthiasm@8: end
matthiasm@8: mcmc_post = normalise(mcmc_post);
matthiasm@8: 
matthiasm@8: