% Bayes_1D_HPD_MR.m % % function [HPD_lo,HPD_hi] = Bayes_1D_HPD_MR( ... % bin_1Dc, bin_1Dp, j_plot_1D, alpha); % % 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_MR.m % % Routine inputs: % bin_1Dc : the bin center positions of each marginal % density histogram % bin_1Dp : the bin probability values of each % marginal density histogram % j_plot_1D: 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/15/2004. ver. 12/15/2004 function [HPD_lo,HPD_hi] = Bayes_1D_HPD_MR( ... bin_1Dc, bin_1Dp, j_plot_1D, alpha); % extract dimensioning parameters N_bins = size(bin_1Dc,2); N_plots = size(bin_1Dc,1); bin_delta = zeros(N_plots,1); for m=1:N_plots bin_delta(m) = bin_1Dc(m,2) - bin_1Dc(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_1D(m); % get # of parameter % set vector of trial probability contour values p_c_max = max(bin_1Dp(m,:)); 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); k_above = find(bin_1Dp(m,:) >= p_c); alpha_vals(n) = bin_delta(m)*sum(bin_1Dp(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_bins 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_1Dp(m,:) >= p_c_HPD); % make plot of HPD plot_HPD = zeros(1,N_bins); plot_HPD(k_above) = 1; figure; plot(bin_1Dc(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_1Dc(m,:)) max(bin_1Dc(m,:)) ... -0.1 1.1]); % find lower and upper bounds of HPD HPD_lo(m) = min(bin_1Dc(m,k_above)); HPD_hi(m) = max(bin_1Dc(m,k_above)); end return;