TuringModels

Over-dispersed Oceanic tool complexity

This is model m12.6 in Statistical Rethinking Edition 1.

  1. Data
  2. Model
  3. Output
  4. Original output

Data

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

Model

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;

Output

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"

Original output

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