Calculate Credible (Highest Posterior Density, HPD) Intervals in Julia using Distributions.jl

cred_interval.jl

using PyPlot
using Distributions

function credible_interval(D::UnivariateDistribution; c=0.95, nx=1000)
	# Discretize over the support
	r = support(D)
	lb,ub = r.lb,r.ub

	# Histogram approximation of area under pdf
	x = linspace(lb,ub,nx)
	f = pdf(D,x)
	a = f * (x[2]-x[1])

	# Sort all the histogram bins
	sp = sortperm(a)
	x = x[sp]
	a = a[sp]

	a_tot = 0.0          # running sum of area
	cred_x = (Float64)[] # inputs to pdf within credible interval
	
	i = length(a)
	while a_tot < c
		a_tot += a[i]
		push!(cred_x,x[i])
		i -= 1
	end

	return sort(cred_x)
end

## example
x = linspace(0,1,200)
B = Beta(30,30)
cred_x = credible_interval(B;c=0.90)

figure()
fill_between(cred_x,pdf(B,cred_x),0,facecolor=(1,0.5,0.5))
plot(x,pdf(B,x))
title("pdf with 90% region highlighted")

Here is the output: