#################################################################################
# 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 <- 10

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

# 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,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)

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

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