m8.1d · 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", "08")
cd(ProjDir)

snippet 5.1

d = CSV.read(rel_path("..", "data", "rugged.csv"), delim=';');
df = convert(DataFrame, d);

dcc = filter(row -> !(ismissing(row[:rgdppc_2000])), df)
dcc[!, :log_gdp] = log.(dcc[!, :rgdppc_2000])
dcc[!, :cont_africa] = Array{Float64}(convert(Array{Int}, dcc[!, :cont_africa]))

170-element Array{Float64,1}:
 1.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 1.0
 ⋮  
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 1.0
 1.0
 1.0

First 5 rows with data

first(dcc[!, [:rugged, :cont_africa, :log_gdp]], 5)

struct m_8_1_model{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_8_1_model)(θ)
    @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, 10), X[2]) # bR = X[2]
    ll += logpdf(Normal(0, 10), X[3]) # bA = X[3]
    ll += logpdf(Normal(0, 10), X[4]) # bAR = X[4]
    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[!, :rugged], dcc[!, :cont_africa], dcc[!, :rugged].*dcc[!, :cont_africa]);
y = convert(Vector{Float64}, dcc[!, :log_gdp])
p = m_8_1_model(y, X);
p((β = [1.0, 2.0, 1.0, 2.0], σ = 1.0))

-2770.7279383721293

Write a function to return properly dimensioned transformation.

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

problem_transformation (generic function with 1 method)

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

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

ForwardDiff AD wrapper for TransformedLogDensity of dimension 5, w/ chunk size 5

Tune and sample.

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

(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([9.259743274414754, -0.21985164837307386, -2.228464113060997, 0.467306005393707, -0.0074551199983883334], -248.49145614398196, 2, DynamicHMC.AdjacentTurn, 0.9851456717458211, 7), NUTS_Transition{Array{Float64,1},Float64}([9.060852311416351, -0.12926117516342084, -1.4262832577978206, 0.273268691801494, -0.07125530520315583], -253.0321129523481, 3, DynamicHMC.DoubledTurn, 0.8136593911137707, 7), NUTS_Transition{Array{Float64,1},Float64}([9.168574797906373, -0.2527242176049811, -1.8462644829890758, 0.4849699493704176, -0.033428582584551336], -250.46251982061, 2, DynamicHMC.DoubledTurn, 1.0, 3), NUTS_Transition{Array{Float64,1},Float64}([9.14578591240642, -0.2764847352810241, -1.90458326145807, 0.39775385153987936, -0.037230603425289006], -251.3368825093078, 2, DynamicHMC.DoubledTurn, 0.8025255825552358, 3), NUTS_Transition{Array{Float64,1},Float64}([9.305514355527977, -0.3571986717708979, -2.048244767995564, 0.5250481370075013, -0.02358184936168269], -251.80855686146623, 2, DynamicHMC.AdjacentTurn, 0.9812527135330881, 7), NUTS_Transition{Array{Float64,1},Float64}([9.085692001327148, -0.08752522163343948, -1.6519203080630656, 0.10046661253010765, 0.019671998156646824], -253.21023297009856, 3, DynamicHMC.DoubledTurn, 0.9836129529969568, 7), NUTS_Transition{Array{Float64,1},Float64}([9.436912464762965, -0.33400375886627326, -2.480477224694146, 0.7546412441651487, -0.11189453902587662], -253.57256651713877, 3, DynamicHMC.DoubledTurn, 0.9138219193745495, 7), NUTS_Transition{Array{Float64,1},Float64}([9.437743065133393, -0.31090048528688524, -1.961323282014375, 0.5092102231499678, -0.08480914614854856], -253.1210849046568, 3, DynamicHMC.DoubledTurn, 1.0, 7), NUTS_Transition{Array{Float64,1},Float64}([9.262588213417143, -0.19983797782351737, -2.0424681081125984, 0.5068240531065921, -0.07905299184144908], -251.60389246943444, 2, DynamicHMC.DoubledTurn, 0.8856921112207438, 3), NUTS_Transition{Array{Float64,1},Float64}([9.210175784658063, -0.21902707089582754, -1.8718149873064465, 0.2808132208328478, -0.025100521101536648], -249.62714694914504, 3, DynamicHMC.DoubledTurn, 0.8624317972270211, 7)  …  NUTS_Transition{Array{Float64,1},Float64}([9.164279382260926, -0.1946018019244994, -2.113107787396484, 0.5280998543729035, -0.07425457807229888], -251.95309202967135, 3, DynamicHMC.DoubledTurn, 0.94449621177122, 7), NUTS_Transition{Array{Float64,1},Float64}([9.073901956304056, -0.116295411403447, -1.7819059585117996, 0.24286127827403425, -0.07861928010988282], -248.56059326992224, 3, DynamicHMC.DoubledTurn, 0.9942622782743638, 7), NUTS_Transition{Array{Float64,1},Float64}([9.353625613640233, -0.2801587463212483, -2.114526733502246, 0.5515941674891409, -0.05099972103823799], -247.9727548914604, 3, DynamicHMC.DoubledTurn, 0.9957903755204078, 7), NUTS_Transition{Array{Float64,1},Float64}([9.142461663635467, -0.1312776147900947, -1.709393329834966, 0.21776125524598494, -0.11438862659227964], -248.54647592321348, 3, DynamicHMC.DoubledTurn, 0.9476201811249287, 7), NUTS_Transition{Array{Float64,1},Float64}([9.317918712883717, -0.19443232125119883, -1.9380784615853854, 0.39383386543048404, -0.09930473737578996], -249.80657913140692, 3, DynamicHMC.DoubledTurn, 0.8947526846706795, 7), NUTS_Transition{Array{Float64,1},Float64}([9.282894750552977, -0.19437092146215548, -1.7905051225482536, 0.42477171117864915, -0.06404352789564008], -250.46645499842296, 2, DynamicHMC.DoubledTurn, 0.8493095040750975, 3), NUTS_Transition{Array{Float64,1},Float64}([9.225044073662097, -0.24412814171285097, -1.9385133309291755, 0.2112374411234422, -0.08263953487620854], -252.04413853371673, 3, DynamicHMC.DoubledTurn, 0.9405452298742186, 7), NUTS_Transition{Array{Float64,1},Float64}([9.25572379879544, -0.284620961866809, -2.1326810212267384, 0.6720180945116022, -0.1606849202797319], -254.56508425913142, 3, DynamicHMC.DoubledTurn, 0.9267528185835889, 7), NUTS_Transition{Array{Float64,1},Float64}([9.122465192964691, -0.13658838790423733, -1.799504933087528, 0.36830836662229094, -0.05876459504610845], -251.23050242282196, 2, DynamicHMC.DoubledTurn, 1.0, 3), NUTS_Transition{Array{Float64,1},Float64}([9.500981765409975, -0.36938394862721974, -1.977153741859942, 0.5064742061730746, -0.09226138321782884], -252.83149739378322, 3, DynamicHMC.DoubledTurn, 0.7307951155710583, 7)], NUTS sampler in 5 dimensions
  stepsize (ϵ) ≈ 0.691
  maximum depth = 10
  Gaussian kinetic energy, √diag(M⁻¹): [0.13663659233386266, 0.08257438397606842, 0.22507115125409952, 0.13370781898943024, 0.05790189497228747]
)

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}:
 (β = [9.259743274414754, -0.21985164837307386, -2.228464113060997, 0.467306005393707], σ = 0.9925726004794355) 
 (β = [9.060852311416351, -0.12926117516342084, -1.4262832577978206, 0.273268691801494], σ = 0.9312241154140981)
 (β = [9.168574797906373, -0.2527242176049811, -1.8462644829890758, 0.4849699493704176], σ = 0.9671239782600249)
 (β = [9.14578591240642, -0.2764847352810241, -1.90458326145807, 0.39775385153987936], σ = 0.963453933952559)   
 (β = [9.305514355527977, -0.3571986717708979, -2.048244767995564, 0.5250481370075013], σ = 0.9766940296142794)

Extract the parameter posterior means: β,

posterior_β = mean(first, posterior)

4-element Array{Float64,1}:
  9.225975760047557  
 -0.20453796246352907
 -1.9499601907821267 
  0.3970206897599782

then σ:

posterior_σ = mean(last, posterior)

0.9464996307998951

Effective sample sizes (of untransformed draws)

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

1×5 Array{Float64,2}:
 1000.0  1000.0  1000.0  1000.0  1000.0

NUTS-specific statistics

NUTS_statistics(chain)

Hamiltonian Monte Carlo sample of length 1000
  acceptance rate mean: 0.91, min/25%/median/75%/max: 0.51 0.86 0.94 0.98 1.0
  termination: AdjacentTurn => 18% DoubledTurn => 82%
  depth: 1 => 0% 2 => 43% 3 => 56% 4 => 0%

Result rethinking

rethinking = "
       mean   sd  5.5% 94.5% n_eff Rhat
a      9.22 0.14  9.00  9.46   282    1
bR    -0.21 0.08 -0.33 -0.08   275    1
bA    -1.94 0.24 -2.33 -1.59   268    1
bAR    0.40 0.14  0.18  0.62   271    1
sigma  0.96 0.05  0.87  1.04   339    1
"

"\n       mean   sd  5.5% 94.5% n_eff Rhat\na      9.22 0.14  9.00  9.46   282    1\nbR    -0.21 0.08 -0.33 -0.08   275    1\nbA    -1.94 0.24 -2.33 -1.59   268    1\nbAR    0.40 0.14  0.18  0.62   271    1\nsigma  0.96 0.05  0.87  1.04   339    1\n"

Summary

[posterior_β, posterior_σ]

2-element Array{Any,1}:
  [9.225975760047557, -0.20453796246352907, -1.9499601907821267, 0.3970206897599782]
 0.9464996307998951

End of 08/m8.1s.jl

This page was generated using Literate.jl.