#Function to generate alpha and beta when mean and standard deviation of the distribution is known

Beta_Parameters <- function(mean, variance) {
  alpha <- ((1 - mean) / variance - 1 / mean) * mean ^ 2
  beta <- alpha * (1 / mean - 1)
  return(Beta_Parameters = list(alpha = alpha, beta = beta))
}

#Generate alpha and beta for beta distributions priors for the three ads (A,B, and C)
a<-Beta_Parameters(0.15,0.0016)
b<-Beta_Parameters(0.16,0.0016)
c<-Beta_Parameters(0.17,0.0016)

#Plot prior beta distibutions for the three ads (A,B, and C)
x=seq(0,0.4,length=10000)
y=dbeta(x,a$alpha,a$beta)
y1=dbeta(x,b$alpha,b$beta)
y2=dbeta(x,c$alpha,c$beta)
plot(x,y, type="l", col="blue",ylim = c(0,12),xlim=c(0,.40),lwd=2,ylab='Density',xlab= expression(theta))
lines(x,y1, col="orange",lwd=2)
lines(x,y2,col="black",lwd=2)
legend(0.1,12,legend = c("A", "B", "C"),col=c("blue","orange","black"), lwd=5, cex=1, horiz = TRUE)

#Simulate data to identify probability of ads A better than B with the prior distibutions
a_simulation=rbeta(10^6,a$alpha,a$beta)
b_simulation=rbeta(10^6,b$alpha,b$beta)
c_simulation=rbeta(10^6,c$alpha,c$beta)
mean(b_simulation>a_simulation)
mean(c_simulation>b_simulation)
mean(c_simulation>a_simulation)

########### 1st hour ##############################
#Results from the campaign in the first hour 

a1=4
an1=58
b1=14
bn1=101
c1=23
cn1=129

#Plot posterior beta distibutions for the three ads (A,B, and C) afterincorporating evidence from the first hour
x=seq(0,0.4,length=10000)
y=dbeta(x,a$alpha+a1,a$beta+an1-a1)
y1=dbeta(x,b$alpha+b1,b$beta+bn1-b1)
y2=dbeta(x,c$alpha+c1,c$beta+cn1-c1)
plot(x,y, type="l", col="blue",ylim = c(0,18),xlim=c(0,.40),lwd=2,ylab='Density',xlab= expression(theta))
lines(x,y1, col="orange",lwd=2)
lines(x,y2,col="black",lwd=2)
legend(0.1,18,legend = c("A", "B", "C"),col=c("blue","orange","black"), lwd=5, cex=1, horiz = TRUE)

#Simulate data to identify probability of ads A better than B with the prior distibutions
a_simulation=rbeta(10^6,a$alpha+a1,a$beta+an1-a1)
b_simulation=rbeta(10^6,b$alpha+b1,b$beta+bn1-b1)
c_simulation=rbeta(10^6,c$alpha+c1,c$beta+cn1-c1)
mean(b_simulation>a_simulation)
mean(c_simulation>b_simulation)
mean(c_simulation>a_simulation)

######### Complete Experiment ##############
#The final results from the campaign after the 5th day
a1=41
an1=554
b1=98
bn1=922
c1=230
cn1=1235

#Plot posterior beta distibutions for the three ads (A,B, and C) afterincorporating the final evidence from the campaign
x=seq(0,0.4,length=10000)
y=dbeta(x,a$alpha+a1,a$beta+an1-a1)
y1=dbeta(x,b$alpha+b1,b$beta+bn1-b1)
y2=dbeta(x,c$alpha+c1,c$beta+cn1-c1)

#Simulate data to identify probability of ads A better than B with the prior distibutions
plot(x,y, type="l", col="blue",ylim = c(0,46),xlim=c(0,0.4),lwd=2,ylab='Density',xlab= expression(theta))
lines(x,y1, col="orange",lwd=2)
lines(x,y2,col="black",lwd=2)
legend(0.085,46.5,legend = c("A", "B", "C"),col=c("blue","orange","black"), lwd=5, cex=1, horiz = TRUE)

a_simulation=rbeta(10^6,a$alpha+a1,a$beta+an1-a1)
b_simulation=rbeta(10^6,b$alpha+b1,b$beta+bn1-b1)
c_simulation=rbeta(10^6,c$alpha+c1,c$beta+cn1-c1)
mean(b_simulation>a_simulation)
mean(c_simulation>b_simulation)
mean(c_simulation>a_simulation)