R/plotmcmc.R

In atmcmc: Automatically Tuned Markov Chain Monte Carlo

#' plots output from `atmcmc'
#' 
#' Produces traceplots and histograms of output from an `atmcmc' run 
#' 
#' Histograms and traceplots are drawn for each coordinate separately. For the sampling phase, traceplots show only one replicative chain (chain that takes the last value of the 2nd adaption phase as the starting value of the sampling phase).
#' @param output output from an `atmcmc' run 
#' @param name name of the project
#' @param multimodal whether to assume target density function is strongly multimodal. Argument takes T or F
#' @param plottype whether to produce traceplots or histograms, or both. Takes argument `trace', `hist', or `all'
#' @param format file format. R graphics, pdf, or png. Takes argument `default', `pdf', or `png'
#' @param m every `m'th iteration is plotted. Has to be a factor of `batchwidth'/2 from the `atmcmc' run
#' @param phase.start first phase to be plotted in traceplots. Takes argument 1, 2, 3, 4, or 5. 1 = 1st adaption phase, 2 = transient phase, 3 = 2nd adaption phase, 4 = 1st half of sampling phase, and 5 = 2nd half of sampling phase
#' @param phase.end last phase to be plotted in traceplots. Takes argument 1, 2, 3, 4, or 5. 1 = 1st adaption phase, 2 = transient phase, 3 = 2nd adaption phase, 4 = 1st half of sampling phase, and 5 = 2nd half of sampling phase
#' @param nrow.trace number of rows for traceplots per page 
#' @param ncol.trace number of columns for traceplots per page 
#' @param nrow.hist number of rows for histograms per page 
#' @param ncol.hist number of columns for histograms per page
#' @return Plots including
#'  \item{`name'_histogram}{histograms of the sample obtained from the `atmcmc' run}
#'  \item{`name'_traceplot}{traceplots from the `atmcmc' run. Every `m'th iteration is plotted}
#' @examples ## see examples in `atmcmc'
#' @export
plotmcmc <- function ( output, name = "MCMC", multimodal=F, plottype='trace', format='default', m=10, phase.start=1, phase.end=5, nrow.trace=3, ncol.trace=1, nrow.hist=3, ncol.hist=3) {
  
  batchwidth=output$batchwidth
  if ((multimodal!=F)&(multimodal!=T)) {stop("Argument multimodal should be either T or F", sep="")}
  if ((plottype!='hist')&(plottype!='trace')&(plottype!='all')) {stop("Argument format should be one of 'trace', 'hist' or 'all'", sep="")}
  if ((format!='default')&(format!='pdf')&(format!='png')) {stop("Argument format should be one of 'default', 'pdf' or 'png'", sep="")}
  if (length(m)!=1){stop("Argument m should be scalar", sep="")}
  if (((batchwidth/2) %% m)!=0){stop("m=",m, " has to be a factor of batchwidth/2=",batchwidth/2, sep="")}
  if (length(phase.start)!=1){stop("Argument phase.start should be scalar", sep="")}
  if (length(phase.end)!=1){stop("Argument phase.end should be scalar", sep="")}
  if ((phase.start!=1)&(phase.start!=2)&(phase.start!=3)&(phase.start!=4)&(phase.start!=5)){stop("Argument phase.start should be one of 1, 2, 3, 4 or 5", sep="")}
  if ((phase.end!=1)&(phase.end!=2)&(phase.end!=3)&(phase.end!=4)&(phase.end!=5)){stop("Argument phase.start should be one of 1, 2, 3, 4 or 5", sep="")}
  if (phase.end<phase.start){stop("Argument phase.end should be greater than or equal to phase.start", sep="")}
  if (length(nrow.hist)!=1){stop("Argument nrow.hist should be scalar", sep="")}
  if (length(ncol.hist)!=1){stop("Argument ncol.hist should be scalar", sep="")}
  if (length(nrow.trace)!=1){stop("Argument nrow.trace should be scalar", sep="")}
  if (length(ncol.trace)!=1){stop("Argument ncol.trace should be scalar", sep="")}
  
  retrieve <- function(argString, argInt) {return( eval(as.name(paste(argString, argInt, sep=""))))}
  
  Xveclistbase = output$Xveclistbase
  dim = output$dim
  Xveclistdim1=output$Xveclistdim1
  if (dim>1){
    for (j in 1:dim){
      expression =paste("Xveclistdim",j,"=output$Xveclistdim",j, sep = "")
      eval(parse(text = expression))
    }
  }
  
  if(multimodal==F){
    adpstop1 = output$sumchain[1]
    transientstop = output$sumchain[2]
    adpstop2 = output$sumchain[3]
    discard = output$sumchain[4]
    chainstop = output$sumchain[5]
    Xveclist=Xveclistdim1[(1:(chainstop-adpstop2)),1]
    if (dim>1){
      for (j in 2:dim){
        expression =paste("Xveclist=cbind(Xveclist,Xveclistdim",j, "[(1:(chainstop-adpstop2)),1])", sep = "")
        eval(parse(text = expression))
      }
    }
    
    #histogram of the sample generated by MCMC run 
    #combines all points from all replicative chains in the sampling phase
    #is drawn for each coordinate seperately 
    if ((plottype=='hist')|(plottype=='all')){
      numtemp=0
      xlab <- bquote(bold(paste("n = ", .((chainstop-discard)*(dim(Xveclistdim1)[2]))))) 
      if (format=='pdf'){
        pdf(paste(name, "_histogram.pdf", sep=""))
        par(mfrow=c(nrow.hist,ncol.hist))
      }
      for (i in 1:dim){
        title <- bquote(bold(paste("Density Plot Coordinate " , . (i)))) 
        if (format=='default'){
          if ((numtemp %% (nrow.hist*ncol.hist)) == 0){
            dev.new(noRStudioGD = TRUE)
            par(mfrow=c(nrow.hist,ncol.hist)) 
          }
        }
        if (format=='png'){
          if ((numtemp %% (nrow.hist*ncol.hist)) == 0){
            png(paste(name, (numtemp/(nrow.hist*ncol.hist)+1), "_histogram.png", sep=""))
            par(mfrow=c(nrow.hist,ncol.hist))
          }
        }
        hist(retrieve("Xveclistdim", i)[((discard-adpstop2+1):(chainstop-adpstop2)),],freq=FALSE,breaks=50,xlab=xlab,ylab=".",main=title)
        numtemp=numtemp+1
        if (format=='png'){
          if (((numtemp>0)&((numtemp %% (nrow.hist*ncol.hist)) == 0))|(i==dim)){
            dev.off()
          }
        }
      }
      if (format =='pdf'){
        dev.off()
      }
    }
    
    #Traceplots of each coordinate (every `m'th)
    if ((plottype=='trace')|(plottype=='all')){
      if (phase.end==1){endpoint=adpstop1}
      if (phase.end==2){endpoint=transientstop}
      if (phase.end==3){endpoint=adpstop2}
      if (phase.end==4){endpoint=discard}
      if (phase.end==5){endpoint=chainstop}
      numtemp=0
      xveclist<-matrix(0,chainstop,1)
      if (format=='pdf'){
        pdf(paste(name, "_traceplot.pdf", sep=""))
        par(mfrow=c(nrow.trace,ncol.trace))         
      }
      for (i in 1:dim){
        xveclist[1:adpstop2]<-Xveclistbase[,i]
        if (dim==1){xveclist[(adpstop2+1):chainstop]<-Xveclist} else{xveclist[(adpstop2+1):chainstop]<-Xveclist[,i]}
        Xveclisttmp <- xveclist[seq(m,length(xveclist),by=m)]
        title <- bquote(bold(paste("Trace Plot Coordinate " , . (i)))) 
        if (format=='default'){
          if ((numtemp %% (nrow.trace*ncol.trace)) == 0){
            dev.new(noRStudioGD = TRUE)
            par(mfrow=c(nrow.trace,ncol.trace)) 
          }
        }
        if (format=='png'){
          if ((numtemp %% (nrow.trace*ncol.trace)) == 0){
            png(paste(name, (numtemp/(nrow.trace*ncol.trace)+1), "_traceplot.png", sep=""))
            par(mfrow=c(nrow.trace,ncol.trace))
          }
        }
        if (phase.start==1){
          plot((seq(m, adpstop1, by = m)), Xveclisttmp[(1:(adpstop1/m))], type="l", main=title, xlab = "iteration", ylab = " ", col="purple", xlim=c(1,(endpoint)),ylim=c(min(Xveclisttmp[1:(endpoint/m)]),max(Xveclisttmp[1:(endpoint/m)])))
          if (phase.end>1){lines((seq((adpstop1), transientstop, by = m)),Xveclisttmp[(adpstop1/m):(transientstop/m)], type="l", col="orange")}
          if (phase.end>2){lines((seq((transientstop), adpstop2, by = m)),Xveclisttmp[(transientstop/m):(adpstop2/m)], type="l", col="red")}
          if (phase.end>3){lines((seq((adpstop2), discard, by = m)),Xveclisttmp[(adpstop2/m):(discard/m)], type="l", col="blue")}
          if (phase.end>4){lines((seq((discard), chainstop, by = m)),Xveclisttmp[(discard/m):(chainstop/m)], type="l", col="green")}
        }
        if (phase.start==2){
          plot((seq((adpstop1+m), transientstop, by = m)), Xveclisttmp[((adpstop1/m+1):(transientstop/m))], type="l", main=title, xlab = "iteration", ylab = " ", col="orange", xlim=c((adpstop1+1),(endpoint)),ylim=c(min(Xveclisttmp[(adpstop1/m+1):(endpoint/m)]),max(Xveclisttmp[(adpstop1/m+1):(endpoint/m)])))
          if (phase.end>2){lines((seq((transientstop), adpstop2, by = m)),Xveclisttmp[(transientstop/m):(adpstop2/m)], type="l", col="red")}
          if (phase.end>3){lines((seq((adpstop2), discard, by = m)),Xveclisttmp[(adpstop2/m):(discard/m)], type="l", col="blue")}
          if (phase.end>4){lines((seq((discard), chainstop, by = m)),Xveclisttmp[(discard/m):(chainstop/m)], type="l", col="green")}
        }
        if (phase.start==3){
          plot((seq((transientstop+m), adpstop2, by = m)), Xveclisttmp[((transientstop/m+1):(adpstop2/m))], type="l", main=title, xlab = "iteration", ylab = " ", col="red", xlim=c((transientstop+1),(endpoint)),ylim=c(min(Xveclisttmp[(transientstop/m+1):(endpoint/m)]),max(Xveclisttmp[(transientstop/m+1):(endpoint/m)])))
          if (phase.end>3){lines((seq((adpstop2), discard, by = m)),Xveclisttmp[(adpstop2/m):(discard/m)], type="l", col="blue")}
          if (phase.end>4){lines((seq((discard), chainstop, by = m)),Xveclisttmp[(discard/m):(chainstop/m)], type="l", col="green")}
        }
        if (phase.start==4){
          plot((seq((adpstop2+m), discard, by = m)), Xveclisttmp[((adpstop2/m+1):(discard/m))], type="l", main=title, xlab = "iteration", ylab = " ", col="blue", xlim=c((adpstop2+1),(endpoint)),ylim=c(min(Xveclisttmp[(adpstop2/m+1):(endpoint/m)]),max(Xveclisttmp[(adpstop2/m+1):(endpoint/m)])))
          if (phase.end>4){lines((seq((discard), chainstop, by = m)),Xveclisttmp[(discard/m):(chainstop/m)], type="l", col="green")}
        }
        if (phase.start==5){
          plot((seq((discard+m), chainstop, by = m)), Xveclisttmp[((discard/m+1):(chainstop/m))], type="l", main=title, xlab = "iteration", ylab = " ", col="green", xlim=c((discard+1),(endpoint)),ylim=c(min(Xveclisttmp[(discard/m+1):(endpoint/m)]),max(Xveclisttmp[(discard/m+1):(endpoint/m)])))
        }
        numtemp=numtemp+1
        if (format=='png'){
          if (((numtemp>0)&((numtemp %% (nrow.trace*ncol.trace)) == 0))|(i==dim)){
            dev.off()
          }
        }
      }
      if (format =='pdf'){
        dev.off()
      }
    }
    
  }
  
  else{
    adpstop1 = output$sumchain[1,]
    transientstop = output$sumchain[2,]
    adpstop2 = output$sumchain[3,]
    adpstop2_max = output$sumchain[4,1]
    discard = output$sumchain[5,1]
    chainstop = output$sumchain[6,1]
    nummodes = output$nummodes
    Xveclist=Xveclistdim1[(1:(chainstop-adpstop2_max)),1]
    if (dim>1){
      for (j in 2:dim){
        expression =paste("Xveclist=cbind(Xveclist, Xveclistdim",j, "[(1:(chainstop-adpstop2_max)),1])", sep = "")
        eval(parse(text = expression))
      }
    }
    
    #histogram of the sample generated by MCMC run 
    #combines all points from all replicative chains in the sampling phase
    #is drawn for each coordinate seperately 
    if ((plottype=='hist')|(plottype=='all')){
      numtemp=0
      xlab <- bquote(bold(paste("n = ", .((chainstop-discard)*(dim(Xveclistdim1)[2]))))) 
      if (format=='pdf'){
        pdf(paste(name, "_histogram.pdf", sep=""))
        par(mfrow=c(nrow.hist,ncol.hist))
      }
      for (i in 1:dim){
        title <- bquote(bold(paste("Density Plot Coordinate " , . (i)))) 
        if (format=='default'){
          if ((numtemp %% (nrow.hist*ncol.hist)) == 0){
            dev.new(noRStudioGD = TRUE)
            par(mfrow=c(nrow.hist,ncol.hist)) 
          }
        }
        if (format=='png'){
          if ((numtemp %% (nrow.hist*ncol.hist)) == 0){
            png(paste(name, (numtemp/(nrow.hist*ncol.hist)+1), "_histogram.png", sep=""))
            par(mfrow=c(nrow.hist,ncol.hist))
          }
        }
        hist(retrieve("Xveclistdim", i)[((discard-adpstop2_max+1):(chainstop-adpstop2_max)),],freq=FALSE,breaks=50,xlab=xlab,ylab=".",main=title)
        numtemp=numtemp+1
        if (format=='png'){
          if (((numtemp>0)&((numtemp %% (nrow.hist*ncol.hist)) == 0))|(i==dim)){
            dev.off()
          }
        }
      }
      if (format =='pdf'){
        dev.off()
      }
    }
    
    
    #Traceplots of each coordinate (every `m'th)
    if ((plottype=='trace')|(plottype=='all')){
      if (phase.end==1){endpoint=adpstop1}
      if (phase.end==2){endpoint=transientstop}
      if (phase.end==3){endpoint=adpstop2}
      if (phase.end==4){endpoint=discard}
      if (phase.end==5){endpoint=chainstop}
      if (phase.start==2){startpoint=adpstop1}
      if (phase.start==3){startpoint=transientstop}
      ymax = ymin = ymax.simple.partial = ymin.simple.partial = matrix(0,nummodes,dim)
      for (j in 1:nummodes){
        xveclist<-matrix(0,chainstop,1)
        for (i in 1:dim){
          if (dim==1){
            xveclist[1:adpstop2[j]]<-Xveclistbase[(1:adpstop2[j]),j]
            xveclist[(adpstop2_max+1):chainstop]<-Xveclist
          } 
          else{
            xveclist[1:adpstop2[j]]<-Xveclistbase[(1:adpstop2[j]),i,j]
            xveclist[(adpstop2_max+1):chainstop]<-Xveclist[,i]
          }
          
          Xveclisttmp <- xveclist[seq(m,length(xveclist),by=m)]
          if (phase.start==1){
            if (phase.end<4){
              ymax[j,i] = max(Xveclisttmp[1:(adpstop2[j]/m)])
              ymin[j,i] = min(Xveclisttmp[1:(adpstop2[j]/m)])
            }
            else{
              ymax[j,i] = max(max(Xveclisttmp[1:(adpstop2[j]/m)]),max(Xveclisttmp[(adpstop2_max/m+1):(endpoint/m)]))
              ymin[j,i] = min(min(Xveclisttmp[1:(adpstop2[j]/m)]),min(Xveclisttmp[(adpstop2_max/m+1):(endpoint/m)]))
            }
          }
          if ((phase.start==2)|(phase.start==3)){
            if (phase.end<4){
              ymax[j,i] = max(Xveclisttmp[(startpoint[j]/m+1):(adpstop2[j]/m)])
              ymin[j,i] = min(Xveclisttmp[(startpoint[j]/m+1):(adpstop2[j]/m)])
            }
            else{
              ymax[j,i] = max(max(Xveclisttmp[(startpoint[j]/m+1):(adpstop2[j]/m)]),max(Xveclisttmp[(adpstop2_max/m+1):(endpoint/m)]))
              ymin[j,i] = min(min(Xveclisttmp[(startpoint[j]/m+1):(adpstop2[j]/m)]),min(Xveclisttmp[(adpstop2_max/m+1):(endpoint/m)]))
            }
          }
          if (phase.start>3){
            ymax[j,i] = max(Xveclisttmp[(adpstop2_max/m+1):(endpoint/m)])
            ymin[j,i] = min(Xveclisttmp[(adpstop2_max/m+1):(endpoint/m)])
          }
        }
      }
      
      numtemp=0
      if (format=='pdf'){
        pdf(paste(name, "_traceplot.pdf", sep=""))
        par(mfrow=c(nrow.trace,ncol.trace))         
      }
      if (phase.start>3){nummodes=1}
      for (i in 1:dim){
        for (j in 1:nummodes){
          xveclist<-matrix(0,chainstop,1)
          if (dim==1){
            xveclist[1:adpstop2[j]]<-Xveclistbase[(1:adpstop2[j]),j]
            xveclist[(adpstop2_max+1):chainstop]<-Xveclist
          } 
          else{
            xveclist[1:adpstop2[j]]<-Xveclistbase[(1:adpstop2[j]),i,j]
            xveclist[(adpstop2_max+1):chainstop]<-Xveclist[,i]
          }
          Xveclisttmp <- xveclist[seq(m,length(xveclist),by=m)]
          if (j==1){X0.sampling=Xveclisttmp[adpstop2[j]/m]}
          if (phase.start<4){title <- bquote(bold(paste("Trace Plot coordinate " , . (i), " start with mode ", . (j))))}
          else {title <- bquote(bold(paste("Trace Plot coordinate " , . (i))))}
          if (format=='default'){
            if ((numtemp %% (nrow.trace*ncol.trace)) == 0){
              dev.new(noRStudioGD = TRUE)
              par(mfrow=c(nrow.trace,ncol.trace)) 
            }
          }
          if (format=='png'){
            if ((numtemp %% (nrow.trace*ncol.trace)) == 0){
              png(paste(name, (numtemp/(nrow.trace*ncol.trace)+1), "_traceplot.png", sep=""))
              par(mfrow=c(nrow.trace,ncol.trace))
            }
          }
          if (phase.start==1){
            plot((seq(m, adpstop1[j], by = m)), Xveclisttmp[(1:(adpstop1[j]/m))], type="l", main=title, xlab = "iteration", ylab = " ", col="purple", xlim=c(1,(max(endpoint))),ylim=c(min(ymin[,i]),max(ymax[,i])))
            if (phase.end>1){lines((seq((adpstop1[j]), transientstop[j], by = m)),Xveclisttmp[(adpstop1[j]/m):(transientstop[j]/m)], type="l", col="orange")}
            if (phase.end>2){lines((seq((transientstop[j]), adpstop2[j], by = m)),Xveclisttmp[(transientstop[j]/m):(adpstop2[j]/m)], type="l", col="red")}
            if (phase.end>3){lines((seq((adpstop2_max), discard, by = m)),c(X0.sampling, Xveclisttmp[(adpstop2_max/m+1):(discard/m)]), type="l", col="blue")}
            if (phase.end>4){lines((seq((discard), chainstop, by = m)),Xveclisttmp[(discard/m):(chainstop/m)], type="l", col="green")}
          }
          if (phase.start==2){
            plot((seq((adpstop1[j]+m), transientstop[j], by = m)), Xveclisttmp[((adpstop1[j]/m+1):(transientstop[j]/m))], type="l", main=title, xlab = "iteration", ylab = " ", col="orange", xlim=c((min(startpoint)),(max(endpoint))),ylim=c(min(ymin[,i]),max(ymax[,i])))
            if (phase.end>2){lines((seq((transientstop[j]), adpstop2[j], by = m)),Xveclisttmp[(transientstop[j]/m):(adpstop2[j]/m)], type="l", col="red")}
            if (phase.end>3){lines((seq((adpstop2_max), discard, by = m)),c(X0.sampling, Xveclisttmp[(adpstop2_max/m+1):(discard/m)]), type="l", col="blue")}
            if (phase.end>4){lines((seq((discard), chainstop, by = m)),Xveclisttmp[(discard/m):(chainstop/m)], type="l", col="green")}
          }
          if (phase.start==3){
            plot((seq((transientstop[j]+m), adpstop2[j], by = m)), Xveclisttmp[((transientstop[j]/m+1):(adpstop2[j]/m))], type="l", main=title, xlab = "iteration", ylab = " ", col="red", xlim=c((min(startpoint)),(max(endpoint))),ylim=c(min(ymin[,i]),max(ymax[,i])))
            if (phase.end>3){lines((seq((adpstop2_max), discard, by = m)),c(X0.sampling, Xveclisttmp[(adpstop2_max/m+1):(discard/m)]), type="l", col="blue")}
            if (phase.end>4){lines((seq((discard), chainstop, by = m)),Xveclisttmp[(discard/m):(chainstop/m)], type="l", col="green")}
          }
          if (phase.start==4){
            plot((seq((adpstop2_max+m), discard, by = m)), Xveclisttmp[((adpstop2_max/m+1):(discard/m))], type="l", main=title, xlab = "iteration", ylab = " ", col="blue", xlim=c(adpstop2_max,endpoint),ylim=c(min(ymin[,i]),max(ymax[,i])))
            if (phase.end>4){lines((seq(discard, chainstop, by = m)),Xveclisttmp[(discard/m):(chainstop/m)], type="l", col="green")}
          }
          if (phase.start==5){
            plot((seq((discard+m), chainstop, by = m)), Xveclisttmp[((discard/m+1):(chainstop/m))], type="l", main=title, xlab = "iteration", ylab = " ", col="green", xlim=c(discard,endpoint),ylim=c(min(ymin[,i]),max(ymax[,i])))
          }
          numtemp=numtemp+1
          if (format=='png'){
            if (((numtemp>0)&((numtemp %% (nrow.trace*ncol.trace)) == 0))|((i==dim)&(j==nummodes))){
              dev.off()
            }
          }
        }
      }
      if (format =='pdf'){
        dev.off()
      }
    }
    
  }
  
}