Need help with highest posterior density interval plot using Makie

I am trying to create the highest posterior density interval plot, but haven’t been able to create the one where the HDPI region is shaded and follows the density curve.

p_grid = range(0,stop=1,length=1000)
prior = fill(1,1000)
likelihood = pdf.(Binomial.(3,p_grid),3)
posterior = likelihood .* prior
posterior = posterior / sum(posterior)
samples = sample(p_grid, Weights(posterior),1000)
hd = hpdi(samples,alpha=0.5)

using Makie
x= samples[hd[1] .< samples .< hd[2]]
a = kde(x,npoints = 496)
lims1, lims2 = searchsortedfirst(x, -1), searchsortedlast(x, 1)
y = pdf(a,x)
fig = Figure()
ax = Axis(fig,xlabel="Proportion Water", ylabel="Density")
lines = Makie.density!(ax,samples, strokewidth=1, color=:transparent)
fillb = Makie. band!(ax,x[lims1:lims2], fill(0,length(x)),y; color = (:dodgerblue, 0.3))
fig[1,1] = ax

The kind of plots I am getting are

k = kde(randn(100))

f, ax, l = lines(k)
band!(k.x[500:1500], zeros(1001), k.density[500:1500])


Thank You :smiley: