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.

# Import library
library(rethinking)

# Select number of data points for grid approximation
N_POINTS <- 1e3

# Define grid
p_grid <- seq(from=0, to=1, length.out=N_POINTS)

# Define prior
prob_p <- rep(1, N_POINTS)

# Compute likelihood at each value in grid (or prob_data)
NUM_WATER <- 8
NUM_TOSSES <- 15
likelihood <- dbinom(NUM_WATER, size=NUM_TOSSES, prob=p_grid)

# Compute product of likelihood and prior
posterior <- likelihood * prob_p

# Standardize the posterior, so it sums to one
posterior <- posterior / sum(posterior)

# Plot the posterior
plot(x=p_grid,
     y=posterior,
     type='l',
     main='Posterior of Globe Tossing Experiment',
     xlab='Probability of Water',
     ylab='Posterior Probability')

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

# Number of samples
N_SAMPLES <- 1e4
set.seed(42)

# Create samples
samples <- sample(p_grid, prob=posterior, size=N_SAMPLES, replace=TRUE)

# Calculate 90% Highest Posterior Density Interval (HPDI) for p
hpdi_val <- HPDI(samples, prob=0.9)
cat('The 90% Highest Posterior Density Interval (HPDI) is between,', hpdi_val[1], 'and', hpdi_val[2])
## The 90% Highest Posterior Density Interval (HPDI) is between, 0.3393393 and 0.7267267

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 outcome
water <- rbinom(N_SAMPLES, size=NUM_TOSSES, prob=samples)
simplehist(water)

# Compute probability
prob_water <- sum(water==NUM_WATER) / N_SAMPLES
cat('The probability of observing 8 water in 15 tosses is', prob_water)
## The probability of observing 8 water in 15 tosses is 0.1464

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 outcome given new parameters
NEW_NUM_TOSSES <- 9
NEW_NUM_WATER <- 6
new_water <- rbinom(N_SAMPLES, size=NEW_NUM_TOSSES, prob=samples)
simplehist(new_water)

# Compute probability
prob_new_water <- sum(new_water==NEW_NUM_WATER) / N_SAMPLES
cat('The probability of observing 6 water in 9 tosses is', prob_new_water)
## The probability of observing 6 water in 9 tosses is 0.1711

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.

# Select number of data points for grid approximation
N_POINTS <- 1e3

# Define grid
p_grid <- seq(from=0, to=1, length.out=N_POINTS)

# Define prior
prior <- ifelse(p_grid < 0.5, 0, 1)

# Compute likelihood at each value in grid (or prob_data)
NUM_WATER <- 8
NUM_TOSSES <- 15
likelihood <- dbinom(NUM_WATER, size=NUM_TOSSES, prob=p_grid)

# Compute product of likelihood and prior
posterior <- likelihood * prior

# Standardize the posterior, so it sums to one
posterior <- posterior / sum(posterior)

# Plot the posterior
plot(x=p_grid,
     y=posterior,
     type='l',
     main='Posterior of Globe Tossing Experiment',
     xlab='Probability of Water',
     ylab='Posterior Probability')

################################################################################
# Number of samples
N_SAMPLES <- 1e4

# Create samples
samples <- sample(p_grid, prob=posterior, size=N_SAMPLES, replace=TRUE)

# Calculate 90% Highest Posterior Density Interval (HPDI) for p
hpdi_val <- HPDI(samples, prob=0.9)
cat('The 90% Highest Posterior Density Interval (HPDI) is between,', hpdi_val[1], 'and', hpdi_val[2], '\n')
## The 90% Highest Posterior Density Interval (HPDI) is between, 0.5005005 and 0.7097097
################################################################################
# Simulate the outcome
water <- rbinom(N_SAMPLES, size=NUM_TOSSES, prob=samples)
simplehist(water)

# Compute probability
prob_water <- sum(water==NUM_WATER) / N_SAMPLES
cat('The probability of observing 8 water in 15 tosses is', prob_water, '\n')
## The probability of observing 8 water in 15 tosses is 0.1631
################################################################################
# Simulate the outcome given new parameters
NEW_NUM_TOSSES <- 9
NEW_NUM_WATER <- 6
new_water <- rbinom(N_SAMPLES, size=NEW_NUM_TOSSES, prob=samples)
simplehist(new_water)

# Compute probability
prob_new_water <- sum(new_water==NEW_NUM_WATER) / N_SAMPLES
cat('The probability of observing 6 water in 9 tosses is', prob_new_water, '\n')
## The probability of observing 6 water in 9 tosses is 0.2254

  1. Plot The new prior resulted in the posterior for probability of water < 0.5 to be 0, while the distribution for probability of water > 0.5 has increased in comparison to the flat posterior.

  2. 90% HPDI A better prior resulted in a narrower range for the 90% HPDI. This reduction in range means we have a higher confident that the real value falls between this new interval. Note that the lower limit has increase from 0.3 to 0.5, while the upper limit has stayed similar.

  3. 8 waters in 15 tosses With this new prior, while the distribution of the outcomes looks similar, note that there are now higher frequency at 10 instead of 8, which is opposite of the outcome from the previous prior.

  4. 6 waters in 9 tosses With this new prior, while the distribution of the outcomes looks similar, note that the distribution now centers around 6 instead of 5, which was what it was when simulating using the previous prior.

Overall, this new prior shows us what happens to our simulations when we eliminate the probability of water < 0.5. Looking at the simulation results, we see that the HPDI range closer to 0.7 than it was before which is encouraging as the true value of p is 0.7.

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?

# Create function to compute number of tosses needed
compute_perc_interval_width <- function(NUM_TOSSES, TRUE_PROB){
  # Select number of data points for grid approximation
  N_POINTS <- 1e3
  
  # Define grid
  p_grid <- seq(from=0, to=1, length.out=N_POINTS)
  
  # Compute likelihood
  likelihood <- dbinom(round(NUM_TOSSES * TRUE_PROB), size=NUM_TOSSES, prob=p_grid)
  
  # Compute product of likelihood and prior
  posterior <- likelihood * prior
  
  # Standardize the posterior, so it sums to one
  posterior <- posterior / sum(posterior)
  
  # Number of samples
  N_SAMPLES <- 1e6
  PROB_LIMIT <- 0.99
  
  # Create samples
  samples <- sample(p_grid, prob=posterior, size=N_SAMPLES, replace=TRUE)
  
  # Get Percent Interval
  perc_interval <- PI(samples, prob=PROB_LIMIT)
  
  # Compute interval width
  names(perc_interval) <- NULL
  diff <- diff(perc_interval)
}

# Find number of tosses that resulted with percent interval width closest to 0.05
# True probability
TRUE_PROB <- 0.7

# Number of tosses to test
num_tosses <- c(2100:2200)

# Create array to store width of percent interval
width_list <- c(1:length(num_tosses))

# Compute width for tosses of interests
for (toss in num_tosses) {
  width_list[toss-(num_tosses[1]-1)] <- compute_perc_interval_width(toss, TRUE_PROB)
}

# Find location of smallest width to 0.05
diff_of_diff <- width_list - 0.05
min_loc <- which.min(diff_of_diff)
num_of_tosses <- num_tosses[1] + min_loc

cat('The number of tosses needed to get to a 99% interval with a PI width of 0.05 and is', num_of_tosses)
## The number of tosses needed to get to a 99% interval with a PI width of 0.05 and is 2142