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.35193742737501355, 0.24690060220202914, 0.05827588239230507, 0.1165517598296827, 0.1748276485360773]
Σ: [0.004969260325419835 -0.00022439878810641406 … 0.0005188996329265915 -0.00023506395820184131; -0.00022439878810642013 0.02586718625550903 … -0.0177998683150072 -0.00019797809456924796; … ; 0.0005188996329265512 -0.017799868315007177 … 0.7439352913289171 -0.43849537809019495; -0.0002350639582018271 -0.00019797809456919592 … -0.43849537809019495 0.2653368500944556]
)

init_params = result_pf.draws[:, 1]

5-element Vector{Float64}:
 -0.24466990408861164
  0.2908259603073542
  0.5016560101883027
 -1.6244919269371485
  1.1758786640700607

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.00496926   -0.000224399  -4.08535e-5   0.0005189  -0.000235064
 -0.000224399   0.0258672    -0.00139475  -0.0177999  -0.000197978
 -4.08535e-5   -0.00139475    0.038944    -0.145762    0.0852677
  0.0005189    -0.0177999    -0.145762     0.743935   -0.438495
 -0.000235064  -0.000197978   0.0852677   -0.438495    0.265337

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.36659943199900774 -0.28377718936871 … -0.3071816879625347 -0.2830738493846023; 0.2661870825077825 0.24835318486678934 … 0.20445559052445647 0.20836936120039706; … ; -0.31259080661506183 -0.492492037066512 … 0.3783830301132299 -0.2089911107465618; 0.38706634080943925 0.4738169392910579 … -0.3845837386045893 0.7702570313906439], tree_statistics = DynamicHMC.TreeStatisticsNUTS[DynamicHMC.TreeStatisticsNUTS(-115.10302911428727, 7, turning at positions -98:29, 0.973083040824819, 127, DynamicHMC.Directions(0x983ba09d)), DynamicHMC.TreeStatisticsNUTS(-116.86198623262455, 5, turning at positions 24:27, 0.5747358587383129, 35, DynamicHMC.Directions(0xc79762b7)), DynamicHMC.TreeStatisticsNUTS(-115.63847644628811, 6, turning at positions 0:63, 0.926819828700998, 63, DynamicHMC.Directions(0x4b9fb8bf)), DynamicHMC.TreeStatisticsNUTS(-119.23912354345026, 6, turning at positions -27:36, 0.8985477306472753, 63, DynamicHMC.Directions(0xae7de6a4)), DynamicHMC.TreeStatisticsNUTS(-122.85372402740389, 7, turning at positions -127:0, 0.9393527209513303, 127, DynamicHMC.Directions(0xfb1c2280)), DynamicHMC.TreeStatisticsNUTS(-120.15686276182612, 6, turning at positions -19:44, 0.9683455954495853, 63, DynamicHMC.Directions(0x535af5ec)), DynamicHMC.TreeStatisticsNUTS(-115.22815608200246, 6, turning at positions -55:8, 0.9996197410486959, 63, DynamicHMC.Directions(0x17063348)), DynamicHMC.TreeStatisticsNUTS(-116.40359785414745, 6, turning at positions -33:-96, 0.9616615002234465, 127, DynamicHMC.Directions(0xf2947d1f)), DynamicHMC.TreeStatisticsNUTS(-117.56190049807923, 6, turning at positions -36:-99, 0.9310943492312304, 127, DynamicHMC.Directions(0xc12bb81c)), DynamicHMC.TreeStatisticsNUTS(-115.84939214078962, 6, turning at positions -24:39, 0.9980454889044957, 63, DynamicHMC.Directions(0x89f83de7))  …  DynamicHMC.TreeStatisticsNUTS(-115.58294095853034, 6, turning at positions -27:-90, 0.968998059164561, 127, DynamicHMC.Directions(0xf61df1a5)), DynamicHMC.TreeStatisticsNUTS(-116.49799114136601, 6, turning at positions -2:61, 0.8637081146656936, 63, DynamicHMC.Directions(0xaf68f0bd)), DynamicHMC.TreeStatisticsNUTS(-114.75370049253777, 6, turning at positions -21:42, 0.9984332464960698, 63, DynamicHMC.Directions(0x05e6b4ea)), DynamicHMC.TreeStatisticsNUTS(-114.96990802146773, 6, turning at positions -24:39, 0.9914907406767843, 63, DynamicHMC.Directions(0x31b49927)), DynamicHMC.TreeStatisticsNUTS(-115.59263026870481, 6, turning at positions 35:98, 0.9977643619342167, 127, DynamicHMC.Directions(0xfca253e2)), DynamicHMC.TreeStatisticsNUTS(-117.60241075847627, 6, turning at positions -59:4, 0.9562898769246392, 63, DynamicHMC.Directions(0x93aca744)), DynamicHMC.TreeStatisticsNUTS(-116.69252351500293, 5, turning at positions -15:16, 0.9681241961630677, 31, DynamicHMC.Directions(0x21013e50)), DynamicHMC.TreeStatisticsNUTS(-118.6962152409409, 6, turning at positions -47:16, 0.8321196449377951, 63, DynamicHMC.Directions(0xbaae5e10)), DynamicHMC.TreeStatisticsNUTS(-116.80535019694102, 6, turning at positions -40:23, 0.9996254767679531, 63, DynamicHMC.Directions(0x5b8f75d7)), DynamicHMC.TreeStatisticsNUTS(-114.68533922884298, 7, turning at positions -2:125, 0.9769575766546554, 127, DynamicHMC.Directions(0x41273f7d))], logdensities = [-112.23632768661639, -112.9233192510133, -114.5666241620301, -118.5919920040288, -118.83714359957229, -114.74136577747774, -113.58747743232472, -114.00195614243628, -115.29056392597111, -113.86048947039374  …  -113.81894694483357, -114.30251322358794, -113.82230903884145, -112.4708070672628, -114.82218526910953, -115.71402447210737, -113.63011802645711, -116.02991259633545, -114.0919337552539, -113.56477595084887], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.07949241285913249, 0.14384813751280245, 0.9256505526570562, 0.9077252522400574, 0.6436850445259535], ϵ = 0.04452302397437675)

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.3701973983811422 -0.37938468433087674 … -0.3996554350739519 -0.3947835577456003; 0.0626277138424248 0.10369820258167403 … 0.0736532268056179 0.33396925553964035; … ; 0.471742526180727 0.48331101021611306 … -0.40748661387556495 0.3114169622957113; -0.1396448564774172 -0.18131730586397318 … 0.3067738840979057 -0.13565734721666473], tree_statistics = DynamicHMC.TreeStatisticsNUTS[DynamicHMC.TreeStatisticsNUTS(-115.43366353805492, 3, turning at positions -7:0, 0.9999999999999999, 7, DynamicHMC.Directions(0x2b5d9d40)), DynamicHMC.TreeStatisticsNUTS(-114.87939777782366, 2, turning at positions -3:0, 0.8204876951566576, 3, DynamicHMC.Directions(0x2d5eb4b8)), DynamicHMC.TreeStatisticsNUTS(-116.6019226143793, 3, turning at positions -2:5, 0.9947201332797962, 7, DynamicHMC.Directions(0x4119b11d)), DynamicHMC.TreeStatisticsNUTS(-115.72336589117255, 3, turning at positions -2:5, 0.9692112155051756, 7, DynamicHMC.Directions(0xdbd2d71d)), DynamicHMC.TreeStatisticsNUTS(-116.32938404958803, 4, turning at positions -8:7, 0.9214534122062351, 15, DynamicHMC.Directions(0xf7f25ac7)), DynamicHMC.TreeStatisticsNUTS(-115.74888075612655, 2, turning at positions -3:0, 0.9385893622347471, 3, DynamicHMC.Directions(0x4313dc0c)), DynamicHMC.TreeStatisticsNUTS(-116.68148854516426, 4, turning at positions -1:14, 0.8903939545161369, 15, DynamicHMC.Directions(0x38df760e)), DynamicHMC.TreeStatisticsNUTS(-119.70044552545107, 3, turning at positions 4:11, 0.9734979778043491, 15, DynamicHMC.Directions(0x2f2832cb)), DynamicHMC.TreeStatisticsNUTS(-119.6780124214395, 4, turning at positions -7:8, 0.7931778423617768, 15, DynamicHMC.Directions(0xd8812c38)), DynamicHMC.TreeStatisticsNUTS(-118.62321156755361, 3, turning at positions -1:6, 0.9143506151972965, 7, DynamicHMC.Directions(0x4b18143e))  …  DynamicHMC.TreeStatisticsNUTS(-118.06595718824167, 3, turning at positions -2:5, 0.9894645817446197, 7, DynamicHMC.Directions(0x837cb65d)), DynamicHMC.TreeStatisticsNUTS(-118.5833300357164, 4, turning at positions -5:10, 0.9753708965275825, 15, DynamicHMC.Directions(0x5f72e20a)), DynamicHMC.TreeStatisticsNUTS(-116.41353812390234, 4, turning at positions 0:15, 0.8449627601205768, 15, DynamicHMC.Directions(0x3e10917f)), DynamicHMC.TreeStatisticsNUTS(-117.77543316387639, 3, turning at positions 0:7, 0.9359497462854052, 7, DynamicHMC.Directions(0xad9da1ff)), DynamicHMC.TreeStatisticsNUTS(-119.38377424785634, 4, turning at positions -13:2, 0.9822910520519328, 15, DynamicHMC.Directions(0x376e5632)), DynamicHMC.TreeStatisticsNUTS(-116.45162282293164, 3, turning at positions -5:2, 0.9999999999999999, 7, DynamicHMC.Directions(0xdf02f81a)), DynamicHMC.TreeStatisticsNUTS(-118.97890714401484, 2, turning at positions 2:5, 0.8877334107124147, 7, DynamicHMC.Directions(0xe7bd721d)), DynamicHMC.TreeStatisticsNUTS(-114.03179711964616, 3, turning at positions -7:0, 0.9879173040235025, 7, DynamicHMC.Directions(0x263b4de0)), DynamicHMC.TreeStatisticsNUTS(-115.61485029760915, 3, turning at positions 0:7, 0.9815744430227087, 7, DynamicHMC.Directions(0x4a599b57)), DynamicHMC.TreeStatisticsNUTS(-114.93415552113994, 3, turning at positions -3:4, 0.9582022118210364, 7, DynamicHMC.Directions(0x02cafd74))], logdensities = [-113.27842823836157, -114.03891791903577, -114.57566063516852, -114.62797820090161, -114.00386998078572, -114.56946600758259, -114.57364673134585, -114.13299150873901, -116.36172786097615, -114.29380682810874  …  -116.86154313606418, -113.74286953285409, -113.93315182123287, -115.8676086324056, -115.6080266613828, -112.73118851749184, -112.89152604297617, -113.46564496167831, -113.83536066357901, -113.02808087470729], κ = Gaussian kinetic energy (Symmetric), √diag(M⁻¹): [0.07381828801778338, 0.12650585411964896, 1.0050163821142786, 0.9363489211997947, 0.6759210750169675], ϵ = 0.35346049048199285)

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.3799727829288582 -0.30702924016268895 … -0.2856044985013047 -0.1491290433873244; 0.15145322124191238 0.20043130892387553 … 0.2556564034965453 0.2824137640610929; … ; 0.5107936956314771 0.8057820265462009 … 0.11094930972338735 -0.8505896073213106; 0.32923957424962325 -0.39150248890147027 … 0.2797559013424403 0.7551581304798214], tree_statistics = DynamicHMC.TreeStatisticsNUTS[DynamicHMC.TreeStatisticsNUTS(-119.01080595577714, 2, turning at positions -1:2, 0.8443761900346621, 3, DynamicHMC.Directions(0x11d21912)), DynamicHMC.TreeStatisticsNUTS(-114.77677422037736, 5, turning at positions 30:33, 0.9395880622792583, 39, DynamicHMC.Directions(0xafe4acb9)), DynamicHMC.TreeStatisticsNUTS(-113.39207936573753, 2, turning at positions -2:1, 0.9945990955389569, 3, DynamicHMC.Directions(0xfde49301)), DynamicHMC.TreeStatisticsNUTS(-114.75539468720345, 3, turning at positions -5:2, 0.8775391555688469, 7, DynamicHMC.Directions(0x4675614a)), DynamicHMC.TreeStatisticsNUTS(-115.71214114931017, 4, turning at positions 11:18, 0.993284407345505, 23, DynamicHMC.Directions(0x4923dd9a)), DynamicHMC.TreeStatisticsNUTS(-114.90838737725214, 5, turning at positions -26:-33, 0.9286186566601875, 39, DynamicHMC.Directions(0xa766d3c6)), DynamicHMC.TreeStatisticsNUTS(-114.69934977795174, 2, turning at positions -2:-5, 0.9833084061950362, 7, DynamicHMC.Directions(0x9f242dea)), DynamicHMC.TreeStatisticsNUTS(-119.38929522325577, 4, turning at positions 19:22, 0.8149816390507921, 23, DynamicHMC.Directions(0xd5d44c9e)), DynamicHMC.TreeStatisticsNUTS(-121.29814638012965, 3, turning at positions -7:0, 0.7715474079148957, 7, DynamicHMC.Directions(0x24a3ace8)), DynamicHMC.TreeStatisticsNUTS(-117.10895312392859, 3, turning at positions -5:2, 0.9065931507969763, 7, DynamicHMC.Directions(0x78937a32))  …  DynamicHMC.TreeStatisticsNUTS(-116.08959750988545, 3, turning at positions 0:7, 0.9999999999999999, 7, DynamicHMC.Directions(0xd251cde7)), DynamicHMC.TreeStatisticsNUTS(-117.57147398520104, 4, turning at positions 10:17, 0.9455738191731414, 31, DynamicHMC.Directions(0x2cdc6cd1)), DynamicHMC.TreeStatisticsNUTS(-120.28353942792185, 3, turning at positions -8:-15, 0.8796816259731949, 15, DynamicHMC.Directions(0xc6deb240)), DynamicHMC.TreeStatisticsNUTS(-115.8085790720673, 3, turning at positions -4:3, 0.9999999999999999, 7, DynamicHMC.Directions(0xf3e2489b)), DynamicHMC.TreeStatisticsNUTS(-114.41415259278949, 3, turning at positions -5:2, 0.9635401716688629, 7, DynamicHMC.Directions(0x8a2b9252)), DynamicHMC.TreeStatisticsNUTS(-115.72751121104815, 3, turning at positions -3:4, 0.9113317424495618, 7, DynamicHMC.Directions(0xf2e027f4)), DynamicHMC.TreeStatisticsNUTS(-124.31151976008302, 2, turning at positions 0:3, 0.472683092923128, 3, DynamicHMC.Directions(0x117ac1f3)), DynamicHMC.TreeStatisticsNUTS(-121.36065299530152, 4, turning at positions -25:-28, 0.9820469213205849, 31, DynamicHMC.Directions(0x7d34f843)), DynamicHMC.TreeStatisticsNUTS(-116.74917403130189, 3, turning at positions -2:5, 0.886082809249989, 7, DynamicHMC.Directions(0xd7108165)), DynamicHMC.TreeStatisticsNUTS(-118.17927529923242, 2, turning at positions -1:-4, 0.7592303402330097, 7, DynamicHMC.Directions(0x92c59a23))], logdensities = [-113.10888536769221, -113.18406730255994, -112.38507962598713, -114.23794286601193, -112.90067030777837, -114.04352447842352, -113.45744293819934, -117.21288585489604, -115.26871517826716, -115.35438484177513  …  -114.49734198890053, -115.54323413421007, -115.54861630006297, -113.52026241104248, -113.81416894270016, -114.46365156632652, -119.3334163281969, -113.74121018510141, -113.44129257640427, -116.8174158044207], κ = Gaussian kinetic energy (WoodburyPDMat), √diag(M⁻¹): [0.07049298068190786, 0.16083278973986936, 0.19734223780493018, 0.862516835388688, 0.5151085808782995], ϵ = 0.70298741679258)

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.3870400041227244, 0.3653591764636401, 0.02919121059605427, -0.13793610101552425, 0.25491921982811944], [-0.3775684849205111, 0.33955742076229106, 0.055851031365678525, -0.15698357876317537, 0.34833150004615276], [-0.3086476937033776, 0.4256508819324544, -0.050895281217571094, -0.5767065865062615, 0.5502748687886092], [-0.23012029091250386, 0.08443817504099631, 0.4247194743090758, -0.8005403236623807, 0.7643153063660049], [-0.38646085253627366, 0.10243208826543998, 0.9204627371391044, 0.4794100380166985, -0.26237276603653786], [-0.4061422639250686, 0.14588518249156052, 0.800837845276798, -0.10052731271190943, 0.1327836203301695], [-0.3741594335620468, 0.2190886343992991, -0.6771423927437, 0.33907845502231115, 0.23521737601177756], [-0.32372736667188395, 0.110570924271601, -1.2664889817179383, 1.4146760190169922, -0.20215704343338858], [-0.28284422206670956, 0.11050295736543814, -1.2342751535181629, 1.4373356384837, -0.1635966323206412], [-0.43395713430280614, 0.41238451883894095, -1.559758906316286, 0.1821008884406956, 0.5916675957789879]  …  [-0.34612829451055793, 0.16306042646902694, -0.3663934691866386, -0.05350778508506867, 0.46308377642062704], [-0.3079248347329725, 0.3185549354587214, 0.1003612220706796, -0.14384677950012706, 0.2632752667182644], [-0.3149180249495505, 0.287609374979178, 0.2126308806063746, 1.3356475365860057, -0.7628288857572288], [-0.37410016192540985, 0.16131994011792472, 0.678196971296346, -0.5014682447014898, 0.4119082788647123], [-0.3728475418674468, 0.33587375020001925, 0.503279932008085, -0.5759195350751393, 0.45035251477762483], [-0.33038543829130484, 0.21339860746020428, -0.2618362765117194, -0.17343096338188152, 0.4942830797758483], [-0.27474899369071093, 0.2622374180423606, -0.4307335115454139, 0.06285962618940352, 0.36332778871665955], [-0.38571750169936264, 0.32609518148297606, 0.40384021376668755, 0.859987234134923, -0.41951306212007483], [-0.26468255666064117, 0.39771985854261754, 1.4224575945828155, 0.10928925288155722, -0.44149807436303046], [-0.30986258764576613, 0.14787284763063488, -0.7072941466576084, -0.6882410681598596, 0.9591122594872731]], NamedTuple[(n_steps = 39, is_accept = true, acceptance_rate = 0.8631127403923833, log_density = -113.16594689590295, hamiltonian_energy = 117.06643593044501, hamiltonian_energy_error = 0.18174264039218713, max_hamiltonian_energy_error = 0.31593163906377697, tree_depth = 5, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.992720202720335, log_density = -112.85685843359632, hamiltonian_energy = 113.72159492107446, hamiltonian_energy_error = -0.010043013384830601, max_hamiltonian_energy_error = -0.22084265470449793, tree_depth = 4, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 47, is_accept = true, acceptance_rate = 1.0, log_density = -113.67694154061755, hamiltonian_energy = 115.71833387036395, hamiltonian_energy_error = -0.08677390455581246, max_hamiltonian_energy_error = -0.2320584225379605, tree_depth = 5, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.5829483596188245, log_density = -114.75120469555239, hamiltonian_energy = 117.26549616647893, hamiltonian_energy_error = 0.2009289653199744, max_hamiltonian_energy_error = 1.4055800530510822, tree_depth = 6, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.9912630136475262, log_density = -113.45178795655734, hamiltonian_energy = 116.17724483028067, hamiltonian_energy_error = -0.034864302553145876, max_hamiltonian_energy_error = -0.13030416338570205, tree_depth = 6, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.9786645263191941, log_density = -113.01644796387686, hamiltonian_energy = 113.87827363150686, hamiltonian_energy_error = -0.03330250268341217, max_hamiltonian_energy_error = 0.09042319008179334, tree_depth = 6, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 127, is_accept = true, acceptance_rate = 0.6950343287894325, log_density = -113.25112048106101, hamiltonian_energy = 114.78748860189123, hamiltonian_energy_error = 0.43284659216661225, max_hamiltonian_energy_error = 0.8020221560324927, tree_depth = 7, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.9693585318765344, log_density = -114.53052456571761, hamiltonian_energy = 115.41173321111529, hamiltonian_energy_error = -0.2679234828458874, max_hamiltonian_energy_error = -0.39294564442913327, tree_depth = 6, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 0.8432629990927004, log_density = -115.22706844627201, hamiltonian_energy = 115.92179954385799, hamiltonian_energy_error = 0.16728773372170735, max_hamiltonian_energy_error = 0.3129744987129328, tree_depth = 3, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 127, is_accept = true, acceptance_rate = 0.9313802691999397, log_density = -115.03447431843941, hamiltonian_energy = 117.60156508233784, hamiltonian_energy_error = -0.284808197987914, max_hamiltonian_energy_error = 0.44562728484154945, tree_depth = 6, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false)  …  (n_steps = 23, is_accept = true, acceptance_rate = 0.8934816662545918, log_density = -112.42379324094142, hamiltonian_energy = 114.86378750005211, hamiltonian_energy_error = -0.10842293188397889, max_hamiltonian_energy_error = 0.29987192465280543, tree_depth = 4, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.9195102989733916, log_density = -112.77777928989003, hamiltonian_energy = 113.76074525762431, hamiltonian_energy_error = 0.14494705432672106, max_hamiltonian_energy_error = 0.19783475052520316, tree_depth = 5, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.9262155154133918, log_density = -114.24931706461157, hamiltonian_energy = 115.91773043019597, hamiltonian_energy_error = 0.18363370189841532, max_hamiltonian_energy_error = 0.2636031661401006, tree_depth = 5, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 127, is_accept = true, acceptance_rate = 0.8835962459439473, log_density = -112.77177865117194, hamiltonian_energy = 115.12777126271264, hamiltonian_energy_error = -0.37895041500780735, max_hamiltonian_energy_error = 0.48009010056996715, tree_depth = 6, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 15, is_accept = true, acceptance_rate = 0.9902922271592095, log_density = -112.72159391682413, hamiltonian_energy = 113.09454898578183, hamiltonian_energy_error = -0.016198525202781866, max_hamiltonian_energy_error = 0.029261784932444357, tree_depth = 4, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.9923849349604309, log_density = -112.29051813932159, hamiltonian_energy = 113.30785605752237, hamiltonian_energy_error = -0.012905394108685186, max_hamiltonian_energy_error = 0.030865666197101405, tree_depth = 6, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 23, is_accept = true, acceptance_rate = 0.9621002363418837, log_density = -112.7879842907841, hamiltonian_energy = 113.56360968912446, hamiltonian_energy_error = 0.04355267582667466, max_hamiltonian_energy_error = 0.08328899135655377, tree_depth = 4, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 47, is_accept = true, acceptance_rate = 0.7336883562903593, log_density = -113.9034079830573, hamiltonian_energy = 117.07508907831954, hamiltonian_energy_error = 0.269678892834861, max_hamiltonian_energy_error = 0.9471576890981339, tree_depth = 5, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.7443145527104333, log_density = -116.58382815534837, hamiltonian_energy = 117.561423058912, hamiltonian_energy_error = 0.46937976189791186, max_hamiltonian_energy_error = 0.9328895990065718, tree_depth = 6, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, is_adapt = false), (n_steps = 63, is_accept = true, acceptance_rate = 0.7802753747650127, log_density = -114.31394714172932, hamiltonian_energy = 121.24863427882245, hamiltonian_energy_error = -0.3051823710187165, max_hamiltonian_energy_error = 2.1426079551829815, tree_depth = 6, numerical_error = false, step_size = 0.04019718899483467, nom_step_size = 0.04019718899483467, 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.48405909119090584, 0.222537263914171, 1.289665351771141, 1.389251655987481, -1.1196126249348963], [-0.18296556779501763, 0.30446317245739285, 1.8572933028128424, -1.2527727888154754, 0.4998341558954196], [-0.49288289419472475, 0.22663514366397997, 1.208561037348258, 1.4642921009481278, -1.1473753407800311], [-0.4565675104148432, 0.20712062269850534, 0.1248710195065324, 1.9379072201572842, -1.0323913180855768], [-0.2442595698820892, 0.3005249646611569, 0.605249240899271, -1.6114207862305503, 1.112321756335048], [-0.20000634840741888, -0.013223239267102643, -0.41326567962282534, -0.7408021889782792, 1.0066283194514], [-0.13439144663700747, 0.06932255008385889, -0.5841912852791198, -0.17605729633444478, 0.6833992125883097], [-0.5277023157352512, 0.4835270090841027, -0.8596863886179792, 0.31440696195245477, 0.21407956177626142], [-0.22639748142732644, -0.08284539435155625, -0.39004293446010996, -0.30780723980581404, 0.7605560752879135], [-0.40834690891764663, 0.014234827392397342, 0.8733181400767179, 1.735981721623575, -1.0897843755505767]  …  [-0.27244908873501705, 0.23870318343413746, -0.08562454769057083, 0.589919124055628, -0.0962771929433612], [-0.27336853513074205, 0.37540894167779093, 1.445632629673298, -0.3385389288096445, -0.022414613005506845], [-0.4056095165626381, 0.3328605638772755, 1.1144371027244417, -0.7720834873378256, 0.36994766027886494], [-0.2985000624214291, 0.31379222201082174, 0.3341025368022972, 0.7272564661840073, -0.34731215210102573], [-0.3213333251469311, 0.13107134929127387, 0.49690609313905787, 0.07642229964573777, 0.16126306104724053], [-0.3227692771794899, 0.29189872085081037, -1.5063155383363094, -0.7108300330239402, 1.247134724223113], [-0.3545507015478428, 0.2716355627021048, -1.5744489211526915, -0.5992463622836145, 1.098706727476086], [-0.35580789146705266, 0.2608830987242595, 1.939108727982254, -0.7914942140224719, 0.08420555110166239], [-0.3870235834758584, 0.3463502769882344, 1.9431186917087715, -0.5152262653916211, -0.1455319692119888], [-0.3352816424829154, 0.3763622478834259, 2.1991381953562517, -1.2489726631046607, 0.2686280173503506]], NamedTuple[(n_steps = 3, is_accept = true, acceptance_rate = 0.4449486646132428, log_density = -117.33635001601859, hamiltonian_energy = 120.82662734184441, hamiltonian_energy_error = 0.5200338004953551, max_hamiltonian_energy_error = 1.0748258246303664, tree_depth = 2, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 0.9953988120695966, log_density = -117.69201511722295, hamiltonian_energy = 118.40854652036467, hamiltonian_energy_error = 0.03273841684088552, max_hamiltonian_energy_error = -0.533315024377174, tree_depth = 2, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.9693624924206471, log_density = -117.77352973573596, hamiltonian_energy = 118.72961139933422, hamiltonian_energy_error = 0.09641456442444962, max_hamiltonian_energy_error = -0.6403366502212293, tree_depth = 2, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 23, is_accept = true, acceptance_rate = 0.8689670312958578, log_density = -115.731360977501, hamiltonian_energy = 120.73230176723011, hamiltonian_energy_error = -0.16641312776928885, max_hamiltonian_energy_error = 0.4650901067366675, tree_depth = 4, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 1.0, log_density = -115.31395769822008, hamiltonian_energy = 115.99639076419686, hamiltonian_energy_error = -0.09092730286459982, max_hamiltonian_energy_error = -0.6661822823494532, tree_depth = 2, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 27, is_accept = true, acceptance_rate = 0.7666833896264693, log_density = -116.48146300315842, hamiltonian_energy = 121.29349064590004, hamiltonian_energy_error = 0.2988381941811298, max_hamiltonian_energy_error = 1.0902644011881932, tree_depth = 4, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 0.9889768250084624, log_density = -117.16961714693166, hamiltonian_energy = 118.72879562127348, hamiltonian_energy_error = 0.08030181825894545, max_hamiltonian_energy_error = -0.4905860810841318, tree_depth = 2, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.924000359581962, log_density = -118.23025485439891, hamiltonian_energy = 120.58033645745441, hamiltonian_energy_error = 0.2587693316191064, max_hamiltonian_energy_error = -0.406799140469019, tree_depth = 2, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 1.0, log_density = -116.15991393434216, hamiltonian_energy = 119.72705257797557, hamiltonian_energy_error = -0.3565545078308503, max_hamiltonian_energy_error = -1.0803549690736958, tree_depth = 2, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.873038048642848, log_density = -116.73912211869315, hamiltonian_energy = 120.1793410455605, hamiltonian_energy_error = 0.044339882525548546, max_hamiltonian_energy_error = -0.7295532864038137, tree_depth = 4, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false)  …  (n_steps = 27, is_accept = true, acceptance_rate = 0.962611964727777, log_density = -112.8507473866307, hamiltonian_energy = 113.26887351770456, hamiltonian_energy_error = 0.08292669872945169, max_hamiltonian_energy_error = 0.08468926870142468, tree_depth = 4, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 27, is_accept = true, acceptance_rate = 0.8968476988167637, log_density = -114.2991158621813, hamiltonian_energy = 116.37644135504343, hamiltonian_energy_error = 0.12866469394882074, max_hamiltonian_energy_error = 0.2849916027413144, tree_depth = 4, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 15, is_accept = true, acceptance_rate = 0.9723736503906427, log_density = -113.56002572255296, hamiltonian_energy = 115.11379371603917, hamiltonian_energy_error = -0.051279164926398835, max_hamiltonian_energy_error = -0.12010974838473487, tree_depth = 4, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 15, is_accept = true, acceptance_rate = 1.0, log_density = -112.84536252773424, hamiltonian_energy = 113.66007116042483, hamiltonian_energy_error = -0.10004432798122309, max_hamiltonian_energy_error = -0.1753949944318549, tree_depth = 4, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.8226533498357673, log_density = -113.2667965382735, hamiltonian_energy = 114.42850864481557, hamiltonian_energy_error = 0.056490070166418604, max_hamiltonian_energy_error = 0.31648454033746987, tree_depth = 2, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 27, is_accept = true, acceptance_rate = 0.8016975616830301, log_density = -114.48125410746998, hamiltonian_energy = 117.78722857434228, hamiltonian_energy_error = -0.05203926016329774, max_hamiltonian_energy_error = 0.6043075444523822, tree_depth = 4, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.8027528393600494, log_density = -116.35130573594581, hamiltonian_energy = 117.17188411698363, hamiltonian_energy_error = 0.3447411691232105, max_hamiltonian_energy_error = 0.3447411691232105, tree_depth = 2, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.8877170282255072, log_density = -115.45710264851415, hamiltonian_energy = 120.18900273710416, hamiltonian_energy_error = -0.18968097224555436, max_hamiltonian_energy_error = 0.368783601899068, tree_depth = 5, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 0.9669967199843337, log_density = -115.84676916687265, hamiltonian_energy = 116.71038647891268, hamiltonian_energy_error = 0.10848223326229345, max_hamiltonian_energy_error = -0.2963428999397877, tree_depth = 2, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, is_adapt = false), (n_steps = 11, is_accept = true, acceptance_rate = 0.9849294383683103, log_density = -116.19609638066909, hamiltonian_energy = 117.81517248390028, hamiltonian_energy_error = 0.014415557355206943, max_hamiltonian_energy_error = -0.2704895167994863, tree_depth = 3, numerical_error = false, step_size = 0.9028587033521537, nom_step_size = 0.9028587033521537, 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.42282133658100696, 0.12146908470801988, -0.7384788409421505, -1.3669256520142943, 1.4764786839280144], [-0.3127739542303199, 0.2762772946173255, -0.5085893629697196, 0.47360865505930705, 0.19200469949374965], [-0.3710389602124542, 0.2516757595167619, -0.20748114789876612, 0.2394742808715117, 0.16886728876303347], [-0.30696283863905094, 0.31315525348786194, 0.18783649212676148, -1.83677222570517, 1.3671742793470707], [-0.3696180336326971, 0.07245462488111301, 0.5542684034389744, 0.8913021651529398, -0.38324526108214907], [-0.3220119561283321, 0.39714242457801, 0.8419584215187171, -0.5901292429878128, 0.27034640278708044], [-0.26908902654761563, 0.4198037824329849, 1.2236431235304481, 0.032701565998691895, -0.24674702771462853], [-0.24512148899816807, 0.20181750051780079, 0.06369626027649664, -0.40871568064592456, 0.5097264779263966], [-0.23896445139142714, 0.2005552848146755, 0.04970411503970722, 0.6115956239719916, -0.12945208358332794], [-0.3723967320870173, 0.43316845861198655, 0.1807593209125947, -0.7066762215068987, 0.5449620756237389]  …  [-0.3915819889097766, 0.5778443630051373, -1.3287614317261012, 0.9381383852081487, -0.10417620619188733], [-0.3278852710715162, 0.10357583515081356, -0.06298564143182009, 1.544605569890877, -0.723843654471936], [-0.317799004103397, 0.025180062587003615, 1.1850936325045356, 0.06682361415262317, -0.05552398103673872], [-0.3250359175239661, -0.07029515502994879, 1.3692728264847598, -1.0131050163733821, 0.6549897300872756], [-0.3334303461779773, 0.22037930804916647, -0.15124181533547468, 0.1835546033343507, 0.19955849104206924], [-0.4359246597232659, 0.1137121533021126, 0.06504527930163478, 0.8457010648241753, -0.24784726074938201], [-0.38667266029332437, -0.0077239451627573105, 1.1915465430021688, -0.5840970120595976, 0.3990494996952355], [-0.2508948533597599, 0.3102876383858316, 0.8109139217407667, 1.195971949946051, -0.8845700204020126], [-0.4271980014648589, 0.2407740441648927, 1.5192290684414307, -1.0606958825752923, 0.5211296907856442], [-0.314898647003706, 0.212773347754305, 1.5985224208947977, -0.9685980654279565, 0.398019371921441]], NamedTuple[(n_steps = 39, is_accept = true, acceptance_rate = 0.8702771760999344, log_density = -115.43067439837797, hamiltonian_energy = 116.83194561384168, hamiltonian_energy_error = 0.2415648320751842, max_hamiltonian_energy_error = 0.3046385631446924, tree_depth = 5, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 23, is_accept = true, acceptance_rate = 0.9183168056131119, log_density = -114.3291392142908, hamiltonian_energy = 117.86752256494718, hamiltonian_energy_error = -0.023136028498925043, max_hamiltonian_energy_error = 0.30091910780893727, tree_depth = 4, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 23, is_accept = true, acceptance_rate = 0.8857205240180513, log_density = -112.18379407237296, hamiltonian_energy = 116.86964948362694, hamiltonian_energy_error = -0.2685539551291356, max_hamiltonian_energy_error = 0.402119322204598, tree_depth = 4, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.517804470290237, log_density = -115.32122852401824, hamiltonian_energy = 118.65060276108247, hamiltonian_energy_error = 0.4780323941931215, max_hamiltonian_energy_error = 1.0186352796932425, tree_depth = 1, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 19, is_accept = true, acceptance_rate = 0.9956288654647422, log_density = -113.94874888318057, hamiltonian_energy = 116.58587130185994, hamiltonian_energy_error = -0.1906750612479584, max_hamiltonian_energy_error = -0.35343215563496244, tree_depth = 4, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 1.0, log_density = -113.73780952755145, hamiltonian_energy = 114.011940255037, hamiltonian_energy_error = -0.059499699158706676, max_hamiltonian_energy_error = -0.2426408477949309, tree_depth = 3, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 23, is_accept = true, acceptance_rate = 0.988536321644165, log_density = -114.1753088877625, hamiltonian_energy = 115.3925962967937, hamiltonian_energy_error = -0.00188947073570489, max_hamiltonian_energy_error = -0.13235910218011782, tree_depth = 4, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 39, is_accept = true, acceptance_rate = 0.9730197054369949, log_density = -113.70013400670591, hamiltonian_energy = 115.58780843856654, hamiltonian_energy_error = 0.0112264599601275, max_hamiltonian_energy_error = 0.17119374641356444, tree_depth = 5, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.9888339892438872, log_density = -113.48037615155722, hamiltonian_energy = 114.98467161616881, hamiltonian_energy_error = -0.010170109786258763, max_hamiltonian_energy_error = -0.12134886585667459, tree_depth = 4, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 23, is_accept = true, acceptance_rate = 0.9237220986247124, log_density = -114.13461884108308, hamiltonian_energy = 115.75771943839928, hamiltonian_energy_error = 0.11709381957302867, max_hamiltonian_energy_error = -0.20026432286238105, tree_depth = 4, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false)  …  (n_steps = 11, is_accept = true, acceptance_rate = 0.8500055187141654, log_density = -116.97254061445095, hamiltonian_energy = 119.04683349092585, hamiltonian_energy_error = 0.3516338021117633, max_hamiltonian_energy_error = 0.3780315381052901, tree_depth = 3, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.9958610961637856, log_density = -115.07283004418265, hamiltonian_energy = 119.00794276726604, hamiltonian_energy_error = -0.2152058290535308, max_hamiltonian_energy_error = -0.722900826442725, tree_depth = 4, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.9839113147639321, log_density = -114.05977455118007, hamiltonian_energy = 116.895710278864, hamiltonian_energy_error = -0.24977670356393844, max_hamiltonian_energy_error = -0.4597311887870035, tree_depth = 5, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 3, is_accept = true, acceptance_rate = 0.8428039231851684, log_density = -116.28279730858782, hamiltonian_energy = 117.37389240668534, hamiltonian_energy_error = 0.3363123847891245, max_hamiltonian_energy_error = 0.3363123847891245, tree_depth = 2, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 35, is_accept = true, acceptance_rate = 0.9908451062498774, log_density = -112.20099191405598, hamiltonian_energy = 116.87839741310675, hamiltonian_energy_error = -0.5663994155873695, max_hamiltonian_energy_error = -0.5663994155873695, tree_depth = 5, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 67, is_accept = true, acceptance_rate = 0.7825224622077603, log_density = -113.73473859506026, hamiltonian_energy = 115.91412735706594, hamiltonian_energy_error = 0.20464846045241813, max_hamiltonian_energy_error = 0.5675568649972149, tree_depth = 6, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 31, is_accept = true, acceptance_rate = 0.8481260446706104, log_density = -114.96459847065874, hamiltonian_energy = 117.52505706102288, hamiltonian_energy_error = 0.019954557341051782, max_hamiltonian_energy_error = 0.6051318760976443, tree_depth = 5, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 15, is_accept = true, acceptance_rate = 0.8985069894888915, log_density = -115.37687391618788, hamiltonian_energy = 118.35170870036292, hamiltonian_energy_error = 0.031157920544686135, max_hamiltonian_energy_error = 0.38311002943127903, tree_depth = 4, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 15, is_accept = true, acceptance_rate = 0.9752546135098723, log_density = -115.22888755533076, hamiltonian_energy = 116.66795581664802, hamiltonian_energy_error = -0.07631251952201978, max_hamiltonian_energy_error = -0.3439683031399312, tree_depth = 3, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false), (n_steps = 7, is_accept = true, acceptance_rate = 0.9999553387663411, log_density = -114.03689064828993, hamiltonian_energy = 116.18095153972405, hamiltonian_energy_error = -0.2028311777465177, max_hamiltonian_energy_error = -0.22646461229369663, tree_depth = 2, numerical_error = false, step_size = 0.7513641510268986, nom_step_size = 0.7513641510268986, is_adapt = false)])