Illustrations of use of NestedSamplers.jl

Here are two illustrations of the use of NestedSamplers.jl v0.5.0 which are inspired by dynesty.

First the environment is setup. Setting a seed will ensure reproducibility of these examples.

using Random

Random.seed!(1225);

Next we need to add the latest released version (v0.5.0) of NestedSamplers package, and also add the Distributions , LinearAlgebra , StatsBase, and MCMCChains packages.

]add NestedSamplers#v0.5.0

]add Distributions LinearAlgebra StatsBase MCMCChains

using NestedSamplers

using Distributions

using LinearAlgebra

using StatsBase: sample, Weights

using MCMCChains: Chains

NestedSamplers v0.5.0 supports five proposal algorithms:

1. Proposals.Uniform() is suitable for low dimensions,

2. Proposals.RWalk() and Proposals.RStagger() are suitable for low to moderate dimensions,

3. Proposals.Slice() and Proposals.RSlice() are suitable for high dimensions.

25-D correlated normal likelihood function

This example demonstrates how Proposals.Slice() and Proposals.RSlice() perform for 25-D correlated multivariate normal distribution.

ndims = 25    # number of dimensions

C = Matrix(1.0I, ndims, ndims)    # covariance is set to be an identity Matrix

C[C.==0] .= 0.4    # off-diagonal terms set to be correlated

Cinv = inv(C)    # the precision Matrix

lnorm = -0.5 * (log(2π)*ndims + log(det(C)))    # ln(normalization) constant

# 25-D correlated multivariate normal log-likelihood function

function logl(x)

    return -0.5 * (x' * (Cinv * x)) + lnorm

end

# prior transform

function prior_transform(u)

    return @. 10.0 * u - 5.0   # unit cube samples `u` are transformed to a uniform prior in -5.0 to 5.0 and 0 elsewhere, for each variable

end

Create the model using the log-likelihood and prior transform functions as arguments for NestedModel.

model = NestedModel(logl, prior_transform)

For the 25 (ndims) parameters, create the sampler with 500 active points. Since this is a unimodal case, the Ellipsoid (single) bound is selected to initialize the sampler.

spl = Nested(ndims, 500, bounds=Bounds.Ellipsoid, proposal=Proposals.Slice())

Next the chains object is constructed.

# by default, dlogz_convergence is used

chain = sample(model, spl; dlogz=0.01, chain_type=Chains)

Another sample is drawn using the random slice sampling proposal algorithm. Random slice sampling differs from slice sampling in the sense that, the former is a standard random implementation where each slice is along a random direction based on the provided axes.

spl = Nested(ndims, 500, bounds=Bounds.Ellipsoid, proposal=Proposals.RSlice())

chain = sample(model, spl; dlogz=0.01, chain_type=Chains)

Noisy Likelihoods

In some estimation and inference problems, the likelihood inline_formula not implemented is either intractable or its computation is numerically infeasible. In such cases, a noisy estimateinline_formula not implemented can be obtained, which is a function of the actual likelihood inline_formula not implemented and hyper-parameters inline_formula not implemented .

Using a 3-D multivariate normal distribution, sampling from a noiseless and a noisy likelihood case is illustrated here. The Proposals.Uniform() algorithm is used for sampling in both the case.

The noiseless case is demonstrated below:

ndims = 3    # number of dimensions

C = Matrix(1.0I, ndims, ndims)    # covariance is set to be an identity Matrix

Cinv = inv(C)    # the precision Matrix

lnorm = -0.5 * (log(2π)*ndims + log(det(C)))    # ln(normalization) constant

# 3-D multivariate normal log-likelihood function

function logl(x)

    return -0.5 * (x' * (Cinv * x)) + lnorm

end

# prior transform which transforms the unit cube samples `u` to be uniform in each dimension from -10.0 to 10.0 and 0 elsewhere

function prior_transform(u)

    return @. 20.0 * u - 10.0

end

model = NestedModel(logl, prior_transform)

spl = Nested(ndims, 100, bounds=Bounds.Ellipsoid, proposal=Proposals.Uniform()) 

# by default, dlogz_convergence is used

chain = sample(model, spl; dlogz=0.2, param_names = ["x", "y", "z"], chain_type=Chains)

A noisy case of the above model is demonstrated below. Noise is introduced by setting a variable noise = 1.0. Hence, noisy estimate of the above model, which is inline_formula not implemented, is set as inline_formula not implemented .

noise = 1.0

# 3-D multivariate normal log-likelihood function

function logl2(x)

    xp = @. rand(Normal(x, noise), 1)

    loglikeli = -0.5 * (xp' * (Cinv * xp)) + lnorm

    scale = -0.5 * (noise^2)    # scaling

    bias_correction = scale * ndims    # bias correction term, introduced to ensure that the noisy likelihood is unbiased relative to the true likelihood

    return loglikeli - bias_correction

end

model = NestedModel(logl2, prior_transform)

spl = Nested(ndims, 100, bounds=Bounds.Ellipsoid, proposal=Proposals.Uniform()) 

chain = sample(model, spl; dlogz=0.2, param_names = ["x", "y", "z"], chain_type=Chains)

References:

1. https://github.com/joshspeagle/dynesty

2. Miles Lucas, Saranjeet Kaur, Hong Ge, & Cameron Pfiffer. (2020, July 22). TuringLang/NestedSamplers.jl: v0.5.0 (Version v0.5.0). Zenodo. http://doi.org/10.5281/zenodo.3955218

3. Speagle, J. S. (2020). dynesty: A dynamic nested sampling package for estimating Bayesian posteriors and evidences. Monthly Notices of the Royal Astronomical Society, 493(3), 3132-3158