BCPA/PlotBCPA

Plotting function for BCPA output.

PlotBCPA <- function (t, x, ws, pp, threshold=10, col.cp= rgb(1,0.5,0,.5)) 
{
  plot(t, x, type = "n", main = "", ylim = c(min(x), max(x) * 1.3))
  x.breaks <- ws$Break
  x.model <- ws$Model
  x.breaks <- x.breaks[x.model > 0]
  x.model <- x.model[x.model > 0]
  mids <- hist(x.breaks, breaks = t, plot = F)$mid
  freq <- hist(x.breaks, breaks = t, plot = F)$count
  if(threshold > max(freq))
  {
      print("Threshold too high!")
      return(NULL)
  }
  
  goodbreaks <- mids[freq > threshold]
  freq <- freq[freq > threshold]
  break.cols <- heat.colors(max(freq))
  break.cols <- break.cols[length(break.cols):1]
  abline(v = goodbreaks, lwd = freq, col = col.cp)
  lines(t, x, col = "darkgrey")
  lines(t, pp$mu.hat, lwd = 2)
  lines(t, pp$mu.hat + pp$s.hat, col = 2, lwd = 1.5)
  lines(t, pp$mu.hat - pp$s.hat, col = 2, lwd = 1.5)
  rho.hat <- pp$rho.hat
  rho.int <- round((rho.hat - min(rho.hat, na.rm = 1)) * 1000) + 
    1
  rho.col <- topo.colors(max(rho.int, na.rm=TRUE))[rho.int]
  rho.cex <- (rho.hat - min(rho.hat, na.rm=TRUE))/(diff(range(rho.hat))) * 
    1.5 + 0.5
  points(t, x, col = rho.col, pch = 19)
  rho.index <- quantile(0:max(rho.int, na.rm=TRUE))
  rho.index[1] <- 1
  legend.cols <- topo.colors(max(rho.int, na.rm=TRUE))[round(rho.index)]
  legend.rhos <- round(quantile((min(rho.hat, na.rm=TRUE) * 1000):(max(rho.hat, na.rm=TRUE) * 
                                                         1000))/1000, 2)
  legend("bottomright", bg = "white", legend = c(expression(hat(rho)), 
                                                 legend.rhos), pch = c(0, rep(19, 5)), ncol = 2, col = c(0, 
                                                                                                         legend.cols), xjust = 0.5, yjust = 0.5, cex = 1.2)
  legend("topright", bg = "white", legend = c(expression(hat(mu)), 
                                              expression(hat(mu) %+-% hat(sigma))), lty = 1, lwd = 2:1, 
         col = 1:2, xjust = 0.5, yjust = 0.5, cex = 1.2)
}