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