# Shading interval in density plot

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))
``````

2 Likes