Initializing HMC with Pathfinder

The MCMC warm-up phase

When using MCMC to draw samples from some target distribution, there is often a lengthy warm-up phase with 2 phases:

converge to the typical set (the region of the target distribution where the bulk of the probability mass is located)
adapt any tunable parameters of the MCMC sampler (optional)

While (1) often happens fairly quickly, (2) usually requires a lengthy exploration of the typical set to iteratively adapt parameters suitable for further exploration. An example is the widely used windowed adaptation scheme of Hamiltonian Monte Carlo (HMC) in Stan [3], where a step size and positive definite metric (aka mass matrix) are adapted. For posteriors with complex geometry, the adaptation phase can require many evaluations of the gradient of the log density function of the target distribution.

Pathfinder can be used to initialize MCMC, and in particular HMC, in 3 ways:

Initialize MCMC from one of Pathfinder's draws (replace phase 1 of the warm-up).
Initialize an HMC metric adaptation from the inverse of the covariance of the multivariate normal approximation (replace the first window of phase 2 of the warm-up).
Use the inverse of the covariance as the metric without further adaptation (replace phase 2 of the warm-up).

This tutorial demonstrates all three approaches with DynamicHMC.jl and AdvancedHMC.jl. Both of these packages have standalone implementations of adaptive HMC (aka NUTS) and can be used independently of any probabilistic programming language (PPL). Both the Turing and Soss PPLs have some DynamicHMC integration, while Turing also integrates with AdvancedHMC.

For demonstration purposes, we'll use the following dummy data:

using LinearAlgebra, Pathfinder, Random, StatsFuns, StatsPlots

Random.seed!(91)

x = 0:0.01:1
y = @. sin(10x) + randn() * 0.2 + x

scatter(x, y; xlabel="x", ylabel="y", legend=false, msw=0, ms=2)

We'll fit this using a simple polynomial regression model:

\[\begin{aligned} \sigma &\sim \text{Half-Normal}(\mu=0, \sigma=1)\\ \alpha, \beta_j &\sim \mathrm{Normal}(\mu=0, \sigma=1)\\ \hat{y}_i &= \alpha + \sum_{j=1}^J x_i^j \beta_j\\ y_i &\sim \mathrm{Normal}(\mu=\hat{y}_i, \sigma=\sigma) \end{aligned}\]

We just need to implement the log-density function of the resulting posterior.

struct RegressionProblem{X,Z,Y}
    x::X
    J::Int
    z::Z
    y::Y
end
function RegressionProblem(x, J, y)
    z = x .* (1:J)'
    return RegressionProblem(x, J, z, y)
end

function (prob::RegressionProblem)(θ)
    (; σ, α, β) = θ
    (; z, y) = prob
    lp = normlogpdf(σ) + logtwo
    lp += normlogpdf(α)
    lp += sum(normlogpdf, β)
    y_hat = muladd(z, β, α)
    lp += sum(eachindex(y_hat, y)) do i
        return normlogpdf(y_hat[i], σ, y[i])
    end
    return lp
end

J = 3
dim = J + 2
model = RegressionProblem(x, J, y)
ndraws = 1_000;

DynamicHMC.jl

To use DynamicHMC, we first need to transform our model to an unconstrained space using TransformVariables.jl and wrap it in a type that implements the LogDensityProblems interface (see DynamicHMC's worked example):

using DynamicHMC, ForwardDiff, LogDensityProblems, LogDensityProblemsAD, TransformVariables
using TransformedLogDensities: TransformedLogDensity

transform = as((σ=asℝ₊, α=asℝ, β=as(Array, J)))
P = TransformedLogDensity(transform, model)
∇P = ADgradient(:ForwardDiff, P)

ForwardDiff AD wrapper for TransformedLogDensity of dimension 5, w/ chunk size 5

Pathfinder can take any object that implements this interface.

result_pf = pathfinder(∇P)

Single-path Pathfinder result
  tries: 1
  draws: 5
  fit iteration: 15 (total: 17)
  fit ELBO: -117.64 ± 0.23
  fit distribution: 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: 5
μ: [-0.35193742737688777, 0.24690060220372106, 0.058275882395669526, 0.11655175981944875, 0.1748276485411877]
Σ: [0.004969588924912717 -0.00022940256988482796 … 0.0005296595888833644 -0.00023978958346961997; -0.0002294025698848301 0.0258505173524506 … -0.017788412328272638 -0.00020020475341149188; … ; 0.0005296595888833648 -0.017788412328272833 … 0.743957303003753 -0.43851071755252763; -0.00023978958346961606 -0.00020020475341137067 … -0.43851071755252674 0.2653446953844584]
)

init_params = result_pf.draws[:, 1]

5-element Vector{Float64}:
 -0.2447069082536234
  0.29076814691738456
  0.501652842232575
 -1.6245650931617113
  1.1759862272969104

inv_metric = result_pf.fit_distribution.Σ

5×5 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}}}}:
  0.00496959   -0.000229403  -4.15468e-5   0.00052966  -0.00023979
 -0.000229403   0.0258505    -0.00138677  -0.0177884   -0.000200205
 -4.15468e-5   -0.00138677    0.0389324   -0.145779     0.0852796
  0.00052966   -0.0177884    -0.145779     0.743957    -0.438511
 -0.00023979   -0.000200205   0.0852796   -0.438511     0.265345

Initializing from Pathfinder's draws

Here we just need to pass one of the draws as the initial point q:

result_dhmc1 = mcmc_with_warmup(
    Random.default_rng(),
    ∇P,
    ndraws;
    initialization=(; q=init_params),
    reporter=NoProgressReport(),
)

(posterior_matrix = [-0.38196010598063457 -0.339200542658622 … -0.21286788208096405 -0.36488500531819285; 0.36891554897200945 0.22414671843018813 … 0.008252884174824759 0.19846084203414682; … ; -0.47010869387580156 0.7562176553747938 … 0.803444035127791 1.277081548859845; 0.7315556542188062 -0.08074612793142652 … 0.1805208793567661 -0.214174308455518], tree_statistics = DynamicHMC.TreeStatisticsNUTS[DynamicHMC.TreeStatisticsNUTS(-118.03684869028952, 6, turning at positions -29:34, 0.9954936974481424, 63, DynamicHMC.Directions(0xf1d2b2a2)), DynamicHMC.TreeStatisticsNUTS(-113.83378691934725, 6, turning at positions -53:-116, 0.9925274901796707, 127, DynamicHMC.Directions(0x506cb68b)), DynamicHMC.TreeStatisticsNUTS(-113.74268928091125, 6, turning at positions -60:3, 0.9920737081431408, 63, DynamicHMC.Directions(0xea352e03)), DynamicHMC.TreeStatisticsNUTS(-116.31959323457387, 6, turning at positions -60:3, 0.8917265229493909, 63, DynamicHMC.Directions(0xc0020b03)), DynamicHMC.TreeStatisticsNUTS(-116.9260123573649, 6, turning at positions -47:16, 0.9995981463025626, 63, DynamicHMC.Directions(0x233a4cd0)), DynamicHMC.TreeStatisticsNUTS(-117.78301511134116, 6, turning at positions -30:33, 0.8380558504632759, 63, DynamicHMC.Directions(0xe857fe21)), DynamicHMC.TreeStatisticsNUTS(-119.51643474768481, 5, turning at positions 18:21, 0.8918111391611169, 39, DynamicHMC.Directions(0x252143ad)), DynamicHMC.TreeStatisticsNUTS(-113.58379360601693, 6, turning at positions -53:10, 0.9902028864996022, 63, DynamicHMC.Directions(0x02050a0a)), DynamicHMC.TreeStatisticsNUTS(-115.92155858915203, 5, turning at positions -26:-57, 0.9937775147348274, 63, DynamicHMC.Directions(0xe7456046)), DynamicHMC.TreeStatisticsNUTS(-118.61953370413833, 7, turning at positions -124:3, 0.9712925390304312, 127, DynamicHMC.Directions(0x0aaaee83))  …  DynamicHMC.TreeStatisticsNUTS(-114.43573210651813, 7, turning at positions -65:62, 0.9777455464011678, 127, DynamicHMC.Directions(0xbc10d23e)), DynamicHMC.TreeStatisticsNUTS(-116.10457800153098, 6, turning at positions -30:-93, 0.9891880230727232, 127, DynamicHMC.Directions(0x9d553d22)), DynamicHMC.TreeStatisticsNUTS(-115.40116405975122, 6, turning at positions -39:24, 0.9885140324287273, 63, DynamicHMC.Directions(0xc7a5ec58)), DynamicHMC.TreeStatisticsNUTS(-114.94950335072556, 6, turning at positions -37:26, 0.9816177987155759, 63, DynamicHMC.Directions(0xe9c32f1a)), DynamicHMC.TreeStatisticsNUTS(-114.33297366694698, 6, turning at positions -28:35, 0.9395855602904865, 63, DynamicHMC.Directions(0x2fb8f923)), DynamicHMC.TreeStatisticsNUTS(-116.87743631729751, 7, turning at positions -19:108, 0.9540734839187606, 127, DynamicHMC.Directions(0x2646376c)), DynamicHMC.TreeStatisticsNUTS(-116.43132032990374, 6, turning at positions -55:8, 0.9999435092629735, 63, DynamicHMC.Directions(0x6c15e108)), DynamicHMC.TreeStatisticsNUTS(-116.61343738735164, 7, turning at positions -42:85, 0.8974009419450286, 127, DynamicHMC.Directions(0xa723b1d5)), DynamicHMC.TreeStatisticsNUTS(-118.34698972585295, 6, turning at positions -37:26, 0.9969469439902109, 63, DynamicHMC.Directions(0x0174891a)), DynamicHMC.TreeStatisticsNUTS(-119.50961446962357, 6, turning at positions -24:-87, 0.8928413548631957, 127, DynamicHMC.Directions(0x79e08aa8))], logdensities = [-113.239214714429, -112.60085172136098, -112.3160543494414, -115.409311001143, -113.59759176072411, -115.97961715947396, -113.1424202042765, -112.80155714263111, -113.97261109470224, -113.97652445727975  …  -113.76896461223564, -113.81365264254624, -113.58187318645764, -112.85835620131422, -113.25155787952356, -115.05762214859536, -113.78118185600012, -115.51528980533614, -116.95289072273472, -114.0796130107218], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.07425342098736351, 0.13670310737189154, 1.0212628693529269, 0.8723270822471226, 0.6335544476337074], ϵ = 0.046014685450374515)

Initializing metric adaptation from Pathfinder's estimate

To start with Pathfinder's inverse metric estimate, we just need to initialize a DynamicHMC.GaussianKineticEnergy object with it as input:

result_dhmc2 = mcmc_with_warmup(
    Random.default_rng(),
    ∇P,
    ndraws;
    initialization=(; q=init_params, κ=GaussianKineticEnergy(inv_metric)),
    warmup_stages=default_warmup_stages(; M=Symmetric),
    reporter=NoProgressReport(),
)

(posterior_matrix = [-0.33742951807169014 -0.3590462978442824 … -0.41277431452716273 -0.3815558298390257; 0.19993685254148025 0.35658596411279286 … 0.20120033694498837 0.10556629959797192; … ; -1.0034307241976048 0.07727641271274083 … 0.4055758497848926 0.3842362364112818; 0.9261601644555799 0.4169213575249797 … -0.006496473240251294 0.5593803918473792], tree_statistics = DynamicHMC.TreeStatisticsNUTS[DynamicHMC.TreeStatisticsNUTS(-114.65666851417791, 4, turning at positions -11:4, 0.9990987834917432, 15, DynamicHMC.Directions(0x1d811054)), DynamicHMC.TreeStatisticsNUTS(-113.74592323811973, 3, turning at positions -7:-14, 0.994871622929985, 15, DynamicHMC.Directions(0x3aac6a51)), DynamicHMC.TreeStatisticsNUTS(-115.81660343079031, 3, turning at positions 0:7, 0.9184541593551013, 7, DynamicHMC.Directions(0xbbda6e1f)), DynamicHMC.TreeStatisticsNUTS(-114.79236425169404, 2, turning at positions 0:3, 0.9231718352801566, 3, DynamicHMC.Directions(0x7f5cf633)), DynamicHMC.TreeStatisticsNUTS(-115.17923846164062, 2, turning at positions -1:2, 0.9846418165222502, 3, DynamicHMC.Directions(0x6889b3ae)), DynamicHMC.TreeStatisticsNUTS(-117.16559963981172, 3, turning at positions -1:6, 0.9985811409829103, 7, DynamicHMC.Directions(0x7deb39c6)), DynamicHMC.TreeStatisticsNUTS(-118.16690972343832, 3, turning at positions 0:7, 0.961751436745328, 7, DynamicHMC.Directions(0x847f6b6f)), DynamicHMC.TreeStatisticsNUTS(-116.850327884055, 3, turning at positions -5:2, 0.990610553048516, 7, DynamicHMC.Directions(0x59cf457a)), DynamicHMC.TreeStatisticsNUTS(-117.55224989465096, 3, turning at positions -7:0, 0.9446816676764309, 7, DynamicHMC.Directions(0x5fcfaf10)), DynamicHMC.TreeStatisticsNUTS(-119.23294500440144, 3, turning at positions -8:-15, 0.9444263030086547, 15, DynamicHMC.Directions(0x08c63cf0))  …  DynamicHMC.TreeStatisticsNUTS(-117.64078844584708, 3, turning at positions -7:0, 0.584291515069383, 7, DynamicHMC.Directions(0x43780d98)), DynamicHMC.TreeStatisticsNUTS(-116.83581403746392, 3, turning at positions -3:4, 0.9986957790493988, 7, DynamicHMC.Directions(0xab513b64)), DynamicHMC.TreeStatisticsNUTS(-113.86582386706856, 3, turning at positions -4:3, 0.9814082609662634, 7, DynamicHMC.Directions(0x5dcd8c4b)), DynamicHMC.TreeStatisticsNUTS(-115.55884693680161, 3, turning at positions -6:1, 0.937216593606017, 7, DynamicHMC.Directions(0x8135dde1)), DynamicHMC.TreeStatisticsNUTS(-114.8474587230229, 3, turning at positions -4:3, 0.9336431099298913, 7, DynamicHMC.Directions(0x14f5bc8b)), DynamicHMC.TreeStatisticsNUTS(-114.68481131011848, 3, turning at positions -7:0, 0.9925439906517933, 7, DynamicHMC.Directions(0xc44e5ca8)), DynamicHMC.TreeStatisticsNUTS(-115.6222608741956, 3, turning at positions -5:2, 0.9817273539216548, 7, DynamicHMC.Directions(0xd64ffc7a)), DynamicHMC.TreeStatisticsNUTS(-116.95132398257807, 3, turning at positions -7:0, 0.9994343568224784, 7, DynamicHMC.Directions(0xb45fcf30)), DynamicHMC.TreeStatisticsNUTS(-115.8645320716711, 3, turning at positions -3:4, 0.9972360669579999, 7, DynamicHMC.Directions(0x6db6e324)), DynamicHMC.TreeStatisticsNUTS(-114.99257016669672, 3, turning at positions -7:0, 0.9568476492000234, 7, DynamicHMC.Directions(0x04029ba8))], logdensities = [-113.05794116026404, -112.78759416920896, -112.50476567105275, -114.44863359489659, -115.05097564400037, -114.39050627522695, -115.07798332168154, -114.38999717421473, -116.51572150607979, -115.45316843985375  …  -116.72021513300471, -113.2002315391702, -112.5466069716753, -113.26454150143202, -113.54503309982547, -113.98554563782218, -115.40447655611962, -115.37241744151707, -112.62186148261598, -113.91196020441166], κ = Gaussian kinetic energy (Symmetric), √diag(M⁻¹): [0.07490535365229321, 0.13406333870827467, 0.888781775826629, 0.7718521461558405, 0.5499358853995192], ϵ = 0.44267467010095046)

We also specified that DynamicHMC should tune a dense Symmetric matrix. However, we could also have requested a Diagonal metric.

Use Pathfinder's metric estimate for sampling

To turn off metric adaptation entirely and use Pathfinder's estimate, we just set the number and size of the metric adaptation windows to 0.

result_dhmc3 = mcmc_with_warmup(
    Random.default_rng(),
    ∇P,
    ndraws;
    initialization=(; q=init_params, κ=GaussianKineticEnergy(inv_metric)),
    warmup_stages=default_warmup_stages(; middle_steps=0, doubling_stages=0),
    reporter=NoProgressReport(),
)

(posterior_matrix = [-0.28219054045100517 -0.38197783026878945 … -0.32046300525273297 -0.27115234824040313; 0.16018468970016705 0.3138886160351999 … 0.34767685580637087 0.19194536938558346; … ; -0.5755585612465779 0.3132810846872743 … -0.18106669739492243 0.5500048723602323; 0.44375726239477437 -0.11673018439312616 … 0.020732377101093546 -0.345372365290185], tree_statistics = DynamicHMC.TreeStatisticsNUTS[DynamicHMC.TreeStatisticsNUTS(-115.05949158993661, 2, turning at positions -3:-6, 0.8993343556295831, 7, DynamicHMC.Directions(0xf4ace331)), DynamicHMC.TreeStatisticsNUTS(-114.42936934449013, 3, turning at positions -3:4, 0.9669777811571946, 7, DynamicHMC.Directions(0x158f7ab4)), DynamicHMC.TreeStatisticsNUTS(-114.45712597732741, 2, turning at positions -3:0, 0.9065055732770694, 3, DynamicHMC.Directions(0x3ef3e86c)), DynamicHMC.TreeStatisticsNUTS(-115.67548379359854, 3, turning at positions -11:-14, 0.9591471691415301, 15, DynamicHMC.Directions(0x48b57331)), DynamicHMC.TreeStatisticsNUTS(-114.71867434172422, 2, turning at positions 2:5, 0.9999999999999999, 7, DynamicHMC.Directions(0x1b5de87d)), DynamicHMC.TreeStatisticsNUTS(-113.86635846149113, 4, turning at positions -4:-11, 1.0, 23, DynamicHMC.Directions(0x8281804c)), DynamicHMC.TreeStatisticsNUTS(-115.20547681499033, 2, turning at positions -3:0, 0.8474531210618542, 3, DynamicHMC.Directions(0x3374c5a8)), DynamicHMC.TreeStatisticsNUTS(-115.869090752354, 4, turning at positions -22:-25, 0.9764999663842461, 27, DynamicHMC.Directions(0x234ab542)), DynamicHMC.TreeStatisticsNUTS(-112.8590960700354, 2, turning at positions 2:5, 0.991553445280433, 7, DynamicHMC.Directions(0x8e0dc01d)), DynamicHMC.TreeStatisticsNUTS(-114.54672766286994, 5, turning at positions -11:-18, 0.9539242027685929, 39, DynamicHMC.Directions(0x5c4287d5))  …  DynamicHMC.TreeStatisticsNUTS(-119.91623068268446, 4, turning at positions -12:-19, 0.9377271436161124, 23, DynamicHMC.Directions(0x886ca684)), DynamicHMC.TreeStatisticsNUTS(-116.4525985016515, 3, turning at positions -3:4, 0.9999999999999999, 7, DynamicHMC.Directions(0x37b648fc)), DynamicHMC.TreeStatisticsNUTS(-119.12167494397156, 3, turning at positions 2:9, 0.9654665302371433, 15, DynamicHMC.Directions(0x1290f149)), DynamicHMC.TreeStatisticsNUTS(-117.11622759329698, 2, turning at positions 2:5, 0.9883790621197815, 7, DynamicHMC.Directions(0xa920d6b5)), DynamicHMC.TreeStatisticsNUTS(-116.00237838968259, 3, turning at positions -1:6, 0.9638578430968537, 7, DynamicHMC.Directions(0xd8aefcee)), DynamicHMC.TreeStatisticsNUTS(-121.50079499317714, 4, turning at positions 12:15, 0.7495020126242702, 19, DynamicHMC.Directions(0xee68e6db)), DynamicHMC.TreeStatisticsNUTS(-117.49565190032361, 3, turning at positions -6:1, 0.9999999999999999, 7, DynamicHMC.Directions(0x1fb64869)), DynamicHMC.TreeStatisticsNUTS(-116.57432958834991, 3, turning at positions 0:7, 0.9922393679813691, 7, DynamicHMC.Directions(0xe75d4b17)), DynamicHMC.TreeStatisticsNUTS(-117.2387033855675, 3, turning at positions 8:15, 0.8496432307002232, 15, DynamicHMC.Directions(0x9c9c8e0f)), DynamicHMC.TreeStatisticsNUTS(-115.26651386508914, 3, turning at positions -4:-11, 0.9281001866080739, 15, DynamicHMC.Directions(0x4f30bb44))], logdensities = [-113.27589148339064, -113.04553596982849, -113.35478675550712, -114.21779011666311, -113.48648539441253, -112.82466076072542, -114.74600750938598, -112.61042810892464, -112.30618125972252, -114.233927032386  …  -116.01800781876888, -115.8650731563811, -115.64868826582087, -114.21177048596996, -114.9773658930544, -115.99685044553503, -114.95522092925599, -114.2491133465807, -113.75971625943382, -113.73095463492326], κ = Gaussian kinetic energy (WoodburyPDMat), √diag(M⁻¹): [0.07049531136829397, 0.16078096078967372, 0.19731291336757925, 0.8625295954364424, 0.5151161960028614], ϵ = 0.7634060459179786)

AdvancedHMC.jl

Similar to Pathfinder and DynamicHMC, AdvancedHMC can also work with a LogDensityProblem:

using AdvancedHMC

nadapts = 500;

Initializing from Pathfinder's draws

metric = DiagEuclideanMetric(dim)
hamiltonian = Hamiltonian(metric, ∇P)
ϵ = find_good_stepsize(hamiltonian, init_params)
integrator = Leapfrog(ϵ)
kernel = HMCKernel(Trajectory{MultinomialTS}(integrator, GeneralisedNoUTurn()))
adaptor = StepSizeAdaptor(0.8, integrator)
samples_ahmc1, stats_ahmc1 = sample(
    hamiltonian,
    kernel,
    init_params,
    ndraws + nadapts,
    adaptor,
    nadapts;
    drop_warmup=true,
    progress=false,
)

([[-0.3458358942284875, 0.26577960545565427, -0.5742436205393942, -0.14231631489795032, 0.4887459004158137], [-0.32064939802679954, 0.317134376619134, 0.5702720008789725, -0.8886819346747703, 0.5737515634225554], [-0.3353112544828748, 0.1273059928292335, -1.4140238355909558, 0.3548200548227511, 0.5163836372377264], [-0.2966816810888126, 0.19824032297417965, 0.1300140090817286, 0.6721035715950026, -0.15155780619248718], [-0.2866011916194203, 0.18793109630406504, 0.1358597702676998, 0.6403484107570575, -0.1810245263810895], [-0.2490274300971668, 0.264910085426941, 0.22579917445388198, 0.6141406968834296, -0.2080800295938152], [-0.29532760437634165, 0.1825819856441852, 0.0031610946869207145, 0.20558649007608426, 0.14326404714140312], [-0.30989702391372165, 0.15954932593760263, 0.6019596987403389, -0.9446706078653209, 0.7580484028533889], [-0.3561329452303376, 0.2388953084760252, -0.023237117952837193, -1.2394241694033452, 1.1461842891418974], [-0.3620508816371867, 0.146972326714781, 0.0461426504047245, -1.069120940080641, 1.0650799268035291]  …  [-0.40482977964001987, 0.35781486979798904, -1.6105050636946039, -0.8371118244313375, 1.4078221077341677], [-0.2791843491254503, 0.12508093398448966, -1.1344042523546811, -0.7482605021697061, 1.1458817449795087], [-0.3082827516905339, 0.06271643421430918, -1.3498852921883846, -0.6058700219379757, 1.271646210432331], [-0.38606199069579816, 0.12367994563249744, -2.143373357781447, 0.24085686890034624, 0.8747732667408151], [-0.34782945118864717, 0.14154094234628112, -0.7619240892689018, 0.9784246122166438, -0.053687919546340696], [-0.41375421848671384, 0.012033358898293382, -0.6197469732100194, 1.0128661553584761, -0.054682092483315986], [-0.4040334552028756, 0.03552451700784978, -0.7025180776630265, 1.0076193909900164, -0.0803785165195383], [-0.29644027048860394, 0.27965660940056547, -0.6744061778769709, 0.9671088565482747, -0.16936350697882485], [-0.3504444054872439, 0.29319762734312954, 1.6737692165364, -0.866712725490039, 0.29021843596237223], [-0.38185169897668847, 0.22996558356122157, 0.2063036952836244, 0.7518771754385023, -0.29735038512747003]], NamedTuple[(n_steps = 55, is_accept = true, acceptance_rate = 0.7974205069874294, log_density = -113.46385668747193, hamiltonian_energy = 116.09995420280134, hamiltonian_energy_error = 0.401760005250253, max_hamiltonian_energy_error = 0.534660569909363, tree_depth = 5, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.9635652005696805, log_density = -114.24987642436467, hamiltonian_energy = 115.02782994363821, hamiltonian_energy_error = 0.0458679346491806, max_hamiltonian_energy_error = -0.5504992083511979, tree_depth = 6, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.9941057356353268, log_density = -114.37260187510694, hamiltonian_energy = 116.74082494783856, hamiltonian_energy_error = -0.15125220890863034, max_hamiltonian_energy_error = -0.5732445253818241, tree_depth = 6, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.9773595347457686, log_density = -113.15912103743759, hamiltonian_energy = 115.08469558338784, hamiltonian_energy_error = -0.1158174631000719, max_hamiltonian_energy_error = -0.429275318503727, tree_depth = 6, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 0.8672074719287919, log_density = -112.7932595886393, hamiltonian_energy = 113.8429063473409, hamiltonian_energy_error = -0.17691492960722144, max_hamiltonian_energy_error = 0.449360313623302, tree_depth = 3, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 127, is_accept = true, acceptance_rate = 0.9817192346482566, log_density = -113.33559547574269, hamiltonian_energy = 114.54246621064934, hamiltonian_energy_error = 0.04847811881107589, max_hamiltonian_energy_error = 0.08424675806101334, tree_depth = 6, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 43, is_accept = true, acceptance_rate = 0.9824091738720322, log_density = -112.62371653696434, hamiltonian_energy = 115.32480122205315, hamiltonian_energy_error = -0.013174715160957362, max_hamiltonian_energy_error = 0.13019862478597588, tree_depth = 5, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 47, is_accept = true, acceptance_rate = 0.8649742217537959, log_density = -113.36034370651352, hamiltonian_energy = 114.84071370511049, hamiltonian_energy_error = -0.049947812682546555, max_hamiltonian_energy_error = 0.3979525722020867, tree_depth = 5, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.8697553731236707, log_density = -113.87455980173206, hamiltonian_energy = 114.90220216092084, hamiltonian_energy_error = 0.13274806911583426, max_hamiltonian_energy_error = 0.3388701461831545, tree_depth = 6, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 39, is_accept = true, acceptance_rate = 0.9436830250136916, log_density = -114.11414872549598, hamiltonian_energy = 115.25154072767128, hamiltonian_energy_error = 0.08217368709141226, max_hamiltonian_energy_error = 0.19847493264245486, tree_depth = 5, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false)  …  (n_steps = 79, is_accept = true, acceptance_rate = 0.8247701243152082, log_density = -118.48836077987899, hamiltonian_energy = 122.26087296336641, hamiltonian_energy_error = 0.7077340772986247, max_hamiltonian_energy_error = -0.8691316675417369, tree_depth = 6, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 35, is_accept = true, acceptance_rate = 0.9969689566144894, log_density = -115.54986200968953, hamiltonian_energy = 119.41804391059554, hamiltonian_energy_error = -0.835098746754511, max_hamiltonian_energy_error = -1.646690932881512, tree_depth = 5, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.7627706932543892, log_density = -115.75910452133604, hamiltonian_energy = 120.49394141247194, hamiltonian_energy_error = -0.4471568869072513, max_hamiltonian_energy_error = 0.9586427297120679, tree_depth = 5, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 127, is_accept = true, acceptance_rate = 0.9378576799871178, log_density = -115.34677391415198, hamiltonian_energy = 117.80997925383086, hamiltonian_energy_error = -0.34543862096307976, max_hamiltonian_energy_error = -0.3907318128318451, tree_depth = 6, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.9993436389457647, log_density = -113.22089648682747, hamiltonian_energy = 116.30457102199803, hamiltonian_energy_error = -0.011035870728036912, max_hamiltonian_energy_error = -0.055944422454587084, tree_depth = 6, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 0.9132353006489439, log_density = -114.99921192354819, hamiltonian_energy = 116.54001544857243, hamiltonian_energy_error = 0.08532063008527757, max_hamiltonian_energy_error = 0.1716370838074539, tree_depth = 3, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 23, is_accept = true, acceptance_rate = 0.7371344030551816, log_density = -114.51869334395836, hamiltonian_energy = 116.3888661132698, hamiltonian_energy_error = 0.06248149447310425, max_hamiltonian_energy_error = 0.6638253717192839, tree_depth = 4, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 0.7535274084147162, log_density = -113.09182503277364, hamiltonian_energy = 115.79036373295997, hamiltonian_energy_error = 0.016889358560888468, max_hamiltonian_energy_error = 0.5879355582158468, tree_depth = 3, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 95, is_accept = true, acceptance_rate = 0.979773167653382, log_density = -114.06186957179256, hamiltonian_energy = 114.634380042959, hamiltonian_energy_error = 0.03797368877302176, max_hamiltonian_energy_error = 0.06908240369901364, tree_depth = 6, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.8548399186844843, log_density = -112.51365492894948, hamiltonian_energy = 115.91535598146493, hamiltonian_energy_error = -0.09095844052112056, max_hamiltonian_energy_error = 0.387630386466185, tree_depth = 6, numerical_error = false, step_size = 0.040131052556286254, nom_step_size = 0.040131052556286254, is_adapt = false)])

Initializing metric adaptation from Pathfinder's estimate

We can't start with Pathfinder's inverse metric estimate directly. Instead we need to first extract its diagonal for a DiagonalEuclideanMetric or make it dense for a DenseEuclideanMetric:

metric = DenseEuclideanMetric(Matrix(inv_metric))
hamiltonian = Hamiltonian(metric, ∇P)
ϵ = find_good_stepsize(hamiltonian, init_params)
integrator = Leapfrog(ϵ)
kernel = HMCKernel(Trajectory{MultinomialTS}(integrator, GeneralisedNoUTurn()))
adaptor = StepSizeAdaptor(0.8, integrator)
samples_ahmc2, stats_ahmc2 = sample(
    hamiltonian,
    kernel,
    init_params,
    ndraws + nadapts,
    adaptor,
    nadapts;
    drop_warmup=true,
    progress=false,
)

([[-0.3744233845664814, 0.20589575597454146, -1.5637138072071357, 1.2081523993503351, 0.03914251273638314], [-0.3553524893935575, 0.2967049958381658, -0.3468336193421929, 0.24997550428225068, 0.2441010584376998], [-0.28529338521470815, 0.29854705167613593, -0.7876729875290569, 0.2026760873918081, 0.37625504018853567], [-0.3934185740495537, 0.2622079501382567, -0.9157835086337673, -0.07994584325695045, 0.6353093493487877], [-0.3194358482817694, 0.2937236268731123, -0.32684330177587045, 1.2199863718494668, -0.4970820642940947], [-0.39881872556096254, 0.30507126964099957, -0.3789715603117769, -0.8326373225840131, 0.9315328040088319], [-0.39881872556096254, 0.30507126964099957, -0.3789715603117769, -0.8326373225840131, 0.9315328040088319], [-0.3494620344317025, 0.2398752533266597, 0.2265251504138964, 0.25256073739494256, 0.02742957295527937], [-0.4465017445496079, 0.10178465596011077, 0.5444121200550913, -0.8463054355930886, 0.7538928139021418], [-0.2582236529108202, 0.2004408538119687, -0.12563747363684835, 1.2421330422542782, -0.5355594841336053]  …  [-0.2955091090931978, 0.1423537221787053, 0.45517142286612156, 0.30863000164833254, 0.0356886077073596], [-0.39512665551019294, 0.31927375773005, 0.511992878885136, 0.34837095680406593, -0.14390228753931594], [-0.486788899399558, 0.201190444158927, 1.0516261333990289, -1.426050213500536, 0.9692884381139577], [-0.2569732288816123, 0.19409032495274164, 0.3046588686295624, 1.498335357523853, -0.7917018390529333], [-0.36147949946771063, 0.005678348588198934, 1.9777254230339065, -1.327279282445068, 0.5587315883648603], [-0.2872368673243825, 0.36862376364826216, 1.6377197909714618, 0.43694263548987766, -0.5938561013490306], [-0.3234081417142246, 0.09828083577451718, 2.005727995106004, -0.5820087086720811, 0.010212859159416478], [-0.3325642254896341, 0.2662482903188057, 1.36394143246208, 0.729225777628138, -0.6496756041925245], [-0.3797800756107558, 0.421614681205582, 1.6015188980381871, -0.7604227891852504, 0.12855843404623657], [-0.3797800756107558, 0.421614681205582, 1.6015188980381871, -0.7604227891852504, 0.12855843404623657]], NamedTuple[(n_steps = 31, is_accept = true, acceptance_rate = 0.8914758116386297, log_density = -114.39731873513912, hamiltonian_energy = 118.6263060791208, hamiltonian_energy_error = -0.14530825502841083, max_hamiltonian_energy_error = 0.4336964590888783, tree_depth = 4, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 19, is_accept = true, acceptance_rate = 0.9474336801708986, log_density = -112.94284113336904, hamiltonian_energy = 115.63440391452775, hamiltonian_energy_error = -0.10504772656281602, max_hamiltonian_energy_error = 0.2922323478468911, tree_depth = 4, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 15, is_accept = true, acceptance_rate = 0.9930480593282239, log_density = -112.94542319064392, hamiltonian_energy = 113.90280026327652, hamiltonian_energy_error = -0.04745941657617436, max_hamiltonian_energy_error = -0.10491875560505548, tree_depth = 3, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.764091653869185, log_density = -112.90570220862561, hamiltonian_energy = 116.57212768947196, hamiltonian_energy_error = -0.14406753636447434, max_hamiltonian_energy_error = 0.7860168589184582, tree_depth = 2, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 23, is_accept = true, acceptance_rate = 0.8547445603045918, log_density = -113.63045277732844, hamiltonian_energy = 114.87673269471419, hamiltonian_energy_error = 0.2951955993817421, max_hamiltonian_energy_error = 0.31984228780606827, tree_depth = 4, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 15, is_accept = true, acceptance_rate = 0.8721946789288401, log_density = -113.2448646768983, hamiltonian_energy = 116.20045998567944, hamiltonian_energy_error = -0.10094710254948325, max_hamiltonian_energy_error = 0.41695830892329866, tree_depth = 3, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.2847054782601081, log_density = -113.2448646768983, hamiltonian_energy = 121.11404474546364, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 2.1689129707209958, tree_depth = 2, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 19, is_accept = true, acceptance_rate = 0.9285130068963812, log_density = -112.11092693453726, hamiltonian_energy = 114.26617611951437, hamiltonian_energy_error = -0.228962906367272, max_hamiltonian_energy_error = -0.22970566504488943, tree_depth = 4, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.5445905582705005, log_density = -114.72585334826675, hamiltonian_energy = 115.8200872396967, hamiltonian_energy_error = 0.5775234245215728, max_hamiltonian_energy_error = 0.6863786342540834, tree_depth = 2, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 11, is_accept = true, acceptance_rate = 0.8812146303701699, log_density = -114.35509850695978, hamiltonian_energy = 117.21523237992466, hamiltonian_energy_error = 0.04349925415800726, max_hamiltonian_energy_error = 0.5191624433220454, tree_depth = 3, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false)  …  (n_steps = 23, is_accept = true, acceptance_rate = 0.7206675049352078, log_density = -114.06077649025522, hamiltonian_energy = 117.46986937263816, hamiltonian_energy_error = 0.18151936800327917, max_hamiltonian_energy_error = 0.5845728568918247, tree_depth = 4, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 0.9689389609456293, log_density = -112.77157204180419, hamiltonian_energy = 114.85711273508495, hamiltonian_energy_error = -0.27437402817565726, max_hamiltonian_energy_error = -0.27437402817565726, tree_depth = 2, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 15, is_accept = true, acceptance_rate = 0.4479737030149065, log_density = -118.38347676454102, hamiltonian_energy = 120.99854625959136, hamiltonian_energy_error = 1.4291924665403997, max_hamiltonian_energy_error = 1.6632397746698473, tree_depth = 3, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 1.0, log_density = -114.4660368980589, hamiltonian_energy = 119.3231995042888, hamiltonian_energy_error = -0.8101691039312158, max_hamiltonian_energy_error = -0.8962431760434271, tree_depth = 2, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 19, is_accept = true, acceptance_rate = 0.6896306143758822, log_density = -117.68961380628784, hamiltonian_energy = 120.20185629870045, hamiltonian_energy_error = 0.3055967818958152, max_hamiltonian_energy_error = 1.0956922057325897, tree_depth = 4, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 1.0, log_density = -114.69260251836529, hamiltonian_energy = 117.29948452207596, hamiltonian_energy_error = -0.7269252762841347, max_hamiltonian_energy_error = -0.7617969801750775, tree_depth = 2, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.5545308784430069, log_density = -115.73663143508188, hamiltonian_energy = 118.89307435537367, hamiltonian_energy_error = 0.35156919277827114, max_hamiltonian_energy_error = 0.9611186709615538, tree_depth = 2, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 15, is_accept = true, acceptance_rate = 0.9843772536159937, log_density = -113.7422182217931, hamiltonian_energy = 117.5032291737333, hamiltonian_energy_error = -0.3318328632178975, max_hamiltonian_energy_error = -0.3384905629236812, tree_depth = 3, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.6300671273608408, log_density = -114.86038780492832, hamiltonian_energy = 116.00242437035044, hamiltonian_energy_error = 0.2878877748910327, max_hamiltonian_energy_error = 0.6076809306458983, tree_depth = 2, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.4922756026289467, log_density = -114.86038780492832, hamiltonian_energy = 120.25668055452755, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 1.3800246702443104, tree_depth = 2, numerical_error = false, step_size = 0.9338144777561775, nom_step_size = 0.9338144777561775, is_adapt = false)])

Use Pathfinder's metric estimate for sampling

To use Pathfinder's metric with no metric adaptation, we need to use Pathfinder's own Pathfinder.RankUpdateEuclideanMetric type, which just wraps our inverse metric estimate for use with AdvancedHMC:

nadapts = 75
metric = Pathfinder.RankUpdateEuclideanMetric(inv_metric)
hamiltonian = Hamiltonian(metric, ∇P)
ϵ = find_good_stepsize(hamiltonian, init_params)
integrator = Leapfrog(ϵ)
kernel = HMCKernel(Trajectory{MultinomialTS}(integrator, GeneralisedNoUTurn()))
adaptor = StepSizeAdaptor(0.8, integrator)
samples_ahmc3, stats_ahmc3 = sample(
    hamiltonian,
    kernel,
    init_params,
    ndraws + nadapts,
    adaptor,
    nadapts;
    drop_warmup=true,
    progress=false,
)

([[-0.43852913717814335, 0.2447463158796643, 1.2798858940705204, -0.9671220517160993, 0.4990732469530218], [-0.38848723385475387, 0.08867550355478691, 1.4755969831043736, -0.9226867446886533, 0.48760673424014944], [-0.40478392293834137, 0.2690102943172251, 1.3751398900725778, -0.7204652273978505, 0.24536483755245403], [-0.40844501276669365, 0.24636691760536186, 1.1660496923519985, 0.5548247637546115, -0.4517258896569402], [-0.3132389762992508, 0.25427029786839284, -1.6242592284472785, 0.04381032295649062, 0.7672731889981518], [-0.4368321381399493, 0.2721641087369217, -1.5497271226079006, 0.006173273077439388, 0.8460588076139506], [-0.3488451259583163, 0.22952119347560976, -1.4757336634803553, -0.6302633329280755, 1.2042512849671543], [-0.34371816778844677, 0.23517520572760672, -1.2723437289097672, 0.4135338852281372, 0.43181004154973485], [-0.3574256423725241, 0.32966138923346044, 1.2440422890562928, 0.6708175037308084, -0.6459502508555482], [-0.3735624703393793, 0.16359464369278887, -0.3877055951685178, -0.5245605026832315, 0.7976922059732151]  …  [-0.4522227408114806, 0.1646682509879295, 0.5999348916095613, -1.7229747342996615, 1.2409557578749484], [-0.3440545732536459, 0.16038680792709575, -0.3239847108211308, 1.9609467315801279, -0.8516658562662185], [-0.3228848703240608, 0.14862053257719632, -0.038804370178773806, 0.0749332471360622, 0.29771459019987523], [-0.3502439262153331, 0.11343579028327032, -0.1729672176805908, -0.4302913235749198, 0.6712163633412604], [-0.30002354979218415, 0.29059641331974495, -0.21521830801863506, -0.7065963698013069, 0.7680263222084696], [-0.34762449978019894, 0.395696671849239, -0.1513582722091446, 1.607548571201352, -0.8399883725137138], [-0.3261141889959478, 0.281234195274528, -0.13801863983824111, -1.101920264238168, 1.0544631702326375], [-0.3938494029089846, 0.3573454265714501, 0.10381714876827786, 0.5912099772942094, -0.22609241899540494], [-0.4189847140336396, 0.30848517591739527, 0.13121511023835503, 1.3837803163921325, -0.7236069478035551], [-0.20934101768931632, 0.27746813685313043, 0.8045399632748923, -0.9683667207781124, 0.5865964547416206]], NamedTuple[(n_steps = 7, is_accept = true, acceptance_rate = 1.0, log_density = -114.28084677397293, hamiltonian_energy = 117.56472879655175, hamiltonian_energy_error = -0.29757807527732894, max_hamiltonian_energy_error = -0.488086544471912, tree_depth = 2, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 0.9083097528230691, log_density = -114.57491393582815, hamiltonian_energy = 116.71061063852595, hamiltonian_energy_error = 0.029744385942748863, max_hamiltonian_energy_error = 0.26799033578876674, tree_depth = 2, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.9821909191090722, log_density = -114.07776481938906, hamiltonian_energy = 115.59134152343762, hamiltonian_energy_error = -0.08684038015725548, max_hamiltonian_energy_error = -0.12997105475714932, tree_depth = 2, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.7901896103205605, log_density = -113.72263597044727, hamiltonian_energy = 117.27017753575824, hamiltonian_energy_error = -0.06657477115000177, max_hamiltonian_energy_error = 0.4038718845884546, tree_depth = 2, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 35, is_accept = true, acceptance_rate = 0.9571538508414204, log_density = -113.86165180344403, hamiltonian_energy = 115.8431845550255, hamiltonian_energy_error = -0.07056832764114063, max_hamiltonian_energy_error = 0.15043635002564315, tree_depth = 5, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.8439005357722799, log_density = -116.38987759314216, hamiltonian_energy = 116.58177518115173, hamiltonian_energy_error = 0.40533268784095355, max_hamiltonian_energy_error = 0.40533268784095355, tree_depth = 2, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 1.0, log_density = -114.10461971109655, hamiltonian_energy = 116.66209363465711, hamiltonian_energy_error = -0.39787838069294423, max_hamiltonian_energy_error = -0.40795371299610395, tree_depth = 2, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 15, is_accept = true, acceptance_rate = 0.9944284413574209, log_density = -113.05173363638174, hamiltonian_energy = 115.6670086723513, hamiltonian_energy_error = -0.054845488755617566, max_hamiltonian_energy_error = -0.05666866538244619, tree_depth = 3, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 35, is_accept = true, acceptance_rate = 0.9396315740143856, log_density = -113.51722466506004, hamiltonian_energy = 114.3724084985679, hamiltonian_energy_error = 0.05811362590483782, max_hamiltonian_energy_error = 0.12827362269891296, tree_depth = 5, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.9802643003434056, log_density = -112.81729113849296, hamiltonian_energy = 116.67259742956277, hamiltonian_energy_error = 0.007787614160946532, max_hamiltonian_energy_error = 0.07255870872089076, tree_depth = 4, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false)  …  (n_steps = 23, is_accept = true, acceptance_rate = 0.8484786889965857, log_density = -115.94340843668695, hamiltonian_energy = 119.69669419635859, hamiltonian_energy_error = 0.19333452785495808, max_hamiltonian_energy_error = 0.37001326650015187, tree_depth = 4, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 27, is_accept = true, acceptance_rate = 0.9194519033460226, log_density = -114.89383183816824, hamiltonian_energy = 118.98501346698508, hamiltonian_energy_error = -0.11765603477850561, max_hamiltonian_energy_error = -0.35739167024519247, tree_depth = 4, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 0.9901581481487075, log_density = -112.46242147032105, hamiltonian_energy = 115.20033056030321, hamiltonian_energy_error = -0.3777853631622605, max_hamiltonian_energy_error = -0.42146080965005694, tree_depth = 2, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.9557205140659444, log_density = -112.87960396828043, hamiltonian_energy = 113.50637719744607, hamiltonian_energy_error = 0.056485746795843284, max_hamiltonian_energy_error = 0.12688563109844608, tree_depth = 4, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 23, is_accept = true, acceptance_rate = 0.9781905652251018, log_density = -113.10381236408672, hamiltonian_energy = 113.70395298984178, hamiltonian_energy_error = 0.019695042225990278, max_hamiltonian_energy_error = -0.08660550908410869, tree_depth = 4, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.8594500083316888, log_density = -114.38443104556667, hamiltonian_energy = 116.14458313783913, hamiltonian_energy_error = 0.19613559472480802, max_hamiltonian_energy_error = 0.34117261103996555, tree_depth = 4, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 15, is_accept = true, acceptance_rate = 0.9689682310477726, log_density = -113.41943252053339, hamiltonian_energy = 115.99631754737538, hamiltonian_energy_error = -0.1469780457453993, max_hamiltonian_energy_error = -0.30716043455124975, tree_depth = 3, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.9980339923350239, log_density = -112.86330704960085, hamiltonian_energy = 114.91770846439579, hamiltonian_energy_error = -0.06983607013066262, max_hamiltonian_energy_error = -0.1599237309989121, tree_depth = 4, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.9199187723964709, log_density = -113.8668902280801, hamiltonian_energy = 114.46728055252301, hamiltonian_energy_error = 0.1353867016573247, max_hamiltonian_energy_error = 0.22299978873310522, tree_depth = 5, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false), (n_steps = 11, is_accept = true, acceptance_rate = 0.829870731193343, log_density = -115.43564364851827, hamiltonian_energy = 118.03207899836775, hamiltonian_energy_error = 0.23812291743651315, max_hamiltonian_energy_error = 0.514060204870944, tree_depth = 3, numerical_error = false, step_size = 0.7618635047601371, nom_step_size = 0.7618635047601371, is_adapt = false)])