Baytes

Baytes.jl is a sampling library to perform Markov Chain and Sequential Monte Carlo proposal steps. It consists of several sub-libraries, such as BaytesMCMC.jl, BaytesFilters.jl, BaytesPMCMC.jl and BaytesSMC.jl, and provides an interface to combine kernels from these libraries.

In order to start, we have to define parameter and an objective function first. Let us use the model initially defined in the ModelWrappers.jl introduction:

using ModelWrappers, Baytes
using Distributions, Random, UnPack
_rng = Random.GLOBAL_RNG

#Create initial model and data
myparameter = (μ = Param(Normal(), 2.0), σ = Param(Gamma(), 3.0))
mymodel = ModelWrapper(myparameter)
data = randn(1000)
#Create objective for both μ and σ and define a target function for it
myobjective = Objective(mymodel, data, (:μ, :σ))
function (objective::Objective{<:ModelWrapper{BaseModel}})(θ::NamedTuple)
	@unpack data = objective
	lprior = Distributions.logpdf(Distributions.Normal(),θ.μ) + Distributions.logpdf(Distributions.Exponential(), θ.σ)
    llik = sum(Distributions.logpdf( Distributions.Normal(θ.μ, θ.σ), data[iter] ) for iter in eachindex(data))
	return lprior + llik
end

Sampling this model is straightforward. For instance, we can jointly estimate μ and σ via NUTS:

trace1, algorithm1 = sample(_rng, mymodel, data, NUTS((:μ, :σ)))

Trace.val contains samples stored as NamedTuple. Trace.diagnostics contains diagnostics that were returned when sampling the corresponding parameter. Trace.info contains useful summary information printing and summarizing the parameter estimates. Algorithm returns the kernel that was used and tuned during the sampling process. Note that a new kernel is initiated for each chain separately.

Alternatively, you might want to sample parameter with different kernels. This is usually useful when working with state space models.

trace2, algorithm2 = sample(_rng, mymodel, data, NUTS((:μ, )), Metropolis( (:σ,) ))

Similar to other Baytes packages, sampling arguments may be tweaked via a Default struct.

sampledefault = SampleDefault(;
    dataformat=Batch(),
    tempering=IterationTempering(Float64, UpdateFalse(), 1.0, 1000),
    chains=4,
    iterations=2000,
    burnin=max(1, Int64(floor(2000/10))),
    thinning = 1,
    safeoutput=false,
    printoutput=true,
    printdefault=PrintDefault(),
    report=ProgressReport(),
)
trace3, algorithm3 = sample(_rng, mymodel, data, NUTS((:μ, :σ)); default = sampledefault)

Note that hyperparameter that are specific to any kernel will have to be assigned in the MCMC constructor itself, i.e.:

mcmc4 = HMC((:μ,); proposal = ConfigProposal(; metric = MDense()), GradientBackend = :ReverseDiff,)
trace4, algorithm4 = sample(_rng, mymodel, data, mcmc4; default = sampledefault)

You can reuse returned algorithm container to sample again with the pre-tuned kernels. In this this case, on can call sample!. The information in trace will assure that sampling is continued with the correct specifications. After sample! is finished, a new trace is returned that contains sampling information.

iterations = 10^3
trace5, algorithm4 = sample!(iterations, _rng, mymodel, data, trace4, algorithm4)

Per default, sample and sample! return summary information of the chain and diagnostics using MCMCDiagnosticTools.jl, unless printdefault = false. For further inference, one can manually convert trace.val into a 3D-array that is compatible with MCMCChains.jl:

burnin = 0
thinning = 1
tagged = Tagged(mymodel, (:μ,:σ) )
tracetransform = TraceTransform(trace5, mymodel, tagged)
array_3dim = Baytes.trace_to_3DArray(trace5, tracetransform)

using MCMCChains
MCMCChains.Chains(array_3dim, trace5.info.sampling.paramnames)

This package is still highly experimental - suggestions and comments are always welcome!