samples.Sigma = samples.Sigma(:,:,burnin:end);
Nsamples = size(samples.mu, 2);

% Plot posterior on mu
[p, grid1, grid2] = ksdensity2d(samples.mu');
figure;
contour(grid1, grid2, p);
hold on
ndx = 1:1:Nsamples;
plot(samples.mu(1,ndx), samples.mu(2,ndx), '.');
line([35 65], [35 65], 'color', 'r', 'linewidth', 3);
%set(gca, 'xlim', mu1range);
%set(gca, 'ylim', mu2range);
xlabel(sprintf('%s', '\mu_1'))
ylabel(sprintf('%s', '\mu_2'))
title('posterior on \mu')
printPmtkFigure('gibbsGaussParamsPostMu')

quantileP