#---------------------------------------------- # 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) source("dati_Bernoulli.R") # file con la sequenza modello = "inf_p_pred.bug" # file con il 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')