function I=credibleregion(f, prob)
if nargin<2, prob=0.95; end

% sort the densities, largest first
sf = flipud(sort(f(:)));

% compute cumulative probabilities
cf = cumsum(sf);

% find where prob is first reached
i = find(cf >= prob*sum(f(:)), 1);
dthreshold = sf(i);

% indicate areas where density exceeds that threshold
I = (f >= dthreshold);