R/diagplots.R

In BNPTSclust: A Bayesian Nonparametric Algorithm for Time Series Clustering

diagplots <-
function(L){
    
    # Function that produces the diagnostic plots to assess the convergence
    # of the Markov Chains generated by either of the functions:
    # "tseriesca", "tseriescm" or "tseriescq". 
    #
    # IN:
    # 
    # L       <- output list from the functions: "tseriesca", "tseriescm"
    #            or "tseriescq".
    #
    # OUT:
    #
    # Trace plots, histograms and ergodic mean plots of the posterior
    # distribution sample from the model parameters. 
    
    sig2epssample <- L$sig2epssample
    sig2alphasample <- L$sig2alphasample
    sig2betasample <- L$sig2betasample
    sig2thesample <- L$sig2thesample
    rhosample <- L$rhosample
    asample <- L$asample
    bsample <- L$bsample
    msample <- L$msample
    
    n <- length(msample)
    
    # Trace plots
    
    par(mfrow = c(2,2))
    plot(1:n,sig2epssample[,1],type = "l",main = "Trace plot of sig2eps",xlab = "Iteration number",ylab = "Simulated value")
    plot(1:n,sig2thesample,type = "l",main = "Trace plot of sig2theta",xlab = "Iteration number",ylab = "Simulated value")
    if(is.null(sig2betasample) == TRUE){
      plot(1:n,sig2alphasample[,1],type = "l",main = "Trace plot of sig2alpha",xlab = "Iteration number",ylab = "Simulated value")
    }else if(is.null(sig2alphasample) == TRUE){
      plot(1:n,sig2betasample[,1],type = "l",main = "Trace plot of sig2beta",xlab = "Iteration number",ylab = "Simulated value")
    }else{
      plot(1:n,sig2alphasample[,1],type = "l",main = "Trace plot of sig2alpha",xlab = "Iteration number",ylab = "Simulated value")
      plot(1:n,sig2betasample[,1],type = "l",main = "Trace plot of sig2beta",xlab = "Iteration number",ylab = "Simulated value")
    }
    
    par(mfrow = c(2,2))
    plot(1:n,rhosample,type = "l",main = "Trace plot of rho",xlab = "Iteration number",ylab = "Simulated value")
    plot(1:n,asample,type = "l",main = "Trace plot of a",xlab = "Iteration number",ylab = "Simulated value")
    plot(1:n,bsample,type = "l",main = "Trace plot of b",xlab = "Iteration number",ylab = "Simulated value")
    plot(1:n,msample,type = "l",main = "m",xlab = "Iteration number",ylab = "Number of groups")
    
    # Histograms
    
    par(mfrow = c(2,2))
    hist(sig2epssample[,1],main = "Hist. of sig2eps",xlab = "Simulated values")
    hist(sig2thesample,main = "Hist. of sig2the",xlab = "Simulated values")
    if(is.null(sig2betasample) == TRUE){
      hist(sig2alphasample[,1],main = "Hist. of sig2alpha",xlab = "Simulated values")
    }else if(is.null(sig2alphasample) == TRUE){
      hist(sig2betasample[,1],main = "Hist. of sig2beta",xlab = "Simulated values")
    }else{
      hist(sig2alphasample[,1],main = "Hist. of sig2alpha",xlab = "Simulated values")
      hist(sig2betasample[,1],main = "Hist. of sig2beta",xlab = "Simulated values")
    }
    
    par(mfrow = c(2,2))
    hist(rhosample,main = "Hist. of rho",xlab = "Simulated values")
    hist(asample,main = "Hist. of a",xlab = "Simulated values")
    hist(bsample,main = "Hist. of b",xlab = "Simulated values")
    hist(msample,main = "Hist. of m",xlab = "Number of groups")
    
    # Ergodic means
    
    par(mfrow = c(2,2))
    plot(1:n,cumsum(sig2epssample[,1])/1:n,type = "l",main = "Ergodic mean of sig2eps",xlab = "Iteration number",ylab = "")
    plot(1:n,cumsum(sig2thesample)/1:n,type = "l",main = "Ergodic mean of sig2the",xlab = "Iteration number",ylab = "")
    if(is.null(sig2betasample) == TRUE){
      plot(1:n,cumsum(sig2alphasample[,1])/1:n,type = "l",main = "Ergodic mean of sig2alpha",xlab = "Iteration number",ylab = "")
    }else if(is.null(sig2alphasample) == TRUE){
      plot(1:n,cumsum(sig2betasample[,1])/1:n,type = "l",main = "Ergodic mean of sig2beta",xlab = "Iteration number",ylab = "")
    }else{
      plot(1:n,cumsum(sig2alphasample[,1])/1:n,type = "l",main = "Ergodic mean of sig2alpha",xlab = "Iteration number",ylab = "")
      plot(1:n,cumsum(sig2betasample[,1])/1:n,type = "l",main = "Ergodic mean of sig2beta",xlab = "Iteration number",ylab = "")
    }
    
    par(mfrow = c(2,2))
    plot(1:n,cumsum(rhosample)/1:n,type = "l",main = "Ergodic mean of rho",xlab = "Iteration number",ylab = "")
    plot(1:n,cumsum(asample)/1:n,type = "l",main = "Ergodic mean of a",xlab = "Iteration number",ylab = "")
    plot(1:n,cumsum(bsample)/1:n,type = "l",main = "Ergodic mean of b",xlab = "Iteration number",ylab = "")
    plot(1:n,cumsum(msample)/1:n,type = "l",main = "Ergodic mean of m",xlab = "Iteration number",ylab = "")
    
  }