Bayesian parametric histogram

histogram.jl

# Analytical dirichlet model with flat priors --- see BDA3 ch23.
###

using Distributions, Plots, StatsBase

"""
    bayes_hist(x)

Bayesian parametric histogram using the Freedman-Diaconis rule for
number of bins.
"""
function bayes_hist(x)
    # Note, our intervals are right inclusive, (a, b]
    a, b = round(minimum(x) - 0.1 * (maximum(x) - minimum(x)), RoundDown, digits=1),
        round(maximum(x) + 0.1 * (maximum(x) - minimum(x)), RoundUp, digits=1)

    width = round(2 * iqr(x) / (length(x) ^ (1/3)), digits=2)
    K = ceil(Int, (b - a) / width)
    bins = range(a, step=width, length=K+1)

    tbl = countmap(findfirst(>(n), bins) - 1 for n in x)
    posteriors = [get(tbl, k, 0) + 1 for k in 1:K]

    heights = posteriors ./ sum(posteriors)
    mid_points = (bins[1:end-1] + bins[2:end]) / 2

    bar(mid_points, heights, linewidth=0.5, legend=false)
end

x = rand(Laplace(-5, 4), 1_000)

p1 = bayes_hist(x)
p2 = histogram(x, normalize=:probability, legend=false)

plot(p1, p2)