Easy
p_grid = seq(0,1,length=1000)
prior=rep(1,1000)
likelyhood = dbinom(6,9,prob=p_grid)
unst.posterior = likelyhood*prior
posterior = unst.posterior/sum(unst.posterior)
set.seed(100)
samples = sample(p_grid, prob=posterior, size=10^4, replace=T)
# 1
mean(samples < 0.2)
## [1] 5e-04
# 2
mean(samples > 0.8)
## [1] 0.1117
# 3
mean(samples > 0.2 & samples < 0.8)
## [1] 0.8878
# 4&5
quantile(samples, probs = c(0.2,0.8))
## 20% 80%
## 0.5195195 0.7567568
# 6
library(rethinking)
HPDI(samples, 0.5)
## |0.5 0.5|
## 0.5655656 0.7537538
PI(samples, 0.5)
## 25% 75%
## 0.5445445 0.7357357
Medium
3M1-2
# 1
p_grid = seq(0,1,length=1000)
prior=rep(1,1000)
likelyhood = dbinom(8,15,prob=p_grid)
unst.posterior = likelyhood*prior
posterior = unst.posterior/sum(unst.posterior)
# 2
set.seed(100)
samples = sample(p_grid, prob=posterior, size=10^4, replace=T)
HPDI(samples, 0.9)
## |0.9 0.9|
## 0.3243243 0.7157157
# 3
sim <- rbinom(10^4, 15, prob = samples)
mean(sim==8)
## [1] 0.1475
# 4
sim <- rbinom(10^4, 9, prob = samples)
mean(sim==6)
## [1] 0.1766
3M5
# 1
p_grid = seq(0,1,length=1000)
prior=rep(1,1000)
prior[p_grid<0.5]=0
likelyhood = dbinom(8,15,prob=p_grid)
unst.posterior = likelyhood*prior
posterior = unst.posterior/sum(unst.posterior)
# 2
set.seed(100)
samples = sample(p_grid, prob=posterior, size=10^4, replace=T)
HPDI(samples, 0.9)
## |0.9 0.9|
## 0.5005005 0.7077077
# 3
sim <- rbinom(10^4, 15, prob = samples)
mean(sim==8)
## [1] 0.1617
sim <- rbinom(10^4, 9, prob = samples)
mean(sim==6)
## [1] 0.2376
Hard
data(homeworkch3)
total.births <- length(birth1) + length(birth2)
boys.born <- sum(birth1 + birth2)
girls.born <- total.births - boys.born
p_grid = seq(0,1,length=1000)
prior=rep(1,1000)
likelyhood = dbinom(boys.born,total.births,prob=p_grid)
unst.posterior = likelyhood*prior
posterior = unst.posterior/sum(unst.posterior)
p_grid[which.max(posterior)]
## [1] 0.5545546
plot(y=posterior, x=p_grid, type="l")

samples = sample(p_grid, prob=posterior, size=10^4, replace=T)
HPDI(samples, c(0.5,0.89,0.97))
## |0.97 |0.89 |0.5 0.5| 0.89| 0.97|
## 0.4824825 0.4994995 0.5305305 0.5765766 0.6106106 0.6326326
# 3
sim <- rbinom(10^4, total.births, prob = samples)
dens(sim)
abline(v=boys.born)

# likelyhood = dbinom(sum(birth1),length(birth1),prob=p_grid)
# unst.posterior = likelyhood*prior
# posterior = unst.posterior/sum(unst.posterior)
# note: don't make a new model buth check the model with parts of the data
p_grid[which.max(posterior)]
## [1] 0.5545546
plot(y=posterior, x=p_grid, type="l")

samples = sample(p_grid, prob=posterior, size=10^4, replace=T)
HPDI(samples, c(0.5,0.89,0.97))
## |0.97 |0.89 |0.5 0.5| 0.89| 0.97|
## 0.4804805 0.5005005 0.5285285 0.5755756 0.6116116 0.6326326
# 3
sim <- rbinom(10^4, length(birth1), prob = samples)
dens(sim)
abline(v=sum(birth1))

boys.following.girls <- birth2[which(birth1==0)]
# likelyhood = dbinom(sum(boys.following.girls),length(boys.following.girls),prob=p_grid)
# unst.posterior = likelyhood*prior
# posterior = unst.posterior/sum(unst.posterior)
#
# p_grid[which.max(posterior)]
#
# plot(y=posterior, x=p_grid, type="l")
samples = sample(p_grid, prob=posterior, size=10^4, replace=T)
HPDI(samples, c(0.5,0.89,0.97))
## |0.97 |0.89 |0.5 0.5| 0.89| 0.97|
## 0.4794795 0.4944945 0.5305305 0.5765766 0.6066066 0.6306306
# 3
sim <- rbinom(10^4, length(boys.following.girls), prob = samples)
dens(sim)
abline(v=sum(birth1))
