Gaussian Shells
This example will explore the classic Gaussian shells model using Models.GaussianShells
.
Setup
For this example, you'll need to add the following packages
julia>]add Distributions MCMCChains Measurements NestedSamplers StatsBase StatsPlots
Define model
using NestedSamplers
model, logz = Models.GaussianShells()
let's take a look at a couple of parameters to see what the likelihood surface looks like
using StatsPlots
x = range(-6, 6, length=1000)
y = range(-2.5, 2.5, length=1000)
loglike = model.prior_transform_and_loglikelihood.loglikelihood
logf = [loglike([xi, yi]) for yi in y, xi in x]
heatmap(
x, y, exp.(logf),
xlims=extrema(x),
ylims=extrema(y),
xlabel="x",
ylabel="y",
)
Sample
using MCMCChains
using StatsBase
# using multi-ellipsoid for bounds
# using default rejection sampler for proposals
sampler = Nested(2, 1000)
chain, state = sample(model, sampler; dlogz=0.05, param_names=["x", "y"])
# resample chain using statistical weights
chain_resampled = sample(chain, Weights(vec(chain[:weights])), length(chain));
Results
chain_resampled
Chains MCMC chain (7154×3×1 Array{Float64, 3}):
Iterations = 1:7154
Number of chains = 1
Samples per chain = 7154
parameters = x, y
internals = weights
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
x -0.0245 3.7464 0.0450 7002.6038 6928.2097 1.0000 ⋯
y 0.0366 1.4344 0.0172 6908.8202 6791.4721 0.9999 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
x -5.5045 -3.4476 -1.3050 3.4685 5.4566
y -2.0577 -1.4026 0.0519 1.4635 2.0473
marginalkde(chain[:x], chain[:y])
plot!(xlims=(-6, 6), ylims=(-2.5, 2.5), sp=2)
plot!(xlims=(-6, 6), sp=1)
plot!(ylims=(-2.5, 2.5), sp=3)
density(chain_resampled)
vline!([-5.5, -1.5, 1.5, 5.5], c=:black, ls=:dash, sp=1)
vline!([-2, 2], c=:black, ls=:dash, sp=2)
using Measurements
logz_est = state.logz ± state.logzerr
diff = logz_est - logz
println("logz: $logz")
println("estimate: $logz_est")
println("diff: $diff")
logz: -1.75
estimate: -1.75 ± 0.051
diff: 0.00045 ± 0.051