ThermodynamicIntegration
Thermodynamic Integration
Thermodynamic integration is a technique from physics to get an accurate estimate of the log evidence. By creating a schedule going from the prior to the posterior and estimating the log likelihood at each step one gets a stable ad robust estimate of the log evidence.
A simple example
A simple package to compute Thermodynamic Integration for computing the evidence in a Bayesian setting.
You need to provide the logprior and the loglikelihood as well as an initial sample:
using Distributions, ThermodynamicIntegration
D = 5
prior = MvNormal(0.5 * ones(D)) # The prior distribution
likelihood = MvNormal(2.0 * ones(D))
logprior(x) = logpdf(prior, x) # The log-prior function
loglikelihood(x) = logpdf(likelihood, x) # The log-likelihood function
alg = ThermInt(n_samples=5000) # We are going to sample 5000 samples at every step
logZ = alg(logprior, loglikelihood, rand(prior)) # Compute the log evidence
# -8.244829688529377
true_logZ = -0.5 * (logdet(cov(prior) + cov(likelihood)) + D * log(2π)) # we compare twith the true value
# -8.211990123364176You can also simply pass a Turing model :
using Turing
@model function gauss(y)
x ~ prior
y ~ MvNormal(x, cov(likelihood))
end
alg = ThermInt(n_samples=5000)
model = gauss(zeros(D))
turing_logZ = alg(model)
# # -8.211990123364176Parallel sampling
The algorithm also works on multiple threads by calling :
alg = ThermInt(n_samples=5000)
logZ = alg(logprior, loglikelihood, rand(prior), TIParallelThreads())Sampling methods
Right now sampling is based on AdvancedHMC.jl, with the ForwardDiff AD backend.
To change the backend to Zygote or ReverseDiff (recommended for variables with large dimensions you can do:
using Zygote # (or ReverseDiff)
ThermoDynamicIntegration.set_adbackend(:Zygote) # (or :ReverseDiff)More samplers will be available in the future.
Further options
You can disactivate the progress by calling progress=false
alg()