ThermodynamicIntegration.jl

Thermodynamic Integration for Turing models and more
Author theogf
Popularity
10 Stars
Updated Last
1 Year Ago
Started In
January 2021

ThermodynamicIntegration

Stable Dev Build Status Coverage Status Code Style: Blue ColPrac: Contributor's Guide on Collaborative Practices for Community Packages

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. 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 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.211990123364176

You 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.211990123364176

Parallel sampling

The 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.

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

Reference

Lartillot, N., & Philippe, H. (2006). Computing Bayes factors using thermodynamic integration