This is model m12.6
in Statistical Rethinking Edition 1.
import CSV
import Random
using DataFrames
using TuringModels: project_root
Random.seed!(1)
path = joinpath(project_root, "data", "Kline.csv")
df = CSV.read(path, DataFrame; delim=';')
df.log_pop = log.(df.population)
df.society = 1:nrow(df)
df
10×7 DataFrame
Row │ culture population contact total_tools mean_TU log_pop society
│ String15 Int64 String7 Int64 Float64 Float64 Int64
─────┼──────────────────────────────────────────────────────────────────────────
1 │ Malekula 1100 low 13 3.2 7.00307 1
2 │ Tikopia 1500 low 22 4.7 7.31322 2
3 │ Santa Cruz 3600 low 24 4.0 8.18869 3
4 │ Yap 4791 high 43 5.0 8.47449 4
5 │ Lau Fiji 7400 high 33 5.0 8.90924 5
6 │ Trobriand 8000 high 19 4.0 8.9872 6
7 │ Chuuk 9200 high 40 3.8 9.12696 7
8 │ Manus 13000 low 28 6.6 9.4727 8
9 │ Tonga 17500 high 55 5.4 9.76996 9
10 │ Hawaii 275000 low 71 6.6 12.5245 10
using Turing
@model function m12_6(total_tools, log_pop, society)
N = length(total_tools)
α ~ Normal(0, 10)
βp ~ Normal(0, 1)
σ_society ~ truncated(Cauchy(0, 1), 0, Inf)
N_society = length(unique(society)) ## 10
α_society ~ filldist(Normal(0, σ_society), N_society)
for i in 1:N
λ = exp(α + α_society[society[i]] + βp*log_pop[i])
total_tools[i] ~ Poisson(λ)
end
end;
chns = sample(
m12_6(df.total_tools, df.log_pop, df.society),
NUTS(0.95),
1000
)
Chains MCMC chain (1000×25×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 15.84 seconds
Compute duration = 15.84 seconds
parameters = α, βp, σ_society, α_society[1], α_society[2], α_society[3], α_society[4], α_society[5], α_society[6], α_society[7], α_society[8], α_society[9], α_society[10]
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ess_per_sec
Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64
α 1.1560 0.7612 0.0241 0.0409 416.4094 0.9991 26.2901
βp 0.2548 0.0852 0.0027 0.0045 415.8743 0.9990 26.2563
σ_society 0.3106 0.1236 0.0039 0.0071 306.0761 1.0006 19.3242
α_society[1] -0.2235 0.2325 0.0074 0.0072 615.6183 0.9995 38.8672
α_society[2] 0.0312 0.2173 0.0069 0.0093 550.3464 0.9992 34.7463
α_society[3] -0.0488 0.1805 0.0057 0.0064 704.1932 0.9991 44.4594
α_society[4] 0.3201 0.1856 0.0059 0.0086 549.3299 0.9993 34.6821
α_society[5] 0.0417 0.1671 0.0053 0.0065 574.2615 0.9995 36.2562
α_society[6] -0.3231 0.2070 0.0065 0.0071 586.2791 0.9992 37.0149
α_society[7] 0.1480 0.1730 0.0055 0.0073 521.0055 1.0005 32.8938
α_society[8] -0.1694 0.1941 0.0061 0.0075 572.8922 0.9992 36.1697
α_society[9] 0.2805 0.1819 0.0058 0.0090 449.0373 0.9994 28.3501
α_society[10] -0.0707 0.3214 0.0102 0.0163 410.3765 0.9990 25.9092
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
α -0.3836 0.6693 1.1643 1.5910 2.7934
βp 0.0708 0.2056 0.2550 0.3089 0.4234
σ_society 0.1286 0.2256 0.2919 0.3764 0.5978
α_society[1] -0.7233 -0.3625 -0.2060 -0.0676 0.1969
α_society[2] -0.3952 -0.1095 0.0301 0.1637 0.4468
α_society[3] -0.4197 -0.1634 -0.0420 0.0685 0.2967
α_society[4] 0.0065 0.1861 0.3065 0.4473 0.6906
α_society[5] -0.2825 -0.0651 0.0444 0.1475 0.3868
α_society[6] -0.7490 -0.4575 -0.3081 -0.1798 0.0525
α_society[7] -0.1516 0.0333 0.1354 0.2629 0.4858
α_society[8] -0.5591 -0.2942 -0.1645 -0.0505 0.2083
α_society[9] -0.0467 0.1561 0.2721 0.3975 0.6590
α_society[10] -0.6826 -0.2476 -0.0741 0.0996 0.6445
using StatsPlots
StatsPlots.plot(chns)
"/home/runner/work/TuringModels.jl/TuringModels.jl/__site/assets/models/over-dispersed-oceanic/code/output/chns.svg"
m12_6rethinking = "
Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
a 1.11 0.75 -0.05 2.24 1256 1
bp 0.26 0.08 0.13 0.38 1276 1
a_society[1] -0.20 0.24 -0.57 0.16 2389 1
a_society[2] 0.04 0.21 -0.29 0.38 2220 1
a_society[3] -0.05 0.19 -0.36 0.25 3018 1
a_society[4] 0.32 0.18 0.01 0.60 2153 1
a_society[5] 0.04 0.18 -0.22 0.33 3196 1
a_society[6] -0.32 0.21 -0.62 0.02 2574 1
a_society[7] 0.14 0.17 -0.13 0.40 2751 1
a_society[8] -0.18 0.19 -0.46 0.12 2952 1
a_society[9] 0.27 0.17 -0.02 0.52 2540 1
a_society[10] -0.10 0.30 -0.52 0.37 1433 1
sigma_society 0.31 0.13 0.11 0.47 1345 1
";