#----------------------------------------------
# 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_r_bck_measured.bug"  # file con il modello

dati <- NULL  # oggetto con i dati
dati$X <- 100   
dati$T <- 10
dati$TB <- 4   # tempo di misura del background
dati$XB <- 20  # numeri di background misurati

# Modo alternativo di fornire definire l'oggetto 'dati'
# [commentato dopo il casino fatto a lezione, 2/2/2024]
# dati=list(X=100, T=10, TB=4, XB=20)

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

update(jm, 100)                  # burn in 

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

print(summary(catena))
plot(catena)

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

layout(matrix(c(1:2),c(2,1)))
# plot a modo nostro
hist(catena.df$r, nc=100, prob=TRUE,
     col='yellow',
     xlab='r', ylab='f(r)', main='Inferenza su r, intensity of process')
hist(catena.df$rB, nc=100, prob=TRUE,
     col='red',
     xlab='rB', ylab='f(rB)', main='Inferenza su rB, intensity of backd process')

layout(1)
      
pausa()
plot(catena.df$rB, catena.df$r, col='blue')