BDA3 Chapter 3 Exercise 3

Here’s my solution to exercise 3, chapter 3, of Gelman’s Bayesian Data Analysis (BDA), 3rd edition. There are solutions to some of the exercises on the book’s webpage.

Suppose we have \(n\) measurements \(y \mid \mu, \sigma \sim \dnorm(\mu, \sigma)\), where the prior \(p(\mu, \log \sigma) \propto 1\) is uniform. The calculations on page 66 show that the marginal posterior distribution of \(\mu\) is \(\mu \mid y \sim \dt(\bar y, s / n)\), where \(s\) is the sample standard deviation. The measurements are as follows.

control <- list(
  n = 32,
  mean = 1.013,
  sd = 0.24
)

treatment <- list(
  n = 36,
  mean = 1.173,
  sd = 0.2
)

The t-distribution in base-R is a standardised t-distribution. For a more general t-distribution (with arbitrary location and scale), we’ll use the LaplacesDemon package.

library(LaplacesDemon)

This allows us to plot the marginal posterior means.

mp <- tibble(value = seq(0, 2, 0.01)) %>% 
  mutate(
    ctrl = dst(value, control$mean, control$sd / sqrt(control$n), control$n - 1),
    trt = dst(value, treatment$mean, treatment$sd / sqrt(treatment$n), treatment$n - 1)
  ) %>% 
  gather(cohort, density, ctrl, trt) 

The 95% credible interval of the marginal posterior means is:

draws <- 10000

difference <- tibble(draw = 1:draws) %>% 
  mutate(
    ctrl = rst(n(), control$mean, control$sd / sqrt(control$n), control$n - 1),
    trt = rst(n(), treatment$mean, treatment$sd / sqrt(treatment$n), treatment$n - 1),
    difference = trt - ctrl
  ) 

ci <- difference$difference %>% 
  quantile(probs = c(0.05, 0.95)) 

ci

        5%        95% 
0.06752309 0.25173035 

We can also plot the distribution of the difference.