Commit | Line | Data |
---|---|---|
3a38473a BA |
1 | #' Plot curves |
2 | #' | |
3 | #' Plot a range of curves in data. | |
4 | #' | |
5 | #' @inheritParams computeError | |
6 | #' @param indices Range of indices (integers or dates) | |
7 | #' | |
8 | #' @export | |
9 | plotCurves <- function(data, indices=seq_len(data$getSize())) | |
10 | { | |
11 | series = data$getSeries(indices) | |
12 | yrange = quantile(series, probs=c(0.025,0.975), na.rm=TRUE) | |
13 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) | |
14 | matplot(series, type="l", ylim=yrange, xlab="Time (hours)", ylab="PM10") | |
15 | } | |
16 | ||
17 | #' Plot error | |
18 | #' | |
19 | #' Draw error graphs, potentially from several runs of \code{computeForecast()}. | |
20 | #' | |
21 | #' @param err Error as returned by \code{computeError()} | |
22 | #' @param cols Colors for each error (default: 1,2,3,...) | |
23 | #' @param agg Aggregation level ("day", "week" or "month") | |
24 | #' | |
25 | #' @seealso \code{\link{plotCurves}}, \code{\link{plotPredReal}}, | |
26 | #' \code{\link{plotSimils}}, \code{\link{plotFbox}}, \code{\link{computeFilaments}}, | |
27 | #' \code{\link{plotFilamentsBox}}, \code{\link{plotRelVar}} | |
28 | #' | |
29 | #' @export | |
30 | plotError <- function(err, cols=seq_along(err), agg="day") | |
31 | { | |
32 | if (!is.null(err$abs)) | |
33 | err = list(err) | |
34 | par(mfrow=c(2,2), mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) | |
35 | L = length(err) | |
36 | ||
37 | yrange = range( sapply(1:L, function(i) err[[i]]$abs$day), na.rm=TRUE ) | |
38 | matplot(sapply( seq_len(L), function(i) err[[i]]$abs$day ), type="l", | |
39 | xlab="Time (hours)", ylab="Mean |y - y_hat|", ylim=yrange, col=cols, lwd=2, lty=1) | |
40 | ||
41 | agg_curves <- sapply( seq_len(L), function(i) { | |
42 | curve <- err[[i]]$abs$indices | |
43 | delta <- if (agg=="day") 1 else if (agg=="week") 7 else if (agg=="month") 30 | |
44 | vapply( seq(1,length(curve),delta), function(i) { | |
45 | mean(curve[i:(i+delta-1)], na.rm=TRUE) | |
46 | }, vector("double",1), USE.NAMES=FALSE ) | |
47 | }) | |
48 | yrange = range(agg_curves, na.rm=TRUE) | |
49 | matplot(agg_curves, type="l", xlab=paste("Time (",agg,"s)", sep=""), | |
50 | ylab="Mean |y - y_hat|", ylim=yrange, col=cols, lwd=2, lty=1) | |
51 | ||
52 | yrange = range( sapply(1:L, function(i) err[[i]]$MAPE$day), na.rm=TRUE ) | |
53 | matplot(sapply( seq_len(L), function(i) err[[i]]$MAPE$day ), type="l", | |
54 | xlab="Time (hours)", ylab="Mean MAPE", ylim=yrange, col=cols, lwd=2, lty=1) | |
55 | ||
56 | agg_curves <- sapply( seq_len(L), function(i) { | |
57 | curve <- err[[i]]$MAPE$indices | |
58 | delta <- if (agg=="day") 1 else if (agg=="week") 7 else if (agg=="month") 30 | |
59 | vapply( seq(1,length(curve),delta), function(i) { | |
60 | mean(curve[i:(i+delta-1)], na.rm=TRUE) | |
61 | }, vector("double",1), USE.NAMES=FALSE ) | |
62 | }) | |
63 | yrange = range(agg_curves, na.rm=TRUE) | |
64 | matplot(agg_curves, type="l", xlab=paste("Time (",agg,"s)", sep=""), | |
65 | ylab="Mean MAPE", ylim=yrange, col=cols, lwd=2, lty=1) | |
66 | } | |
67 | ||
68 | #' Plot measured / predicted | |
69 | #' | |
70 | #' Plot measured curve (in black) and predicted curve (in blue). | |
71 | #' | |
72 | #' @inheritParams computeError | |
73 | #' @param index Index in forecasts (integer or date) | |
74 | #' | |
75 | #' @export | |
76 | plotPredReal <- function(data, pred, index) | |
77 | { | |
78 | prediction = pred$getForecast(index) | |
79 | measure = data$getSerie( pred$getIndexInData(index) )[1:length(pred$getForecast(1))] | |
80 | ||
81 | # Remove the common part, where prediction == measure | |
82 | dot_mark <- ifelse(prediction[1]==measure[1], | |
83 | which.max(seq_along(prediction)[prediction==measure]), 0) | |
84 | prediction = prediction[(dot_mark+1):length(prediction)] | |
85 | measure = measure[(dot_mark+1):length(measure)] | |
86 | ||
87 | yrange = range(measure, prediction) | |
88 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=3) | |
89 | plot(measure, type="l", ylim=yrange, xlab="Time (hours)", ylab="PM10") | |
90 | par(new=TRUE) | |
91 | plot(prediction, type="l", col="#0000FF", ylim=yrange, xlab="", ylab="") | |
92 | } | |
93 | ||
94 | #' Plot similarities | |
95 | #' | |
96 | #' Plot histogram of similarities (weights), for 'Neighbors' method. | |
97 | #' | |
98 | #' @inheritParams computeError | |
99 | #' @param index Index in forecasts (integer or date) | |
100 | #' | |
101 | #' @export | |
102 | plotSimils <- function(pred, index) | |
103 | { | |
104 | weights = pred$getParams(index)$weights | |
105 | if (is.null(weights)) | |
106 | stop("plotSimils only works on 'Neighbors' forecasts") | |
107 | par(mfrow=c(1,2), mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) | |
108 | small_weights = weights[ weights < 1/length(weights) ] | |
109 | large_weights = weights[ weights >= 1/length(weights) ] | |
110 | hist(small_weights, nclass=25, main="", xlab="Weight < 1/N", ylab="Count") | |
111 | hist(large_weights, nclass=25, main="", xlab="Weight >= 1/N", ylab="Count") | |
112 | } | |
113 | ||
114 | #' Functional boxplot | |
115 | #' | |
116 | #' Draw the functional boxplot on the left, and bivariate plot on the right. | |
117 | #' | |
118 | #' @inheritParams computeError | |
119 | #' @inheritParams plotCurves | |
120 | #' | |
121 | #' @export | |
122 | plotFbox <- function(data, indices=seq_len(data$getSize())) | |
123 | { | |
124 | if (!requireNamespace("rainbow", quietly=TRUE)) | |
125 | stop("Functional boxplot requires the rainbow package") | |
126 | ||
127 | series_matrix = data$getSeries(indices) | |
128 | # Remove series with NAs | |
129 | no_NAs_indices = sapply( 1:ncol(series_matrix), | |
130 | function(i) all(!is.na(series_matrix[,i])) ) | |
131 | series_matrix = series_matrix[,no_NAs_indices] | |
132 | ||
133 | series_fds = rainbow::fds(seq_len(nrow(series_matrix)), series_matrix) | |
134 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) | |
135 | rainbow::fboxplot(series_fds, "functional", "hdr", xlab="Time (hours)", ylab="PM10", | |
136 | plotlegend=FALSE, lwd=2) | |
137 | rainbow::fboxplot(series_fds, "bivariate", "hdr", plotlegend=FALSE) | |
138 | } | |
139 | ||
140 | #' Compute filaments | |
141 | #' | |
142 | #' Obtain similar days in the past, and (optionally) plot them -- as black as distances | |
143 | #' are small. | |
144 | #' | |
145 | #' @inheritParams computeError | |
146 | #' @param index Index in forecast (integer or date) | |
147 | #' @param limit Number of neighbors to consider | |
148 | #' @param plot Should the result be plotted? | |
149 | #' | |
150 | #' @return A list with | |
151 | #' \itemize{ | |
152 | #' \item index : index of the current serie ('today') | |
153 | #' \item neighb_indices : indices of its neighbors | |
154 | #' \item colors : colors of neighbors curves (shades of gray) | |
155 | #' } | |
156 | #' | |
157 | #' @export | |
158 | computeFilaments <- function(data, pred, index, limit=60, plot=TRUE) | |
159 | { | |
160 | weights <- pred$getParams(index)$weights | |
161 | if (is.null(weights) || is.na(pred$getParams(index)$weights[1])) | |
162 | stop("computeFilaments requires a serie without NAs") | |
163 | ||
164 | nn <- min(limit, length(weights)) | |
165 | sorted_dists = sort(-log(weights), index.return=TRUE) | |
166 | # Compute colors for each neighbor (from darkest to lightest), if weights differ | |
167 | if ( any( weights != weights[1] ) ) | |
168 | { | |
169 | min_dist = min(sorted_dists$x[1:nn]) | |
170 | max_dist = max(sorted_dists$x[1:nn]) | |
171 | color_values = floor(19.5*(sorted_dists$x[1:nn]-min_dist)/(max_dist-min_dist)) + 1 | |
172 | colors = gray.colors(20,0.1,0.9)[color_values] #TODO: 20 == magic number | |
173 | } | |
174 | else | |
175 | colors <- rep(colors()[17], length(weights)) | |
176 | ||
177 | if (plot) | |
178 | { | |
179 | # Complete series with (past and present) tomorrows | |
180 | ref_serie = c( data$getCenteredSerie( pred$getIndexInData(index)-1 ), | |
181 | data$getCenteredSerie( pred$getIndexInData(index) ) ) | |
182 | centered_series = rbind( | |
183 | data$getCenteredSeries( pred$getParams(index)$indices-1 ), | |
184 | data$getCenteredSeries( pred$getParams(index)$indices ) ) | |
185 | yrange = range( ref_serie, | |
186 | quantile(centered_series, probs=c(0.025,0.975), na.rm=TRUE) ) | |
187 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=2) | |
188 | for (i in nn:1) | |
189 | { | |
190 | plot(centered_series[,sorted_dists$ix[i]], ylim=yrange, type="l", col=colors[i], | |
191 | xlab=ifelse(i==1,"Time (hours)",""), ylab=ifelse(i==1,"Centered PM10","")) | |
192 | par(new=TRUE) | |
193 | } | |
194 | # Also plot ref curve, in red | |
195 | plot(ref_serie, ylim=yrange, type="l", col="#FF0000", xlab="", ylab="") | |
196 | dot_mark <- 0.5 + which.max( pred$getForecast(1) == | |
197 | data$getSerie( pred$getIndexInData(1) )[1:length(pred$getForecast(1))] ) | |
198 | abline(v=24+dot_mark, lty=2, col=colors()[56], lwd=1) | |
199 | } | |
200 | ||
201 | list( | |
202 | "index"=pred$getIndexInData(index), | |
203 | "neighb_indices"=pred$getParams(index)$indices[sorted_dists$ix[1:nn]], | |
204 | "colors"=colors) | |
205 | } | |
206 | ||
207 | #' Functional boxplot on filaments | |
208 | #' | |
209 | #' Draw the functional boxplot on filaments obtained by \code{computeFilaments()}. | |
210 | #' | |
211 | #' @inheritParams computeError | |
212 | #' @param fil Output of \code{computeFilaments} | |
213 | #' @param predict_from First predicted time step | |
214 | #' | |
215 | #' @export | |
216 | plotFilamentsBox = function(data, fil, predict_from) | |
217 | { | |
218 | if (!requireNamespace("rainbow", quietly=TRUE)) | |
219 | stop("Functional boxplot requires the rainbow package") | |
220 | ||
221 | series_matrix = rbind( | |
222 | data$getSeries(fil$neighb_indices-1), data$getSeries(fil$neighb_indices) ) | |
223 | series_fds = rainbow::fds(seq_len(nrow(series_matrix)), series_matrix) | |
224 | ||
225 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) | |
226 | rainbow::fboxplot(series_fds, "functional", "hdr", xlab="Time (hours)", ylab="PM10", | |
227 | plotlegend=FALSE, lwd=2) | |
228 | ||
229 | # "Magic": http://stackoverflow.com/questions/13842560/get-xlim-from-a-plot-in-r | |
230 | usr <- par("usr") | |
231 | yr <- (usr[4] - usr[3]) / 27 | |
232 | par(new=TRUE) | |
233 | plot(c(data$getSerie(fil$index-1),data$getSerie(fil$index)), type="l", lwd=2, lty=2, | |
234 | ylim=c(usr[3] + yr, usr[4] - yr), xlab="", ylab="") | |
235 | abline(v=24+predict_from-0.5, lty=2, col=colors()[56]) | |
236 | } | |
237 | ||
238 | #' Plot relative conditional variability / absolute variability | |
239 | #' | |
240 | #' Draw the relative conditional variability / absolute variability based on filaments | |
241 | #' obtained by \code{computeFilaments()}. | |
242 | #' | |
243 | #' @inheritParams computeError | |
244 | #' @inheritParams plotFilamentsBox | |
245 | #' | |
246 | #' @export | |
247 | plotRelVar = function(data, fil, predict_from) | |
248 | { | |
249 | ref_var = c( apply(data$getSeries(fil$neighb_indices-1),1,sd), | |
250 | apply(data$getSeries(fil$neighb_indices),1,sd) ) | |
251 | tdays = .getNoNA2(data, 2, fil$index) | |
252 | global_var = c( | |
253 | apply(data$getSeries(tdays-1),1,sd), | |
254 | apply(data$getSeries(tdays),1,sd) ) | |
255 | ||
256 | yrange = range(ref_var, global_var) | |
257 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) | |
258 | plot(ref_var, type="l", col=1, lwd=3, ylim=yrange, | |
259 | xlab="Time (hours)", ylab="Standard deviation") | |
260 | par(new=TRUE) | |
261 | plot(global_var, type="l", col=2, lwd=3, ylim=yrange, xlab="", ylab="") | |
262 | abline(v=24+predict_from-0.5, lty=2, col=colors()[56]) | |
263 | } |