m10.02d · DynamicHMCModels.jl

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

using DynamicHMCModels

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])
first(df, 5)

struct m_10_02d_model{TY <: AbstractVector, TX <: AbstractMatrix}
    "Observations."
    y::TY
    "Covariates"
    X::TX
    "Number of observations"
    N::Int
end

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

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

Instantiate the model with data and inits.

N = size(df, 1)
X = hcat(ones(Int64, N), df[!, :prosoc_left]);
y = df[!, :pulled_left]
p = m_10_02d_model(y, X, N);
θ = (β = [1.0, 2.0],)
p(θ)

-487.6540051035728

Write a function to return properly dimensioned transformation.

problem_transformation(p::m_10_02d_model) =
  as( (β = as(Array, size(p.X, 2)), ) )

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);

Tune and sample.

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

(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([-0.17191950474219544, 0.5649306585708268], -347.9692565959096, 2, DynamicHMC.DoubledTurn, 0.9061181975703501, 3), NUTS_Transition{Array{Float64,1},Float64}([-0.16702295656560337, 0.59457210995657], -347.70503238045376, 1, DynamicHMC.DoubledTurn, 1.0, 1), NUTS_Transition{Array{Float64,1},Float64}([-0.12703120945317609, 0.5785153899671025], -347.0346839234007, 1, DynamicHMC.DoubledTurn, 1.0, 1), NUTS_Transition{Array{Float64,1},Float64}([-0.13261124413868342, 0.5287583247339599], -347.04209513947745, 3, DynamicHMC.DoubledTurn, 0.9738731987575703, 7), NUTS_Transition{Array{Float64,1},Float64}([0.006827837045709759, 0.44817150940813116], -347.4608537956937, 2, DynamicHMC.DoubledTurn, 0.9746047472249431, 3), NUTS_Transition{Array{Float64,1},Float64}([0.12114398119084357, 0.4930433102036059], -345.86443215944223, 2, DynamicHMC.DoubledTurn, 0.9768498487499375, 3), NUTS_Transition{Array{Float64,1},Float64}([0.1027431877571647, 0.27473394932395306], -346.41907576071054, 2, DynamicHMC.DoubledTurn, 0.877421135114389, 3), NUTS_Transition{Array{Float64,1},Float64}([0.11019730525451413, 0.302912645688865], -346.31955977593566, 1, DynamicHMC.DoubledTurn, 1.0, 1), NUTS_Transition{Array{Float64,1},Float64}([0.11333975375325317, 0.57081369675486], -346.3908200788499, 3, DynamicHMC.DoubledTurn, 0.9865729162002049, 7), NUTS_Transition{Array{Float64,1},Float64}([0.04292221934981454, 0.6559306396544702], -345.21337415800735, 2, DynamicHMC.DoubledTurn, 0.9922128363646054, 3)  …  NUTS_Transition{Array{Float64,1},Float64}([0.12808716048776947, 0.4752504175295959], -345.408849208841, 2, DynamicHMC.DoubledTurn, 0.9925057675581502, 3), NUTS_Transition{Array{Float64,1},Float64}([0.17501535468023074, 0.38238422008731476], -345.4982876315726, 3, DynamicHMC.DoubledTurn, 0.9626966920255848, 7), NUTS_Transition{Array{Float64,1},Float64}([0.03633875561997968, 0.5872154010396271], -345.6994461771976, 3, DynamicHMC.DoubledTurn, 0.9763116018466506, 7), NUTS_Transition{Array{Float64,1},Float64}([0.08306932227097359, 0.5969757485666433], -345.1376091146886, 3, DynamicHMC.AdjacentTurn, 0.9650939640952693, 15), NUTS_Transition{Array{Float64,1},Float64}([0.1757995511354754, 0.5397664742242448], -345.7137579404565, 2, DynamicHMC.DoubledTurn, 0.9379627094933953, 3), NUTS_Transition{Array{Float64,1},Float64}([-0.02542833180500391, 0.4937411381255956], -345.8164449663965, 3, DynamicHMC.AdjacentTurn, 0.9991799296940971, 15), NUTS_Transition{Array{Float64,1},Float64}([0.2180848959437626, 0.5270169855887885], -346.44089136129094, 3, DynamicHMC.DoubledTurn, 0.9516169656256057, 7), NUTS_Transition{Array{Float64,1},Float64}([0.1371905758749442, 0.6554361164294937], -346.3939903568805, 3, DynamicHMC.DoubledTurn, 1.0, 7), NUTS_Transition{Array{Float64,1},Float64}([0.15205296028743825, 0.5390172930778476], -346.0670430516673, 2, DynamicHMC.DoubledTurn, 1.0, 3), NUTS_Transition{Array{Float64,1},Float64}([-0.09547579182478785, 0.6806433144307732], -345.7952641878946, 2, DynamicHMC.DoubledTurn, 0.9693028918534483, 3)], NUTS sampler in 2 dimensions
  stepsize (ϵ) ≈ 0.766
  maximum depth = 10
  Gaussian kinetic energy, √diag(M⁻¹): [0.12543267446126483, 0.176071720019298]
)

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}}},1}:
 (β = [-0.17191950474219544, 0.5649306585708268],) 
 (β = [-0.16702295656560337, 0.59457210995657],)   
 (β = [-0.12703120945317609, 0.5785153899671025],) 
 (β = [-0.13261124413868342, 0.5287583247339599],) 
 (β = [0.006827837045709759, 0.44817150940813116],)

Extract the parameter posterior means: β,

posterior_β = mean(first, posterior)

2-element Array{Float64,1}:
 0.04404074787996636
 0.5645292603234728

Effective sample sizes (of untransformed draws)

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

1×2 Array{Float64,2}:
 661.607  731.257

NUTS-specific statistics

NUTS_statistics(chain)

Hamiltonian Monte Carlo sample of length 1000
  acceptance rate mean: 0.94, min/25%/median/75%/max: 0.41 0.91 0.97 1.0 1.0
  termination: AdjacentTurn => 20% DoubledTurn => 80%
  depth: 1 => 15% 2 => 52% 3 => 30% 4 => 3%

CmdStan result

m_10_2s_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.05103234 0.12579086 0.0019889282 0.0035186307 1000
bp 0.55711212 0.18074275 0.0028577937 0.0040160451 1000

Quantiles:
       2.5%        25.0%       50.0%      75.0%      97.5%
 a -0.19755400 -0.029431425 0.05024655 0.12978825 0.30087758
bp  0.20803447  0.433720250 0.55340400 0.67960975 0.91466915
";

"\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.05103234 0.12579086 0.0019889282 0.0035186307 1000\nbp 0.55711212 0.18074275 0.0028577937 0.0040160451 1000\n\nQuantiles:\n       2.5%        25.0%       50.0%      75.0%      97.5%\n a -0.19755400 -0.029431425 0.05024655 0.12978825 0.30087758\nbp  0.20803447  0.433720250 0.55340400 0.67960975 0.91466915\n"

Extract the parameter posterior means: β,

posterior_β = mean(first, posterior)

2-element Array{Float64,1}:
 0.04404074787996636
 0.5645292603234728

End of 10/m10.02d.jl

This page was generated using Literate.jl.