BDA3 Chapter 2 Exercise 11

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

Assume the sampling distribution is \(\dcauchy(y \mid \theta, 1)\) with uniform prior \(p(\theta) \propto 1\) on \([0, 100]\). Given observations \(y\), we can approximate the posterior for \(\theta\) by dividing the interval \([0, 100]\) into partitions of length \(\frac{1}{m}\). The unnormalised posterior for \(\theta\) on this grid is then computed as follows.

# observations
y <- c(43, 44, 45, 46.5, 47.5) 

# grid granularity
m <- 100

# L(θ) := p(y | θ)
likelihood <- function(theta) 
  y %>% 
    map(dcauchy, theta, 1) %>% 
    reduce(prod)

# unnormalised posterior grid
posterior_unnorm <- tibble(theta = seq(0, 100, 1 / m)) %>% 
  mutate(density = map(theta, likelihood) %>% unlist())

We can approximate the normalising constant by summing the approximate area on each partition. Each partition has width \(\frac{1}{m}\) and approximate height given by the density, so the approximate area is the multiple of the two.

# grid approx to area under curve
normalising_constant <- posterior_unnorm %>% 
  summarise(sum(density) / m) %>% 
  pull()

# normalised posterior grid
posterior <- posterior_unnorm %>% 
  mutate(density = density / normalising_constant)

normalising_constant

[1] 3.418359e-05

Let’s zoom in on the region \([40, 50]\) where most of the density lies.

Sampling from this posterior yields a histogram with a similar shape.

posterior_draws <- posterior %>% 
  sample_n(1000, replace = TRUE, weight = density) %>% 
  select(theta)

We can draw from the posterior predictive distribution by first drawing \(\tilde\theta\) from the posterior of \(\theta\), then drawing \(\tilde y\) from \(\dcauchy(\tilde\theta, 1)\). The tails of the posterior predictive distribution are much wider than for \(\theta\) so we plot this histogram on the interval \([10, 90]\) (although there are a few observations outside this interval).

posterior_predictive <- posterior_draws %>% 
  mutate(pp = rcauchy(n(), theta, 1))