% Bayes_1D_HPD_SR.m
% This routine returns desired (1-alpha)
% Highest Probability Density (HPD) regions
% given the results of the 1-D marginal
% posterior density calculator,
% Bayes_MCMC_1Dmarginal_SR.m
%
% Routine inputs:
% bin_c : the bin center positions of each marginal
%         density histogram
% bin_p : the bin probability values of each
%         marginal density histogram
% j_plot: the numbers of each parameter whose
%         marginal density has been calculated
%
% alpha : for 95% CI, use alpha = 0.05. The HPD is
%         returned that has (1-alpha) fraction of
%         total marginal density
% Routine outputs:
% HPD_lo : the minimum values in each HPD
% HPD_hi : the maximum values in each HPD
%
% K. J. Beers. MIT ChE. 12/13/2004.
% ver. 7/20/2005, increased # of contour values

function [HPD_lo,HPD_hi] = Bayes_1D_HPD_SR( ...
    bin_c, bin_p, j_plot, alpha);

% extract dimensioning parameters
N_bins = size(bin_c,2);
N_contours = 10*N_bins;
N_plots = size(bin_c,1);
bin_delta = zeros(N_plots,1);
for m=1:N_plots
    bin_delta(m) = bin_c(m,2) - bin_c(m,1);
end

% initialize the output values
HPD_lo = zeros(N_plots,1);
HPD_hi = zeros(N_plots,1);

% for each parameter to be considered
for m = 1:N_plots
    j = j_plot(m);  % get # of parameter
    % set vector of trial probability contour values
    p_c_max = max(bin_p(m,:));
    p_c_vals = linspace(0,p_c_max,N_contours);
    alpha_vals = zeros(size(p_c_vals));
    % for each contour value, compute total density
    % within corresponding HPD
    for n=1:N_contours
        p_c = p_c_vals(n);
        k_above = find(bin_p(m,:) >= p_c);
        alpha_vals(n) = bin_delta(m)*sum(bin_p(m,k_above));
    end
    % make plot of alpha vs. p_c
    figure;
    plot(p_c_vals,alpha_vals);
    xlabel('p_c');
    ylabel('\alpha');
    title(['HPD contour levels for \theta(', ...
            int2str(j), ')']);
    % find p_c value for desired HPD
    for n=1:N_contours
        if(alpha_vals(n) <= (1-alpha))
            p_c_HPD = p_c_vals(n);
            break;
        end
    end
    % compute HPD for this contour value
    k_above = find(bin_p(m,:) >= p_c_HPD);
    % make plot of HPD
    plot_HPD = zeros(1,N_bins);
    plot_HPD(k_above) = 1;
    figure;
    plot(bin_c(m,:),plot_HPD);
    xlabel(['\theta(', int2str(j), ')']);
    ylabel('HPD region');
    title(['1-D HPD for \theta(', int2str(j), ...
            '), \alpha = ', num2str(alpha), ...
            ', p_c = ', num2str(p_c_HPD)]); 
    axis([min(bin_c(m,:)) max(bin_c(m,:))  ...
        -0.1   1.1]);
    % find lower and upper bounds of HPD
    HPD_lo(m) = min(bin_c(m,k_above));
    HPD_hi(m) = max(bin_c(m,k_above));
    
end



return;