function [lowerBound,upperBound] = hpd(draws,p,nBins)
%{

Given draws from a posterior find the p-percent highest posterior
density assuming the density is unimodal. Sort the draws into nBins

Dick Startz
October 2016
%}
% get counts, then cumulate highest until they add up to p

[Values,edges] = histcounts(draws,nBins);
highestBinIndex = find(Values == max(Values),1,'first');
Values = Values/length(draws);

lowerLimit = highestBinIndex;
upperLimit = highestBinIndex;
cumPercentage = max(Values);

% now look to the right and the left
while (cumPercentage <= p)
    if lowerLimit>1 && upperLimit<nBins
        if Values(lowerLimit-1) > Values(upperLimit+1)
            lowerLimit = lowerLimit - 1;
            cumPercentage = cumPercentage + Values(lowerLimit);
        else
            upperLimit = upperLimit + 1;
            cumPercentage = cumPercentage + Values(upperLimit);
        end
    elseif lowerLimit == 1 %can't go lower
        upperLimit = upperLimit + 1;
        cumPercentage = cumPercentage + Values(upperLimit);
    elseif upperLimit == nBins %can't go higher
        lowerLimit = lowerLimit - 1;
        cumPercentage = cumPercentage + Values(lowerLimit);
    else
        error('HPD is being stupid'); % this shpu;ld never happen
    end
end
lowerBound = edges(lowerLimit);
if upperLimit < nBins
    upperBound = edges(upperLimit+1);
else
    upperBound = edges(upperLimit);
end

end

