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);
0.3s

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
62.3s
]add Distributions LinearAlgebra StatsBase MCMCChains
1.0s
using NestedSamplers
using Distributions
using LinearAlgebra
using StatsBase: sample, Weights
using MCMCChains: Chains
52.8s

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
2.5s
prior_transform (generic function with 1 method)

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

model = NestedModel(logl, prior_transform)
0.6s
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())
6.9s
Nested(ndims=25, nactive=500, enlarge=1.25, update_interval=56250) bounds=Ellipsoid{Float64}(ndims=25) proposal=NestedSamplers.Proposals.Slice slices: Int64 5 scale: Float64 1.0 logz=-1.0e300 log_vol=-6.215607931755531 H=0.0

Next the chains object is constructed.

# by default, dlogz_convergence is used
chain = sample(model, spl; dlogz=0.01, chain_type=Chains)
47.4s

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())
0.6s
Nested(ndims=25, nactive=500, enlarge=1.25, update_interval=5000) bounds=Ellipsoid{Float64}(ndims=25) proposal=NestedSamplers.Proposals.RSlice slices: Int64 5 scale: Float64 1.0 logz=-1.0e300 log_vol=-6.215607931755531 H=0.0
chain = sample(model, spl; dlogz=0.01, chain_type=Chains)
3.4s

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
0.2s
prior_transform (generic function with 1 method)
model = NestedModel(logl, prior_transform)
0.2s
NestedModel(logl, prior_transform)
spl = Nested(ndims, 100, bounds=Bounds.Ellipsoid, proposal=Proposals.Uniform()) 
0.8s
Nested(ndims=3, nactive=100, enlarge=1.25, update_interval=150) bounds=Ellipsoid{Float64}(ndims=3) proposal=NestedSamplers.Proposals.Uniform logz=-1.0e300 log_vol=-4.610166019324897 H=0.0
# by default, dlogz_convergence is used
chain = sample(model, spl; dlogz=0.2, param_names = ["x", "y", "z"], chain_type=Chains)
2.2s

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
0.5s
logl2 (generic function with 1 method)
model = NestedModel(logl2, prior_transform)
0.2s
NestedModel(logl2, prior_transform)
spl = Nested(ndims, 100, bounds=Bounds.Ellipsoid, proposal=Proposals.Uniform()) 
0.2s
Nested(ndims=3, nactive=100, enlarge=1.25, update_interval=150) bounds=Ellipsoid{Float64}(ndims=3) proposal=NestedSamplers.Proposals.Uniform logz=-1.0e300 log_vol=-4.610166019324897 H=0.0
chain = sample(model, spl; dlogz=0.2, param_names = ["x", "y", "z"], chain_type=Chains)
1.5s

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 Society493(3), 3132-3158

Runtimes (1)