#---------------------------------------------------------------------
# inferenza di lambda avendo osservato n conteggi
# e assumendo una prior uniforme
# 
# nota: questo script a scopo didattico e funziona solo per n 'piccoli'
# 
# versione bis: 
# - f.lambda contiene già a denominatore n!
# - ma è mantenuto, come controprova, il calcolo di 'norm'    
#
#  GdA giugno 2023
#----------------------------------------------------------------------

# f.lambda 
# - Se chiamata senza 'norm' (default 1)  ritorna la funzione
#   non normalizzato (a meno che norm non valga proprio 1)
# - Se chiamata con norm e pow=0 (default) calcola la pdf
# - il parametro pow serve a calcolare valore attesi di potenze di lambda 
f.lambda <- function(l, n, norm=1, pow=0)  l^pow * exp(-l)*l^n/factorial(n) / norm

# per ogni n osservato fa l'inferenza e il plot (se richiesto)
inf.lambda <- function(n, l, plot=FALSE) {
    norm <- integrate(f.lambda, 0, Inf, n=n)$value
    # plot
    if(plot) points(l, f.lambda(l,n,norm), ty='l', col='blue')
    # 'riassunti'
    E.lambda   <- integrate(f.lambda, 0, Inf, n=n, norm=norm, pow=1)$value
    E.lambda2  <- integrate(f.lambda, 0, Inf, n=n, norm=norm, pow=2)$value
    Var.lambda <- E.lambda2 - E.lambda
    sigma.lambda  <- sqrt(Var.lambda)
    cat(sprintf("n = %d, norm: %.1f,  E[lambda] = %.1f,  Var[...] = %.1f, sigma[...] = %.1f\n",
                n, norm, E.lambda, Var.lambda, sigma.lambda ))
}

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

l.max <- 12
plot = TRUE
if (plot) {
    l <- seq(0, l.max, len=501)
    plot(NA, NA, ty='l', xaxs = "i", yaxs = "i" , xlim=range(l), ylim=c(0,1),         
         xlab='lambda', ylab='f(lambda)')
}

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

Sys.sleep(1)
cat("\n Alcuni valori di probabilità di lambda:\n")
# esempi di probabilità che lambda sia in un certo intervallo
n = 0; l.min=0; l.max=3
cat(sprintf("n=%d =>  P(%.0f <= lambda <= %.0f) = %.3f\n",
            n, l.min, l.max, prob.lambda(l.min, l.max, n) ))
n = 2; l.min=0; l.max=2
cat(sprintf("n=%d =>  P(%.0f <= lambda <= %.0f) = %.3f\n",
            n, l.min, l.max, prob.lambda(l.min, l.max, n) ))
n = 2; l.min=2; l.max=Inf
cat(sprintf("n=%d =>  P(%.0f <= lambda <= %.0f) = %.3f\n",
            n, l.min, l.max, prob.lambda(l.min, l.max, n) ))

n = 5; l.min=4; l.max=6
cat(sprintf("n=%d =>  P(%.0f <= lambda <= %.0f) = %.3f\n",
            n, l.min, l.max, prob.lambda(l.min, l.max, n) ))