% Bayes_2D_HPD_SR.m
% This routine returns desired (1-alpha)
% Highest Probability Density (HPD) regions
% given the results of the 2-D marginal
% posterior density calculator,
% Bayes_MCMC_2Dmarginal_SR.m
%
% Routine inputs:
% bin_ic = the center coordinates of each bin in the
%          i_plot histogram
% bin_jc = the center coordinates of each bin in the
%          j_plot histogram
% bin_p = the marginal probability values for each bin,
%         where bin_p(m,n) contains the marginal density
%         at theta(i_plot) = bin_ic(m),
%         theta(j_plot_ = bin_jc(n).
% i_plot = the label of the first parameter
% j_plot = the label of the second parameter
% alpha : for 95% CI, use alpha = 0.05. The HPD is
%         returned that has (1-alpha) fraction of
%         total marginal density
% Routine outputs:
% HPD_2D = contains a 1 for each bin if it is in the HPD, else
%       it contains a zero
%
% K. J. Beers. MIT ChE. 12/13/2004. ver. 12/13/2004

function HPD_2D = Bayes_2D_HPD_SR(bin_2Dic, bin_2Djc, bin_2Dp, ...
    i_plot, j_plot, alpha);

% extract dimension parameters
N_bins = length(bin_2Dic);

% set area of each bin
bin_delta = zeros(2,1);
bin_delta(1) = bin_2Dic(2) - bin_2Dic(1);
bin_delta(2) = bin_2Djc(2) - bin_2Djc(1);
bin_area = bin_delta(1)*bin_delta(2);

% initialize HPD output structure
HPD_2D = zeros(N_bins,N_bins);

% set vector of trial probability contour values
p_c_max = max(max(bin_2Dp));
p_c_vals = linspace(0,p_c_max,N_bins);
alpha_vals = zeros(size(p_c_vals));

% for each contour value, compute total density
% within corresponding HPD
for n=1:N_bins
    p_c = p_c_vals(n);
    % compute net fraction of probability above this
    % contour value
    for ki=1:N_bins
        for kj=1:N_bins
            if(bin_2Dp(ki,kj) >= p_c)
                alpha_vals(n) = alpha_vals(n) + ...
                    bin_2Dp(ki,kj);
            end
        end
    end
    alpha_vals(n) = alpha_vals(n)*bin_area;
end
% make plot of alpha vs. p_c
figure;
plot(p_c_vals,alpha_vals);
xlabel('p_c');
ylabel('\alpha');
title(['2-D HPD contour levels for p(\theta(', ...
        int2str(i_plot), '), \theta(', ...
        int2str(j_plot), ') | y )']);
% find p_c value for desired HPD
p_c_HPD = 0;
for n=1:N_bins
    if(alpha_vals(n) <= (1-alpha))
        p_c_HPD = p_c_vals(n);
        break;
    end
end
% compute HPD for this contour value
for ki=1:N_bins
    for kj=1:N_bins
        if(bin_2Dp(ki,kj) >= p_c_HPD)
            HPD_2D(ki,kj) = 1;
        end
    end
end

% make contoura plot of HPD
figure;
contour(bin_2Dic,bin_2Djc,HPD_2D,1);
xlabel(['\theta(', int2str(i_plot), ')']);
ylabel(['\theta(', int2str(j_plot), ')']);
title(['2-D HPD for p(\theta(', ...
        int2str(i_plot), '), \theta(', ...
        int2str(j_plot), ') | y )']);
figure;
contourf(bin_2Dic,bin_2Djc,HPD_2D,1);
xlabel(['\theta(', int2str(i_plot), ')']);
ylabel(['\theta(', int2str(j_plot), ')']);
title(['2-D HPD for p(\theta(', ...
        int2str(i_plot), '), \theta(', ...
        int2str(j_plot), ') | y )']);

return;