Commit bfc11bc4 authored by aho's avatar aho
Browse files

Update to master

parents 19ef45fa 742ca4eb
......@@ -36,7 +36,9 @@ Imports:
plyr,
ncdf4,
NbClust,
multiApply (>= 2.1.1)
multiApply (>= 2.1.1),
SpecsVerification (>= 0.5.0),
easyNCDF
Suggests:
easyVerification,
testthat
......
......@@ -4,6 +4,8 @@ export(ACC)
export(AMV)
export(AnimateMap)
export(Ano)
export(Ano_CrossValid)
export(BrierScore)
export(CDORemap)
export(Clim)
export(Cluster)
......@@ -21,8 +23,12 @@ export(ConfigRemoveEntry)
export(ConfigShowDefinitions)
export(ConfigShowSimilarEntries)
export(ConfigShowTable)
export(Consist_Trend)
export(Corr)
export(EOF)
export(Eno)
export(EuroAtlanticTC)
export(Filter)
export(GMST)
export(GSAT)
export(Histo2Hindcast)
......@@ -30,31 +36,41 @@ export(InsertDim)
export(LeapYear)
export(Load)
export(MeanDims)
export(NAO)
export(Persistence)
export(PlotACC)
export(PlotAno)
export(PlotBoxWhisker)
export(PlotClim)
export(PlotEquiMap)
export(PlotLayout)
export(PlotMatrix)
export(PlotSection)
export(PlotStereoMap)
export(PlotVsLTime)
export(ProbBins)
export(ProjectField)
export(REOF)
export(RMS)
export(RMSSS)
export(RandomWalkTest)
export(RatioRMS)
export(RatioSDRMS)
export(Regression)
export(Reorder)
export(SPOD)
export(Season)
export(Smoothing)
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)
......@@ -70,6 +86,7 @@ importFrom(ClimProjDiags,Subset)
importFrom(ClimProjDiags,WeightedMean)
importFrom(abind,abind)
importFrom(abind,adrop)
importFrom(easyNCDF,ArrayToNc)
importFrom(grDevices,bmp)
importFrom(grDevices,col2rgb)
importFrom(grDevices,colorRampPalette)
......@@ -105,5 +122,7 @@ importFrom(stats,quantile)
importFrom(stats,rnorm)
importFrom(stats,sd)
importFrom(stats,setNames)
importFrom(stats,spectrum)
importFrom(stats,ts)
importFrom(stats,varimax)
importFrom(stats,window)
#'Compute anomalies in cross-validation mode
#'
#'Compute the anomalies from the arrays of the experimental and observational
#'data output by subtracting the climatologies computed with a leave-one-out
#'cross validation technique and a per-pair method (Garcia-Serrano and
#'Doblas-Reyes, CD, 2012).
#'Per-pair climatology means that only the start dates covered by the
#'whole experiments/observational datasets will be used. In other words, the
#'startdates which do not all have values along 'dat_dim' dimension of both
#'the 'exp' and 'obs' are excluded when computing the climatologies.
#'
#'@param exp A named numeric array of experimental data, with at least
#' dimensions 'time_dim' and 'dat_dim'.
#'@param obs A named numeric array of observational data, same dimensions as
#' parameter 'exp' except along 'dat_dim'.
#'@param time_dim A character string indicating the name of the time dimension.
#' The default value is 'sdate'.
#'@param dat_dim A character vector indicating the name of the dataset and
#' member dimensions. When calculating the climatology, if data at one
#' startdate (i.e., 'time_dim') is not complete along 'dat_dim', this startdate
#' along 'dat_dim' will be discarded. The default value is
#' "c('dataset', 'member')".
#'@param memb_dim A character string indicating the name of the member
#' dimension. Only used when parameter 'memb' is FALSE. It must be one element
#' in 'dat_dim'. The default value is 'member'.
#'@param memb A logical value indicating whether to subtract the climatology
#' based on the individual members (TRUE) or the ensemble mean over all
#' members (FALSE) when calculating the anomalies. 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
#'A list of 2:
#'\item{$exp}{
#' A numeric array with the same dimensions as 'exp'. The dimension order may
#' change.
#'}
#'\item{$obs}{
#' A numeric array with the same dimensions as 'obs'.The dimension order may
#' change.
#'}
#'
#'@examples
#'# Load sample data as in Load() example:
#'example(Load)
#'anomalies <- Ano_CrossValid(sampleData$mod, sampleData$obs)
#'\dontrun{
#'PlotAno(anomalies$exp, anomalies$obs, startDates,
#' toptitle = paste('anomalies'), ytitle = c('K', 'K', 'K'),
#' legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano_crossvalid.eps')
#'}
#'@import multiApply
#'@importFrom ClimProjDiags Subset
#'@export
Ano_CrossValid <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'),
memb_dim = 'member', memb = 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 numeric array.")
}
if (is.null(dim(exp)) | is.null(dim(obs))) {
stop(paste0("Parameter 'exp' and 'obs' must have at least dimensions ",
"time_dim and dat_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(!all(names(dim(exp)) %in% names(dim(obs))) |
!all(names(dim(obs)) %in% names(dim(exp)))) {
stop("Parameter 'exp' and 'obs' must have same dimension name.")
}
## 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.")
}
## dat_dim
if (!is.character(dat_dim)) {
stop("Parameter 'dat_dim' must be a character vector.")
}
if (!all(dat_dim %in% names(dim(exp))) | !all(dat_dim %in% names(dim(obs)))) {
stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.")
}
## memb
if (!is.logical(memb) | length(memb) > 1) {
stop("Parameter 'memb' must be one logical value.")
}
## memb_dim
if (!memb) {
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)) | !memb_dim %in% names(dim(obs))) {
stop("Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension.")
}
if (!memb_dim %in% dat_dim) {
stop("Parameter 'memb_dim' must be one element in parameter 'dat_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.")
}
}
## exp and obs (2)
name_exp <- sort(names(dim(exp)))
name_obs <- sort(names(dim(obs)))
for (i in 1:length(dat_dim)) {
name_exp <- name_exp[-which(name_exp == dat_dim[i])]
name_obs <- name_obs[-which(name_obs == dat_dim[i])]
}
if(!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
stop(paste0("Parameter 'exp' and 'obs' must have the same length of ",
"all dimensions except 'dat_dim'."))
}
###############################
# Sort dimension
name_exp <- names(dim(exp))
name_obs <- names(dim(obs))
order_obs <- match(name_exp, name_obs)
if (any(order_obs != sort(order_obs))) {
obs <- Reorder(obs, order_obs)
}
#-----------------------------------
# Per-paired method: If any sdate along dat_dim is NA, turn all sdate points along dat_dim into NA.
pos <- rep(0, length(dat_dim)) # dat_dim: [dataset, member]
for (i in 1:length(dat_dim)) {
pos[i] <- which(names(dim(obs)) == dat_dim[i])
}
outrows_exp <- MeanDims(exp, pos, na.rm = FALSE) +
MeanDims(obs, pos, na.rm = FALSE)
outrows_obs <- outrows_exp
for (i in 1:length(pos)) {
outrows_exp <- InsertDim(outrows_exp, pos[i], dim(exp)[pos[i]])
outrows_obs <- InsertDim(outrows_obs, pos[i], dim(obs)[pos[i]])
}
exp_for_clim <- exp
obs_for_clim <- obs
exp_for_clim[which(is.na(outrows_exp))] <- NA
obs_for_clim[which(is.na(outrows_obs))] <- NA
#-----------------------------------
res <- Apply(list(exp, obs, exp_for_clim, obs_for_clim),
target_dims = c(time_dim, dat_dim),
fun = .Ano_CrossValid,
memb_dim = memb_dim, memb = memb,
ncores = ncores)
return(res)
}
.Ano_CrossValid <- function(exp, obs, exp_for_clim, obs_for_clim,
memb_dim = 'member', memb = TRUE, ncores = NULL) {
# exp: [sdate, dat_dim, memb_dim]
# obs: [sdate, dat_dim, memb_dim]
ano_exp_list <- vector('list', length = dim(exp)[1]) #length: [sdate]
ano_obs_list <- vector('list', length = dim(obs)[1])
for (tt in 1:dim(exp)[1]) { #[sdate]
# calculate clim
exp_sub <- ClimProjDiags::Subset(exp_for_clim, 1, c(1:dim(exp)[1])[-tt])
obs_sub <- ClimProjDiags::Subset(obs_for_clim, 1, c(1:dim(obs)[1])[-tt])
clim_exp <- apply(exp_sub, c(1:length(dim(exp)))[-1], mean, na.rm = TRUE) # average out time_dim -> [dat, memb]
clim_obs <- apply(obs_sub, c(1:length(dim(obs)))[-1], mean, na.rm = TRUE)
# ensemble mean
if (!memb) {
if (is.null(dim(clim_exp)) | length(dim(clim_exp)) == 1) { #dim: [member]
clim_exp <- mean(clim_exp, na.rm = TRUE) # a number
clim_obs <- mean(clim_obs, na.rm = TRUE)
} else {
pos <- which(names(dim(clim_exp)) == memb_dim)
pos <- c(1:length(dim(clim_exp)))[-pos]
dim_name <- names(dim(clim_exp))
dim_exp_ori <- dim(clim_exp)
dim_obs_ori <- dim(clim_obs)
clim_exp <- apply(clim_exp, pos, mean, na.rm = TRUE)
clim_obs <- apply(clim_obs, pos, mean, na.rm = TRUE)
if (is.null(names(dim(as.array(clim_exp))))) {
clim_exp <- as.array(clim_exp)
clim_obs <- as.array(clim_obs)
names(dim(clim_exp)) <- dim_name[pos]
names(dim(clim_obs)) <- dim_name[pos]
}
# Expand it back
clim_exp_tmp <- array(clim_exp, dim = c(dim_exp_ori[pos], dim_exp_ori[-pos]))
clim_obs_tmp <- array(clim_obs, dim = c(dim_obs_ori[pos], dim_obs_ori[-pos]))
# Reorder it back to dim(clim_exp)
tmp <- match(dim_exp_ori, dim(clim_exp_tmp))
if (any(tmp != sort(tmp))) {
clim_exp <- Reorder(clim_exp_tmp, tmp)
clim_obs <- Reorder(clim_obs_tmp, tmp)
} else {
clim_exp <- clim_exp_tmp
clim_obs <- clim_obs_tmp
}
}
}
# calculate ano
ano_exp_list[[tt]] <- ClimProjDiags::Subset(exp, 1, tt, drop = 'selected') - clim_exp
ano_obs_list[[tt]] <- ClimProjDiags::Subset(obs, 1, tt, drop = 'selected') - clim_obs
}
ano_exp <- array(unlist(ano_exp_list), dim = c(dim(exp)[-1], dim(exp)[1]))
ano_exp <- Reorder(ano_exp, c(length(dim(exp)), 1:(length(dim(exp)) - 1)))
ano_obs <- array(unlist(ano_obs_list), dim = c(dim(obs)[-1], dim(obs)[1]))
ano_obs <- Reorder(ano_obs, c(length(dim(obs)), 1:(length(dim(obs)) - 1)))
return(list(exp = ano_exp, obs = ano_obs))
}
#'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()---------