bern-with-target.R

library(cmdstanr)

m1 <- cmdstan_model("bern-with-target.stan") # compile the model
data <- list(N = 5, y = c(0, 1, 1, 0, 1), a = 2, b = 7)

f1 <- m1$sample(       # for other options to sample, help(sample)
  data = data,         # pass data as a list, match the vars name to Stan
  seed = 123,          # to reproduce results, Stan does not rely on R's seed
  chains = 4,          # total chains, the more, the better
  parallel_chains = 4, # for multi-processor CPUs
  iter_warmup = 500,   # number of draws for warmup (per chain)
  iter_sampling = 1000  # number of draws for samples (per chain)
)
f1$summary()

m2 <- cmdstan_model("bern.stan") # compile the model
f2 <- m2$sample(       # for other options to sample, help(sample)
  data = data,         # pass data as a list, match the vars name to Stan
  seed = 123,          # to reproduce results, Stan does not rely on R's seed
  chains = 4,          # total chains, the more, the better
  parallel_chains = 4, # for multi-processor CPUs
  iter_warmup = 500,   # number of draws for warmup (per chain)
  iter_sampling = 1000  # number of draws for samples (per chain)
)
f2$summary()

library(tidybayes)
library(ggplot2)

theta1 <- spread_draws(f1, theta)
theta2 <- spread_draws(f2, theta)
d <- data.frame(theta1 = theta1$theta, theta2 = theta2$theta)

p <- ggplot(aes(x = theta1), data = d)
p + geom_density() + geom_density(aes(x = theta2), color = 'red')