Set up grid approximate for the easy questions

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 )
plot(p_grid, posterior, type = "b")

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

You take the sum of the samples < 0.2 and divide by the number of samples taken from the posterior.

sum(samples < 0.2)/1e4
## [1] 4e-04

4e-04 in decimal form is 0.0004.

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

sum(samples > 0.8)/1e4
## [1] 0.1116

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

sum(samples > 0.2 & samples < 0.8)/1e4
## [1] 0.888

89% of the posterior probability lies between p = 0.2 and p = 0.8.

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

quantile(samples, 0.2)
##       20% 
## 0.5185185

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

quantile( samples, 0.8)
##       80% 
## 0.7557558

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

HPDI(samples, prob = 0.66)
##     |0.66     0.66| 
## 0.5085085 0.7737738

This interval captures the parameters with narrowest posterior probability: 0.24 in width.

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

PI(samples, prob = 0.66)
##       17%       83% 
## 0.5025025 0.7697698

This interval assigns 17% of the probability mass above and below the interval.The width of this interval, 0.32, is wider than with the HDPI method.

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.

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

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

set.seed(100)
samples <- sample( p_grid , prob=posterior , size=1e4 , replace=TRUE )

HPDI(samples, prob = 0.90)
##      |0.9      0.9| 
## 0.3343343 0.7217217

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?

w <- rbinom(1e4 , size=15 , prob=samples )
mean(w == 8)
## [1] 0.1499
sum(w == 8)/1e4
## [1] 0.1499

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

w <- rbinom(1e4 , size=9 , prob=samples )
sum(w == 6)/1e4
## [1] 0.1842

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 <- ifelse( p_grid < 0.5 , 0 , 1)
likelihood <- dbinom(8 , size=15 , prob=p_grid )
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)

plot( posterior ~ p_grid , type="l" )

set.seed(100)
samples <- sample( p_grid , prob=posterior , size=1e4 , replace=TRUE )

HPDI(samples, prob = 0.90)
##      |0.9      0.9| 
## 0.5005005 0.7097097
w <- rbinom(1e4 , size=15 , prob=samples )
sum(w == 8)/1e4
## [1] 0.163
w <- rbinom(1e4 , size=9 , prob=samples )
sum(w == 6)/1e4
## [1] 0.2353

HDPI is narrower now. Probability of observing 8 W in 15 tosses has increased from 0.15 to 0.16, and observing 6 waters has increased from 0.18 to 0.22.