GitHub

Logistic regression

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, StatsFuns, LogExpFunctions

then diagnostic and benchmark tools,

using MCMCDiagnosticTools, BenchmarkTools

and use ForwardDiff for AD since the dimensions is small.

import ForwardDiff

"""
Logistic regression.

For each draw, ``logit(Pr(yᵢ == 1)) ∼ Xᵢ β``. Uses a `β ∼ Normal(0, σ)` prior.

`X` is supposed to include the `1`s for the intercept.
"""
struct LogisticRegression{Ty, TX, Tσ}
    y::Ty
    X::TX
    σ::Tσ
end

function (problem::LogisticRegression)(θ)
    @unpack y, X, σ = problem
    @unpack β = θ
    ℓ_y = mapreduce((y, x) -> logpdf(Bernoulli(logistic(dot(x, β))), y), +, y, eachrow(X))
    ℓ_β =  loglikelihood(Normal(0, σ), β)
    ℓ_y + ℓ_β
end

Make up parameters, generate data using random draws.

N = 1000
β = [1.0, 2.0]
X = hcat(ones(N), randn(N))
y = rand.(Bernoulli.(logistic.(X*β)));

Create a problem, apply a transformation, then use automatic differentiation.

p = LogisticRegression(y, X, 10.0)   # data and (vague) priors
t = as((β = as(Array, length(β)), )) # identity transformation, just to get the dimension
P = TransformedLogDensity(t, p)      # transformed
∇P = ADgradient(:ForwardDiff, P)
ForwardDiff AD wrapper for TransformedLogDensity of dimension 2, w/ chunk size 2

Benchmark

@btime p((β = $β,))
-439.4945572222126

Sample using NUTS, random starting point.

results = map(_ -> mcmc_with_warmup(Random.default_rng(), ∇P, 1000), 1:5)
5-element Vector{@NamedTuple{posterior_matrix::Matrix{Float64}, tree_statistics::Vector{DynamicHMC.TreeStatisticsNUTS}, logdensities::Vector{Float64}, κ::DynamicHMC.GaussianKineticEnergy{LinearAlgebra.Diagonal{Float64, Vector{Float64}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, ϵ::Float64}}:
 (posterior_matrix = [0.87944345589263 0.8647223862783247 … 1.0662082730492644 0.8166771393744177; 2.1792645205916976 2.048174003240724 … 1.9222289145384048 2.0508585725935564], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-439.5720826933488, 1, turning at positions -1:0, 1.0, 1, DynamicHMC.Directions(0xbfd63614)), DynamicHMC.TreeStatisticsNUTS(-439.1280421297494, 2, turning at positions 0:3, 0.9930816234743974, 3, DynamicHMC.Directions(0xd5235b47)), DynamicHMC.TreeStatisticsNUTS(-440.0479515590839, 2, turning at positions 0:3, 0.8498859475662591, 3, DynamicHMC.Directions(0xe5b4eecb)), DynamicHMC.TreeStatisticsNUTS(-439.9399502771731, 2, turning at positions -2:1, 0.9999999999999999, 3, DynamicHMC.Directions(0x36b7fc39)), DynamicHMC.TreeStatisticsNUTS(-438.9757235250866, 2, turning at positions -1:2, 0.9742710719212725, 3, DynamicHMC.Directions(0xfa4ff90e)), DynamicHMC.TreeStatisticsNUTS(-439.51240870213695, 2, turning at positions 0:3, 0.9105688014208008, 3, DynamicHMC.Directions(0x9e3e049f)), DynamicHMC.TreeStatisticsNUTS(-441.1890317011903, 1, turning at positions 1:2, 0.8123941151238045, 3, DynamicHMC.Directions(0x6c992f72)), DynamicHMC.TreeStatisticsNUTS(-439.52361878747143, 2, turning at positions 0:3, 0.9780312462827458, 3, DynamicHMC.Directions(0xe7436ca7)), DynamicHMC.TreeStatisticsNUTS(-438.83299822066743, 2, turning at positions 0:3, 0.9643317811521065, 3, DynamicHMC.Directions(0x11a34677)), DynamicHMC.TreeStatisticsNUTS(-439.0545297318966, 2, turning at positions -3:-6, 0.9347463578779759, 7, DynamicHMC.Directions(0xbf9c5e11))  …  DynamicHMC.TreeStatisticsNUTS(-439.30016212795493, 2, turning at positions 0:3, 0.9332706879258764, 3, DynamicHMC.Directions(0x5726f283)), DynamicHMC.TreeStatisticsNUTS(-439.22971593611237, 1, turning at positions 2:3, 0.9539591429735544, 3, DynamicHMC.Directions(0x24f07bc3)), DynamicHMC.TreeStatisticsNUTS(-443.684173568597, 2, turning at positions -3:0, 0.5253326087266174, 3, DynamicHMC.Directions(0xbef14c8c)), DynamicHMC.TreeStatisticsNUTS(-440.5588899558829, 1, turning at positions 0:1, 0.9558415105458695, 1, DynamicHMC.Directions(0x583503c5)), DynamicHMC.TreeStatisticsNUTS(-440.4715988045996, 2, turning at positions 0:3, 0.9999999999999999, 3, DynamicHMC.Directions(0xf14fea57)), DynamicHMC.TreeStatisticsNUTS(-439.43001988764416, 1, turning at positions -1:0, 1.0, 1, DynamicHMC.Directions(0x10f0762a)), DynamicHMC.TreeStatisticsNUTS(-440.03228413469014, 2, turning at positions -3:0, 0.9695449591394388, 3, DynamicHMC.Directions(0xbc337434)), DynamicHMC.TreeStatisticsNUTS(-439.181626275365, 2, turning at positions -5:-6, 0.9958683394565677, 7, DynamicHMC.Directions(0x50121de1)), DynamicHMC.TreeStatisticsNUTS(-443.05580439725816, 2, turning at positions -1:2, 0.5235760334819489, 3, DynamicHMC.Directions(0xe66d2afe)), DynamicHMC.TreeStatisticsNUTS(-441.5709972366418, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0xbf47e112))], logdensities = [-438.96437538844685, -438.47792646763844, -439.7621459591536, -438.8044333905295, -438.47046075653907, -439.27230214553146, -439.31812332071286, -438.46614801020263, -438.5250677424096, -438.66174082942575  …  -438.4923492096084, -438.6725697343373, -440.1999668010675, -440.3747031552267, -439.3943447304377, -439.2712458053219, -439.04249934218933, -438.9655333309843, -442.02040049349665, -438.8582790209976], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.10269349012060792, 0.13817477050665383], ϵ = 0.6795076666825658)
 (posterior_matrix = [0.8640682431611753 0.7378876889878866 … 1.0305722456016844 1.0233339873339644; 2.0668478540781345 1.9265729043857363 … 2.2387136753088197 2.2821121546133116], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-438.74243422046544, 1, turning at positions -1:-2, 0.9666206657811003, 3, DynamicHMC.Directions(0x7d3ddff9)), DynamicHMC.TreeStatisticsNUTS(-440.65791818256423, 2, turning at positions -1:2, 0.7465777082347822, 3, DynamicHMC.Directions(0x2c7ece56)), DynamicHMC.TreeStatisticsNUTS(-440.5513906640512, 2, turning at positions -2:1, 0.9552064675212524, 3, DynamicHMC.Directions(0xe9b28d4d)), DynamicHMC.TreeStatisticsNUTS(-441.0290929652785, 2, turning at positions -1:-2, 0.6932236621142238, 5, DynamicHMC.Directions(0x14133033)), DynamicHMC.TreeStatisticsNUTS(-440.0083646888829, 2, turning at positions 0:3, 0.9999999999999999, 3, DynamicHMC.Directions(0x399b857b)), DynamicHMC.TreeStatisticsNUTS(-439.45138549641246, 2, turning at positions -2:1, 0.9999999999999999, 3, DynamicHMC.Directions(0x6b681be5)), DynamicHMC.TreeStatisticsNUTS(-438.9952321703886, 3, turning at positions 0:7, 0.9962758518440101, 7, DynamicHMC.Directions(0x739c8b0f)), DynamicHMC.TreeStatisticsNUTS(-439.0942697062829, 2, turning at positions -2:1, 0.956290016705421, 3, DynamicHMC.Directions(0x3a5bf82d)), DynamicHMC.TreeStatisticsNUTS(-439.6620114649722, 2, turning at positions 1:4, 0.9426746903765401, 7, DynamicHMC.Directions(0xe5de1c14)), DynamicHMC.TreeStatisticsNUTS(-442.69740509801596, 2, turning at positions -1:2, 0.5793973720560807, 3, DynamicHMC.Directions(0x9cc4b406))  …  DynamicHMC.TreeStatisticsNUTS(-442.3260494463705, 1, turning at positions 1:2, 0.6893941356687696, 3, DynamicHMC.Directions(0x44c771ba)), DynamicHMC.TreeStatisticsNUTS(-441.3530788325306, 2, turning at positions -2:1, 0.9999999999999999, 3, DynamicHMC.Directions(0xb2d2ab51)), DynamicHMC.TreeStatisticsNUTS(-440.4891564394362, 2, turning at positions 1:4, 0.9167130212025149, 7, DynamicHMC.Directions(0xfb9a3ab4)), DynamicHMC.TreeStatisticsNUTS(-439.643175175741, 2, turning at positions 3:6, 0.8649399022326769, 7, DynamicHMC.Directions(0x7e64c8ee)), DynamicHMC.TreeStatisticsNUTS(-439.22810906177943, 2, turning at positions -3:0, 0.9964071133513469, 3, DynamicHMC.Directions(0x48e5e1bc)), DynamicHMC.TreeStatisticsNUTS(-439.40364612533415, 2, turning at positions -2:1, 0.9980946023550964, 3, DynamicHMC.Directions(0xc19ff1bd)), DynamicHMC.TreeStatisticsNUTS(-438.8805825496269, 2, turning at positions -1:2, 0.9907682827763925, 3, DynamicHMC.Directions(0xdf3bd2e6)), DynamicHMC.TreeStatisticsNUTS(-439.7872073327701, 1, turning at positions -2:-3, 0.8967417843884345, 3, DynamicHMC.Directions(0x604e0298)), DynamicHMC.TreeStatisticsNUTS(-440.01308325096124, 3, turning at positions -5:2, 0.9328826240923458, 7, DynamicHMC.Directions(0x247ec1c2)), DynamicHMC.TreeStatisticsNUTS(-440.0144407488083, 1, turning at positions -1:0, 0.976999334372214, 1, DynamicHMC.Directions(0x48bcc790))], logdensities = [-438.4997574044444, -440.0868516591243, -438.94828715420743, -439.72718041579316, -439.4483142261621, -438.82423177930536, -438.7723447765679, -438.914869021018, -439.4665095825823, -439.4665095825823  …  -441.8027141074687, -439.81498094777106, -438.81943049709184, -438.88760045956377, -439.0713529338895, -438.8733712337936, -438.4817715129646, -438.7853632724955, -439.6681227629271, -439.92589068207246], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.09642450264393365, 0.1298589244570766], ϵ = 0.8344144495979853)
 (posterior_matrix = [0.9445151247804454 0.9752192081349853 … 0.8866411781618694 0.8732884007314617; 2.0537684729784296 1.8221381361259725 … 1.864201088085394 2.197745639570347], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-439.5299905013875, 2, turning at positions -1:2, 0.9793438635962252, 3, DynamicHMC.Directions(0x8ce05712)), DynamicHMC.TreeStatisticsNUTS(-442.593047629159, 2, turning at positions -1:2, 0.5337389039080945, 3, DynamicHMC.Directions(0x826de0fe)), DynamicHMC.TreeStatisticsNUTS(-441.3386938414052, 1, turning at positions -1:0, 1.0, 1, DynamicHMC.Directions(0x2ac128ee)), DynamicHMC.TreeStatisticsNUTS(-440.28275520823576, 3, turning at positions -1:6, 0.9874619594057025, 7, DynamicHMC.Directions(0xb16b7286)), DynamicHMC.TreeStatisticsNUTS(-438.84843476036957, 3, turning at positions 0:7, 0.9727019897444551, 7, DynamicHMC.Directions(0xc7e51817)), DynamicHMC.TreeStatisticsNUTS(-438.83916364659024, 1, turning at positions -1:0, 0.9727530414966901, 1, DynamicHMC.Directions(0x48319506)), DynamicHMC.TreeStatisticsNUTS(-439.2546691886395, 1, turning at positions -2:-3, 0.9817928926106814, 3, DynamicHMC.Directions(0x7b397e40)), DynamicHMC.TreeStatisticsNUTS(-439.45639860222485, 2, turning at positions -2:1, 0.867216273139534, 3, DynamicHMC.Directions(0x2c155e19)), DynamicHMC.TreeStatisticsNUTS(-439.10891210153585, 3, turning at positions -3:4, 0.9999999999999999, 7, DynamicHMC.Directions(0x90556adc)), DynamicHMC.TreeStatisticsNUTS(-439.592508123346, 2, turning at positions -1:2, 0.8496853363731455, 3, DynamicHMC.Directions(0xad68e192))  …  DynamicHMC.TreeStatisticsNUTS(-440.3103283860569, 2, turning at positions -1:2, 0.9888617854938744, 3, DynamicHMC.Directions(0x35552bca)), DynamicHMC.TreeStatisticsNUTS(-440.15306474201543, 1, turning at positions 0:1, 1.0, 1, DynamicHMC.Directions(0xa861328f)), DynamicHMC.TreeStatisticsNUTS(-439.80952734608593, 2, turning at positions 1:4, 0.9999999999999999, 7, DynamicHMC.Directions(0xfc188054)), DynamicHMC.TreeStatisticsNUTS(-439.2988219256166, 2, turning at positions -1:2, 0.9786658237272771, 3, DynamicHMC.Directions(0x82db81c6)), DynamicHMC.TreeStatisticsNUTS(-440.69721583799566, 2, turning at positions 0:3, 0.7298085347316716, 3, DynamicHMC.Directions(0xe11cc1e3)), DynamicHMC.TreeStatisticsNUTS(-439.10114054510353, 2, turning at positions -2:1, 0.994404858533286, 3, DynamicHMC.Directions(0x7815d281)), DynamicHMC.TreeStatisticsNUTS(-440.354381572245, 1, turning at positions -1:-2, 0.8087841974709847, 3, DynamicHMC.Directions(0x0289813d)), DynamicHMC.TreeStatisticsNUTS(-439.8056276453669, 2, turning at positions 0:3, 0.9811812679990076, 3, DynamicHMC.Directions(0x9ce2233f)), DynamicHMC.TreeStatisticsNUTS(-439.9000582594927, 1, turning at positions -1:0, 0.9816424460519466, 1, DynamicHMC.Directions(0x2751b316)), DynamicHMC.TreeStatisticsNUTS(-440.3429108878749, 3, turning at positions -3:4, 0.9852768239241428, 7, DynamicHMC.Directions(0xc118bdf4))], logdensities = [-438.56052953214476, -441.6314663836549, -439.88459018118016, -438.47639077590105, -438.6716675686295, -438.79975316040066, -438.5042256680139, -439.0471157002816, -438.6631440477166, -439.0950480675291  …  -440.0142538738821, -439.7794702526984, -439.13618588394866, -438.5014635871413, -438.6849719028434, -438.9338461419015, -439.1896074825734, -439.5810917941901, -439.7018922241374, -439.1780192196609], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.09134309451958587, 0.12401005922373814], ϵ = 0.7867914533841917)
 (posterior_matrix = [0.8770533030311569 0.8944738332748822 … 1.0142682102445264 1.0341461883672682; 2.0508541196420143 1.976462937304785 … 2.170280103678981 2.202184408528402], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-438.56637950897965, 2, turning at positions 2:5, 0.9921164393904939, 7, DynamicHMC.Directions(0x383ed015)), DynamicHMC.TreeStatisticsNUTS(-438.8591354663403, 1, turning at positions -2:-3, 0.959345491617852, 3, DynamicHMC.Directions(0x79b9c424)), DynamicHMC.TreeStatisticsNUTS(-439.388825352863, 2, turning at positions -3:0, 0.9352039450618067, 3, DynamicHMC.Directions(0xb8fb7a90)), DynamicHMC.TreeStatisticsNUTS(-439.79981237309477, 2, turning at positions 3:4, 0.9899269944205991, 7, DynamicHMC.Directions(0xc5bf5b14)), DynamicHMC.TreeStatisticsNUTS(-439.39465150751556, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0x00a469e2)), DynamicHMC.TreeStatisticsNUTS(-439.3019949256806, 2, turning at positions -2:-3, 0.9637046795753706, 5, DynamicHMC.Directions(0x19453212)), DynamicHMC.TreeStatisticsNUTS(-439.58696422936725, 1, turning at positions -1:-2, 0.9668388172665862, 3, DynamicHMC.Directions(0x2fb893f1)), DynamicHMC.TreeStatisticsNUTS(-439.0844809960064, 2, turning at positions -1:-4, 0.9632308183639054, 7, DynamicHMC.Directions(0x0246e22b)), DynamicHMC.TreeStatisticsNUTS(-439.9243232317953, 3, turning at positions -4:3, 0.9238029926246354, 7, DynamicHMC.Directions(0xc93803ab)), DynamicHMC.TreeStatisticsNUTS(-439.2023300609611, 2, turning at positions -1:2, 0.976307727032979, 3, DynamicHMC.Directions(0x22614b86))  …  DynamicHMC.TreeStatisticsNUTS(-443.0252366367844, 1, turning at positions -1:-2, 0.7378571860830658, 3, DynamicHMC.Directions(0x66d743a5)), DynamicHMC.TreeStatisticsNUTS(-443.3946477573647, 2, turning at positions 2:5, 0.7084149394759124, 7, DynamicHMC.Directions(0x67fcc52d)), DynamicHMC.TreeStatisticsNUTS(-440.6767460227307, 3, turning at positions 0:7, 0.9705428511721131, 7, DynamicHMC.Directions(0x05da1497)), DynamicHMC.TreeStatisticsNUTS(-442.0113956269375, 2, turning at positions -3:0, 0.8070759910239489, 3, DynamicHMC.Directions(0xebd7da8c)), DynamicHMC.TreeStatisticsNUTS(-440.06732067224084, 3, turning at positions -3:4, 0.9999999999999999, 7, DynamicHMC.Directions(0xa287ede4)), DynamicHMC.TreeStatisticsNUTS(-440.23850373748155, 2, turning at positions -4:-5, 0.9809166137985825, 5, DynamicHMC.Directions(0x03be9e70)), DynamicHMC.TreeStatisticsNUTS(-439.96681307019907, 2, turning at positions 4:5, 0.9640642578842529, 5, DynamicHMC.Directions(0x23011b97)), DynamicHMC.TreeStatisticsNUTS(-439.6159657086531, 2, turning at positions -3:0, 0.9441508217570408, 3, DynamicHMC.Directions(0xf20ad118)), DynamicHMC.TreeStatisticsNUTS(-439.54797139086014, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0xbc06af42)), DynamicHMC.TreeStatisticsNUTS(-439.55922932987687, 2, turning at positions -2:1, 0.9927477330330753, 3, DynamicHMC.Directions(0x07e14561))], logdensities = [-438.43439950081944, -438.64265688870984, -439.07897384098726, -439.4492301472862, -438.45112408619286, -438.98709105936433, -438.44551733863767, -438.9467167982689, -438.8646605614268, -438.7973678645117  …  -439.3167498364597, -439.3167498364597, -439.58555430000445, -439.69823918777274, -439.7166061749489, -439.0048214189683, -438.76082886285394, -439.5925616457904, -439.19692501497707, -439.53286356702995], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.09408912042114634, 0.13005077526225964], ϵ = 0.6556617844181579)
 (posterior_matrix = [0.8054780999680442 0.7389910250263252 … 0.7516762857269293 0.9709467729410717; 1.8763113953813577 2.0070258793297104 … 1.8507775796630035 2.305077964934351], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-440.7877788517282, 2, turning at positions 3:6, 0.93871416682262, 7, DynamicHMC.Directions(0x435d04de)), DynamicHMC.TreeStatisticsNUTS(-440.49124669426214, 2, turning at positions -3:0, 0.95875843621694, 3, DynamicHMC.Directions(0xeab6c6e0)), DynamicHMC.TreeStatisticsNUTS(-442.0914296958468, 2, turning at positions -2:1, 0.9515263233804129, 3, DynamicHMC.Directions(0x000af9b1)), DynamicHMC.TreeStatisticsNUTS(-441.7884922802155, 2, turning at positions 0:3, 0.9791758118996051, 3, DynamicHMC.Directions(0x8f0a5387)), DynamicHMC.TreeStatisticsNUTS(-440.8701339006456, 2, turning at positions -1:2, 0.9947009075869913, 3, DynamicHMC.Directions(0xab2326a2)), DynamicHMC.TreeStatisticsNUTS(-439.53317556010944, 2, turning at positions -4:-7, 0.9627866155067897, 7, DynamicHMC.Directions(0x6833b640)), DynamicHMC.TreeStatisticsNUTS(-441.8778177989409, 2, turning at positions 0:3, 0.7305299331493784, 3, DynamicHMC.Directions(0x39c4183b)), DynamicHMC.TreeStatisticsNUTS(-440.45885953385243, 1, turning at positions -1:-2, 0.7415469305431422, 3, DynamicHMC.Directions(0x1d249bd9)), DynamicHMC.TreeStatisticsNUTS(-443.4877340216282, 3, turning at positions -1:6, 0.6673486348657259, 7, DynamicHMC.Directions(0xe52274ae)), DynamicHMC.TreeStatisticsNUTS(-442.10537374159026, 3, turning at positions -7:0, 0.9655062532150518, 7, DynamicHMC.Directions(0xa9e51de8))  …  DynamicHMC.TreeStatisticsNUTS(-438.76312837460256, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0x2ec2e46e)), DynamicHMC.TreeStatisticsNUTS(-438.9668474741652, 2, turning at positions 0:3, 0.9811657464208879, 3, DynamicHMC.Directions(0x2bf362cf)), DynamicHMC.TreeStatisticsNUTS(-439.30929676137947, 2, turning at positions -1:2, 0.965471941389168, 3, DynamicHMC.Directions(0xf5af9a86)), DynamicHMC.TreeStatisticsNUTS(-441.724050097528, 2, turning at positions -1:2, 0.7308331866409911, 3, DynamicHMC.Directions(0xa9572376)), DynamicHMC.TreeStatisticsNUTS(-440.63754137228693, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0xdb8b05ae)), DynamicHMC.TreeStatisticsNUTS(-440.77249712753223, 2, turning at positions -3:0, 0.9178845759856107, 3, DynamicHMC.Directions(0xdbc22948)), DynamicHMC.TreeStatisticsNUTS(-439.13021720383904, 2, turning at positions -2:1, 0.9853354033439574, 3, DynamicHMC.Directions(0x97f674b5)), DynamicHMC.TreeStatisticsNUTS(-439.0541275832939, 3, turning at positions -7:0, 0.9999999999999999, 7, DynamicHMC.Directions(0x1c287908)), DynamicHMC.TreeStatisticsNUTS(-441.5506189239825, 2, turning at positions 1:4, 0.8860505500104517, 7, DynamicHMC.Directions(0xea50d614)), DynamicHMC.TreeStatisticsNUTS(-442.30985164733494, 2, turning at positions 1:4, 0.9008561121158785, 7, DynamicHMC.Directions(0xb3943964))], logdensities = [-439.5760174942234, -440.0097354842661, -441.125088042166, -440.3665599675572, -438.76310869130725, -438.984493382159, -438.47686827881046, -438.47686827881046, -440.96864943285027, -439.878576596822  …  -438.65884605113746, -438.87860470565374, -438.97627263312495, -440.62144729690056, -439.61846860377824, -438.8491061298699, -439.12017118530423, -438.72392714501655, -440.36082181816926, -439.95222259586814], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.10147942769212548, 0.1547542005668675], ϵ = 0.579334893818833)

Extract the posterior. (Here the transformation was not really necessary).

β_posterior = first.(transform.(t, eachcol(pool_posterior_matrices(results))))
5000-element Vector{Vector{Float64}}:
 [0.87944345589263, 2.1792645205916976]
 [0.8647223862783247, 2.048174003240724]
 [0.7576969020742198, 2.04872210995903]
 [0.8361628405388281, 2.0959118764729183]
 [0.8817761602469009, 2.0149925354649496]
 [0.9104967267297177, 2.2327501754428543]
 [1.0085625237238838, 2.0444133726832785]
 [0.9160855823377348, 2.034649712129222]
 [0.9024301848360111, 2.1249342942684537]
 [0.9487021998405073, 2.0319001284803884]
 ⋮
 [0.845753931911677, 1.9392259081660268]
 [0.8737430286119123, 1.924181321877326]
 [1.0003045586582562, 1.9062854706963746]
 [0.9482112174681953, 1.915100093291728]
 [0.8199984098452495, 2.0616608827710796]
 [0.8039386318117738, 2.0795893614441283]
 [0.9239548071095479, 2.170728588354794]
 [0.7516762857269293, 1.8507775796630035]
 [0.9709467729410717, 2.305077964934351]

Check that we recover the parameters.

mean(β_posterior)
2-element Vector{Float64}:
 0.9030311052635918
 2.071722336197655

Quantiles

qs = [0.05, 0.25, 0.5, 0.75, 0.95]
quantile(first.(β_posterior), qs)

quantile(last.(β_posterior), qs)
5-element Vector{Float64}:
 1.8597942447342857
 1.9811928364650786
 2.07113752454796
 2.161025549208615
 2.288842316129517

Check that mixing is good.

ess, R̂ = ess_rhat(stack_posterior_matrices(results))
(ess = [2772.55834574209, 2424.1360124192865], rhat = [1.0007108964141362, 1.003441227345371])

This page was generated using Literate.jl.