Statistical Rethinking (Code)
Chapter 8
March, 2017
Code from Statistical Rethinking modified by R Pruim is shown below. Differences to the oringal include:
- a preference for putting data into containers (data frames, mostly), rather than working with lose vectors.
- use of
ggplot2(viaggformula) rather than base graphics - use of
tidyversefor data transformation - better (in my opinion) naming conventions
R Code 8.1
The book does this as a chunk of code. It is much nicer to wrap this in a function and to return a data frame that makes plotting nicer.
KingMarkov <- function(
num_steps = 1e5,
population = 1:10,
start = 10) {
num_islands <- length(population)
position <- rep(NA, num_steps) # trick to pre-alocate memory
current <- start
for (i in 1:num_steps) {
# record current position
position[i] <- current
# flip coin to generate proposal
proposal <- 1 + (current + sample(c(-2, 0), size = 1) ) %% 10
# move?
prob_move <- population[proposal] / population[current]
current <- ifelse(runif(1) < prob_move, proposal, current)
}
data_frame(
step = 1:num_steps,
position = position
)
}In the short-run, any particular path looks like a jumble, and the desired probabilities are not acheived.
set.seed(100)
KM <- KingMarkov(100)
gf_point(position ~ step, data = KM)
tally( ~ position, data = KM)## position
## 4 5 6 7 8 9 10
## 1 3 9 28 28 15 16
With a larger number of steps, things are still not perfect, but the are much closer to the target probabilities.
set.seed(123)
KM <- KingMarkov(1e5)
gf_histogram( ~ position, data = KM, binwidth = 1)
KM %>%
group_by(position) %>%
summarise(n = n(), prop = n / 1e5) %>%
mutate(target = (1:10)/(sum(1:10)))## # A tibble: 10 × 4
## position n prop target
## <dbl> <int> <dbl> <dbl>
## 1 1 1813 0.01813 0.01818182
## 2 2 3658 0.03658 0.03636364
## 3 3 5316 0.05316 0.05454545
## 4 4 7193 0.07193 0.07272727
## 5 5 9118 0.09118 0.09090909
## 6 6 11030 0.11030 0.10909091
## 7 7 12889 0.12889 0.12727273
## 8 8 14772 0.14772 0.14545455
## 9 9 16426 0.16426 0.16363636
## 10 10 17785 0.17785 0.18181818
Alternatively, we can propose any island instead of just neighboring islands.
KingMarkov2 <- function(
num_steps = 1e5,
population = 1:10,
start = 10) {
num_islands <- length(population)
position <- rep(NA, num_steps) # trick to pre-alocate memory
current <- start
for (i in 1:num_steps) {
# record current position
position[i] <- current
# propose any one of the other islands
proposal <- sample(setdiff(1:num_islands, current), 1)
# move?
prob_move <- population[proposal] / population[current]
current <- ifelse(runif(1) < prob_move, proposal, current)
}
data_frame(
step = 1:num_steps,
position = position
)
}The mixing is faster now.
set.seed(100)
KM2 <- KingMarkov2(100)
gf_point(position ~ step, data = KM2)
tally( ~ position, data = KM2)## position
## 1 2 3 4 5 6 7 8 9 10
## 1 4 5 7 12 12 9 15 11 24
With a larger number of steps, things are still not perfect, but the are much closer to the target probabilities.
set.seed(123)
KM2 <- KingMarkov2(1e5)
gf_histogram( ~ position, data = KM2, binwidth = 1)
KM2 %>%
group_by(position) %>%
summarise(n = n(), prop = n / 1e5) %>%
mutate(target = (1:10)/(sum(1:10)))## # A tibble: 10 × 4
## position n prop target
## <dbl> <int> <dbl> <dbl>
## 1 1 1798 0.01798 0.01818182
## 2 2 3578 0.03578 0.03636364
## 3 3 5587 0.05587 0.05454545
## 4 4 7158 0.07158 0.07272727
## 5 5 9146 0.09146 0.09090909
## 6 6 10768 0.10768 0.10909091
## 7 7 12917 0.12917 0.12727273
## 8 8 14467 0.14467 0.14545455
## 9 9 16315 0.16315 0.16363636
## 10 10 18266 0.18266 0.18181818
R Code 8.2
library(rethinking)
data(rugged)
Nations <- rugged %>%
mutate(log_gdp = log(rgdppc_2000)) %>%
filter(!is.na(rgdppc_2000))R Code 8.3
m8.1 <- map(
alist(
log_gdp ~ dnorm(mu, sigma),
mu <- a + bR * rugged + bA * cont_africa + bAR * rugged * cont_africa,
a ~ dnorm(0, 100),
bR ~ dnorm(0, 10),
bA ~ dnorm(0, 10),
bAR ~ dnorm(0, 10),
sigma ~ dunif(0, 10)
),
data = Nations
)
precis(m8.1)## Mean StdDev 5.5% 94.5%
## a 9.22 0.14 9.00 9.44
## bR -0.20 0.08 -0.32 -0.08
## bA -1.95 0.22 -2.31 -1.59
## bAR 0.39 0.13 0.19 0.60
## sigma 0.93 0.05 0.85 1.01R Code 8.4
Restrict to 3 variables
inspect(Nations) # Note: several variables have missing values##
## categorical variables:
## name class levels n missing
## 1 isocode factor 234 170 0
## 2 country factor 234 170 0
## distribution
## 1 AGO (0.6%), ALB (0.6%), ARE (0.6%) ...
## 2 Albania (0.6%), Algeria (0.6%) ...
##
## quantitative variables:
## name class min Q1 median
## 1 isonum integer 8.000000 212.500000 417.500000
## 2 rugged numeric 0.003000 0.442250 0.979500
## 3 rugged_popw numeric 0.003000 0.283000 0.579500
## 4 rugged_slope numeric 0.005000 1.287750 2.792500
## 5 rugged_lsd numeric 0.001000 0.118500 0.284500
## 6 rugged_pc numeric 0.000000 1.954250 12.117500
## 7 land_area numeric 3.000000 2800.500000 14148.000000
## 8 lat numeric -41.806000 2.225000 17.238000
## 9 lon numeric -174.847000 -11.578000 19.452500
## 10 soil numeric 0.000000 18.617250 40.438000
## 11 desert numeric 0.000000 0.000000 0.000000
## 12 tropical numeric 0.000000 0.000000 1.617000
## 13 dist_coast numeric 0.000000 0.039000 0.164500
## 14 near_coast numeric 0.000000 8.887290 33.900950
## 15 gemstones integer 0.000000 0.000000 0.000000
## 16 rgdppc_2000 numeric 466.647000 1880.829250 5314.742000
## 17 rgdppc_1950_m numeric 289.153000 660.201500 1410.050500
## 18 rgdppc_1975_m numeric 512.136000 1047.390500 2518.934000
## 19 rgdppc_2000_m numeric 216.950000 1352.019000 3637.062000
## 20 rgdppc_1950_2000_m numeric 469.724000 1007.087250 2423.392000
## 21 q_rule_law numeric -2.050000 -0.807000 -0.203000
## 22 cont_africa integer 0.000000 0.000000 0.000000
## 23 cont_asia integer 0.000000 0.000000 0.000000
## 24 cont_europe integer 0.000000 0.000000 0.000000
## 25 cont_oceania integer 0.000000 0.000000 0.000000
## 26 cont_north_america integer 0.000000 0.000000 0.000000
## 27 cont_south_america integer 0.000000 0.000000 0.000000
## 28 legor_gbr integer 0.000000 0.000000 0.000000
## 29 legor_fra integer 0.000000 0.000000 0.000000
## 30 legor_soc integer 0.000000 0.000000 0.000000
## 31 legor_deu integer 0.000000 0.000000 0.000000
## 32 legor_sca integer 0.000000 0.000000 0.000000
## 33 colony_esp integer 0.000000 0.000000 0.000000
## 34 colony_gbr integer 0.000000 0.000000 0.000000
## 35 colony_fra integer 0.000000 0.000000 0.000000
## 36 colony_prt integer 0.000000 0.000000 0.000000
## 37 colony_oeu integer 0.000000 0.000000 0.000000
## 38 africa_region_n integer 0.000000 0.000000 0.000000
## 39 africa_region_s integer 0.000000 0.000000 0.000000
## 40 africa_region_w integer 0.000000 0.000000 0.000000
## 41 africa_region_e integer 0.000000 0.000000 0.000000
## 42 africa_region_c integer 0.000000 0.000000 0.000000
## 43 slave_exports numeric 0.000000 0.000000 0.000000
## 44 dist_slavemkt_atlantic numeric 3.647000 5.121000 5.712000
## 45 dist_slavemkt_indian numeric 0.032000 2.682000 7.643000
## 46 dist_slavemkt_saharan numeric 0.310000 2.642000 3.353000
## 47 dist_slavemkt_redsea numeric 0.064000 2.293000 3.515000
## 48 pop_1400 integer 0.000000 78895.750000 348220.000000
## 49 european_descent numeric 0.000000 0.000000 3.450000
## 50 log_gdp numeric 6.145573 7.539468 8.578233
## Q3 max mean sd n missing
## 1 6.452500e+02 8.940000e+02 4.293941e+02 2.550246e+02 170 0
## 2 1.957250e+00 6.202000e+00 1.333182e+00 1.168467e+00 170 0
## 3 9.882500e-01 4.165000e+00 7.352353e-01 5.944169e-01 170 0
## 4 5.575500e+00 1.759500e+01 3.813453e+00 3.382790e+00 170 0
## 5 5.565000e-01 1.559000e+00 3.793471e-01 3.280890e-01 170 0
## 6 3.133000e+01 9.024900e+01 1.878796e+01 2.026075e+01 170 0
## 7 5.445675e+04 1.638134e+06 7.325901e+04 1.954638e+05 170 0
## 8 4.007150e+01 6.499000e+01 1.856515e+01 2.503034e+01 170 0
## 9 4.484400e+01 1.714780e+02 1.517724e+01 6.562387e+01 170 0
## 10 5.623850e+01 1.000000e+02 3.956522e+01 2.557699e+01 170 0
## 11 0.000000e+00 7.728000e+01 3.423624e+00 1.108766e+01 170 0
## 12 9.770950e+01 1.000000e+02 4.054217e+01 4.534055e+01 170 0
## 13 4.085000e-01 2.206000e+00 3.180941e-01 4.177448e-01 170 0
## 14 9.612520e+01 1.000000e+02 4.560283e+01 4.030518e+01 170 0
## 15 0.000000e+00 2.641540e+05 6.513724e+03 3.292907e+04 170 0
## 16 1.310015e+04 5.779209e+04 9.094893e+03 9.699991e+03 170 0
## 17 2.413561e+03 2.887814e+04 2.392489e+03 3.370413e+03 124 46
## 18 6.451961e+03 2.546500e+04 4.630049e+03 4.825339e+03 124 46
## 19 8.219698e+03 2.844885e+04 6.675195e+03 7.049194e+03 145 25
## 20 6.096618e+03 2.011999e+04 4.466931e+03 4.632625e+03 124 46
## 21 6.430000e-01 2.023000e+00 -3.151479e-02 9.554936e-01 169 1
## 22 1.000000e+00 1.000000e+00 2.882353e-01 4.542793e-01 170 0
## 23 0.000000e+00 1.000000e+00 2.352941e-01 4.254356e-01 170 0
## 24 0.000000e+00 1.000000e+00 2.176471e-01 4.138652e-01 170 0
## 25 0.000000e+00 1.000000e+00 5.882353e-02 2.359892e-01 170 0
## 26 0.000000e+00 1.000000e+00 1.352941e-01 3.430479e-01 170 0
## 27 0.000000e+00 1.000000e+00 6.470588e-02 2.467329e-01 170 0
## 28 1.000000e+00 1.000000e+00 3.176471e-01 4.669368e-01 170 0
## 29 1.000000e+00 1.000000e+00 4.411765e-01 4.979946e-01 170 0
## 30 0.000000e+00 1.000000e+00 1.823529e-01 3.872759e-01 170 0
## 31 0.000000e+00 1.000000e+00 2.941176e-02 1.694569e-01 170 0
## 32 0.000000e+00 1.000000e+00 2.941176e-02 1.694569e-01 170 0
## 33 0.000000e+00 1.000000e+00 1.176471e-01 3.231416e-01 170 0
## 34 1.000000e+00 1.000000e+00 3.352941e-01 4.734878e-01 170 0
## 35 0.000000e+00 1.000000e+00 1.647059e-01 3.720107e-01 170 0
## 36 0.000000e+00 1.000000e+00 3.529412e-02 1.850673e-01 170 0
## 37 0.000000e+00 1.000000e+00 2.352941e-02 1.520254e-01 170 0
## 38 0.000000e+00 1.000000e+00 2.352941e-02 1.520254e-01 170 0
## 39 0.000000e+00 1.000000e+00 5.294118e-02 2.245776e-01 170 0
## 40 0.000000e+00 1.000000e+00 8.823529e-02 2.844747e-01 170 0
## 41 0.000000e+00 1.000000e+00 7.647059e-02 2.665347e-01 170 0
## 42 0.000000e+00 1.000000e+00 4.705882e-02 2.123903e-01 170 0
## 43 0.000000e+00 3.610000e+06 9.189638e+04 3.746049e+05 170 0
## 44 1.013100e+01 1.639300e+01 7.496224e+00 3.349753e+00 49 121
## 45 9.300000e+00 1.583300e+01 6.660061e+00 4.079243e+00 49 121
## 46 4.821000e+00 6.637000e+00 3.539510e+00 1.578397e+00 49 121
## 47 4.423000e+00 6.465000e+00 3.412878e+00 1.502717e+00 49 121
## 48 1.191290e+06 8.094376e+07 2.037433e+06 8.723193e+06 168 2
## 49 8.990525e+01 1.000000e+02 3.405936e+01 4.218259e+01 152 18
## 50 9.480244e+00 1.096461e+01 8.517117e+00 1.166495e+00 170 0Nations.trim <- Nations %>% select(log_gdp, rugged, cont_africa)
inspect(Nations.trim) # No missing values among these three variables##
## quantitative variables:
## name class min Q1 median Q3 max
## 1 log_gdp numeric 6.145573 7.539468 8.578233 9.480244 10.96461
## 2 rugged numeric 0.003000 0.442250 0.979500 1.957250 6.20200
## 3 cont_africa integer 0.000000 0.000000 0.000000 1.000000 1.00000
## mean sd n missing
## 1 8.5171175 1.1664946 170 0
## 2 1.3331824 1.1684674 170 0
## 3 0.2882353 0.4542793 170 0R Code 8.4
Fit with Stan:
m8.1stan <- map2stan(
alist(
log_gdp ~ dnorm(mu, sigma),
mu <- a + bR * rugged + bA * cont_africa + bAR * rugged * cont_africa,
a ~ dnorm(0, 100),
bR ~ dnorm(0, 10),
bA ~ dnorm(0, 10),
bAR ~ dnorm(0, 10),
sigma ~ dcauchy(0, 2)
),
data = Nations.trim,
refresh = 0
)R Code 8.6
show(m8.1stan)## map2stan model fit
## 1000 samples from 1 chain
##
## Formula:
## log_gdp ~ dnorm(mu, sigma)
## mu <- a + bR * rugged + bA * cont_africa + bAR * rugged * cont_africa
## a ~ dnorm(0, 100)
## bR ~ dnorm(0, 10)
## bA ~ dnorm(0, 10)
## bAR ~ dnorm(0, 10)
## sigma ~ dcauchy(0, 2)
##
## Log-likelihood at expected values: -229.45
## Deviance: 458.9
## DIC: 468.52
## Effective number of parameters (pD): 4.81
##
## WAIC (SE): 468.9 (14.8)
## pWAIC: 4.87precis(m8.1stan)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 9.24 0.14 9.02 9.45 376 1
## bR -0.21 0.08 -0.32 -0.08 449 1
## bA -1.97 0.21 -2.29 -1.62 493 1
## bAR 0.41 0.13 0.22 0.62 493 1
## sigma 0.95 0.05 0.86 1.03 468 1plot(m8.1stan)
R Code 8.7
m8.1stan_4chains <- map2stan(m8.1stan, chains = 4, cores = 4)##
## SAMPLING FOR MODEL 'log_gdp ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 6.6e-05 seconds (Sampling)
## 7e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]precis(m8.1stan_4chains)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 9.22 0.14 9.00 9.45 1762 1
## bR -0.20 0.08 -0.32 -0.07 1678 1
## bA -1.95 0.23 -2.35 -1.61 1796 1
## bAR 0.39 0.13 0.19 0.60 1803 1
## sigma 0.95 0.05 0.86 1.03 3032 1plot(m8.1stan_4chains)
R Code 8.8
m8.1s.post <- extract.samples(m8.1stan)
str(m8.1s.post)## List of 5
## $ a : num [1:1000(1d)] 9.3 9.17 9.37 9.05 9.21 ...
## $ bR : num [1:1000(1d)] -0.224 -0.167 -0.307 -0.113 -0.198 ...
## $ bA : num [1:1000(1d)] -1.98 -2.04 -2.21 -1.63 -2.09 ...
## $ bAR : num [1:1000(1d)] 0.445 0.363 0.586 0.193 0.516 ...
## $ sigma: num [1:1000(1d)] 0.959 0.881 0.92 0.912 0.961 ...str(m8.1s.post %>% as.data.frame())## 'data.frame': 1000 obs. of 5 variables:
## $ a : num [1:1000(1d)] 9.3 9.17 9.37 9.05 9.21 ...
## $ bR : num [1:1000(1d)] -0.224 -0.167 -0.307 -0.113 -0.198 ...
## $ bA : num [1:1000(1d)] -1.98 -2.04 -2.21 -1.63 -2.09 ...
## $ bAR : num [1:1000(1d)] 0.445 0.363 0.586 0.193 0.516 ...
## $ sigma: num [1:1000(1d)] 0.959 0.881 0.92 0.912 0.961 ...gf_dens(~ bR, data = m8.1s.post %>% as.data.frame())
R Code 8.9
The rethinking package includes methods for pairs().
pairs(m8.1s.post)
R Code 8.10
pairs(m8.1stan)
R Code 8.11
show(m8.1stan)## map2stan model fit
## 1000 samples from 1 chain
##
## Formula:
## log_gdp ~ dnorm(mu, sigma)
## mu <- a + bR * rugged + bA * cont_africa + bAR * rugged * cont_africa
## a ~ dnorm(0, 100)
## bR ~ dnorm(0, 10)
## bA ~ dnorm(0, 10)
## bAR ~ dnorm(0, 10)
## sigma ~ dcauchy(0, 2)
##
## Log-likelihood at expected values: -229.45
## Deviance: 458.9
## DIC: 468.52
## Effective number of parameters (pD): 4.81
##
## WAIC (SE): 468.9 (14.8)
## pWAIC: 4.87R Code 8.12
plot(m8.1stan)
Want to see raw Stan code? Here it is.
stancode(m8.1stan)## data{
## int<lower=1> N;
## real log_gdp[N];
## real rugged[N];
## int cont_africa[N];
## }
## parameters{
## real a;
## real bR;
## real bA;
## real bAR;
## real<lower=0> sigma;
## }
## model{
## vector[N] mu;
## sigma ~ cauchy( 0 , 2 );
## bAR ~ normal( 0 , 10 );
## bA ~ normal( 0 , 10 );
## bR ~ normal( 0 , 10 );
## a ~ normal( 0 , 100 );
## for ( i in 1:N ) {
## mu[i] = a + bR * rugged[i] + bA * cont_africa[i] + bAR * rugged[i] * cont_africa[i];
## }
## log_gdp ~ normal( mu , sigma );
## }
## generated quantities{
## vector[N] mu;
## real dev;
## dev = 0;
## for ( i in 1:N ) {
## mu[i] = a + bR * rugged[i] + bA * cont_africa[i] + bAR * rugged[i] * cont_africa[i];
## }
## dev = dev + (-2)*normal_lpdf( log_gdp | mu , sigma );
## }R Code 8.13
y <- c(-1, 1)
m8.2 <- map2stan(
alist(y ~ dnorm(mu, sigma),
mu <- alpha),
data = list(y = y),
start = list(alpha = 0, sigma = 1),
chains = 2, cores = 2,
iter = 4000, warmup = 1000
)## Warning: There were 556 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: There were 2 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low## Warning: Examine the pairs() plot to diagnose sampling problems##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 3.9e-05 seconds (Sampling)
## 4.3e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]## Warning in map2stan(alist(y ~ dnorm(mu, sigma), mu <- alpha), data = list(y = y), : There were 556 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.R Code 8.14
precis(m8.2)## Warning in precis(m8.2): There were 556 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 1174967 3438009 -1771670.19 8980565 14 1.1
## sigma 12053381 108539021 350.03 15853582 515 1.0plot(m8.2)
stancode(m8.2)## data{
## int<lower=1> N;
## real y[N];
## }
## parameters{
## real alpha;
## real<lower=0> sigma;
## }
## model{
## vector[N] mu;
## for ( i in 1:N ) {
## mu[i] = alpha;
## }
## y ~ normal( mu , sigma );
## }
## generated quantities{
## vector[N] mu;
## real dev;
## dev = 0;
## for ( i in 1:N ) {
## mu[i] = alpha;
## }
## dev = dev + (-2)*normal_lpdf( y | mu , sigma );
## }
R Code 8.15
m8.3 <- map2stan(
alist(
y ~ dnorm(mu, sigma),
mu <- alpha,
alpha ~ dnorm(1, 10),
sigma ~ dcauchy(0, 1)
),
data = list(y = y),
start = list(alpha = 0, sigma = 1),
chains = 2,
iter = 4000,
warmup = 1000
)##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 1, Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 1, Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 1, Iteration: 1001 / 4000 [ 25%] (Sampling)
## Chain 1, Iteration: 1400 / 4000 [ 35%] (Sampling)
## Chain 1, Iteration: 1800 / 4000 [ 45%] (Sampling)
## Chain 1, Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 1, Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 1, Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 1, Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 1, Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 1, Iteration: 4000 / 4000 [100%] (Sampling)
## Elapsed Time: 0.01494 seconds (Warm-up)
## 0.039659 seconds (Sampling)
## 0.054599 seconds (Total)
##
##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 4000 [ 0%] (Warmup)
## Chain 2, Iteration: 400 / 4000 [ 10%] (Warmup)
## Chain 2, Iteration: 800 / 4000 [ 20%] (Warmup)
## Chain 2, Iteration: 1001 / 4000 [ 25%] (Sampling)
## Chain 2, Iteration: 1400 / 4000 [ 35%] (Sampling)
## Chain 2, Iteration: 1800 / 4000 [ 45%] (Sampling)
## Chain 2, Iteration: 2200 / 4000 [ 55%] (Sampling)
## Chain 2, Iteration: 2600 / 4000 [ 65%] (Sampling)
## Chain 2, Iteration: 3000 / 4000 [ 75%] (Sampling)
## Chain 2, Iteration: 3400 / 4000 [ 85%] (Sampling)
## Chain 2, Iteration: 3800 / 4000 [ 95%] (Sampling)
## Chain 2, Iteration: 4000 / 4000 [100%] (Sampling)
## Elapsed Time: 0.015877 seconds (Warm-up)
## 0.042117 seconds (Sampling)
## 0.057994 seconds (Total)
##
##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 1.9e-05 seconds (Sampling)
## 2.2e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]precis(m8.3)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha 0.08 1.81 -2.29 2.47 1181 1
## sigma 2.07 2.04 0.45 3.61 934 1R Code 8.16
D <-
data_frame(
i = 1:1e4,
y = rcauchy(1e4, 0, 5),
mu = cumsum(y) / i
)
gf_line(mu ~ i, data = D)
R Code 8.17
y <- rnorm(100, mean = 0, sd = 1)R Code 8.18
m8.4 <- map2stan(
alist(y ~ dnorm(mu, sigma),
mu <- a1 + a2,
sigma ~ dcauchy(0, 1)),
data = list(y = y),
start = list(a1 = 0, a2 = 0, sigma = 1),
chains = 2, cores = 2,
iter = 4000,
warmup = 1000
)##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 3e-05 seconds (Sampling)
## 3.3e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]precis(m8.4)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a1 1639.19 802.36 484.04 2913.38 5 1.49
## a2 -1639.08 802.36 -2913.22 -483.92 5 1.49
## sigma 0.98 0.06 0.89 1.07 39 1.03R Code 8.19
m8.5 <- map2stan(
alist(
y ~ dnorm(mu, sigma),
mu <- a1 + a2,
a1 ~ dnorm(0, 10),
a2 ~ dnorm(0, 10),
sigma ~ dcauchy(0, 1)
),
data = list(y = y),
start = list(a1 = 0, a2 = 0, sigma = 1),
chains = 2, cores = 2,
iter = 4000, warmup = 1000
)##
## SAMPLING FOR MODEL 'y ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 3.2e-05 seconds (Sampling)
## 3.5e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 600 / 6000 ]
[ 1200 / 6000 ]
[ 1800 / 6000 ]
[ 2400 / 6000 ]
[ 3000 / 6000 ]
[ 3600 / 6000 ]
[ 4200 / 6000 ]
[ 4800 / 6000 ]
[ 5400 / 6000 ]
[ 6000 / 6000 ]precis(m8.5)## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a1 0.11 7.10 -11.15 11.41 1473 1
## a2 0.00 7.10 -10.88 11.70 1474 1
## sigma 0.97 0.07 0.86 1.08 2214 1R Code 8.20
mp <-
map2stan(
alist(a ~ dnorm(0, 1), b ~ dcauchy(0, 1)),
data = list(y = rep(1,10)),
start = list(a = 0, b = 0),
iter = 1e4, warmup = 100,
WAIC = FALSE
)##
## SAMPLING FOR MODEL 'a ~ dnorm(0, 1)' NOW (CHAIN 1).
## WARNING: The initial buffer, adaptation window, and terminal buffer
## overflow the total number of warmup iterations.
## Defaulting to a 15%/75%/10% partition,
## init_buffer = 15
## adapt_window = 75
## term_buffer = 10
##
##
## Chain 1, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1, Iteration: 101 / 10000 [ 1%] (Sampling)
## Chain 1, Iteration: 1100 / 10000 [ 11%] (Sampling)
## Chain 1, Iteration: 2100 / 10000 [ 21%] (Sampling)
## Chain 1, Iteration: 3100 / 10000 [ 31%] (Sampling)
## Chain 1, Iteration: 4100 / 10000 [ 41%] (Sampling)
## Chain 1, Iteration: 5100 / 10000 [ 51%] (Sampling)
## Chain 1, Iteration: 6100 / 10000 [ 61%] (Sampling)
## Chain 1, Iteration: 7100 / 10000 [ 71%] (Sampling)
## Chain 1, Iteration: 8100 / 10000 [ 81%] (Sampling)
## Chain 1, Iteration: 9100 / 10000 [ 91%] (Sampling)
## Chain 1, Iteration: 10000 / 10000 [100%] (Sampling)
## Elapsed Time: 0.001631 seconds (Warm-up)
## 0.163983 seconds (Sampling)
## 0.165614 seconds (Total)
##
##
## SAMPLING FOR MODEL 'a ~ dnorm(0, 1)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 4e-06 seconds (Warm-up)
## 2.4e-05 seconds (Sampling)
## 2.8e-05 seconds (Total)R Code 8.21
Simulating height and leg length.
N <- 100 # number of individuals
Legs <-
data_frame(
height = rnorm(N, 10, 2), # sim total height of each
leg_prop = runif(N, 0.4, 0.5), # leg as proportion of height
leg_left = leg_prop * height + rnorm(N, 0, 0.02),
leg_right = leg_prop * height + rnorm(N, 0, 0.02)
)R Code 8.22
m5.8s <- map2stan(
alist(
height ~ dnorm(mu, sigma),
mu <- a + bl * leg_left + br * leg_right,
a ~ dnorm(10, 100),
bl ~ dnorm(2, 10),
br ~ dnorm(2, 10),
sigma ~ dcauchy(0, 1)
),
data = Legs,
chains = 4, cores = 4,
start = list(
a = 10, bl = 0,
br = 0, sigma = 1
)
)##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 3e-06 seconds (Warm-up)
## 4e-05 seconds (Sampling)
## 4.3e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]R Code 8
m5.8s2 <- map2stan(
alist(
height ~ dnorm(mu, sigma),
mu <- a + bl * leg_left + br * leg_right,
a ~ dnorm(10, 100),
bl ~ dnorm(2, 10),
br ~ dnorm(2, 10) & T[0, ],
sigma ~ dcauchy(0, 1)
),
data = Legs,
chains = 4, cores = 4,
start = list(
a = 10, bl = 0,
br = 0, sigma = 1
)
)## Warning: There were 2307 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: Examine the pairs() plot to diagnose sampling problems##
## SAMPLING FOR MODEL 'height ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
## performed for num_warmup < 20
##
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## Elapsed Time: 5e-06 seconds (Warm-up)
## 4.1e-05 seconds (Sampling)
## 4.6e-05 seconds (Total)## Computing WAIC## Constructing posterior predictions## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]## Warning in map2stan(alist(height ~ dnorm(mu, sigma), mu <- a + bl * leg_left + : There were 2307 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.