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)$.

Download this tutorial as a Jupyter notebook

Specification of the Log Posterior Density

When using the AdaptEllipticalSliceSampler.jl package, the main component that the user must specify is the log density of the target distribution. We will construct a Julia function to evaluate the log density (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
using MCMCChains

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)

Chains MCMC chain (10000×3×1 Array{Float64, 3}):

Iterations        = 1:1:10000
Number of chains  = 1
Samples per chain = 10000
Wall duration     = 4.07 seconds
Compute duration  = 4.07 seconds
parameters        = param_1, param_2
internals         = lp

Use `describe(chains)` for summary statistics and quantiles.

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). Using the Turing ecosystem we can easily view the results as follows:

# Discard the first 2500 samples to burn-in
results = results[2500:end, :, :]
# Note `describe()` requires MCMCChains
describe(results)

Chains MCMC chain (7501×3×1 Array{Float64, 3}):

Iterations        = 2500:1:10000
Number of chains  = 1
Samples per chain = 7501
Wall duration     = 4.07 seconds
Compute duration  = 4.07 seconds
parameters        = param_1, param_2
internals         = lp

Summary Statistics

  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

     param_1    0.2188    0.9036    0.0714   247.4827   158.7572    1.0231     ⋯
     param_2    0.5947    0.9022    0.0513   257.4737   171.6587    1.0236     ⋯

                                                                1 column omitted

Quantiles

  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

     param_1   -2.1318   -0.1824    0.5050    0.8964    1.1304
     param_2   -0.9310   -0.1439    0.5702    1.3027    2.2065

Note that in scenarios where sampling is quick, most of the Wall duration or Compute duration will be due to compilation costs (recall AGESS using the same call as before, and you should see a significant drop in computational costs). Additionally, we can view the MCMC samples by directly extracting the MCMC samples and plotting them ourselves.

# Note the 3rd dimension is used when you run multiple chains in parallel
# Here we only run one chain
scatter(results[:,1,1], results[:,2,1], legend = false)

Additionally, we can look at the plot of the log posterior density evaluated at each iteration of the Markov chain:

# Note this is results[:lp] can also be called by results[:,3,1]
plot(results[:lp], 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)

Chains MCMC chain (500000×3×1 Array{Float64, 3}):

Iterations        = 1:1:500000
Number of chains  = 1
Samples per chain = 500000
Wall duration     = 4.74 seconds
Compute duration  = 4.74 seconds
parameters        = param_1, param_2
internals         = lp

Use `describe(chains)` for summary statistics and quantiles.

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 assess the results as follows:

# Discard the first 125000 samples to burn-in
results_TB = results_TB[125001:500000,:,:]
describe(results_TB)

Chains MCMC chain (375000×3×1 Array{Float64, 3}):

Iterations        = 125001:1:500000
Number of chains  = 1
Samples per chain = 375000
Wall duration     = 4.74 seconds
Compute duration  = 4.74 seconds
parameters        = param_1, param_2
internals         = lp

Summary Statistics

  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ⋯
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64   ⋯

     param_1    1.1687    4.8624    0.0667   5347.0396   4975.3641    1.0003   ⋯
     param_2   -0.0840    2.5327    0.0339   6251.3444   9330.1457    1.0000   ⋯

                                                                1 column omitted

Quantiles

  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

     param_1   -5.8425   -4.5964    4.6891    4.9457    6.0798
     param_2   -3.1191   -2.1490   -1.5198    2.7675    3.9313

Similarly to before, we can view the MCMC samples. This time we will call them by the variable names, which are by default $\text{param_1},\dots,\text{param_P}$ unless specified through the param_names argument when calling AGESS.

# Note the 3rd dimension is used when you run multiple chains in parallel
# Here we only run one chain
scatter(results_TB[:param_1], results_TB[:param_2], legend = false)

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

plot(results_TB[:lp], 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].

# Free the sampled chains before the next tutorial's page is built.
results = results_TB = nothing
GC.gc()

This page was generated using Literate.jl.