rm(list = ls())

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.

#install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
#options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
#install.packages('rethinking',type='source')
library(rethinking)
p_grid_tossing <- seq(from = 0, to = 1, length.out = 1000)

#Every prior is posterior of some other inference

prior_tossing <- rep(1, 1000)

likelihood_tossing <- dbinom(8, size = 15, prob = p_grid_tossing)

# Every posterior is a prior for next observation. Bayesian estimate is always posterior distribution
# over parameters, Pr(parameters|data)

posterior_tossing <- likelihood_tossing * prior_tossing
posterior_tossing <- posterior_tossing / sum(posterior_tossing)

plot(x = p_grid_tossing, y = posterior_tossing, type = "b", main = "Posterior distribution for Globe Tossing Data", xlab = "Grid Tossing", ylab = "Posetrior Tossing")

#plot(x = p_grid_tossing, y = posterior_tossing, type = "l", main = "Posterior distribution for Globe Tossing Data", xlab = "Grid Tossing", ylab = "Posetrior Tossing")

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

samples_10k = sample(p_grid_tossing, prob = posterior_tossing, size = 10000, replace = T)

# HPDI - These functions compute highest posterior density (HPDI) and percentile (PI) intervals, using samples from a posterior density or simulated outcomes.

HPDI(samples_10k, prob = .90)
##      |0.9      0.9| 
## 0.3413413 0.7317317
plot(x = p_grid_tossing, y = posterior_tossing, type = "h", main = "Posterior distribution for Globe Tossing Data", xlab = "Grid Tossing", ylab = "Posetrior Tossing")

#hist(HPDI, xlim = "Frequency", main = "Highest Posterior Density", breaks = #10, col = "black")

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?

obs_probability = rbinom(10000, size = 15, prob = samples_10k)

mean(obs_probability == 8)
## [1] 0.1507
hist(obs_probability, xlab = "Probability", ylab = "Frequency", main = "Observation Probability", breaks = 10, col = "lightgray")

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

obs_probability = rbinom(10000, size = 9, prob = samples_10k)

mean(obs_probability == 6)
## [1] 0.1759
hist(obs_probability, xlab = "Probability", ylab = "Frequency", main = "Observation Probability", breaks = 10, col = "lightgray")

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.

p_grid_tossing <- seq(from = 0, to = 1, length.out = 1000)

#Every prior is posterior of some other inference

prior_tossing <- ifelse(p_grid_tossing < 0.5, 0, 1)

likelihood_tossing <- dbinom(8, size = 15, prob = p_grid_tossing)

# Every posterior is a prior for next observation. Bayesian estimate is always posterior distribution
# over parameters, Pr(parameters|data)

posterior_tossing <- likelihood_tossing * prior_tossing
posterior_tossing <- posterior_tossing / sum(posterior_tossing)

plot(x = p_grid_tossing, y = posterior_tossing, type = "b", main = "Posterior distribution for Globe Tossing Data", xlab = "Grid Tossing", ylab = "Posetrior Tossing")

The type of the curve that was produced is very similar to the one that was produced with original posterior distribution of global data. The structure has changed significantly with prior condition. The posterior probability has flattened out i.e. has become 0 whenever p < 0.5 condition is satisfied. The rest of the curve shifts to the original bell shape when p < 0.5 condition is not met.

samples_10k = sample(p_grid_tossing, prob = posterior_tossing, size = 10000, replace = T)

# HPDI - These functions compute highest posterior density (HPDI) and percentile (PI) intervals, using samples from a posterior density or simulated outcomes.

HPDI(samples_10k, prob = .90)
##      |0.9      0.9| 
## 0.5005005 0.7097097
plot(x = p_grid_tossing, y = posterior_tossing, type = "h", main = "Posterior distribution for Globe Tossing Data", xlab = "Grid Tossing", ylab = "Posetrior Tossing")

In the original HPDI calculation, after changing sample value to 10000, we get HDPI values as

 |0.9      0.9|

0.3393393 0.7217217

but when we change p value condition for prior, we get HDPI value as

 |0.9      0.9|

0.5005005 0.7117117

This signifies that change in prior values have significant effect on the lower range of HDPI as compared to upper range.

obs_probability = rbinom(10000, size = 15, prob = samples_10k)

mean(obs_probability == 8)
## [1] 0.1576
hist(obs_probability, xlab = "Probability", ylab = "Frequency", main = "Observation Probability", breaks = 10, col = "lightgray")

When we construct a posterior predictive check for the new model, we can see that the histogram distribution is skewed as compared to the originsl model. The probability of observing 8 water in 15 tosses shifted to 15.6% from 15.8% as compared to the original probability.

obs_probability = rbinom(10000, size = 9, prob = samples_10k)

mean(obs_probability == 6)
## [1] 0.23
hist(obs_probability, xlab = "Probability", ylab = "Frequency", main = "Observation Probability", breaks = 10, col = "lightgray")

When we calculated the probability of observing 6 water in 9 tosses, we get probability around 18.21% but when we apply new prior condition and calculate the probability, we see significant change in histogram shape and probability value. The probability value shifted to 22.82% and histrogram is no more symmetric losing it's bell curved shape. To summrize, the new p-value condition has infact improved the probability.

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?

p_grid_tossing <- seq(from = 0, to = 1, length.out = 1000)

#Every prior is posterior of some other inference

prior_tossing <- rep(1, 1000)

prob_tossing = 0.05

obs_probability = rbinom(1000, size = 9, prob = prob_tossing)

likelihood_tossing <- dbinom(obs_probability, size = 15, prob = p_grid_tossing)

# Every posterior is a prior for next observation. Bayesian estimate is always posterior distribution
# over parameters, Pr(parameters|data)

posterior_tossing <- likelihood_tossing * prior_tossing
posterior_tossing <- posterior_tossing / sum(posterior_tossing)

plot(x = p_grid_tossing, y = posterior_tossing, type = "b", main = "Posterior distribution for Globe Tossing Data", xlab = "Grid Tossing", ylab = "Posetrior Tossing")