#-----------------------------------------------------------------------
# inferenza di p della binomiale avendo osservato x successi in n prove
# e assumendo una prior uniforme
# 
# nota: questo script a scopo didattico e funziona solo per n 'piccoli'
# 
# Nota questo script è stato sulla traccia di inf_lambda_bis.R:
#  - la funzione (da normalizzare) ha la stessa espressione matematica
#     della binomiale (con diverso significato di p!)
#   - viene anche riportata, per confronto, la funzione
#     di probabilità Beta con parametri legati a n e x
#     (analoga della Gamma per la poissoniana)
#
#  GdA giugno 2023
#----------------------------------------------------------------------

f.p <- function(p, n, x, norm=1, pow=0)  {
    return( p^pow * factorial(n)/(factorial(x)*factorial(n-x)) *            
            p^x * (1-p)^(n-x)/ norm )    
}
# per ogni n e x l'inferenza e il plot (se richiesto)
inf.p <- function(n, x, p, plot=FALSE) {
    norm <- integrate(f.p, 0, 1, n=n, x=x)$value
    print(norm)
    if(plot) {
        points(p, f.p(p,n,x,norm), ty='l', col='blue')   # pdf calcolata 'a mano'
        points(p, dbeta(p, x+1, (n-x)+1), cex=0.8, col='red')
    }
    
    # 'riassunti'
    E.p   <- integrate(f.p, 0, 1, n=n, x=x, norm=norm, pow=1)$value
    E.p2  <- integrate(f.p, 0, 1, n=n, x=x, norm=norm, pow=2)$value
    Var.p <- E.p2 - E.p^2
    sigma.p  <- sqrt(Var.p)
    cat(sprintf("n=%d, x=%d: norm=%.4f, 1/norm=%.1f, E[p]=%.3f, Var[]=%.3f, sigma[]=%.3f\n",
                n, x, norm, 1/norm, E.p, Var.p, sigma.p ))
}

# calcolo della probabilità che p sia in un certo intervallo
prob.p <- function(p.min, p.max, n, x) {
    norm <- integrate(f.p, 0, 1, n=n, x=x)$value
    return( integrate(f.p, p.min, p.max, n=n, x=x, norm=norm)$value )    
}


n <- 10   # nr di prove

# plot vuoto, se richiesto 
plot = TRUE
if (plot) {
    p <- seq(0, 1, len=101)
     titolo = sprintf("inferenza di p con n=%d e x da %d a %d", n, 0, n)
    plot(NA, NA, ty='l', xaxs = "i", yaxs = "i" , xlim=range(p), ylim=c(0,n+1),         
         xlab='p', ylab='f (p)', main=titolo)
}

# inferenze per n prove e numero di successi da 0 a n
for(x in 0:n) {
    inf.p(n, x, p, plot)
    abline(v=x/n, lty=2, col='cyan')
    Sys.sleep(1)
}
cat("\n Quindi, qual'è l'espressione di f(p|n,x) con 'prior piatta'?\n\n")

Sys.sleep(1)
cat("\n Alcuni valori di probabilità di p\n")
x=3; p.min=0; p.max=0.3
cat(sprintf("n=%d, x=%d =>  P(%.2f <= p <= %.2f) = %.3f\n",
            n, x, p.min, p.max, prob.p(p.min, p.max, n, x) ))
cat(sprintf("n=%d, x=%d =>  P(%.2f <= p <= %.2f) = %.3f  [usando pbeta()]\n\n", n, x,
            p.min, p.max, pbeta(p.max, x+1, (n-x)+1) - pbeta(p.min, x+1, (n-x)+1)))

x=5; p.min=0; p.max=0.5
cat(sprintf("n=%d, x=%d =>  P(%.2f <= p <= %.2f) = %.3f\n",
            n, x, p.min, p.max, prob.p(p.min, p.max, n, x) ))
cat(sprintf("n=%d, x=%d =>  P(%.2f <= p <= %.2f) = %.3f  [usando pbeta()]\n\n", n, x,
            p.min, p.max, pbeta(p.max, x+1, (n-x)+1) - pbeta(p.min, x+1, (n-x)+1)))

x=5; p.min=0.2; p.max=0.8
cat(sprintf("n=%d, x=%d =>  P(%.2f <= p <= %.2f) = %.3f\n",
            n, x, p.min, p.max, prob.p(p.min, p.max, n, x) ))
cat(sprintf("n=%d, x=%d =>  P(%.2f <= p <= %.2f) = %.3f  [usando pbeta()]\n\n", n, x,
            p.min, p.max, pbeta(p.max, x+1, (n-x)+1) - pbeta(p.min, x+1, (n-x)+1)))