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 )

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

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

To get what proportion of the values in samples are less than 0.2, I use logical indexing, sume the logicals, then divide by the total number of values.

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

0.1% of the posterior probability lies below p=0.2.

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

By the same logic, I change less than 0.2 to greater than 0.8.

mean( samples > 0.8 )
## [1] 0.1116

Around 12% of the posterior probability lies above p=0.8.

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

mean( samples > 0.2 & samples < 0.8 )
## [1] 0.888

Around 88% 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, probs = 0.2 )
##       20% 
## 0.5185185

20% of the posterior probability is below – 0.5185185

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

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

20% of the posterior probability is above– 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

66% of the posterior probability is between 0.51 and 0.77

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

To get the interval that has equal probability mass in each tail, we need to find the percentile interval (page 56)

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

The values of 0.50<p<0.77 contain 66% of the posterior probability with equal posterior probability both below and above this interval.

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 = 1e3)
prior <- rep(1, 1e3)
likelihood <- dbinom(8, size = 15, prob = p_grid)
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)
plot(x = p_grid, y = posterior, type = "l")

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)
HPDI(samples, prob = 0.9 )
##      |0.9      0.9| 
## 0.3293293 0.7167167

The 90% HPDI is thus between p=0.33 and p=0.71.

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)
simplehist(w)

The probability of observing 8 water in 15 tosses is equal to

mean( w == 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.

w <- rbinom(1e4, size = 9, prob = samples)
simplehist(w)

mean(w == 6)
## [1] 0.1751

The probability of observing of 6 water in 9 tosses is 17.6%.

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 = 1e3)
prior <- ifelse(p_grid < .5, 0, 1)
likelihood <- dbinom(8, size = 15, prob = p_grid)
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)
plot(x = p_grid, y = posterior, type = "l")

With the new, informed prior, the probability mass below 0.5 is zero. The posterior with the flat prior still puts quite a bit of probability mass below 0.5:

mean( samples < 0.5 )
## [1] 0.3957

We can compare how much probability mass the two different models put on the interval between 0.6 and 0.8 (this includes the true value 0.7):

mean( samples > 0.6 & samples < 0.8 )
## [1] 0.2787

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?

I want to try some values by using true value of p for the proportion of water to generate some “observed” data.

compute_pi_width <- function(N, true_p) {
  likelihood <- dbinom( round(N*true_p), size=N, prob=p_grid )
  posterior <- likelihood * prior
  posterior <- posterior / sum(posterior)
  samples <- sample(p_grid, prob=posterior, size=1e4, replace=TRUE )
  interval <- PI(samples, prob=0.99)
  names(interval) <- NULL
  diff( interval )
}

true_p <- 0.7
N <- 10
compute_pi_width(N, true_p)
## [1] 0.4274274

With 10 observations, we’re not even close yet.

N <- 100
compute_pi_width(N, true_p)
## [1] 0.2302352

Still some way to go.

N <- 1000
compute_pi_width(N, true_p)
## [1] 0.07507508

With 1000 observations, we’re getting quite close to an PI of 0.05 width. With N = 2200, we then get a 95%-interval that is 0.05 wide.

N <- 2200
compute_pi_width(N, true_p)
## [1] 0.05105105
N <- 150
compute_pi_width(N, true_p = 0.995)
## [1] 0.04704705

Seemingly, using a more informed prior would also influence how many observations are needed.