#Set up the plotting region
def.par <- par(no.readonly = TRUE)
nf <- layout(matrix(c(1,2),2,1,byrow=TRUE), heights=c(2,1),widths=c(1))
layout.show(nf)

if(n.modes == 1) p.max <- max(dnorm(x,true.mean,true.sd))
if(n.modes == 2) p.max <- max(p.1*dnorm(x,true.mean.1,true.sd.1)+p.2*dnorm(x,true.mean.2,true.sd.2))

par(ask=T,mar=c(6,2,2,2))
for(k in 1:K){
      theta.proposed <- rnorm(1,theta.old,proposal.sd)
      
      if(n.modes == 1){
                 plot(x,dnorm(x,true.mean,true.sd),type='l',yaxt='n',bty='n',xaxs='r',xpd=F,xlab="theta",ylab="",main="Posterior Mean of theta")
                 points(theta.old,dnorm(theta.old,true.mean,true.sd),col="red",pch=19,cex=1.5)
                 points(theta.proposed,dnorm(theta.proposed,true.mean,true.sd),col="green",pch=19,cex=1.5) 
                 arrows(theta.old,dnorm(theta.old,true.mean,true.sd),theta.proposed,dnorm(theta.proposed,true.mean,true.sd),lwd=2,length=.1)
                 alpha <- calculate.alpha.unimodal(theta.proposed,theta.old,true.mean,true.sd)
                 }
      if(n.modes == 2){
                 plot(x,p.1*dnorm(x,true.mean.1,true.sd.1)+p.2*dnorm(x,true.mean.2,true.sd.2),type='l',yaxt='n',bty='n',xaxs='r',xpd=F,xlab="theta",ylab="",main="Posterior Mean of theta")
                 points(theta.old,p.1*dnorm(theta.old,true.mean.1,true.sd.1)+p.2*dnorm(theta.old,true.mean.2,true.sd.2),col="red",pch=19,cex=1.5)
                 points(theta.proposed,p.1*dnorm(theta.proposed,true.mean.1,true.sd.1)+p.2*dnorm(theta.proposed,true.mean.2,true.sd.2),col="green",pch=19,cex=1.5)  
                 arrows(theta.old,p.1*dnorm(theta.old,true.mean.1,true.sd.1)+p.2*dnorm(theta.old,true.mean.2,true.sd.2),theta.proposed,p.1*dnorm(theta.proposed,true.mean.1,true.sd.1)+p.2*dnorm(theta.proposed,true.mean.2,true.sd.2),lwd=2,length=.1)
                 alpha <- calculate.alpha.bimodal(theta.proposed,theta.old,p.1,p.2,true.mean.1,true.sd.1,true.mean.2,true.sd.2)
                 }
      text(10,p.max,labels=paste("alpha = ",round(alpha,2)),pos=1)
      if(runif(1,0,1) < alpha){
                      accept[k] <- 1
                      theta.old <- theta.proposed
                      text(10,p.max-.02,labels="Accept",col="green",cex=2,pos=1)
                      }else{ 
                            accept[k] <- 0
                            text(10,p.max-.02,labels="Reject",col="red",cex=2,pos=1)
                      }
      theta.samples[k] <- theta.old
      if(k > k.hist){
       hist(theta.samples[1:k],prob=T,add=T,border='lightgrey')
      }
                
      plot(1:k,theta.samples[1:k],type='l',xlim=c(1,K),ylim=c(0,10),xlab="",ylab="0",sub=(paste("Acceptance Rate = ",round(mean(accept[1:k]),2))),main="Trace plot of theta samples")
      
    
}

par(def.par)   #reset to default