motif()
theta_seq(0,1,length=5000)
par(mfrow=c(4,5))   #examine beta densities
for(i in c(1,5,10,20)) {
  for(j in c(1,5,10,20,100)) plot(theta,dbeta(theta,i,j))
}

par(mfrow=c(2,3))   #sample posteriors
for(i in 1:6) {
  y_rbinom(1,10*i,.2)
  plot(theta,dbeta(theta,y+1,10*i-y+1))
}
#repeat

par(mfrow=c(1,1))   #quantiles of the beta
plot(theta,dbeta(theta,2,10))
abline(v=qbeta(0.5,2,10),col=2)
abline(v=qbeta(0.025,2,10),col=2)
abline(v=qbeta(0.975,2,10),col=2)

x_runif(5000)
y_qbeta(x,2,10)     #sampling from the CDF
hist(y,prob=T,nclass=25)
lines(theta,dbeta(theta,2,10))
lines(theta[theta<max(y)],dbeta(theta[theta<max(y)],2,10))

z_sort(y)           #compare simulated and theoretical quantiles
z[5000*c(0.025,0.975)]
qbeta(c(0.025,0.975),2,10)

c(mean(y),var(y))
c(2/(2+10),2*10/((2+10)^2*(2+10+1)))

c(306/651,300/705)

plot(theta,dbeta(theta,301,406),type="l",col=1)
lines(theta,dbeta(theta,307,346),col=2)

d_rbeta(5000,301,406)-rbeta(5000,307,346)
hist(d,nclass=30,prob=T)
summary(d)
sum(d>0)/5000
mean(d>0)
quantile(d,prob=c(0.05,0.95))