Load packages
library(cmdstanr)
options(mc.cores = 1)
library(posterior)
library(loo)
library(tidyr)
library(dplyr)
options(pillar.neg=FALSE)
library(ggplot2)
library(gridExtra)
library(bayesplot)
theme_set(bayesplot::theme_default(base_family = "sans"))
library(rprojroot)
root<-has_file(".BDA_R_demos_root")$make_fix_file()
SEED <- 48927 # set random seed for reproducability
This notebook contains several examples of how to use Stan in R with cmdstanr. This notebook assumes basic knowledge of Bayesian inference and MCMC. The Stan models are stored in separate .stan-files. The examples are related to Bayesian data analysis course.
Toy data with sequence of failures (0) and successes (1). We would like to learn about the unknown probability of success.
data_bern <- list(N = 10, y = c(1, 1, 1, 0, 1, 1, 1, 0, 1, 0))
Bernoulli model with a proper Beta(1,1) (uniform) prior
code_bern <- root("demos_rstan", "bern.stan")
writeLines(readLines(code_bern))
// Bernoulli model
data {
int<lower=0> N; // number of observations
int<lower=0,upper=1> y[N]; // vector of binary observations
}
parameters {
real<lower=0,upper=1> theta; // probability of success
}
model {
// model block creates the log density to be sampled
theta ~ beta(1, 1); // prior
y ~ bernoulli(theta); // observation model / likelihood
// the notation using ~ is syntactic sugar for
// target += beta_lpdf(theta | 1, 1); // lpdf for continuous theta
// target += bernoulli_lpmf(y | theta); // lpmf for discrete y
// target is the log density to be sampled
//
// y is an array of integers and
// y ~ bernoulli(theta);
// is equivalent to
// for (i in 1:N) {
// y[i] ~ bernoulli(theta);
// }
// which is equivalent to
// for (i in 1:N) {
// target += bernoulli_lpmf(y[i] | theta);
// }
}
Sample form the posterior and show the summary
mod_bern <- cmdstan_model(stan_file = code_bern)
fit_bern <- mod_bern$sample(data = data_bern, seed = SEED)
fit_bern$summary()
Plot the histogram of the posterior draws
draws <- fit_bern$draws()
mcmc_hist(draws, pars='theta')