Running Pathfinder on Turing.jl models

This tutorial demonstrates how Turing can be used with Pathfinder.

We'll demonstrate with a regression example.

using AdvancedHMC, Pathfinder, Random, Turing
Random.seed!(39)

@model function regress(x, y)
    α ~ Normal()
    β ~ Normal()
    σ ~ truncated(Normal(); lower=0)
    y ~ product_distribution(Normal.(α .+ β .* x, σ))
end
x = 0:0.1:10
y = @. 2x + 1.5 + randn() * 0.2

model = regress(collect(x), y)
n_chains = 8

8

For convenience, pathfinder and multipathfinder can take Turing models as inputs and produce MCMCChains.Chains objects as outputs.

result_single = pathfinder(model; ndraws=1_000)

Single-path Pathfinder result
  tries: 1
  draws: 1000
  fit iteration: 22 (total: 27)
  fit ELBO: -1.95 ± 0.03
  fit distribution: Distributions.MvNormal{Float64, Pathfinder.WoodburyPDMat{Float64, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}, Matrix{Float64}, Pathfinder.WoodburyPDFactorization{Float64, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}, LinearAlgebra.UpperTriangular{Float64, Matrix{Float64}}}}, Vector{Float64}}(
dim: 3
μ: [1.4759589631397645, 2.006026941373265, -1.5478356855806579]
Σ: [0.0018286217896122648 -0.0002601667183752173 -0.0005873937036806257; -0.0002601667183752181 5.100356929144418e-5 5.198418656871739e-5; -0.0005873937036806262 5.198418656871667e-5 0.0064369057163803384]
)

result_multi = multipathfinder(model, 1_000; nruns=n_chains)

Multi-path Pathfinder result
  runs: 8
  draws: 1000
  Pareto shape diagnostic: 0.79 (bad)

Here, the Pareto shape diagnostic indicates that it is likely safe to use these draws to compute posterior estimates.

When passed a DynamicPPL.Model, Pathfinder also gives access to the posterior draws in a familiar Chains object.

result_multi.draws_transformed

Chains MCMC chain (1000×3×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
parameters        = α, β, σ

Use `describe(chains)` for summary statistics and quantiles.

We can also use these posterior draws to initialize MCMC sampling.

init_params = collect.(eachrow(result_multi.draws_transformed.value[1:n_chains, :, 1]))

8-element Vector{Vector{Float64}}:
 [1.4156874588659516, 2.016094841223518, 0.24468823022975472]
 [1.4867712189092908, 2.009537160953487, 0.22457994999991854]
 [1.4691467301037398, 2.0111813084487418, 0.22204577368872686]
 [1.4225743419997652, 2.0150513003204957, 0.19291870250704835]
 [1.375596352156451, 2.0205385232497464, 0.23071888349800143]
 [1.436418141335313, 2.0118829758483967, 0.20678691584596665]
 [1.4770215057770222, 2.0025924142409632, 0.2236454231227252]
 [1.4483836340938132, 2.015921280636908, 0.21591143831903215]

chns = sample(model, Turing.NUTS(), MCMCThreads(), 1_000, n_chains; init_params, progress=false)

Chains MCMC chain (1000×15×8 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 8
Samples per chain = 1000
Wall duration     = 9.23 seconds
Compute duration  = 7.27 seconds
parameters        = α, β, σ
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size

Use `describe(chains)` for summary statistics and quantiles.

We can use Pathfinder's estimate of the metric and only perform enough warm-up to tune the step size.

inv_metric = result_multi.pathfinder_results[1].fit_distribution.Σ
metric = Pathfinder.RankUpdateEuclideanMetric(inv_metric)
kernel = HMCKernel(Trajectory{MultinomialTS}(Leapfrog(0.0), GeneralisedNoUTurn()))
adaptor = StepSizeAdaptor(0.8, 1.0)  # adapt only the step size
nuts = AdvancedHMC.HMCSampler(kernel, metric, adaptor)

n_adapts = 50
n_draws = 1_000
chns = sample(
    model,
    externalsampler(nuts),
    MCMCThreads(),
    n_draws + n_adapts,
    n_chains;
    n_adapts,
    init_params,
    progress=false,
)[n_adapts + 1:end, :, :]  # drop warm-up draws

Chains MCMC chain (1000×16×8 Array{Float64, 3}):

Iterations        = 51:1:1050
Number of chains  = 8
Samples per chain = 1000
Wall duration     = 5.76 seconds
Compute duration  = 5.34 seconds
parameters        = α, β, σ
internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt

Use `describe(chains)` for summary statistics and quantiles.

See Initializing HMC with Pathfinder for further examples.