Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

# SR2 Chapter 3 Hard

Here’s my solutions to the hard exercises in chapter 3 of McElreath’s Statistical Rethinking, 2nd edition.

$$\DeclareMathOperator{\dbinomial}{Binomial} \DeclareMathOperator{\dbernoulli}{Bernoulli} \DeclareMathOperator{\dpoisson}{Poisson} \DeclareMathOperator{\dnormal}{Normal} \DeclareMathOperator{\dt}{t} \DeclareMathOperator{\dcauchy}{Cauchy} \DeclareMathOperator{\dexponential}{Exp} \DeclareMathOperator{\duniform}{Uniform} \DeclareMathOperator{\dgamma}{Gamma} \DeclareMathOperator{\dinvpamma}{Invpamma} \DeclareMathOperator{\invlogit}{InvLogit} \DeclareMathOperator{\logit}{Logit} \DeclareMathOperator{\ddirichlet}{Dirichlet} \DeclareMathOperator{\dbeta}{Beta}$$

Let’s first put the data into a tibble for easier manipulation later.

data(homeworkch3)

df <- tibble(birth1 = birth1, birth2 = birth2) %>%
mutate(birth = row_number())
The first few rows of the data.
birth1 birth2 birth
1 0 1
0 1 2
0 0 3
0 1 4
1 0 5
1 1 6

## 3H1

Let’s check we have the correct total cound and the correct number of boys.

h1_counts <- df %>%
gather(order, gender, -birth) %>%
summarise(boys = sum(gender), births = n())

Now we can grid approximate the posterior as before.

granularity <- 1000

h1_grid <- tibble(p = seq(0, 1, length.out = granularity)) %>%
mutate(prior = 1)

h1_posterior <- h1_grid %>%
mutate(
likelihood = dbinom(h1_counts$boys, h1_counts$births, p),
posterior = prior * likelihood,
posterior = posterior / sum(posterior)
)

The maximum a posteriori (MAP) value is the value of p that maximises the posterior.

h1_map <- h1_posterior %>%
slice(which.max(posterior)) %>%
pull(p)

h1_map
 0.5545546

## 3H2

We draw samples with weight equalt to the posterior. We then apply the HPDI function to these samples, each time with a different width.

h2_samples <- h1_posterior %>%
sample_n(10000, replace = TRUE, weight = posterior) %>%
pull(p)

h2_hpdi <- h2_samples %>%
crossing(prob = c(0.5, 0.89, 0.97)) %>%
group_by(prob) %>%
group_map(HPDI)

h2_hpdi
[]
|0.5      0.5|
0.4574575 0.5735736

[]
|0.89     0.89|
0.4534535 0.6606607

[]
|0.97     0.97|
0.4294294 0.6616617 

## 3H3

The posterior predictive samples are possible observations according to our posterior.

h3_posterior_predictive <- rbinom(10000, 200, h2_samples)

The number of observed births is very close to the MAP of the posterior predictive distribution, suggesting we have a decent fit.

## 3H4

Our data are from birth pairs and so far we didn’t make any distinction between the first and second births. To test this assumption, we can perform a posterior predictive check as in 3H3, but this time for first births.

h4_posterior_predictive <- rbinom(10000, 100, h2_samples)

The fit doesn’t look quite as good for first births as it did for all births together. It also doesn’t look bad since there is still a fair bit of probability mass around the observed number of first birth boys.

## 3H5

As the final posterior predictive check, let’s check the number of boys born after a girl.

h5_counts <- df %>%
filter(birth1 == 0) %>%
summarise(boys = sum(birth2), births = n())

h5_posterior_predictive <- rbinom(10000, h5_counts\$births, h2_samples)