Title: ‘Assignment #2’ Author: “Cem Soylu” Date: “2022-08-16” Output: html_document

Chapter 3 - Sampling the Imaginary

This chapter introduced the basic procedures for manipulating posterior distributions. Our fundamental tool is samples of parameter values drawn from the posterior distribution. These samples can be used to produce intervals, point estimates, posterior predictive checks, as well as other kinds of simulations. Posterior predictive checks combine uncertainty about parameters, as described by the posterior distribution, with uncertainty about outcomes, as described by the assumed likelihood function. These checks are useful for verifying that your software worked correctly. They are also useful for prospecting for ways in which your models are inadequate.

Place each answer inside the code chunk (grey box). The code chunks should contain a text response or a code that completes/answers the question or activity requested. Make sure to include plots if the question requests them.

Finally, upon completion, name your final output .html file as: YourName_ANLY505-Year-Semester.html and publish the assignment to your R Pubs account and submit the link to Canvas. Each question is worth 5 points

Questions

3-1. Suppose the globe tossing data had turned out to be 8 water in 15 tosses. Construct the posterior distribution, using grid approximation. Plot the posterior. Use the same flat prior as before.

library(dplyr)
library(ggplot2)

# using grid approximation
p_grid <- seq(from = 0, to = 1, length.out = 1000)
prior <- rep(1, 1000)

# tossing data had turned out to be 8 water in 15 tosses
likelihood <- dbinom(8, size = 15, prob = p_grid)

# calculating the posterior
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)

# Plot the posterior
tibble(p = p_grid, posterior = posterior) %>%
  ggplot(aes(x = p, y = posterior)) +
  geom_point() +
  geom_line() +
  labs(x = "Proportion Water/Tosses (p)", y = "Posterior Density")

3-2. Draw 10,000 samples from the grid approximation from above. Then use the samples to calculate the 90% HPDI for p.

library(rethinking)

set.seed(2022)

#  Draw 10,000 samples from the grid approximation from above

samples <- sample(p_grid, prob = posterior, size = 10000, replace = TRUE)

# Use the samples to calculate the 90% HPDI for p

HPDI(samples, prob = 0.9)

##      |0.9      0.9| 
## 0.3413413 0.7357357

# 90% HPDI for p is (0.341, 0.736)

3-3. Construct a posterior predictive check for this model and data. This means simulate the distribution of samples, averaging over the posterior uncertainty in p. What is the probability of observing 8 water in 15 tosses?

# simulate the distribution of samples
set.seed(2022)

water <- rbinom(10000, size = 15, prob = samples)

# calculating the probability of observing 8 water in 15 tosses
mean(water == 8)

## [1] 0.1474

# the probability of observing 8 water in 15 tosses is 14.74%

3-4. Using the posterior distribution constructed from the new (8/15) data, now calculate the probability of observing 6 water in 9 tosses.

# simulate the distribution of samples
set.seed(2022)

water <- rbinom(10000, size = 9, prob = samples)

# calculating the probability of observing 6 water in 9 tosses

mean(water == 6)

## [1] 0.1789

# the probability of observing 6 water in 9 tosses is 17.89%

3-5. Start over at 3-1, but now use a prior that is zero below p = 0.5 and a constant above p = 0.5. This corresponds to prior information that a majority of the Earth’s surface is water. Repeat each problem above and compare the inferences. Plot the posterior. What difference does the better prior make? If it helps, compare inferences (using both priors) to the true value p = 0.7.

# using grid approximation
p_grid <- seq(from = 0, to = 1, length.out = 1000)
prior <- case_when(p_grid < 0.5 ~ 0L, TRUE ~ 1L)

# simulate the distribution of samples
set.seed(2022)

likelihood <- dbinom(8, size = 15, prob = p_grid)

# calculating the posterior
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)

# Repeat 3.1 - Plot the posterior

tibble(p = p_grid, posterior = posterior) %>%
  ggplot(aes(x = p, y = posterior)) +
  geom_point() +
  geom_line() +
  labs(x = "Proportion Water/Tosses (p)", y = "Posterior Density")

# Repeat 3.2 - calculating the 90% HPDI for p

samples <- sample(p_grid, prob = posterior, size = 10000, replace = TRUE)
HPDI(samples, prob = 0.9)

##      |0.9      0.9| 
## 0.5005005 0.7157157

# Answer
#     |0.9      0.9| 
#0.5005005 0.7157157 


# Repeat 3.3 - calculating the probability of observing 8 water in 15 tosses
water <- rbinom(10000, size = 15, prob = samples)
mean(water == 8)

## [1] 0.1566

# Answer
# probability = 0.1566


# Repeat 3.4 - calculating the probability of observing 6 water in 9 tosses
water <- rbinom(10000, size = 9, prob = samples)
mean(water == 6)

## [1] 0.2358

# Answer
# probability = 0.2358

3-6. Suppose you want to estimate the Earth’s proportion of water very precisely. Specifically, you want the 99% percentile interval of the posterior distribution of p to be only 0.05 wide. This means the distance between the upper and lower bound of the interval should be 0.05. How many times will you have to toss the globe to do this?

library(purrr)
library(tidyr)

single_sim <- function(tosses, prior_type = c("uniform", "step")) {
  prior_type <- match.arg(prior_type)
  obs <- rbinom(1, size = tosses, prob = 0.7)
  
  p_grid <- seq(from = 0, to = 1, length.out = 1000)
  prior <- rep(1, 1000)
  if (prior_type == "step") prior[1:500] <- 0
  
  likelihood <- dbinom(obs, size = tosses, prob = p_grid)
  posterior <- likelihood * prior
  posterior <- posterior / sum(posterior)
  
  samples <- sample(p_grid, prob = posterior, size = 10000, replace = TRUE)
  
  # 99% percentile interval of the posterior distribution of p to be only 0.05 wide
  
  interval <- PI(samples, prob = 0.99)
  width <- interval[2] - interval[1]
}


single_cond <- function(tosses, prior_type, reps = 100) {
  tibble(tosses = tosses,
         prior_type = prior_type,
         width = map_dbl(seq_len(reps), ~single_sim(tosses = tosses,
                                                    prior_type = prior_type)))
}

simulation <- crossing(tosses = seq(1000, 5000, by = 100),
                       prior_type = c("uniform", "step")) %>%
  pmap_dfr(single_cond, reps = 100) %>%
  group_by(tosses, prior_type) %>%
  summarize(avg_width = mean(width), .groups = "drop") %>%
  mutate(prior_type = case_when(prior_type == "uniform" ~ "Uniform Prior",
                                prior_type == "step" ~ "Step Prior"),
         prior_type = factor(prior_type, levels = c("Uniform Prior",
                                                    "Step Prior")))

ggplot(simulation, aes(x = tosses, y = avg_width)) +
  facet_wrap(~prior_type, nrow = 1) +
  geom_point() +
  geom_line() +
  labs(x = "Tosses", y = "Average Interval Width") +
  theme(panel.spacing.x = unit(2, "lines"))

# This figure below shows that the the average interval width is approximately 0.074 given 100 simulations, the true proportion p is 0.7, and if we toss the globe 1,000 times. 
# To get an interval width of .05 or smaller, we would need around 2,300 tosses.