This is the first Stan model in Statistical Rethinking Edition 1 (page 249). In Edition 2 (page 242)
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 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);
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"
"""
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
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