# 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