#class16-prep --checking prob. intervals
#----------
a = 11
b = 12
pbeta(.65,a,b) - pbeta(.3,a,b)
pbeta(.57,a,b) - pbeta(.42,a,b)

# We make the plot in RStudio and export it to class16-beta.pdf
x = seq(0,1,.01)
y = dbeta(x,a,b)
plot(x,y, type='l', col='green', xlab=expression(theta), ylab='', lwd=3, cex.axis=1.3, cex.lab=2)
lines(c(.42,.57), c(.5,.5), col='red', lw=4)
lines(c(.3,.65), c(.4,.4), col='magenta', lw=4)

#-----------
# Severe respiratory failure tables
#------------
prior.flat = rep(1/16,16)
prior.flat = matrix(prior.flat, nrow=4, ncol=4)
x = c(18,18,32,32)
prior.informed = matrix(x, nrow=4, ncol=4)
for (j in 1:4)
  prior.informed[j,] = x

theta.c = c(1,3,5,7)/8
theta.e = c(1,3,5,7)/8
ntreated.c = 10
ncured.c = 6
ntreated.e = 29
ncured.e = 28
lik.c = (theta.c^ncured.c)*(1-theta.c)^(ntreated.c-ncured.c)
lik.e = (theta.e^ncured.e)*(1-theta.e)^(ntreated.e-ncured.e)
lik.joint = matrix(nrow=4,ncol=4)
for (i in 1:4){
    for (j in 1:4){
      lik.joint[i,j] = lik.c[i]*lik.e[j]
    }
}
print(lik.e)
print(lik.c)
unnormpost.flat = prior.flat*lik.joint
post.flat = unnormpost.flat/sum(unnormpost.flat)
unnormpost.informed = prior.informed*lik.joint
post.informed = unnormpost.informed/sum(unnormpost.informed)

print(prior.flat)
print(prior.informed)
print(lik.joint)
print(post.flat)
print(post.informed)

probE_betterthan_C.flat = 0
probE_point6_better.flat = 0
probE_betterthan_C.informed = 0
probE_point6_better.informed = 0
for (i in 1:4) {
  for (j in 1:4) {
        if (theta.e[j] > theta.c[i]){
          probE_betterthan_C.flat = probE_betterthan_C.flat + post.flat[i,j]
          probE_betterthan_C.informed = probE_betterthan_C.informed + post.informed[i,j]
        }   
        if (theta.e[j] - theta.c[i] >= .6){
           probE_point6_better.flat = probE_point6_better.flat + post.flat[i,j]
           probE_point6_better.informed = probE_point6_better.informed + post.informed[i,j]
        }
  }
}
print(probE_betterthan_C.flat)
print(probE_point6_better.flat)
print(probE_betterthan_C.informed)
print(probE_point6_better.informed)