===== 2010 Sep 20 =====
<code>
geweke<-function(x, frac1 = 0.2, frac2 = 0.5){

	xstart <- c(start(x), end(x) - frac2 * (end(x) - start(x)))[c(1,3)]
   xend <- c(start(x) + frac1 * (end(x) - start(x)), end(x))[c(1,3)]

   y.variance <- y.mean <- vector("list", 2)
   for (i in 1:2) {
        y <- window(x, start = xstart[i], end = xend[i])
        y.mean[[i]] <- apply(as.matrix(y), 2, mean)
        y.variance[[i]] <- apply(as.matrix(y), 2, var)
    }
    z <- (y.mean[[1]] - y.mean[[2]])/sqrt(y.variance[[1]] + y.variance[[2]])
    z    
	}

HPD <- function(obj, prob = 0.95, ...)
{
    obj <- as.matrix(obj)
    vals <- apply(obj, 2, sort)
    if (!is.matrix(vals)) stop("obj must have nsamp > 1")
    nsamp <- nrow(vals)
    npar <- ncol(vals)
    gap <- max(1, min(nsamp - 1, round(nsamp * prob)))
    init <- 1:(nsamp - gap)
    inds <- apply(vals[init + gap, ,drop=FALSE] - vals[init, ,drop=FALSE],
                  2, which.min)
    ans <- cbind(vals[cbind(inds, 1:npar)],
                 vals[cbind(inds + gap, 1:npar)])
    ans
}


res<-NULL

stem<-'Tdata'

for (i in 1001:2000){
  codafile<-paste("/pscratch/zzhang4/Private/student-grm/CODA/TdataR/TdataCODAN300I",i,".txt",sep='')
  
  if (file.exists(codafile)){
     codasamp<-read.table(codafile)

     par<-matrix(c(codasamp[,2]),ncol=7)

     gew.max<- max(geweke(par))

     m.par<-apply(par,2,mean)
     sd.par<-apply(par,2,sd)
     p025.par<-apply(par,2,quantile,.025)
     p975.par<-apply(par,2,quantile,.975)
     
     hpd.ci<-c(t(HPD(par)))
     
     par.i<-c(m.par,sd.par,p025.par,p975.par,hpd.ci, gew.max,i)
     
     res<-rbind(res,par.i)
  }
  file.res<-paste('res-Tdata-T-N-300-Wishart.txt', sep='')

}

  write.table(res, file.res) 
</code>

  * T model results

<code>
## process the results
## 1) convergence diagnostic

## geweke function
geweke<-function(x, frac1 = 0.2, frac2 = 0.5){

	xstart <- c(start(x), end(x) - frac2 * (end(x) - start(x)))[c(1,3)]
   xend <- c(start(x) + frac1 * (end(x) - start(x)), end(x))[c(1,3)]

   y.variance <- y.mean <- vector("list", 2)
   for (i in 1:2) {
        y <- window(x, start = xstart[i], end = xend[i])
        y.mean[[i]] <- apply(as.matrix(y), 2, mean)
        y.variance[[i]] <- apply(as.matrix(y), 2, var)
    }
    z <- (y.mean[[1]] - y.mean[[2]])/sqrt(y.variance[[1]] + y.variance[[2]])
    z    
	}

HPD <- function(obj, prob = 0.95, ...)
{
    obj <- as.matrix(obj)
    vals <- apply(obj, 2, sort)
    if (!is.matrix(vals)) stop("obj must have nsamp > 1")
    nsamp <- nrow(vals)
    npar <- ncol(vals)
    gap <- max(1, min(nsamp - 1, round(nsamp * prob)))
    init <- 1:(nsamp - gap)
    inds <- apply(vals[init + gap, ,drop=FALSE] - vals[init, ,drop=FALSE],
                  2, which.min)
    ans <- cbind(vals[cbind(inds, 1:npar)],
                 vals[cbind(inds + gap, 1:npar)])
    ans
}


res<-NULL

stem<-'Tdata'

for (i in 1001:2000){
  codafile<-paste(stem,"CODAN500I",i,".txt",sep='')
  
  if (file.exists(codafile)){
     codasamp<-read.table(codafile)

     par<-matrix(c(codasamp[,2]),ncol=7)

     gew.max<- max(geweke(par))

     m.par<-apply(par,2,mean)
     sd.par<-apply(par,2,sd)
     p025.par<-apply(par,2,quantile,.025)
     p975.par<-apply(par,2,quantile,.975)
     
     hpd.ci<-c(t(HPD(par)))
     
     par.i<-c(m.par,sd.par,p025.par,p975.par,hpd.ci, gew.max,i)
     
     res<-rbind(res,par.i)
  }
  file.res<-paste('res-',stem,'.txt', sep='')

}

  write.table(res, file.res) 

</code>

====  Normal model results ====

<code>
## process the results
## 1) convergence diagnostic

## geweke function
geweke<-function(x, frac1 = 0.2, frac2 = 0.5){

	xstart <- c(start(x), end(x) - frac2 * (end(x) - start(x)))[c(1,3)]
   xend <- c(start(x) + frac1 * (end(x) - start(x)), end(x))[c(1,3)]

   y.variance <- y.mean <- vector("list", 2)
   for (i in 1:2) {
        y <- window(x, start = xstart[i], end = xend[i])
        y.mean[[i]] <- apply(as.matrix(y), 2, mean)
        y.variance[[i]] <- apply(as.matrix(y), 2, var)
    }
    z <- (y.mean[[1]] - y.mean[[2]])/sqrt(y.variance[[1]] + y.variance[[2]])
    z    
	}

HPD <- function(obj, prob = 0.95, ...)
{
    obj <- as.matrix(obj)
    vals <- apply(obj, 2, sort)
    if (!is.matrix(vals)) stop("obj must have nsamp > 1")
    nsamp <- nrow(vals)
    npar <- ncol(vals)
    gap <- max(1, min(nsamp - 1, round(nsamp * prob)))
    init <- 1:(nsamp - gap)
    inds <- apply(vals[init + gap, ,drop=FALSE] - vals[init, ,drop=FALSE],
                  2, which.min)
    ans <- cbind(vals[cbind(inds, 1:npar)],
                 vals[cbind(inds + gap, 1:npar)])
    ans
}


res<-NULL

stem<-'Tdata'

for (i in 1001:2000){
  
codafile<-paste("/pscratch/zzhang4/Private/student-grm/CODA/Tdata/TdataCODAN400I",i,".txt",sep=''
)
  
  if (file.exists(codafile)){
     codasamp<-read.table(codafile)

     par<-matrix(c(codasamp[,2]),ncol=6)

     gew.max<- max(geweke(par))

     m.par<-apply(par,2,mean)
     sd.par<-apply(par,2,sd)
     p025.par<-apply(par,2,quantile,.025)
     p975.par<-apply(par,2,quantile,.975)
     
     hpd.ci<-c(t(HPD(par)))
     
     par.i<-c(m.par,sd.par,p025.par,p975.par,hpd.ci, gew.max,i)
     
     res<-rbind(res,par.i)
  }
  file.res<-paste('res-Tdata-N-N-400-Uniform.txt', sep='')

}

  write.table(res, file.res) 

</code>