projects / aggexp.git / blame
| Commit | Line | Data |
|---|---|---|
| a961f8a1 BA |
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 | } |