Below are the solutions to these exercises on “MCMC Using STAN – An Introduction with The RSTAN Package: Exercises.”
#### Run the folowing lines before doing the exercises # ================================================================ ## Vector of packages to use in this report pkgs <- c("rstan", "faraway", "tidyverse", "parallel") ## Install packages From CRAN if they are not installed. for (p in pkgs) { if (p %in% rownames(installed.packages()) == FALSE) { install.packages(p, repos = 'http://cran.us.r-project.org') } } ## Load the packages for (p in pkgs) suppressPackageStartupMessages(library(p, quietly = T, character.only = T)) ## Personalized theme for ggplot (just for aesthetics of the plots) "theme_piss" <- function(size_p = 18, size_c = 14, size_l = 12, theme = theme_bw(), ...) { text <- function(size, ...) element_text(size = size, colour = "#33666C", face = "bold", ...) theme + theme( plot.title = text(size_p, hjust = 0.5), axis.title = text(size_c), legend.title = text(size_l) ) + theme( legend.background = element_rect(colour = "black"), legend.key = element_rect(fill = "white"), ... ) } # ===============================================================
############### # # # Exercise 1 # # # ############### rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores()) ############### # # # Exercise 2 # # # ############### data(pulp, package = "faraway") str(pulp)
## 'data.frame': 20 obs. of 2 variables: ## $ bright : num 59.8 60 60.8 60.8 59.8 59.8 60.2 60.4 59.9 60 ... ## $ operator: Factor w/ 4 levels "a","b","c","d": 1 1 1 1 1 2 2 2 2 2 ...
# or glimpse(pulp)
## Observations: 20 ## Variables: 2 ## $ bright <dbl> 59.8, 60.0, 60.8, 60.8, 59.8, 59.8, 60.2, 60.4, 59.9,... ## $ operator <fctr> a, a, a, a, a, b, b, b, b, b, c, c, c, c, c, d, d, d...
# Plot ggplot(pulp, aes(x = operator, y = bright)) + geom_jitter(width = 0.1) + theme_piss()
# Yes, a model with random operator would be preferable. ############### # # # Exercise 3 # # # ############### stan_model <- " data { int<lower=0> N; int<lower=0> J; vector[N] response; int<lower=1,upper=J> predictor[N]; } parameters { vector[J] eta; real mu; real<lower=0> sigma_alpha; real<lower=0> sigma_epsilon; } transformed parameters { vector[J] a; vector[N] yhat; a = mu + sigma_alpha * eta; for (i in 1:N) yhat[i] = a[predictor[i]]; } model { eta ~ normal(0, 1); response ~ normal(yhat, sigma_epsilon); } " ############### # # # Exercise 4 # # # ############### data_list <- list( N = nrow(pulp), J = nlevels(pulp$operator), response = pulp$bright, predictor = as.numeric(pulp$operator) ) ############### # # # Exercise 5 # # # ############### # Reproducible example for the whole script set.seed(1234) # Compiling and producing posterior samples from the model. stan_fit <- stan( model_code = stan_model, data = data_list, chains = 4, iter = 1000, warmup = 200, cores = parallel::detectCores(), # Use all you available cores seed = 1234 )
############### # # # Exercise 6 # # # ############### stan_fit2 <- stan( model_code = stan_model, data = data_list, chains = 4, iter = 1000, warmup = 200, cores = parallel::detectCores(), # Use all you available cores seed = 1234, control = list(adapt_delta = 0.95, max_treedepth = 20) ) # Only 1 divergent transition is left after warmp-up ! # matrix of output plots from the stan object pairs(stan_fit2, pars = c("mu", "sigma_alpha", "sigma_epsilon", "a"))
############### # # # Exercise 7 # # # ############### traceplot(stan_fit2, pars = c("mu", "sigma_alpha", "sigma_epsilon", "a"))
############### # # # Exercise 8 # # # ############### sampler_par <- get_sampler_params(stan_fit2, inc_warmup = TRUE) # a. Chains by chains lapply(sampler_par, summary, digits = 2)
## [[1]] ## accept_stat__ stepsize__ treedepth__ n_leapfrog__ ## Min. :0.00 Min. :3.5e-04 Min. : 0.0 Min. : 1 ## 1st Qu.:0.96 1st Qu.:1.2e-01 1st Qu.: 4.0 1st Qu.: 15 ## Median :0.99 Median :1.2e-01 Median : 5.0 Median : 31 ## Mean :0.93 Mean :1.2e-01 Mean : 4.5 Mean : 129 ## 3rd Qu.:1.00 3rd Qu.:1.2e-01 3rd Qu.: 5.0 3rd Qu.: 31 ## Max. :1.00 Max. :1.1e+01 Max. :15.0 Max. :37119 ## divergent__ energy__ ## Min. :0.000 Min. :-1.3e+01 ## 1st Qu.:0.000 1st Qu.:-6.6e+00 ## Median :0.000 Median :-4.3e+00 ## Mean :0.007 Mean : 1.5e+03 ## 3rd Qu.:0.000 3rd Qu.:-1.4e+00 ## Max. :1.000 Max. : 1.1e+06 ## ## [[2]] ## accept_stat__ stepsize__ treedepth__ n_leapfrog__ ## Min. :0.00 Min. : 0.0015 Min. : 0.0 Min. : 1 ## 1st Qu.:0.95 1st Qu.: 0.1747 1st Qu.: 4.0 1st Qu.: 15 ## Median :0.98 Median : 0.1747 Median : 4.0 Median : 31 ## Mean :0.94 Mean : 0.1726 Mean : 4.4 Mean : 49 ## 3rd Qu.:1.00 3rd Qu.: 0.1747 3rd Qu.: 5.0 3rd Qu.: 31 ## Max. :1.00 Max. :10.9517 Max. :13.0 Max. :8959 ## divergent__ energy__ ## Min. :0.000 Min. :-12.0 ## 1st Qu.:0.000 1st Qu.: -7.1 ## Median :0.000 Median : -5.0 ## Mean :0.007 Mean : -2.5 ## 3rd Qu.:0.000 3rd Qu.: -2.6 ## Max. :1.000 Max. :552.5 ## ## [[3]] ## accept_stat__ stepsize__ treedepth__ n_leapfrog__ ## Min. :0.00 Min. : 0.0004 Min. : 0.0 Min. : 1 ## 1st Qu.:0.94 1st Qu.: 0.0176 1st Qu.: 6.0 1st Qu.: 127 ## Median :0.98 Median : 0.0176 Median : 7.0 Median : 191 ## Mean :0.95 Mean : 0.0280 Mean : 6.9 Mean : 371 ## 3rd Qu.:1.00 3rd Qu.: 0.0176 3rd Qu.: 8.0 3rd Qu.: 255 ## Max. :1.00 Max. :10.9517 Max. :15.0 Max. :49979 ## divergent__ energy__ ## Min. :0.000 Min. : -15.1 ## 1st Qu.:0.000 1st Qu.: -7.2 ## Median :0.000 Median : -5.2 ## Mean :0.013 Mean : 115.3 ## 3rd Qu.:0.000 3rd Qu.: -2.5 ## Max. :1.000 Max. :37196.9 ## ## [[4]] ## accept_stat__ stepsize__ treedepth__ n_leapfrog__ ## Min. :0.00 Min. :5.3e-04 Min. : 0.0 Min. : 1 ## 1st Qu.:0.94 1st Qu.:1.6e-01 1st Qu.: 4.0 1st Qu.: 15 ## Median :0.98 Median :1.6e-01 Median : 4.0 Median : 31 ## Mean :0.94 Mean :1.5e-01 Mean : 4.6 Mean : 166 ## 3rd Qu.:1.00 3rd Qu.:1.6e-01 3rd Qu.: 5.0 3rd Qu.: 31 ## Max. :1.00 Max. :1.1e+01 Max. :15.0 Max. :39135 ## divergent__ energy__ ## Min. :0.000 Min. :-1.4e+01 ## 1st Qu.:0.000 1st Qu.:-6.8e+00 ## Median :0.000 Median :-4.7e+00 ## Mean :0.019 Mean : 1.6e+03 ## 3rd Qu.:0.000 3rd Qu.:-2.0e+00 ## Max. :1.000 Max. : 1.3e+06
# b. The whole aggregated chain summary(do.call(rbind, sampler_par), digits = 2)
## accept_stat__ stepsize__ treedepth__ n_leapfrog__ ## Min. :0.00 Min. :3.5e-04 Min. : 0.0 Min. : 1 ## 1st Qu.:0.95 1st Qu.:1.8e-02 1st Qu.: 4.0 1st Qu.: 15 ## Median :0.98 Median :1.2e-01 Median : 5.0 Median : 31 ## Mean :0.94 Mean :1.2e-01 Mean : 5.1 Mean : 179 ## 3rd Qu.:1.00 3rd Qu.:1.6e-01 3rd Qu.: 6.0 3rd Qu.: 79 ## Max. :1.00 Max. :1.1e+01 Max. :15.0 Max. :49979 ## divergent__ energy__ ## Min. :0.000 Min. :-1.5e+01 ## 1st Qu.:0.000 1st Qu.:-6.9e+00 ## Median :0.000 Median :-4.8e+00 ## Mean :0.011 Mean : 8.1e+02 ## 3rd Qu.:0.000 3rd Qu.:-2.2e+00 ## Max. :1.000 Max. : 1.3e+06
############### # # # Exercise 9 # # # ############### print(stan_fit2, pars = c("mu", "sigma_alpha", "sigma_epsilon", "a"))
## Inference for Stan model: 316aef3d3154fd7bedd9618e29c60ef5. ## 4 chains, each with iter=1000; warmup=200; thin=1; ## post-warmup draws per chain=800, total post-warmup draws=3200. ## ## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat ## mu 60.41 0.01 0.33 59.76 60.27 60.41 60.55 61.14 697 1.01 ## sigma_alpha 0.51 0.02 0.43 0.05 0.25 0.38 0.63 1.70 486 1.01 ## sigma_epsilon 0.36 0.00 0.07 0.25 0.31 0.35 0.40 0.53 1092 1.00 ## a[1] 60.28 0.00 0.15 59.98 60.18 60.28 60.38 60.58 3200 1.00 ## a[2] 60.14 0.00 0.16 59.83 60.03 60.13 60.25 60.48 1611 1.00 ## a[3] 60.57 0.00 0.16 60.25 60.46 60.57 60.67 60.87 2105 1.00 ## a[4] 60.62 0.00 0.16 60.29 60.51 60.62 60.72 60.92 1655 1.00 ## ## Samples were drawn using NUTS(diag_e) at Thu Jan 25 16:29:50 2018. ## 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).
# Rhat is very close to 1 for any parameters and hence we have support that convergence did happen.
Related exercise sets:
Hacking statistics or: How I Learned to Stop Worrying About Calculus and Love Stats Exercises (Part-2)
Scripting Loops exercises Vol.2
Using factor variables like a pro (Part-2)
Parallel Computing Exercises: Foreach and DoParallel (Part-2)
Hacking statistics or: How I Learned to Stop Worrying About Calculus and Love Stats Exercises (Part-6)
Leave a Reply