BCPA/PlotBCPA

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

Latest revision as of 05:23, 18 December 2012

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)
}
Personal tools