m10.04d2 · DynamicHMCModels.jl

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

using DynamicHMCModels, ForwardDiff, Flux, ReverseDiff
gr(size=(400,400))

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[[:actor, :pulled_left, :prosoc_left, :condition]], 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), ) )

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

P = TransformedLogDensity(problem_transformation(p), p)

TransformedLogDensity of dimension 9

For stress testing

stresstest = false

#ad = :Flux
ad = :ForwardDiff
#ad = :ReverseDiff

if stresstest
  ∇P = ADgradient(:ForwardDiff, P);
  LogDensityProblems.stresstest(p, N=1000, scale=1.0)
else
  ∇P = LogDensityRejectErrors(ADgradient(ad, P));
end

LogDensityRejectErrors{InvalidLogDensityException,LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :α),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}},Main.ex-m10.04d2.m_10_04d_model{Array{Int64,1},Array{Int64,2},Array{Int64,1}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :α),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}},Main.ex-m10.04d2.m_10_04d_model{Array{Int64,1},Array{Int64,2},Array{Int64,1}}}},Float64},Float64,9,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :α),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}},Main.ex-m10.04d2.m_10_04d_model{Array{Int64,1},Array{Int64,2},Array{Int64,1}}}},Float64},Float64,9},1}}}}(ForwardDiff AD wrapper for TransformedLogDensity of dimension 9, w/ chunk size 9)

Run N chains

posterior = Vector{Array{NamedTuple{(:β, :α),Tuple{Array{Float64,1},
  Array{Float64,1}}},1}}(undef, 4)

for i in 1:4
  chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000);
  posterior[i] = TransformVariables.transform.(Ref(problem_transformation(p)),
    get_position.(chain));
end

MCMC, adapting ϵ (75 steps)
0.014 s/step ...done
MCMC, adapting ϵ (25 steps)
0.015 s/step ...done
MCMC, adapting ϵ (50 steps)
0.013 s/step ...done
MCMC, adapting ϵ (100 steps)
0.0099 s/step ...done
MCMC, adapting ϵ (200 steps)
step 158 (of 200), 0.0063 s/step
0.0063 s/step ...done
MCMC, adapting ϵ (400 steps)
step 208 (of 400), 0.0049 s/step
0.0044 s/step ...done
MCMC, adapting ϵ (50 steps)
0.0062 s/step ...done
MCMC (1000 steps)
step 215 (of 1000), 0.0047 s/step
step 446 (of 1000), 0.0045 s/step
step 675 (of 1000), 0.0045 s/step
step 892 (of 1000), 0.0045 s/step
0.0045 s/step ...done
MCMC, adapting ϵ (75 steps)
0.012 s/step ...done
MCMC, adapting ϵ (25 steps)
0.011 s/step ...done
MCMC, adapting ϵ (50 steps)
0.011 s/step ...done
MCMC, adapting ϵ (100 steps)
0.0082 s/step ...done
MCMC, adapting ϵ (200 steps)
step 149 (of 200), 0.0067 s/step
0.0063 s/step ...done
MCMC, adapting ϵ (400 steps)
step 187 (of 400), 0.0054 s/step
0.0046 s/step ...done
MCMC, adapting ϵ (50 steps)
0.0056 s/step ...done
MCMC (1000 steps)
step 220 (of 1000), 0.0046 s/step
step 437 (of 1000), 0.0046 s/step
step 650 (of 1000), 0.0046 s/step
step 858 (of 1000), 0.0047 s/step
0.0047 s/step ...done
MCMC, adapting ϵ (75 steps)
0.019 s/step ...done
MCMC, adapting ϵ (25 steps)
0.01 s/step ...done
MCMC, adapting ϵ (50 steps)
0.014 s/step ...done
MCMC, adapting ϵ (100 steps)
0.0063 s/step ...done
MCMC, adapting ϵ (200 steps)
step 155 (of 200), 0.0065 s/step
0.006 s/step ...done
MCMC, adapting ϵ (400 steps)
step 212 (of 400), 0.0047 s/step
0.0045 s/step ...done
MCMC, adapting ϵ (50 steps)
0.0046 s/step ...done
MCMC (1000 steps)
step 209 (of 1000), 0.0048 s/step
step 423 (of 1000), 0.0048 s/step
step 639 (of 1000), 0.0047 s/step
step 863 (of 1000), 0.0047 s/step
0.0047 s/step ...done
MCMC, adapting ϵ (75 steps)
0.017 s/step ...done
MCMC, adapting ϵ (25 steps)
0.0098 s/step ...done
MCMC, adapting ϵ (50 steps)
0.01 s/step ...done
MCMC, adapting ϵ (100 steps)
0.008 s/step ...done
MCMC, adapting ϵ (200 steps)
step 165 (of 200), 0.0061 s/step
0.006 s/step ...done
MCMC, adapting ϵ (400 steps)
step 233 (of 400), 0.0043 s/step
0.0044 s/step ...done
MCMC, adapting ϵ (50 steps)
0.0048 s/step ...done
MCMC (1000 steps)
step 206 (of 1000), 0.0049 s/step
step 420 (of 1000), 0.0048 s/step
step 624 (of 1000), 0.0048 s/step
step 840 (of 1000), 0.0048 s/step
0.0047 s/step ...done

Result rethinking

rethinking = "
      mean   sd  5.5% 94.5% n_eff Rhat
a[1] -0.74 0.27 -1.19 -0.31  2899    1
a[2] 10.77 5.20  4.60 20.45  1916    1
a[3] -1.05 0.28 -1.50 -0.62  3146    1
a[4] -1.05 0.28 -1.50 -0.61  3525    1
a[5] -0.73 0.28 -1.17 -0.28  3637    1
a[6]  0.22 0.27 -0.21  0.67  3496    1
a[7]  1.82 0.41  1.21  2.50  3202    1
bp    0.83 0.27  0.42  1.27  2070    1
bpC  -0.13 0.31 -0.62  0.34  3430    1
";

"\n      mean   sd  5.5% 94.5% n_eff Rhat\na[1] -0.74 0.27 -1.19 -0.31  2899    1\na[2] 10.77 5.20  4.60 20.45  1916    1\na[3] -1.05 0.28 -1.50 -0.62  3146    1\na[4] -1.05 0.28 -1.50 -0.61  3525    1\na[5] -0.73 0.28 -1.17 -0.28  3637    1\na[6]  0.22 0.27 -0.21  0.67  3496    1\na[7]  1.82 0.41  1.21  2.50  3202    1\nbp    0.83 0.27  0.42  1.27  2070    1\nbpC  -0.13 0.31 -0.62  0.34  3430    1\n"

Set varable names

parameter_names = ["bp", "bpC"]
pooled_parameter_names = ["a[$i]" for i in 1:7]

7-element Array{String,1}:
 "a[1]"
 "a[2]"
 "a[3]"
 "a[4]"
 "a[5]"
 "a[6]"
 "a[7]"

Create a3d

a3d = Array{Float64, 3}(undef, 1000, 9, 4);
for j in 1:4
  for i in 1:1000
    a3d[i, 1:2, j] = values(posterior[j][i][1])
    a3d[i, 3:9, j] = values(posterior[j][i][2])
  end
end

chns = MCMCChains.Chains(a3d,
  vcat(parameter_names, pooled_parameter_names),
  Dict(
    :parameters => parameter_names,
    :pooled => pooled_parameter_names
  )
);

Object of type Chains, with data of type 1000×9×4 Array{Float64,3}

Log evidence      = 0.0
Iterations        = 1:1000
Thinning interval = 1
Chains            = 1, 2, 3, 4
Samples per chain = 1000
pooled            = a[1], a[2], a[3], a[4], a[5], a[6], a[7]
parameters        = bp, bpC

parameters
     Mean    SD   Naive SE  MCSE   ESS
 bp 1.4813 3.5252   0.0557 0.0754 1000
bpC 0.4098 0.2437   0.0039 0.0034 1000

Describe the chain

describe(chns)

Log evidence      = 0.0
Iterations        = 1:1000
Thinning interval = 1
Chains            = 1, 2, 3, 4
Samples per chain = 1000
parameters        = bp, bpC

Empirical Posterior Estimates
──────────────────────────────────────────
parameters
     Mean    SD   Naive SE  MCSE   ESS
 bp 1.4813 3.5252   0.0557 0.0754 1000
bpC 0.4098 0.2437   0.0039 0.0034 1000

Quantiles
──────────────────────────────────────────
parameters
      2.5%    25.0%   50.0%  75.0%  97.5%
 bp -10.5565 -0.8979 1.4917 3.8829 12.1537
bpC  -0.4207  0.2435 0.4107 0.5756  1.5195

Describe the chain

describe(chns, section=:pooled)

Log evidence      = 0.0
Iterations        = 1:1000
Thinning interval = 1
Chains            = 1, 2, 3, 4
Samples per chain = 1000
pooled            = a[1], a[2], a[3], a[4], a[5], a[6], a[7]

Empirical Posterior Estimates
─────────────────────────────────────────────
pooled
       Mean    SD   Naive SE  MCSE   ESS
a[1] -1.9205 3.5272   0.0558 0.0750 1000
a[2] 10.0882 6.1479   0.0972 0.1462 1000
a[3] -2.2276 3.5309   0.0558 0.0756 1000
a[4] -2.2246 3.5311   0.0558 0.0761 1000
a[5] -1.9258 3.5353   0.0559 0.0743 1000
a[6] -0.9982 3.5373   0.0559 0.0757 1000
a[7]  0.5566 3.5314   0.0558 0.0771 1000

Quantiles
─────────────────────────────────────────────
pooled
       2.5%    25.0%   50.0%   75.0%   97.5%
a[1] -12.9120 -4.3099 -1.9193  0.5005 10.3960
a[2]  -5.6394  5.6250  9.5843 13.9126 32.8906
a[3] -13.0093 -4.6430 -2.2524  0.1682  9.4858
a[4] -13.3029 -4.5972 -2.2408  0.1639 10.1828
a[5] -12.6149 -4.3398 -1.9261  0.4659 10.4256
a[6] -11.5104 -3.3880 -1.0156  1.4154 11.1344
a[7] -10.7024 -1.8301  0.5457  2.9491 12.8552

Plot the chain parameters

plot(chns)

Plot the chain pooled parameters

plot(chns, section=:pooled)

End of m10.04d1.jl

This page was generated using Literate.jl.