MCMC Using STAN – Visualization With The Shinystan Package: Solutions

#### Run the folowing lines before doing the exercises 
# ================================================================

## Vector of packages to use in this report
pkgs <- c("shinystan", "rstan", "rstanarm", "bayesplot", "evd")
## 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))

# ===============================================================

###############
#             #
# Exercise 1  #
#             #
###############

launch_shinystan_demo(demo_name = "eight_schools")



###############
#             #
# Exercise 2  #
#             #
###############

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# Reproducible example
set.seed(123)

n            <- 100
mu_sample    <- 10
sigma_sample <- 3

sample_gumb <- rgev(
  n     = n,
  loc   = mu_sample,
  scale = sigma_sample,
  shape = 0
)


###############
#             #
# Exercise 3  #
#             #
###############

stan_model <- '
data{
    int<lower=0> n;
    real y[n];
}

parameters{
  real mu;
  real<lower=0> sigma;
}

model{
  mu    ~ normal(mean(y), 5);
  sigma ~ normal(sd(y), 5);
  y     ~ gumbel(mu,sigma);
}

generated quantities{
  real y_pred;
  y_pred = gumbel_rng(mu, sigma);
}
'


###############
#             #
# Exercise 4  #
#             #
###############

data_list <- list(
  n = n,
  y = sample_gumb
)


###############
#             #
# 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 your available cores
  seed       = 1234
)

# You can always quickly check the traceplots of the two parameters
# traceplot(stan_fit, pars = c("mu", "sigma"))


###############
#             #
# Exercise 6  #
#             #
###############

launch_shinystan(stan_fit, rstudio = TRUE)


###############
#             #
# Exercise 7  #
#             #
###############

# Transform the stan object into a shinstan object and change its name. 
shiny_stan_fit <- as.shinystan(stan_fit, model_name = "Gumbel distribution", ppd = TRUE)

# Run the app 
launch_shinystan(shiny_stan_fit, rstudio = FALSE)



###############
#             #
# Exercise 8  #
#             #
###############

###############
#             #
# Exercise 9  #
#             #
###############

a.

b.

c.



###############
#             #
# Exercise 10 #
#             #
###############

y_pred <- extract(stan_fit, 'y_pred')$y_pred[1:1000] %>% as.numeric()

launch_shinystan(shiny_stan_fit, rstudio = FALSE)