% Plot
figure(2);
clf;

% Choose colormap.
cmap = colormap(hot);
%cmap = [linspace(1,0,100)' linspace(0,0,100)' linspace(0,0,100)'; linspace(0,0,100)' linspace(0,0,100)' linspace(0,1,100)'];
colormap(cmap);

% find combinations of Nf,Nr that we've tested
indices = find(mask == 1);

% one sigma
c = 1;
ci = cis(c); % expected fraction that should fall in this confidence interval

% find scale limits
X = squeeze(Pfrs_asymptotic(:,:,c));
Y = squeeze(Pfrs_bayesian(:,:,c));
%max_dev = max(max([X(indices)-1 Y(indices)-1 1-X(indices) 1-Y(indices)]))

subplot(2,3,1);
show_deviation(X, [0 1], mask, cmap);
title(sprintf('ABAR rel fraction for %.2f CI', ci));

subplot(2,3,4);
show_deviation(Y, [0 1], mask, cmap);
title(sprintf('BBAR rel fraction for %.2f CI', ci));

% TWO SIGMA
c = 2;
ci = cis(c); % expected fraction that should fall in this confidence interval

% find scale limits
X = squeeze(Pfrs_asymptotic(:,:,c));
Y = squeeze(Pfrs_bayesian(:,:,c));
%max_dev = max(max([X(indices)-1 Y(indices)-1 1-X(indices) 1-Y(indices)]))

subplot(2,3,2);
show_deviation(X, [0 1], mask, cmap);
title(sprintf('ABAR rel fraction for %.2f CI', ci));

subplot(2,3,5);
show_deviation(Y, [0 1], mask, cmap);
title(sprintf('BBAR rel fraction for %.2f CI', ci));

% BIAS

% find scale limits
X = squeeze(BAR_ML_bias_fr / abs(true_df));
Y = squeeze(BAR_mean_bias_fr / abs(true_df));
max_dev = max(max([X(indices)-1 Y(indices)-1 1-X(indices) 1-Y(indices)]))
max_dev = 2;

subplot(2,3,3);
show_deviation(X, [-max_dev +max_dev], mask, cmap);
title(sprintf('ABAR rel bias of ML'));

subplot(2,3,6);
show_deviation(Y, [-max_dev +max_dev], mask, cmap);
title(sprintf('BBAR rel bias of posterior mean'));

% print
print -depsc bar-data.eps
unix('epstopdf bar-data.eps');