Banana Distribution

The Banana distribution^{[1] [2]} is a classical example of a non-convex and non-elliptical two-dimensional distribution. While it is a relatively low-dimensional target distribution, gradient-based samplers often have trouble drawing samples from the target distribution due to areas of high curvature^{[3] [4]}. The Banana distribution on $\boldsymbol{\theta} := (\theta_1, \theta_2)$ can be specified in the following hierarchical way:

\[Y_i \sim \mathcal{N}\left((\theta_1 - \mu_1) + (\theta_2 - \mu_2)^2, 1\right)\;\;\;\;\; \theta_1 \sim \mathcal{N}(0, 4)\;\;\;\;\; \theta_2 \sim \mathcal{N}(0.5, 4),\]

where $Y_i$ is generated from $\mathcal{N}(1, 1)$.

Specification of the Log Posterior Density

When using the AdaptEllipticalSliceSampler.jl package, the main component that the user must specify is the log pdf of the target distribution. We will construct a Julia function to evaluate the log pdf (up to a constant) of the posterior distribution of $\boldsymbol{\theta}$ based on the hierarchical representation above.

import Random
using AdaptEllipticalSliceSampler
using Distributions
using Plots
using LinearAlgebra

Random.seed!(123)

function log_posterior(theta::AbstractVector{<:AbstractFloat}, Y::Vector{<:AbstractFloat},
                       mu::Vector{<:AbstractFloat})
    mean = (theta[1] - mu[1]) + (theta[2] - mu[2])^2
    ## Normal distribution for likelihood
    lpdf = -0.5 * norm(Y .- mean)^2
    ## Prior distributions
    lpdf += logpdf(Normal(0, 2), theta[1]) + logpdf(Normal(0.5, 2), theta[2])

    return lpdf
end

log_posterior (generic function with 1 method)

We can see that the function log_posterior takes 3 arguments: theta (parameters of interest), Y (data), and mu the mean parameters.

Sampling via AGESS

To generate samples from the posterior distribution, we will use adaptive generalized elliptical slice sampling^[3] (AGESS). We will also specify the prior mean $\mu_0$ and prior scale $\Sigma_0$ parameters.

### Generate our data (100 observations)
Y = randn(100) .+ 1

### random mean component of our model
mu = randn(2) * 3

### Sample Via AGESS
### 10000 MCMC samples (N_MCMC = 10000)
### Target distribution dimension (P = 2)
μ_0 = [0.0, 0.5]
Σ_0 = diagm(4*ones(2))
results = AGESS(t -> log_posterior(t, Y, mu), 10000, 2; μ_0 = μ_0, Σ_0 = Σ_0)

AdaptEllipticalSliceSampler.MCMC_output([0.0 0.0; -0.05651690658244081 -0.3997655799791846; … ; 1.1539346083255393 0.3238943077018044; 0.9136356054699446 0.8545593822040143], [-76.16808908656144, -49.61565986487586, -47.433538761845845, -47.427276626757795, -48.150718079827584, -48.30909513344463, -46.400491505036484, -46.42910715521959, -46.551665449056344, -46.362137111488366  …  -47.33794337482202, -47.378936205988325, -48.66355428717292, -47.89853939557331, -46.48261193539375, -46.282653552090046, -46.21195868123721, -46.20307032202008, -47.385508720699406, -46.19116980710264], AGESS_MCMC_params(Main.var"#2#3"(), [0.0, 0.5], [4.0 0.0; 0.0 4.0], true, 6.0, 0.5, 0.05, 0.5, 0.05, 10000, 2), [1.2421769420006796 -0.5164845728052629; -0.5164845728052629 1.1242535241127352], [-0.13264061038800024, 0.7398460480695371])

We can see that we use an anonymous function to specify the likelihood, where theta (our parameters of interest) are the only parameters of interest. Here, Y and mu are the globally specified parameters (not-changing). We can visualize the results as follows:

## Discard the first 2500 samples to burn-in
scatter(results.samps[2501:10000,1], results.samps[2501:10000,2], alpha = 0.4,
        legend = false)

We can also view the log posterior density evaluated at each iteration of the Markov chain:

plot(results.l_pdf[1:10000], legend = false)

Twin Banana Distribution

The Twin Banana^[3] expands the banana distribution to a target distribution with two isolated bananas. Similarly, we can construct through a similar hierarchical representation:

\[Y_i \sim \mathcal{N}\left(0.1(\theta_1 - \mu_1)^2 - 0.5(\theta_2 - \mu_2)^4 - 10(\theta_1 - \mu_1)(\theta_2 - \mu_2), 100\right),\]

\[\theta_1 \sim \mathcal{N}(0, 4)\;\;\;\;\; \theta_2 \sim \mathcal{N}(0, 4),\]

where $Y_i$ is generated from $\mathcal{N}(100, 100)$.

Specification of the Log Posterior Density

We will construct a Julia function to evaluate the log pdf (up to a constant) of the posterior distribution of $\boldsymbol{\theta}$ based on the hierarchical representation of the Twin Banana distribution.

function TB_log_posterior(theta::AbstractVector{<:AbstractFloat}, Y::Vector{<:AbstractFloat},
                          mu::Vector{<:AbstractFloat})
    m = (0.1 * (theta[1] - mu[1])^2  -  0.5 * (theta[2] - mu[2])^4 -
            10 * (theta[1]- mu[1]) * (theta[2]- mu[2]))
    ## Normal distribution for likelihood
    lpdf = -(0.5 / 100) * norm(Y .- m)^2
    ## Prior distributions
    lpdf += logpdf(Normal(0, 2), theta[1]) + logpdf(Normal(0, 2), theta[2])

    return lpdf
end

TB_log_posterior (generic function with 1 method)

Again, we can see that the function TB_log_posterior takes 3 arguments: theta (parameters of interest), Y (data), and mu the mean parameters.

Sampling via AGESS

To generate samples from the posterior distribution, we will use adaptive generalized elliptical slice sampling^[3] (AGESS). We will specify the prior scale $\Sigma_0$ parameter, however since the prior mean is centered at 0, we do not have to explicitly specify $\mu_0$.

Random.seed!(123)
### Generate our data (100 observations)
Y_TB = randn(100) * 10 .+ 100

### random mean component of our model
mu_TB = randn(2) .* 3

### Sample Via AGESS
### 500000 MCMC samples (N_MCMC = 500000)
### Target distribution dimension (P = 2)
Σ_0 = diagm(ones(2) .* 4.0)
results_TB = AGESS(t -> TB_log_posterior(t, Y_TB, mu_TB), 500000, 2; Σ_0 = Σ_0)

AdaptEllipticalSliceSampler.MCMC_output([0.0 0.0; 0.04197061357566199 0.29796767033645155; … ; 4.660194314103457 -2.3458176740072947; 4.630823403531412 -2.329004171046704], [-5045.382519814066, -5007.468682915566, -4899.634110809638, -4885.355839703444, -4746.096383072505, -4580.3938817299495, -4558.393261095823, -4551.340452247577, -4326.186966148448, -4089.9122883119344  …  -51.50402254498055, -50.582396665504646, -49.90643028035699, -50.057582303918124, -49.90263985007265, -50.239512353055375, -49.84027577126728, -49.97438685016611, -49.47295840482789, -49.84565795492394], AGESS_MCMC_params(Main.var"#5#6"(), [0.0, 0.0], [4.0 0.0; 0.0 4.0], true, 6.0, 0.5, 0.05, 0.5, 0.05, 500000, 2), [23.368216848417898 -11.935752557741019; -11.935752557741019 6.675591374813882], [1.3573921857229818, -0.10617800585922516])

We can see that we use an anonymous function to specify the posterior pdf, where theta (our parameters of interest) are the only parameters of interest. Here, Y and mu are the globally specified parameters (not-changing). We can visualize the results as follows:

## Discard the first 125000 samples to burn-in
scatter(results_TB.samps[125001:500000,1], results_TB.samps[125001:500000,2], alpha = 0.4,
        legend = false)

Again, we can also view the log posterior density evaluated at each iteration of the Markov chain:

plot(results_TB.l_pdf[1:500000], legend = false)

Conclusion

Sampling from these challenging, yet low-dimensional distributions is fairly simple using adaptive generalized elliptical slice sampling^[3]. Using this package, sampling from these target distributions is as simple as efficiently constructing a function that evaluates the log posterior density. Comparisons between AGESS and other standard sampling schemes can be found in the main manuscript^[3].