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. Problems are labeled Easy (E), Medium (M), and Hard(H).

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

Simulate the posterior distribution

p_grid <- seq( from=0 , to=1 , length.out=1000 )
prior <- rep( 1 , 1000 )
likelihood <- dbinom( 6 , size=9 , prob=p_grid )
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)
set.seed(100)
samples <- sample( p_grid , prob=posterior , size=1e4 , replace=TRUE )

Questions

Use the values in samples to answer the questions that follow.

3E1. How much posterior probability lies below p = 0.2?

#using samples to calculate posterior probability below p = 0.2
mean(samples < 0.2)
## [1] 4e-04

3E2. How much posterior probability lies above p = 0.8?

#using samples to calculate posterior probability above p = 0.8
mean(samples > 0.8)
## [1] 0.1116

3E3. How much posterior probability lies between p = 0.2 and p = 0.8?

#using samples to calculate posterior probability below p = 0.2 and above p = 0.8
mean(samples > 0.2 & samples < 0.8)
## [1] 0.888

3E4. 20% of the posterior probability lies below which value of p?

#20% of the posterior probability lies below following p value
quantile(samples, 0.2)
##       20% 
## 0.5185185

3E5. 20% of the posterior probability lies above which value of p?

#20% of the posterior probability lies above following p value
quantile(samples, 0.8)
##       80% 
## 0.7557558

3E6. Which values of p contain the narrowest interval equal to 66% of the posterior probability?

#loading required library
library(rethinking)
## Loading required package: rstan
## Loading required package: StanHeaders
## Loading required package: ggplot2
## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
## Loading required package: parallel
## rethinking (Version 2.13)
## 
## Attaching package: 'rethinking'
## The following object is masked from 'package:stats':
## 
##     rstudent
#Following values of p contain the narrowest interval equal to 66% of the posterior probability- Using Highest Posterior Density Intervals (HPDI) function
HPDI(samples, prob = 0.66)
##     |0.66     0.66| 
## 0.5085085 0.7737738

3E7. Which values of p contain 66% of the posterior probability, assuming equal posterior probability both below and above the interval?

#Using Percentile Intervals function
PI(samples, prob = 0.66)
##       17%       83% 
## 0.5025025 0.7697698

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

likelihood <- dbinom( 8 , size=15 , prob=p_grid )
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)

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

samples <- sample( p_grid , prob=posterior , size=1e4 , replace=TRUE )
#90% HPDI for p
HPDI(samples, prob = 0.9)
##      |0.9      0.9| 
## 0.3293293 0.7167167

3M3. 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?

nw <- rbinom(1e4, size = 15, prob = samples)
#probability of observing 8 water in 15 tosses
mean(nw == 8)
## [1] 0.1444

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

nw <- rbinom(1e4, size = 9, prob = samples)
#probability of observing 8 water in 15 tosses
mean(nw == 6)
## [1] 0.1751

3M5. Start over at 3M1, 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. What difference does the better prior make? If it helps, compare inferences (using both priors) to the true value p = 0.7.

p_grid <- seq( from=0 , to=1 , length.out=1000 )
prior_2 <- ifelse(p_grid < 0.5, 0, 1 )
likelihood <- dbinom( 8 , size= 15 , prob=p_grid )
posterior <- likelihood * prior_2
posterior <- posterior / sum(posterior)
#Plot
plot(p_grid, posterior, type = "b")

samples_2 <- sample( p_grid , prob=posterior , size=1e4 , replace=TRUE )

mean(samples_2 < 0.5)
## [1] 0
mean(samples_2 > 0.5)
## [1] 1
mean(samples_2 > 0.5 & samples_2 < 0.8)
## [1] 0.9889
PI(samples_2, prob = 0.66)
##       17%       83% 
## 0.5315315 0.6796797
HPDI(samples_2, prob = 0.66)
##     |0.66     0.66| 
## 0.5005005 0.6306306
nw <- rbinom(1e4, size = 15, prob = samples_2)
nw <- rbinom(1e4, size = 9, prob = samples_2)

3M6. 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?

#PI Width function 
pi_width <- function(N, true_p) #N- number of observations and true_p- true value of p
{
prob_1 <- dbinom(round(N*true_p), size = N, prob = p_grid)
posterior <- prob_1 * prior
posterior <- posterior / sum(posterior)
samples <- sample( p_grid , prob=posterior , size=1e4 , replace=TRUE )  
#99% percentile interval of the posterior distribution
int <- PI(samples, prob = 0.99)
names(int) <- NULL
diff(int)
}

#Using PI Width function with different number of observations for the true p value of 0.7
#When N = 500 and true_p probability = 0.7
true_p <- 0.7
N <- 500
pi_width(N, true_p)
## [1] 0.1061061
#When N = 1000 and true_p probability = 0.7
true_p <- 0.7
N <- 1000
pi_width(N, true_p)
## [1] 0.07507508
#When N = 1500 and true_p probability = 0.7
true_p <- 0.7
N <- 1500
pi_width(N, true_p)
## [1] 0.06006006
#When N = 2000 and true_p probability = 0.7
true_p <- 0.7
N <- 2000
pi_width(N, true_p)
## [1] 0.05405405
#Seems like N = 2000 and p value 0.7 is giving the distance between upper and lower bound interval to be 0.05