# 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)