#---------------------------------
#  Beta prior, given mu and sigma
#
#  GdA 15 feb 2025
#---------------------------------

# function pausa
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

mu    = 0.95
sigma = 0.05
#sigma = 0.025

# Wikipedia (https://en.wikipedia.org/wiki/Beta_distribution)
r <- mu * ( mu*(1-mu)/sigma^2 - 1 )
s <- r * (1-mu) / mu
# OK!

# r=96; s=6

print(r)
print(s)

# s <- max(1.1, s)

cat(sprintf("\n =>  r = %4f,  s = %4f \n", r, s))

x <- seq(0,1,len=501) 
plot(x, dbeta(x, r, s), ty='l', lwd=3, col='cyan', xlim=c(0.8,1))
grid()
pausa()

mu.r <- r /(r+s)
sigma.r <- sqrt( r*s / ((r+s+1)*(r+s)^2) ) 

cat(sprintf("\n Check:  =>  mu = %3f,  sigma = %3f\n\n",
            mu.r, sigma.r) )

n <- 100000
p <- rbeta(n, r, s)
hist(p, n =100, col='cyan', prob=TRUE)
print(summary(p))
cat(sprintf("\n std=  %.4f\n", sd(p)))