Custom explorers

Pigeons have several built-in explorer kernels such as AutoMALA and a SliceSampler. However when the state space is neither the reals nor the integers, or for performance reasons, it may be necessary to create custom exploration MCMC kernels.

Creating a new explorer

We show how to create a new explorer, for pedagogy, a simple independence Metropolis algorithm, applied to our familiar unidentifiable toy example, based on Julia black-box implementation. 

struct MyIndependenceSampler
    which_parameter_index::Int
end
function Pigeons.step!(explorer::MyIndependenceSampler, replica, shared)
    state = replica.state
    rng = replica.rng
    i = explorer.which_parameter_index
    # Note: the log_potential is an InterpolatedLogPotential between the target and reference
    log_potential = Pigeons.find_log_potential(replica, shared.tempering, shared)
    log_pr_before = log_potential(state)
    # propose
    state_before = state[i]
    state[i] = rand(rng)
    log_pr_after = log_potential(state)
    # accept-reject step
    accept_ratio = exp(log_pr_after - log_pr_before)
    if accept_ratio < 1 && rand(rng) > accept_ratio
        # reject: revert the move we just proposed
        state[i] = state_before
    end # (nothing to do if accept, we work in-place)
end

Creating combinations of explorers

To alternate between two explorers, use Compose: for example continuing on our example, we want to alternate between sampling the two parameters of our model:

pt = pigeons(
        target = MyLogPotential(100, 50),
        reference = MyLogPotential(0, 0),
        explorer = Compose(MyIndependenceSampler(1), MyIndependenceSampler(2))
    )
┌ Info: Neither traces, disk, nor online recorders included.
│    You may not have access to your samples (unless you are using a custom recorder, or maybe you just want log(Z)).
└    To add recorders, use e.g. pigeons(target = ..., record = [traces; record_default()])
┌ Warning: It looks like sample_iid!() is not implemented for a
│ reference_log_potential of type Main.MyLogPotential.
│ Instead, using step!().
└ @ Pigeons ~/work/Pigeons.jl/Pigeons.jl/src/targets/target.jl:51
────────────────────────────────────────────────────────────────────────────
  scans        Λ        time(s)    allc(B)  log(Z₁/Z₀)   min(α)     mean(α)
────────── ────────── ────────── ────────── ────────── ────────── ──────────
        2       2.49     0.0467   2.37e+06      -15.4    5.9e-08      0.724
        4       0.85   0.000107   1.56e+04      -5.56      0.574      0.906
        8       1.54   0.000159   2.77e+04      -4.68      0.237      0.829
       16        1.7   0.000281   5.79e+04      -4.82      0.648      0.811
       32       1.27   0.000479   7.94e+04      -5.42      0.669      0.858
       64       1.44   0.000801   8.12e+04      -4.77       0.73       0.84
      128        1.5    0.00137   7.79e+04      -5.14      0.756      0.833
      256       1.57    0.00255   7.79e+04      -5.08      0.782      0.826
      512       1.46     0.0049   7.79e+04      -4.94      0.801      0.838
 1.02e+03        1.5    0.00958   7.79e+04      -4.96      0.808      0.834
────────────────────────────────────────────────────────────────────────────

Adaptation

We assume the following model for MCMC explorer adaptation:

  1. during each PT round, statistics are collected distributively,
  2. at the end of each round, the statistics are reduced and shared, and the explorers are given an opportunity to update based on these statistics.

To control (1), use explorer_recorder_builders. For example, AutoMALA requests online statistics to be computed on uncontrainted parameters to perform pre-conditioning. By default, explorer_recorder_builders returns an empty list.

To control (2), use adapt_explorer which is fed the reduced statistics. By default, adapt_explorer is a no-op.