#--------------------------------------------------
# Inferenza del numero di processi di Bernoulli
# dal numero di successi e dall'ipotesi su p
#
# GdA marzo 2011 (rivisto marzo 2012)
#--------------------------------------------------
# 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 = "inf_n.bug"  # file con il modello


dati <- NULL          # oggetto con i dati
dati$p  <- 0.2       # probabilita' nota a priori 
dati$s  <- as.integer(sum(X))     # numero di successi

inits <- NULL
inits$n <- max(dati$s, as.integer(dati$s/dati$p))   

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

update(jm, 100)                  # burn in 

catena <- coda.samples(jm, c("n"), 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
ty <- table(catena.df$n)
barplot(ty/sum(ty), col='red', xlab='n', ylab='f(n)',
        main='Inferenza sul numero di prove')