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')
}
