TuringModels

Africa

This is the first Stan model in Statistical Rethinking Edition 1 (page 249). In Edition 2 (page 242)

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

Data

import CSV
import Random

using DataFrames
using Turing
using TuringModels

Random.seed!(1)

data_path = joinpath(TuringModels.project_root, "data", "rugged.csv")
df = CSV.read(data_path, DataFrame)

df.log_gdp = log.(df.rgdppc_2000)
dropmissing!(df)

df = select(df, :log_gdp, :rugged, :cont_africa);

DataFrame df is shown Section df.

Model

@model function model_fn(y, x₁, x₂)
  σ ~ truncated(Cauchy(0, 2), 0, Inf)
  βAR ~ Normal(0, 10)
  βR ~ Normal(0, 10)
  βA ~ Normal(0, 10)
  α ~ Normal(0, 100)

  μ = α .+ βR * x₁ .+ βA * x₂ .+ βAR * x₁ .* x₂
  y ~ MvNormal(μ, σ)
end

model = model_fn(df.log_gdp, df.rugged, df.cont_africa);

Output

chns = sample(model, NUTS(), 1000)
Chains MCMC chain (1000×17×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 24.95 seconds
Compute duration  = 24.95 seconds
parameters        = α, σ, βA, βAR, βR
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

           σ    0.9144    0.1016     0.0032    0.0033   630.1360    0.9994       25.2570
         βAR    0.3836    7.0259     0.2222    0.3365   425.2160    0.9991       17.0434
          βR   -0.2557    7.0306     0.2223    0.3372   425.2452    0.9991       17.0446
          βA    0.3246   10.4599     0.3308    0.6616   250.2059    1.0001       10.0287
           α    7.0292   10.4585     0.3307    0.6620   250.2245    1.0001       10.0294

Quantiles
  parameters       2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol    Float64   Float64   Float64   Float64   Float64

           σ     0.7437    0.8418    0.9029    0.9761    1.1439
         βAR   -13.8296   -4.0967    0.5810    5.1416   14.3395
          βR   -14.3662   -5.0192   -0.4538    4.1932   14.0195
          βA   -20.2170   -6.7925    0.3892    7.7400   20.8784
           α   -13.5815   -0.4123    6.9437   14.1505   27.6747

using StatsPlots

StatsPlots.plot(chns)

Original output

"""
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  9.22360053 0.139119116 0.0021996664 0.0034632816 1000
   bR -0.20196346 0.076106388 0.0012033477 0.0018370185 1000
   bA -1.94430980 0.227080488 0.0035904578 0.0057840746 1000
  bAR  0.39071684 0.131889143 0.0020853505 0.0032749642 1000
sigma  0.95036370 0.052161768 0.0008247500 0.0009204073 1000

Quantiles:
          2.5%       25.0%       50.0%      75.0%        97.5%
    a  8.95307475  9.12719750  9.2237750  9.31974000  9.490234250
   bR -0.35217930 -0.25334425 -0.2012855 -0.15124725 -0.054216855
   bA -2.39010825 -2.09894500 -1.9432550 -1.78643000 -1.513974250
  bAR  0.13496995  0.30095575  0.3916590  0.47887625  0.650244475
sigma  0.85376115  0.91363250  0.9484920  0.98405750  1.058573750
""";

df

df
45×3 DataFrame
 Row │ log_gdp  rugged   cont_africa
     │ Float64  Float64  Int64
─────┼───────────────────────────────
   1 │ 7.49261    0.858            1
   2 │ 6.43238    1.78             1
   3 │ 6.86612    0.141            1
   4 │ 6.90617    0.236            1
   5 │ 8.9493     0.181            1
   6 │ 7.04582    0.197            1
   7 │ 7.36241    0.224            1
   8 │ 7.54046    0.515            1
   9 │ 6.50528    0.443            1
  10 │ 6.86422    0.152            1
  11 │ 7.47886    3.328            1
  12 │ 8.48861    2.367            1
  13 │ 8.59747    0.51             1
  14 │ 8.18822    0.723            1
  15 │ 8.7191     0.218            1
  16 │ 7.54582    0.228            1
  17 │ 7.58873    0.74             1
  18 │ 7.45483    0.353            1
  19 │ 6.66246    0.491            1
  20 │ 9.62836    0.559            1
  21 │ 6.92631    0.669            1
  22 │ 7.66037    6.202            1
  23 │ 8.16446    2.413            1
  24 │ 6.71568    1.169            1
  25 │ 6.65878    0.147            1
  26 │ 6.77636    0.612            1
  27 │ 7.45542    0.115            1
  28 │ 9.17186    0.949            1
  29 │ 6.37375    1.027            1
  30 │ 8.70911    0.913            1
  31 │ 6.55556    0.178            1
  32 │ 6.78286    0.312            1
  33 │ 6.94618    3.309            1
  34 │ 7.42814    0.442            1
  35 │ 7.26898    0.244            1
  36 │ 6.14557    0.498            1
  37 │ 8.45136    3.063            1
  38 │ 6.71071    0.419            1
  39 │ 7.27021    0.28             1
  40 │ 8.74058    0.726            1
  41 │ 6.25734    0.677            1
  42 │ 7.1257     0.913            1
  43 │ 9.1505     1.761            1
  44 │ 6.65158    0.533            1
  45 │ 7.82373    1.194            1