function [low,up] = klaus_hpdi(x,a,B)
%
%
% This generates a highest posterior density interval for draws in vector
% "x".
% The posterior of x, i.e. p(x|y), must be unimodal.
% a is a number between 0 and 1, and indicates the 100(1-a)% desired "level of confidence"
% for the interval. Usually, a will be 0.05. 
% B is the number of bins you want to have for the underlying histogram
% that captures empirical frequencies. The higher this number, the more
% accurate, but slower, is the algorithm.

int=sort(x)';%sort x from smallest to largest
l=length(int); %total number of draws
bw=(max(int) - min(int))/(B-1);  %set binwidth s.t. we get B equally sized bins
% one might increase that number for higher accuracy
t= min(int):bw:max(int);
h=histc(int,t)/l;%produces empirical frequency for each bin

lt=length(t); %number of bins
startmat=[];
stopmat=[];


for i=1:lt
    int=0;
    j=i;
    while j<lt
        int=int+h(j); %add up probabilities until we hit (1-a)
        if int>(1-a)
            startx=min(x)+(i-1)*bw;
            stopx=min(x)+j*bw;
            startmat=[startmat startx]; 
            stopmat=[stopmat stopx];%j*bw indicates the lower bound of the jth bin,
            % this is a consevative choice for an upper bound
            break %exit while-loop
        end
        j=j+1;
    end
end


% now we need to find the smallest of these intervals
rangemat=stopmat-startmat;
all=[rangemat' startmat' stopmat'];
all=sortrows(all,1);

low=all(1,2); % lower x for HDPI
up=all(1,3); % upper x for HDPI






   

