Description

A synthetic target introduced in Neal (2003) to benchmark algorithms on situations where the local curvature of the target density varies from one part of the space to another. Specifically, the shape of the target (see pair plots below) is such that when $y < 0$, the posterior is very narrow while for $y > 0$ it is wide.

The funnel is formally defined as follows:

\[\begin{aligned} y &\sim \text{Normal}(0, 3) \\ x_i &\sim \text{Normal}(0, \exp(y / \text{scale})) \;\text{ for }i \in \{1, \dots, d\} \\ \end{aligned}\]

Here we use the values:

While the example is artificial, it is useful since it combines certain features present in many real challenging targets (varying local curvature), while having known moments for the difficult dimension to explore, $y$ (since $y$ is marginally a normal distribution with known moments).

Stan implementation

Pair plot

Diagonal entries show estimates of the marginal densities as well as the (0.16, 0.5, 0.84) quantiles (dotted lines). Off-diagonal entries show estimates of the pairwise densities.

Movie linked below (🍿) superimposes 100 iterations of MCMC.

Trace plots

Moments

Cumulative traces

For each iteration $i$, shows the running average up to $i$, $\frac{1}{i} \sum_{n = 1}^{i} x_n$.

Local communication barrier

When the global communication barrier is large, many chains may be required to obtain tempered restarts.

The local communication barrier can be used to visualize the cause of a high global communication barrier. For example, if there is a sharp peak close to a reference constructed from the prior, it may be useful to switch to a variational approximation.

GCB estimation progress

Estimate of the Global Communication Barrier (GCB) as a function of the adaptation round.

The global communication barrier can be used to set the number of chains. The theoretical framework of Syed et al., 2021 yields that under simplifying assumptions, it is optimal to set the number of chains (the argument n_chains in pigeons()) to roughly 2Λ.

Last round estimate: $1.48023218918939$

Evidence estimation progress

Estimate of the log normalization (computed using the stepping stone estimator) as a function of the adaptation round.

Last round estimate: $0.022091653206380624$

Round trips

Number of tempered restarts as a function of the adaptation round.

A tempered restart happens when a sample from the reference percolates to the target. When the reference supports iid sampling, tempered restarts can enable large jumps in the state space.

Swaps plot

Pigeons summary

Pigeons inputs

Reproducibility

run(`git clone https://github.com/Julia-Tempering/InferenceReport.jl`)
cd("InferenceReport.jl")
run(`git checkout f26f12430f67a2c8c64988c8c8c1d80696a7d7fd`)

using Pkg 
Pkg.activate(".")
Pkg.instantiate()
 

using Pigeons
inputs = Inputs(; target = Pigeons.stan_funnel(2), variational = GaussianReference(first_tuning_round = 5), n_rounds = 10, record = [traces; round_trip; record_default()])
 

pt = pigeons(inputs)