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. You can find a good reference for the method in the paper "Computing Bayes Factors Using Thermodynamic Integration" Additionally I wrote a short blog post about it.
For a different way of computing the evidence integral see also my BayesianQuadrature package.
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(Diagonal(0.5 * ones(D))) # The prior distribution
likelihood = MvNormal(Diagonal(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.211990123364176The algorithm also works on multiple threads by calling :
alg = ThermInt(n_samples=5000)
logZ = alg(logprior, loglikelihood, rand(prior), TIThreads())or on multiple processes:
alg = ThermInt(n_samples=5000)
logZ = alg(logprior, loglikelihood, rand(prior), TIDistributed())Note that you need to load ThermodynamicIntegration and other necessary external packages on your additional processes via @everywhere.
Right now sampling is based on AdvancedHMC.jl, with the ForwardDiff.jl 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.
You can disactivate the progress by calling progress=false
Lartillot, N., & Philippe, H. (2006). Computing Bayes factors using thermodynamic integration