#---------------------------------------------------------
# Inferenza di p del processo di Bernoulli (-> binomiale)
# dati il nr di prove e di successi
#
# GdA maggio 2025
#---------------------------------------------------------

dati <- as.numeric( unlist(
    strsplit( readline(prompt="dai nr prove e successi separati da virgola\n"), ',') ) )

n <- dati[1]
x <- dati[2]

cat(sprintf("\nDati:  %d prove con %d successi\n", n, x))

# prior deella Beta
ri <- 1  # alpha iniziale
si <- 1  # beta iniziale
# posterior 
r <- ri + x
s <- si + (n-x)
cat(sprintf("\n     ->  parametri Beta():  r = %d,  s = %d \n", r, s))
# summary di f(p)
E.p    <- r / (r+s)
var.p  <- r*s / ( (r+s+1)*(r+s)^2 )
sd.p   <- sqrt(var.p)
cat(sprintf("\n     ->  E(p) = %.3f, sd(p) = %.3f\n", E.p, sd.p))
pp <- seq(0,1, len=101)
plot(pp, dbeta(pp, r, s), ty='l', col='blue',
     xlab='p', ylab='f(p)', xaxs = "i", yaxs = "i")