Hi,
Per the title of the post, trying to figure out how to shade a density plot interval using StatsPlots
I don’t really have any usable example code. A couple of days playing with this and none the wiser.
Trying to reproduce Statistical Rethinking Chapter 3 plots.
Tried the formula outlined here which uses calls to plot(...) but couldn’t get it to work for density(...). No clue how to approach this problem using StatsPlots.
The goal,
thx,
It will work the same, i.e. you can do (roughly):
density(something)
plot!(x, y, fillrange = zero(x))
but you will need the y values for the range of x you want to fill, so something like pdf.(distribution, x)
Very big thanks to you Nils. That set me on the right path.
FWIW I was able to replicate the Statistical Rethinking plots after considerable homework.
The big missing piece for me was that denstiy() plots were based on kernel density calculations. That was novel enough. Enter KernelDensity.jl… Then to make it all work using Plots.jl.
It’s no wonder the graphic appears in an introductory chapter of SR without explanation.
Here is the code, and the plot itself.
using Distributions
using KernelDensity
using StatsPlots
p_grid = collect(range(0,1,10^3))
prior = ones(10^3)
likelihood = pdf.(Binomial.(9, p_grid), 6)
posterior = likelihood .* prior
posterior = posterior ./ sum(posterior)
samples = sample( p_grid, Weights(posterior), 10^4)
ksamp = kde(samples)
x = collect(ksamp.x)
y = ksamp.density
xi = findall( i -> i<0.5, x)
d11 = density(samples, xlabel="proportion water (p)", ylabel="density", label="kde", xticks=0.0:0.25:1, legend=:topleft)
plot!(x[xi], y[xi], fillrange=zeros(1), fc=:blues, label="p<0.5")
xi = findall(i -> (0.5<i) & (i<0.75), x)
d12 = density(samples, xlabel="proportion water (p)", ylabel="density", label="kde", xticks=0.0:0.25:1, legend=:topleft)
plot!(x[xi], y[xi], fillrange=zeros(1), fc=:blues, label="0.5<p<0.75")
qa, qb = quantile(samp, (0.0, 0.8))
xi = findall(i -> (qa<=i) & (i<=qb), x)
d21 = density(samples, xlabel="proportion water (p)", ylabel="density", label="kde", xticks=0.0:0.25:1, legend=:topleft)
plot!(x[xi], y[xi], fillrange=zeros(1), fc=:blues, label="lower 80%")
qa, qb = quantile(samp, (0.1, 0.9))
xi = findall(i -> (qa<=i) & (i<=qb), x)
d22 = density(samples, xlabel="proportion water (p)", ylabel="density", label="kde", xticks=0.0:0.25:1, legend=:topleft)
plot!(x[xi], y[xi], fillrange=zeros(1), fc=:blues, label="middle 80%")
plot(d11, d12, d21, d22, layout=(2, 2), size=(800,700))
