m5.6d · DynamicHMCModels.jl

Load Julia packages (libraries) needed for the snippets in chapter 0

using DynamicHMCModels

CmdStan uses a tmp directory to store the output of cmdstan

ProjDir = rel_path_d("..", "scripts", "05")
cd(ProjDir)

Read the milk data

wd = CSV.read(rel_path("..", "data", "milk.csv"), delim=';')
df = convert(DataFrame, wd);
dcc = filter(row -> !(row[:neocortex_perc] == "NA"), df)
dcc[:kcal_per_g] = convert(Vector{Float64}, dcc[:kcal_per_g])
dcc[:log_mass] = log.(convert(Vector{Float64}, dcc[:mass]))

17-element Array{Float64,1}:
  0.6678293725756554
  1.6582280766035324
  1.6808279085207734
  0.9202827531436925
 -0.3856624808119846
 -2.120263536200091
 -0.7550225842780328
 -1.1394342831883648
  0.4382549309311553
  1.1755733298042381
  2.509599262378372
  1.6808279085207734
  3.5689691574413787
  4.374876130645041
  3.70721041079866
  3.4998353515591547
  4.006423680849631

Show first 5 rows

first(dcc[[3, 7, 9]], 5)

5 rows × 3 columns

	kcal_per_g	mass	log_mass
	Float64	Float64⍰	Float64
1	0.49	1.95	0.667829
2	0.47	5.25	1.65823
3	0.56	5.37	1.68083
4	0.89	2.51	0.920283
5	0.92	0.68	-0.385662

Define the model struct

struct m_5_6{TY <: AbstractVector, TX <: AbstractMatrix}
    "Observations."
    y::TY
    "Covariates"
    X::TX
end

Make the type callable with the parameters as a single argument.

function (problem::m_5_6)(θ)
    @unpack y, X, = problem   # extract the data
    @unpack β, σ = θ            # works on the named tuple too
    ll = 0.0
    ll += logpdf(Normal(0, 100), X[1]) # a = X[1]
    ll += logpdf(Normal(0, 1), X[2]) # b1 = X[2]
    ll += logpdf(TDist(1.0), σ)
    ll += loglikelihood(Normal(0, σ), y .- X*β)
    ll
end

Instantiate the model with data and inits.

N = size(dcc, 1)
X = hcat(ones(N), dcc[:log_mass]);
y = dcc[:kcal_per_g]
p = m_5_6(y, X);
p((β = [1.0, 2.0], σ = 1.0))

-242.8761844035513

Write a function to return properly dimensioned transformation.

problem_transformation(p::m_5_6) =
    as((β = as(Array, size(p.X, 2)), σ = asℝ₊))

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

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

LogDensityRejectErrors{InvalidLogDensityException,LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ShiftedExp{true,Float64}}},Main.ex-m5.6d.m_5_6{Array{Float64,1},Array{Float64,2}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ShiftedExp{true,Float64}}},Main.ex-m5.6d.m_5_6{Array{Float64,1},Array{Float64,2}}}},Float64},Float64,3,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ShiftedExp{true,Float64}}},Main.ex-m5.6d.m_5_6{Array{Float64,1},Array{Float64,2}}}},Float64},Float64,3},1}}}}(ForwardDiff AD wrapper for TransformedLogDensity of dimension 3, w/ chunk size 3)

Tune and sample.

chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000);

MCMC, adapting ϵ (75 steps)
8.2e-5 s/step ...done
MCMC, adapting ϵ (25 steps)
7.6e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
0.00049 s/step ...done
MCMC, adapting ϵ (100 steps)
3.6e-5 s/step ...done
MCMC, adapting ϵ (200 steps)
2.6e-5 s/step ...done
MCMC, adapting ϵ (400 steps)
2.6e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
3.9e-5 s/step ...done
MCMC (1000 steps)
4.3e-5 s/step ...done
(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([0.7660167498350918, -0.03067551385655505, -2.0142316537717546], -7.261585381393523, 3, DynamicHMC.DoubledTurn, 0.9327334134105877, 7), NUTS_Transition{Array{Float64,1},Float64}([0.6679867364835417, -0.03732587763263839, -1.6460952049542756], -6.855413645619663, 3, DynamicHMC.DoubledTurn, 0.9987616953363678, 7), NUTS_Transition{Array{Float64,1},Float64}([0.6936497891860188, -0.06088410612602795, -1.7128365738556521], -4.539634972815555, 3, DynamicHMC.DoubledTurn, 0.9819894656629481, 7), NUTS_Transition{Array{Float64,1},Float64}([0.7344745424115171, -0.008136137163578975, -1.9105027415111635], -6.917262260869699, 3, DynamicHMC.DoubledTurn, 0.896229212685492, 7), NUTS_Transition{Array{Float64,1},Float64}([0.7325104075129135, -0.011176448247363361, -1.820877212614075], -5.112432888293861, 2, DynamicHMC.DoubledTurn, 1.0, 3), NUTS_Transition{Array{Float64,1},Float64}([0.6973225512645181, -0.008123504244726534, -1.7017384511698943], -5.146131956317651, 3, DynamicHMC.DoubledTurn, 0.9816105133441929, 7), NUTS_Transition{Array{Float64,1},Float64}([0.7135565014372451, -0.04696549325498229, -1.8490354214382658], -3.729812137569709, 3, DynamicHMC.DoubledTurn, 0.9931012881839074, 7), NUTS_Transition{Array{Float64,1},Float64}([0.6592009197600381, -0.04797563140187742, -1.930682037531562], -5.585380335913473, 3, DynamicHMC.DoubledTurn, 0.8473860166173888, 7), NUTS_Transition{Array{Float64,1},Float64}([0.7793439620466398, -0.03291042032686272, -1.3797442880302175], -6.847093983755506, 3, DynamicHMC.DoubledTurn, 0.977417190616957, 7), NUTS_Transition{Array{Float64,1},Float64}([0.560453528472016, 0.00854288639627449, -1.8138000193588577], -8.12411925030281, 3, DynamicHMC.DoubledTurn, 0.9561368879742211, 7)  …  NUTS_Transition{Array{Float64,1},Float64}([0.6557220178319205, -0.0058719182984590486, -1.7219197064506484], -3.727388700159157, 2, DynamicHMC.AdjacentTurn, 0.9890118950496625, 7), NUTS_Transition{Array{Float64,1},Float64}([0.696261699898732, -0.041976765208038884, -1.8227118021104554], -3.6475677532710145, 2, DynamicHMC.AdjacentTurn, 0.9963415496908918, 7), NUTS_Transition{Array{Float64,1},Float64}([0.6661546809471731, -0.043159482935956985, -1.9501721756973431], -4.804586327183684, 2, DynamicHMC.DoubledTurn, 0.8430696504027972, 3), NUTS_Transition{Array{Float64,1},Float64}([0.6414621270669906, -0.025670021965331243, -1.8444996019439515], -5.078711497552484, 3, DynamicHMC.DoubledTurn, 1.0, 7), NUTS_Transition{Array{Float64,1},Float64}([0.7197605314234589, -0.026937298847983987, -1.5447473128795484], -4.515216286313192, 3, DynamicHMC.DoubledTurn, 0.9964767409912343, 7), NUTS_Transition{Array{Float64,1},Float64}([0.7161346544519858, -0.018942061688723245, -1.6370736729646806], -3.9558006038196343, 2, DynamicHMC.AdjacentTurn, 0.995002356585693, 7), NUTS_Transition{Array{Float64,1},Float64}([0.6655196958880732, -0.053358273382678886, -1.80086109633995], -8.141312774332352, 3, DynamicHMC.DoubledTurn, 0.8036202958747404, 7), NUTS_Transition{Array{Float64,1},Float64}([0.6737865280631026, -0.003028607177065464, -1.7211433185887997], -5.362071152494633, 3, DynamicHMC.DoubledTurn, 0.9996749862385598, 7), NUTS_Transition{Array{Float64,1},Float64}([0.7352826534630607, -0.05918728501645078, -1.8694444072010938], -4.738230240075646, 3, DynamicHMC.DoubledTurn, 0.9856482973496364, 7), NUTS_Transition{Array{Float64,1},Float64}([0.7067548558507162, -0.010727040429838, -1.8323383005944658], -5.56280254081228, 3, DynamicHMC.DoubledTurn, 0.9659791030949462, 7)], NUTS sampler in 3 dimensions
  stepsize (ϵ) ≈ 0.602
  maximum depth = 10
  Gaussian kinetic energy, √diag(M⁻¹): [0.058868933149553, 0.02384099774198216, 0.17401953523976757]
)

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

posterior = TransformVariables.transform.(Ref(problem_transformation(p)), get_position.(chain));
posterior[1:5]

5-element Array{NamedTuple{(:β, :σ),Tuple{Array{Float64,1},Float64}},1}:
 (β = [0.7660167498350918, -0.03067551385655505], σ = 0.13342287895834978)
 (β = [0.6679867364835417, -0.03732587763263839], σ = 0.19280129019326991)
 (β = [0.6936497891860188, -0.06088410612602795], σ = 0.18035348038765234)
 (β = [0.7344745424115171, -0.008136137163578975], σ = 0.14800595914858294)
 (β = [0.7325104075129135, -0.011176448247363361], σ = 0.1618836822226995)

Extract the parameter posterior means: β,

posterior_β = mean(first, posterior)

2-element Array{Float64,1}:
  0.7048584329283026
 -0.032422053825325516

then σ:

posterior_σ = mean(last, posterior)

0.18142338883757908

Effective sample sizes (of untransformed draws)

ess = mapslices(effective_sample_size,
                get_position_matrix(chain); dims = 1)

NUTS-specific statistics

NUTS_statistics(chain)

Hamiltonian Monte Carlo sample of length 1000
  acceptance rate mean: 0.95, min/25%/median/75%/max: 0.55 0.93 0.97 0.99 1.0
  termination: AdjacentTurn => 16% DoubledTurn => 84%
  depth: 1 => 3% 2 => 38% 3 => 59%

cmdstan result

cmdstan_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
          Mean         SD        Naive SE       MCSE      ESS
    a  0.70472876 0.057040655 0.00090189195 0.0011398893 1000
   bm -0.03150330 0.023642759 0.00037382484 0.0004712342 1000
sigma  0.18378372 0.039212805 0.00062000888 0.0011395979 1000

Quantiles:
          2.5%       25.0%       50.0%        75.0%       97.5%
    a  0.59112968  0.66848775  0.70444950  0.741410500 0.81915225
   bm -0.07729257 -0.04708425 -0.03104865 -0.015942925 0.01424901
sigma  0.12638780  0.15605950  0.17800600  0.204319250 0.27993590
";

"\nIterations = 1:1000\nThinning interval = 1\nChains = 1,2,3,4\nSamples per chain = 1000\n\nEmpirical Posterior Estimates:\n          Mean         SD        Naive SE       MCSE      ESS\n    a  0.70472876 0.057040655 0.00090189195 0.0011398893 1000\n   bm -0.03150330 0.023642759 0.00037382484 0.0004712342 1000\nsigma  0.18378372 0.039212805 0.00062000888 0.0011395979 1000\n\nQuantiles:\n          2.5%       25.0%       50.0%        75.0%       97.5%\n    a  0.59112968  0.66848775  0.70444950  0.741410500 0.81915225\n   bm -0.07729257 -0.04708425 -0.03104865 -0.015942925 0.01424901\nsigma  0.12638780  0.15605950  0.17800600  0.204319250 0.27993590\n"

Extract the parameter posterior means: [β, σ],

[posterior_β, posterior_σ]

2-element Array{Any,1}:
  [0.7048584329283026, -0.032422053825325516]
 0.18142338883757908

End of 05/5.6d.jl

This page was generated using Literate.jl.