#---------------------------------------------------------------------
#  Inferenza di p avendo osservato x successi in n
#  e previsione del numero di successi in una seguente misura
#  (ricordiamo che siamo partiti da una prior in p uniforme)   
#
#  Usa la distribuzione beta (vedi inf_p.R)
# 
# Nota: la maggior parte del codice dello script è per
#       fare plot, calcolare i 'riassunti' e stamparli.
#       Le istruzioni chiave sono messe in evidenza mediante 
#       i commenti '# <<---'
#
#  GdA giugno 2023
#----------------------------------------------------------------------

# nr. prove effettuate e nr si successi osservati
n <- 10
x <- 3

# per fare due plot nella finestra grafica
par(mfrow= c(2,1) )

# random 
N <- 1000000

# estraiamo i valori di p mediante la sua p.d.f. (-> Beta(x+1, (n-x)+1) )
rp <- rbeta(N, x+1, (n-x)+1) # <<---

titolo <- sprintf("Inferenza di p della binomiale da n = %d e x = %d", n, x)
hist(rp, nc=100, col='cyan', prob=TRUE, xlab='p', ylab='f ( p )',
     xlim=c(0, 1), main=titolo)
p <- seq(0, 1, len=101)
points(p, dbeta(p, x+1, (n-x)+1) , ty='l', lwd=2, col='blue')
abline(v=x/n, lty=2, col='red', lwd=2)
testo <- sprintf("E[p] = %.3f \n\nsigma(p) = %.3f",
                 mean(rp), sd(rp) )
text(3/4, 2.0, testo, cex=1.5 )

# predizione del nr di successi in una misura futura ('xf')
xf <- rbinom(N, n, rp)  # <<---

titolo <- sprintf("Probabilità del nr di successi su %d prove, avendone precedentemente osservati %d su %d", n, x, n)
barplot(table(xf)/N, col='magenta', xlab='x', ylab='f(x)',  main=titolo )
E.xf  <- mean(xf)
sd.xf <- sd(xf) 
testo <- sprintf("E[xf] = %.1f \n\nsigma(xf) = %.2f",
                 E.xf, sd.xf )
text(10, 0.15, testo, cex=1.5 )

# resettiamo a un solo plot nella finestra grafica 
par(mfrow= c(1,1) )