Setup

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

1 Introduction

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.

2 Bernoulli model

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')