# questo non e' un vero script da eseguire, ma una serie di
# comandi da eseguire in successione

library(coda)

syst <-read.coda("bugs.out","bugs.ind",quiet=TRUE)
summary(syst)

layout(matrix(c(1:6),c(2,3)))
traceplot(syst)

plot(syst)   # (si autoridefinisce il layout)

syst.m <- as.matrix(syst)
mu1 <- syst.m[,1]
mu2 <- syst.m[,2]
sigma1 <- syst.m[,3]
sigma2 <- syst.m[,4]
y1 <- syst.m[,5]
y2 <- syst.m[,6]

layout(matrix(c(1:2),c(1,2)))
hist(mu1,nc=100,freq=F)
hist(mu2,nc=100,freq=F)
hist(sigma1,nc=100,freq=F)
hist(sigma2,nc=100,freq=F)
hist(y1,nc=100,freq=F)
hist(y2,nc=100,freq=F)
layout(1)

mysum <- function(x) {
    print(summary(x))
    media <- mean(x)
    stdev <- sd(x)
    print( sprintf("media = %f;   stddev = %f", media, stdev), quote=F )
}
corr <- function(x1,x2) {
    x1.m <- mean(x1)
    x2.m <- mean(x2)
    x1.sd <- sd(x1)
    x2.sd <- sd(x2)
    cov.x1.x2 <- mean(x1*x2)-x1.m*x2.m
    rho.x1.x2 <- cov.x1.x2 / (x1.sd * x2.sd)
    print( sprintf("cov = %f;  rho = %f", cov.x1.x2, rho.x1.x2), quote=F )
}
# in realta', esistono anche le funzioni di R cov() e cor()
# che fanno la stessa cosa

mysum(mu1)
mysum(mu2)
corr(mu1,mu2)
# nota: 
#  sqrt(cov.mu1.mu2) # viene 0.3, come deve.

mysum(sigma1)
mysum(sigma2)
corr(sigma1,sigma2)

mysum(y1)
mysum(y2)