#---------------------------------------------------------------------
#  Inferenza di lambda avendo osservato n conteggi
#  e previsione del numero di conteggi in una seguente misura
#
#  Nota: il CODICE IMPORTANTE è quello segnalato da '# <<---'
#  (tutto il resto è per stampare/graficare i risultati)

#  GdA maggio 2025
#----------------------------------------------------------------------

# 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("nr di conteggi osservati:  %d", n)
hist(rl, nc=100, col='cyan', prob=TRUE, xlab='lambda', ylab='f ( lambda )',
     main=titolo)

# sovrapponiamo anche la curva della gamma
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(4/5*max(rl), 4/5*max(dgamma(l, n+1,1)), testo, cex=1.5 )


# predizione del nr di conteggi in una misura futura ('nf')
# si noti come il secondo parametro di rpois() è il vettore 'rl'
# => ogni elemento di 'nf' è estratto con il corrispondente elemento di 'rl' 
nf <- rpois(N, rl)  # <<---

barplot(table(nf)/N, col='magenta', xlab='n', ylab='f(n)', xlim=c(0, max(nf)) )
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 la finestra grafica a un singolo plot
par(mfrow= c(1,1) )