[ANN] AdaptEllipticalSliceSampler.jl: Adaptive Generalized Elliptical Slice Sampling

Package Announcements
, , ,
1

I’d like to announce AdaptEllipticalSliceSampler.jl — a Julia implementation of an adaptive generalized elliptical slice sampler (AGESS) from a recent paper that we have been working on.

Elliptical slice sampling (ESS) is a common gradient-free MCMC algorithm that can adapt to localized features in the target distribution, while scaling efficiently as the dimension of the target distribution grows. However, that efficiency quickly degrades as the mismatch between the prior (used to generate the ellipses of ESS) and the target distribution grows—a common situation when using diffuse priors, in high-dimensional settings, or when the likelihood is highly informative.

AGESS directly addresses this by adapting the ellipse-defining distribution to the geometry of the target distribution. The resulting ergodic adaptive algorithm scales well with dimension, while retaining the gradient-free, mode traversing features of ESS. While AGESS can be used to sample from fairly general continuous target distributions, we have found the following settings particularly promising:

  1. Non-differentiable target distributions (ruling out the use of gradient-based methods such as HMC)
  2. Multimodal target distributions with strong inter-parameter dependencies
  3. Elliptically contoured target distributions

Usage
To use AGESS, simply write a function that evaluates your log posterior density, then call the AGESS function. Tutorials covering a wide range of scenarios (Banana distributions, regression settings, Bayesian neural networks, deep GP surrogate modeling, constrained inference, and factor models) are available in the documentation. A simple example of Bayesian linear regression can be found in the following code:

using AdaptEllipticalSliceSampler
using Distributions, LinearAlgebra, Random

function generate_data(N::T, D::T) where {T<:Integer}
    β = randn(D) * (2 * log(D))^(1.0 / 4)
    x = randn(N, D)
    y = zeros(Float64, N)
    for i in 1:N
        y[i] = randn() * 0.5 + dot(x[i,:], β)
    end

    return β, x, y
end

## Generate data with 1000 observations and 10 covariates
D = 10
β, X, y = generate_data(1000, D)

## Function evaluation log posterior density
function log_posterior(Param::AbstractVector{Y}, X::AbstractMatrix{Y},
                       y::AbstractVector{Y}) where {Y<:AbstractFloat}
    P = length(Param)
    @views lpdf = -0.5 * (1 / exp(Param[P])) * norm(X * Param[1:P-1] - y)^2 -
            (0.5 * length(y) * Param[P])
    @views lpdf += -0.5 * norm(Param[1:P-1])^2      # std normal prior on coefficients
    lpdf += -1 * Param[P] - (1 / exp(Param[P]))     # IG(1,1) prior on variance parameter 
    return lpdf
end

P = D + 1            # dimension of target distribution
n_MCMC = 10000       # number of MCMC iterations

## Run AGESS
results = AGESS(Param -> log_posterior(Param, X, y), n_MCMC, P)

Any feedback would be greatly appreciated!

Links:

3 Likes