m12.6d · DynamicHMCModels.jl

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(θ)

-182.19705519520178

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ℝ₊) )

problem_transformation (generic function with 1 method)

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);

ForwardDiff AD wrapper for TransformedLogDensity of dimension 13, w/ chunk size 7

Tune and sample.

chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000);

(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([1.2203755612727307, 0.2378165738228704, -0.1237179505928725, 0.11971607222613388, -0.09349239119345844, 0.4749656902667122, 0.06597443540569618, -0.3380957213705037, 0.27339373001931766, 0.012068134741907544, 0.29767652250978444, 0.14878470401886668, -1.3702669414946793], -46.819766857637276, 4, DynamicHMC.DoubledTurn, 0.9137705829793278, 15), NUTS_Transition{Array{Float64,1},Float64}([0.5538654425823543, 0.34991508881851197, -0.49811724919833117, -0.10744907741730066, -0.07401932418409736, -0.020166009851455842, -0.13824806488420466, -0.6387847805565272, -0.19764554139418591, -0.7223177768917445, 0.05570455664277531, -0.8339835586651703, -0.799935757050696], -43.63845830705484, 4, DynamicHMC.DoubledTurn, 0.7537659932077398, 15), NUTS_Transition{Array{Float64,1},Float64}([2.0373566501323603, 0.1386912452201454, -0.5012498284149148, 0.022919737763497094, -0.043114416310875064, 0.5849577882986725, 0.20299285868159506, -0.36385722211363236, 0.4266868281815913, 0.16012573053705398, 0.49487985984057575, 0.46953549477064416, -0.41079065128637804], -46.5585312456529, 4, DynamicHMC.DoubledTurn, 0.9559785276158455, 15), NUTS_Transition{Array{Float64,1},Float64}([-0.27996241177503967, 0.41276506148827247, 0.25565801573845864, 0.5623159894141463, 0.13262597038342389, 0.425512418897236, -0.013552135375441976, -0.5032493792020711, 0.391445214643828, -0.24926126717838615, 0.22616933867067426, -0.555979099800734, -1.5313412714705037], -48.990698749099195, 4, DynamicHMC.DoubledTurn, 0.8737723225318734, 15), NUTS_Transition{Array{Float64,1},Float64}([2.7031331837400403, 0.08332294644401146, -0.9839480018480282, -0.3070982477072085, -0.15806428966654823, 0.38620962649378765, -0.03368677829655739, -0.3817598233224514, 0.4493955360871747, -0.18759340058301988, 0.5815769251255268, 0.6445326288799329, -0.6837090072869402], -52.31336608568331, 4, DynamicHMC.DoubledTurn, 0.9059151889332643, 15), NUTS_Transition{Array{Float64,1},Float64}([-0.4621485038404696, 0.42697337726123846, 0.45062358051190593, 0.41405990707363716, 0.17770577471496957, 0.21589737398278072, 0.057515761527746616, -0.28366484875132614, 0.12897098138374055, -0.2229282904141837, 0.08697400272127552, -0.7260967862927867, -1.419248612589498], -49.718737500963634, 4, DynamicHMC.DoubledTurn, 0.9724767547309744, 15), NUTS_Transition{Array{Float64,1},Float64}([-0.9912693310454507, 0.492792766996423, 0.24779551506992706, 0.5674364110070784, 0.15234133282533024, 0.37657495045629935, 0.191072740090469, -0.3220438433146633, 0.3097061807763117, -0.2568696370687405, 0.08522694151335075, -0.7619059264781495, -1.0076496399569563], -50.84869657693219, 4, DynamicHMC.DoubledTurn, 0.9117270806514945, 15), NUTS_Transition{Array{Float64,1},Float64}([1.5026289463317428, 0.19435833769811167, -0.5085518856370949, -0.3483528690347095, -0.01757483989140505, 0.7178034783632348, 0.10777218927575058, -0.01729196660251937, 0.2273349191369352, 0.05080192278415599, 0.47180822450445403, 0.19666503007650826, -1.259667893084798], -47.43491949003435, 4, DynamicHMC.DoubledTurn, 0.9251064753966337, 15), NUTS_Transition{Array{Float64,1},Float64}([0.43896085547054997, 0.34808111871523356, -0.139917302946324, 0.34773805784726375, -0.10663885092365522, -0.024654988553841017, 0.01583168329855785, -0.7460908197143858, 0.2868551889877035, -0.299824985716154, 0.134820060037955, -0.36194179421184114, -0.7251256615746962], -50.815869691547064, 4, DynamicHMC.DoubledTurn, 0.9975770921783201, 15), NUTS_Transition{Array{Float64,1},Float64}([-1.487029079394939, 0.5221684068644677, 0.4022761960187978, 0.739269307757662, 0.1118950126911627, 0.7657977990487655, 0.5528791057027023, -0.7701110915229118, 0.4511947845924382, -0.16014002815734826, 0.27612870448809973, -0.8662197754366091, -0.380634620784727], -56.36596730800085, 4, DynamicHMC.DoubledTurn, 0.9824450253399158, 15)  …  NUTS_Transition{Array{Float64,1},Float64}([0.0017112073194228972, 0.3845543836324173, 0.033548122093712024, 0.1577106621475498, -0.022878225690235766, 0.17486657282085327, 0.17054532894076063, -0.3138693986098904, 0.009136975074122318, -0.3075412419367686, 0.18315727989589758, -0.3724323705777571, -1.4126106991039442], -41.63246213002792, 4, DynamicHMC.DoubledTurn, 0.9993139857418915, 15), NUTS_Transition{Array{Float64,1},Float64}([1.7510494235671414, 0.18119092093391265, -0.4214303713120266, 0.09434593133254607, -0.1592785271578073, 0.2983054568090144, 0.013728383080516632, -0.360397692524499, 0.31309427846556265, -0.091875256140752, 0.24552560451173153, 0.2686498348424631, -1.6057045244578525], -43.02249475882323, 4, DynamicHMC.DoubledTurn, 0.8649985948534377, 15), NUTS_Transition{Array{Float64,1},Float64}([0.9963021520452982, 0.2839745567542142, -0.07966727177794153, 0.1455391000239462, -0.1517927208630831, 0.148581308580176, -0.06954249925328368, -0.3313176985578806, 0.32112634606481055, -0.1813812003746495, 0.10855530407435343, -0.27617343410893413, -1.5617089930711543], -43.42481268892823, 4, DynamicHMC.DoubledTurn, 0.9691629051587682, 15), NUTS_Transition{Array{Float64,1},Float64}([2.1685257638632884, 0.15722555252034318, -0.6846262572403995, -0.015570293993886326, -0.2778026925576971, 0.28426551221358454, 0.1325252693894859, -0.36271289479565755, -0.0001733768980721126, -0.2650006730220744, 0.26179840578287644, 0.11693168582176937, -1.1246241456590247], -42.69969192997442, 4, DynamicHMC.DoubledTurn, 0.9592226794012404, 15), NUTS_Transition{Array{Float64,1},Float64}([-0.14145191217481165, 0.3832831546325597, 0.019424139144858826, 0.18054697732624583, 0.31885746312825014, 0.481353302604173, -0.1181659418123487, -0.36795232606667894, 0.43417930915402203, -0.30631902758412, 0.3880616835421962, -0.4023929706520109, -1.278532289776621], -44.746866698898586, 4, DynamicHMC.DoubledTurn, 0.9853759116477465, 15), NUTS_Transition{Array{Float64,1},Float64}([1.8863139086785945, 0.185691381703495, -0.29673711385312346, 0.12453674367980784, -0.20833807940123475, 0.17072940191618047, 0.32762200119630785, -0.5900387497977263, -0.18825095315548113, -0.004848859015249259, 0.16101842187966173, 0.10126781189636982, -1.2648453791323322], -45.97515668261367, 4, DynamicHMC.DoubledTurn, 0.9951199932783178, 15), NUTS_Transition{Array{Float64,1},Float64}([1.5277156856088179, 0.20886046255425975, -0.2738777950584005, 0.27157309831515486, -0.07010518791882075, 0.2697103234324432, 0.2439572958768483, -0.48743369905940437, -0.10267647542769169, -0.06462996345199831, 0.27335021231562023, 0.06790991696075531, -1.313621177218367], -44.37359962965381, 4, DynamicHMC.DoubledTurn, 0.9997360572158585, 15), NUTS_Transition{Array{Float64,1},Float64}([0.6885123024587521, 0.30915154514503734, 0.033037228358797396, -0.30082197373457636, -0.20693730467896052, 0.0977850186452572, 0.0166112462772203, -0.07500821482969042, 0.1247295765026037, -0.11075291637817458, 0.24896109483953213, -0.30734381411514217, -1.7158915164883295], -44.66135560834853, 4, DynamicHMC.DoubledTurn, 0.9908199428794114, 15), NUTS_Transition{Array{Float64,1},Float64}([2.093640748276437, 0.1430759894215697, -0.4511422438735716, 0.2961169563426361, 0.009189391531225496, 0.6339332334332928, -0.07330559514464403, -0.5215695658615613, 0.1368791696443336, -0.035140002659830215, 0.3545574765346653, 0.2651164687764882, -0.9345784098687602], -46.390642578031745, 4, DynamicHMC.DoubledTurn, 0.8408279594041793, 15), NUTS_Transition{Array{Float64,1},Float64}([2.170059414155504, 0.13927258971937095, -0.1768119528287272, 0.003336835042214714, -0.2342277001424952, 0.016418437062444877, -0.19820851411971382, -0.18900868923613362, -0.007693703492201892, -0.4431141152524658, 0.5566666249319283, 0.376047020138992, -1.051882777761266], -51.20599207856391, 4, DynamicHMC.DoubledTurn, 0.9350746374054926, 15)], NUTS sampler in 13 dimensions
  stepsize (ϵ) ≈ 0.245
  maximum depth = 10
  Gaussian kinetic energy, √diag(M⁻¹): [0.6852245484733178, 0.07825352538288516, 0.2584909523664461, 0.222209124234028, 0.17327859587726244, 0.17634863508500134, 0.1737020276433239, 0.18319716224420812, 0.15566395073422917, 0.1868482759791672, 0.16577290359297922, 0.28212099598652646, 0.4080327398486221]
)

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{NamedTuple{(:β, :α, :σ),Tuple{Array{Float64,1},Array{Float64,1},Float64}},1}:
 (β = [1.2203755612727307, 0.2378165738228704], α = [-0.1237179505928725, 0.11971607222613388, -0.09349239119345844, 0.4749656902667122, 0.06597443540569618, -0.3380957213705037, 0.27339373001931766, 0.012068134741907544, 0.29767652250978444, 0.14878470401886668], σ = 0.2540391369139499)        
 (β = [0.5538654425823543, 0.34991508881851197], α = [-0.49811724919833117, -0.10744907741730066, -0.07401932418409736, -0.020166009851455842, -0.13824806488420466, -0.6387847805565272, -0.19764554139418591, -0.7223177768917445, 0.05570455664277531, -0.8339835586651703], σ = 0.44935783126232953)
 (β = [2.0373566501323603, 0.1386912452201454], α = [-0.5012498284149148, 0.022919737763497094, -0.043114416310875064, 0.5849577882986725, 0.20299285868159506, -0.36385722211363236, 0.4266868281815913, 0.16012573053705398, 0.49487985984057575, 0.46953549477064416], σ = 0.6631257415913456)       
 (β = [-0.27996241177503967, 0.41276506148827247], α = [0.25565801573845864, 0.5623159894141463, 0.13262597038342389, 0.425512418897236, -0.013552135375441976, -0.5032493792020711, 0.391445214643828, -0.24926126717838615, 0.22616933867067426, -0.555979099800734], σ = 0.21624542889086107)        
 (β = [2.7031331837400403, 0.08332294644401146], α = [-0.9839480018480282, -0.3070982477072085, -0.15806428966654823, 0.38620962649378765, -0.03368677829655739, -0.3817598233224514, 0.4493955360871747, -0.18759340058301988, 0.5815769251255268, 0.6445326288799329], σ = 0.5047414266428223)

Extract the parameter posterior means.

posterior_β = mean(posterior[i].β for i in 1:length(posterior))
posterior_α = mean(posterior[i].α for i in 1:length(posterior))
posterior_σ = mean(posterior[i].σ for i in 1:length(posterior))

0.3162411564068822

Effective sample sizes (of untransformed draws)

ess = mapslices(effective_sample_size, get_position_matrix(chain); dims = 1)
ess

1×13 Array{Float64,2}:
 865.901  860.172  1000.0  988.114  …  567.313  554.357  801.263  209.454

NUTS-specific statistics

NUTS_statistics(chain)

Hamiltonian Monte Carlo sample of length 1000
  acceptance rate mean: 0.93, min/25%/median/75%/max: 0.01 0.9 0.97 0.99 1.0
  termination: AdjacentTurn => 8% DoubledTurn => 92%
  depth: 2 => 0% 3 => 6% 4 => 93% 5 => 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{Array{Float64,1},1}:
 [1.1138301195940192, 0.25943353128574054]                                                                                                                                                                           
 [-0.1978887526076462, 0.03617271804735259, -0.044314618373848584, 0.3324531582922331, 0.04626053391099113, -0.3121525657507764, 0.1590760067945609, -0.16418496183383896, 0.27710295869553603, -0.08568801581464212]
 [0.3162411564068822]

End of m12.6d.jl

This page was generated using Literate.jl.