#-------------------------------------------------------------------
# script per generare ed analizzare campioni gaussiani
# affetti da errori sistematici sia di scala che di zero
#   
#  G. D'Agostini, gennaio 2011
#-------------------------------------------------------------------
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

library(rjags)

# simula i dati --------------------------------------------------
# tre valori veri
n1      <- 100   
mu1.T   <- 10.0  
n2      <- 50   
mu2.T   <- 15.0 
n3      <- 20   
mu3.T   <- 20.0 
sigma   <- 2    
# sistematiche
sigma.z <- 0.15
sigma.f <- 0.015

z <- rnorm(1, 0, sigma.z)
f <- rnorm(1, 1, sigma.f)
cat(sprintf("\nvalori generati delle sistematiche: z = %f;  f = %f\n\n",z, f))

dati <- NULL
dati$X1 <- rnorm(n1, f * mu1.T + z, sigma * f)
dati$X2 <- rnorm(n2, f * mu2.T + z, sigma * f)
dati$X3 <- rnorm(n3, f * mu3.T + z, sigma * f)
dati$sigma.z  <- sigma.z
dati$sigma.f  <- sigma.f
  
cat( sprintf("Campione 1: n1 = %d; media = %f, sigma = %f (sigma/sqrt(n) = %f)\n",
             n1, mean(dati$X1), sd(dati$X1), sd(dati$X1)/sqrt(n1) ))
cat( sprintf("Campione 2: n2 = %d; media = %f, sigma = %f (sigma/sqrt(n) = %f)\n",
             n2, mean(dati$X2), sd(dati$X2), sd(dati$X2)/sqrt(n2) ))
cat( sprintf("Campione 1: n1 = %d; media = %f, sigma = %f (sigma/sqrt(n) = %f)\n\n",
             n3, mean(dati$X3), sd(dati$X3), sd(dati$X3)/sqrt(n2) ))

# usa rjags -----------------------------------------------------
jm <- jags.model("norm-syst.bug", dati)

update(jm, 1000)

samples <- coda.samples(jm, c("mu1","mu2","mu3","sigma", "Delta.mu.21", "Delta.mu.32"), n.iter=50000)

# analizza la catena --------------------------------------------
print(summary(samples))
plot(samples)

# cambiamo formato
samples.df <- as.data.frame(as.mcmc(samples))
print(cor(samples.df))

pausa()
layout(matrix(c(1:3),c(3,3)))
hist(samples.df$mu1, nc=100, col='red', xlim=c(8,23))
hist(samples.df$mu2, nc=100, col='red', xlim=c(8,23))
hist(samples.df$mu3, nc=100, col='red', xlim=c(8,23))
layout(1)

pausa()
layout(matrix(c(1:3),c(3,3)))
hist(samples.df$sigma, nc=100, col='blue', xlim=c(1.5,2.5))
hist(samples.df$Delta.mu.21, nc=100, col='yellow', xlim=c(3,7))
hist(samples.df$Delta.mu.32, nc=100, col='yellow', xlim=c(3,7))
layout(1)


# etc. etc.