pkg/R/z_plot.R

   1 #' @include z_plotHelper.R
   2
   3 #' @title Plot forecasts/observations
   4 #'
   5 #' @description Plot the measures at one station versus all experts forecasts.
   6 #'
   7 #' @param r Output of \code{\link{runAlgorithm}}.
   8 #' @param station Name or index of the station to consider. Default: the first one
   9 #' @param interval Time interval for the plot. Default: all time range.
  10 #' @param experts Subset of experts for the plot. Default: all experts.
  11 #' @param ... Additional arguments to be passed to graphics::plot method.
  12 #'
  13 #' @export
  14 plotCurves = function(r, station=1, interval=1:(nrow(r$data)/length(r$stations)), experts=r$experts, cols=rainbow(length(experts)), ...)
  15 {
  16     if (is.character(station))
  17         station = match(station, r$stations)
  18     if (is.numeric(experts))
  19         experts = r$experts[experts]
  20
  21     XY = subset(r$data[interval,], subset = (Station == station), select = c(experts,"Measure"))
  22     indices = getNoNAindices(XY)
  23     XY = XY[indices,]
  24     X = as.matrix(XY[,names(XY) %in% experts])
  25     Y = XY[,"Measure"]
  26
  27     yRange = range(XY)
  28     par(mar=c(5,4.5,1,1), cex=1.5)
  29     for (i in 1:length(experts))
  30     {
  31         plot(X[,i],ylim=yRange,type="l",lty="dotted",col=cols[i],xlab="",ylab="",xaxt="n",yaxt="n", lwd=2, ...)
  32         par(new=TRUE)
  33     }
  34     plot(Y, type="l", ylim=yRange, xlab="", ylab="", lwd=2, cex.axis=1.5, ...)
  35     title(xlab="Time",ylab="Forecasts / Measures", cex.lab=1.6)
  36     legend("topright", title="Historical PM10",lwd=c(2,1),lty=c("solid","dotted"),horiz=TRUE,legend=c("Measures","Forecasts"))
  37 }
  38
  39 #' @title Plot error
  40 #'
  41 #' @description Plot the absolute error over time at one station.
  42 #'
  43 #' @param r Output of \code{\link{runAlgorithm}}.
  44 #' @param station Name or index of the station to consider. Default: the first one
  45 #' @param start First index to consider (too much variability in early errors)
  46 #' @param noNA TRUE to show only errors associated with full lines (old behavior)
  47 #' @param ... Additional arguments to be passed to graphics::plot method.
  48 #'
  49 #' @export
  50 plotError = function(r, station=1, start=1, noNA=TRUE, ...)
  51 {
  52     if (is.character(station))
  53         station = match(station, r$stations)
  54
  55     XY = subset(r$data, subset = (Station == station), select = c(r$experts,"Measure","Prediction"))
  56     Y = XY[,"Measure"]
  57     hatY = XY[,"Prediction"]
  58     indices = !is.na(Y) & !is.na(hatY)
  59     if (noNA)
  60     {
  61         X = XY[,names(XY) %in% r$experts]
  62         indices = indices & getNoNAindices(X)
  63     }
  64     Y = Y[indices]
  65     hatY = hatY[indices]
  66
  67     error = abs(Y - hatY)
  68     par(mar=c(5,4.5,1,1), cex=1.5)
  69     plot(error, type="l", xaxt="n", xlab="Time",ylab="L1 error", cex.lab=1.6, cex.axis=1.5, ...)
  70     axis(side=1, at=(seq(from=start,to=length(Y),by=30) - start), labels=seq(from=start,to=length(Y),by=30), cex.axis=1.5)
  71 }
  72
  73 #' @title Plot regret
  74 #'
  75 #' @description Plot the regret over time at one station.
  76 #'
  77 #' @param r Output of \code{\link{runAlgorithm}}.
  78 #' @param vs Linear weights to compare with. Can be obtained by the \code{getBestXXX} methods, or by any other mean.
  79 #' @param station Name or index of the station to consider. Default: the first one
  80 #' @param start First index to consider (too much variability in early errors)
  81 #' @param ... Additional arguments to be passed to graphics::plot method.
  82 #'
  83 #' @export
  84 plotRegret = function(r, vs, station=1, start=1, ...)
  85 {
  86     if (is.character(station))
  87         station = match(station, r$stations)
  88
  89     XY = subset(r$data, subset = (Station == station), select = c(r$experts,"Measure","Prediction"))
  90     X = XY[,names(XY) %in% r$experts]
  91     Y = XY[,"Measure"]
  92     hatY = XY[,"Prediction"]
  93
  94     indices = !is.na(Y) & !is.na(hatY) & getNoNAindices(X)
  95     X = as.matrix(X[indices,])
  96     Y = Y[indices]
  97     hatY = hatY[indices]
  98
  99     error2 = abs(Y - hatY)^2
 100     vsError2 = abs(Y - X %*% vs)^2
 101     cumErr2 = cumsum(error2) / seq_along(error2)
 102     cumVsErr2 = cumsum(vsError2) / seq_along(vsError2)
 103     regret = cumErr2 - cumVsErr2
 104
 105     par(mar=c(5,4.5,1,1), cex=1.5)
 106     plot(regret, type="l", xaxt="n", xlab="Time", ylab="Regret", cex.lab=1.6, cex.axis=1.5, ...)
 107     abline(a=0., b=0., col=2)
 108     axis(side=1, at=(seq(from=start,to=length(Y),by=30) - start), labels=seq(from=start,to=length(Y),by=30), cex.axis=1.5)
 109 }
 110
 111 #' @title Plot predicted/expected cloud
 112 #'
 113 #' @description Plot the cloud of forecasts/observations + statistical indicators.
 114 #'
 115 #' @param r Output of \code{\link{runAlgorithm}}.
 116 #' @param thresh Threshold to consider for alerts (usually 30 or 50)
 117 #' @param hintThresh thresholds to draw on the plot to help visualization. Often \code{c(30,50,80)}
 118 #' @param station Name or index of the station to consider. Default: the first one
 119 #' @param noNA TRUE to show only errors associated with full lines (old behavior)
 120 #' @param ... Additional arguments to be passed to graphics::plot method.
 121 #'
 122 #' @export
 123 plotCloud = function(r, thresh=30, hintThresh=c(30,50,80), station=1, noNA=TRUE, ...)
 124 {
 125     if (is.character(station))
 126         station = match(station, r$stations)
 127
 128     XY = subset(r$data, subset = (Station == station), select = c(r$experts,"Measure","Prediction"))
 129     Y = XY[,"Measure"]
 130     hatY = XY[,"Prediction"]
 131     indices = !is.na(Y) & !is.na(hatY)
 132     if (noNA)
 133     {
 134         X = XY[,names(XY) %in% r$experts]
 135         indices = indices & getNoNAindices(X)
 136     }
 137     Y = Y[indices]
 138     hatY = hatY[indices]
 139
 140     indics = getIndicators(r, thresh, station, noNA)
 141
 142     par(mar=c(5,5,3,2), cex=1.5)
 143     plot(Y, hatY, xlab="Measured PM10", ylab="Predicted PM10",
 144         cex.lab=1.6, cex.axis=1.5, xlim=c(0,120), ylim=c(0,120), ...)
 145     abline(0,1,h=hintThresh,v=hintThresh,col=2,lwd=2)
 146 #   legend("topleft",legend=c(paste("EV   ",indics$EV),paste("RMSE ",indics$RMSE)),cex=1.2)
 147     legend("topleft",legend=paste("RMSE ",indics$RMSE))
 148     legend("bottomright",legend=c(paste("TS ",indics$TS)))
 149 }