====== Process the results ======
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] )/N
cvg1.l<-sum(temp[,3]>par[i])/N
cvg1.u<-sum(temp[,4]par[i] )/N
cvg2.l<-sum(temp[,5]>par[i])/N
cvg2.u<-sum(temp[,6]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] )/N
cvg1.l<-sum(temp[,3]>par[i])/N
cvg1.u<-sum(temp[,4]par[i] )/N
cvg2.l<-sum(temp[,5]>par[i])/N
cvg2.u<-sum(temp[,6]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 =====
## 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 =====
## 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')
}
}
}