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)
Plot the chain pooled parameters
plot(chns, section=:pooled)
End of m10.04d1.jl
This page was generated using Literate.jl.