diff --git a/R/ProbBins.R b/R/ProbBins.R index c6e4fec19ff4832ac1d8b799df00a9fbe4dc24da..4aee08171eca6db67d659befe728ba42740338f9 100644 --- a/R/ProbBins.R +++ b/R/ProbBins.R @@ -1,9 +1,8 @@ -ProbBins <- function(ano, fcyr, thr, quantile=TRUE, posdates = 3, - posdim = 2, compPeriod= "Full period") { +ProbBins <- function(ano, fcyr = 'all', thr, quantile = TRUE, posdates = 3, + posdim = 2, compPeriod = "Full period") { # define dimensions #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ nbdim <- length(dim(ano)) - nfcyr<-length(fcyr) if (nbdim < 7){ ano <- Enlarge(ano, 7) @@ -13,78 +12,101 @@ ProbBins <- function(ano, fcyr, thr, quantile=TRUE, posdates = 3, posdim <- setdiff(posdim, posdates) nbpos <- length(posdim) #permute dimensions in ano - ano <- aperm(ano, c(posdates, posdim, setdiff(seq(1,7,1), c(posdates, posdim)))) + if (posdates != 1 || posdim != 2) { + dimnames_backup <- names(dim(ano)) + perm <- c(posdates, posdim, (1:7)[-c(posdates, posdim)]) + ano <- aperm(ano, perm) + names(dim(ano)) <- dimnames_backup[perm] + posdates <- 1 + posdim <- 2 + } dimano <- dim(ano) - nsdates=dimano[1] + nsdates <- dimano[1] #calculate the number of elements on posdim dimension in ano - nmemb=1 + nmemb <- 1 if (nbpos > 0){ for (idim in 2:(nbpos+1)){ - nmemb=nmemb*dimano[idim] + nmemb <- nmemb*dimano[idim] } } - - # separate forecast and hindcast - #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - fore <- array(ano[fcyr, , , , , , ], dim = c(nfcyr, - dimano[2:7])) - # the members and startdates are combined in one dimension - sample_fore <- array(fore, dim=c(nfcyr*nmemb, dimano[(nbpos+2):7])) - - if(compPeriod=="Full period"){ - hind <- ano - sample <- array(hind, dim=c(nsdates*nmemb, dimano[(nbpos+2):7])) - } - - if (compPeriod=="Without fcyr"){ - hind <- array(ano[-fcyr, , , , , , ], dim = c(nsdates-nfcyr, - dimano[2:7])) - sample <- array(hind, dim=c((nsdates-nfcyr)*nmemb, dimano[(nbpos+2):7])) - } - - #quantiles for each grid point and experiment - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - if (quantile==TRUE){ - qum <- apply(sample, seq(2,7-nbpos,1), FUN=quantile,probs=thr,na.rm=TRUE,names=FALSE,type=8) - }else{ - qum<-array(thr,dim=c(length(thr), dimano[(nbpos+2):7])) + + if (fcyr == 'all') { + fcyr <- 1:nsdates } + nfcyr <- length(fcyr) - # This function assign the values to a category which is limited by the thresholds - # It provides binary information - - counts <- function (dat, nbthr){ - thr <- dat[1:nbthr] - data <- dat[nbthr+1:(length(dat)-nbthr)] - prob <- array(NA, dim=c(nbthr+1, length(dat)-nbthr)) - prob[1,]=1*(data <= thr[1]) - if(nbthr!=1){ - for (ithr in 2:(nbthr)){ - prob[ithr,]=1*((data > thr[ithr - 1]) & (data <= thr[ithr])) + if (compPeriod == "Cross-validation") { + result <- NULL + for (iyr in fcyr) { + if (is.null(result)) { + result <- ProbBins(ano, iyr, thr, quantile, + posdates, posdim, "Without fcyr") + } else { + dimnames_backup <- names(dim(result)) + result <- abind(result, ProbBins(ano, iyr, thr, quantile, + posdates, posdim, "Without fcyr"), + along = 2) + names(dim(result)) <- dimnames_backup } } - prob[nbthr+1,]=1*(data > thr[nbthr]) - return(prob) - } - - # The thresholds and anomalies are combined to use apply - data <- abind(qum, sample_fore, along = 1) - - # PBF:Probabilistic bins of a forecast. - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # This array contains zeros and ones that indicate the category where your forecast is. - - PBF <- array(apply(data, seq(2,7-nbpos,1), FUN=counts, nbthr=length(thr)), - dim=c(length(thr)+1, nfcyr, nmemb, dimano[(nbpos+2):nbdim])) + return(result) + } else if (compPeriod %in% c("Full period", "Without fcyr")) { + # separate forecast and hindcast + #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + fore <- array(ano[fcyr, , , , , , ], dim = c(nfcyr, + dimano[2:7])) + # the members and startdates are combined in one dimension + sample_fore <- array(fore, dim=c(nfcyr*nmemb, dimano[(nbpos+2):7])) + + if(compPeriod == "Full period") { + hind <- ano + sample <- array(hind, dim=c(nsdates*nmemb, dimano[(nbpos+2):7])) + } else if (compPeriod == "Without fcyr") { + hind <- array(ano[-fcyr, , , , , , ], dim = c(nsdates-nfcyr, + dimano[2:7])) + sample <- array(hind, dim=c((nsdates-nfcyr)*nmemb, dimano[(nbpos+2):7])) + } - return(PBF) + #quantiles for each grid point and experiment + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + if (quantile==TRUE){ + qum <- apply(sample, seq(2,7-nbpos,1), FUN=quantile,probs=thr,na.rm=TRUE,names=FALSE,type=8) + }else{ + qum<-array(thr,dim=c(length(thr), dimano[(nbpos+2):7])) + } + # This function assign the values to a category which is limited by the thresholds + # It provides binary information + + counts <- function (dat, nbthr){ + thr <- dat[1:nbthr] + data <- dat[nbthr+1:(length(dat)-nbthr)] + prob <- array(NA, dim=c(nbthr+1, length(dat)-nbthr)) + prob[1,]=1*(data <= thr[1]) + if(nbthr!=1){ + for (ithr in 2:(nbthr)){ + prob[ithr,]=1*((data > thr[ithr - 1]) & (data <= thr[ithr])) + } + } + prob[nbthr+1,]=1*(data > thr[nbthr]) + return(prob) + } + + # The thresholds and anomalies are combined to use apply + data <- abind(qum, sample_fore, along = 1) + + # PBF:Probabilistic bins of a forecast. + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # This array contains zeros and ones that indicate the category where your forecast is. - if (compPeriod == "cross-validation"){ - for (iyr in 1:fcyr){ - ProbBins(ano,fcyr=iyr,thr=thr,posdates=posdates,posdim=posdim) - } + PBF <- array(apply(data, seq(2,7-nbpos,1), FUN=counts, nbthr=length(thr)), + dim=c(length(thr)+1, nfcyr, nmemb, dimano[(nbpos+2):nbdim])) + + names(dim(PBF)) <- c('bin', 'sdate', 'member', names(dim(ano))[(nbpos+2):nbdim]) + return(PBF) + } else { + stop("Parameter 'compPeriod' must be one of 'Full period', 'Without fcyr' or 'Cross-validation'.") } } diff --git a/man/ProbBins.Rd b/man/ProbBins.Rd index 96ec05767efa099e6678380311beff046b6b704b..a1c5ac5005d9393b061c029cd655d447435975d2 100644 --- a/man/ProbBins.Rd +++ b/man/ProbBins.Rd @@ -7,7 +7,7 @@ Computes probabilistic information of a forecast relative to a threshold or a qu Compute probabilistic bins of a set of forecast years ('fcyr') relative to the forecast climatology over the whole period of anomalies, optionally excluding the selected forecast years ('fcyr') or the forecast year for which the probabilistic bins are being computed (see 'compPeriod'). } \usage{ -ProbBins(ano, fcyr, thr, quantile = TRUE, posdates = 3, posdim = 2, +ProbBins(ano, fcyr = 'all', thr, quantile = TRUE, posdates = 3, posdim = 2, compPeriod = "Full period") } \arguments{ @@ -16,8 +16,8 @@ Array of anomalies from Ano().\cr Must be of dimension (nexp/nobs, nmemb, nsdates, nleadtime, nlat, nlon) } \item{fcyr}{ -Indices of the forecast years of the anomalies of which to compute the probabilistic bins.\cr -Ex: c(1:5), c(1, 4) or 4. +Indices of the forecast years of the anomalies which to compute the probabilistic bins for, or 'all' to compute the bins for all the years.\cr +Ex: c(1:5), c(1, 4), 4 or 'all'. } \item{thr}{ Values used as thresholds to bin the anomalies. @@ -33,12 +33,12 @@ Position of the dimension in \code{ano} that corresponds to the start dates (def Position of the dimension in \code{ano} which will be combined with 'posdates' to compute the quantiles (default = 2, ensemble members). } \item{compPeriod}{ -Three options: "Full period"/"Without fcyr"/"cross-validation" (The probabilities are computed with the terciles based on ano/ano with all 'fcyr's removed/cross-validation). The default is "Full period". +Three options: "Full period"/"Without fcyr"/"Cross-validation" (The probabilities are computed with the terciles based on ano/ano with all 'fcyr's removed/cross-validation). The default is "Full period". } } \value{ Matrix with probabilistic information and dimensions:\cr -c(length('thr'+1), nfcyr, nmemb/nparam, nmod/nexp/nobs, nltime, nlat, nlon)\cr +c(length('thr') + 1, length(fcyr), nmemb/nparam, nmod/nexp/nobs, nltime, nlat, nlon)\cr The values along the first dimension take values 0 or 1 depending on which of the 'thr'+1 cathegories the forecast/observation at the corresponding grid point, time step, member and starting date belongs to. } \examples{ @@ -74,6 +74,7 @@ PB <- ProbBins(ano_exp, fcyr = 3, thr = c(1/3, 2/3), quantile = TRUE, posdates = \author{ History:\cr 1.0 - 2013 (F.Lienert) - Original code\cr -2.0 - 2014-03 (N. Gonzalez and V.Torralba, \email{veronica.torralba@ic3.cat}) - Debugging +2.0 - 2014-03 (N. Gonzalez and V. Torralba, \email{veronica.torralba at bsc.es}) - Debugging +2.1 - 2017-02 (V. Torralba and N. Manubens, \email{veronica.torralba at bsc.es}) - Fix bug with cross-validation } \keyword{datagen}