#---------------------------------------------------------------------
#  Inferenza di lambda avendo osservato n conteggi
#  e previsione del numero di conteggi in una seguente misura
#  (ricordiamo che siamo partiti da una prior in lambda uniforme)   
#
#  Usa la distribuzione gamma (vedi inf_lambda_gamma.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. conteggi osservati
n <- 10

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

# random 
N <- 1000000

# estraiamo i valori di lambda mediante la sua p.d.f. (-> Gamma(n+1, 1) )
rl <- rgamma(N, n+1, 1) # <<---

titolo <- sprintf("Inferenza di lambda di Poisson avendo osservato n = %d", n)
hist(rl, nc=100, col='cyan', prob=TRUE, xlab='lambda', ylab='f ( lambda )',
     xlim=c(0, max(rl)), main=titolo)
l <- seq(0, max(rl), len=101)
points(l, dgamma(l, n+1,1) , ty='l', lwd=2, col='blue')
testo <- sprintf("E[Lambda] = %.1f \n\nsigma(Lambda) = %.1f",
                 mean(rl), sd(rl) )
text(3/4*max(rl), 4/5*max(dgamma(l, n+1,1)), testo, cex=1.5 )

# predizione del nr di conteggi in una misura futura ('nf')
nf <- rpois(N, rl)  # <<---

titolo <- sprintf("Previsione del nr di conteggi, avendo precedentemente osservato n = %d", n)
# l'istogramma è ingannevole! 
#hnf <- hist(nf+0.1, nc=100, col='magenta', prob=TRUE, xlab='n', ylab='f(n)',
#     xlim=c(0, max(rl)), main=titolo )
barplot(table(nf)/N, col='magenta', xlab='n', ylab='f(n)',
        xlim=c(0, max(nf)), main=titolo )
E.nf  <- mean(nf)
sd.nf <- sd(nf) 
testo <- sprintf("E[n] = %.1f \n\nsigma(n) = %.1f",
                 E.nf, sd.nf )
text(E.nf+4*sd.nf, 4/5*max(table(nf)/N), testo, cex=1.5 )

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