Edit on GitHub

Linear regression

We estimate simple linear regression model with a half-T prior. First, we load the packages we use.

First, we import DynamicHMC and related libraries,

using TransformVariables, LogDensityProblems, LogDensityProblemsAD,
    DynamicHMC, TransformedLogDensities

then some packages that help code the log posterior,

using Parameters, Statistics, Random, Distributions, LinearAlgebra

then diagnostic and benchmark tools,

using MCMCDiagnosticTools, DynamicHMC.Diagnostics, BenchmarkTools

and use ForwardDiff for AD since the dimensions is small.

import ForwardDiff

Then define a structure to hold the data: observables, covariates, and the degrees of freedom for the prior.

"""
Linear regression model ``y ∼ Xβ + ϵ``, where ``ϵ ∼ N(0, σ²)`` IID.

Weakly informative prior for `β`, half-T for `σ`.
"""
struct LinearRegressionProblem{TY <: AbstractVector, TX <: AbstractMatrix,
                               Tν <: Real}
    "Observations."
    y::TY
    "Covariates"
    X::TX
    "Degrees of freedom for prior."
    ν::Tν
end
Main.LinearRegressionProblem

Then make the type callable with the parameters as a single argument.

function (problem::LinearRegressionProblem)(θ)
    @unpack y, X, ν = problem   # extract the data
    @unpack β, σ = θ            # works on the named tuple too
    ϵ_distribution = Normal(0, σ) # the error term
    ℓ_error = mapreduce((y, x) -> logpdf(ϵ_distribution, y - dot(x, β)), +,
                        y, eachrow(X))    # likelihood for error
    ℓ_σ = logpdf(TDist(ν), σ)             # prior for σ
    ℓ_β = loglikelihood(Normal(0, 10), β) # prior for β
    ℓ_error + ℓ_σ + ℓ_β
end

Make up random data and test the function runs.

N = 100
X = hcat(ones(N), randn(N, 2));
β = [1.0, 2.0, -1.0]
σ = 0.5
y = X*β .+ randn(N) .* σ;
p = LinearRegressionProblem(y, X, 1.0);
p((β = β, σ = σ))
-80.8792205436739

It is usually a good idea to benchmark and optimize your log posterior code at this stage. Above, we have carefully optimized allocations away using mapreduce.

@btime p((β = $β, σ = $σ))
-80.8792205436739

For this problem, we write a function to return the transformation (as it varies with the number of covariates).

function problem_transformation(p::LinearRegressionProblem)
    as((β = as(Array, size(p.X, 2)), σ = asℝ₊))
end
problem_transformation (generic function with 1 method)

Wrap the problem with a transformation, then use ForwardDiff for the gradient.

t = problem_transformation(p)
P = TransformedLogDensity(t, p)
∇P = ADgradient(:ForwardDiff, P);

Finally, we sample from the posterior. chain holds the chain (positions and diagnostic information), while the second returned value is the tuned sampler which would allow continuation of sampling.

results = map(_ -> mcmc_with_warmup(Random.default_rng(), ∇P, 1000), 1:5)
5-element Vector{NamedTuple{(:posterior_matrix, :tree_statistics, :κ, :ϵ), Tuple{Matrix{Float64}, Vector{DynamicHMC.TreeStatisticsNUTS}, DynamicHMC.GaussianKineticEnergy{LinearAlgebra.Diagonal{Float64, Vector{Float64}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, Float64}}}:
 (posterior_matrix = [1.096513631693949 1.005432358056836 … 0.9897438374387043 1.1205882841049597; 2.065053869125014 1.963571277538607 … 2.0004266431667967 2.020704279492226; -0.8633331732347446 -0.9615230116981112 … -0.939815856324936 -0.8764063062325042; -0.7008155487627915 -0.8781898305753181 … -0.7699622993067623 -0.7020763943161614], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-83.44052411154092, 3, turning at positions -5:2, 0.9999999999999999, 7, DynamicHMC.Directions(0xfbd5e392)), DynamicHMC.TreeStatisticsNUTS(-89.2654024641494, 3, turning at positions -1:6, 0.6583697445485608, 7, DynamicHMC.Directions(0xf2fed866)), DynamicHMC.TreeStatisticsNUTS(-82.82932344410146, 2, turning at positions -3:0, 0.9999999999999999, 3, DynamicHMC.Directions(0xa6fa3adc)), DynamicHMC.TreeStatisticsNUTS(-86.08793772899207, 3, turning at positions -2:5, 0.6687522265268054, 7, DynamicHMC.Directions(0x0d138fe5)), DynamicHMC.TreeStatisticsNUTS(-90.56291280454724, 3, turning at positions -2:5, 0.8832820205169692, 7, DynamicHMC.Directions(0x01c61dbd)), DynamicHMC.TreeStatisticsNUTS(-92.9739956196894, 2, turning at positions -1:2, 0.9659079220250376, 3, DynamicHMC.Directions(0xc880beb6)), DynamicHMC.TreeStatisticsNUTS(-90.65563880966245, 2, turning at positions 4:7, 0.9787658995559471, 7, DynamicHMC.Directions(0x00e9e9d7)), DynamicHMC.TreeStatisticsNUTS(-83.03952148679143, 3, turning at positions -4:3, 0.9972777854387125, 7, DynamicHMC.Directions(0x89d6395b)), DynamicHMC.TreeStatisticsNUTS(-82.63615134400172, 3, turning at positions -2:5, 0.9565005337539793, 7, DynamicHMC.Directions(0xc6f7f745)), DynamicHMC.TreeStatisticsNUTS(-83.01227610560737, 3, turning at positions -5:2, 0.9828511180658708, 7, DynamicHMC.Directions(0x56b8e5c2))  …  DynamicHMC.TreeStatisticsNUTS(-82.10751546811547, 2, turning at positions -4:-7, 0.9986073319086053, 7, DynamicHMC.Directions(0x5596d350)), DynamicHMC.TreeStatisticsNUTS(-82.95333372795375, 3, turning at positions -3:4, 0.8841908190137998, 7, DynamicHMC.Directions(0x290b6094)), DynamicHMC.TreeStatisticsNUTS(-83.95212841499952, 3, turning at positions -1:6, 0.8847310176169642, 7, DynamicHMC.Directions(0x3947f8c6)), DynamicHMC.TreeStatisticsNUTS(-86.24407063218368, 1, turning at positions 1:2, 0.6362531400184562, 3, DynamicHMC.Directions(0x65f1b89e)), DynamicHMC.TreeStatisticsNUTS(-80.70655881247357, 2, turning at positions 0:3, 0.9999999999999999, 3, DynamicHMC.Directions(0xc00cc9df)), DynamicHMC.TreeStatisticsNUTS(-82.87900762269834, 2, turning at positions -3:0, 0.778429740481827, 3, DynamicHMC.Directions(0x344570f4)), DynamicHMC.TreeStatisticsNUTS(-83.49027563434053, 3, turning at positions -3:4, 0.9961298312889305, 7, DynamicHMC.Directions(0x30fef4cc)), DynamicHMC.TreeStatisticsNUTS(-83.44113080566527, 3, turning at positions -2:5, 0.9656277610817152, 7, DynamicHMC.Directions(0x8c60cd75)), DynamicHMC.TreeStatisticsNUTS(-85.19508217962905, 2, turning at positions -1:2, 0.6608878960162168, 3, DynamicHMC.Directions(0xd60f6696)), DynamicHMC.TreeStatisticsNUTS(-82.62518471909202, 3, turning at positions -6:1, 0.9219181241328621, 7, DynamicHMC.Directions(0x7a4a9159))], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.05672761320829226, 0.055413416659354044, 0.04673416731864061, 0.07046874512482806], ϵ = 0.6814165465004662)
 (posterior_matrix = [0.9990449055903169 0.9699695623781616 … 1.0996206122106684 1.0056779014320816; 1.9936635746259794 1.8962906896139355 … 1.9556001489566541 2.0880749558851095; -0.9590055826332764 -0.9575250495508824 … -0.9426424607077304 -0.8948395284852693; -0.7500292831727061 -0.771014586931457 … -0.6822539197522672 -0.7394293741419053], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-82.72708668516144, 2, turning at positions -1:2, 0.9961741985515847, 3, DynamicHMC.Directions(0xf78ed7be)), DynamicHMC.TreeStatisticsNUTS(-83.30361933266195, 2, turning at positions 0:3, 0.7926129579664893, 3, DynamicHMC.Directions(0xb7add4cb)), DynamicHMC.TreeStatisticsNUTS(-82.56843579208754, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0x83ad5842)), DynamicHMC.TreeStatisticsNUTS(-82.42012604102301, 2, turning at positions 0:3, 0.7823069348567899, 3, DynamicHMC.Directions(0xec1cff63)), DynamicHMC.TreeStatisticsNUTS(-82.20595562006295, 2, turning at positions -1:2, 0.9899235984755017, 3, DynamicHMC.Directions(0xee409b02)), DynamicHMC.TreeStatisticsNUTS(-82.20048534625658, 3, turning at positions -4:3, 0.9785792001064878, 7, DynamicHMC.Directions(0x037d7a8b)), DynamicHMC.TreeStatisticsNUTS(-81.19587340224848, 2, turning at positions -3:0, 0.9397382465690766, 3, DynamicHMC.Directions(0x6877dca0)), DynamicHMC.TreeStatisticsNUTS(-81.22552895139394, 2, turning at positions -2:1, 0.8684573252714007, 3, DynamicHMC.Directions(0x4c285db1)), DynamicHMC.TreeStatisticsNUTS(-82.39132596754757, 2, turning at positions -1:2, 0.8717089678603468, 3, DynamicHMC.Directions(0x3b387542)), DynamicHMC.TreeStatisticsNUTS(-81.02246085000125, 2, turning at positions -2:1, 0.9735787066629961, 3, DynamicHMC.Directions(0xd652bcf9))  …  DynamicHMC.TreeStatisticsNUTS(-80.89084646668465, 3, turning at positions -3:4, 0.99763832736686, 7, DynamicHMC.Directions(0xd5c517ec)), DynamicHMC.TreeStatisticsNUTS(-82.57162586374807, 2, turning at positions 4:7, 0.8419769246354613, 7, DynamicHMC.Directions(0x9724e9c7)), DynamicHMC.TreeStatisticsNUTS(-81.47753766826563, 2, turning at positions -2:1, 0.9999999999999999, 3, DynamicHMC.Directions(0xab1bcf49)), DynamicHMC.TreeStatisticsNUTS(-81.72376049382743, 3, turning at positions -2:5, 0.9381999219058398, 7, DynamicHMC.Directions(0xa159e44d)), DynamicHMC.TreeStatisticsNUTS(-84.83814257852754, 2, turning at positions 0:3, 0.693656455307726, 3, DynamicHMC.Directions(0xdb1c1c13)), DynamicHMC.TreeStatisticsNUTS(-83.69973343295808, 2, turning at positions -2:1, 0.7591412311471534, 3, DynamicHMC.Directions(0xdab81b81)), DynamicHMC.TreeStatisticsNUTS(-83.57821497871487, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0xe5c58cc6)), DynamicHMC.TreeStatisticsNUTS(-83.74182574723473, 2, turning at positions -3:-6, 0.8681449025428439, 7, DynamicHMC.Directions(0xd3af3e99)), DynamicHMC.TreeStatisticsNUTS(-84.3833134130398, 2, turning at positions 0:3, 0.8270732535047808, 3, DynamicHMC.Directions(0x5978b2ab)), DynamicHMC.TreeStatisticsNUTS(-84.13737412666471, 2, turning at positions -3:0, 0.8859584844362871, 3, DynamicHMC.Directions(0x818d4490))], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.05167172601472521, 0.04875284909476561, 0.04919857919006001, 0.0709088444803029], ϵ = 0.8440297318165975)
 (posterior_matrix = [1.0292932705651583 1.017730768751647 … 1.02482922735323 1.007609230850035; 2.0264294216762795 1.9783038515066032 … 1.9534064673461464 1.9921411876406772; -0.948348145730263 -0.8868759532884808 … -0.9509428637791082 -0.9465863920513283; -0.8015655721153506 -0.686605010751253 … -0.6663604162418404 -0.8635942605521775], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-83.08872763158323, 2, turning at positions -3:0, 0.8440633616762518, 3, DynamicHMC.Directions(0x369d395c)), DynamicHMC.TreeStatisticsNUTS(-82.0920718656376, 3, turning at positions -7:0, 0.9550975160038694, 7, DynamicHMC.Directions(0x318c8340)), DynamicHMC.TreeStatisticsNUTS(-81.34866241486283, 3, turning at positions -4:3, 0.9674369938737396, 7, DynamicHMC.Directions(0xa1fc4a5b)), DynamicHMC.TreeStatisticsNUTS(-84.16030932304663, 2, turning at positions 0:3, 0.7275120806783195, 3, DynamicHMC.Directions(0xd44e7f0f)), DynamicHMC.TreeStatisticsNUTS(-86.29635664513819, 3, turning at positions -6:1, 0.9682320675178453, 7, DynamicHMC.Directions(0xd1016071)), DynamicHMC.TreeStatisticsNUTS(-84.52594297990017, 2, turning at positions -1:2, 0.9585629214179372, 3, DynamicHMC.Directions(0x856eafea)), DynamicHMC.TreeStatisticsNUTS(-83.77425900478148, 3, turning at positions -7:0, 0.9543713911654229, 7, DynamicHMC.Directions(0xb2032f78)), DynamicHMC.TreeStatisticsNUTS(-82.64018338816598, 2, turning at positions -3:0, 0.9781880844815559, 3, DynamicHMC.Directions(0x61890350)), DynamicHMC.TreeStatisticsNUTS(-83.09796669502298, 2, turning at positions -2:1, 0.8808630720566462, 3, DynamicHMC.Directions(0x59e825ad)), DynamicHMC.TreeStatisticsNUTS(-85.28855799765203, 2, turning at positions -3:-6, 0.9328630747572337, 7, DynamicHMC.Directions(0x45baa939))  …  DynamicHMC.TreeStatisticsNUTS(-83.66943807989209, 2, turning at positions -2:1, 0.9703142657539656, 3, DynamicHMC.Directions(0x1b288955)), DynamicHMC.TreeStatisticsNUTS(-84.18981571866482, 3, turning at positions 0:7, 0.9338529260005338, 7, DynamicHMC.Directions(0x833d821f)), DynamicHMC.TreeStatisticsNUTS(-84.82484025415546, 3, turning at positions -4:3, 0.977964290533226, 7, DynamicHMC.Directions(0x51895b4b)), DynamicHMC.TreeStatisticsNUTS(-83.0603509925001, 2, turning at positions 3:6, 0.9961771120927043, 7, DynamicHMC.Directions(0x7932c23e)), DynamicHMC.TreeStatisticsNUTS(-81.91445704314717, 3, turning at positions -5:2, 0.9954570641263516, 7, DynamicHMC.Directions(0x0676a04a)), DynamicHMC.TreeStatisticsNUTS(-85.91361843438726, 2, turning at positions 0:3, 0.7045878504660862, 3, DynamicHMC.Directions(0xb1731633)), DynamicHMC.TreeStatisticsNUTS(-85.44800343062855, 2, turning at positions 0:3, 0.7792024093377744, 3, DynamicHMC.Directions(0xe16a5efb)), DynamicHMC.TreeStatisticsNUTS(-83.88504910294624, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0x8f65faca)), DynamicHMC.TreeStatisticsNUTS(-82.49225782353258, 2, turning at positions -3:0, 0.9842163196862829, 3, DynamicHMC.Directions(0xbbc517c8)), DynamicHMC.TreeStatisticsNUTS(-84.56706758215529, 2, turning at positions -3:0, 0.7991273159698001, 3, DynamicHMC.Directions(0xc1e83b50))], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.0498848342840449, 0.055275739592553705, 0.04555527409987744, 0.06538395780435943], ϵ = 0.7775089074693241)
 (posterior_matrix = [1.0376101101279498 1.0494493556765083 … 1.0360028708885867 1.1139557641389424; 1.966713077395872 2.0461393507988244 … 2.0224101133517824 1.985187078592518; -0.9077413011637416 -0.9372064135850116 … -0.941273413865512 -0.9286698063577096; -0.5744075268619704 -0.8102500938005742 … -0.707348153157809 -0.7108624173934107], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-83.17314647196136, 3, turning at positions -2:5, 0.997115211750315, 7, DynamicHMC.Directions(0xa16f4955)), DynamicHMC.TreeStatisticsNUTS(-82.96475602550021, 3, turning at positions -1:6, 0.9999999999999999, 7, DynamicHMC.Directions(0x66d28256)), DynamicHMC.TreeStatisticsNUTS(-82.77494378939105, 3, turning at positions -6:1, 0.9807040507614033, 7, DynamicHMC.Directions(0xd9434e91)), DynamicHMC.TreeStatisticsNUTS(-83.27948623153928, 3, turning at positions -7:0, 0.9376233990071174, 7, DynamicHMC.Directions(0xe96b1b80)), DynamicHMC.TreeStatisticsNUTS(-82.0749947248405, 2, turning at positions 0:3, 0.9999999999999999, 3, DynamicHMC.Directions(0x4af862d7)), DynamicHMC.TreeStatisticsNUTS(-81.42268027461435, 3, turning at positions -4:3, 0.9710265466284077, 7, DynamicHMC.Directions(0x8795b7e3)), DynamicHMC.TreeStatisticsNUTS(-82.87671515503212, 3, turning at positions -5:2, 0.8085991186233129, 7, DynamicHMC.Directions(0xdefe742a)), DynamicHMC.TreeStatisticsNUTS(-84.40016247937288, 2, turning at positions -1:2, 0.807917177008465, 3, DynamicHMC.Directions(0xe84b352a)), DynamicHMC.TreeStatisticsNUTS(-85.51191607960763, 3, turning at positions -1:6, 0.9851170424534581, 7, DynamicHMC.Directions(0x6cda9e7e)), DynamicHMC.TreeStatisticsNUTS(-82.81608831584043, 3, turning at positions -4:3, 0.9956544559744962, 7, DynamicHMC.Directions(0xb60b0e63))  …  DynamicHMC.TreeStatisticsNUTS(-86.02817157207421, 3, turning at positions -3:4, 0.9999999999999999, 7, DynamicHMC.Directions(0x555ce6bc)), DynamicHMC.TreeStatisticsNUTS(-83.46859303697599, 3, turning at positions 0:7, 0.9999999999999999, 7, DynamicHMC.Directions(0xbaeda647)), DynamicHMC.TreeStatisticsNUTS(-84.02204197474796, 2, turning at positions 0:3, 0.8573179417900482, 3, DynamicHMC.Directions(0xd952c483)), DynamicHMC.TreeStatisticsNUTS(-83.31635171678433, 2, turning at positions 2:5, 0.9415061546963418, 7, DynamicHMC.Directions(0xc6662fc5)), DynamicHMC.TreeStatisticsNUTS(-84.19104294667912, 3, turning at positions -7:0, 0.9614458180701939, 7, DynamicHMC.Directions(0x731ac990)), DynamicHMC.TreeStatisticsNUTS(-83.96611780907529, 3, turning at positions -2:5, 0.9434910045911477, 7, DynamicHMC.Directions(0xa13234c5)), DynamicHMC.TreeStatisticsNUTS(-85.54705920809789, 2, turning at positions 0:3, 0.8688290169701157, 3, DynamicHMC.Directions(0x375b3c9b)), DynamicHMC.TreeStatisticsNUTS(-82.09453886190158, 2, turning at positions -3:0, 0.9848540491057381, 3, DynamicHMC.Directions(0x76285d2c)), DynamicHMC.TreeStatisticsNUTS(-81.63960793279392, 2, turning at positions 0:3, 0.989537942204366, 3, DynamicHMC.Directions(0x3e9dc383)), DynamicHMC.TreeStatisticsNUTS(-82.86293456632164, 2, turning at positions 3:6, 0.8092997824997257, 7, DynamicHMC.Directions(0x44039406))], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.052263667781537376, 0.05342935440821451, 0.04586240677339016, 0.07401538642017806], ϵ = 0.7385699558173652)
 (posterior_matrix = [1.0060287393626166 1.0053706352077223 … 1.1457983363458275 1.123533907025646; 1.9922848794943462 2.067157958608109 … 2.0065811628598933 2.0069572725776132; -0.948550124324446 -0.9810222998674596 … -0.937211277710722 -0.9049661071455365; -0.7152774168304687 -0.6601909406036002 … -0.6151116485632289 -0.7307225543404727], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-81.81138471647641, 2, turning at positions 0:3, 0.9999999999999999, 3, DynamicHMC.Directions(0x2e02e823)), DynamicHMC.TreeStatisticsNUTS(-82.70418416424995, 2, turning at positions 0:3, 0.80252568277739, 3, DynamicHMC.Directions(0x7fc9f65f)), DynamicHMC.TreeStatisticsNUTS(-86.47971423276407, 3, turning at positions -2:5, 0.8276713926951008, 7, DynamicHMC.Directions(0xb43ca0f5)), DynamicHMC.TreeStatisticsNUTS(-82.81160118253366, 3, turning at positions -6:1, 0.9596643836179293, 7, DynamicHMC.Directions(0x0af08c09)), DynamicHMC.TreeStatisticsNUTS(-82.10384293733144, 3, turning at positions -1:6, 0.9705315161467627, 7, DynamicHMC.Directions(0xe7194b7e)), DynamicHMC.TreeStatisticsNUTS(-85.23679374937416, 2, turning at positions -3:0, 0.793626928769017, 3, DynamicHMC.Directions(0x0545cc44)), DynamicHMC.TreeStatisticsNUTS(-83.2518046193763, 3, turning at positions -6:1, 0.9988949657690342, 7, DynamicHMC.Directions(0xe3fca0b1)), DynamicHMC.TreeStatisticsNUTS(-82.35915866551566, 2, turning at positions 0:3, 0.991890420142556, 3, DynamicHMC.Directions(0x4e045b1f)), DynamicHMC.TreeStatisticsNUTS(-82.9575964811603, 3, turning at positions -2:5, 0.8767604043071426, 7, DynamicHMC.Directions(0x60db1255)), DynamicHMC.TreeStatisticsNUTS(-83.06107751077523, 3, turning at positions -7:0, 0.9550457349805519, 7, DynamicHMC.Directions(0x5a5fe540))  …  DynamicHMC.TreeStatisticsNUTS(-85.50328563497345, 2, turning at positions 4:7, 0.9675966912270433, 7, DynamicHMC.Directions(0x67f8dc17)), DynamicHMC.TreeStatisticsNUTS(-88.97568197687548, 3, turning at positions -3:4, 0.6096998601872536, 7, DynamicHMC.Directions(0x9c16ab2c)), DynamicHMC.TreeStatisticsNUTS(-87.72553760007887, 2, turning at positions -1:2, 0.9845542876993355, 3, DynamicHMC.Directions(0x90f061d6)), DynamicHMC.TreeStatisticsNUTS(-87.60053910785531, 3, turning at positions -6:1, 0.7546601905244638, 7, DynamicHMC.Directions(0x98cf4db1)), DynamicHMC.TreeStatisticsNUTS(-86.46277325044548, 3, turning at positions -2:5, 0.9999999999999999, 7, DynamicHMC.Directions(0x3ea6712d)), DynamicHMC.TreeStatisticsNUTS(-85.09272772677495, 3, turning at positions -1:6, 0.9999999999999999, 7, DynamicHMC.Directions(0xd0f25dae)), DynamicHMC.TreeStatisticsNUTS(-84.95225675266333, 3, turning at positions -6:1, 0.97421943314652, 7, DynamicHMC.Directions(0x7311b511)), DynamicHMC.TreeStatisticsNUTS(-83.50221084495841, 3, turning at positions 0:7, 0.9522831769445365, 7, DynamicHMC.Directions(0xc759d16f)), DynamicHMC.TreeStatisticsNUTS(-84.58700155179103, 2, turning at positions -3:0, 0.7754689311517867, 3, DynamicHMC.Directions(0x7c2d4fe8)), DynamicHMC.TreeStatisticsNUTS(-84.95269109397431, 2, turning at positions -2:-5, 0.9367005146787176, 7, DynamicHMC.Directions(0xc93293e2))], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.04570681741743813, 0.05350375475305488, 0.04796949273623679, 0.07387906259282682], ϵ = 0.748069911936487)

We use the transformation to obtain the posterior from the chain.

posterior = transform.(t, eachcol(pool_posterior_matrices(results)));

Extract the parameter posterior means: β,

posterior_β = mean(first, posterior)
3-element Vector{Float64}:
  1.032862334981682
  1.9798614581093363
 -0.9243007940491812

then σ:

posterior_σ = mean(last, posterior)
0.4908138055588217

Effective sample sizes (of untransformed draws)

ess, R̂ = ess_rhat(stack_posterior_matrices(results))
(ess = [5159.6081335923445, 4894.600925259005, 4607.091718331742, 4411.055483737393], rhat = [1.0000879793590638, 1.0016250901879526, 1.0003141971578404, 1.000224526873462])

summarize NUTS-specific statistics of all chains

summarize_tree_statistics(mapreduce(x -> x.tree_statistics, vcat, results))
Hamiltonian Monte Carlo sample of length 5000
  acceptance rate mean: 0.9, 5/25/50/75/95%: 0.67 0.85 0.94 0.99 1.0
  termination: divergence => 0%, max_depth => 0%, turning => 100%
  depth: 0 => 0%, 1 => 2%, 2 => 58%, 3 => 41%

This page was generated using Literate.jl.