Adding functionalities for analyses.

NAMESPACE

@@ -1,6 +1,8 @@
 # Generated by roxygen2: do not edit by hand
 
 export(crps)
+export(crpsContribClasses)
 export(crpsNbAnalogs)
 export(optimalNbAnalog)
+export(optimalNbAnalogOfAnalogs)
 export(parseNcOutputs)

R/evalNbAnalogs.R

@@ -10,7 +10,7 @@
 #'
 #' @examples
 #' \dontrun{
-#' data <- atmoswing::parseNcOutputs('optimizer-outputs/1/results', 1, 'validation')
+#' data <- atmoswing::parseNcOutputs('optim/1/results', 1, 'validation')
 #' res <- atmoswing::crpsNbAnalogs(data, 9)
 #' }
 #' 
@@ -57,21 +57,21 @@ crpsNbAnalogs <- function(A, filter.size = 9) {
 #'
 #' @param crps.matrix Results of atmoswing::crpsNbAnalogs.
 #'
-#' @return Vector of the optimal number of analogues for each target date..
+#' @return Vector of the optimal number of analogues for each target date.
 #'
 #' @examples
 #' \dontrun{
-#' data <- atmoswing::parseNcOutputs('optimizer-outputs/1/results', 1, 'validation')
-#' res.crps <- atmoswing::crpsNbAnalogs(data, 9)
-#' res.nb <- atmoswing::optimalNbAnalog(res.crps$crps.smoothed, 9)
+#' data <- atmoswing::parseNcOutputs('optim/1/results', 1, 'validation')
+#' res.crps <- atmoswing::crpsNbAnalogs(data)
+#' res.nb <- atmoswing::optimalNbAnalog(res.crps$crps.smoothed)
 #' }
 #' 
 #' @export
 #' 
 optimalNbAnalog <- function(crps.matrix) {
 
   nb.min.first <- apply(crps.matrix, 1, which.min)
-  nb.min.med <- apply(crps.matrix, 1, function(x) median(which(x==min(x, na.rm = TRUE))))
+  nb.min.med <- apply(crps.matrix, 1, function(x) stats::median(which(x==min(x, na.rm = TRUE))))
   nb.min.last <- apply(crps.matrix, 1, function(x) max(which(x==min(x, na.rm = TRUE))))
 
   # Pack results
@@ -81,3 +81,39 @@ optimalNbAnalog <- function(crps.matrix) {
 
   return(res)
 }
+
+
+#' Assign the optimal number of analogues to every analogue date
+#'
+#' Assign the optimal number of analogues to every analogue date. This only
+#' works on the calibration period, where the optimal number of analogues can
+#' be assessed for the non independent dates.
+#'
+#' @param A Results of AtmoSwing as parsed by atmoswing::parseNcOutputs.
+#' @param target.a.nb Results of atmoswing::optimalNbAnalog
+#'
+#' @return Matrix of the optimal number of analogues for each analogue date.
+#'
+#' @examples
+#' \dontrun{
+#' data <- atmoswing::parseNcOutputs('optim/1/results', 1, 'validation')
+#' res.crps <- atmoswing::crpsNbAnalogs(data)
+#' target.a.nb <- atmoswing::optimalNbAnalog(res.crps$crps.smoothed)
+#' res <- atmoswing::optimalNbAnalogOfAnalogs(data, target.a.nb$nb.min.med)
+#' }
+#' 
+#' @export
+#' 
+optimalNbAnalogOfAnalogs <- function(A, target.a.nb) {
+
+  analog.a.nb <- matrix(NA, nrow = nrow(A$analog.dates.MJD), ncol = ncol(A$analog.dates.MJD))
+
+  for (i in 1:length(A$target.dates.MJD)){
+    targ.date <- A$target.dates.MJD[i]
+    targ.nb.analog <- target.a.nb[i]
+
+    analog.a.nb[A$analog.dates.MJD==targ.date] <- targ.nb.analog
+  }
+
+  return(analog.a.nb)
+}

R/performanceScore.R

@@ -59,4 +59,53 @@ crps <- function(x, x0, a=0.44, b=0.12) {
   }
 
   return(res)
-} 
+}
+
+#' Get the contribution to the CRPS per predictand classes.
+#'
+#' Get the contribution to the CRPS for different predictand classes.
+#'
+#' @param crps.vect A vector of the CRPS value for every target day.
+#' @param obs.vect A vector of the observed value (same length as crps.vect).
+#'
+#' @return Contribution to the CRPS (%).
+#'
+#' @examples
+#' \dontrun{
+#' data <- atmoswing::parseNcOutputs('optim/1/results', 1, 'validation')
+#' crps.matrix <- atmoswing::crpsNbAnalogs(data)
+#' res <- atmoswing::crpsContribClasses(crps.matrix$crps[, 30], data$target.values.raw)
+#' barplot(res$labels)
+#' title('Contribution to the CRPS (not cumulative)')
+#' }
+#' 
+#' @export
+#' 
+crpsContribClasses <- function(crps.vect, obs.vect) {
+
+  contrib <- matrix(NA, ncol = 7)
+
+  contrib[1] <- sum(crps.vect[obs.vect == 0])
+  contrib[2] <- sum(crps.vect[obs.vect > 0
+                              & obs.vect <= stats::quantile(obs.vect, probs = 0.8)])
+  contrib[3] <- sum(crps.vect[obs.vect > stats::quantile(obs.vect, probs = 0.8)
+                              & obs.vect <= stats::quantile(obs.vect, probs = 0.9)])
+  contrib[4] <- sum(crps.vect[obs.vect > stats::quantile(obs.vect, probs = 0.9)
+                              & obs.vect <= stats::quantile(obs.vect, probs = 0.95)])
+  contrib[5] <- sum(crps.vect[obs.vect > stats::quantile(obs.vect, probs = 0.95)
+                              & obs.vect <= stats::quantile(obs.vect, probs = 0.97)])
+  contrib[6] <- sum(crps.vect[obs.vect > stats::quantile(obs.vect, probs = 0.97)
+                              & obs.vect <= stats::quantile(obs.vect, probs = 0.99)])
+  contrib[7] <- sum(crps.vect[obs.vect > stats::quantile(obs.vect, probs = 0.99)])
+
+  contrib <- 100 * contrib / sum(contrib)
+
+  labels <- c("P = 0", "0 < P <= q.8", "q.8 < P <= q.9", "q.9 < P <= q.95",
+              "q.95 < P <= q.97", "q.97 < P <= q.99", "P > q.99")
+
+  # Pack results
+  res <- list(contrib = contrib,
+              labels = labels)
+
+  return(res)
+}

tests/testthat/testOptimalNbAnalog.R

@@ -1,7 +1,7 @@
 library(atmoswing)
 context("Optimal nb analogues")
 
-test_that("Optimal nb of analogues number is correct", {
+test_that("Optimal nb of analogues is correct", {
   A <- atmoswing::parseNcOutputs(file.path('test_files', 'optim', '1', 'results'),
                                  1, 'calibration')
 

tests/testthat/testOptimalNbAnalogOfAnalogs.R

@@ -0,0 +1,17 @@
+library(atmoswing)
+context("Optimal nb analogues of analogues")
+
+test_that("Optimal nb of analogues of analogues is correct", {
+  A <- atmoswing::parseNcOutputs(file.path('test_files', 'optim', '1', 'results'),
+                                 1, 'calibration')
+
+  crps.nb.analogs <- atmoswing::crpsNbAnalogs(A, 5)
+
+  target.a.nb <- atmoswing::optimalNbAnalog(crps.nb.analogs$crps.smoothed)
+
+  res <- atmoswing::optimalNbAnalogOfAnalogs(A, target.a.nb$nb.min.med)
+
+  expect_equal(res[5, 2], 3)
+  expect_equal(res[6, 5], 5.5)
+  expect_equal(res[24, 10], 8)
+})