Africa

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

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     = 30.71 seconds
Compute duration  = 30.71 seconds
parameters        = σ, βAR, βR, βA, α
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.9124    0.1038     0.0033    0.0041   564.1423    0.9991       18.3718
         βAR   -0.4022    6.7703     0.2141    0.3039   364.0892    0.9994       11.8569
          βR    0.5288    6.7725     0.2142    0.3031   365.8438    0.9994       11.9140
          βA    0.3043    9.7038     0.3069    0.4403   360.1388    0.9991       11.7282
           α    7.0447    9.7034     0.3068    0.4383   361.5474    0.9991       11.7741

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

           σ     0.7276    0.8384    0.9056    0.9736    1.1367
         βAR   -13.7031   -5.0066   -0.4931    4.1581   12.3653
          βR   -12.2815   -4.0114    0.5909    5.1958   13.9103
          βA   -19.2804   -6.0694   -0.1506    6.9126   19.6269
           α   -12.2354    0.4536    7.4986   13.4632   26.5998

using StatsPlots

StatsPlots.plot(chns)

"/home/runner/work/TuringModels.jl/TuringModels.jl/__site/assets/models/africa/code/output/chns.svg"

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