m10.04d · DynamicHMCModels.jl

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

using DynamicHMCModels, Random
Random.seed!(12345)

MersenneTwister(UInt32[0x00003039], Random.DSFMT.DSFMT_state(Int32[-870096391, 1072918504, -1812426662, 1073255081, -733866021, 1073404543, 807620846, 1073368448, 1919433844, 1072852359  …  -362113007, 1073100625, -166402106, 1073460158, -1907020342, 721295190, -750225566, -1300227565, 382, 0]), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000  …  0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000], 1002, 0)

CmdStan uses a tmp directory to store the output of cmdstan

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

snippet 10.4

d = CSV.read(rel_path("..", "data", "chimpanzees.csv"), delim=';');
df = convert(DataFrame, d);
df[!, :pulled_left] = convert(Array{Int64}, df[!, :pulled_left])
df[!, :prosoc_left] = convert(Array{Int64}, df[!, :prosoc_left])
df[!, :condition] = convert(Array{Int64}, df[!, :condition])
df[!, :actor] = convert(Array{Int64}, df[!, :actor])
first(df, 5)

struct m_10_04d_model{TY <: AbstractVector, TX <: AbstractMatrix,
  TA <: AbstractVector}
    "Observations."
    y::TY
    "Covariates"
    X::TX
    "Actors"
    A::TA
    "Number of observations"
    N::Int
    "Number of unique actors"
    N_actors::Int
end

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

function (problem::m_10_04d_model)(θ)
    @unpack y, X, A, N, N_actors = problem   # extract the data
    @unpack β, α = θ  # works on the named tuple too
    ll = 0.0
    ll += sum(logpdf.(Normal(0, 10), β)) # bp & bpC
    ll += sum(logpdf.(Normal(0, 10), α)) # alpha[1:7]
    ll += sum(
      [loglikelihood(Binomial(1, logistic(α[A[i]] + dot(X[i, :], β))), [y[i]]) for i in 1:N]
    )
    ll
end

Instantiate the model with data and inits.

N = size(df, 1)
N_actors = length(unique(df[!, :actor]))
X = hcat(ones(Int64, N), df[!, :prosoc_left] .* df[!, :condition]);
A = df[!, :actor]
y = df[!, :pulled_left]
p = m_10_04d_model(y, X, A, N, N_actors);
θ = (β = [1.0, 0.0], α = [-1.0, 10.0, -1.0, -1.0, -1.0, 0.0, 2.0])
p(θ)

-305.21943396408915

Write a function to return properly dimensioned transformation.

problem_transformation(p::m_10_04d_model) =
    as( (β = as(Array, size(p.X, 2)), α = as(Array, p.N_actors), ) )

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 = LogDensityRejectErrors(ADgradient(:ForwardDiff, P));
#∇P = ADgradient(:ForwardDiff, P);

#import Zygote
#∇P = LogDensityRejectErrors(ADgradient(:Zygote, P));
#∇P = ADgradient(:Zygote, P);

Tune and sample.

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

(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([4.805608442141108, 0.2389779757501151, -5.4228120431192135, 21.102838682057858, -5.087605475240926, -5.339686571980073, -5.239935978535041, -4.187801436108223, -2.196616929378975], -304.5866229204268, 3, DynamicHMC.AdjacentTurn, 0.768997189893174, 15), NUTS_Transition{Array{Float64,1},Float64}([3.1779855333330698, 0.5965490813806572, -3.332679123628399, 22.4757391092144, -4.052102062224823, -4.302760278141705, -3.5144792187419234, -2.9698061668305353, -0.9186563627898878], -302.46355876847804, 4, DynamicHMC.DoubledTurn, 0.8817693965854141, 15), NUTS_Transition{Array{Float64,1},Float64}([3.6751239531827453, 0.5510236952460794, -4.006795241697499, 17.439170211065864, -4.428966426100576, -4.43289953739397, -3.92856451258628, -3.195481244391613, -1.7185600988784437], -298.25202804076116, 3, DynamicHMC.DoubledTurn, 0.9024392971710403, 7), NUTS_Transition{Array{Float64,1},Float64}([6.811143293480398, 0.2662578097605355, -7.316318256645592, 2.853053665655831, -7.543343290655876, -7.55305997384782, -7.380809100313019, -6.3777111648238165, -4.189699859996347], -296.9038935516846, 4, DynamicHMC.DoubledTurn, 0.9649599897794795, 15), NUTS_Transition{Array{Float64,1},Float64}([-4.431487462956057, 0.45030287799194935, 4.0152969847455955, 12.178060714384907, 3.4452039433839743, 3.6231689167024097, 3.6010109315830148, 4.893472780944509, 5.919666770432145], -296.3356338592542, 4, DynamicHMC.DoubledTurn, 0.9015772416871323, 15), NUTS_Transition{Array{Float64,1},Float64}([4.490182263508208, 0.5065554671206052, -4.943540488632991, 1.464391745910337, -5.313004477982811, -5.339099247194863, -4.876204044081245, -4.085591926342948, -2.4815996696069513], -299.9606063570173, 4, DynamicHMC.DoubledTurn, 0.7631075656010746, 15), NUTS_Transition{Array{Float64,1},Float64}([-0.964622553717225, 0.6066199123840172, 0.3282799508332197, 7.089783469589104, 0.38435984583946103, -0.36474326179547734, -0.05488830200426276, 1.6449033827036064, 3.45100414506115], -300.5782982146635, 3, DynamicHMC.AdjacentTurn, 0.7343085342331114, 15), NUTS_Transition{Array{Float64,1},Float64}([5.100615421804797, 0.39911408670542514, -5.479766108474879, 2.6198228762535187, -6.113476907468543, -5.5512517371641295, -5.0840043247655515, -4.777868627555236, -3.867943388203419], -301.58513763203723, 3, DynamicHMC.AdjacentTurn, 0.9590716763113012, 15), NUTS_Transition{Array{Float64,1},Float64}([3.738219859209382, 0.46733090462215915, -4.342028850044987, 4.533039261561904, -4.250271499770893, -4.851800724124278, -4.424552097090727, -3.2242796648746923, -0.7296647298145693], -299.6963608541688, 4, DynamicHMC.DoubledTurn, 0.9991673955334325, 15), NUTS_Transition{Array{Float64,1},Float64}([4.461766427126224, 0.34309188192690426, -5.015079921874218, 8.368636952189947, -5.109147378700952, -5.350184126009973, -5.1293676024118655, -3.74435858164871, -2.0491135548590287], -299.05636522849744, 4, DynamicHMC.DoubledTurn, 0.9360238627867953, 15)  …  NUTS_Transition{Array{Float64,1},Float64}([-0.38642953192546503, 1.1786523993623372, 0.10454509940070844, 8.28434312999926, -1.011814266408297, -1.056539432900534, -0.10979839312805073, 0.7920875489791117, 2.3570894867914065], -303.75939517144064, 3, DynamicHMC.DoubledTurn, 0.8980544634833463, 7), NUTS_Transition{Array{Float64,1},Float64}([-0.3950342295081849, 0.11350458830841759, 0.23020014136049557, 9.339016118519055, -0.5462965486944091, -0.21691415593067623, -0.07260202415332803, 1.1420424682229546, 2.6346296106390747], -302.2375571761599, 3, DynamicHMC.DoubledTurn, 1.0, 7), NUTS_Transition{Array{Float64,1},Float64}([5.246588022226619, 0.7388497489270798, -5.9509920607123945, 1.1387728250740052, -5.932090558819487, -5.818674770889562, -5.357543023894947, -4.9232544451233125, -3.3702366088479394], -299.40148987508763, 4, DynamicHMC.AdjacentTurn, 0.7834332614572548, 31), NUTS_Transition{Array{Float64,1},Float64}([5.24182170097153, 0.44659319112185725, -5.431041796005965, 2.2329701259037504, -5.756562171285025, -6.315445582997034, -5.555054569821388, -4.877133299511366, -3.1880045618671122], -298.7624038787372, 4, DynamicHMC.DoubledTurn, 0.9876529936578515, 15), NUTS_Transition{Array{Float64,1},Float64}([-2.5081063195101674, 0.27603443159958174, 2.0077630791930976, 12.863973883427384, 1.4038617470916857, 1.7291996369249203, 2.2799062898601705, 2.9870501976620067, 4.456040596858176], -296.5455787425553, 4, DynamicHMC.DoubledTurn, 0.9535432786885711, 15), NUTS_Transition{Array{Float64,1},Float64}([-3.40468654247805, 0.4995366681058604, 3.0759957714473725, 10.416880742658936, 2.1825928159467174, 2.813259682506419, 3.2445807579904447, 4.001406741046448, 5.140347961549395], -297.2487837927424, 3, DynamicHMC.DoubledTurn, 0.9806031086916019, 7), NUTS_Transition{Array{Float64,1},Float64}([-3.267741254973158, 0.40250542489973684, 2.860557172672515, 7.570394337779512, 2.2447285817808895, 2.719509543539649, 3.2281315856912176, 3.979508307501073, 5.1697342510054325], -299.67884873001975, 3, DynamicHMC.DoubledTurn, 0.32392284807998084, 7), NUTS_Transition{Array{Float64,1},Float64}([0.8008515963540074, 0.37470646358002, -1.0325534835836687, 12.238286271593084, -1.3043883778214844, -1.3991125180948512, -1.4568155348859873, -0.3527517530158185, 1.261384309984194], -301.5646408582166, 3, DynamicHMC.DoubledTurn, 0.8947219857853242, 7), NUTS_Transition{Array{Float64,1},Float64}([2.9229350071055156, 0.42383194376865246, -3.5816331293748593, 13.030658178389741, -3.6864428101037525, -3.6134893706820437, -3.0055392653713953, -2.217182038438321, -0.9724604086053438], -301.7700168853612, 3, DynamicHMC.AdjacentTurn, 0.8493194331162169, 15), NUTS_Transition{Array{Float64,1},Float64}([2.593753038372206, 0.5075955776475296, -3.202451958419349, 13.271795091219946, -3.350486137806249, -3.3771453057809726, -2.721798338999288, -1.8985968509889792, -0.5927563278095931], -294.7972872518422, 3, DynamicHMC.DoubledTurn, 0.9999218039793589, 7)], NUTS sampler in 9 dimensions
  stepsize (ϵ) ≈ 0.227
  maximum depth = 10
  Gaussian kinetic energy, √diag(M⁻¹): [3.3502608205751976, 0.4627608063451792, 3.3564056225007555, 6.276846524952666, 3.3567570922034013, 3.3456824123322613, 3.3434606552779105, 3.348013666585637, 3.3748551468304293]
)

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},Array{Float64,1}}},1}:
 (β = [4.805608442141108, 0.2389779757501151], α = [-5.4228120431192135, 21.102838682057858, -5.087605475240926, -5.339686571980073, -5.239935978535041, -4.187801436108223, -2.196616929378975]) 
 (β = [3.1779855333330698, 0.5965490813806572], α = [-3.332679123628399, 22.4757391092144, -4.052102062224823, -4.302760278141705, -3.5144792187419234, -2.9698061668305353, -0.9186563627898878])
 (β = [3.6751239531827453, 0.5510236952460794], α = [-4.006795241697499, 17.439170211065864, -4.428966426100576, -4.43289953739397, -3.92856451258628, -3.195481244391613, -1.7185600988784437])  
 (β = [6.811143293480398, 0.2662578097605355], α = [-7.316318256645592, 2.853053665655831, -7.543343290655876, -7.55305997384782, -7.380809100313019, -6.3777111648238165, -4.189699859996347])   
 (β = [-4.431487462956057, 0.45030287799194935], α = [4.0152969847455955, 12.178060714384907, 3.4452039433839743, 3.6231689167024097, 3.6010109315830148, 4.893472780944509, 5.919666770432145])

Extract the parameter posterior means: β,

posterior_β = mean(first, posterior)

2-element Array{Float64,1}:
 1.451237768661092 
 0.4233537383099833

Extract the parameter posterior means: β,

posterior_α = mean(last, posterior)

7-element Array{Float64,1}:
 -1.9000576019113276
 10.145744904561283 
 -2.194854907777116 
 -2.183382830567691 
 -1.9048052110036864
 -0.971437677726311 
  0.5935301675225381

Effective sample sizes (of untransformed draws)

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

1×9 Array{Float64,2}:
 393.489  1000.0  399.712  261.483  …  394.444  403.401  395.341  396.188

NUTS-specific statistics

NUTS_statistics(chain)

Hamiltonian Monte Carlo sample of length 1000
  acceptance rate mean: 0.91, min/25%/median/75%/max: 0.0 0.89 0.96 0.99 1.0
  termination: AdjacentDivergent => 2% AdjacentTurn => 32% DoubledTurn => 66%
  depth: 0 => 0% 1 => 0% 2 => 2% 3 => 58% 4 => 39% 5 => 1% 6 => 0%

Result rethinking

rethinking = "
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.1 -0.74503184 0.26613979 0.0042080396 0.0060183398 1000
a.2 10.77955494 5.32538998 0.0842018089 0.1269148045 1000
a.3 -1.04982353 0.28535997 0.0045119373 0.0049074219 1000
a.4 -1.04898135 0.28129307 0.0044476339 0.0056325117 1000
a.5 -0.74390933 0.26949936 0.0042611590 0.0052178124 1000
a.6  0.21599365 0.26307574 0.0041595927 0.0045153523 1000
a.7  1.81090866 0.39318577 0.0062168129 0.0071483527 1000
 bp  0.83979926 0.26284676 0.0041559722 0.0059795826 1000
bpC -0.12913322 0.29935741 0.0047332562 0.0049519863 1000
";

"\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\na.1 -0.74503184 0.26613979 0.0042080396 0.0060183398 1000\na.2 10.77955494 5.32538998 0.0842018089 0.1269148045 1000\na.3 -1.04982353 0.28535997 0.0045119373 0.0049074219 1000\na.4 -1.04898135 0.28129307 0.0044476339 0.0056325117 1000\na.5 -0.74390933 0.26949936 0.0042611590 0.0052178124 1000\na.6  0.21599365 0.26307574 0.0041595927 0.0045153523 1000\na.7  1.81090866 0.39318577 0.0062168129 0.0071483527 1000\n bp  0.83979926 0.26284676 0.0041559722 0.0059795826 1000\nbpC -0.12913322 0.29935741 0.0047332562 0.0049519863 1000\n"

Means of draws

[posterior_β, posterior_α]

2-element Array{Array{Float64,1},1}:
 [1.451237768661092, 0.4233537383099833]                                                                                                       
 [-1.9000576019113276, 10.145744904561283, -2.194854907777116, -2.183382830567691, -1.9048052110036864, -0.971437677726311, 0.5935301675225381]

End of 10/m10.04d.jl

This page was generated using Literate.jl.