f522f0fdc0df4edef9b2cd7769b8ece45a63e0d9
[aggexp.git] / pkg / R / z_plotHelper.R
1 #' @include z_runAlgorithm.R
2
3 #' @title Get best expert index
4 #'
5 #' @description Return the weights corresponding to the best expert (...0,1,0...)
6 #'
7 #' @param r Output of \code{\link{runAlgorithm}}
8 #'
9 #' @export
10 getBestExpert = function(r)
11 {
12 X = as.matrix(r$data[,names(r$data) %in% r$experts])
13 Y = r$data[,"Measure"]
14
15 bestIndex = which.min(colMeans(abs(X - Y)^2, na.rm=TRUE))
16 res = rep(0.0, length(r$experts))
17 res[bestIndex] = 1.0
18 return (res)
19 }
20
21 #' @title Get best convex combination
22 #'
23 #' @description Return the weights p minimizing the quadratic error ||X*p-Y||^2 under convexity contraint.
24 #'
25 #' @param r Output of \code{\link{runAlgorithm}}
26 #'
27 #' @export
28 getBestConvexCombination = function(r)
29 {
30 X = r$data[,r$experts]
31 Y = as.double(r$data[,"Measure"])
32 indices = getNoNAindices(X) & !is.na(Y)
33 X = as.matrix(X[indices,])
34 Y = Y[indices]
35
36 K = length(r$experts)
37 return (constrOptim(theta=rep(1.0/K,K),
38 method="Nelder-Mead", #TODO: others not better... why?
39 f=function(p){return(sum((X%*%p-Y)^2))},
40 grad=NULL, #function(p){return(2.*t(X)%*%(X%*%p-Y))},
41 ui=rbind(rep(1.,K),rep(-1.,K),diag(K)), ci=c(0.99999,-1.00001, rep(0.,K)),
42 control=list(ndeps=1e-3,maxit=10000))$par)
43 }
44
45 #' @title Get best linear combination
46 #'
47 #' @description Return the weights u minimizing the quadratic error ||r$X*u-r$Y||^2
48 #'
49 #' @param r Output of \code{\link{runAlgorithm}}
50 #'
51 #' @export
52 getBestLinearCombination = function(r)
53 {
54 X = r$data[,r$experts]
55 Y = r$data[,"Measure"]
56 indices = getNoNAindices(X) & !is.na(Y)
57 X = as.matrix(X[indices,])
58 Y = Y[indices]
59
60 return (mpPsInv(X) %*% Y)
61 }
62
63 #' @title Get statistical indicators
64 #'
65 #' @description Return respectively the TS, FA, MA, RMSE, EV indicators in a list.
66 #'
67 #' @param r Output of \code{\link{runAlgorithm}}
68 #' @param thresh Threshold to compute alerts indicators.
69 #' @param station Name or index of the station to consider. Default: the first one
70 #' @param noNA TRUE to show only errors associated with full lines (old behavior)
71 #'
72 #' @export
73 getIndicators = function(r, thresh, station=1, noNA=TRUE)
74 {
75 if (is.character(station))
76 station = match(station, r$stations)
77
78 #TODO: duplicated block (same in plotCloud())
79 XY = subset(r$data, subset = (Station == station), select = c(r$experts,"Measure","Prediction"))
80 Y = XY[,"Measure"]
81 hatY = XY[,"Prediction"]
82 indices = !is.na(Y) & !is.na(hatY)
83 if (noNA)
84 {
85 X = XY[,names(XY) %in% r$experts]
86 indices = indices & getNoNAindices(X)
87 }
88 Y = Y[indices]
89 hatY = hatY[indices]
90
91 RMSE = round(sqrt(sum((Y - hatY)^2) / length(Y)),2)
92 EV = round(1 - var(Y-hatY) / var(Y), 2)
93 A = sum(hatY >= thresh & Y >= thresh, na.rm=TRUE) #right alarm
94 B = sum(hatY >= thresh & Y < thresh, na.rm=TRUE) #false alarm
95 C = sum(hatY < thresh & Y >= thresh, na.rm=TRUE) #missed alert
96 TS = round(A/(A+B+C),2)
97 FA = B/(A+B)
98 MA = C/(A+C)
99 return (list("TS"=TS, "FA"=FA, "MA"=MA, "RMSE"=RMSE, "EV"=EV))
100 }