NestedSamplers.jl

Implementations of single- and multi-ellipsoidal nested sampling algorithms in pure Julia. We implement the AbstractMCMC.jl interface, allowing straightforward sampling from a variety of statistical models.

This package was heavily influenced by nestle, dynesty, and NestedSampling.jl.

Citing

If you use this library, or a derivative of it, in your work, please consider citing it. This code is built off a multitude of academic works, which have been noted in the docstrings where appropriate. These references, along with references for the more general calculations, can all be found in CITATION.bib

Installation

To use the nested samplers first install this library

julia> ]add NestedSamplers

Background

For statistical background and a more in-depth introduction to nested sampling, I recommend the dynesty documentation. In short, nested sampling is a technique for simultaneously estimating the Bayesian evidence and the posterior distribution (according to Bayes' theorem) from nested iso-likelihood shells. These shells allow a quadrature estimate of the integral for the Bayesian evidence, which we can use for model selection, as well as the statistical weights for the underlying "live" points, which is where we get our posterior samples from!

Usage

The samplers are built using the AbstractMCMC.jl interface. To use it, we need to create a NestedModel.

using Random
using AbstractMCMC
AbstractMCMC.setprogress!(false)
Random.seed!(8452);

[ Info: progress logging is disabled globally

using Distributions
using LinearAlgebra
using NestedSamplers
using LogExpFunctions: logaddexp

# Gaussian mixture model
σ = 0.1
μ1 = ones(2)
μ2 = -ones(2)
inv_σ = diagm(0 => fill(1 / σ^2, 2))

function logl(x)
    dx1 = x .- μ1
    dx2 = x .- μ2
    f1 = -dx1' * (inv_σ * dx1) / 2
    f2 = -dx2' * (inv_σ * dx2) / 2
    return logaddexp(f1, f2)
end
priors = [
    Uniform(-5, 5),
    Uniform(-5, 5)
]
# or equivalently
prior_transform(X) = 10 .* X .- 5
# create the model
# or model = NestedModel(logl, prior_transform)
model = NestedModel(logl, priors);

now, we set up our sampling using StatsBase.

Important: the state of the sampler is returned in addition to the chain by sample.

using StatsBase: sample, Weights

# create our sampler
# 2 parameters, 1000 active points, multi-ellipsoid. See docstring
spl = Nested(2, 1000)
# by default, uses dlogz for convergence. Set the keyword args here
# currently Chains and Array are support chain_types
chain, state = sample(model, spl; dlogz=0.2, param_names=["x", "y"])
# optionally resample the chain using the weights
chain_res = sample(chain, Weights(vec(chain["weights"])), length(chain));

Chains MCMC chain (9342×3×1 Array{Float64, 3}):

Iterations        = 1:9342
Number of chains  = 1
Samples per chain = 9342
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.0854    0.9986    0.0103   9397.6626   9401.2976    1.0000   ⋯
           y   -0.0866    1.0036    0.0103   9397.5975   9279.0105    1.0000   ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

           x   -1.1730   -1.0062   -0.8499    0.9870    1.1582
           y   -1.1713   -1.0140   -0.8576    0.9892    1.1538

let's take a look at the resampled posteriors

using StatsPlots
density(chain_res)
# analytical posterior maxima
vline!([-1, 1], c=:black, ls=:dash, subplot=1)
vline!([-1, 1], c=:black, ls=:dash, subplot=2)

and compare our estimate of the Bayesian (log-)evidence to the analytical value

analytic_logz = log(4π * σ^2 / 100)
# within 2-sigma
@assert isapprox(analytic_logz, state.logz, atol=2state.logzerr)

Contributions and Support

Primary Author: Miles Lucas (@mileslucas)

Contributions are always welcome! In general, contributions should follow ColPrac. Take a look at the issues for ideas of open problems! To discuss ideas or plan contributions, open a discussion.