#Generate posterior plots and histograms

#Generate a postcript file
postscript(file='moreh4.ps',horizontal=F,width=8.0,height=10)

#One figure on one page
par(mfrow=c(1,1), las=1.0)

#Read in the output from c/c++ programs
taodata<-scan("moreh4.dat", what=list(talpha=0, tbeta=0, tao=0, salpha=0, sbeta=0, sigma=0))

m<-999
x<-seq(1, m)
r<-300
phistogram<-rep(NA, r)
y1<-seq(1:r)/15
y2<-seq(1:m)/50

#Create histograms fot tao^2
for (i in (1:r)) {
  phistogram[i]<-length(taodata$tao[taodata$tao>(i-1)/15 & taodata$tao<=i/15])*15/m
  }

pdensity<-rep(NA, m)

#Create mean posterior histograms for tao^2
for (i in (1:m)) {
  pdensity[i]<-mean(dgamma(taodata$tbeta/(i/50),taodata$talpha)*taodata$tbeta)/(i/50)/(i/50)
  }

phistsigma<-rep(NA, r)
y3<-seq(1:r)/50  
y4<-seq(1:m)/150

#Create histograms for sigma^2
for (i in (1:r)) {
  phistsigma[i]<-length(taodata$sigma[taodata$sigma>(i-1)/50 &
  taodata$sigma<=i/50])*50/m
  }
 
pdensigma<-rep(NA, m)

#Create mean posterior histograms for sigma^2
for (i in (1:m)) {
  pdensigma[i]<-mean(dgamma(taodata$sbeta/(i/150),taodata$salpha)*taodata$sbeta)/(i/150)/(i/150)
  }

#plot(x, taodata$tao, type="l", ylab=" ", xlab=" ", ylim=c(0, 80), cex=0.7)

#Generate plots with no x-axis and y-axis labels but scales are twice as big as normal
plot(y1, phistogram, type="l", xlab =" ", ylab=" ", xlim=c(0, 20), ylim=c(0, 0.4), cex=2)
plot(y2, pdensity, type="l", xlab =" ", ylab=" ", xlim=c(0,20), cex=2)

#plot(x, taodata$sigma, type="l", ylab=" ", xlab=" ", ylim=c(0, 6), cex=0.7)
plot(y3, phistsigma, type="l", xlab =" ", ylab=" ", xlim=c(0, 7), ylim=c(0,1.4), cex=2)  
plot(y4, pdensigma, type="l", xlab =" ", ylab=" ", xlim=c(0, 7), cex=2)