# un po' di utili istruzioni (anche ridondanti) per analizzare
# una catena risultante da un mcmc

library(coda)

# leggiamo la catena
gauss.msy <-read.coda("CODAchain1.txt","CODAindex.txt",quiet=TRUE)

# cos'e' gauss.msy?
str(gauss.msy)
is.mcmc(gauss.msy)
dimnames(gauss.msy)[[2]]

# a proposito: per vedere il contenuto di CODAindex.txt da R:
system("cat CODAindex.txt")
# e il risultato puo' essere messo in una lista
s <- system("cat CODAindex.txt", intern=TRUE)
s

# un'occhiata generale
summary(gauss.msy)
plot(gauss.msy)

# anche 
summary(gauss.msy[,1])
summary(gauss.msy[,2])
summary(gauss.msy[,3])

# una riga alla volta:
gauss.msy[1,]

# oppure
mean(gauss.msy[,1])
sd(gauss.msy[,1])

# istogrammi
layout(matrix(c(1:3),c(1,3)))
hist(gauss.msy[,1], nc=100)
hist(gauss.msy[,2], nc=100)
hist(gauss.msy[,3], nc=100)
layout(1)

# covarianza e coeff di corr:
cov(gauss.msy)
cor(gauss.msy)

# inoltre
sqrt(diag(cov(gauss.msy)))

# altro modo per guardare la catena
gauss.df <- as.dataframe(gauss.msy)

mean(gauss.df)
sd(gauss.df)

layout(matrix(c(1:3),c(1,3)))
hist(gauss.df$mu, nc=100)
hist(gauss.df$sigma, nc=100)
hist(gauss.df$y, nc=100)
layout(1)

cor(gauss.df)
var(gauss.df)


# cor() e var() si possono fare anche su gauss.msy (vedi sopra)
# ottenendo gli stessi risultati, ma mean() e sd() danno cose strane...