| title | author | date | output | ||||
|---|---|---|---|---|---|---|---|
Chapter3 |
Julin N Maloof |
4/18/2019 |
|
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 )hist(samples)
How much posterior probability lies below p = 0.2?
sum(samples < .2) / length(samples)## [1] 5e-04
How much posterior probability lies above p = 0.8?
sum(samples > 0.8) / length(samples)## [1] 0.1117
How much posterior probability lies between p = 0.2 and p = 0.8?
sum(samples >0.2 & samples < 0.8) / length(samples)## [1] 0.8878
20% of the posterior probability lies below which value of p?
quantile(samples, 0.2)## 20%
## 0.5195195
20% of the posterior probability lies above which value of p?
quantile(samples, 0.8)## 80%
## 0.7567568
Which values of p contain the narrowest interval equal to 66% of the posterior probability?
HPDI(samples, 0.66)## |0.66 0.66|
## 0.5205205 0.7847848
Which values of p contain 66% of the posterior probability, assuming equal posterior probability both below and above the interval?
PI(samples, 0.66)## 17% 83%
## 0.5005005 0.7687688
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(0,1, length.out = 1000)
prior <- rep(1,length(p_grid))
likelihood <- dbinom(8, 15, prob = p_grid)
plot(likelihood)
unstd.post <- prior*likelihood
post <- unstd.post / sum(unstd.post)
plot(post)
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, size = 10000, replace = TRUE, prob = post)
hist(samples)
HPDI(samples, 0.9)## |0.9 0.9|
## 0.3383383 0.7317317
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?
predictions <- rbinom(10000, 15, prob = samples)
sum(predictions==8) / length(predictions)## [1] 0.1428
Using the posterior distribution constructed from the new (8/15) data, now calculate the probability of observing 6 water in 9 tosses.
predictions9 <- rbinom(10000, 9, prob = samples)
sum(predictions9==6) / length(predictions9)## [1] 0.1695
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(0,1, length.out = 1000)
prior <- rep(1,length(p_grid)) %>% ifelse(p_grid < .5, 0, .)
likelihood <- dbinom(8, 15, prob = p_grid)
plot(likelihood)
unstd.post <- prior*likelihood
post <- unstd.post / sum(unstd.post)
plot(post)
Draw 10,000 samples from the grid approximation from above. Then use the samples to calculate the 90% HPDI for p.
samples2 <- sample(p_grid, size = 10000, replace = TRUE, prob = post)
hist(samples2)
HPDI(samples2, 0.9)## |0.9 0.9|
## 0.5005005 0.7097097
Confidence intervals are tighter
3H1. Using grid approximation, compute the posterior distribution for the probability of a birth being a boy. Assume a uniform prior probability. Which parameter value maximizes the posterior probability?
library(rethinking)
data(homeworkch3)
sum(birth1) + sum(birth2)## [1] 111
p_grid <- seq(0,1, length.out = 1000)
prior <- rep(1,1000)
likelihood <- dbinom(111, size=200, prob=p_grid)
posterior <- prior*likelihood
posterior <- posterior / sum(posterior)
plot(p_grid, posterior)
p_grid[which.max(posterior)]## [1] 0.5545546
Using the sample function, draw 10,000 random parameter values from the posterior distribution you calculated above. Use these samples to estimate the 50%, 89%, and 97% highest posterior density intervals.
sample <- sample(p_grid, size=10000, replace = TRUE, prob=posterior)
HPDI(sample,c(.5,.89,.97))## |0.97 |0.89 |0.5 0.5| 0.89| 0.97|
## 0.4784785 0.4964965 0.5255255 0.5735736 0.6076076 0.6276276
Use rbinom to simulate 10,000 replicates of 200 births. You should end up with 10,000 num- bers, each one a count of boys out of 200 births. Compare the distribution of predicted numbers of boys to the actual count in the data (111 boys out of 200 births). There are many good ways to visualize the simulations, but the dens command (part of the rethinking package) is probably the easiest way in this case. Does it look like the model fits the data well? That is, does the distribution of predictions include the actual observation as a central, likely outcome?
sims <- rbinom(10000, size = 200, prob = sample)
hist(sims)
dens(sims)
mean(sims)## [1] 110.7069
median(sims)## [1] 111
sum(sims==111) / length(sims)## [1] 0.0398
sum(sims > 109 & sims < 113) / length(sims)## [1] 0.1168
##3H4. Now compare 10,000 counts of boys from 100 simulated first borns only to the number of boys in the first births, birth1. How does the model look in this light?
sims <- rbinom(10000, size = 100, prob = sample)
hist(sims)
dens(sims)
cat("mean: ", mean(sims))## mean: 55.3854
cat("median: ", median(sims))## median: 55
cat("observed: ",sum(birth1))## observed: 51
sum(sims==51) / length(sims)## [1] 0.0528
sum(sims==55) / length(sims)## [1] 0.0639
HPDI(sims,.95)## |0.95 0.95|
## 43 66
clearly not great
3H5. The model assumes that sex of first and second births are independent. To check this assump- tion, focus now on second births that followed female first borns. Compare 10,000 simulated counts of boys to only those second births that followed girls. To do this correctly, you need to count the number of first borns who were girls and simulate that many births, 10,000 times. Compare the counts of boys in your simulations to the actual observed count of boys following girls. How does the model look in this light? Any guesses what is going on in these data?
#number of female first borns:
sum(birth1==0)## [1] 49
#number of boys after female first borns
sum(birth1==0 & birth2==1)## [1] 39
sims <- rbinom(10000, size = 49, prob = sample)
hist(sims)
dens(sims)
cat("mean: ", mean(sims))## mean: 27.2215
cat("median: ", median(sims))## median: 27
cat("observed: ",sum(birth1))## observed: 51
sum(sims==39) / length(sims)## [1] 9e-04
HPDI(sims,.95)## |0.95 0.95|
## 19 34
sum(sims>=39) / length(sims)## [1] 0.0014
model predicts we should see 27 boys out of 49, but we observed 39. This is outside of the 95% confidnce interval.