#Plot some summary graphs and print some summary statistics
#excluding the specified burnin, and possibly 'thinning' the sample

plot.mcmc = function (samples, ylab=samplename, burnin=0, thin=1,mu.target=0, sig.target=1) {

  #default value for ylab is just the name of the vector 'samples'
  samplename=paste(deparse(substitute(samples)))
  #samples is either an n-vector (from Metropolis)
  #or a nx2 matrix (from Gibbs)
  if (is.vector(samples))  {
    p=1
    n=length(samples)
    #make this a matrix so we can treat both in the same way
    samples=matrix(samples,n,p)
  }
  else {
    p=dim(samples)[2]
    n=dim(samples)[1]
  }

  #i selects which samples to summarise
  i=seq(burnin+1,n,by=thin)

  #Set up plot area (depends on how many samples we have to plot) 
  par(mfrow=c(p,3))
  for (j in 1:p) {
    #produce plots for each sample
    plot(c(min(i),max(i)),
         c(min(min(samples[i,j]),-3),max(max(samples[i,j]),3)),
         xlab="iterations",ylab=ylab[j],type="n")

    lines(i,samples[i,j],col="red")
    lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
    lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")
    plot(density(samples[i,j]),col="red",main=paste("Density of ",ylab[j]))
    curve(dnorm(x,mu.target[j],sqrt(sig.target[j])),add=T,col="blue",lty=2)
    acf(samples[i,j],main=paste("Autocorrelation for ",ylab[j]))

    cat("\n",ylab[j],":\n")
    cat("\tSamples:\t",length(i),"\n")
    cat("\tMean:\t\t",mean(samples[i,j]),"\n")
    cat("\tVariance:\t",var(samples[i,j]),"\n")
    cat("\t> 1.96:\t\t",mean(samples[i,j]>1.96)*100,"%\n")
  }
  if (p>1) cat("\nCorrelation:",cor(samples[i,1],samples[i,2]),"\n")
  invisible(samples[i,])
}