m12.6d1
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)
0 250 500 750 1000 -2 0 2 4 a Iteration Sample value -2 0 2 4 6 0.0 0.1 0.2 0.3 0.4 0.5 0.6 a Sample value Density 0 250 500 750 1000 -0.2 0.0 0.2 0.4 0.6 bp Iteration Sample value -0.2 0.0 0.2 0.4 0.6 0 1 2 3 4 5 6 bp Sample value Density 0 250 500 750 1000 0.0 0.2 0.4 0.6 0.8 1.0 1.2 sigma_society Iteration Sample value 0.0 0.3 0.6 0.9 1.2 0 1 2 3 4 sigma_society Sample value Density

Plot the chain pooled parameters

plot(chns, section=:pooled)
0 250 500 750 1000 -1.5 -1.0 -0.5 0.0 0.5 a_society[1] Iteration Sample value -1.5 -1.0 -0.5 0.0 0.5 1.0 0.0 0.5 1.0 1.5 a_society[1] Sample value Density 0 250 500 750 1000 -0.5 0.0 0.5 1.0 a_society[2] Iteration Sample value -1.0 -0.5 0.0 0.5 1.0 0.0 0.5 1.0 1.5 2.0 a_society[2] Sample value Density 0 250 500 750 1000 -0.75 -0.50 -0.25 0.00 0.25 0.50 a_society[3] Iteration Sample value -1.0 -0.5 0.0 0.5 0.0 0.5 1.0 1.5 2.0 a_society[3] Sample value Density 0 250 500 750 1000 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 a_society[4] Iteration Sample value -0.5 0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5 2.0 a_society[4] Sample value Density 0 250 500 750 1000 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 a_society[5] Iteration Sample value -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 0.0 0.5 1.0 1.5 2.0 2.5 a_society[5] Sample value Density 0 250 500 750 1000 -1.00 -0.75 -0.50 -0.25 0.00 0.25 a_society[6] Iteration Sample value -1.0 -0.5 0.0 0.5 0.0 0.5 1.0 1.5 2.0 a_society[6] Sample value Density 0 250 500 750 1000 -0.50 -0.25 0.00 0.25 0.50 0.75 a_society[7] Iteration Sample value -0.5 0.0 0.5 1.0 0.0 0.5 1.0 1.5 2.0 2.5 a_society[7] Sample value Density 0 250 500 750 1000 -0.9 -0.6 -0.3 0.0 0.3 a_society[8] Iteration Sample value -1.0 -0.5 0.0 0.5 0.0 0.5 1.0 1.5 2.0 a_society[8] Sample value Density 0 250 500 750 1000 -0.3 0.0 0.3 0.6 0.9 a_society[9] Iteration Sample value -0.5 0.0 0.5 1.0 0.0 0.5 1.0 1.5 2.0 a_society[9] Sample value Density 0 250 500 750 1000 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 a_society[10] Iteration Sample value -1 0 1 2 0.0 0.5 1.0 1.5 a_society[10] Sample value Density

End of m12.6d.jl

This page was generated using Literate.jl.