Table of Contents

Process the results

Click to display ⇲

Click to hide ⇱ 

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<-'N50T5'

for (i in 1:1000){
  codafile<-paste("/pscratch/zzhang4/Private/LCSM/results/",stem,'/CODA-',i,".txt",sep='')
  
  if (file.exists(codafile)){
     codasamp<-read.table(codafile)

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

     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)
     cat(i, "\n")
  }
  file.res<-paste(stem,'.txt', sep='')
}
  write.table(res, file.res)

## final process of results
resname<-'N50T5.txt'
res<-read.table(resname,header=T,row.names=NULL)

res<-res[res[,62]<2, ]
   
N<-dim(res)[1]
par.name<-c('mu1','mu2','e2','v2','vls11','vls12','vls22','beta0','beta1','beta2')
par<-c(3.5,10,1,.1,2,0,.1,-.2,-.2,.1)
 
for (i in 1:10){
 ind<-c(i+1, i+11, i+21, i+31, 2*i+41, 2*i+42)
 temp<-res[,ind]
 
 avg<-mean(temp[,1])
 
 if (par[i]==0){
     bias<-(avg-par[i])*100
   }else{
     bias<-(avg-par[i])/par[i]*100
   }
 
 esd<-sd(temp[,1])
 msd<-mean(temp[,2])
 
 cvg1<-sum( temp[,3]<par[i] & temp[,4]>par[i] )/N
 cvg1.l<-sum(temp[,3]>par[i])/N
 cvg1.u<-sum(temp[,4]<par[i])/N
 
 cvg2<-sum( temp[,5]<par[i] & temp[,6]>par[i] )/N
 cvg2.l<-sum(temp[,5]>par[i])/N
 cvg2.u<-sum(temp[,6]<par[i])/N
 
 power1<-sum( temp[,3]>0 | temp[,4]<0 )/N
 power1.l<-sum(temp[,3]>0)/N
 power1.u<-sum(temp[,4]<0)/N
 
 power2<-sum( temp[,5]>0 | temp[,6]<0 )/N
 power2.l<-sum(temp[,5]>0)/N
 power2.u<-sum(temp[,6]<0)/N
 
 
 cat(par.name[i], c(avg, bias, msd, esd, cvg1, cvg1.l, cvg1.u, power1, power1.l, power1.u, cvg2, cvg2.l, cvg2.u, power2, power2.l, power2.u, N), '\n')
 }


####################
## final process of results
## 1
   2  3  4  5  6  7  8   9  10   
   v1 v2 v3 v4 v5 v6 v7 v8 v9 (par)
   11 12 13 14 15 16 17 18 19
   v1 v2 v3 v4 v5 v6 v7 v8 v9 (sd)
   20 21 22 23 24 25 26 27 28 
   v1 v2 v3 v4 v5 v6 v7 v8 v9 (p025)
   29 30 31 32 33 34 35 36 37
   v1 v2 v3 v4 v5 v6 v7 v8 v9 (p975)
   38 39 40 41 42 43 44 45 46 47   

resname<-'N50T4.txt'
res<-read.table(resname,header=T,row.names=NULL)

res<-res[res$V55<2, ]
   
N<-dim(res)[1]
P<-9
   
par.name<-c('mu1','mu2','e2','v2','vls11','vls12','vls22','beta0','beta1','beta2')
par<-c(10,-.2,1,.01,2,0,.05,-.2,.1)
 
for (i in 1:P){
 ind<-c(i+1, i+P+1, i+P*2+1, i+P*3+1, 2*i+P*4, 2*i+P*4+1)
 temp<-res[,ind]
 
 avg<-mean(temp[,1])
 
 if (par[i]==0){
     bias<-(avg-par[i])*100
   }else{
     bias<-(avg-par[i])/par[i]*100
   }
 
 esd<-sd(temp[,1])
 msd<-mean(temp[,2])
 
 cvg1<-sum( temp[,3]<par[i] & temp[,4]>par[i] )/N
 cvg1.l<-sum(temp[,3]>par[i])/N
 cvg1.u<-sum(temp[,4]<par[i])/N
 
 cvg2<-sum( temp[,5]<par[i] & temp[,6]>par[i] )/N
 cvg2.l<-sum(temp[,5]>par[i])/N
 cvg2.u<-sum(temp[,6]<par[i])/N
 
 power1<-sum( temp[,3]>0 | temp[,4]<0 )/N
 power1.l<-sum(temp[,3]>0)/N
 power1.u<-sum(temp[,4]<0)/N
 
 power2<-sum( temp[,5]>0 | temp[,6]<0 )/N
 power2.l<-sum(temp[,5]>0)/N
 power2.u<-sum(temp[,6]<0)/N
 
 
 cat(par.name[i], c(avg, bias, msd, esd, cvg1, cvg1.l, cvg1.u, power1, power1.l, power1.u, cvg2, cvg2.l, cvg2.u, power2, power2.l, power2.u, N), '\n')
 }

No split

Click to display ⇲

Click to hide ⇱ 

## Generate scripts to run the simulation

for (N in c(50,100)){
	for (T in c(4,5)){
		## create the analysis folder
		ana.folder<-paste('/pscratch/zzhang4/Private/LCSM/analysis/N',N,'T',T,sep='')
		dir.create(ana.folder)		
		setwd(ana.folder)
		
		## copy files to the folder
		file.copy('/pscratch/zzhang4/Private/LCSM/model/ranmodel.txt', 'model.txt', overwrite = TRUE)
		file.copy('/pscratch/zzhang4/Private/LCSM/model/raninit.txt', 'inits.txt', overwrite = TRUE)
		file.copy('/pscratch/zzhang4/Private/LCSM/model/ranscript.txt', 'script.txt', overwrite = TRUE)
		file.copy('/pscratch/zzhang4/Private/LCSM/model/ransub.sh', paste('N',N,'T',T,'.sh',sep=''), overwrite = TRUE)
		
		repfolder<-paste("sed 's/NSIM/",N,"/g' /pscratch/zzhang4/Private/LCSM/model/runsim.R > ransim0.R", sep='')
		system(repfolder)
		
		repfolder<-paste("sed 's/TSIM/",T,"/g' ransim0.R > runsim.R", sep='')
		system(repfolder)
		unlink('ransim0.R')
		ana.folder<-paste('/pscratch/zzhang4/Private/LCSM/results/N',N,'T',T,sep='')
		dir.create(ana.folder)		
	}
}


## submit jobs

for (N in c(50,100)){
	for (T in c(4,5)){
		## create the analysis folder
		ana.folder<-paste('/pscratch/zzhang4/Private/LCSM/analysis/N',N,'T',T,sep='')	
		setwd(ana.folder)
		
		system('qsub *.sh')	
	}
}

Split the analysis

Click to display ⇲

Click to hide ⇱ 

## Generate simulation scripts
for (N in c(50,100)){
	for (T in c(4,5)){
		for (K in 1:10){
			## create the analysis folder
			ana.folder<-paste('/pscratch/zzhang4/Private/LCSM/analysis/N',N,'T',T,'K',K,sep='')
			dir.create(ana.folder)		
			setwd(ana.folder)
		
			## copy files to the folder
			file.copy('/pscratch/zzhang4/Private/LCSM/model/ranmodel.txt', 'model.txt', overwrite = TRUE)
			file.copy('/pscratch/zzhang4/Private/LCSM/model/raninit.txt', 'inits.txt', overwrite = TRUE)
			file.copy('/pscratch/zzhang4/Private/LCSM/model/ranscript.txt', 'script.txt', overwrite = TRUE)
			file.copy('/pscratch/zzhang4/Private/LCSM/model/ransub.sh', paste('N',N,'T',T,'K',K,'.sh',sep=''), overwrite = TRUE)
		
			repfolder<-paste("sed 's/NSIM/",N,"/g' /pscratch/zzhang4/Private/LCSM/model/runsim.R > runsim0.R", sep='')
			system(repfolder)
		
			repfolder<-paste("sed 's/TSIM/",T,"/g' runsim0.R > runsim1.R", sep='')
			system(repfolder)
			
			repfolder<-paste("sed 's/startN/",(K-1)*100+1,"/g' runsim1.R > runsim2.R", sep='')
			system(repfolder)
			
			repfolder<-paste("sed 's/endN/",K*100,"/g' runsim2.R > runsim3.R", sep='')
			system(repfolder)
			
			repfolder<-paste("sed 's/KValue/",K,"/g' runsim3.R > runsim.R", sep='')
			system(repfolder)
			
			unlink('runsim0.R')
			unlink('runsim1.R')
			unlink('runsim2.R')
			unlink('runsim3.R')
		}
		ana.folder<-paste('/pscratch/zzhang4/Private/LCSM/results/N',N,'T',T,sep='')
		dir.create(ana.folder)		
	}
}


for (N in c(50,100)){
	for (T in c(4,5)){
	for (K in 1:10){
		## create the analysis folder
		ana.folder<-paste('/pscratch/zzhang4/Private/LCSM/analysis/N',N,'T',T,'K',K,sep='')	
		setwd(ana.folder)
		
		system('qsub *.sh')	
		}
	}
}