Status

Estimated and checked against the book:

Libraries

library(tidyverse)
library(tidybayes)
library(bayesplot)
library(rstan)
library(patchwork)
options(mc.cores = parallel::detectCores())

Section 9.4: Easy HMC – ulam

rugged <- read.csv('data/rugged.csv', sep = ';') %>%
          mutate(log_gdp = log(rgdppc_2000),
                 log_gdp_std = log_gdp / mean(log_gdp, na.rm = TRUE),
                 rugged_std = rugged / max(rugged, na.rm = TRUE),
                 region = ifelse(cont_africa == 1, 'Africa', 'Not Africa')) %>%
          select(country, log_gdp_std, rugged_std, region) %>%
          drop_na
stan_data <- compose_data(rugged)

stan_program <- '
data {
  int<lower=1> n;        // number of observations
  vector[n] log_gdp_std; // outcome
  vector[n] rugged_std;  // regressor
  int region[n];            // africa indicator
}
parameters {
  real<lower=0> sigma;
  vector[2] a;
  vector[2] b;
}
model {
  vector[n] mu;
  for (i in 1:n) {
    mu[i] = a[region[i]] + b[region[i]] * (rugged_std[i] - 0.215);
  }
  a ~ normal(1, 0.1);
  b ~ normal(0, 0.3);
  sigma ~ exponential(1);
  log_gdp_std ~ normal(mu, sigma);
}
'

m9.1 <- stan(model_code = stan_program, data = stan_data)
m9.1
## Inference for Stan model: b01cd93ce6fd8681789fca5db2e21b98.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##         mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## sigma   0.11    0.00 0.01   0.10   0.11   0.11   0.12   0.12  3988    1
## a[1]    0.89    0.00 0.02   0.86   0.88   0.89   0.91   0.93  3931    1
## a[2]    1.04    0.00 0.01   1.02   1.04   1.04   1.05   1.06  4456    1
## b[1]    0.17    0.00 0.09  -0.02   0.11   0.17   0.23   0.34  3965    1
## b[2]   -0.18    0.00 0.07  -0.32  -0.23  -0.18  -0.13  -0.03  3969    1
## lp__  285.08    0.04 1.65 281.06 284.19 285.42 286.33 287.25  1849    1
## 
## Samples were drawn using NUTS(diag_e) at Sat Aug  1 13:47:34 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
mcmc_trace(m9.1)