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

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

shape <- 1

i <- 1
repeat{
shape <- shape + shape.step
rate <- (shape - 1)/ mode
scale <- 1/rate
percent <- pgamma(x.up,shape, scale)
if(percent >= p.up){
shape.final <- shape
break
}
i <- i + 1
}
rate.final <- (shape.final - 1)/mode
scale.final <- 1/rate.final
final <- c(shape.final, rate.final, scale.final)
return(final)
}