#----------------------------------------------
# Inferenza di p di Bernoulli da una sequenza
# + previsione di successi in ulteriori prove
#
# GdA marzo 2011
#----------------------------------------------
# function pausa
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

library(rjags)

modello = "inf_lambda_pred_N.bug"  # file con il modello

dati <- NULL  # oggetto con i dati
dati$X  <- c(23, 12, 18, 18, 33, 25, 17, 18, 24, 12, 21, 18, 12, 25, 18, 17)
# in alternativs
# dati$X <- rpois(16, 20)

jm <- jags.model(modello, dati)  # definisce il modello

update(jm, 100)                  # burn in 

catena <- coda.samples(jm, c("lambda","Y"), n.iter=10000)  # sampling 

print(summary(catena))
plot(catena)

pausa()
# trasformiamo il formato della catena
catena.df <- as.data.frame( as.mcmc(catena) )

# plot a modo nostro
layout(matrix(1:2,2,1)) 
hist(catena.df$lambda, nc=100, prob=TRUE,
     #xlim=c(40,170),
     col='yellow',
     xlab='lambda', ylab='f(lambda)', main='Inferenza su lambda')

ty <- table(catena.df$Y)
barplot(ty/sum(ty), col='red', xlab='Y', ylab='f(Y)',
#         xlim=c(40,170),
        main=sprintf('Numero di conteggi futuri', dati$n1))
layout(1)

pausa()
# plot di correlazione
print(cor(catena.df$lambda, catena.df$Y))
plot(catena.df$lambda, catena.df$Y,
     xlab='lambda', ylab='Y', main='Correlazione y-p')


# Risultato analitico ------------------------------------
# consideriamo il numero di conteggi totali
X.tot <- sum(dati$X)
n.obs <- length(dati$X)

# inferenza su lambda.tot (nel tempo totale di misura) 
E.lambda.tot     <- X.tot + 1
sigma.lambda.tot <- sqrt(X.tot + 1)

# risultato scalato al tempo della singola misura
E.lambda     <- E.lambda.tot / n.obs
sigma.lambda <- sigma.lambda.tot / n.obs 

cat("\nRisultato analitico:\n")
cat( sprintf("E[lambda] = %.2f\n", E.lambda) )
cat( sprintf("sigma[lambda] = %.2f\n", sigma.lambda) )

# previsione di conteggi in misura futura:
E.X.f  <- E.lambda

# incertezza sui  conteggi in misura futura:
# teniamo conto (in modo approssimativo della convoluzione su lambda incerto!)
sigma.X.f <- sqrt( E.X.f + sigma.lambda^2)

cat( sprintf("\n\nE[X_fut] = %.2f\n", E.X.f) )
cat( sprintf("sigma[X_fut] = %.2f  ( > sqrt(E[X_fut] =  %.2f ))\n",
             sigma.X.f, sqrt(E.X.f)) )