number <- 5
skip <- 25
alag <- 100

library(MASS)
library(ts)

ps.options(horizontal=FALSE,onefile=FALSE)
ps.options(pagecentre=TRUE)
ps.options(paper="letter")
ps.options(width=7.0,height=10.0)
  
list <- c("01","02","03","04","05","06","07","08","09","10","11","12","13")

list <- c("01")

for (j in list) {

  run <- paste("0",j,"000",sep="")

  filename <- paste("pi",run,".dat",sep="")
  tmp <- scan(filename)
  rows <- tmp[1]
  cols <- tmp[2]

  pi <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)

  if (number > 0) {
    for (i in (1:number) ) {
      if (i < 10) {
        run <- paste("0",j,"00",i,sep="")
      } else {
        run <- paste("0",j,"0",i,sep="")
      }
      filename <- paste("pi",run,".dat",sep="")
      tmp <- scan(filename)
      rows <- tmp[1]
      cols <- tmp[2]
      pii <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
      pi <- cbind(pi,pii);
    }
  }

  idx <- seq(1,ncol(pi),skip);
  pi <-pi[idx];
}
  
for (j in list) { 

  run <- paste("0",j,"000",sep="") 
  
  filename <- paste("theta",run,".dat",sep="") 
  tmp <- scan(filename)
  rows <- tmp[1]
  cols <- tmp[2]
  
  post <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
  
  if (number > 0) {
    for (i in (1:number) ) {
      if (i < 10) {
        run <- paste("0",j,"00",i,sep="")
      } else {
        run <- paste("0",j,"0",i,sep="") 
      }
      filename <- paste("theta",run,".dat",sep="")
      tmp <- scan(filename)
      rows <- tmp[1]
      cols <- tmp[2]
      posti <- matrix(tmp[3:(2+rows*cols)],nrow=rows,ncol=cols)
      post <- cbind(post,posti);
    }
  }

  idx <- seq(1,ncol(post),skip);
  post <- post[,idx];

  post <- rbind(post,pi)

  rows <- rows + 1

  rm(tmp)

  tpost <- t(post);

  theta <- c("b1","b2","B11","B21","B12","B22","R11","R12","R22","pi")
  
  run <- paste("0",j,sep="") 
  
  dimnames(tpost) <- list(NULL,theta);
  tpost <- as.data.frame(tpost)
  print("  ")
  print(mean(tpost))
  print("  ")
  print(cor(tpost))
  
  print("  ")
  print("Minima and Maxima")
  for (j in (1:rows)) {
    print(paste(theta[j],"  ", min(tpost[,j]),"  ", max(tpost[,j]),sep=""))
  }
  print("  ")

  filename <- paste("thser",run,".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)
  {
    plot(idx,post[j,],type="n",ylab="")
      lines(idx,post[j,],lty="solid")
      title(theta[j])
  }
  dev.off()

  filename <- paste("thden",run,".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(post[j,],adjust=4)
    x <- den$x
    y <- den$y
    plot(x,y,type="n",ylab="",yaxt="n")
      lines(x,y,lty="solid")
      title(theta[j])
  }
  dev.off()

  rm(den,x,y)

  rles1 <- rows-1

  filename <- paste("thaut",run,".eps",sep="")
  postscript(file=filename)
  par(mfrow=c(rles1,1),mar=c(2.5,4,1.5,2)+0.1) # mar=c(bottom,left,top,right)
  for (j in 1:rles1)
  {
    aut <- acf(tpost[,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")
      title(theta[j])
  }
  dev.off()

  rm(aut,y,aidx)
 
  filename <- paste("thcor",run,".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(tpost[,1:rles1])
  dev.off()

# filename <- paste("thbiv",run,".eps",sep="") 
# postscript(file=filename) 
# par(mfrow=c(rows,rows),mar=c(0.0,0.0,0.0,0.0)+0.1)
# for (j in 1:rows)
# {  
#   for (i in 1:rows) 
#   {  
#      f2 <- kde2d(x=post[i,],y=post[j,],n=70)
#      persp(f2)
#   }
# }
# dev.off() 

}