How to use tempering with custom sampler

[itsdfish]

Hello,

Thank you for putting together this package. I have a use case in which tempring is needed, but the random walk algorithm does not scale very well. I am interested in using tempering with another sampler, such as differential evolution MCMC. The README provides an example of using a custom sampler, but it was not clear to me how to incorporate tempering. Would you be willing to provide guidance or a simple example please?

Thank you in advance.

[mvihola]

Thanks for your question!

The package has been developed so that the random-walk adaptations can be used modularly within a custom sampler (which involves also other types of updates). The tempering, in contrast, has not been currently implemented to accomodate custom samplers easily.

It is still possible, here is a small example:

using AdaptiveMCMC

# Define custom data type for you own sampler
struct NaiveSampler <: AdaptState
end

# Define a specialised arwm_step! function, which does the update 
# (and possible adaptation) on a single level, on log-target of the form "beta*log_p + log_pr".
# This is a simple scalar random-walk sampler with Unif[-1,1] proposals.
import AdaptiveMCMC: arwm_step!, PVals
function arwm_step!(r::RWMState, s::NaiveSampler,  p::PVals{FT},
    beta::FT, log_p::Function, log_pr::Function, k::Int, do_adapt::Bool) :: Int where {
    FT <: AbstractFloat }

    # Form the proposal
    r.y[] = r.x[] + rand()*2-1 

    # Calculate log-density values & accept probabillity:
    p_y = log_p(r.y); p_pr_y = log_pr(r.y)
    dp = beta*(p_y - p.x) + p_pr_y - p.pr_x
    rxy = exp(dp); alpha = (rxy<one(FT)) ? rxy : one(FT)
    if rand(r.rng) <= alpha
        # Accept: interchange current/proposed states & log-target values
        accept!(r)
        # The function needs also to record the accepted state log_p and log_pr values to p:
        p.x = p_y; p.pr_x = p_pr_y
        acc = 1
    else
        acc = 0
    end
    acc # This is mostly for reporting acceptance statistics...
end

# Try the algorithm with two levels, targetting standard normal:
L = 2
S = [NaiveSampler() for i=1:L]
log_p(x) = -.5*sum(x.^2)
out = adaptive_rwm([0.0], log_p, 100_000; algorithm=S, L=L)

[itsdfish]

Much appreciated. Thank you!