betamu <- function (a,b) {
  a/(a+b)
   }

  betavar <- function(a,b) {
  a*b/((a+b)^2 * (a+b+1))
   }

####priors
  a<-c(7,14,28,56)
   b<-c(3,6,12,24)

round(betamu(a,b),5)

 round(sqrt(betavar(a,b)),5)

  round(pbeta(.5,a,b),5)


######posterior
  a<-c(92,99,113,141)
  b<-c(18,21,27,39)
   round(betamu(a,b),5)
  round(sqrt(betavar(a,b)),5)

####### Plotting priors and posteriors

par(mfrow=c(2,1))
 

plot(density(rbeta(10000,56,24)), lty=1, main = "Prior distributions")
lines(density(rbeta(10000,28,12)), lty=2, col='green')
 lines(density(rbeta(10000,14,6)), lty=3, col='blue')
 lines(density(rbeta(10000,7,3)), lty=4, col='red')

plot(density(rbeta(10000,141,39)))
 lines(density(rbeta(10000,113,27)), lty=2, col='green')
 lines(density(rbeta(10000,99,21)), lty=3, col='blue')
 lines(density(rbeta(10000,92,18)), lty=4, col='red')