#----------------------------------------------
# Esempio di fit con JAGS/rjags
# (soliti dati della molla)
#
# GdA marzo 2011
#----------------------------------------------
# function pausa
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

library(rjags)

modello <- "molla.bug"

# dati
dati <- NULL
dati$m.molla <- 0.063
dati$m.disco <- 0.0789
dati$l       <- c(0.014, 0.032, 0.049, 0.066, 0.085, 0.103, 0.119, 0.137)
dati$t       <- c(0.501, 0.557, 0.624, 0.678, 0.728, 0.779, 0.813, 0.863)
dati$Nmin    <- 3
dati$pi2     <- pi^2

# inits
inits <- NULL
inits$c.l   <- 0
inits$m.l   <- 0
inits$tau.l <- 1000
inits$tau.t <- 1000
inits$m.t   <- 1

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

update(jm, 1000)

samples <- coda.samples(jm, c("c.l","m.l","sigma.l","c.t","m.t","sigma.t","k", "g"), n.iter=10000)

print(summary(samples))
plot(samples)

pausa()
# trasformiamo formato
samples.df <- as.data.frame(as.mcmc(samples))
print(cor(molla.df))

layout(matrix(c(1:2),c(2,1)))
hist(molla.df$g, nc=100, col='blue', xlab='g', prob=TRUE, main="inferenza g")
hist(molla.df$k, nc=100, col='red', xlab='k', prob=TRUE, main="inferenza k")
layout(1)