#Function to 
# 1. compute summary statistics
# 2. make trace plot of the last 1,000 sampled generations
# 3. make histograms of the posterior density from multiple chains
#for the output from MrBayes

#Arguments are
# data1-5 - names of matrices containing the data read in from MrBayes - must be created in advance
# string1 - in quotes, the name of data set being analyzed
# string2 - in quotes, the name of the variable being analyzed (these are used for axis labels)
# column - which column in data1-5 will be used in the analysis

#Copy files trace.ps and post_dens.ps to other file names before further analysis


diag <- function(data1,data2,data3,data4,data5,string1,string2,column,ylimit) {

#first compute means and variances after burn-in

cat("\n\nSummary Information for the",string1,"Dataset --",string2,"\n\n")
cat("Run     Mean   Variance    Min      Max \n")
cat(" 1  ",mean(data1[101:10001,column]),var(data1[101:10001,column]),min(data1[101:10001,column]),max(data1[101:10001,column]),"\n")
cat(" 2  ",mean(data2[101:10001,column]),var(data2[101:10001,column]),min(data2[101:10001,column]),max(data2[101:10001,column]),"\n")
cat(" 3  ",mean(data3[101:10001,column]),var(data3[101:10001,column]),min(data3[101:10001,column]),max(data3[101:10001,column]),"\n")
cat(" 4  ",mean(data4[101:10001,column]),var(data4[101:10001,column]),min(data4[101:10001,column]),max(data4[101:10001,column]),"\n")
cat(" 5  ",mean(data5[101:10001,column]),var(data5[101:10001,column]),min(data5[101:10001,column]),max(data5[101:10001,column]),"\n\n")

min1 = min(min(data1[9001:10001,column]),min(data2[9001:10001,column]),min(data3[9001:10001,column]),min(data4[9001:10001,column]),min(data5[9001:10001,column]))
max1 = max(max(data1[9001:10001,column]),max(data2[9001:10001,column]),max(data3[9001:10001,column]),max(data4[9001:10001,column]),max(data5[9001:10001,column]))
min2 = min(min(data1[101:10001,column]),min(data2[101:10001,column]),min(data3[101:10001,column]),min(data4[101:10001,column]),min(data5[101:10001,column]))
max2 = max(max(data1[101:10001,column]),max(data2[101:10001,column]),max(data3[101:10001,column]),max(data4[101:10001,column]),max(data5[101:10001,column]))
incr=(max2-min2)/30

# now make plots

#6-panel plots
postscript(file="plots.ps",horizontal=T)
par(mfrow=c(3,2))

#trace plot for last 1,000 collected generations
plot(c(9000,10000),c(min1,max1),type='n',xlab="Generation x 1,000",ylab=paste(string2),main=paste(string1),cex.lab=1.5,cex.axis=1.5,cex.main=1.8)
points(data1[9001:10001,1]/1000,data1[9001:10001,column],pch=1,col=1)
points(data2[9001:10001,1]/1000,data2[9001:10001,column],pch=2,col=2)
points(data3[9001:10001,1]/1000,data3[9001:10001,column],pch=3,col=3)
points(data4[9001:10001,1]/1000,data4[9001:10001,column],pch=4,col=4)
points(data5[9001:10001,1]/1000,data5[9001:10001,column],pch=5,col=5)
#legend(9001,min1+(max1-min1)/5,legend=c("Chain 1","Chain 2","Chain 3","Chain 4","Chain 5"),col=c(1:5),pch=c(1:5))

#histograms of posterior densities from the five runs
hist(data1[101:10001,column],breaks=c(min2,min2+incr,min2+2*incr,min2+3*incr,min2+4*incr,min2+5*incr,min2+6*incr,min2+7*incr,min2+8*incr,min2+9*incr,min2+10*incr,min2+11*incr,min2+12*incr,min2+13*incr,min2+14*incr,min2+15*incr,min2+16*incr,min2+17*incr,min2+18*incr,min2+19*incr,min2+20*incr,min2+21*incr,min2+22*incr,min2+23*incr,min2+24*incr,min2+25*incr,min2+26*incr,min2+27*incr,min2+28*incr,min2+29*incr,min2+30*incr),xlim=c(min2,max2),ylim=c(0,ylimit),xlab=paste(string2),main=paste(string1," - Chain 1"))
hist(data2[101:10001,column],breaks=c(min2,min2+incr,min2+2*incr,min2+3*incr,min2+4*incr,min2+5*incr,min2+6*incr,min2+7*incr,min2+8*incr,min2+9*incr,min2+10*incr,min2+11*incr,min2+12*incr,min2+13*incr,min2+14*incr,min2+15*incr,min2+16*incr,min2+17*incr,min2+18*incr,min2+19*incr,min2+20*incr,min2+21*incr,min2+22*incr,min2+23*incr,min2+24*incr,min2+25*incr,min2+26*incr,min2+27*incr,min2+28*incr,min2+29*incr,min2+30*incr),xlim=c(min2,max2),ylim=c(0,ylimit),xlab=paste(string2),main=paste(string1," - Chain 2"))
hist(data3[101:10001,column],breaks=c(min2,min2+incr,min2+2*incr,min2+3*incr,min2+4*incr,min2+5*incr,min2+6*incr,min2+7*incr,min2+8*incr,min2+9*incr,min2+10*incr,min2+11*incr,min2+12*incr,min2+13*incr,min2+14*incr,min2+15*incr,min2+16*incr,min2+17*incr,min2+18*incr,min2+19*incr,min2+20*incr,min2+21*incr,min2+22*incr,min2+23*incr,min2+24*incr,min2+25*incr,min2+26*incr,min2+27*incr,min2+28*incr,min2+29*incr,min2+30*incr),xlim=c(min2,max2),ylim=c(0,ylimit),xlab=paste(string2),main=paste(string1," - Chain 3"))
hist(data4[101:10001,column],breaks=c(min2,min2+incr,min2+2*incr,min2+3*incr,min2+4*incr,min2+5*incr,min2+6*incr,min2+7*incr,min2+8*incr,min2+9*incr,min2+10*incr,min2+11*incr,min2+12*incr,min2+13*incr,min2+14*incr,min2+15*incr,min2+16*incr,min2+17*incr,min2+18*incr,min2+19*incr,min2+20*incr,min2+21*incr,min2+22*incr,min2+23*incr,min2+24*incr,min2+25*incr,min2+26*incr,min2+27*incr,min2+28*incr,min2+29*incr,min2+30*incr),xlim=c(min2,max2),ylim=c(0,ylimit),xlab=paste(string2),main=paste(string1," - Chain 4"))
hist(data5[101:10001,column],breaks=c(min2,min2+incr,min2+2*incr,min2+3*incr,min2+4*incr,min2+5*incr,min2+6*incr,min2+7*incr,min2+8*incr,min2+9*incr,min2+10*incr,min2+11*incr,min2+12*incr,min2+13*incr,min2+14*incr,min2+15*incr,min2+16*incr,min2+17*incr,min2+18*incr,min2+19*incr,min2+20*incr,min2+21*incr,min2+22*incr,min2+23*incr,min2+24*incr,min2+25*incr,min2+26*incr,min2+27*incr,min2+28*incr,min2+29*incr,min2+30*incr),xlim=c(min2,max2),ylim=c(0,ylimit),xlab=paste(string2),main=paste(string1," - Chain 5"))
dev.off()

#Gelman & Rubin PSRF
data=cbind(data1[101:10001,column],data2[101:10001,column],data3[101:10001,column],data4[101:10001,column],data5[101:10001,column])
#cat("\n Gelman & Rubin PSRF:\n\n   50%      95%\n -------- --------\n",gpar(data,keep.all=T)$confshrink,"\n\n");
cat("\n\nGelman \& Rubin PSRF: ",gpar(data,keep.all=T)$confshrink,"\n\n");

}