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((β = $β,))

-442.0958140546533

Sample using NUTS, random starting point.

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 = [0.9057322257501476 0.8752679781282163 … 1.081989392236934 0.9659275918008315; 2.0698856647660255 2.087665605543101 … 2.340439330491711 2.2381915172343634], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-443.3528567183514, 2, turning at positions 0:3, 0.9445557575988827, 3, DynamicHMC.Directions(0x5182015f)), DynamicHMC.TreeStatisticsNUTS(-441.5725186159271, 3, turning at positions -5:2, 0.9480620033673984, 7, DynamicHMC.Directions(0xc85f9732)), DynamicHMC.TreeStatisticsNUTS(-441.7320447795951, 2, turning at positions -2:-5, 0.988118190522935, 7, DynamicHMC.Directions(0x2ed3903a)), DynamicHMC.TreeStatisticsNUTS(-444.2064327233359, 2, turning at positions 1:4, 0.8767335115252227, 7, DynamicHMC.Directions(0x44a947dc)), DynamicHMC.TreeStatisticsNUTS(-448.2071145757616, 2, turning at positions 0:3, 0.7488347970349228, 3, DynamicHMC.Directions(0x74f24abb)), DynamicHMC.TreeStatisticsNUTS(-447.3634116562838, 1, turning at positions 0:1, 1.0, 1, DynamicHMC.Directions(0xec0b0775)), DynamicHMC.TreeStatisticsNUTS(-446.5274667130504, 2, turning at positions -1:2, 0.9419992194563586, 3, DynamicHMC.Directions(0xf526c21e)), DynamicHMC.TreeStatisticsNUTS(-446.37789410143046, 2, turning at positions -3:0, 0.9395219018589276, 3, DynamicHMC.Directions(0x8803cfc8)), DynamicHMC.TreeStatisticsNUTS(-444.4717075689284, 3, turning at positions -3:4, 0.9999999999999999, 7, DynamicHMC.Directions(0xec59361c)), DynamicHMC.TreeStatisticsNUTS(-442.5212534151549, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0xf331e566))  …  DynamicHMC.TreeStatisticsNUTS(-444.6569715953046, 1, turning at positions 0:1, 1.0, 1, DynamicHMC.Directions(0x97f9e5c7)), DynamicHMC.TreeStatisticsNUTS(-444.5780982336841, 2, turning at positions -4:-5, 0.984122274725701, 5, DynamicHMC.Directions(0x4753e3c8)), DynamicHMC.TreeStatisticsNUTS(-442.02399370147555, 2, turning at positions -1:2, 0.9922188287564989, 3, DynamicHMC.Directions(0xcc748372)), DynamicHMC.TreeStatisticsNUTS(-441.8239880917822, 3, turning at positions -1:6, 0.9937090530455414, 7, DynamicHMC.Directions(0x724cb7f6)), DynamicHMC.TreeStatisticsNUTS(-441.2463884970516, 2, turning at positions -2:1, 0.9950741052725235, 3, DynamicHMC.Directions(0xd371a7c1)), DynamicHMC.TreeStatisticsNUTS(-442.086152405855, 2, turning at positions -3:0, 0.9533984620216637, 3, DynamicHMC.Directions(0xcf2b1bac)), DynamicHMC.TreeStatisticsNUTS(-442.0752523452796, 2, turning at positions -1:2, 0.9995437988520789, 3, DynamicHMC.Directions(0xf959955a)), DynamicHMC.TreeStatisticsNUTS(-444.51767779384454, 2, turning at positions -3:0, 0.8103095026269881, 3, DynamicHMC.Directions(0x69a19204)), DynamicHMC.TreeStatisticsNUTS(-444.88667014785193, 2, turning at positions -3:0, 0.9999999999999999, 3, DynamicHMC.Directions(0xab99eb64)), DynamicHMC.TreeStatisticsNUTS(-445.68844264411405, 3, turning at positions 0:7, 0.9491656736376227, 7, DynamicHMC.Directions(0x35595f77))], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.09859876161419211, 0.1364409718719911], ϵ = 0.6649562272646156)
 (posterior_matrix = [0.8680389315238533 0.9664918910126306 … 0.9655759747225012 0.854318350943319; 2.0745742741377855 2.0795199583839286 … 2.0080975391058846 2.036809331962731], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-441.37236695288294, 2, turning at positions 5:6, 0.9995407148165915, 7, DynamicHMC.Directions(0x70dd045e)), DynamicHMC.TreeStatisticsNUTS(-441.42375252578125, 2, turning at positions -3:0, 0.9919669305076559, 3, DynamicHMC.Directions(0xe769bca4)), DynamicHMC.TreeStatisticsNUTS(-443.12215863771576, 2, turning at positions -2:-5, 0.9190122858612002, 7, DynamicHMC.Directions(0x76c5f462)), DynamicHMC.TreeStatisticsNUTS(-443.68193144993614, 3, turning at positions -2:5, 0.9035353238187245, 7, DynamicHMC.Directions(0xbb3cf07d)), DynamicHMC.TreeStatisticsNUTS(-444.4932652737158, 2, turning at positions -2:1, 0.9552951515905713, 3, DynamicHMC.Directions(0x3b227481)), DynamicHMC.TreeStatisticsNUTS(-445.11994644002203, 3, turning at positions -4:3, 0.9828357024991398, 7, DynamicHMC.Directions(0x7049a4fb)), DynamicHMC.TreeStatisticsNUTS(-442.9902592024522, 2, turning at positions -3:0, 0.9875734026513524, 3, DynamicHMC.Directions(0xb1eec538)), DynamicHMC.TreeStatisticsNUTS(-446.1693414503824, 2, turning at positions 2:5, 0.754882384097319, 7, DynamicHMC.Directions(0x5546dc3d)), DynamicHMC.TreeStatisticsNUTS(-445.6267111968051, 2, turning at positions -3:0, 0.9214545148566867, 3, DynamicHMC.Directions(0xcf8cdf54)), DynamicHMC.TreeStatisticsNUTS(-446.7134979526038, 3, turning at positions -4:3, 0.8299499167138817, 7, DynamicHMC.Directions(0x0defcd03))  …  DynamicHMC.TreeStatisticsNUTS(-442.4555606214261, 2, turning at positions -2:1, 0.9544693876495488, 3, DynamicHMC.Directions(0x31b43939)), DynamicHMC.TreeStatisticsNUTS(-442.3740804977069, 2, turning at positions -3:-6, 0.9364269567361699, 7, DynamicHMC.Directions(0x07339c09)), DynamicHMC.TreeStatisticsNUTS(-443.34401365782986, 2, turning at positions 0:3, 0.8096203205897968, 3, DynamicHMC.Directions(0xe33a272f)), DynamicHMC.TreeStatisticsNUTS(-443.48497342081873, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0xdb785be2)), DynamicHMC.TreeStatisticsNUTS(-442.0381069465855, 3, turning at positions -7:0, 0.9921880883844194, 7, DynamicHMC.Directions(0x69a09c00)), DynamicHMC.TreeStatisticsNUTS(-443.52789837960785, 1, turning at positions 2:3, 0.8912095963359657, 3, DynamicHMC.Directions(0xf9fa9d23)), DynamicHMC.TreeStatisticsNUTS(-442.21751069467683, 1, turning at positions 0:1, 1.0, 1, DynamicHMC.Directions(0xc6f0c1bd)), DynamicHMC.TreeStatisticsNUTS(-442.1512115483032, 3, turning at positions -4:3, 0.9974897177515113, 7, DynamicHMC.Directions(0x35a6d55b)), DynamicHMC.TreeStatisticsNUTS(-442.84191981254264, 1, turning at positions 1:2, 0.8929360575709236, 3, DynamicHMC.Directions(0xe7cd4a92)), DynamicHMC.TreeStatisticsNUTS(-441.7398168025174, 2, turning at positions 4:7, 0.9955648068416462, 7, DynamicHMC.Directions(0x3d1c676f))], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.09444138639632087, 0.14191598872483382], ϵ = 0.7259540380709076)
 (posterior_matrix = [0.9300005743967242 1.0469861723794402 … 0.9995717115022427 0.9166429882059818; 2.1332074966789634 2.2760111600973976 … 1.979633480318152 2.121782736485274], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-444.32129990269414, 2, turning at positions -2:1, 0.9999999999999999, 3, DynamicHMC.Directions(0x570c5c59)), DynamicHMC.TreeStatisticsNUTS(-442.4998709888647, 2, turning at positions -3:0, 0.9024303396394538, 3, DynamicHMC.Directions(0x71d85400)), DynamicHMC.TreeStatisticsNUTS(-444.10107319654065, 3, turning at positions -7:0, 0.8934132870263864, 7, DynamicHMC.Directions(0xcad3ed38)), DynamicHMC.TreeStatisticsNUTS(-444.53460602675847, 2, turning at positions -1:2, 0.7334697982860136, 3, DynamicHMC.Directions(0x6473c8d2)), DynamicHMC.TreeStatisticsNUTS(-442.4375388439867, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0xb579a316)), DynamicHMC.TreeStatisticsNUTS(-442.7923386684332, 2, turning at positions -3:0, 0.9264636405869505, 3, DynamicHMC.Directions(0x1b583428)), DynamicHMC.TreeStatisticsNUTS(-442.41047741992253, 2, turning at positions -2:1, 0.9895275525369986, 3, DynamicHMC.Directions(0xd0492d9d)), DynamicHMC.TreeStatisticsNUTS(-441.8367253317915, 1, turning at positions 0:1, 0.9418813860719479, 1, DynamicHMC.Directions(0x92b46f6f)), DynamicHMC.TreeStatisticsNUTS(-442.1375189371335, 3, turning at positions 0:7, 0.9962231840938022, 7, DynamicHMC.Directions(0x8c689a4f)), DynamicHMC.TreeStatisticsNUTS(-442.4126192711258, 2, turning at positions -2:1, 0.9109099202109029, 3, DynamicHMC.Directions(0x92e1e645))  …  DynamicHMC.TreeStatisticsNUTS(-443.23208072729636, 2, turning at positions -2:-5, 0.9722059577990623, 7, DynamicHMC.Directions(0x082e8fea)), DynamicHMC.TreeStatisticsNUTS(-443.182341631735, 2, turning at positions -3:0, 0.9604741821091735, 3, DynamicHMC.Directions(0x877d17d4)), DynamicHMC.TreeStatisticsNUTS(-444.8212267639885, 2, turning at positions -3:0, 0.6447940135481102, 3, DynamicHMC.Directions(0x7d7e6548)), DynamicHMC.TreeStatisticsNUTS(-441.7767212313912, 1, turning at positions 2:3, 0.9999999999999999, 3, DynamicHMC.Directions(0xd0d6e89b)), DynamicHMC.TreeStatisticsNUTS(-441.46752458324903, 2, turning at positions 0:3, 0.9841060996753638, 3, DynamicHMC.Directions(0x162d6a8b)), DynamicHMC.TreeStatisticsNUTS(-442.14902419701457, 2, turning at positions 0:3, 0.8630806885044894, 3, DynamicHMC.Directions(0x157d274f)), DynamicHMC.TreeStatisticsNUTS(-441.93622244085407, 2, turning at positions 0:3, 0.87930819743465, 3, DynamicHMC.Directions(0x774fbda7)), DynamicHMC.TreeStatisticsNUTS(-441.60217301382494, 2, turning at positions -2:1, 0.9863206574623171, 3, DynamicHMC.Directions(0xa77294d1)), DynamicHMC.TreeStatisticsNUTS(-442.9271414999545, 2, turning at positions -3:0, 0.7640950549142134, 3, DynamicHMC.Directions(0x98612aa4)), DynamicHMC.TreeStatisticsNUTS(-442.41517078821255, 1, turning at positions 2:3, 0.9651911949218475, 3, DynamicHMC.Directions(0x8514eb17))], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.08925736349698302, 0.13477170039858566], ϵ = 0.8094384027458571)
 (posterior_matrix = [1.0030030352565988 1.0030030352565988 … 0.9867499302097263 1.0359124883679902; 2.371792974245721 2.371792974245721 … 2.188550450962671 2.2447280798094913], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-443.64048447238906, 2, turning at positions -2:1, 0.9698900620100616, 3, DynamicHMC.Directions(0xb89e94e1)), DynamicHMC.TreeStatisticsNUTS(-444.35224563456603, 1, turning at positions 0:1, 0.8431875642366928, 1, DynamicHMC.Directions(0x96f20d6b)), DynamicHMC.TreeStatisticsNUTS(-443.2853078951459, 2, turning at positions 0:3, 0.9995984901112641, 3, DynamicHMC.Directions(0xe13e0b4f)), DynamicHMC.TreeStatisticsNUTS(-442.9288443544148, 3, turning at positions -1:6, 0.9634019309335874, 7, DynamicHMC.Directions(0x2834c10e)), DynamicHMC.TreeStatisticsNUTS(-442.6770945903042, 1, turning at positions 0:1, 0.9631774962153405, 1, DynamicHMC.Directions(0x0a1314ff)), DynamicHMC.TreeStatisticsNUTS(-443.842099703705, 2, turning at positions -4:-5, 0.8812748397684533, 7, DynamicHMC.Directions(0xa65c85fa)), DynamicHMC.TreeStatisticsNUTS(-442.56401682059885, 2, turning at positions -1:2, 0.9999999999999999, 3, DynamicHMC.Directions(0x78dfbdb2)), DynamicHMC.TreeStatisticsNUTS(-445.4041369321881, 3, turning at positions -2:5, 0.6920656479222848, 7, DynamicHMC.Directions(0x8fef85b5)), DynamicHMC.TreeStatisticsNUTS(-444.2727911605053, 2, turning at positions -4:-7, 0.9562963351738432, 7, DynamicHMC.Directions(0x77af2020)), DynamicHMC.TreeStatisticsNUTS(-444.98983146025563, 2, turning at positions 0:3, 0.8333393850847409, 3, DynamicHMC.Directions(0xe8bc130f))  …  DynamicHMC.TreeStatisticsNUTS(-444.7089587665511, 2, turning at positions -3:0, 0.9114398461845945, 3, DynamicHMC.Directions(0x82abb6c4)), DynamicHMC.TreeStatisticsNUTS(-443.7272173515674, 2, turning at positions -1:2, 0.7768076043796553, 3, DynamicHMC.Directions(0x2bf00376)), DynamicHMC.TreeStatisticsNUTS(-443.33658914353066, 1, turning at positions -1:0, 1.0, 1, DynamicHMC.Directions(0xb180046c)), DynamicHMC.TreeStatisticsNUTS(-444.09102563333516, 3, turning at positions -2:5, 0.9832285972598562, 7, DynamicHMC.Directions(0xac4785ad)), DynamicHMC.TreeStatisticsNUTS(-441.42512738719324, 2, turning at positions -1:2, 0.9913589631483116, 3, DynamicHMC.Directions(0xcd0ec4e6)), DynamicHMC.TreeStatisticsNUTS(-442.9231512829026, 2, turning at positions -1:2, 0.8380180535019042, 3, DynamicHMC.Directions(0x414a2d5a)), DynamicHMC.TreeStatisticsNUTS(-442.84532714068513, 2, turning at positions -3:0, 0.9999999999999999, 3, DynamicHMC.Directions(0xcacab790)), DynamicHMC.TreeStatisticsNUTS(-443.83774902198036, 2, turning at positions 1:4, 0.8935284985176232, 7, DynamicHMC.Directions(0x54be3c44)), DynamicHMC.TreeStatisticsNUTS(-442.4057760022732, 2, turning at positions -2:1, 0.9780363787853594, 3, DynamicHMC.Directions(0x52fb0eb1)), DynamicHMC.TreeStatisticsNUTS(-442.56917751473884, 2, turning at positions -3:0, 0.9409434307900381, 3, DynamicHMC.Directions(0x352209fc))], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.09490765609315233, 0.1317832735244582], ϵ = 0.688749434496165)
 (posterior_matrix = [1.023965195505918 0.9612844060676489 … 1.038612798013524 0.8372371094184493; 2.096776974658813 2.2347616882219263 … 2.0979135839585576 2.1676600553195975], tree_statistics = [DynamicHMC.TreeStatisticsNUTS(-441.93408846099106, 2, turning at positions -2:1, 0.9229687873646237, 3, DynamicHMC.Directions(0x89d5fd69)), DynamicHMC.TreeStatisticsNUTS(-442.67537180773945, 1, turning at positions -2:-3, 0.9237945296409088, 3, DynamicHMC.Directions(0xea894d9c)), DynamicHMC.TreeStatisticsNUTS(-441.87359697226066, 3, turning at positions -3:4, 0.994454049694762, 7, DynamicHMC.Directions(0x0d7ca824)), DynamicHMC.TreeStatisticsNUTS(-442.4125465647762, 1, turning at positions -1:-2, 0.8243943721897288, 3, DynamicHMC.Directions(0x78bf2cc5)), DynamicHMC.TreeStatisticsNUTS(-444.45523627044963, 2, turning at positions -2:-5, 0.9166369557331793, 7, DynamicHMC.Directions(0x3a5ea212)), DynamicHMC.TreeStatisticsNUTS(-442.59804154443725, 3, turning at positions -3:4, 0.9999999999999999, 7, DynamicHMC.Directions(0x3ba86dd4)), DynamicHMC.TreeStatisticsNUTS(-442.24555406456886, 3, turning at positions -2:5, 0.9713780592304351, 7, DynamicHMC.Directions(0x01ea3bcd)), DynamicHMC.TreeStatisticsNUTS(-441.7901645880536, 2, turning at positions -2:1, 0.9477135200482004, 3, DynamicHMC.Directions(0xfdd24bb1)), DynamicHMC.TreeStatisticsNUTS(-442.77563602452506, 2, turning at positions -2:-5, 0.9551041493409107, 7, DynamicHMC.Directions(0x967d6f22)), DynamicHMC.TreeStatisticsNUTS(-443.1330212355244, 2, turning at positions -1:2, 0.9922072970486213, 3, DynamicHMC.Directions(0xa76c6a02))  …  DynamicHMC.TreeStatisticsNUTS(-442.7721119293837, 2, turning at positions -3:0, 0.8927040123809284, 3, DynamicHMC.Directions(0x3671e674)), DynamicHMC.TreeStatisticsNUTS(-443.9719495032516, 2, turning at positions -3:0, 0.9468481945660057, 3, DynamicHMC.Directions(0x126e9710)), DynamicHMC.TreeStatisticsNUTS(-444.0624305160022, 2, turning at positions -3:0, 0.8937550830777706, 3, DynamicHMC.Directions(0x39c5a8c0)), DynamicHMC.TreeStatisticsNUTS(-443.08932499266155, 2, turning at positions 2:3, 0.9021751069999207, 5, DynamicHMC.Directions(0x0d1ef27d)), DynamicHMC.TreeStatisticsNUTS(-443.7040080919868, 3, turning at positions -6:1, 0.8912249850997995, 7, DynamicHMC.Directions(0x0977abd9)), DynamicHMC.TreeStatisticsNUTS(-442.00014884707815, 2, turning at positions -1:2, 0.9859256161677346, 3, DynamicHMC.Directions(0x8baa2062)), DynamicHMC.TreeStatisticsNUTS(-442.9864407716316, 3, turning at positions -1:6, 0.9339052447501205, 7, DynamicHMC.Directions(0xc0bd397e)), DynamicHMC.TreeStatisticsNUTS(-442.75422933100975, 2, turning at positions -3:0, 0.9767157931063353, 3, DynamicHMC.Directions(0x2e57a7a8)), DynamicHMC.TreeStatisticsNUTS(-442.5302427710723, 2, turning at positions -2:1, 0.9999999999999999, 3, DynamicHMC.Directions(0x499be545)), DynamicHMC.TreeStatisticsNUTS(-442.4837120262183, 2, turning at positions -3:0, 0.9782734102479482, 3, DynamicHMC.Directions(0x91caeb94))], κ = Gaussian kinetic energy (Diagonal), √diag(M⁻¹): [0.09806430050584929, 0.1236349114171122], ϵ = 0.7578842425880905)

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.9057322257501476, 2.0698856647660255]
 [0.8752679781282163, 2.087665605543101]
 [0.8402597039008076, 2.029435420448242]
 [1.1131455345514611, 2.3509021824182565]
 [1.2731405154201603, 2.3360131755631364]
 [1.2007524256286601, 2.346046862036971]
 [1.1095775456852723, 2.4554393527487015]
 [1.1725024345443675, 2.21938389276195]
 [1.0241001157755028, 2.2948145423489406]
 [0.984523536186026, 2.2728866805870336]
 ⋮
 [0.9363047422995505, 1.8691745924506784]
 [0.9924990026586501, 2.1172847228672955]
 [0.935119163629553, 1.9586176676293778]
 [0.9649240686619207, 1.9922454819049382]
 [0.9726998146069341, 1.9834644171619502]
 [1.0584331502362228, 2.138866889724439]
 [0.8935991288162878, 2.258989695960221]
 [1.038612798013524, 2.0979135839585576]
 [0.8372371094184493, 2.1676600553195975]

Check that we recover the parameters.

mean(β_posterior)

2-element Vector{Float64}:
 0.9279946260356703
 2.0975251777947475

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.8797198000821063
 2.0028218185172713
 2.092802233488305
 2.189963324353938
 2.331852696914637

Check that mixing is good.

ess, R̂ = ess_rhat(stack_posterior_matrices(results))

(ess = [2179.132637809816, 2062.1423194196213], rhat = [1.0011882272307435, 1.0011149334356764])

This page was generated using Literate.jl.