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 }