MTK( samples.mu(2,:) - samples.mu(1,:), [0.025 0.5 0.975])
mean(samples.mu(2,:) > samples.mu(1,:))


% Compute posterior predictive
Xpred = zeros(2,Nsamples);
for i=1:Nsamples
  x = gaussSample(samples.mu(:,i), samples.Sigma(:,:,i), 1);
  Xpred(:,i) = colvec(x);
end

% Now plot posterior predictive
[p, xgrid1, xgrid2] = ksdensity2d(Xpred');
figure;
contour(xgrid1, xgrid2, p)
hold on
plot(Xpred(1,:), Xpred(2,:), '.');
plot(data(:,1), data(:,2), 'x', 'color', 'r', ...
  'markersize', 8, 'linewidth', 2);
line(