using DynamicHMCModels
ProjDir = rel_path_d("..", "scripts", "12")
df = CSV.read(rel_path( "..", "data", "Kline.csv"), delim=';');
size(df) # Should be 10x5
(10, 5)
New col logpop, set log() for population data
df[:logpop] = map((x) -> log(x), df[:population]);
df[:society] = 1:10;
first(df[[:total_tools, :logpop, :society]], 5)
struct m_12_06d_model{TY <: AbstractVector, TX <: AbstractMatrix,
TS <: AbstractVector}
"Observations (total_tools)."
y::TY
"Covariates (logpop)"
X::TX
"Society"
S::TS
"Number of observations (10)"
N::Int
"Number of societies (also 10)"
N_societies::Int
end
Make the type callable with the parameters as a single argument.
function (problem::m_12_06d_model)(θ)
@unpack y, X, S, N, N_societies = problem # extract the data
@unpack β, α, σ = θ # β : a, bp, α : a_society
ll = 0.0
ll += logpdf(Cauchy(0, 1), σ)
ll += sum(logpdf.(Normal(0, σ), α)) # α[1:10]
ll += logpdf.(Normal(0, 10), β[1]) # a
ll += logpdf.(Normal(0, 1), β[2]) # a
ll += sum(
[loglikelihood(Poisson(exp(α[S[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_societies = length(unique(df[:society]))
X = hcat(ones(Int64, N), df[:logpop]);
S = df[:society]
y = df[:total_tools]
p = m_12_06d_model(y, X, S, N, N_societies);
θ = (β = [1.0, 0.25], α = rand(Normal(0, 1), N_societies), σ = 0.2)
p(θ)
-748.6658051475423
Write a function to return properly dimensioned transformation.
problem_transformation(p::m_12_06d_model) =
as( (β = as(Array, size(p.X, 2)), α = as(Array, p.N_societies), σ = asℝ₊) )
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.
posterior = Vector{Array{NamedTuple{(:β, :α, :σ),Tuple{Array{Float64,1},
Array{Float64,1},Float64}},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.0025 s/step ...done
MCMC, adapting ϵ (25 steps)
0.0033 s/step ...done
MCMC, adapting ϵ (50 steps)
0.0038 s/step ...done
MCMC, adapting ϵ (100 steps)
0.0012 s/step ...done
MCMC, adapting ϵ (200 steps)
0.00091 s/step ...done
MCMC, adapting ϵ (400 steps)
0.00053 s/step ...done
MCMC, adapting ϵ (50 steps)
0.00068 s/step ...done
MCMC (1000 steps)
0.00043 s/step ...done
MCMC, adapting ϵ (75 steps)
0.0037 s/step ...done
MCMC, adapting ϵ (25 steps)
0.0032 s/step ...done
MCMC, adapting ϵ (50 steps)
0.0027 s/step ...done
MCMC, adapting ϵ (100 steps)
0.0013 s/step ...done
MCMC, adapting ϵ (200 steps)
0.00089 s/step ...done
MCMC, adapting ϵ (400 steps)
0.00053 s/step ...done
MCMC, adapting ϵ (50 steps)
0.00061 s/step ...done
MCMC (1000 steps)
0.0004 s/step ...done
MCMC, adapting ϵ (75 steps)
0.0034 s/step ...done
MCMC, adapting ϵ (25 steps)
0.0038 s/step ...done
MCMC, adapting ϵ (50 steps)
0.0032 s/step ...done
MCMC, adapting ϵ (100 steps)
0.0016 s/step ...done
MCMC, adapting ϵ (200 steps)
0.00089 s/step ...done
MCMC, adapting ϵ (400 steps)
0.00058 s/step ...done
MCMC, adapting ϵ (50 steps)
0.00068 s/step ...done
MCMC (1000 steps)
0.00041 s/step ...done
MCMC, adapting ϵ (75 steps)
0.0028 s/step ...done
MCMC, adapting ϵ (25 steps)
0.0036 s/step ...done
MCMC, adapting ϵ (50 steps)
0.0029 s/step ...done
MCMC, adapting ϵ (100 steps)
0.0015 s/step ...done
MCMC, adapting ϵ (200 steps)
0.00093 s/step ...done
MCMC, adapting ϵ (400 steps)
0.00054 s/step ...done
MCMC, adapting ϵ (50 steps)
0.00073 s/step ...done
MCMC (1000 steps)
0.00042 s/step ...done
CmdStan result
m_12_6_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 1.076167468 0.7704872560 0.01218247319 0.0210530022 1000.000000
bp 0.263056273 0.0823415805 0.00130193470 0.0022645077 1000.000000
a_society.1 -0.191723568 0.2421382537 0.00382854195 0.0060563054 1000.000000
a_society.2 0.054569029 0.2278506876 0.00360263570 0.0051693148 1000.000000
a_society.3 -0.035935050 0.1926364647 0.00304584994 0.0039948433 1000.000000
a_society.4 0.334355037 0.1929971201 0.00305155241 0.0063871707 913.029080
a_society.5 0.049747513 0.1801287716 0.00284808595 0.0043631095 1000.000000
a_society.6 -0.311903245 0.2096126337 0.00331426674 0.0053000536 1000.000000
a_society.7 0.148637507 0.1744680594 0.00275858223 0.0047660246 1000.000000
a_society.8 -0.164567976 0.1821341074 0.00287979309 0.0034297298 1000.000000
a_society.9 0.277066965 0.1758237250 0.00278001719 0.0055844175 991.286501
a_society.10 -0.094149204 0.2846206232 0.00450024719 0.0080735022 1000.000000
sigma_society 0.310352849 0.1374834682 0.00217380450 0.0057325226 575.187461
";
"\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 1.076167468 0.7704872560 0.01218247319 0.0210530022 1000.000000\n bp 0.263056273 0.0823415805 0.00130193470 0.0022645077 1000.000000\n a_society.1 -0.191723568 0.2421382537 0.00382854195 0.0060563054 1000.000000\n a_society.2 0.054569029 0.2278506876 0.00360263570 0.0051693148 1000.000000\n a_society.3 -0.035935050 0.1926364647 0.00304584994 0.0039948433 1000.000000\n a_society.4 0.334355037 0.1929971201 0.00305155241 0.0063871707 913.029080\n a_society.5 0.049747513 0.1801287716 0.00284808595 0.0043631095 1000.000000\n a_society.6 -0.311903245 0.2096126337 0.00331426674 0.0053000536 1000.000000\n a_society.7 0.148637507 0.1744680594 0.00275858223 0.0047660246 1000.000000\n a_society.8 -0.164567976 0.1821341074 0.00287979309 0.0034297298 1000.000000\n a_society.9 0.277066965 0.1758237250 0.00278001719 0.0055844175 991.286501\n a_society.10 -0.094149204 0.2846206232 0.00450024719 0.0080735022 1000.000000\nsigma_society 0.310352849 0.1374834682 0.00217380450 0.0057325226 575.187461\n"
Set varable names
parameter_names = ["a", "bp", "sigma_society"]
pooled_parameter_names = ["a_society[$i]" for i in 1:10]
10-element Array{String,1}:
"a_society[1]"
"a_society[2]"
"a_society[3]"
"a_society[4]"
"a_society[5]"
"a_society[6]"
"a_society[7]"
"a_society[8]"
"a_society[9]"
"a_society[10]"
Create a3d
a3d = Array{Float64, 3}(undef, 1000, 13, 4);
for j in 1:4
for i in 1:1000
a3d[i, 1:2, j] = values(posterior[j][i][1])
a3d[i, 3, j] = values(posterior[j][i][3])
a3d[i, 4:13, 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×13×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_society[1], a_society[2], a_society[3], a_society[4], a_society[5], a_society[6], a_society[7], a_society[8], a_society[9], a_society[10]
parameters = a, bp, sigma_society
parameters
Mean SD Naive SE MCSE ESS
a 1.0941 0.7381 0.0117 0.0127 1000
bp 0.2617 0.0804 0.0013 0.0014 1000
sigma_society 0.3070 0.1285 0.0020 0.0029 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 = a, bp, sigma_society
Empirical Posterior Estimates
─────────────────────────────────────────────────
parameters
Mean SD Naive SE MCSE ESS
a 1.0941 0.7381 0.0117 0.0127 1000
bp 0.2617 0.0804 0.0013 0.0014 1000
sigma_society 0.3070 0.1285 0.0020 0.0029 1000
Quantiles
─────────────────────────────────────────────────
parameters
2.5% 25.0% 50.0% 75.0% 97.5%
a -2.1502 0.6666 1.1030 1.5344 5.1436
bp -0.2017 0.2142 0.2601 0.3088 0.6134
sigma_society 0.0361 0.2195 0.2864 0.3715 1.2879
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_society[1], a_society[2], a_society[3], a_society[4], a_society[5], a_society[6], a_society[7], a_society[8], a_society[9], a_society[10]
Empirical Posterior Estimates
────────────────────────────────────────────────────
pooled
Mean SD Naive SE MCSE ESS
a_society[1] -0.1998 0.2358 0.0037 0.0040 1000
a_society[2] 0.0402 0.2204 0.0035 0.0029 1000
a_society[3] -0.0434 0.1945 0.0031 0.0031 1000
a_society[4] 0.3281 0.1821 0.0029 0.0031 1000
a_society[5] 0.0472 0.1726 0.0027 0.0027 1000
a_society[6] -0.3145 0.1989 0.0031 0.0030 1000
a_society[7] 0.1440 0.1758 0.0028 0.0026 1000
a_society[8] -0.1719 0.1834 0.0029 0.0032 1000
a_society[9] 0.2739 0.1731 0.0027 0.0027 1000
a_society[10] -0.0978 0.2907 0.0046 0.0050 1000
Quantiles
────────────────────────────────────────────────────
pooled
2.5% 25.0% 50.0% 75.0% 97.5%
a_society[1] -1.4364 -0.3387 -0.1825 -0.0386 0.7476
a_society[2] -0.8518 -0.0958 0.0373 0.1723 1.2128
a_society[3] -0.8876 -0.1685 -0.0444 0.0841 0.6828
a_society[4] -0.3246 0.2024 0.3217 0.4420 1.3397
a_society[5] -0.6857 -0.0612 0.0407 0.1578 0.7351
a_society[6] -1.1907 -0.4388 -0.3014 -0.1748 0.3130
a_society[7] -0.5555 0.0267 0.1428 0.2522 0.9446
a_society[8] -1.0249 -0.2839 -0.1636 -0.0504 0.4897
a_society[9] -0.3994 0.1533 0.2703 0.3840 1.0694
a_society[10] -1.4620 -0.2587 -0.0874 0.0688 1.7537
Plot the chain parameters
plot(chns)
Plot the chain pooled parameters
plot(chns, section=:pooled)
End of m12.6d.jl
This page was generated using Literate.jl.