#----------------------------------------------
# Inferenza di p di Bernoulli da una sequenza
# + previsione di successi in ulteriori prove
# [variante con esempio di come si puo' inserire il modello
#  nel file R, in modo da avere tutto sott'occhio
#  (ovviamente va bene quando il codice per descrivere
#   il modello e' corto)]
#
# GdA marzo 2011
#----------------------------------------------
# function pausa
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

library(rjags)

source("dati_Bernoulli.R")             # file con la sequenza
modello = "temp_model.bug"  # file temporaneo con il modello

#------------  modello  ----------------------------
write("
model {
   for (i in 1:n) {    # nota:  JAGS ammette anche for (in 1:length(X))
     X[i] ~ dbern(p);
   }
   p ~ dbeta(1,1);
   y ~ dbin(p,n1);
}
", modello)
#----------------------------------------------------


dati <- NULL  # oggetto con i dati
dati$X  <- X
dati$n  <- length(X)  # facendo cosi' lo si puo' cambiare
                      # per eventuali inferenze parziali
dati$n1 <- 10

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

update(jm, 100)                  # burn in 

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

print(summary(catena))
plot(catena)

pausa()
# trasformiamo il formato della catena
catena.df <- as.data.frame( as.mcmc(catena) )
cat(sprintf("\n Matrice di correlazione: \n"))
print(cor(catena.df))

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

ty <- table(catena.df$y)
barplot(ty/sum(ty), col='red', xlab='y', ylab='f(y)',
        main=sprintf('Numero di successi in %d prove future', dati$n1))
layout(1)

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