R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > invisible(options(echo = TRUE)) > prefix <- "sv" > first <- 1 > last <- 15 > stride <- 10 > alag <- 100 > library(MASS) > library(ts) Warning message: package 'ts' has been merged into 'stats' > > ps.options(horizontal=FALSE,onefile=FALSE) > ps.options(pagecentre=TRUE) > ps.options(paper="letter") > ps.options(width=7.0,height=10.0) > > for (i in (first:last)) { + if (i < 10) { + filename <- paste(prefix,".pi.00",i,".dat",sep="") + } else { + filename <- paste(prefix,".pi.0",i,".dat",sep="") + } + tmp <- scan(filename) + rows <- tmp[1] + cols <- tmp[2] + pii <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols) + if (i == first) { + pi <- pii + } else { + pi <- cbind(pi,pii) + } + } Read 6002 items Read 6002 items Read 6002 items Read 6002 items Read 6002 items Read 6002 items Read 6002 items Read 6002 items Read 6002 items Read 6002 items Read 6002 items Read 6002 items Read 6002 items Read 6002 items Read 6002 items > > idx <- seq(1,ncol(pi),stride) > pi <- pi[1,idx] > if (is.vector(pi)) pi <- matrix(pi,nrow=1,ncol=length(pi)) > > for (i in (first:last)) { + if (i < 10) { + filename <- paste(prefix,".rho.00",i,".dat",sep="") + } else { + filename <- paste(prefix,".rho.0",i,".dat",sep="") + } + tmp <- scan(filename) + rows <- tmp[1] + cols <- tmp[2] + rhoi <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols) + if (i == first) { + rho <- rhoi + print(dim(rho)) + } else { + rho <- cbind(rho,rhoi) + print(dim(rho)) + } + } Read 12002 items [1] 6 2000 Read 12002 items [1] 6 4000 Read 12002 items [1] 6 6000 Read 12002 items [1] 6 8000 Read 12002 items [1] 6 10000 Read 12002 items [1] 6 12000 Read 12002 items [1] 6 14000 Read 12002 items [1] 6 16000 Read 12002 items [1] 6 18000 Read 12002 items [1] 6 20000 Read 12002 items [1] 6 22000 Read 12002 items [1] 6 24000 Read 12002 items [1] 6 26000 Read 12002 items [1] 6 28000 Read 12002 items [1] 6 30000 > > idx <- seq(1,ncol(rho),stride) > rho <-rho[,idx] > if (is.vector(rho)) rho <- matrix(rho,nrow=1,ncol=length(rho)) > > for (i in (first:last)) { + if (i < 10) { + filename <- paste(prefix,".stats.00",i,".dat",sep="") + } else { + filename <- paste(prefix,".stats.0",i,".dat",sep="") + } + tmp <- scan(filename) + rows <- tmp[1] + cols <- tmp[2] + statsi <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols) + if (i == first) { + stats <- statsi + print(dim(stats)) + } else { + stats <- cbind(stats,statsi) + print(dim(stats)) + } + } Read 16002 items [1] 8 2000 Read 16002 items [1] 8 4000 Read 16002 items [1] 8 6000 Read 16002 items [1] 8 8000 Read 16002 items [1] 8 10000 Read 16002 items [1] 8 12000 Read 16002 items [1] 8 14000 Read 16002 items [1] 8 16000 Read 16002 items [1] 8 18000 Read 16002 items [1] 8 20000 Read 16002 items [1] 8 22000 Read 16002 items [1] 8 24000 Read 16002 items [1] 8 26000 Read 16002 items [1] 8 28000 Read 16002 items [1] 8 30000 > > idx <- seq(1,ncol(stats),stride) > stats <-stats[,idx] > if (is.vector(stats)) stats <- matrix(stats,nrow=1,ncol=length(stats)) > > rows <- nrow(rho) > > filename <- paste(prefix,".rho.chain.eps",sep="") > postscript(file=filename) > par(mfrow=c(rows+1,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(bottom,left,top,right) > for (j in 1:rows) { + plot(idx,rho[j,],type="n",ylab="") + lines(idx,rho[j,],lty="solid") + } > plot(idx,pi,type="n",ylab="") > lines(idx,pi,lty="solid") > dev.off() null device 1 > > trho <- t(rho) > trho <- as.data.frame(trho) > > if (nrow(rho) > 1) { + filename <- paste(prefix,".rho.pairs.eps",sep="") + postscript(file=filename) + par(mfrow=c(1,1),mar=c(4.5,4.0,1.5,1)+0.1,pty="s",xaxt="n",yaxt="n") + pairs(trho) + dev.off() + } null device 1 > > filename <- paste(prefix,".rho.acf.eps",sep="") > postscript(file=filename) > par(mfrow=c(rows,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(bottom,left,top,right) > for (j in 1:rows) { + aut <- acf(trho[,j],lag.max=alag,type="correlation",plot=F,demean=T) + y <- aut$acf + aidx <- 0:alag + plot(c(0,alag),c(0,1),type="n",ylab="") + lines(aidx,y,lty="solid") + } > dev.off() null device 1 > > filename <- paste(prefix,".rho.density.eps",sep="") > postscript(file=filename) > par(mfrow=c(rows,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(bottom,left,top,right) > for (j in 1:rows) { + den <- density(rho[j,],adjust=4) + plot(den$x,den$y,type="n",ylab="",yaxt="n") + lines(den$x,den$y,lty="solid") + } > dev.off() null device 1 > > > if (nrow(rho) > 1) { + corr <- cor(trho) + print(corr) + for (j in (1:rows)) { + for (i in (1:rows)) { + if (abs(corr[i,j]) < 0.6) {corr[i,j] <- 0.0} + } + } + print(corr) + } V1 V2 V3 V4 V5 V6 V1 1.00000000 0.028284722 -0.299880563 -0.17910825 0.134289247 0.051560030 V2 0.02828472 1.000000000 0.009935869 0.00762289 -0.170198578 -0.045095888 V3 -0.29988056 0.009935869 1.000000000 0.79180321 -0.602978957 0.008646002 V4 -0.17910825 0.007622890 0.791803209 1.00000000 -0.736369337 -0.048875707 V5 0.13428925 -0.170198578 -0.602978957 -0.73636934 1.000000000 -0.009771126 V6 0.05156003 -0.045095888 0.008646002 -0.04887571 -0.009771126 1.000000000 V1 V2 V3 V4 V5 V6 V1 1 0 0.0000000 0.0000000 0.0000000 0 V2 0 1 0.0000000 0.0000000 0.0000000 0 V3 0 0 1.0000000 0.7918032 -0.6029790 0 V4 0 0 0.7918032 1.0000000 -0.7363693 0 V5 0 0 -0.6029790 -0.7363693 1.0000000 0 V6 0 0 0.0000000 0.0000000 0.0000000 1 > proc.time() [1] 5.62 0.60 12.34 0.00 0.00 >