#---------------------------------------------------------------------
#  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
np <- 10
x <- 3

# per fare due plot nella finestra grafica
par(mfrow= c(3,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, (np-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 )

n = 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, np)
barplot(table(xf)/N, col='magenta', xlab='xf', ylab='f (xf) ',  main=titolo )
E.xf  <- mean(xf)
sd.xf <- sd(xf) 
testo <- sprintf("E[xf] = %.1f \n\nsigma(xf) = %.2f",
                 E.xf, sd.xf )
text(6, 0.22, testo, cex=2.5 )

n = 20 
# 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, np)
barplot(table(xf)/N, col='magenta', xlab='xf', ylab='f (xf)',  main=titolo )
E.xf  <- mean(xf)
sd.xf <- sd(xf) 
testo <- sprintf("E[xf] = %.1f \n\nsigma(xf) = %.2f",
                 E.xf, sd.xf )
text(18, 0.07, testo, cex=2.5 )

n = 40
# 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, np)
barplot(table(xf)/N, col='magenta', xlab='xf', ylab='f (xf)',  main=titolo )
E.xf  <- mean(xf)
sd.xf <- sd(xf) 
testo <- sprintf("E[xf] = %.1f \n\nsigma(xf) = %.2f",
                 E.xf, sd.xf )
text(35, 0.04, testo, cex=2.5 )


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