R/Eno.R

In s2dverification: Set of Common Tools for Forecast Verification

#'Computes Effective Sample Size With Classical Method
#'
#'Computes the effective number of independent values along the posdim 
#'dimension of a matrix.\cr
#'This effective number of independent observations can be used in 
#'statistical/inference tests.\cr
#'Based on eno function from Caio Coelho from rclim.txt.
#'
#'@param obs Matrix of any number of dimensions up to 10.
#'@param posdim Dimension along which to compute the effective sample size.
#'
#'@return Same dimensions as var but without the posdim dimension.
#'
#'@keywords datagen
#'@author History:\cr
#'0.1  -  2011-05  (V. Guemas)  -  Original code\cr
#'1.0  -  2013-09  (N. Manubens)  -  Formatting to R CRAN  
#'@examples
#'# See examples on Load() to understand the first lines in this example
#'  \dontrun{
#'data_path <- system.file('sample_data', package = 's2dverification')
#'exp <- list(
#'         name = 'experiment',
#'         path = file.path(data_path, 'model/$EXP_NAME$/monthly_mean',
#'                          '$VAR_NAME$_3hourly/$VAR_NAME$_$START_DATES$.nc')
#'       )
#'obs <- list(
#'         name = 'observation',
#'         path = file.path(data_path, 'observation/$OBS_NAME$/monthly_mean',
#'                          '$VAR_NAME$/$VAR_NAME$_$YEAR$$MONTH$.nc')
#'       )
#'# Now we are ready to use Load().
#'startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
#'sampleData <- Load('tos', list(exp), list(obs), startDates,
#'                   leadtimemin = 1, leadtimemax = 4, output = 'lonlat',
#'                   latmin = 27, latmax = 48, lonmin = -12, lonmax = 40)
#'  }
#'  \dontshow{
#'startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
#'sampleData <- s2dverification:::.LoadSampleData('tos', c('experiment'), 
#'                                                c('observation'), startDates,
#'                                                output = 'lonlat', 
#'                                                latmin = 27, latmax = 48, 
#'                                                lonmin = -12, lonmax = 40) 
#'  }
#'sampleData$mod <- Season(sampleData$mod, 4, 11, 1, 12)
#'eno <- Eno(sampleData$mod[1, 1, , 1, , ], 1)
#'PlotEquiMap(eno, sampleData$lon, sampleData$lat)
#'
#'@importFrom stats acf na.pass
#'@export
Eno <- function(obs, posdim) {
  dimsvar <- dim(obs)
  if (is.null(dimsvar)) {
    dimsvar <- length(obs)
  }
  enlobs <- Enlarge(obs, 10)
  outdim <- c(dimsvar, array(1, dim = (10 - length(dimsvar))))
  posaperm <- 1:10
  posaperm[posdim] <- 1
  posaperm[1] <- posdim
  enlobs <- aperm(enlobs, posaperm)
  dimsaperm <- outdim[posaperm]
  #
  #  Loop on all dimensions to compute effective number of observations 
  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  #
  enleno <- array(dim = c(1, dimsaperm[2:10]))
  for (j2 in 1:dimsaperm[2]) {
    for (j3 in 1:dimsaperm[3]) {
      for (j4 in 1:dimsaperm[4]) {
        for (j5 in 1:dimsaperm[5]) {
          for (j6 in 1:dimsaperm[6]) {
            for (j7 in 1:dimsaperm[7]) {
              for (j8 in 1:dimsaperm[8]) {
                for (j9 in 1:dimsaperm[9]) {
                  for (j10 in 1:dimsaperm[10]) {
                    tmp <- enlobs[, j2, j3, j4, j5, j6, j7, j8, j9, j10]
                    if (length(sort(tmp)) > 1 ) {
                      n <- length(sort(tmp))
                      a <- acf(tmp, lag.max = n - 1, plot = FALSE, 
                               na.action = na.pass)$acf[2:n, 1, 1]
                      s <- 0
                      for (k in 1:(n - 1)) {
                        s <- s + (((n - k) / n) * a[k])
                      }
                      enleno[1, j2, j3, j4, j5, j6, j7, j8, j9, 
                             j10] <- min(n / (1 + (2 * s)), n)
                    }
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  #
  #  Back to the original dimensions
  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  #
  #dimsvar <- dimsvar[-posdim]
  if (length(dimsvar) == 1) {
    dimsvar <- 1
  } else {
    dimsvar <- dimsvar[-posdim]
  }
  effnumobs <- array(dim = dimsvar)
  effnumobs[] <- aperm(enleno, posaperm)
  #
  #  Outputs
  # ~~~~~~~~~
  #
  effnumobs
}