m10.04d1

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.04d1.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.04d1.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.04d1.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 single chains

chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 3000);
posterior = TransformVariables.transform.(Ref(problem_transformation(p)),
  get_position.(chain));

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
";
MCMC, adapting ϵ (75 steps)
0.016 s/step ...done
MCMC, adapting ϵ (25 steps)
0.011 s/step ...done
MCMC, adapting ϵ (50 steps)
0.0079 s/step ...done
MCMC, adapting ϵ (100 steps)
0.0086 s/step ...done
MCMC, adapting ϵ (200 steps)
step 159 (of 200), 0.0063 s/step
0.006 s/step ...done
MCMC, adapting ϵ (400 steps)
step 218 (of 400), 0.0046 s/step
0.0042 s/step ...done
MCMC, adapting ϵ (50 steps)
0.0041 s/step ...done
MCMC (3000 steps)
step 238 (of 3000), 0.0042 s/step
step 472 (of 3000), 0.0042 s/step
step 711 (of 3000), 0.0042 s/step
step 968 (of 3000), 0.0041 s/step
step 1216 (of 3000), 0.0041 s/step
step 1455 (of 3000), 0.0041 s/step
step 1704 (of 3000), 0.0041 s/step
step 1949 (of 3000), 0.0041 s/step
step 2206 (of 3000), 0.0041 s/step
step 2449 (of 3000), 0.0041 s/step
step 2699 (of 3000), 0.0041 s/step
step 2941 (of 3000), 0.0041 s/step
0.0041 s/step ...done
"\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, this will be automated using θ

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, 3000, 9, 1);
for i in 1:3000
  a3d[i, 1:2, 1] = values(posterior[i][1])
  a3d[i, 3:9, 1] = values(posterior[i][2])
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 3000×9×1 Array{Float64,3}

Log evidence      = 0.0
Iterations        = 1:3000
Thinning interval = 1
Chains            = 1
Samples per chain = 3000
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.3766 3.5989   0.0657 0.0955 1419.1074
bpC 0.4173 0.2460   0.0045 0.0032 3000.0000

Describe the chain

describe(chns)
Log evidence      = 0.0
Iterations        = 1:3000
Thinning interval = 1
Chains            = 1
Samples per chain = 3000
parameters        = bp, bpC

Empirical Posterior Estimates
───────────────────────────────────────────
parameters
     Mean    SD   Naive SE  MCSE     ESS
 bp 1.3766 3.5989   0.0657 0.0955 1419.1074
bpC 0.4173 0.2460   0.0045 0.0032 3000.0000

Quantiles
───────────────────────────────────────────
parameters
      2.5%    25.0%   50.0%  75.0%  97.5%
 bp -11.4096 -1.2158 1.3415 3.8854 13.5759
bpC  -0.3573  0.2553 0.4150 0.5822  1.2518

Describe the chain

describe(chns, section=:pooled)
Log evidence      = 0.0
Iterations        = 1:3000
Thinning interval = 1
Chains            = 1
Samples per chain = 3000
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.8183 3.6074   0.0659 0.0944 1460.0257
a[2] 10.5963 6.2122   0.1134 0.2258  756.5888
a[3] -2.1281 3.6077   0.0659 0.0965 1397.2852
a[4] -2.1300 3.6071   0.0659 0.0957 1422.0506
a[5] -1.8223 3.6090   0.0659 0.0964 1401.1982
a[6] -0.9018 3.6051   0.0658 0.0947 1449.3533
a[7]  0.6846 3.6236   0.0662 0.0983 1358.3387

Quantiles
─────────────────────────────────────────────
pooled
       2.5%    25.0%   50.0%   75.0%   97.5%
a[1] -14.1852 -4.3162 -1.8036  0.8123 11.4490
a[2]  -5.6534  6.3179  9.9603 14.2412 32.8029
a[3] -14.2744 -4.6472 -2.0859  0.4351 10.8671
a[4] -14.1196 -4.6234 -2.0990  0.4991 10.8758
a[5] -13.8834 -4.3306 -1.7750  0.7966 11.0854
a[6] -13.3491 -3.4026 -0.8532  1.6987 11.7333
a[7] -11.5728 -1.7685  0.7161  3.3611 14.1383

Plot the chain parameters

plot(chns)
0 1000 2000 3000 -10 -5 0 5 10 bp Iteration Sample value -10 -5 0 5 10 15 0.000 0.025 0.050 0.075 0.100 bp Sample value Density 0 1000 2000 3000 -0.25 0.00 0.25 0.50 0.75 1.00 1.25 bpC Iteration Sample value -0.5 0.0 0.5 1.0 0.0 0.5 1.0 1.5 bpC Sample value Density

Plot the chain pooled parameters

plot(chns, section=:pooled)
0 1000 2000 3000 -10 -5 0 5 10 a[1] Iteration Sample value -15 -10 -5 0 5 10 0.000 0.025 0.050 0.075 0.100 a[1] Sample value Density 0 1000 2000 3000 0 10 20 30 a[2] Iteration Sample value -10 0 10 20 30 0.00 0.02 0.04 0.06 a[2] Sample value Density 0 1000 2000 3000 -15 -10 -5 0 5 10 a[3] Iteration Sample value -15 -10 -5 0 5 10 0.000 0.025 0.050 0.075 0.100 a[3] Sample value Density 0 1000 2000 3000 -10 -5 0 5 10 a[4] Iteration Sample value -15 -10 -5 0 5 10 0.000 0.025 0.050 0.075 0.100 a[4] Sample value Density 0 1000 2000 3000 -10 -5 0 5 10 a[5] Iteration Sample value -15 -10 -5 0 5 10 0.000 0.025 0.050 0.075 0.100 a[5] Sample value Density 0 1000 2000 3000 -10 -5 0 5 10 a[6] Iteration Sample value -15 -10 -5 0 5 10 15 0.000 0.025 0.050 0.075 0.100 a[6] Sample value Density 0 1000 2000 3000 -10 -5 0 5 10 a[7] Iteration Sample value -15 -10 -5 0 5 10 15 0.000 0.025 0.050 0.075 0.100 a[7] Sample value Density

End of m10.04d1.jl

This page was generated using Literate.jl.