RStan in CoCalc (Kernel: "R-Project")

# switch to PNG images, because the default of SVG produces very large files
options(jupyter.plot_mimetypes = c('image/png'))

pkgbuild::has_build_tools(debug = TRUE)

packageDescription("rstan")$Version

require(rstan)

# make sure to only use one core!
options(mc.cores = 1)

rstan_options(auto_write = TRUE)

schools8 <- "
data {
  int<lower=0> J;         // number of schools 
  real y[J];              // estimated treatment effects
  real<lower=0> sigma[J]; // standard error of effect estimates 
}
parameters {
  real mu;                // population treatment effect
  real<lower=0> tau;      // standard deviation in treatment effects
  vector[J] eta;          // unscaled deviation from mu by school
}
transformed parameters {
  vector[J] theta = mu + tau * eta;        // school treatment effects
}
model {
  target += normal_lpdf(eta | 0, 1);       // prior log-density
  target += normal_lpdf(y | theta, sigma); // log-likelihood
}"

schools_dat <- list(J = 8,
                    y = c(28, 8, -3, 7, -1, 1, 18, 12),
                    sigma = c(15, 10, 16, 11, 9, 11, 10, 18)
                   )

fit <- stan(model_code = schools8, data = schools_dat)

print(fit)

plot(fit)

pairs(fit, pars = c("mu", "tau", "lp__"))

### return an array of three dimensions: iterations, chains, parameters 
a <- extract(fit, permuted = FALSE) 
a[1:10]

### use S3 functions on stanfit objects
a2 <- as.array(fit)
m <- as.matrix(fit)
d <- as.data.frame(fit)
d[1:10]