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[!, :society] = 1:10;
df[!, :logpop] = map((x) -> log(x), df[!, :population]);
#df[!, :total_tools] = convert(Vector{Int64}, df[!, :total_tools])
first(df[!, [:total_tools, :logpop, :society]], 5)| total_tools | logpop | society | |
|---|---|---|---|
| Int64 | Float64 | Int64 | |
| 1 | 13 | 7.00307 | 1 |
| 2 | 22 | 7.31322 | 2 |
| 3 | 24 | 8.18869 | 3 |
| 4 | 43 | 8.47449 | 4 |
| 5 | 33 | 8.90924 | 5 |
Define problem data structure
struct m_12_06d{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)(θ)
@unpack y, X, S, N, N_societies = problem # extract the data
@unpack β, α, s = trans(θ) # β : a, bp, α : a_society, s
σ = s[1]^2
ll = 0.0
ll += logpdf(Cauchy(0, 1), σ) # sigma
ll += sum(logpdf.(Normal(0, σ), α)) # α[1:10]
ll += logpdf.(Normal(0, 10), β[1]) # a
ll += logpdf.(Normal(0, 1), β[2]) # bp
ll += sum(
[loglikelihood(Poisson(exp(α[S[i]] + dot(X[i, :], β))), [y[i]]) for i in 1:N]
)
endInstantiate 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];
γ = (β = [1.0, 0.25], α = rand(Normal(0, 1), N_societies), s = [0.2]);
p = m_12_06d(y, X, S, N, N_societies);Main.ex-m12.6d2.m_12_06d{CSV.Column{Int64,Int64},Array{Float64,2},Array{Int64,1}}([13, 22, 24, 43, 33, 19, 40, 28, 55, 71], [1.0 7.003065458786462; 1.0 7.313220387090301; … ; 1.0 9.769956159911606; 1.0 12.524526376648708], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 10, 10)Function convert from a single vector of parms to parks NamedTuple
trans = as((β = as(Array, 2), α = as(Array, 10), s = as(Array, 1)));TransformVariables.TransformTuple{NamedTuple{(:β, :α, :s),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ArrayTransform{TransformVariables.Identity,1}}}}((β = TransformVariables.ArrayTransform{TransformVariables.Identity,1}(TransformVariables.Identity(), (2,)), α = TransformVariables.ArrayTransform{TransformVariables.Identity,1}(TransformVariables.Identity(), (10,)), s = TransformVariables.ArrayTransform{TransformVariables.Identity,1}(TransformVariables.Identity(), (1,))), 13)Define input parameter vector
θ = inverse(trans, γ);
p(θ)-4158.936566437799Maximumaposterior
using Optim
x0 = θ;
lower = vcat([0.0, 0.0], -3ones(10), [0.0]);
upper = vcat([2.0, 1.0], 3ones(10), [5.0]);
ll(x) = -p(x);
inner_optimizer = GradientDescent()
res = optimize(ll, lower, upper, x0, Fminbox(inner_optimizer));
res * Status: failure
* Candidate solution
Minimizer: [1.07e+00, 2.66e-01, -1.96e-13, ...]
Minimum: -1.382377e+02
* Found with
Algorithm: Fminbox with Gradient Descent
Initial Point: [1.00e+00, 2.50e-01, -2.84e-01, ...]
* Convergence measures
|x - x'| = 2.01e-11 ≰ 0.0e+00
|x - x'|/|x'| = 1.83e-11 ≰ 0.0e+00
|f(x) - f(x')| = 8.74e-06 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 6.32e-08 ≰ 0.0e+00
|g(x)| = 2.98e+05 ≰ 1.0e-08
* Work counters
Iterations: 1000
f(x) calls: 933209405
∇f(x) calls: 933209405
Minimum gives MAP estimate:
Optim.minimizer(res)13-element Array{Float64,1}:
1.0698724873889005
0.26584195691346685
-1.9591099444923367e-13
3.7224270107018933e-14
-7.820582057420148e-14
4.5995889939024e-13
2.3506923659964085e-14
-4.4735080528122405e-13
1.8301589273103083e-13
-3.141580522265198e-13
4.448076544246152e-13
-2.0434093379376505e-12
6.723611221204989e-5 Write a function to return properly dimensioned transformation.
problem_transformation(p::m_12_06d) =
as( Vector, length(θ) )problem_transformation (generic function with 1 method)Wrap the problem with a transformation, then use ForwardDiff for the gradient.
P = TransformedLogDensity(problem_transformation(p), p)
∇P = LogDensityRejectErrors(ADgradient(:ForwardDiff, P));
#∇P = ADgradient(:ForwardDiff, P);Tune and sample.
chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 4000);(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([1.1253772383240437, 0.2622669279848795, 0.12220828838470743, -0.1369849997766338, 0.02271722912313025, 0.49454421718002867, 0.15341007499761383, -0.05991801907704529, 0.2384526685667614, -0.32467509184525245, 0.21486652403184658, 0.03631578358564913, -0.3827019937005194], -43.33414511577093, 4, DynamicHMC.AdjacentTurn, 0.8564060825108861, 23), NUTS_Transition{Array{Float64,1},Float64}([1.4259718961173498, 0.22738864681711948, -0.15879319565078936, -0.011617306363702652, -0.09403823223438239, 0.18913989177938156, -0.044102485189420025, -0.37557920610303447, 0.048000652553465244, -0.296828840347033, 0.0797573803924331, -0.0813660084647233, -0.5032575715241459], -44.20218755796961, 3, DynamicHMC.DoubledTurn, 0.9482901259585671, 7), NUTS_Transition{Array{Float64,1},Float64}([1.4896572470500349, 0.232544582476378, -0.13263337832673855, -0.03320904446541051, -0.15565347332116952, 0.12616334869529966, -0.16879530401467072, -0.32439652378280576, 0.13086403156261034, -0.2740054411023427, 0.18129734000097128, -0.12150760812066619, -0.5089619013242039], -39.44932755061596, 3, DynamicHMC.AdjacentTurn, 0.9989414643301266, 15), NUTS_Transition{Array{Float64,1},Float64}([0.2935812564862246, 0.32891430656041887, -0.3570150568598451, 0.2245955372364483, 0.3898509932638105, 0.46697636082484634, 0.11662169252784022, -0.38187818673753215, 0.1028822192158862, 0.05596048996103552, 0.5182588578367701, -0.09269333906589824, -0.5678476488254981], -43.22887883243467, 4, DynamicHMC.DoubledTurn, 0.9550514930629946, 15), NUTS_Transition{Array{Float64,1},Float64}([0.45876599189062295, 0.3214209715834159, -0.49306107270893945, 0.15599430263141906, 0.38610774141389387, 0.39569782733058084, 0.06405276789559705, -0.3305935951425482, 0.16574099685100482, 0.08569572862637582, 0.5136785072378698, -0.16054232211081132, -0.5650065908745067], -41.935824311692514, 4, DynamicHMC.DoubledTurn, 0.9995266410672927, 15), NUTS_Transition{Array{Float64,1},Float64}([2.0525233638339273, 0.13941185461299274, 0.19138814235126764, -0.15841036471911374, -0.10260334274179256, 0.34425864838495224, 0.279459969250341, -0.18356396117062082, 0.27277891823162126, -0.09770492564771191, 0.3598556447929344, 0.49578494832033015, -0.6235806665542389], -46.50777773404141, 4, DynamicHMC.DoubledTurn, 0.9857453968684678, 15), NUTS_Transition{Array{Float64,1},Float64}([2.1347874557143247, 0.15667876360787036, -0.11941139675006868, -0.21270558609649412, -0.44625255670217057, 0.15392410232765041, -0.08796587007249489, -0.3194655185289527, 0.012300825910831356, -0.19650263870047166, 0.2743878079839476, 0.2521199280343992, -0.5750778778119698], -51.21788680542076, 4, DynamicHMC.AdjacentTurn, 0.9825905659121666, 31), NUTS_Transition{Array{Float64,1},Float64}([1.889524126069748, 0.17582666015974255, -0.6418550812520067, -0.31386925034888746, -0.1990094457550085, 0.3020706211248802, -0.18255409749033374, -0.38278398303323036, 0.3383861767309641, -0.3847259771863365, 0.3275207975015407, 0.08612539804444633, -0.7142240922422392], -46.7947235653561, 4, DynamicHMC.DoubledTurn, 0.9024595049863732, 15), NUTS_Transition{Array{Float64,1},Float64}([2.1592487768836675, 0.13248319566944627, -0.5021824920983733, -0.08806064502204565, -0.3450836632053417, 0.3178856406533123, 0.3114492276985155, -0.3767930477659342, 0.42646323303807876, -0.12290222123175112, 0.43198655814702286, 0.35861059380536126, -0.46820239534667785], -44.92037185258332, 4, DynamicHMC.DoubledTurn, 0.9903746993061503, 15), NUTS_Transition{Array{Float64,1},Float64}([1.6997158233183465, 0.19706861060129774, -0.8353373701963367, -0.2657313164746481, -0.2773280392824272, 0.42297718443948884, 0.05666219032468828, -0.3223378544655153, 0.18300274154615104, -0.14509881300487187, 0.26928820512804397, 0.06638446331130997, -0.6402516724251436], -45.86404225840802, 4, DynamicHMC.DoubledTurn, 1.0, 15) … NUTS_Transition{Array{Float64,1},Float64}([1.965491787559912, 0.1708354491277778, -0.18894579078175042, -0.24732527870810156, -0.346327452620221, 0.32047402214816395, -0.15886235157396483, -0.5991244051730145, 0.24650943733953037, -0.08781457245317072, 0.22943389829732735, 0.09521095040734276, -0.7284455975683647], -47.406498275122985, 4, DynamicHMC.DoubledTurn, 0.9876365589166842, 15), NUTS_Transition{Array{Float64,1},Float64}([0.6648436045347013, 0.3129801789746841, -0.15155311107839833, 0.04943632143919062, 0.04184952016801375, 0.12786938802107614, 0.1734797137166275, -0.05864450543853605, 0.14473435574336632, -0.16940574441754713, 0.16531340583553492, -0.2238536661109593, -0.2953925530808919], -43.66921504866631, 4, DynamicHMC.DoubledTurn, 0.9942075243290778, 15), NUTS_Transition{Array{Float64,1},Float64}([0.8540819236811621, 0.2930773264673961, -0.15644951331737772, 0.04481376829282129, -0.09912649117641417, 0.23645788886887492, 0.16648258824827292, 0.007036433475064349, 0.09659456445921956, -0.035312048065214784, 0.26583113256704705, -0.13777605749398597, -0.35753597951241106], -42.73494778930227, 4, DynamicHMC.DoubledTurn, 0.9204090476487221, 15), NUTS_Transition{Array{Float64,1},Float64}([0.9792681281874871, 0.26661604183431703, -0.20807597960371227, 0.017448804000165616, -0.11444165518983061, 0.2553106132404711, 0.18253501747017267, 0.02889585464015832, 0.07846908608516752, 0.009190893327119358, 0.27062655315844264, -0.06597974090165275, -0.3545489708250496], -37.824339520550524, 4, DynamicHMC.DoubledTurn, 0.8837281749700167, 15), NUTS_Transition{Array{Float64,1},Float64}([0.7710594282720857, 0.29272299185958345, -0.01704233590432865, -0.1678213306388679, 0.014697885984390764, 0.07181112720053205, -0.020319268098541357, -0.2637964622119228, 0.15350275055251825, -0.24813723189678422, 0.22691973901962995, -0.20757878514259054, -0.41414935588636576], -40.313630616524705, 4, DynamicHMC.DoubledTurn, 0.8603740873164705, 15), NUTS_Transition{Array{Float64,1},Float64}([1.7656885388773764, 0.20017820212876766, -0.3172017120274748, 0.1369970160125108, -0.15468072246908135, 0.3439750969017568, 0.08660168226901324, -0.2017792265279895, 0.07315069810478717, -0.03966781080004931, 0.09772007646162174, 0.06583599269387361, -0.45004916956257923], -39.05573970699302, 4, DynamicHMC.DoubledTurn, 0.9884476347365314, 15), NUTS_Transition{Array{Float64,1},Float64}([1.0508858065208253, 0.2551912330387846, -0.010036796413919608, -0.062100691479728964, 0.07887388170108041, 0.0048897101635304654, -0.024641216170098047, -0.1917768026942107, 0.14504894176457422, -0.23208707847702692, 0.27447137728194676, 0.07464139201714852, -0.35587378214048393], -41.114239011906285, 4, DynamicHMC.DoubledTurn, 0.9877007198718311, 15), NUTS_Transition{Array{Float64,1},Float64}([0.9293494882177882, 0.2897701609158993, -0.18934156281338343, 0.09435997381564987, -0.052596687392960306, 0.3490142437479168, 0.2041428725898484, -0.14794211354903083, 0.07493647874941099, 0.030228181227733797, 0.14222934161255207, -0.29234044898029227, -0.44384879883469597], -39.59473204192665, 4, DynamicHMC.DoubledTurn, 0.9933523175088667, 15), NUTS_Transition{Array{Float64,1},Float64}([1.0406584420948213, 0.2625456156782788, -0.029365417570141335, 0.005821130383794754, -0.017427077605450536, 0.10997591500663077, -0.12163818413762773, -0.2942768080970383, 0.08071363799452377, -0.24998770633276668, 0.3053937866083294, -0.0034306309830703166, -0.41813011859515087], -40.55664465080242, 4, DynamicHMC.DoubledTurn, 0.9947585242046167, 15), NUTS_Transition{Array{Float64,1},Float64}([1.1195202454699016, 0.25857681623700773, -0.4018349788274355, -0.07981549275926866, -0.03512918521246376, 0.38903625325041175, 0.13455603264645116, -0.2884989189962894, 0.01929844262055091, -0.005344296628497919, 0.09154342174175516, -0.21693147424767545, -0.459971425724632], -38.85804020783523, 4, DynamicHMC.DoubledTurn, 0.9885363575243429, 15)], NUTS sampler in 13 dimensions
stepsize (ϵ) ≈ 0.234
maximum depth = 10
Gaussian kinetic energy, √diag(M⁻¹): [0.6614139767691929, 0.07550157762722404, 0.22732720840801104, 0.20407603410337072, 0.1679879265631661, 0.1909443172837496, 0.1799793741449921, 0.19600597361669927, 0.15798805736263583, 0.18210098250850126, 0.1723505726533973, 0.2645453222227464, 0.10852782149118682]
)We use the transformation to obtain the posterior from the chain.
posterior = TransformVariables.transform.(Ref(problem_transformation(p)),
get_position.(chain));
posterior[1:5]5-element Array{Array{Float64,1},1}:
[1.1253772383240437, 0.2622669279848795, 0.12220828838470743, -0.1369849997766338, 0.02271722912313025, 0.49454421718002867, 0.15341007499761383, -0.05991801907704529, 0.2384526685667614, -0.32467509184525245, 0.21486652403184658, 0.03631578358564913, -0.3827019937005194]
[1.4259718961173498, 0.22738864681711948, -0.15879319565078936, -0.011617306363702652, -0.09403823223438239, 0.18913989177938156, -0.044102485189420025, -0.37557920610303447, 0.048000652553465244, -0.296828840347033, 0.0797573803924331, -0.0813660084647233, -0.5032575715241459]
[1.4896572470500349, 0.232544582476378, -0.13263337832673855, -0.03320904446541051, -0.15565347332116952, 0.12616334869529966, -0.16879530401467072, -0.32439652378280576, 0.13086403156261034, -0.2740054411023427, 0.18129734000097128, -0.12150760812066619, -0.5089619013242039]
[0.2935812564862246, 0.32891430656041887, -0.3570150568598451, 0.2245955372364483, 0.3898509932638105, 0.46697636082484634, 0.11662169252784022, -0.38187818673753215, 0.1028822192158862, 0.05596048996103552, 0.5182588578367701, -0.09269333906589824, -0.5678476488254981]
[0.45876599189062295, 0.3214209715834159, -0.49306107270893945, 0.15599430263141906, 0.38610774141389387, 0.39569782733058084, 0.06405276789559705, -0.3305935951425482, 0.16574099685100482, 0.08569572862637582, 0.5136785072378698, -0.16054232211081132, -0.5650065908745067] Extract the parameter posterior means.
posterior_β = mean(trans(posterior[i]).β for i in 1:length(posterior))
posterior_α = mean(trans(posterior[i]).α for i in 1:length(posterior))
posterior_σ = mean(trans(posterior[i]).s for i in 1:length(posterior))[1]^20.27249165372476997Effective sample sizes (of untransformed draws)
ess = mapslices(effective_sample_size, get_position_matrix(chain); dims = 1)
ess1×13 Array{Float64,2}:
3296.24 3301.37 3971.67 3626.47 … 2827.36 1926.44 3197.35 488.611NUTS-specific statistics
NUTS_statistics(chain)Hamiltonian Monte Carlo sample of length 4000
acceptance rate mean: 0.93, min/25%/median/75%/max: 0.0 0.91 0.97 0.99 1.0
termination: AdjacentDivergent => 0% AdjacentTurn => 7% DoubledTurn => 93%
depth: 0 => 0% 1 => 0% 2 => 1% 3 => 6% 4 => 92% 5 => 1% 6 => 0%
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"Show means
[posterior_β, posterior_α, posterior_σ]3-element Array{Any,1}:
[1.11643807294995, 0.25946232932934726]
[-0.1897505644554755, 0.040701885090751795, -0.04180146915130906, 0.31265832957357764, 0.036828532743372715, -0.29682419939231225, 0.13967380955792097, -0.16330897270956143, 0.264447984204219, -0.08758409001758367]
0.27249165372476997 End of m12.6d1.jl
This page was generated using Literate.jl.