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 = shape, scale = 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)
}

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)
}

normal.bestfit <- function(mode, p.up, x.up, sd.step = 0.0001){
    # This function returns the mean and standard deviation parameters 
    # for a normal prior given the mode and upper p percentile.
    
    p.up <- round(p.up, digits = 3)

    mean <- mode
    sd <- 0

    i <- 1
    repeat{
        sd <- sd + sd.step
        percent <- pnorm(x.up,mean,sd)
        if(percent <= p.up){
            sd.final <- sd
            break
        }
        i <- i + 1
    }
    mean.final <- mean
    final <- c(mean.final, sd.final)
    return(final)
}

#################################################################################
# mean and sd for beta and gamma

bm <- function(a,b) a/(a+b)
bs <- function(a,b) sqrt( (a*b)/((a+b)^2*(a+b+1)))

gm <- function(a,b) a/b
gs <- function(a,b) sqrt(a/b^2) 


#################################################################################
# Example of gamma.bestfit

mode <- 4

p <- 0.975

x.up <- 5

final <- gamma.bestfit(mode,p,x.up)
final

# Check

shape <- final[1]
scale <- final[3]

prior.p <- pgamma(x.up, shape, scale)
prior.p

x <- seq(0,(gm(shape,1/scale)+5*gs(shape,1/scale)),0.01)
y <- dgamma(x,shape = shape,scale = scale)
plot(x,y,type="l")
    

#################################################################################
# Example for beta.bestfit to left

mode <- 0.01

p <- 0.95

x.up <- 0.02

final <- beta.bestfit(mode,p,x.up)
final

# Check

a <- final[1]
a
b <- final[2]
b

prior.p <- pbeta(x.up,a,b)
prior.p

x <- seq(0,min((bm(a,b)+6*bs(a,b)),1),0.0001)
y <- dbeta(x,a,b)
plot(x,y,type="l")

#################################################################################
# Example for beta.bestfit to right

mode <- 0.99

p <- 0.05

x.up <- 0.98

final <- beta.bestfit(mode,p,x.up)

# Check

a <- final[1]
a
b <- final[2]
b

prior.p <- pbeta(x.up,a,b)
prior.p

x <- seq(max((bm(a,b)-6*bs(a,b)),0),1,0.0001)
y <- dbeta(x,a,b)
plot(x,y,type="l")

#################################################################################
# Example for normal.bestfit to left

mode <- 0.2

p <- 0.95

x.up <- 0.201

final <- normal.bestfit(mode,p,x.up)
final

# Check

mu <- final[1]
mu
sigma <- final[2]
sigma

prior.p <- pnorm(x.up,mu,sigma)
prior.p

x <- seq((mu-5*sigma),(mu+5*sigma),0.0001)
y <- dnorm(x, mu, sigma)
plot(x,y,type="l")