bayesian_HW4.R

## 4.1
theta1 <- rbeta(5000,58,44)
theta2 <- rbeta(5000,31,21)
mean(theta1 < theta2)

## [1] 0.6248

## 4.3
A = cbind(12,9,12,14,13,13,15,8,15,6)
B = cbind(11,11,10,9,9,8,7,10,6,8,8,9,7)
S = 10000
thetaA <- rgamma(S, 120+sum(A),10+length(A))
thetaB <- rgamma(S, 12+sum(B), 1+length(B))
mean(thetaB < thetaA)

## [1] 0.9947

n0 <- seq(1,200,10)
mean1 <- rep(NA,20)
for(n0 in seq(1,200,10)){
  thetaB = rgamma(S, 12*n0 + sum(B), n0+length(B))
  mean1[which(seq(1,200,10)==n0)] = mean(thetaB<thetaA)
}
mean1

##  [1] 0.9951 0.9481 0.8805 0.8106 0.7618 0.7194 0.6800 0.6565 0.6274 0.6101
## [11] 0.5960 0.5814 0.5686 0.5576 0.5522 0.5430 0.5312 0.5299 0.5199 0.5148