Exercise 3.2 of ”Bayesian Workflow”

exercise_3.2.R

library(cmdstanr)
library(posterior)
library(ggplot2)

model_file <- file.path("models", "exercise_3.2.stan")

stan_data <- list(N = 171, W = 111, q = 0.5)

mod <- cmdstan_model(model_file, pedantic = TRUE)
fit <- mod$sample(data = stan_data, seed = 1)
fit$summary()

draws <- fit$draws() |>
  as_draws_df()

freq <- cut(draws$p, breaks = seq(0, 1, 0.01)) |>
    table()
p <- seq(0, 0.99, 0.01) + 0.005

df <- data.frame(p = p, dens = c(freq))
ggplot(df, aes(x = p, y = dens)) +
  geom_col() +
  scale_y_continuous(labe = \(x) x / nrow(draws) / 0.01) +
  labs(title = "Prior: Beta(1, 1)",
       x = "Proportion of honest participants",
       y = "Posterior probability") +
  theme_gray(base_size = 8, base_family = "Helvetica")
ggsave("ex32-1.png", width = 800, height = 600, units = "px")

exercise_3.2.stan

data {
  int<lower=0> N;           // Num. participants
  int<lower=0, upper=N> W;  // Num. participants claim the prize
  real<lower=0, upper=1> q; // Prob. win the prize
}

parameters {
  real<lower=0, upper=1> p;  // Prob. honest
}

model {
  vector[W + 1] lp;
  for (n in (N - W):N) { // Num. honest participants
    lp[W - (N - n) + 1] = binomial_lpmf(n | N, p)
                         + binomial_lpmf(W - (N - n) | n, q);
  }
  target += log_sum_exp(lp);

  // Flat prior
  target += beta_lpdf(p | 1, 1);
}