[aggexp.git] / pkg / R / z_plot.R

#' @include z_plotHelper.R

#' @title Plot forecasts/observations
#'
#' @description Plot the measures at one station versus all experts forecasts.
#'
#' @param r Output of \code{\link{runAlgorithm}}.
#' @param station Name or index of the station to consider. Default: the first one
#' @param interval Time interval for the plot. Default: all time range.
#' @param experts Subset of experts for the plot. Default: all experts.
#' @param ... Additional arguments to be passed to graphics::plot method.
#'
#' @export
plotCurves = function(r, station=1, interval=1:(nrow(r$data)/length(r$stations)), experts=r$experts, cols=rainbow(length(experts)), ...)
{
	if (is.character(station))
		station = match(station, r$stations)
	if (is.numeric(experts))
		experts = r$experts[experts]

	XY = subset(r$data[interval,], subset = (Station == station), select = c(experts,"Measure"))
	indices = getNoNAindices(XY)
	XY = XY[indices,]
	X = as.matrix(XY[,names(XY) %in% experts])
	Y = XY[,"Measure"]

	yRange = range(XY)
	par(mar=c(5,4.5,1,1), cex=1.5)
	for (i in 1:length(experts))
	{
		plot(X[,i],ylim=yRange,type="l",lty="dotted",col=cols[i],xlab="",ylab="",xaxt="n",yaxt="n", lwd=2, ...)
		par(new=TRUE)
	}
	plot(Y, type="l", ylim=yRange, xlab="", ylab="", lwd=2, cex.axis=1.5, ...)
	title(xlab="Time",ylab="Forecasts / Measures", cex.lab=1.6)
	legend("topright", title="Historical PM10",lwd=c(2,1),lty=c("solid","dotted"),horiz=TRUE,legend=c("Measures","Forecasts"))
}

#' @title Plot error
#'
#' @description Plot the absolute error over time at one station.
#'
#' @param r Output of \code{\link{runAlgorithm}}.
#' @param station Name or index of the station to consider. Default: the first one
#' @param start First index to consider (too much variability in early errors)
#' @param noNA TRUE to show only errors associated with full lines (old behavior)
#' @param ... Additional arguments to be passed to graphics::plot method.
#'
#' @export
plotError = function(r, station=1, start=1, noNA=TRUE, ...)
{
	if (is.character(station))
		station = match(station, r$stations)

	XY = subset(r$data, subset = (Station == station), select = c(r$experts,"Measure","Prediction"))
	Y = XY[,"Measure"]
	hatY = XY[,"Prediction"]
	indices = !is.na(Y) & !is.na(hatY)
	if (noNA)
	{
		X = XY[,names(XY) %in% r$experts]
		indices = indices & getNoNAindices(X)
	}
	Y = Y[indices]
	hatY = hatY[indices]

	error = abs(Y - hatY)
	par(mar=c(5,4.5,1,1), cex=1.5)
	plot(error, type="l", xaxt="n", xlab="Time",ylab="L1 error", cex.lab=1.6, cex.axis=1.5, ...)
	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)
}

#' @title Plot regret
#'
#' @description Plot the regret over time at one station.
#'
#' @param r Output of \code{\link{runAlgorithm}}.
#' @param vs Linear weights to compare with. Can be obtained by the \code{getBestXXX} methods, or by any other mean.
#' @param station Name or index of the station to consider. Default: the first one
#' @param start First index to consider (too much variability in early errors)
#' @param ... Additional arguments to be passed to graphics::plot method.
#'
#' @export
plotRegret = function(r, vs, station=1, start=1, ...)
{
	if (is.character(station))
		station = match(station, r$stations)

	XY = subset(r$data, subset = (Station == station), select = c(r$experts,"Measure","Prediction"))
	X = XY[,names(XY) %in% r$experts]
	Y = XY[,"Measure"]
	hatY = XY[,"Prediction"]

	indices = !is.na(Y) & !is.na(hatY) & getNoNAindices(X)
	X = as.matrix(X[indices,])
	Y = Y[indices]
	hatY = hatY[indices]

	error2 = abs(Y - hatY)^2
	vsError2 = abs(Y - X %*% vs)^2
	cumErr2 = cumsum(error2) / seq_along(error2)
	cumVsErr2 = cumsum(vsError2) / seq_along(vsError2)
	regret = cumErr2 - cumVsErr2

	par(mar=c(5,4.5,1,1), cex=1.5)
	plot(regret, type="l", xaxt="n", xlab="Time", ylab="Regret", cex.lab=1.6, cex.axis=1.5, ...)
	abline(a=0., b=0., col=2)
	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)
}

#' @title Plot predicted/expected cloud
#'
#' @description Plot the cloud of forecasts/observations + statistical indicators.
#'
#' @param r Output of \code{\link{runAlgorithm}}.
#' @param thresh Threshold to consider for alerts (usually 30 or 50)
#' @param hintThresh thresholds to draw on the plot to help visualization. Often \code{c(30,50,80)}
#' @param station Name or index of the station to consider. Default: the first one
#' @param noNA TRUE to show only errors associated with full lines (old behavior)
#' @param ... Additional arguments to be passed to graphics::plot method.
#'
#' @export
plotCloud = function(r, thresh=30, hintThresh=c(30,50,80), station=1, noNA=TRUE, ...)
{
	if (is.character(station))
		station = match(station, r$stations)

	XY = subset(r$data, subset = (Station == station), select = c(r$experts,"Measure","Prediction"))
	Y = XY[,"Measure"]
	hatY = XY[,"Prediction"]
	indices = !is.na(Y) & !is.na(hatY)
	if (noNA)
	{
		X = XY[,names(XY) %in% r$experts]
		indices = indices & getNoNAindices(X)
	}
	Y = Y[indices]
	hatY = hatY[indices]

	indics = getIndicators(r, thresh, station, noNA)

	par(mar=c(5,5,3,2), cex=1.5)
	plot(Y, hatY, xlab="Measured PM10", ylab="Predicted PM10",
		cex.lab=1.6, cex.axis=1.5, xlim=c(0,120), ylim=c(0,120), ...)
	abline(0,1,h=hintThresh,v=hintThresh,col=2,lwd=2)
#	legend("topleft",legend=c(paste("EV   ",indics$EV),paste("RMSE ",indics$RMSE)),cex=1.2)
	legend("topleft",legend=paste("RMSE ",indics$RMSE))
	legend("bottomright",legend=c(paste("TS ",indics$TS)))
}
Commit	Line	Data
a961f8a1 BA	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	}