beta.bestfit <- function(mode, p.up, x.up, a.step = 0.005){
# This function returns the a and b parameters 
# for a beta prior given the mode and upper p percentile.
# Since a mode is assumed a > 1 and b > 1 is assumed.

p.up <- round(p.up, digits = 3)

if(mode <= 0.5){

a <- 1

i <- 1
repeat{
a <- a + a.step
b <- (a*(1-mode)-1+2*mode)/mode
percent <- pbeta(x.up,a, b)
if(percent >= p.up){
a.final <- a
break
}
i <- i + 1
}
b.final <- (a.final*(1-mode)-1+2*mode)/mode
}else{ 
mode <- 1 - mode # fit the beta for a small a and then use symmetry
p.up <- 1 - p.up
x.up <- 1 - x.up

a <- 1

i <- 1
repeat{
a <- a + a.step
b <- (a*(1-mode)-1+2*mode)/mode
percent <- pbeta(x.up,a, b)
if(percent >= p.up){
a.final <- a
break
}
i <- i + 1
}
b.final <- (a.final*(1-mode)-1+2*mode)/mode

temp <- a.final # switch parameters
a.final <- b.final
b.final <- temp

}

final <- c(a.final, b.final)
return(final)
}