Bayesian Statistics Exercises in R: 20 Practice Problems

prior_a <- 2
prior_b <- 2
successes <- 14
failures  <- 100 - successes
ex_1_1 <- c(prior_a + successes, prior_b + failures)
ex_1_1
#> [1] 16 88

Explanation: When the prior on a proportion is Beta(a, b) and the data are binomial with k successes out of n trials, conjugacy gives a closed-form Beta(a+k, b+n-k) posterior. Beta(2, 2) is a soft prior centered at 0.5 with an equivalent weight of about 2 successes plus 2 failures. The resulting Beta(16, 88) posterior shifts strongly toward 0.16 because the data dominate the prior at n=100. This identity is the foundation for streaming A/B updates without refitting.

post_a <- 16
post_b <- 88
post_mean <- post_a / (post_a + post_b)
ci <- qbeta(c(0.025, 0.975), post_a, post_b)
ex_1_2 <- c(mean = post_mean, lower = ci[1], upper = ci[2])
ex_1_2
#>      mean     lower     upper
#> 0.1538462 0.0935023 0.2256617

Explanation: The mean of Beta(a, b) is a / (a + b), and equal-tailed quantiles come straight from qbeta(). The 95% credible interval is a direct probability statement about the parameter: there is a 95% posterior probability the true conversion rate lies between 9.4% and 22.6%, unlike a frequentist confidence interval which is a statement about the procedure. A highest-density interval would be slightly narrower because the Beta posterior is right-skewed.

prior_shape <- 2
prior_rate  <- 1
total_calls <- 120
hours       <- 8
post_shape <- prior_shape + total_calls
post_rate  <- prior_rate + hours
post_mean  <- post_shape / post_rate
ex_1_3 <- c(shape = post_shape, rate = post_rate, mean = post_mean)
ex_1_3
#>      shape       rate       mean
#> 122.000000   9.000000  13.555556

Explanation: For Poisson observations with a Gamma(alpha, beta) prior on the rate, the posterior is Gamma(alpha + total counts, beta + total exposure time), so the Gamma is conjugate to the Poisson and updating is just adding sufficient statistics. The posterior mean of 122/9 (about 13.56 calls per hour) is a precision-weighted average of the prior mean (2) and the empirical rate (15), tilted toward the data because the prior carries only 1 hour of equivalent weight versus 8 hours of observation.

p_grid <- seq(0, 1, length.out = 1000)
prior  <- rep(1, 1000)
lik    <- dbinom(7, size = 10, prob = p_grid)
post   <- (lik * prior) / sum(lik * prior)
ex_1_4 <- sum(p_grid * post)
ex_1_4
#> [1] 0.6666007

Explanation: Grid approximation discretizes the parameter space, evaluates prior times likelihood at each grid point, and normalizes the weights to sum to one. The posterior mean is then a weighted average of the grid values. The exact Beta(8, 4) posterior under a uniform prior has mean 8/12 = 0.6667, so the 1,000-point grid agrees to four decimals. Grid approximation breaks down past three or four parameters because the grid size grows exponentially with dimensionality, which is why MCMC takes over.

ex_2_1 <- rstanarm::stan_glm(
  mpg ~ wt,
  data    = mtcars,
  family  = gaussian(),
  iter    = 2000,
  chains  = 4,
  refresh = 0,
  seed    = 1
)
print(ex_2_1)
#> stan_glm
#>  family:       gaussian [identity]
#>  formula:      mpg ~ wt
#>             Median MAD_SD
#> (Intercept) 37.3   1.9
#> wt          -5.3   0.6
#> sigma        3.2   0.4

Explanation: stan_glm() is the rstanarm wrapper around Stan that fits the same families as base glm() but returns full posterior draws rather than maximum likelihood estimates. Default rstanarm priors are weakly informative normals centered at zero with scales auto-adjusted to the data, so for moderate sample sizes the posterior is essentially data-driven. The printed output shows posterior medians and median absolute deviations rather than point estimates plus standard errors. Set refresh = 0 to suppress chain progress messages in batch scripts.

ex_2_2 <- rstanarm::stan_glm(
  mpg ~ wt,
  data            = mtcars,
  family          = gaussian(),
  prior           = rstanarm::normal(location = -5, scale = 2),
  prior_intercept = rstanarm::normal(location = 20, scale = 2.5),
  iter            = 2000,
  chains          = 4,
  refresh         = 0,
  seed            = 1
)
rstanarm::prior_summary(ex_2_2)

Explanation: The prior argument controls the prior on slope coefficients while prior_intercept controls the prior on the centered intercept; rstanarm internally centers predictors before sampling for efficiency, so the intercept prior is on that centered intercept, not the raw one. rstanarm::normal() is the prior constructor and accepts vector arguments when you want a different normal per coefficient. Always confirm with prior_summary() because rstanarm sometimes rescales user priors when autoscale = TRUE.

ex_2_3 <- rstanarm::stan_glm(
  am ~ mpg + hp,
  data    = mtcars,
  family  = binomial(link = "logit"),
  iter    = 2000,
  chains  = 4,
  refresh = 0,
  seed    = 1
)
print(ex_2_3)

Explanation: For a binomial family, stan_glm() samples directly from the posterior of the log-odds coefficients, so the reported medians are interpretable as posterior median log-odds-ratios per unit change in the predictor. The MAD_SD column reports median absolute deviation, a robust scale measure that is more reliable than the posterior standard deviation when the posterior is heavy-tailed. With only 32 observations the intercept posterior is very wide; the Bayesian framing carries that uncertainty forward into predictions instead of pretending it does not exist.

ex_2_4 <- rstanarm::posterior_interval(ex_2_1, prob = 0.90)
ex_2_4
#>                   5%        95%
#> (Intercept)  34.1438   40.5113
#> wt           -6.3084   -4.3110
#> sigma         2.6051    3.9523

Explanation: posterior_interval() extracts equal-tailed credible intervals from the stored MCMC draws. The default prob = 0.9 is the rstanarm convention because the authors argue 95% intervals overstate certainty in noisy regression scenarios. For a highest-density interval instead of equal-tailed bounds, use bayestestR::hdi(ex_2_1, ci = 0.9), which gives narrower bounds when the posterior is skewed (commonly for variance components or rare-event probabilities).

ex_3_1 <- brms::brm(
  mpg ~ wt,
  data    = mtcars,
  prior   = brms::prior(normal(-4, 1), class = "b", coef = "wt"),
  chains  = 4,
  iter    = 2000,
  refresh = 0,
  seed    = 1
)
summary(ex_3_1)

Explanation: brms compiles a Stan program from the R formula and runs Hamiltonian Monte Carlo. Priors are declared with brms::prior() keyed by class (b for coefficients, Intercept for the centered intercept, sigma for residual SD). When the prior conflicts with the data, the posterior compromises between them weighted by their precision; here Bayes pulls the OLS slope of about -5.3 toward the prior mean of -4 because n = 32 carries only moderate evidence. Run a prior sensitivity sweep before publishing.

ex_3_2 <- brms::brm(
  mpg ~ wt + (1 | cyl),
  data    = mtcars,
  chains  = 4,
  iter    = 2000,
  refresh = 0,
  seed    = 1
)
summary(ex_3_2)

Explanation: The (1 | cyl) term tells brms to estimate one varying intercept per cylinder group, drawn from a shared normal prior whose standard deviation sd(Intercept) is itself estimated from the data. This is partial pooling: a group with few cars is shrunk toward the global mean, while a group with many cars stays close to its own data. Use ranef(ex_3_2) to inspect the per-group intercept deviations; with only three groups the partial-pooling shrinkage is strong because the group-level variance is poorly identified.

fit_tight <- brms::brm(
  mpg ~ wt, data = mtcars,
  prior   = brms::prior(normal(0, 0.1), class = "b"),
  chains  = 4, iter = 2000, refresh = 0, seed = 1
)
fit_wide <- brms::brm(
  mpg ~ wt, data = mtcars,
  prior   = brms::prior(normal(0, 100), class = "b"),
  chains  = 4, iter = 2000, refresh = 0, seed = 1
)
ex_3_3 <- c(
  tight = median(brms::as_draws_df(fit_tight)$b_wt),
  wide  = median(brms::as_draws_df(fit_wide)$b_wt)
)
ex_3_3
#>       tight        wide
#> -0.19030  -5.34218

Explanation: With only 32 observations the data information is finite, so a tight prior centered at zero with scale 0.1 dominates and shrinks the slope toward zero, while the diffuse normal(0, 100) prior recovers essentially the OLS slope. This is the textbook bias-variance trade-off of prior choice: tight priors reduce posterior variance but introduce bias toward the prior mean. Always run a prior sensitivity analysis like this before publishing Bayesian results in any regulated reporting setting.

ex_3_4 <- brms::prior_summary(ex_3_2)
ex_3_4

Explanation: brms picks weakly informative defaults from a Student-t family with 3 degrees of freedom and location and scale derived from the response data; this is more robust than a normal default because the heavier tails accommodate occasional large coefficients. Coefficients (class = "b") default to flat (improper) priors. Always run prior_summary() before publishing because regulated reporting environments like FDA submissions or EMA reviews require that priors be explicitly disclosed rather than left at unstated defaults.

draws <- posterior::as_draws(ex_2_1)
ex_4_1 <- max(posterior::summarise_draws(draws, "rhat")$rhat)
ex_4_1
#> [1] 1.001

Explanation: R-hat (the potential scale reduction factor) compares the within-chain variance to the between-chain variance; values below 1.01 indicate the chains have mixed well and are sampling from the same target distribution. If max(rhat) > 1.05 the chains have not converged and you should run longer, reparameterize, or fix divergences before trusting the posterior. The posterior package is the modern unified interface and supersedes rstan::summary() and coda::gelman.diag() for diagnostics on any MCMC output.

draws <- posterior::as_draws_df(ex_2_1)
ex_4_2 <- posterior::ess_bulk(draws$wt)
ex_4_2
#> [1] 3984.521

Explanation: Bulk ESS is the effective sample size for estimating the bulk of the posterior (mean, median) and corrects for autocorrelation between successive MCMC draws. The rule of thumb is ess_bulk >= 400 for stable median estimates; values below that suggest you should run more iterations or reparameterize. Tail ESS (posterior::ess_tail) is a separate diagnostic for the 5% and 95% quantiles and can be lower than bulk ESS when the posterior has heavy tails; always report both for any published quantity.

ex_4_3 <- bayesplot::mcmc_trace(ex_2_1, pars = "wt")
ex_4_3

Explanation: A healthy trace plot looks like overlapping fuzzy caterpillars with no visible drift, structural breaks, or stuck patches. If one chain wanders off to a different region of parameter space, you have multimodality or a label-switching problem; if all chains crawl slowly with high autocorrelation, increase iter or reparameterize (e.g., switch to non-centered parameterization for hierarchical models). bayesplot is part of the Stan ecosystem and accepts both rstanarm and brms fits directly without conversion.

y_rep <- rstanarm::posterior_predict(ex_2_1, draws = 50)
ex_4_4 <- bayesplot::ppc_dens_overlay(y = mtcars$mpg, yrep = y_rep)
ex_4_4

Explanation: A posterior predictive check simulates new outcomes from the fitted model and compares them to the actual data. If the model captures the data-generating process, the replicates should overlap the observed density. Systematic gaps (e.g., the observed has two modes but the replicates have one) reveal model mis-specification such as an unmodeled subgroup. ppc_dens_overlay() is the workhorse PPC for continuous responses; ppc_bars() is the discrete-count equivalent and ppc_stat() checks calibration of any summary statistic.

ex_5_1 <- rstanarm::posterior_interval(ex_2_1, prob = 0.95, pars = "wt")
ex_5_1
#>       2.5%     97.5%
#> wt -6.5089  -4.1207

Explanation: posterior_interval() returns equal-tailed quantile-based intervals from the stored MCMC draws. The default prob is 0.9 in rstanarm; pass prob = 0.95 to match conventional frequentist reporting when the audience expects 95% intervals. The interpretation is fully Bayesian: there is a 95% posterior probability the true slope lies in [-6.51, -4.12], conditional on the model and prior. This is a direct statement about the parameter, not the procedure.

ex_5_2 <- bayestestR::hdi(ex_2_1, ci = 0.95)
ex_5_2

Explanation: An HDI is the narrowest interval that contains 95% of the posterior density, so every value inside the HDI has higher posterior density than any value outside. For symmetric unimodal posteriors the HDI and equal-tailed interval coincide, but for skewed posteriors (variance parameters, rare-event probabilities) the HDI is more informative and sometimes drops the bound at zero entirely. bayestestR also exposes eti() for equal-tailed intervals and a flexible ci() wrapper that switches methods.

preds <- rstanarm::posterior_predict(
  ex_2_1,
  newdata = data.frame(wt = 3.0),
  draws   = 4000
)
ex_5_3 <- c(
  mean = mean(preds),
  q05  = unname(quantile(preds, 0.05)),
  q95  = unname(quantile(preds, 0.95))
)
ex_5_3
#>     mean      q05      q95
#> 21.45    16.27    26.61

Explanation: A posterior predictive distribution propagates both parameter uncertainty (from the posterior over the coefficients) and residual uncertainty (from sigma) to give the full predictive distribution of a new observation. The 90% predictive interval is much wider than the 90% credible interval for the mean response because it incorporates residual scatter. Always quote predictive intervals (not coefficient intervals) when forecasting a single new outcome rather than an average.

fit_wt_cyl <- rstanarm::stan_glm(
  mpg ~ wt + cyl, data = mtcars,
  family  = gaussian(),
  iter    = 2000, chains = 4, refresh = 0, seed = 1
)
ex_5_4 <- loo::loo_compare(
  loo::loo(ex_2_1),
  loo::loo(fit_wt_cyl)
)
ex_5_4
#>            elpd_diff se_diff
#> fit_wt_cyl   0.0       0.0
#> ex_2_1      -3.6       2.5

Explanation: loo::loo() computes Pareto-smoothed importance-sampling leave-one-out cross-validation (PSIS-LOO), which estimates out-of-sample expected log predictive density (elpd_loo) without refitting the model n times. loo_compare() reports differences in elpd_loo with the best model on top. A difference smaller than 2-3 times se_diff is not strong evidence either way; here -3.6 with se_diff = 2.5 is only weak evidence that adding cyl improves out-of-sample fit.