motif()
normal.1_scan("mcmc.normal.1.out")
plot(normal.1,pch=".")

normal.1_normal.1[1001:10000]
hist(normal.1)
qqnorm(normal.1)


normal.2_scan("mcmc.normal.2.out")
plot(normal.2,pch=".")
normal.2_normal.2[1001:10000]
hist(normal.2)
qqnorm(normal.2)

qqplot(normal.1,normal.2)
lines(normal.1,normal.1)

normal.3_scan("mcmc.normal.3.out")
plot(normal.3,pch=".")
normal.3_normal.3[1001:10000]
hist(normal.3)
qqnorm(normal.3)
qqplot(normal.2,normal.3)


normal.4_read.table("mcmc.normal.4.out")
par(mfrow=c(2,1))
plot(normal.4[,1],pch=".",main="mu")
plot(normal.4[,2],pch=".",main="psi")
normal.4_normal.4[1001:10000,]
par(mfrow=c(1,1))
hist(normal.4[,1])
hist(normal.4[,2])
qqnorm(normal.4[,1])
qqplot(normal.1,normal.4[,1])


#for the VA data
va.data_read.table("va92.orig.dat")
va_read.table("mcmc.va.out")
par(mfrow=c(2,1))
plot(va[,1],pch=".",main="alpha")
plot(va[,2],pch=".",main="beta")

va_va[1001:5000,]
par(mfrow=c(3,1))
plot(va[,1],pch=".",main="alpha")
plot(va[,2],pch=".",main="beta")
plot(va[,1]/(va[,1]+va[,2]),pch=".",main="alpha/beta")

#individual hospitals
va.mean_apply(va,2,mean)
options(object.size=5e8)
par(mfrow=c(3,3))
for(i in 3:11) {
  plot(va[,i],pch=".")
  abline(h=va.mean[i],col=2)
  abline(h=va.data[i-2,1]/va.data[i-2,3],col=6)
}