exercises-02.Rmd

# Chapter 2 Exercises

2021-08-19

## Setup

```{r setup,  message=FALSE, warning=FALSE}
knitr::opts_chunk$set(echo = TRUE, comment = "#>", dpi = 300)

rfiles <- list.files(here::here("src"), full.names = TRUE, pattern = "R$")
for (rfile in rfiles) {
  source(rfile)
}

library(glue)
library(tidyverse)
```

> Complete questions 2.1-2.5, 2.8, 2.9, 2.14, 2.17, and 2.22.

## Question 1

**Posterior inference: suppose you have a $\text{Beta}(4,4)$ prior distribution on the probability $\theta$ that a count will yield a 'head' when spun.**
**The coin is spun 10 times and 'heads' appear *fewer* than 3 times.**
**Calculate the exact posterior density for $\theta$ and sketch it.**

prior: $\text{Beta}(4,4)$

data: $y = 0 \text{ or } 1 \text{ or } 2$

- if $y=2$: $p(\theta | y=2) = \text{Beta}(4+2, 4+8)$
- if $y=1$: $p(\theta | y=1) = \text{Beta}(4+1, 4+9)$
- if $y=0$: $p(\theta | y=0) = \text{Beta}(4+0, 4+10)$

$p(\theta|y) = \frac{1}{3} \text{Beta}(6, 12) + \frac{1}{3} \text{Beta}(5, 13) + \frac{1}{3} \text{Beta}4, 14)$

```{r}
theta <- seq(0, 1, 0.01)
prob_density <- (dbeta(theta, 6, 12) + dbeta(theta, 5, 13) + dbeta(theta, 4, 14)) / 3
plot_dist(theta, prob_density, xlab = "theta", ylab = "probability")
```

## Question 2

**Predictive distributions: consider two coins $C_1$ and $C_2$ with the following characteristics: $\Pr(\text{heads} | C_1) = 0.6$ and $\Pr(\text{heads} | C_2) = 0.4$.**
**Choose one of the coins at random and spin it.**
**Given that the first two spins are tails, what is the expectation of the number of additional spins until a heads?**

Find the probability of each coin given the data and use those as "weights" for the expected number of spins to get heads.

$$
p(C_1|y) = \frac{p(C_1) p(y|C_1)}{p(y)} \\
p(C_1) = \frac{1}{2} \\
p(y|C_1) = (1-0.6)^2 = 0.4^2 = \frac{16}{100} \\
p(y) =  \frac{1}{2} \frac{16}{100} + \frac{1}{2} \frac{36}{100} \\
p(C_1|y) = \frac{\frac{1}{2} \frac{16}{100}}{\frac{1}{2} \frac{16}{100} + \frac{1}{2} \frac{36}{100}} = \frac{8}{26}
$$
Same cacuation for $p(C_2|y)$ resulting in $p(C_2|y) = \frac{18}{26}$.

Expected number $n$ of coin spins until get heads given the probability of getting heads $\theta$:

$$
\text{E}(n|\theta) = 1 \theta + 2(1-\theta)\theta + 3(1 - \theta)^2 \theta + \dots = \frac{1}{\theta}
$$

Thus

$$
\begin{aligned}
\text{E}(n|y) &= p(C_1|y) \text{E}(n|C_1,y) + p(C_2|y) \text{E}(n|C_2,y) \\
&= \frac{8}{26} \frac{1}{0.6} + \frac{18}{26} \frac{1}{0.4} \\
&= 2.24
\end{aligned}
$$

## Question 3

**Predictive distributions: let $y$ be the number of 6's in 1000 rolls of a fair die.**

**a) Sketch the approximate distribution of $y$ based on the normal approximation.**

mean: $\text{E}(y) = \frac{1}{6} 1000 = \frac{500}{3}$

std. dev: $\text{sd}(y) = \sqrt{\frac{1}{6} \frac{5}{6} 1000}$

```{r}
mu <- 500 / 3
sigma <- sqrt(1000 * 5 / (6 * 6))
y <- seq(100, 250, 1)
likelihood <- dnorm(y, mu, sigma)
plot_dist(y, likelihood, xlab = "y", ylab = "likelihood")
```

**b) Using the normal distribution table, give approximate 5%, 25%, 50%, 75%, and 95% points for the distribution of $y$.**

| percentile |     z |                                  formula | value |
|------------|-------|------------------------------------------|-------|
|         5% | -1.65 | $\text{E}(y) - 1.65 \times \text{sd}(y)$ | 147.2 |
|        25% | -0.67 | $\text{E}(y) - 1.65 \times \text{sd}(y)$ | 158.8 |
|        50% |     0 |                            $\text{E}(y)$ | 166.7 |
|        75% |  0.67 | $\text{E}(y) + 1.65 \times \text{sd}(y)$ | 174.6 |
|        95% |  1.65 | $\text{E}(y) + 1.65 \times \text{sd}(y)$ | 186.1 |

```{r}
p <- map_chr(c(5, 25, 50, 75, 95), ~ glue("{.x}%"))
q <- c(147.2, 158.8, 166.7, 174.6, 186.1)
pqs <- tibble(p = p, q = q)

d <- tibble(y = y, likelihood = likelihood)
ggplot(d, aes(x = y, y = likelihood)) +
  geom_line() +
  geom_vline(aes(xintercept = q, color = p), data = pqs, show.legend = FALSE) +
  geom_text(aes(x = q, label = p), data = pqs, y = 0.036, show.legend = FALSE) +
  scale_x_continuous(expand = expansion()) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  theme_bw()
```

## Question 4

**Predictive distributions: let $y$ be the number of 6's in 1000 rolls of a die that may not be fair.**
**Let $\theta$ be the probability that the die lands on 6 with the following priors for different values of $\theta$:**

$$
\begin{aligned}
\Pr(\theta = \frac{1}{12} &= 0.25) \\
\Pr(\theta = \frac{1}{6} &= 0.5) \\
\Pr(\theta = \frac{1}{4} &= 0.25) \\
\end{aligned}
$$

**a) Using the normal approximation for the conditional distributions $p(y|\theta)$, sketch the prior predictive for $y$.**

prior predictive: $p(\tilde{y}|M) = \int p(\tilde{y}|\theta, n, M) p(\theta|M)$

for this exercise:
$$
\begin{aligned}
p(\tilde{y}|M) &= p(\tilde{y} | \theta=\frac{1}{12}, n=1000, M) p(\theta=\frac{1}{12}) + \dots \\
&= N(\frac{1}{12} 1000, \sqrt{1000 \frac{1}{12} (1 - \frac{1}{12})}) \times \frac{1}{4} + \dots \\
\end{aligned}
$$
```{r}
six_likelihood <- function(y, theta, n = 1000, prior) {
  mu <- n * theta
  sigma <- sqrt(n * theta * (1 - theta))
  dnorm(y, mu, sigma) * prior
}

y <- seq(30, 310)
separate_likelihoods <- map2(
  c(1 / 12, 1 / 6, 1 / 4),
  c(0.25, 0.5, 0.25),
  ~ six_likelihood(y, .x, n = 1000, prior = .y)
)

combined_likelihood <- accumulate(separate_likelihoods, ~ .x + .y)[[3]]

plot_dist(y, combined_likelihood, "y", "likelihood")
```

---

```{r}
sessionInfo()
```