diff --git a/DESCRIPTION b/DESCRIPTION
index 9a5a77142ea2554d0164a89dbcf9224b1ff7d99c..ecf2fa02bd9fbef6bcb562d2a7a9275a1aa1f066 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -37,6 +37,7 @@ Imports:
     ncdf4,
     NbClust,
     multiApply (>= 2.1.1),
+    SpecsVerification (>= 0.5.0),
     easyNCDF
 Suggests:
     easyVerification,
diff --git a/NAMESPACE b/NAMESPACE
index 3b097ba5778c78dca252519bd7b86269ed230046..0e2aec97d063b1b2138af5c1b637b5f9f54a1b13 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,6 +5,7 @@ export(AMV)
 export(AnimateMap)
 export(Ano)
 export(Ano_CrossValid)
+export(BrierScore)
 export(CDORemap)
 export(Clim)
 export(Cluster)
@@ -60,10 +61,12 @@ export(Spectrum)
 export(TPI)
 export(ToyModel)
 export(Trend)
+export(UltimateBrier)
 export(clim.colors)
 export(clim.palette)
 import(GEOmap)
 import(NbClust)
+import(SpecsVerification)
 import(bigmemory)
 import(geomapdata)
 import(graphics)
diff --git a/R/BrierScore.R b/R/BrierScore.R
new file mode 100644
index 0000000000000000000000000000000000000000..1363f6155ff53a7c73cb8afd3d630f05fceb343a
--- /dev/null
+++ b/R/BrierScore.R
@@ -0,0 +1,370 @@
+#'Compute Brier score, its decomposition, and Brier skill score
+#'
+#'Compute the Brier score (BS) and the components of its standard decompostion
+#'with the two within-bin components described in Stephenson et al., (2008). It
+#'also returns the bias-corrected decomposition of the BS (Ferro and Fricker, 
+#'2012). BSS has the climatology as the reference forecast. 
+#'
+#'@param exp A vector or a numeric array with named dimensions. It should be 
+#'  the predicted probabilities which are within the range [0, 1] if memb_dim
+#'  doesn't exist. If it has memb_dim, the value should be 0 or 1, and the 
+#'  predicted probabilities will be computed by ensemble mean. The dimensions 
+#'  must at least have 'time_dim'. 
+#'  range [0, 1].
+#'@param obs A numeric array with named dimensions of the binary observations 
+#'  (0 or 1). The dimension must be the same as 'exp' except memb_dim, which is
+#'  optional. If it has 'memb_dim', then the length must be 1. The length of 
+#'  'dat_dim' can be different from 'exp' if it has.
+#'@param thresholds A numeric vector used to bin the forecasts. The default 
+#' value is \code{seq(0.1, 0.9, 0.1)}, which means that the bins are 
+#'  \code{[0, 0.1), [0.1, 0.2), ... [0.9, 1]}.
+#'@param time_dim A character string indicating the name of dimension along  
+#'  which Brier score is computed. The default value is 'sdate'.
+#'@param dat_dim A character string indicating the name of dataset dimension in
+#'  'exp' and 'obs'. The length of this dimension can be different between 
+#'  'exp' and 'obs'. The default value is NULL.
+#'@param memb_dim A character string of the name of the member dimension in 
+#'  'exp' (and 'obs', optional). The function will do the ensemble mean 
+#'  over this dimension. If there is no member dimension, set NULL. The default
+#'  value is NULL.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return 
+#'A list that contains:
+#'\item{$rel}{standard reliability}
+#'\item{$res}{standard resolution}
+#'\item{$unc}{standard uncertainty}  
+#'\item{$bs}{Brier score}
+#'\item{$bs_check_res}{rel - res + unc}
+#'\item{$bss_res}{res - rel / unc}
+#'\item{$gres}{generalized resolution}
+#'\item{$bs_check_gres}{rel - gres + unc}
+#'\item{$bss_gres}{gres - rel / unc}
+#'\item{$rel_bias_corrected}{bias - corrected rel}
+#'\item{$gres_bias_corrected}{bias - corrected gres}
+#'\item{$unc_bias_corrected}{bias - corrected unc}
+#'\item{$bss_bias_corrected}{gres_bias_corrected - rel_bias_corrected / unc_bias_corrected}
+#'\item{$nk}{number of forecast in each bin}
+#'\item{$fkbar}{average probability of each bin}
+#'\item{$okbar}{relative frequency that the observed event occurred}
+#'The data type and dimensions of the items depend on if the input 'exp' and 
+#''obs' are:\cr
+#'(a) Vectors\cr
+#'(b) Arrays with 'dat_dim' specified\cr
+#'(c) Arrays with no 'dat_dim' specified\cr
+#'Items 'rel', 'res', 'unc', 'bs', 'bs_check_res', 'bss_res', 'gres', 
+#''bs_check_gres', 'bss_gres', 'rel_bias_corrected', 'gres_bias_corrected', 
+#''unc_bias_corrected', and 'bss_bias_corrected' are (a) a number (b) an array
+#'with dimensions c(nexp, nobs, all the rest dimensions in 'exp' and 'obs' 
+#'expect 'time_dim' and 'memb_dim') (c) an array with dimensions of
+#''exp' and 'obs' except 'time_dim' and 'memb_dim'\cr
+#'Items 'nk', 'fkbar', and 'okbar' are (a) a vector of length of bin number 
+#'determined by 'threshold' (b) an array with dimensions c(nexp, nobs, 
+#'no. of bins, all the rest dimensions in 'exp' and 'obs' expect 'time_dim' and
+#''memb_dim') (c) an array with dimensions c(no. of bin, all the rest dimensions
+#'in 'exp' and 'obs' expect 'time_dim' and 'memb_dim')
+#' 
+#'@references
+#'Wilks (2006) Statistical Methods in the Atmospheric Sciences.\cr
+#'Stephenson et al. (2008). Two extra components in the Brier score decomposition. 
+#'  Weather and Forecasting, 23: 752-757.\cr
+#'Ferro and Fricker (2012). A bias-corrected decomposition of the BS. 
+#'  Quarterly Journal of the Royal Meteorological Society, DOI: 10.1002/qj.1924.
+#'
+#'@examples
+#'# Inputs are vectors
+#'exp <- runif(10)
+#'obs <- round(exp)
+#'x <- BrierScore(exp, obs)
+#'
+#'# Inputs are arrays
+#'example(Load)
+#'clim <- Clim(sampleData$mod, sampleData$obs)
+#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+#'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+#'bins_ano_exp <- ProbBins(ano_exp, thr = c(1/3, 2/3))
+#'bins_ano_obs <- ProbBins(ano_obs, thr = c(1/3, 2/3))
+#'res <- BrierScore(bins_ano_exp, MeanDims(bins_ano_obs, 'member'), memb_dim = 'member') 
+#'
+#'@import multiApply
+#'@export
+BrierScore <- function(exp, obs, thresholds = seq(0.1, 0.9, 0.1), time_dim = 'sdate',
+                       dat_dim = NULL, memb_dim = NULL, ncores = NULL) {
+
+  # Check inputs
+  ## exp and obs (1)
+  if (is.null(exp) | is.null(obs)) {
+    stop("Parameter 'exp' and 'obs' cannot be NULL.")
+  }
+  if (!is.numeric(exp) | !is.numeric(obs)) {
+    stop("Parameter 'exp' and 'obs' must be a numeric vector or a numeric array.")
+  }
+  if (is.null(dim(exp))) {  #is vector
+    dim(exp) <- c(length(exp))
+    names(dim(exp)) <- time_dim
+  }
+  if (is.null(dim(obs))) {  #is vector
+    dim(obs) <- c(length(obs))
+    names(dim(obs)) <- time_dim
+  }
+  if(any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) |
+     any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) {
+    stop("Parameter 'exp' and 'obs' must have dimension names.")
+  }
+  if (any(!obs %in% c(0, 1))) {
+    stop("Parameter 'obs' must be binary events (0 or 1).")
+  }
+  ## thresholds
+  if (!is.numeric(thresholds) | !is.vector(thresholds)) {
+    stop("Parameter 'thresholds' must be a numeric vector.")
+  }
+  if (any(thresholds <= 0 | thresholds >= 1)) {
+    stop("Parameter 'thresholds' must be between 0 and 1 as the bin-breaks.")
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) > 1) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) {
+    stop("Parameter 'time_dim' is not found in 'exp' and 'obs' dimension.")
+  }
+  ## dat_dim
+  if (!is.null(dat_dim)) {
+    if (!is.character(dat_dim) | length(dat_dim) > 1) {
+      stop("Parameter 'dat_dim' must be a character string.")
+    }
+    if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) {
+      stop("Parameter 'dat_dim' is not found in 'exp' and 'obs' dimension.")
+    }
+  }
+  ## memb_dim
+  if (!is.null(memb_dim)) {
+    if (!is.character(memb_dim) | length(memb_dim) > 1) {
+      stop("Parameter 'memb_dim' must be a character string.")
+    }
+    if (!memb_dim %in% names(dim(exp))) {
+      stop("Parameter 'memb_dim' is not found in 'exp' dimension.")
+    }
+    if (memb_dim %in% names(dim(obs))) {
+      if (dim(obs)[memb_dim] != 1) {
+        stop("The length of parameter 'memb_dim' in 'obs' must be 1.")
+      }
+    }
+  }
+  ## exp and obs (2)
+  if (is.null(memb_dim)) {
+    if (max(exp) > 1 | min(exp) < 0) {
+      stop("Parameter 'exp' must be within [0, 1] range.")
+    }
+  } else {
+    if (any(!exp %in% c(0, 1))) {
+      stop("Parameter 'exp' must be 0 or 1 if it has memb_dim.")
+    }
+  }
+  name_exp <- sort(names(dim(exp)))
+  name_obs <- sort(names(dim(obs)))
+  if (!is.null(memb_dim)) {
+    name_exp <- name_exp[-which(name_exp == memb_dim)]
+    if (memb_dim %in% name_obs) {
+    name_obs <- name_obs[-which(name_obs == memb_dim)]
+    }
+  }
+  if (!is.null(dat_dim)) {
+    name_exp <- name_exp[-which(name_exp == dat_dim)]
+    name_obs <- name_obs[-which(name_obs == dat_dim)]
+  }
+  if (any(name_exp != name_obs)) {
+    stop(paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
+                "of all the dimensions expect 'dat_dim' and 'memb_dim'."))
+  }
+  if (!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
+    stop(paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
+                "of all the dimensions expect 'dat_dim' and 'memb_dim'."))
+  }
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+      length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  } 
+
+  ###############################
+  # Calculate Brier score
+
+  ## ensemble mean
+  if (!is.null(memb_dim)) {
+    exp <- MeanDims(exp, memb_dim)
+    if (memb_dim %in% names(dim(obs))) {
+      obs <- MeanDims(obs, memb_dim)
+    }
+  }
+
+  if (is.null(dat_dim)) {
+    res <- Apply(list(exp, obs), 
+                 target_dims = list(c(time_dim), 
+                                    c(time_dim)),
+                 fun = .BrierScore, 
+                 thresholds = thresholds,
+                 ncores = ncores)
+  } else {
+    res <- Apply(list(exp, obs),
+                 target_dims = list(c(time_dim, dat_dim),
+                                    c(time_dim, dat_dim)),
+                 fun = .BrierScore,
+                 thresholds = thresholds,
+                 ncores = ncores)
+  }
+
+ return(res)
+}
+
+.BrierScore <- function(exp, obs, thresholds = seq(0.1, 0.9, 0.1)) {
+
+  # exp: [sdate] or [sdate, nexp]
+  # obs: [sdate] or [sdate, nobs]
+  if (length(dim(exp)) == 2) {
+    nexp <- as.numeric(dim(exp)[2])
+    nobs <- as.numeric(dim(obs)[2])
+    exp_ori <- exp
+    obs_ori <- obs
+    # Create empty arrays
+    arr_rel <- arr_res <- arr_unc <- arr_bs <- arr_bs_check_res <- arr_bss_res <- 
+      arr_gres <- arr_bs_check_gres <- arr_bss_gres <- arr_rel_bias_corrected <- 
+      arr_gres_bias_corrected  <- arr_unc_bias_corrected <- arr_bss_bias_corrected <- 
+      array(dim = c(nexp = nexp, nobs = nobs))
+    arr_nk <- arr_fkbar <- arr_okbar <- 
+      array(dim = c(nexp = nexp, nobs = nobs, bin = length(thresholds) + 1))
+
+  } else {
+    nexp <- 1
+    nobs <- 1
+  }
+
+  for (n_exp in 1:nexp) {
+    for (n_obs in 1:nobs) {
+      if (exists('exp_ori')) {
+        exp <- exp_ori[, n_exp]
+        obs <- obs_ori[, n_obs]
+      }
+      n <- length(exp)
+      nbins <- length(thresholds) + 1  # Number of bins
+      bins <- vector('list', nbins) #as.list(paste("bin", 1:nbins, sep = ""))
+      for (i in 1:nbins) {
+        if (i == 1) {
+          bins[[i]] <- list(which(exp >= 0 & exp < thresholds[i]))
+        } else if (i == nbins) {
+          bins[[i]] <- list(which(exp >= thresholds[i - 1] & exp <= 1))
+        } else {
+          bins[[i]] <- list(which(exp >= thresholds[i - 1] & exp < thresholds[i]))
+        }
+      }
+    
+      fkbar <- okbar <- nk <- array(0, dim = nbins)
+      for (i in 1:nbins) {
+        nk[i] <- length(bins[[i]][[1]])
+        fkbar[i] <- sum(exp[bins[[i]][[1]]]) / nk[i]
+        okbar[i] <- sum(obs[bins[[i]][[1]]]) / nk[i]
+      }
+        
+    #-----in old .BrierScore()---------
+    #    fkbar[fkbar == Inf] <- 0 
+    #    okbar[is.nan(okbar)] <- 0
+    #----------------------------------
+    
+      obar <- sum(obs) / length(obs)
+      relsum <- ressum <- term1 <- term2 <- 0
+      for (i in 1:nbins) {
+        if (nk[i] > 0) {
+          relsum <- relsum + nk[i] * (fkbar[i] - okbar[i])^2
+          ressum <- ressum + nk[i] * (okbar[i] - obar)^2
+          for (j in 1:nk[i]) {
+            term1 <- term1 + (exp[bins[[i]][[1]][j]] - fkbar[i])^2
+            term2 <- term2 + (exp[bins[[i]][[1]][j]] - fkbar[i]) * (obs[bins[[i]][[1]][j]] - okbar[i])
+          }
+        }
+      }
+      rel <- relsum / n
+      res <- ressum / n
+      unc <- obar * (1 - obar)
+      bs <- sum((exp - obs)^2) / n
+      bs_check_res <- rel - res + unc
+      bss_res <- (res - rel) / unc
+      gres <- res - term1 * (1 / n) + term2 * (2 / n)   # Generalized resolution
+      bs_check_gres <- rel - gres + unc                 # BS using GRES
+      bss_gres <- (gres - rel) / unc                    # BSS using GRES
+        
+      
+      # Estimating the bias-corrected components of the BS 
+      term3 <- array(0, nbins)
+      for (i in 1:nbins) {
+        term3[i] <- (nk[i] / (nk[i] - 1)) * okbar[i] * (1 - okbar[i])
+      }
+      term_a <- sum(term3,  na.rm = T) / n
+      term_b <- (obar * (1 - obar)) / (n - 1)
+      rel_bias_corrected <- rel - term_a
+      gres_bias_corrected <- gres - term_a + term_b
+      if (rel_bias_corrected < 0 || gres_bias_corrected < 0) {
+        rel_bias_corrected2 <- max(rel_bias_corrected, rel_bias_corrected - gres_bias_corrected, 0)
+        gres_bias_corrected2 <- max(gres_bias_corrected, gres_bias_corrected - rel_bias_corrected, 0)
+        rel_bias_corrected <- rel_bias_corrected2
+        gres_bias_corrected <- gres_bias_corrected2
+      }
+      unc_bias_corrected <- unc + term_b
+      bss_bias_corrected <- (gres_bias_corrected - rel_bias_corrected) / unc_bias_corrected
+      
+      #if (round(bs, 8) == round(bs_check_gres, 8) & round(bs_check_gres, 8) == round((rel_bias_corrected - gres_bias_corrected + unc_bias_corrected), 8)) {
+      #  cat("No error found  \ n")
+      #  cat("BS = REL - GRES + UNC = REL_lessbias - GRES_lessbias + UNC_lessbias  \ n")
+      #}
+       
+      # Add name for nk, fkbar, okbar
+      names(dim(nk)) <- 'bin'
+      names(dim(fkbar)) <- 'bin'
+      names(dim(okbar)) <- 'bin'
+    
+      if (exists('exp_ori')) {
+        arr_rel[n_exp, n_obs] <- rel
+        arr_res[n_exp, n_obs] <- res
+        arr_unc[n_exp, n_obs] <- unc
+        arr_bs[n_exp, n_obs] <- bs
+        arr_bs_check_res[n_exp, n_obs] <- bs_check_res
+        arr_bss_res[n_exp, n_obs] <- bss_res
+        arr_gres[n_exp, n_obs] <- gres
+        arr_bs_check_gres[n_exp, n_obs] <- bs_check_gres
+        arr_bss_gres[n_exp, n_obs] <- bss_gres
+        arr_rel_bias_corrected[n_exp, n_obs] <- rel_bias_corrected
+        arr_gres_bias_corrected[n_exp, n_obs] <- gres_bias_corrected
+        arr_unc_bias_corrected[n_exp, n_obs] <- unc_bias_corrected
+        arr_bss_bias_corrected[n_exp, n_obs] <- bss_bias_corrected
+        arr_nk[n_exp, n_obs, ] <- nk
+        arr_fkbar[n_exp, n_obs, ] <- fkbar
+        arr_okbar[n_exp, n_obs, ] <- okbar
+      }
+
+    }
+  }
+
+  if (exists('exp_ori')) {
+    res_list <- list(rel = arr_rel, res = arr_res, unc = arr_unc, bs = arr_bs, 
+                     bs_check_res = arr_bs_check_res, bss_res = arr_bss_res, 
+                     gres = arr_gres, bs_check_gres = arr_bs_check_gres,
+                     bss_gres = arr_bss_gres, rel_bias_corrected = arr_rel_bias_corrected,
+                     gres_bias_corrected = arr_gres_bias_corrected,
+                     unc_bias_corrected = arr_unc_bias_corrected,
+                     bss_bias_corrected = arr_bss_bias_corrected, nk = arr_nk, 
+                     fkbar = arr_fkbar, okbar = arr_okbar) #bins = list(bins), 
+  } else {
+
+    res_list <- list(rel = rel, res = res, unc = unc, bs = bs, bs_check_res = bs_check_res,
+                     bss_res = bss_res, gres = gres, bs_check_gres = bs_check_gres,
+                     bss_gres = bss_gres, rel_bias_corrected = rel_bias_corrected,
+                     gres_bias_corrected = gres_bias_corrected,
+                     unc_bias_corrected = unc_bias_corrected,
+                     bss_bias_corrected = bss_bias_corrected, nk = nk, fkbar = fkbar,
+                     okbar = okbar) #bins = list(bins), 
+  }
+
+  return(invisible(res_list))
+}
diff --git a/R/UltimateBrier.R b/R/UltimateBrier.R
new file mode 100644
index 0000000000000000000000000000000000000000..aeaddcdcadcdd9d25e3a283d1ec8d24534f969d0
--- /dev/null
+++ b/R/UltimateBrier.R
@@ -0,0 +1,318 @@
+#'Compute Brier scores
+#'
+#'Interface to compute probabilistic scores (Brier Score, Brier Skill Score) 
+#'from the forecast and observational data anomalies. It provides six types
+#'to choose.
+#'
+#'@param exp A numeric array of forecast anomalies with named dimensions that
+#'  at least include 'dat_dim', 'memb_dim', and 'time_dim'. It can be provided
+#'  by \code{Ano()}. 
+#'@param obs A numeric array of observational reference anomalies with named
+#'  dimensions that at least include 'dat_dim' and 'time_dim'. If it has 
+#'  'memb_dim', the length  must be 1. The dimensions should be consistent with
+#'  'exp' except 'dat_dim' and 'memb_dim'. It can be provided by \code{Ano()}. 
+#'@param dat_dim A character string indicating the name of the dataset 
+#'  dimension in 'exp' and 'obs'. The default value is 'dataset'.
+#'@param memb_dim A character string indicating the name of the member 
+#'  dimension in 'exp' (and 'obs') for ensemble mean calculation. The default
+#'  value is 'member'.
+#'@param time_dim A character string indicating the dimension along which to 
+#'  compute the probabilistic scores. The default value is 'sdate'.
+#'@param quantile A logical value to decide whether a quantile (TRUE) or a 
+#'  threshold (FALSE) is used to estimate the forecast and observed 
+#'  probabilities. If 'type' is 'FairEnsembleBS' or 'FairEnsembleBSS', it must
+#'  be TRUE. The default value is TRUE.
+#'@param thr A numeric vector to be used in probability calculation (for 'BS', 
+#'  'FairStartDatesBS', 'BSS', and 'FairStartDatesBSS') and binary event 
+#'  judgement (for 'FairEnsembleBS' and 'FairEnsembleBSS'). It is as
+#'  quantiles if 'quantile' is TRUE or as thresholds if 'quantile' is FALSE. 
+#'  The default value is \code{c(0.05, 0.95)} for 'quantile = TRUE'.
+#'@param type A character string of the desired score type. It can be the 
+#'  following values:
+#'\itemize{
+#'  \item{'BS': Simple Brier Score. Use SpecsVerification::BrierDecomp inside.}
+#'  \item{'FairEnsembleBS': Corrected Brier Score computed across ensemble 
+#'    members. Use SpecsVerification::FairBrier inside.}
+#'  \item{'FairStartDatesBS': Corrected Brier Score computed across starting 
+#'    dates. Use s2dv:::.BrierScore inside.}
+#'  \item{'BSS': Simple Brier Skill Score. Use s2dv:::.BrierScore inside.}
+#'  \item{'FairEnsembleBSS': Corrected Brier Skill Score computed across 
+#'    ensemble members. Use SpecsVerification::FairBrierSs inside.}
+#'  \item{'FairStartDatesBSS': Corrected Brier Skill Score computed across 
+#'    starting dates. Use s2dv:::.BrierScore inside.}
+#'}
+#'  The default value is 'BS'.
+#'@param decomposition A logical value to determine whether the decomposition 
+#'  of the Brier Score should be provided (TRUE) or not (FALSE). It is only 
+#'  used when 'type' is 'BS' or 'FairStartDatesBS'. The default value is TRUE.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return 
+#'If 'type' is 'BS' or 'FairStartDatesBS' and 'decomposition' is TRUE, the 
+#'output is a list of 4 arrays (see details below.) In other cases, the output 
+#'is an array of Brier scores or Brier skill scores. All the arrays have the 
+#'same dimensions:
+#'c(nexp, nobs, no. of bins, the rest dimensions of 'exp' except 'time_dim' and
+#''memb_dim'). 'nexp' and 'nobs' is the length of dataset dimension in 'exp'
+#'and 'obs' respectively.\cr
+#'The list of 4 includes: 
+#'  \itemize{
+#'    \item{$bs: Brier Score}
+#'    \item{$rel: Reliability component}
+#'    \item{$res: Resolution component}
+#'    \item{$unc: Uncertainty component}
+#'  }
+#'
+#'@examples
+#'  \dontshow{
+#'startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
+#'sampleData <- s2dv:::.LoadSampleData('tos', c('experiment'),
+#'                                     c('observation'), startDates,
+#'                                     leadtimemin = 1,
+#'                                     leadtimemax = 4,
+#'                                     output = 'lonlat',
+#'                                     latmin = 27, latmax = 48,
+#'                                     lonmin = -12, lonmax = 40)
+#'  }
+#'sampleData$mod <- Season(sampleData$mod, monini = 11, moninf = 12, monsup = 2)
+#'sampleData$obs <- Season(sampleData$obs, monini = 11, moninf = 12, monsup = 2)
+#'clim <- Clim(sampleData$mod, sampleData$obs)
+#'exp <- Ano(sampleData$mod, clim$clim_exp)
+#'obs <- Ano(sampleData$obs, clim$clim_obs)
+#'bs <- UltimateBrier(exp, obs)
+#'bss <- UltimateBrier(exp, obs, type = 'BSS')
+#'
+#'@import SpecsVerification plyr multiApply
+#'@export
+UltimateBrier <- function(exp, obs, dat_dim = 'dataset', memb_dim = 'member', time_dim = 'sdate', 
+                          quantile = TRUE, thr = c(5/100, 95/100), type = 'BS', 
+                          decomposition = TRUE, ncores = NULL) {
+
+  # Check inputs 
+  ## exp and obs (1)
+  if (is.null(exp) | is.null(obs)) {
+    stop("Parameter 'exp' and 'obs' cannot be NULL.")
+  }
+  if (!is.numeric(exp) | !is.numeric(obs)) {
+    stop("Parameter 'exp' and 'obs' must be a vector or a numeric array.")
+  }
+  if(any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) |
+     any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) {
+    stop("Parameter 'exp' and 'obs' must have dimension names.")
+  }
+  ## dat_dim
+  if (!is.character(dat_dim) | length(dat_dim) > 1) {
+    stop("Parameter 'dat_dim' must be a character string.")
+  }
+  if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) {
+    stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.")
+  }
+  ## memb_dim
+  if (!is.character(memb_dim) | length(memb_dim) > 1) {
+    stop("Parameter 'memb_dim' must be a character string.")
+  }
+  if (!memb_dim %in% names(dim(exp))) {
+    stop("Parameter 'memb_dim' is not found in 'exp' dimension.")
+  }
+  if (!memb_dim %in% names(dim(obs))) {
+    # Insert memb_dim into obs for the ease of later calculation
+    obs <- InsertDim(obs, posdim = 2, lendim = 1, name = memb_dim)
+  } else if (dim(obs)[memb_dim] != 1) {
+    stop("The length of parameter 'memb_dim' in 'obs' must be 1.") 
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) > 1) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) {
+    stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
+  }
+  ## exp and obs (2)
+  name_exp <- sort(names(dim(exp)))
+  name_obs <- sort(names(dim(obs)))
+  name_exp <- name_exp[-which(name_exp == memb_dim)]
+  name_obs <- name_obs[-which(name_obs == memb_dim)]
+  name_exp <- name_exp[-which(name_exp == dat_dim)]
+  name_obs <- name_obs[-which(name_obs == dat_dim)]
+  if (any(name_exp != name_obs)) {
+    stop(paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
+                "of all the dimensions expect 'dat_dim' and 'memb_dim'."))
+  }
+  if (!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
+    stop(paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
+                "of all the dimensions expect 'dat_dim' and 'memb_dim'."))
+  }
+  ## quantile
+  if (!is.logical(quantile) | length(quantile) > 1) {
+    stop("Parameter 'quantile' must be one logical value.")
+  }
+  ## thr
+  if (!is.numeric(thr) | !is.vector(thr)) {
+    stop("Parameter 'thr' must be a numeric vector.")
+  }
+  if (quantile) {
+    if (!all(thr < 1 & thr > 0)) {
+      stop("Parameter 'thr' must be between 0 and 1 when quantile is TRUE.")
+    }
+  }
+  if (!quantile & (type %in% c('FairEnsembleBSS', 'FairEnsembleBS'))) {
+    stop("Parameter 'quantile' must be TRUE if 'type' is 'FairEnsembleBSS' or 'FairEnsembleBS'.")
+  }
+  ## type
+  if (!(type %in% c("BS", "BSS", "FairEnsembleBS", "FairEnsembleBSS", "FairStartDatesBS", "FairStartDatesBSS"))) {
+    stop("Parameter 'type' must be one of 'BS', 'BSS', 'FairEnsembleBS', 'FairEnsembleBSS', 'FairStartDatesBS' or 'FairStartDatesBSS'.")
+  }
+  ## decomposition
+  if (!is.logical(decomposition) | length(decomposition) > 1) {
+    stop("Parameter 'decomposition' must be one logical value.")
+  }
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+      length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  } 
+
+  ###############################
+  # Calculate UltimateBrier 
+
+  if (type %in% c('FairEnsembleBSS', 'FairEnsembleBS')) {
+    res <- Apply(list(exp, obs),
+                 target_dims = list(c(time_dim, dat_dim, memb_dim),
+                                    c(time_dim, dat_dim, memb_dim)),
+                 fun = .UltimateBrier,
+                 thr = thr, type = type,
+                 decomposition = decomposition,
+                 ncores = ncores)$output1
+
+  } else {
+    # Calculate probablities by ProbBins() and ensemble mean first.
+    # The first dim will become 'bin' and memb_dim is gone.
+    exp <- MeanDims(
+             ProbBins(exp, thr = thr, time_dim = time_dim, memb_dim = memb_dim,
+                      quantile = quantile, ncores = ncores),
+             memb_dim)
+    obs <- MeanDims(
+             ProbBins(obs, thr = thr, time_dim = time_dim, memb_dim = memb_dim,
+                      quantile = quantile, ncores = ncores),
+             memb_dim)
+
+    res <- Apply(list(exp, obs),
+                 target_dims = list(c(time_dim, dat_dim),
+                                    c(time_dim, dat_dim)),
+                 fun = .UltimateBrier,
+                 thr = thr, type = type,
+                 decomposition = decomposition,
+                 ncores = ncores)
+
+    if (type %in% c('BSS', 'FairStartDatesBSS')) {
+      res <- res$output1
+    } else if (!decomposition) {
+      res <- res$bs
+    }
+  }
+
+  return(res)
+}
+
+.UltimateBrier <- function(exp, obs, thr = c(5/100, 95/100), type = 'BS',
+                           decomposition = TRUE) {
+  # If exp and obs are probablistics
+  # exp: [sdate, nexp]
+  # obs: [sdate, nobs]
+  # If exp and obs are anomalies
+  # exp: [sdate, nexp, memb]
+  # obs: [sdate, nobs, memb]
+  
+  #NOTE: 'thr' is used in 'FairEnsembleBSS' and 'FairEnsembleBS'. But if quantile = F and 
+  #      thr is real value, does it work?
+  if (type == 'FairEnsembleBSS') {
+    size_ens_ref <- prod(dim(obs)[c(1, 3)])
+    res <- array(dim = c(nexp = as.numeric(dim(exp)[2]), 
+                         nobs = as.numeric(dim(obs)[2]), 
+                         bin = length(thr) + 1))
+    for (n_exp in 1:dim(exp)[2]) {
+      for (n_obs in 1:dim(obs)[2]) {
+        ens_ref <- matrix(obs[, n_obs, 1], size_ens_ref, size_ens_ref, byrow = TRUE)
+        for (n_thr in 1:length(c(thr, 1))) {
+          #NOTE: FairBreirSs is deprecated now. Should change to SkillScore (according to 
+          #      SpecsVerification's documentation)
+          res[n_exp, n_obs, n_thr] <- SpecsVerification::FairBrierSs(exp[, n_exp, ] > c(thr, 1)[n_thr],
+                                        ens_ref > c(thr, 1)[n_thr],
+                                        obs[, n_obs, 1] > c(thr, 1)[n_thr])['skillscore']
+        }
+      }
+    }
+
+  } else if (type == 'FairEnsembleBS') {
+    #NOTE: The calculation in s2dverification::UltimateBrier is wrong. In the final stage,
+    #      the function calculates like "take(result, 3, 1) - take(result, 3, 2) + take(result, 3, 3)",
+    #      but the 3rd dim of result is 'bins' instead of decomposition. 'FairEnsembleBS' does
+    #      not have decomposition.
+    #      The calculation is fixed here.
+    res <- array(dim = c(nexp = as.numeric(dim(exp)[2]), 
+                         nobs = as.numeric(dim(obs)[2]), 
+                         bin = length(thr) + 1))
+    for (n_exp in 1:dim(exp)[2]) {
+      for (n_obs in 1:dim(obs)[2]) {
+        for (n_thr in 1:length(c(thr, 1))) {
+          fb <- SpecsVerification::FairBrier(ens = exp[, n_exp, ] > c(thr, 1)[n_thr], 
+                                             obs = obs[, n_obs, 1] > c(thr, 1)[n_thr])
+          res[n_exp, n_obs, n_thr] <- mean(fb, na.rm = T)
+        }
+      }
+    }
+#    tmp <- res[, , 1] - res[, , 2] + res[, , 3]
+#    res <- array(tmp, dim = c(nexp = as.numeric(dim(exp)[2]), nobs = as.numeric(dim(obs)[2])))
+
+  } else if (type == 'BS') {
+    comp <- array(dim = c(nexp = as.numeric(dim(exp)[2]),
+                          nobs = as.numeric(dim(obs)[2]),
+                          comp = 3))
+    for (n_exp in 1:dim(exp)[2]) {
+      for (n_obs in 1:dim(obs)[2]) {
+        #NOTE: Parameter 'bins' is default.
+        comp[n_exp, n_obs, ] <- SpecsVerification::BrierDecomp(p = exp[, n_exp], 
+                                                               y = obs[, n_obs])[1, ]
+      }
+    }
+    if (decomposition) {
+      rel <- comp[, , 1]
+      dim(rel) <- c(nexp = as.numeric(dim(exp)[2]), nobs = as.numeric(dim(obs)[2]))
+      res <- comp[, , 2]
+      dim(res) <- c(nexp = as.numeric(dim(exp)[2]), nobs = as.numeric(dim(obs)[2]))
+      unc <- comp[, , 3]
+      dim(unc) <- c(nexp = as.numeric(dim(exp)[2]), nobs = as.numeric(dim(obs)[2]))
+      bs <- rel - res + unc
+      dim(bs) <- c(nexp = as.numeric(dim(exp)[2]), nobs = as.numeric(dim(obs)[2]))
+      res <- list(bs = bs, rel = rel, res = res, unc = unc)
+    } else {
+      bs <- comp[, , 1] - comp[, , 2] + comp[, , 3]
+      dim(bs) <- c(nexp = as.numeric(dim(exp)[2]), nobs = as.numeric(dim(obs)[2]))
+      res <- list(bs = bs)
+    }
+
+  } else if (type == 'FairStartDatesBS') {
+    #NOTE: parameter 'thresholds' is not specified.
+    res <- .BrierScore(exp = exp, obs = obs)
+    if (decomposition) {
+      res <- list(bs = res$bs, rel = res$rel, res = res$res, unc = res$unc)
+    } else {
+      res <- list(bs = res$bs)
+    }
+
+  } else if (type == 'BSS') {
+    #NOTE: parameter 'thresholds' is not specified.
+    res <- .BrierScore(exp = exp, obs = obs)$bss_res
+
+  } else if (type == 'FairStartDatesBSS') {
+    #NOTE: parameter 'thresholds' is not specified.
+    res <- .BrierScore(exp = exp, obs = obs)$bss_gres
+  }
+  
+  return(res)
+
+}
+
diff --git a/man/BrierScore.Rd b/man/BrierScore.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..9271a2adb8fca99ed95483a36b46545e219cd47c
--- /dev/null
+++ b/man/BrierScore.Rd
@@ -0,0 +1,112 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/BrierScore.R
+\name{BrierScore}
+\alias{BrierScore}
+\title{Compute Brier score, its decomposition, and Brier skill score}
+\usage{
+BrierScore(
+  exp,
+  obs,
+  thresholds = seq(0.1, 0.9, 0.1),
+  time_dim = "sdate",
+  dat_dim = NULL,
+  memb_dim = NULL,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A vector or a numeric array with named dimensions. It should be 
+the predicted probabilities which are within the range [0, 1] if memb_dim
+doesn't exist. If it has memb_dim, the value should be 0 or 1, and the 
+predicted probabilities will be computed by ensemble mean. The dimensions 
+must at least have 'time_dim'. 
+range [0, 1].}
+
+\item{obs}{A numeric array with named dimensions of the binary observations 
+(0 or 1). The dimension must be the same as 'exp' except memb_dim, which is
+optional. If it has 'memb_dim', then the length must be 1. The length of 
+'dat_dim' can be different from 'exp' if it has.}
+
+\item{thresholds}{A numeric vector used to bin the forecasts. The default 
+value is \code{seq(0.1, 0.9, 0.1)}, which means that the bins are 
+ \code{[0, 0.1), [0.1, 0.2), ... [0.9, 1]}.}
+
+\item{time_dim}{A character string indicating the name of dimension along
+which Brier score is computed. The default value is 'sdate'.}
+
+\item{dat_dim}{A character string indicating the name of dataset dimension in
+'exp' and 'obs'. The length of this dimension can be different between 
+'exp' and 'obs'. The default value is NULL.}
+
+\item{memb_dim}{A character string of the name of the member dimension in 
+'exp' (and 'obs', optional). The function will do the ensemble mean 
+over this dimension. If there is no member dimension, set NULL. The default
+value is NULL.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+A list that contains:
+\item{$rel}{standard reliability}
+\item{$res}{standard resolution}
+\item{$unc}{standard uncertainty}  
+\item{$bs}{Brier score}
+\item{$bs_check_res}{rel - res + unc}
+\item{$bss_res}{res - rel / unc}
+\item{$gres}{generalized resolution}
+\item{$bs_check_gres}{rel - gres + unc}
+\item{$bss_gres}{gres - rel / unc}
+\item{$rel_bias_corrected}{bias - corrected rel}
+\item{$gres_bias_corrected}{bias - corrected gres}
+\item{$unc_bias_corrected}{bias - corrected unc}
+\item{$bss_bias_corrected}{gres_bias_corrected - rel_bias_corrected / unc_bias_corrected}
+\item{$nk}{number of forecast in each bin}
+\item{$fkbar}{average probability of each bin}
+\item{$okbar}{relative frequency that the observed event occurred}
+The data type and dimensions of the items depend on if the input 'exp' and 
+'obs' are:\cr
+(a) Vectors\cr
+(b) Arrays with 'dat_dim' specified\cr
+(c) Arrays with no 'dat_dim' specified\cr
+Items 'rel', 'res', 'unc', 'bs', 'bs_check_res', 'bss_res', 'gres', 
+'bs_check_gres', 'bss_gres', 'rel_bias_corrected', 'gres_bias_corrected', 
+'unc_bias_corrected', and 'bss_bias_corrected' are (a) a number (b) an array
+with dimensions c(nexp, nobs, all the rest dimensions in 'exp' and 'obs' 
+expect 'time_dim' and 'memb_dim') (c) an array with dimensions of
+'exp' and 'obs' except 'time_dim' and 'memb_dim'\cr
+Items 'nk', 'fkbar', and 'okbar' are (a) a vector of length of bin number 
+determined by 'threshold' (b) an array with dimensions c(nexp, nobs, 
+no. of bins, all the rest dimensions in 'exp' and 'obs' expect 'time_dim' and
+'memb_dim') (c) an array with dimensions c(no. of bin, all the rest dimensions
+in 'exp' and 'obs' expect 'time_dim' and 'memb_dim')
+}
+\description{
+Compute the Brier score (BS) and the components of its standard decompostion
+with the two within-bin components described in Stephenson et al., (2008). It
+also returns the bias-corrected decomposition of the BS (Ferro and Fricker, 
+2012). BSS has the climatology as the reference forecast.
+}
+\examples{
+# Inputs are vectors
+exp <- runif(10)
+obs <- round(exp)
+x <- BrierScore(exp, obs)
+
+# Inputs are arrays
+example(Load)
+clim <- Clim(sampleData$mod, sampleData$obs)
+ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+bins_ano_exp <- ProbBins(ano_exp, thr = c(1/3, 2/3))
+bins_ano_obs <- ProbBins(ano_obs, thr = c(1/3, 2/3))
+res <- BrierScore(bins_ano_exp, MeanDims(bins_ano_obs, 'member'), memb_dim = 'member') 
+
+}
+\references{
+Wilks (2006) Statistical Methods in the Atmospheric Sciences.\cr
+Stephenson et al. (2008). Two extra components in the Brier score decomposition. 
+ Weather and Forecasting, 23: 752-757.\cr
+Ferro and Fricker (2012). A bias-corrected decomposition of the BS. 
+ Quarterly Journal of the Royal Meteorological Society, DOI: 10.1002/qj.1924.
+}
diff --git a/man/UltimateBrier.Rd b/man/UltimateBrier.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..2cad133c23d620eeee9e0a57d7466550da9478cf
--- /dev/null
+++ b/man/UltimateBrier.Rd
@@ -0,0 +1,114 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/UltimateBrier.R
+\name{UltimateBrier}
+\alias{UltimateBrier}
+\title{Compute Brier scores}
+\usage{
+UltimateBrier(
+  exp,
+  obs,
+  dat_dim = "dataset",
+  memb_dim = "member",
+  time_dim = "sdate",
+  quantile = TRUE,
+  thr = c(5/100, 95/100),
+  type = "BS",
+  decomposition = TRUE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A numeric array of forecast anomalies with named dimensions that
+at least include 'dat_dim', 'memb_dim', and 'time_dim'. It can be provided
+by \code{Ano()}.}
+
+\item{obs}{A numeric array of observational reference anomalies with named
+dimensions that at least include 'dat_dim' and 'time_dim'. If it has 
+'memb_dim', the length  must be 1. The dimensions should be consistent with
+'exp' except 'dat_dim' and 'memb_dim'. It can be provided by \code{Ano()}.}
+
+\item{dat_dim}{A character string indicating the name of the dataset 
+dimension in 'exp' and 'obs'. The default value is 'dataset'.}
+
+\item{memb_dim}{A character string indicating the name of the member 
+dimension in 'exp' (and 'obs') for ensemble mean calculation. The default
+value is 'member'.}
+
+\item{time_dim}{A character string indicating the dimension along which to 
+compute the probabilistic scores. The default value is 'sdate'.}
+
+\item{quantile}{A logical value to decide whether a quantile (TRUE) or a 
+threshold (FALSE) is used to estimate the forecast and observed 
+probabilities. If 'type' is 'FairEnsembleBS' or 'FairEnsembleBSS', it must
+be TRUE. The default value is TRUE.}
+
+\item{thr}{A numeric vector to be used in probability calculation (for 'BS', 
+'FairStartDatesBS', 'BSS', and 'FairStartDatesBSS') and binary event 
+judgement (for 'FairEnsembleBS' and 'FairEnsembleBSS'). It is as
+quantiles if 'quantile' is TRUE or as thresholds if 'quantile' is FALSE. 
+The default value is \code{c(0.05, 0.95)} for 'quantile = TRUE'.}
+
+\item{type}{A character string of the desired score type. It can be the 
+ following values:
+\itemize{
+ \item{'BS': Simple Brier Score. Use SpecsVerification::BrierDecomp inside.}
+ \item{'FairEnsembleBS': Corrected Brier Score computed across ensemble 
+   members. Use SpecsVerification::FairBrier inside.}
+ \item{'FairStartDatesBS': Corrected Brier Score computed across starting 
+   dates. Use s2dv:::.BrierScore inside.}
+ \item{'BSS': Simple Brier Skill Score. Use s2dv:::.BrierScore inside.}
+ \item{'FairEnsembleBSS': Corrected Brier Skill Score computed across 
+   ensemble members. Use SpecsVerification::FairBrierSs inside.}
+ \item{'FairStartDatesBSS': Corrected Brier Skill Score computed across 
+   starting dates. Use s2dv:::.BrierScore inside.}
+}
+ The default value is 'BS'.}
+
+\item{decomposition}{A logical value to determine whether the decomposition 
+of the Brier Score should be provided (TRUE) or not (FALSE). It is only 
+used when 'type' is 'BS' or 'FairStartDatesBS'. The default value is TRUE.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+If 'type' is 'BS' or 'FairStartDatesBS' and 'decomposition' is TRUE, the 
+output is a list of 4 arrays (see details below.) In other cases, the output 
+is an array of Brier scores or Brier skill scores. All the arrays have the 
+same dimensions:
+c(nexp, nobs, no. of bins, the rest dimensions of 'exp' except 'time_dim' and
+'memb_dim'). 'nexp' and 'nobs' is the length of dataset dimension in 'exp'
+and 'obs' respectively.\cr
+The list of 4 includes: 
+ \itemize{
+   \item{$bs: Brier Score}
+   \item{$rel: Reliability component}
+   \item{$res: Resolution component}
+   \item{$unc: Uncertainty component}
+ }
+}
+\description{
+Interface to compute probabilistic scores (Brier Score, Brier Skill Score) 
+from the forecast and observational data anomalies. It provides six types
+to choose.
+}
+\examples{
+ \dontshow{
+startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
+sampleData <- s2dv:::.LoadSampleData('tos', c('experiment'),
+                                    c('observation'), startDates,
+                                    leadtimemin = 1,
+                                    leadtimemax = 4,
+                                    output = 'lonlat',
+                                    latmin = 27, latmax = 48,
+                                    lonmin = -12, lonmax = 40)
+ }
+sampleData$mod <- Season(sampleData$mod, monini = 11, moninf = 12, monsup = 2)
+sampleData$obs <- Season(sampleData$obs, monini = 11, moninf = 12, monsup = 2)
+clim <- Clim(sampleData$mod, sampleData$obs)
+exp <- Ano(sampleData$mod, clim$clim_exp)
+obs <- Ano(sampleData$obs, clim$clim_obs)
+bs <- UltimateBrier(exp, obs)
+bss <- UltimateBrier(exp, obs, type = 'BSS')
+
+}
diff --git a/tests/testthat/test-BrierScore.R b/tests/testthat/test-BrierScore.R
new file mode 100644
index 0000000000000000000000000000000000000000..5668a089972352b29db2d553acb7134a3ede2203
--- /dev/null
+++ b/tests/testthat/test-BrierScore.R
@@ -0,0 +1,269 @@
+context("s2dv::BrierScore tests")
+
+##############################################
+# dat1
+set.seed(1)
+exp1 <- array(runif(10), dim = c(dataset = 1, sdate = 5, ftime = 2))
+set.seed(2)
+obs1 <- array(round(runif(10)), dim = c(dataset = 1, sdate = 5, ftime = 2))
+
+# dat2
+set.seed(1)
+exp2 <- runif(10)
+set.seed(2)
+obs2 <- round(runif(10))
+
+# dat3
+set.seed(1)
+exp3 <- array(sample(c(0, 1), 60, replace = T),
+              dim = c(dataset = 2, member = 3, sdate = 5, ftime = 2))
+set.seed(2)
+obs3 <- array(sample(c(0, 1), 10, replace = T),
+              dim = c(dataset = 1, sdate = 5, ftime = 2))
+
+##############################################
+test_that("1. Input checks", {
+  # exp and obs
+  expect_error(
+  BrierScore(exp1, c()),
+  "Parameter 'exp' and 'obs' cannot be NULL."
+  )
+  expect_error(
+  BrierScore(c('b'), obs1),
+  "Parameter 'exp' and 'obs' must be a numeric vector or a numeric array."
+  )
+  expect_error(
+  BrierScore(array(1:10, dim = c(2, 5)), array(1:4, dim = c(2, 2))),
+  "Parameter 'exp' and 'obs' must have dimension names."
+  )
+  expect_error(
+  BrierScore(exp = 1:10, obs = obs1),
+  "Parameter 'exp' must be within \\[0, 1\\] range."
+  )
+  expect_error(
+  BrierScore(exp = array(exp1, dim = dim(exp3)), obs = obs3, memb_dim = 'member'),
+  "Parameter 'exp' must be 0 or 1 if it has memb_dim."
+  )
+  expect_error(
+  BrierScore(exp1, runif(10)),
+  "Parameter 'obs' must be binary events \\(0 or 1\\)."
+  )
+  expect_error( 
+  BrierScore(exp3, obs3),
+  paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
+         "of all the dimensions expect 'dat_dim' and 'memb_dim'.")
+  )
+  # thresholds
+  expect_error(  
+  BrierScore(exp2, obs2, thresholds = TRUE),
+  "Parameter 'thresholds' must be a numeric vector."
+  )
+  expect_error(
+  BrierScore(exp2, obs2, thresholds = seq(0, 1, length.out = 4)),
+  "Parameter 'thresholds' must be between 0 and 1 as the bin-breaks."
+  )
+  # time_dim
+  expect_error(
+  BrierScore(exp2, obs2, time_dim = 1),
+  "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+  BrierScore(exp1, obs1, memb_dim = 'member', time_dim = 'time'),
+  "Parameter 'time_dim' is not found in 'exp' and 'obs' dimension."
+  )
+  # dat_dim
+  expect_error(
+  BrierScore(exp1, obs1, dat_dim = 2),
+  "Parameter 'dat_dim' must be a character string."
+  )
+  expect_error(
+  BrierScore(exp1, obs1, dat_dim = 'dat'),
+  "Parameter 'dat_dim' is not found in 'exp' and 'obs' dimension."
+  )
+  # memb_dim
+  expect_error(
+  BrierScore(exp2, obs2, memb_dim = 2),
+  "Parameter 'memb_dim' must be a character string."
+  )
+  expect_error(
+  BrierScore(exp2, obs2, memb_dim = 'ensemble'),
+  "Parameter 'memb_dim' is not found in 'exp' dimension."
+  )
+  expect_error(
+  BrierScore(exp3, array(1, dim = c(dataset = 1, member = 3, sdate = 5, ftime = 2)), memb_dim = 'member'),
+  "The length of parameter 'memb_dim' in 'obs' must be 1."
+  )
+  # ncores
+  expect_error(
+  BrierScore(exp2, obs2, ncores = 0),
+  "Parameter 'ncores' must be a positive integer."
+  )
+
+})
+
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+expect_equal(
+length(BrierScore(exp1, obs1)),
+16
+)
+expect_equal(
+names(BrierScore(exp1, obs1)),
+c('rel', 'res', 'unc', 'bs', 'bs_check_res', 'bss_res', 'gres',
+  'bs_check_gres', 'bss_gres', 'rel_bias_corrected', 'gres_bias_corrected',
+  'unc_bias_corrected', 'bss_bias_corrected', 'nk', 'fkbar', 'okbar')
+)
+expect_equal(
+dim(BrierScore(exp1, obs1)$rel),
+c(dataset = 1, ftime = 2)
+)
+expect_equal(
+BrierScore(exp1, obs1)$rel[1, ],
+c(0.3086934, 0.3650011),
+tolerance = 0.0001
+)
+expect_equal(
+BrierScore(exp1, obs1)$res[1, ],
+c(0.14, 0.14),
+tolerance = 0.0001
+)
+expect_equal(
+BrierScore(exp1, obs1)$bs[1, ],
+c(0.4218661, 0.4587647),
+tolerance = 0.0001
+)
+expect_equal(
+dim(BrierScore(exp1, obs1)$okbar),
+c(bin = 10, dataset = 1, ftime = 2)
+)
+expect_equal(
+BrierScore(exp1, obs1)$okbar[, 1, 1],
+c(NaN, NaN, 0.5,1.0, NaN, 1.0, NaN, NaN, NaN, 0.0)
+)
+expect_equal(
+BrierScore(exp1, obs1)$fkbar[, 1, 1],
+c(NaN, NaN, 0.2335953, 0.3721239, NaN, 0.5728534, NaN, NaN, NaN, 0.9082078),
+tolerance = 0.0001
+)
+expect_equal(
+BrierScore(exp1, obs1)$nk[, 1, 1],
+c(0, 0, 2, 1, 0, 1, 0, 0, 0, 1)
+)
+
+expect_equal(
+dim(BrierScore(exp1, obs1, dat_dim = 'dataset')$rel),
+c(nexp = 1, nobs = 1, ftime = 2)
+)
+expect_equal(
+dim(BrierScore(exp1, obs1, dat_dim = 'dataset')$nk),
+c(nexp = 1, nobs = 1, bin = 10, ftime = 2)
+)
+expect_equal(
+as.vector(BrierScore(exp1, obs1, dat_dim = 'dataset')$nk),
+as.vector(BrierScore(exp1, obs1)$nk)
+)
+expect_equal(
+as.vector(BrierScore(exp1, obs1, dat_dim = 'dataset')$bs),
+as.vector(BrierScore(exp1, obs1)$bs)
+)
+
+})
+
+##############################################
+test_that("3. Output checks: dat2", {
+expect_equal(
+length(BrierScore(exp2, obs2)),
+16
+)
+expect_equal(
+dim(BrierScore(exp2, obs2))$bss,
+NULL
+)
+expect_equal(
+length(BrierScore(exp2, obs2)$bss_res),
+1
+)
+expect_equal(
+dim(BrierScore(exp2, obs2))$nk,
+NULL
+) 
+expect_equal(
+length(BrierScore(exp2, obs2)$nk),
+10
+)
+expect_equal(
+BrierScore(exp2, obs2)$bs,
+0.4403154,
+tolerance = 0.0001
+)
+expect_equal(
+BrierScore(exp2, obs2)$gres_bias_corrected,
+0.06313199,
+tolerance = 0.0001
+)
+expect_equal(
+BrierScore(exp2, obs2)$bss_bias_corrected,
+-0.6511828,
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(BrierScore(exp2, obs2)$nk),
+c(1, 0, 2, 1, 0, 1, 2, 0, 1, 2)
+)
+
+})
+##############################################
+test_that("4. Output checks: dat3", {
+
+expect_equal(
+length(BrierScore(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')),
+16
+)
+expect_equal(
+names(BrierScore(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')),
+c('rel', 'res', 'unc', 'bs', 'bs_check_res', 'bss_res', 'gres',
+  'bs_check_gres', 'bss_gres', 'rel_bias_corrected', 'gres_bias_corrected',
+  'unc_bias_corrected', 'bss_bias_corrected', 'nk', 'fkbar', 'okbar')
+)
+expect_equal(
+dim(BrierScore(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$rel),
+c(nexp = 2, nobs = 1, ftime = 2)
+)
+expect_equal(
+BrierScore(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$rel[, 1, 2],
+c(0.3555556, 0.2222222),
+tolerance = 0.0001
+)
+expect_equal(
+BrierScore(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$res[1, 1, ],
+c(0.0000000, 0.1066667),
+tolerance = 0.0001
+)
+expect_equal(
+BrierScore(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$bs[2, 1, ],
+c(0.3555556, 0.4222222),
+tolerance = 0.0001
+)
+expect_equal(
+dim(BrierScore(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$okbar),
+c(nexp = 2, nobs = 1, bin = 10, ftime = 2)
+)
+expect_equal(
+BrierScore(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$okbar[, 1, 1, 1],
+c(NaN, 1)
+)
+expect_equal(
+BrierScore(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', thresholds = 1:2/3)$fkbar[2, 1, , 1],
+c(0.0000000, 0.3333333, 0.6666667),
+tolerance = 0.0001
+)
+expect_equal(
+BrierScore(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', thresholds = 1:2/3)$nk[1, 1, , 1],
+c(0, 5, 0)
+)
+
+})
+
+
diff --git a/tests/testthat/test-UltimateBrier.R b/tests/testthat/test-UltimateBrier.R
new file mode 100644
index 0000000000000000000000000000000000000000..654aad020d5d7df65981316ab1918b3cece94c51
--- /dev/null
+++ b/tests/testthat/test-UltimateBrier.R
@@ -0,0 +1,240 @@
+context("s2dv::UltimateBrier tests")
+
+##############################################
+# dat1
+set.seed(1)
+exp1 <- array(rnorm(30), dim = c(dataset = 1, member = 3, sdate = 5, ftime = 2))
+set.seed(2)
+obs1 <- array(round(rnorm(10)), dim = c(dataset = 1, sdate = 5, ftime = 2))
+
+
+##############################################
+test_that("1. Input checks", {
+  # exp and obs
+  expect_error(
+  UltimateBrier(exp1, c()),
+  "Parameter 'exp' and 'obs' cannot be NULL."
+  )
+  expect_error(
+  UltimateBrier(c('b'), obs1),
+  "Parameter 'exp' and 'obs' must be a vector or a numeric array."
+  )
+  expect_error(
+  UltimateBrier(array(1:10, dim = c(2, 5)), array(1:4, dim = c(2, 2))),
+  "Parameter 'exp' and 'obs' must have dimension names."
+  )
+  # dat_dim
+  expect_error(
+  UltimateBrier(exp1, obs1, dat_dim = 2),
+  "Parameter 'dat_dim' must be a character string."
+  )
+  expect_error(
+  UltimateBrier(exp1, obs1, dat_dim = 'dat'),
+  "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  # memb_dim
+  expect_error(
+  UltimateBrier(exp1, obs1, memb_dim = 2),
+  "Parameter 'memb_dim' must be a character string."
+  )
+  expect_error(
+  UltimateBrier(exp1, obs1, memb_dim = 'ensemble'),
+  "Parameter 'memb_dim' is not found in 'exp' dimension."
+  )
+  expect_error(
+  UltimateBrier(exp1, array(1:10, dim = c(dataset = 1, member = 2, sdate = 5, ftime = 2))),
+  "The length of parameter 'memb_dim' in 'obs' must be 1."
+  )
+  # time_dim
+  expect_error(
+  UltimateBrier(exp1, obs1, time_dim = 1),
+  "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+  UltimateBrier(exp1, obs1, memb_dim = 'member', time_dim = 'time'),
+  "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  # exp and obs (2)
+  expect_error(
+  UltimateBrier(exp1, array(1:10, dim = c(dataset = 1, sdate = 6, ftime = 2))),
+  paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
+                "of all the dimensions expect 'dat_dim' and 'memb_dim'.")
+  )
+  expect_error(
+  UltimateBrier(exp1, array(1:10, dim = c(dataset = 1, sdate = 5, time = 2))),
+  paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
+                "of all the dimensions expect 'dat_dim' and 'memb_dim'.")
+  )
+  # quantile
+  expect_error(
+  UltimateBrier(exp1, obs1, quantile = c(0.05, 0.95)),
+  "Parameter 'quantile' must be one logical value."
+  )
+  expect_error(
+  UltimateBrier(exp1, obs1, quantile = FALSE, thr = 1:3, type = 'FairEnsembleBS'),
+  "Parameter 'quantile' must be TRUE if 'type' is 'FairEnsembleBSS' or 'FairEnsembleBS'."
+  )
+  # thr
+  expect_error(
+  UltimateBrier(exp1, obs1, thr = TRUE),
+  "Parameter 'thr' must be a numeric vector."
+  )
+  expect_error(
+  UltimateBrier(exp1, obs1, quantile = TRUE, thr = 1:3),
+  "Parameter 'thr' must be between 0 and 1 when quantile is TRUE."
+  )
+  # type
+  expect_error(
+  UltimateBrier(exp1, obs1, type = 'UltimateBrier'),
+  "Parameter 'type' must be one of 'BS', 'BSS', 'FairEnsembleBS', 'FairEnsembleBSS', 'FairStartDatesBS' or 'FairStartDatesBSS'."
+  )
+  # decomposition
+  expect_error(
+  UltimateBrier(exp1, obs1, decomposition = 1),
+  "Parameter 'decomposition' must be one logical value."
+  )
+  # ncores
+  expect_error(
+  UltimateBrier(exp1, obs1, ncores = 0),
+  "Parameter 'ncores' must be a positive integer."
+  )
+
+
+})
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+# 'BS'
+expect_equal(
+is.list(UltimateBrier(exp1, obs1)),
+TRUE
+)
+expect_equal(
+names(UltimateBrier(exp1, obs1)),
+c('bs', 'rel', 'res', 'unc')
+)
+expect_equal(
+is.list(UltimateBrier(exp1, obs1, decomposition = FALSE)),
+FALSE
+)
+expect_equal(
+dim(UltimateBrier(exp1, obs1, decomposition = FALSE)),
+c(nexp = 1, nobs = 1, bin = 3, ftime = 2)
+)
+expect_equal(
+dim(UltimateBrier(exp1, obs1, decomposition = FALSE, thr = c(0.25, 0.5, 0.75))),
+c(nexp = 1, nobs = 1, bin = 4, ftime = 2)
+)
+expect_equal(
+UltimateBrier(exp1, obs1)$bs,
+UltimateBrier(exp1, obs1, decomposition = FALSE)
+)
+expect_equal(
+as.vector(UltimateBrier(exp1, obs1)$bs),
+c(0.42222222, 0.44444444, 0.02222222, 0.48888889, 0.37777778, 0.02222222),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(UltimateBrier(exp1, obs1)$rel),
+c(0.22222222, 0.31111111, 0.02222222, 0.28888889, 0.24444444, 0.02222222),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(UltimateBrier(exp1, obs1)$res),
+c(0.0400000, 0.1066667, 0.0000000, 0.0400000, 0.1066667, 0.0000000),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(UltimateBrier(exp1, obs1)$unc),
+c(0.24, 0.24, 0.00, 0.24, 0.24, 0.00),
+tolerance = 0.0001
+)
+
+# 'BSS'
+expect_equal(
+dim(UltimateBrier(exp1, obs1, type = 'BSS')),
+c(nexp = 1, nobs = 1, bin = 3, ftime = 2)
+)
+expect_equal(
+as.vector(UltimateBrier(exp1, obs1, type = 'BSS')),
+c(-0.7592593, -0.8518519, -Inf, -1.0370370, -0.5740741, -Inf),
+tolerance = 0.0001
+)
+
+# 'FairStartDatesBS'
+expect_equal(
+is.list(UltimateBrier(exp1, obs1, type = 'FairStartDatesBS')),
+TRUE
+)
+expect_equal(
+names(UltimateBrier(exp1, obs1, type = 'FairStartDatesBS')),
+c('bs', 'rel', 'res', 'unc')
+)
+expect_equal(
+is.list(UltimateBrier(exp1, obs1, decomposition = FALSE, type = 'FairStartDatesBS')),
+FALSE
+)
+expect_equal(
+dim(UltimateBrier(exp1, obs1, decomposition = FALSE, type = 'FairStartDatesBS')),
+c(nexp = 1, nobs = 1, bin = 3, ftime = 2)
+)
+expect_equal(
+UltimateBrier(exp1, obs1, type = 'FairStartDatesBS')$bs,
+UltimateBrier(exp1, obs1, decomposition = FALSE, type = 'FairStartDatesBS')
+)
+expect_equal(
+as.vector(UltimateBrier(exp1, obs1, type = 'FairStartDatesBS')$bs),
+c(0.42222222, 0.44444444, 0.02222222, 0.48888889, 0.37777778, 0.02222222),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(UltimateBrier(exp1, obs1, type = 'FairStartDatesBS')$rel),
+c(0.22222222, 0.31111111, 0.02222222, 0.28888889, 0.24444444, 0.02222222),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(UltimateBrier(exp1, obs1, type = 'FairStartDatesBS')$res),
+c(0.0400000, 0.1066667, 0.0000000, 0.0400000, 0.1066667, 0.0000000),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(UltimateBrier(exp1, obs1, type = 'FairStartDatesBS')$unc),
+c(0.24, 0.24, 0.00, 0.24, 0.24, 0.00),
+tolerance = 0.0001
+)
+
+# 'FairStartDatesBSS'
+expect_equal(
+dim(UltimateBrier(exp1, obs1, type = 'FairStartDatesBSS')),
+c(nexp = 1, nobs = 1, bin = 3, ftime = 2)
+)
+expect_equal(
+as.vector(UltimateBrier(exp1, obs1, type = 'FairStartDatesBSS')),
+c(-0.7592593, -0.8518519, -Inf, -1.0370370, -0.5740741, -Inf),
+tolerance = 0.0001
+)
+# 'FairEnsembleBS'
+expect_equal(
+dim(UltimateBrier(exp1, obs1, type = 'FairEnsembleBS')),
+c(nexp = 1, nobs = 1, bin = 3, ftime = 2)
+)
+expect_equal(
+as.vector(UltimateBrier(exp1, obs1, type = 'FairEnsembleBS')),
+c(0.1333333, 0.2000000, 0.2000000, 0.1333333, 0.4000000, 0.2000000),
+tolerance = 0.0001
+)
+# 'FairEnsembleBSS'
+expect_equal(
+dim(UltimateBrier(exp1, obs1, type = 'FairEnsembleBSS')),
+c(nexp = 1, nobs = 1, bin = 3, ftime = 2)
+)
+expect_equal(
+as.vector(UltimateBrier(exp1, obs1, type = 'FairEnsembleBSS')),
+c(-0.1111111, -0.6666667, -0.6666667, 0.2592593, -1.2222222, -0.6666667),
+tolerance = 0.0001
+)
+
+})
+
+