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