
## Simulated t data
y<-read.table('data-generated-normal-coda.txt')

ymat<-array(y[, 2], dim=c(1000,5, 500))

t(ymat[2,,])->temp

plot(1:5, temp[1,1:5], ylim=c(30, 150), type='l')
for (i in 2:200){
	lines(1:5, temp[i,1:5])
	}

## write data into txt files

for (i in 1001:2000){
 t(ymat[i-1000,,])->temp
 filename<-paste('./Tdata/TdataN500I',i,'.txt', sep='')
 write.table(temp, filename, row.names=F, col.names=F)
}

source("http://nd.psychstat.org/_media/classdata/rtobugs.r")

for (i in 1001:2000){
 t(ymat[i-1000,,])->temp
 filename<-paste('./Tdata/TdataN500I',i,'bugs.txt', sep='')
 rtobugs(list(N=500, y=temp), filename)
}


## Simulated t data
y<-read.table('N:\\Private\\zzy\\research\\Bayesian analysis of longitudinal data\\Student t distribution\\Simulation\\data-generated-normal-coda.txt')

ymat<-array(y[, 2], dim=c(1000,5, 500))


## write data into txt files

for (i in 1001:2000){
 t(ymat[i-1000,,])->temp
 filename<-paste('TdataN500I',i,'.txt', sep='')
 write.table(temp, filename, row.names=F, col.names=F)
}

source("http://nd.psychstat.org/_media/classdata/rtobugs.r")

for (i in 1001:2000){
 t(ymat[i-1000,,])->temp
 filename<-paste('TdataN500I',i,'bugs.txt', sep='')
 rtobugs(list(y=temp), filename)
}


for (i in 1001:2000){
 filename<-paste('TdataN500I',i,'.txt', sep='')
 data<-read.table(filename)
 temp<-data[1:100,]
 outname<-paste('c:/zzy/student-t/temp/TdataN500I',i,'.txt', sep='')
 write.table(temp, outname, row.names=F, col.names=F)  

}


## process results from Mplus analysis
res<-scan('c:\\mplus\\all_results.txt')

results<-matrix(res, byrow=T, nrow=1000)
n<-dim(results)[1]
t.est<-c(1.5,6,1,2,0,.5)

avg<-apply(results[,1:6],2,mean)
bias<-avg-t.est
sd.est<-apply(results[,1:6],2,sd)
msd<-apply(results[,7:12],2,mean)
cvg<-rep(0,6)
power<-rep(0,6)
for (i in 1:6){
  cvg[i]<-sum( t.est[i]>results[,i]-1.96*results[,i+6] & t.est[i]0 | results[,i]+1.96*results[,i+6]<0 )/n
}

cbind(t.est, avg, bias,sd.est, msd, cvg, power)[c(2,3,4,6,5,1),]