#This first section of the code allows you to enter various parameters for the
#distribution and see what it looks like.

#Enter the values for the distribution
ml <- 4.5 #Most likely value (Make sure it's between min.1 and max.1)
conf <- 10 #Confidence
min.1 <- 3 #Minimum
max.1 <- 5 #Maximum

#Rescale and calculate the parameters for the beta distribution
mlRescale <- (ml - min.1) / (max.1 - min.1);
confRescale <- conf + 3

a <- mlRescale * (confRescale - 3) + 1
b <- confRescale - a - 1

#Plot the result
theta <- seq(.01, .99, len=400)
plot(min.1 + (max.1 - min.1) * theta, dbeta(theta, a, b), type = 'l',
  ylim = c(0, max(dbeta(theta, a, b))),
  xlab = "Control Strength", ylab = "Probability Density")
segments(c(min.1, max.1), c(0, 0), c(min.1 + .01, min.1 + (max.1 - min.1) * .99),
  c(dbeta(.01, a, b), dbeta(.99, a, b)))

################################################################################
#Generate pictures of many different control strength distributions

genCSpic <- function(ml.in, conf.in, min.in, max.in, printTF = T) {

  mlRescale <- (ml.in - min.in) / (max.in - min.in);
  confRescale <- conf.in + 3

  a <- mlRescale * (confRescale - 3) + 1
  b <- confRescale - a - 1

  #Plot the result
  theta <- seq(.01, .99, len=400)
  if (printTF)
    jpeg(paste("Figures", as.character(conf.in), "-m", as.character(ml.in), ".jpg", sep = ""))
  plot(min.in + (max.in - min.in) * theta, dbeta(theta, a, b) / (max.in - min.in), type = 'l',
    ylim = c(0, max(dbeta(theta, a, b) / (max.in - min.in))),
    xlab = "Control Strength", ylab = "Probability Density")
  segments(c(min.in, max.in), c(0, 0),
    c(min.in + (max.in - min.in) * .01, min.in + (max.in - min.in) * .99),
    c(dbeta(.01, a, b) / (max.in - min.in), dbeta(.99, a, b) / (max.in - min.in)))
  if (printTF)
    dev.off()

}

for (conf in c(0, 5, 10, 25, 50, 80)) {
  for (ml in c(45, 50, 55)) {
    genCSpic(ml, conf, 40, 60)
  }
}

#Look at a few pictures on the screen.
par(mfrow = c(2, 2))
for (conf in c(2, 5, 20, 50)) {
  genCSpic(50, conf, 40, 60, F)
}

################################################################################
#Find the percentage of the distribution that falls in the middle 20%

ml <- .5
min.1 <- 0
max.1 <- 1

count <- 1
pctstore <- rep(0, 8) #This holds the answer.
for (conf in c(0, 2, 5, 10, 20, 50, 80, 100)) {
  #Rescale and calculate the parameters for the beta distribution
  mlRescale <- (ml - min.1) / (max.1 - min.1);
  confRescale <- conf + 3

  a <- mlRescale * (confRescale - 3) + 1
  b <- confRescale - a - 1

  pctstore[count] <- pbeta(.6, a, b) - pbeta(.4, a, b)
  count <- count + 1
}