diff --git a/R/BrierScore.R b/R/BrierScore.R index c64434ed5ea97810f72b03b90ef18535943462a4..e1b01bcf7afc64e6e163a20c0131797dba2fb144 100644 --- a/R/BrierScore.R +++ b/R/BrierScore.R @@ -14,14 +14,14 @@ BrierScore <- function(obs, pred, thresholds = seq(0, 1, 0.1)) { bins[[i]] <- list(which(pred >= thresholds[i] & pred < thresholds[i + 1])) } } - + fkbar <- okbar <- nk <- array(0, dim = nbins) for (i in 1:nbins) { nk[i] <- length(bins[[i]][[1]]) fkbar[i] <- sum(pred[bins[[i]][[1]]]) / nk[i] okbar[i] <- sum(obs[bins[[i]][[1]]]) / nk[i] } - + obar <- sum(obs) / length(obs) relsum <- ressum <- term1 <- term2 <- 0 for (i in 1:nbins) { @@ -43,7 +43,7 @@ BrierScore <- function(obs, pred, thresholds = seq(0, 1, 0.1)) { 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 # @@ -63,12 +63,100 @@ BrierScore <- function(obs, pred, thresholds = seq(0, 1, 0.1)) { } 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") #} - + invisible(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 = bins, pred = pred, obs = obs)) } } + + +.BrierScore <- function(exp, obs, thresholds = seq(0, 1, 0.1)) { + if (max(exp) > 1 || min(exp) < 0) { + stop("Parameter 'exp' contains predictions outside [0,1] range. Are you certain this is a probability forecast?") + } else if (max(obs) != 1 && min(obs) != 0) { + .message("Binary events in 'obs' must be either 0 or 1. Are you certain this is a binary event?") + } else { + nbins <- length(thresholds) - 1 # Number of bins + n <- dim(exp)[1] # Number of observations + ens_mean <- rowMeans(exp, na.rm = TRUE) + n.ens <- seq(1, dim(exp)[2], 1) # Number of ensemble members + bins <- as.list(paste("bin", 1:nbins, sep = "")) + for (i in 1:nbins) { + if (i == nbins) { + bins[[i]] <- list(which(ens_mean >= thresholds[i] & ens_mean <= thresholds[i + 1])) + } else { + bins[[i]] <- list(which(ens_mean >= thresholds[i] & ens_mean < thresholds[i + 1])) + } + } + + fkbar <- okbar <- nk <- array(0, dim = nbins) + for (i in 1:nbins) { + nk[i] <- length(bins[[i]][[1]]) + fkbar[i] <- sum(ens_mean[bins[[i]][[1]]]) / nk[i] + okbar[i] <- sum(obs[bins[[i]][[1]]]) / nk[i] + } + + fkbar[fkbar == Inf] <- 0 + okbar[is.nan(okbar)] <- 0 + obar <- sum(obs) / length(obs) + relsum <- ressum <- relsum1 <- ressum1 <- term1 <- term1a <- term2 <- term2a <- 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 + (ens_mean[bins[[i]][[1]][j]] - fkbar[i]) ^ 2 + term2 <- term2 + (ens_mean[bins[[i]][[1]][j]] - fkbar[i]) * (obs[bins[[i]][[1]][j]] - okbar[i]) + } + } + } + } + rel <- relsum / n + res <- ressum / n + unc <- obar * (1 - obar) + #bs <- apply(ens, MARGIN = 2, FUN = function(x) sum((x - obs)^2) / n) + bs <- sum((rowMeans(exp, na.rm = T) - 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") + #} + + invisible(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)) +} diff --git a/R/Corr.R b/R/Corr.R index e23895b7696e8d19b2e66c4b8b6f6e7093bbaa37..10858c29494a38681c96381946144e228fd46101 100644 --- a/R/Corr.R +++ b/R/Corr.R @@ -5,8 +5,8 @@ Corr <- function(var_exp, var_obs, posloop = 1, poscor = 2, compROW = NULL, # Remove data along compROW dim if there is at least one NA between limits # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # - if (is.null(compROW) == FALSE) { - if (is.null(limits) == TRUE) { + if (!is.null(compROW)) { + if (is.null(limits)) { limits <- c(1, dim(var_obs)[compROW]) } outrows <- (is.na(Mean1Dim(var_obs, compROW, narm = FALSE, limits))) @@ -42,7 +42,7 @@ Corr <- function(var_exp, var_obs, posloop = 1, poscor = 2, compROW = NULL, var_obs <- aperm(var_obs, posaperm) dimsaperm <- dim(var_exp) # - + # Check the siglev arguments: # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (siglev > 1 || siglev < 0) { @@ -74,9 +74,9 @@ Corr <- function(var_exp, var_obs, posloop = 1, poscor = 2, compROW = NULL, for (j9 in 1:dimsaperm[9]) { for (j10 in 1:dimsaperm[10]) { tmp1 <- var_exp[jexp, , j3, j4, j5, j6, j7, j8, j9, - j10] + j10] tmp2 <- var_obs[jobs, , j3, j4, j5, j6, j7, j8, j9, - j10] + j10] if (any(!is.na(tmp1)) && sum(!is.na(tmp2)) > 2) { toto <- cor(tmp1, tmp2, use = "pairwise.complete.obs", method = method) CORR[jexp, jobs, dim_val, j3, j4, j5, j6, j7, j8, j9, j10] <- toto @@ -92,17 +92,17 @@ Corr <- function(var_exp, var_obs, posloop = 1, poscor = 2, compROW = NULL, #t <- qt(0.95, eno - 2) t <- qt(siglev, eno - 2) CORR[jexp, jobs, dim_pval, j3, j4, j5, j6, j7, j8, j9, - j10] <- sqrt((t * t) / ((t * t) + eno - 2)) + j10] <- sqrt((t * t) / ((t * t) + eno - 2)) } if (conf) { CORR[jexp, jobs, 1, j3, j4, j5, j6, j7, j8, j9, - j10] <- tanh(atanh(toto) + qnorm(conf_high) / sqrt( - #j10] <- tanh(atanh(toto) + qnorm(0.975) / sqrt( - eno - 3)) + j10] <- tanh(atanh(toto) + qnorm(conf_high) / sqrt( + #j10] <- tanh(atanh(toto) + qnorm(0.975) / sqrt( + eno - 3)) CORR[jexp, jobs, 3, j3, j4, j5, j6, j7, j8, j9, - j10] <- tanh(atanh(toto) + qnorm(conf_low) / sqrt( - #j10] <- tanh(atanh(toto) + qnorm(0.025) / sqrt( - eno - 3)) + j10] <- tanh(atanh(toto) + qnorm(conf_low) / sqrt( + #j10] <- tanh(atanh(toto) + qnorm(0.025) / sqrt( + eno - 3)) } } } @@ -123,3 +123,44 @@ Corr <- function(var_exp, var_obs, posloop = 1, poscor = 2, compROW = NULL, # CORR } + +.Corr <- function(exp, obs, siglev = 0.95, method = 'pearson', + conf = TRUE, pval = TRUE) { + + # Check 'method' + if (!(method %in% c("kendall", "spearman", "pearson"))) { + stop("Parameter 'method' must be one of 'kendall', 'spearman' or 'pearson'.") + # Check 'siglev' + if (siglev > 1 || siglev < 0) { + stop("Parameter 'siglev' must be higher than O and lower than 1.") + } + } + + p_val <- NULL + conflow <- NULL + confhigh <- NULL + ens_mean <- rowMeans(exp) + CORR <- cor(obs, ens_mean, use = "pairwise.complete.obs", method = method) + if (pval || conf) { + if (method == "kendall" || method == "spearman") { + eno <- Eno(rank(obs), 1) + } else if (method == "pearson") { + eno <- Eno(obs, 1) + } + } + if (pval && (method == "pearson")) { + t <- CORR * sqrt((eno - 2) / (1 - (CORR ^ 2))) + p_val <- 1 - pt(t, eno - 2) + } + if (conf && (method == "pearson")) { + conf_low <- (1 - siglev) / 2 + conf_high <- 1 - conf_low + conf_int <- c(tanh(atanh(CORR) + qnorm(conf_low) / sqrt(eno - 3)), + tanh(atanh(CORR) + qnorm(conf_high) / sqrt(eno - 3))) + conf_int <- conf_int[!is.na(CORR)] + conflow <- conf_int[1] + confhigh <- conf_int[2] + } + # Output + invisible(result <- list(corr = CORR, p_val = p_val, conf_low = conflow, conf_high = confhigh)) +} diff --git a/R/PlotEquiMap.R b/R/PlotEquiMap.R index dea625df467a81c26fb406e476e3b3036eaf306e..d8bf07e255299a3ed4f0aebaef9302f6a5036c60 100644 --- a/R/PlotEquiMap.R +++ b/R/PlotEquiMap.R @@ -135,7 +135,6 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL, if (!.IsColor(colNA)) { stop("Parameter 'colNA' must be a valid colour identifier.") } - # Check: brks, cols, subsampleg, bar_limits, color_fun, bar_extra_labels, draw_bar_ticks # draw_separators, triangle_ends_scale, label_scale, units, units_scale, # bar_label_digits diff --git a/R/ProbBins.R b/R/ProbBins.R index c6e4fec19ff4832ac1d8b799df00a9fbe4dc24da..fef21045aa5f71b807852442f8316e85c0aa314e 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,103 @@ 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])) - } - # 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 (length(fcyr) == 1) { + if (fcyr == 'all') { + fcyr <- 1:nsdates } - 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. + nfcyr <- length(fcyr) - 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])) + 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 + } + } + 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/R/RMS.R b/R/RMS.R index 1a63b48fdbf4eac911db3c8f47d2ade99bf8a5a9..b45f8cab0f8813b8f3499348bac2b5d02a6b9a08 100644 --- a/R/RMS.R +++ b/R/RMS.R @@ -56,9 +56,9 @@ RMS <- function(var_exp, var_obs, posloop = 1, posRMS = 2, compROW = NULL, for (jobs in 1:nobs) { dif <- array(dim = dimsaperm[-1]) dif[, , , , , , , , ] <- permvarexp[jexp, , , , , , , , - , ] - permvarobs[jobs, , , , , , , , , ] + , ] - permvarobs[jobs, , , , , , , , , ] enlrms[jexp, jobs, dim_rms, , , , , , , , ] <- Mean1Dim(dif ** 2, 1, - narm = TRUE) ** 0.5 + narm = TRUE) ** 0.5 if (conf) { eno <- Eno(dif, 1) for (j3 in 1:dimsaperm[3]){ @@ -70,17 +70,17 @@ RMS <- function(var_exp, var_obs, posloop = 1, posRMS = 2, compROW = NULL, for (j9 in 1:dimsaperm[9]){ for (j10 in 1:dimsaperm[10]){ ndat <- length(sort(dif[, j3, j4, j5, j6, j7, j8, j9, - j10])) + j10])) enlrms[jexp, jobs, 1, j3, j4, j5, j6, j7, j8, j9, j10] <- (eno[j3, j4, j5, j6, j7, j8, j9, - j10] * enlrms[jexp, jobs, 2, j3, j4, j5, j6, j7, - j8, j9, j10] ** 2 / qchisq(conf_high, eno[j3, j4, j5, - j6, j7, j8, j9, j10] - 1)) ** 0.5 + j10] * enlrms[jexp, jobs, 2, j3, j4, j5, j6, j7, + j8, j9, j10] ** 2 / qchisq(conf_high, eno[j3, j4, j5, + j6, j7, j8, j9, j10] - 1)) ** 0.5 enlrms[jexp, jobs, 3, j3, j4, j5, j6, j7, j8, j9, j10] <- (eno[j3, j4, j5, j6, j7, j8, j9, - j10] * enlrms[jexp, jobs, 2, j3, j4, j5, j6, j7, - j8, j9, j10] ** 2 / qchisq(conf_low, eno[j3, j4, j5, - j6, j7, j8, j9, j10] - 1)) ** 0.5 + j10] * enlrms[jexp, jobs, 2, j3, j4, j5, j6, j7, + j8, j9, j10] ** 2 / qchisq(conf_low, eno[j3, j4, j5, + j6, j7, j8, j9, j10] - 1)) ** 0.5 } } } @@ -92,7 +92,7 @@ RMS <- function(var_exp, var_obs, posloop = 1, posRMS = 2, compROW = NULL, } } } - + dim(enlrms) <- c(nexp, nobs, nvals, dimsvar[-c(posloop, posRMS)]) # # Output @@ -100,3 +100,30 @@ RMS <- function(var_exp, var_obs, posloop = 1, posRMS = 2, compROW = NULL, # enlrms } + +.RMS <- function(exp, obs, siglev = 0.95, conf = TRUE) { + # + # RMS & its confidence interval computation + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # + if (conf) { + conf_low <- (1 - siglev) / 2 + conf_high <- 1 - conf_low + } + dif <- rowMeans(exp) - obs + enlrms <- mean(dif ** 2, na.rm = TRUE) ** 0.5 + if (conf) { + eno <- Eno(dif, 1) + ndat <- length(sort(dif)) + conf_int <- c((eno * enlrms ** 2 / qchisq(conf_high, eno - 1)) ** 0.5, + (eno * enlrms ** 2 / qchisq(conf_low, eno - 1)) ** 0.5) + names(conf_int) <- c("conf_low","conf_high") + } else { + conf_int <- c() + names(conf_int) <- c() + } + + results <- c(enlrms, conf_int) + names(results) <- c("rms", names(conf_int)) + return(results) +} diff --git a/R/RMSSS.R b/R/RMSSS.R index 585a398b8ba0b770c75e1280dff847181b0921db..42622477f39ab291b4c179ec9799dacc1e94787f 100644 --- a/R/RMSSS.R +++ b/R/RMSSS.R @@ -33,13 +33,13 @@ RMSSS <- function(var_exp, var_obs, posloop = 1, posRMS = 2, pval = TRUE) { nvals <- 1 } enlRMSSS <- array(dim = c(nexp, nobs, nvals, dimsaperm[3:10])) - + for (jexp in 1:nexp) { for (jobs in 1:nobs) { dif1 <- array(dim = dimsaperm[-1]) dif2 <- array(dim = dimsaperm[-1]) dif1[, , , , , , , , ] <- permvarexp[jexp, , , , , , , , - , ] - permvarobs[jobs, , , , , , , , , ] + , ] - permvarobs[jobs, , , , , , , , , ] dif2[, , , , , , , , ] <- permvarobs[jobs, , , , , , , , , ] rms1 <- Mean1Dim(dif1 ** 2, 1, narm = TRUE) ** 0.5 rms2 <- Mean1Dim(dif2 ** 2, 1, narm = TRUE) ** 0.5 @@ -63,7 +63,7 @@ RMSSS <- function(var_exp, var_obs, posloop = 1, posRMS = 2, pval = TRUE) { if (is.na(l1) == FALSE & is.na(l2) == FALSE & l1 > 2 & l2 > 2) { enlRMSSS[jexp, jobs, 2, j3, j4, j5, j6, j7, j8, j9, j10] <- 1 - pf(F[j3, j4, j5, j6, j7, j8, j9, - j10], l1 - 1, l2 - 1) + j10], l1 - 1, l2 - 1) } else { enlRMSSS[jexp, jobs, 1, j3, j4, j5, j6, j7, j8, j9, j10] <- NA @@ -79,7 +79,7 @@ RMSSS <- function(var_exp, var_obs, posloop = 1, posRMS = 2, pval = TRUE) { } } } - + dim(enlRMSSS) <- c(nexp, nobs, nvals, dimsvar[-c(posloop, posRMS)]) # # Output @@ -87,3 +87,28 @@ RMSSS <- function(var_exp, var_obs, posloop = 1, posRMS = 2, pval = TRUE) { # enlRMSSS } + +.RMSSS <- function(exp, obs, pval = TRUE) { + dif2 <- obs + dif1 <- rowMeans(exp) - obs + rms1 <- mean(dif1 ** 2, na.rm = TRUE) ** 0.5 + rms2 <- mean(dif2 ** 2, na.rm = TRUE) ** 0.5 + rms2[abs(rms2) <= abs(rms2) / 1000] <- abs(rms2) / 1000 + rmsss <- c(1 - (rms1 / rms2)) + if (pval == TRUE) { + eno1 <- Eno(dif1, 1) + eno2 <- Eno(dif2, 1) + F.stat <- (eno2 * (rms2) ** 2 / (eno2 - 1)) / (eno1 * (rms1) ** 2 / (eno1 - 1)) + if (is.na(eno1) == FALSE && is.na(eno2) == FALSE && eno1 > 2 && eno2 > 2) { + p_val <- 1 - pf(F.stat, eno1 - 1, eno2 - 1) + } + } else { + p_val <- NA + } + + # + # Output + # ~~~~~~~~ + # + list(rmsss = rmsss, p_val = p_val) +} diff --git a/R/RatioRMS.R b/R/RatioRMS.R index a97ec02f8fb59be249cb5b8cef54f9a8a4aa9cba..564f9658d5f2a1566372d8b9d98ca56fe13ac7b5 100644 --- a/R/RatioRMS.R +++ b/R/RatioRMS.R @@ -55,8 +55,8 @@ RatioRMS <- function(var_exp1, var_exp2, var_obs, posRMS = 1, pval = TRUE) { l2 <- eno2[j2, j3, j4, j5, j6, j7, j8, j9, j10] if (!is.na(l1) && !is.na(l2) && l1 > 2 && l2 > 2) { enlratiorms[2, j2, j3, j4, j5, j6, j7, j8, j9, - j10] <- (1 - pf(F[j2, j3, j4, j5, j6, j7, j8, j9, j10], - l1 - 1, l2 - 1)) * 2 + j10] <- (1 - pf(F[j2, j3, j4, j5, j6, j7, j8, j9, j10], + l1 - 1, l2 - 1)) * 2 } else { enlratiorms[1, j2, j3, j4, j5, j6, j7, j8, j9, j10] <- NA } @@ -84,3 +84,29 @@ RatioRMS <- function(var_exp1, var_exp2, var_obs, posRMS = 1, pval = TRUE) { # enlratiorms } + +.RatioRMS <- function(exp, exp_ref, obs, pval = TRUE) { + # + # RMS ratio and its pvalue computation + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # + dif1 <- rowMeans(exp, na.rm = TRUE) - obs + dif2 <- rowMeans(exp_ref, na.rm = TRUE) - obs + rms1 <- mean(dif1 ** 2, na.rm = TRUE) ** 0.5 + rms2 <- mean(dif2 ** 2, na.rm = TRUE) ** 0.5 + rms2[which(abs(rms2) <= (max(abs(rms2), na.rm = TRUE) / 1000))] <- max( + abs(rms2), na.rm = TRUE) / 1000 + enlratiorms <- (rms1 / rms2) + if (pval) { + eno1 <- Eno(dif1, 1) + eno2 <- Eno(dif2, 1) + F <- (eno1 * (rms1) ** 2 / (eno1 - 1)) / (eno2 * (rms2) ** 2 / (eno2 - 1)) + F[which(F < 1)] <- 1 / F[which(F < 1)] + if (!is.na(eno1) && !is.na(eno2) && eno1 > 2 && eno2 > 2) { + p_val <- (1 - pf(F, eno1 - 1, eno2 - 1)) * 2 + } + } + + # Output + list(ratiorms = enlratiorms, p_val = p_val) +} diff --git a/R/RatioSDRMS.R b/R/RatioSDRMS.R index d76a832ba49426ed46568ec98dca871ccecdb052..2bce490dda043a35b3dfc5fd7a9a6308f63873e0 100644 --- a/R/RatioSDRMS.R +++ b/R/RatioSDRMS.R @@ -40,7 +40,7 @@ RatioSDRMS <- function(var_exp, var_obs, pval = TRUE) { } dif <- var_exp - InsertDim(ensmeanexp, 2, dimexp[2]) std <- apply(array(dif, dim = c(dimexp[1], dimexp[2] * dimexp[3], - dimrms[3:6])), c(1, 3, 4, 5, 6), sd, na.rm = TRUE) + dimrms[3:6])), c(1, 3, 4, 5, 6), sd, na.rm = TRUE) enosd <- apply(Eno(dif, 3), c(1, 3, 4, 5, 6), sum, na.rm = TRUE) rms <- array(dim = dimrms) enlratiormssd <- array(dim = dimenlratiormssd) @@ -50,7 +50,7 @@ RatioSDRMS <- function(var_exp, var_obs, pval = TRUE) { rms[jexp,jobs, , , , ] <- Mean1Dim(dif ** 2, 1, narm = TRUE) ** 0.5 enorms <- array(Eno(dif, 1), dim = dimrms[3:6]) enlratiormssd[jexp, jobs, 1, , , , ] <- std[jexp, , , , ] / rms[jexp, - jobs, , , , ] + jobs, , , , ] if (pval) { for (jltime in 1:dimrms[3]) { for (jlev in 1:dimrms[4]) { @@ -59,10 +59,10 @@ RatioSDRMS <- function(var_exp, var_obs, pval = TRUE) { l1 <- enosd[jexp, jltime, jlev, jlat, jlon] l2 <- enorms[jltime, jlev, jlat, jlon] F <- (enosd[jexp, jltime, jlev, jlat, jlon] * (std[jexp, jltime, - jlev, jlat, jlon]) ** 2 / (enosd[jexp, jltime, jlev, jlat, - jlon] - 1)) / (enorms[jltime, jlev, jlat, jlon] * (rms[jexp, - jobs, jltime, jlev, jlat, jlon]) ** 2 / (enorms[jltime, - jlev, jlat, jlon] - 1)) + jlev, jlat, jlon]) ** 2 / (enosd[jexp, jltime, jlev, jlat, + jlon] - 1)) / (enorms[jltime, jlev, jlat, jlon] * (rms[jexp, + jobs, jltime, jlev, jlat, jlon]) ** 2 / (enorms[jltime, + jlev, jlat, jlon] - 1)) if (!is.na(F) && !is.na(l1) && !is.na(l2) && l1 > 2 && l2 > 2) { enlratiormssd[jexp, jobs, 2, jltime, jlev, jlat, jlon] <- 1 - pf(F, l1 - 1, l2 - 1) @@ -83,3 +83,36 @@ RatioSDRMS <- function(var_exp, var_obs, pval = TRUE) { # enlratiormssd } + +.RatioSDRMS <- function(exp, obs, pval = TRUE) { + + # + # Ratio RMSE / SD and its significance level + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # + ensmean <- rowMeans(exp) + enosd <- Eno(ensmean, 1) + dif <- exp - ensmean + std <- sd(dif, na.rm = TRUE) + dif <- ensmean - obs + rms <- mean(dif ** 2, na.rm = TRUE) ** 0.5 + enorms <- Eno(dif, 1) + enlratiormssd <- std / rms + p_val <- 0 + if (pval) { + l1 <- enosd + l2 <- enorms + F <- (enosd * std ** 2 / (enosd - 1)) / (enorms * (rms) ** 2 / (enorms - 1)) + if (!is.na(F) && !is.na(l1) && !is.na(l2) && l1 > 2 && l2 > 2) { + p_val <- 1 - pf(F, l1 - 1, l2 - 1) + } + else { + p_val <- NA + } + } + # + # Output + # ~~~~~~~~ + # + list(ratio = enlratiormssd, p_val = p_val) +} diff --git a/R/Trend.R b/R/Trend.R index 799841601a807be7b1450f05920763d78b800a35..e035df2ab27a6dd71302661db6c8875bd10f92d4 100644 --- a/R/Trend.R +++ b/R/Trend.R @@ -52,7 +52,7 @@ Trend <- function(var, posTR = 2, interval = 1, siglev = 0.95, conf = TRUE) { j10] <- confint(lm.out, level = siglev)[2, 1:2] } enldetrend[is.na(tmp) == FALSE, j2, j3, j4, j5, j6, j7, j8, - j9, j10] <- tmp[is.na(tmp) == FALSE] - lm.out$fitted.values + j9, j10] <- tmp[is.na(tmp) == FALSE] - lm.out$fitted.values } } } @@ -78,3 +78,24 @@ Trend <- function(var, posTR = 2, interval = 1, siglev = 0.95, conf = TRUE) { # invisible(list(trend = enltrend, detrended = enldetrend)) } + +.Trend <- function(exp, interval = 1, siglev = 0.95, conf = TRUE) { + + ensmean <- rowMeans(exp, na.rm = TRUE) + + if (any(!is.na(ensmean))) { + mon <- seq(ensmean) * interval + lm.out <- lm(ensmean ~ mon, na.action = na.omit) + trend <- c(lm.out$coefficients[2], lm.out$coefficients[1]) + if (conf) { + conf.int <- confint(lm.out, level = siglev)[2, 1:2] + } + detrend <- ensmean[is.na(ensmean) == FALSE] - lm.out$fitted.values + } + + # + # Outputs + # ~~~~~~~~~ + # + invisible(list(trend = trend, conf.int = conf.int, detrended = detrend)) +} diff --git a/man/ACC.Rd b/man/ACC.Rd index ccb69e71d65881c345c37c6d02ad6d7a8c76dabd..d4d66aec9b832e2a387e415e5fccd96cf65e1863 100644 --- a/man/ACC.Rd +++ b/man/ACC.Rd @@ -1,13 +1,12 @@ \name{ACC} \alias{ACC} \title{ -Computes Anomaly Correlation Coefficient (Spatial Correlation) +Computes Anomaly Correlation Coefficient } \description{ -Matrix var_exp & var_obs should have dimensions (nexp/nobs, nsdates, nltimes, nlat, nlon) or (nexp/nobs, nsdates, nmember, nltimes, nlat, nlon).\cr -ACC computes the Anomaly Correlation Coefficient for the ensemble mean of each jexp in 1:nexp and each jobs in 1:nobs which gives nexp x nobs ACC for each startdate and each leadtime.\cr -A domain can be selected by providing the list of longitudes/latitudes (lon/lat) of the grid together with the corner of the domain:\cr - lonlatbox = c(lonmin, lonmax, latmin, latmax) +Calculates the Anomaly Correlation Coefficient for the ensemble mean of each model and the corresponding references for each startdate and each leadtime.\cr +The domain of interest can be specified by providing the list of longitudes/latitudes (lon/lat) of the grid together with the corners of the domain:\cr +lonlatbox = c(lonmin, lonmax, latmin, latmax) } \usage{ ACC(var_exp, var_obs, lon = NULL, lat = NULL, lonlatbox = NULL, @@ -15,13 +14,13 @@ ACC(var_exp, var_obs, lon = NULL, lat = NULL, lonlatbox = NULL, } \arguments{ \item{var_exp}{ -Matrix of experimental anomalies with dimensions:\cr +Array of experimental anomalies with dimensions:\cr c(nexp, nsdates, nltimes, nlat, nlon)\cr or\cr c(nexp, nsdates, nmembers, nltimes, nlat, nlon)\cr } \item{var_obs}{ -Matrix of observational anomalies, same dimensions as var_exp except along the first dimension and the second if it corresponds to the member dimension. +Array of observational anomalies, same dimensions as var_exp except along the first dimension and the second if it corresponds to the member dimension. } \item{lon}{ Array of longitudes of the var_exp/var_obs grids, optional. @@ -41,22 +40,23 @@ TRUE/FALSE: confidence intervals and significance level provided or not. To guarantee the statistical robustness of the result, make sure that your experiments/oservations/startdates/leadtimes always have the same number of members. } \item{siglev}{ -Desired confidence level of the computed confidence intervals. +The confidence level for the computed confidence intervals. } } \value{ \item{ACC}{ -If conf set as TRUE, Matrix with dimensions:\cr +If \code{conf = TRUE}, array with dimensions:\cr c(nexp, nobs, nsdates, nleadtimes, 4) \cr The fifth dimension of length 4 corresponds to the lower limit of the \code{siglev}\% confidence interval, the ACC, the upper limit of the \code{siglev}\% confidence interval and the \code{siglev}\% significance level.\cr -If conf set as FALSE, Anomaly Correlation Coefficient with dimensions:\cr +If \code{conf = FALSE}, the array of the Anomaly Correlation Coefficient has dimensions:\cr c(nexp, nobs, nsdates, nleadtimes). } \item{MACC}{ -Mean Anomaly Correlation Coefficient with dimensions:\cr +The array of the Mean Anomaly Correlation Coefficient with dimensions:\cr c(nexp, nobs, nleadtimes) } } + \examples{ # See ?Load for explanations on the first part of this example. \dontrun{ @@ -92,14 +92,17 @@ ano_obs <- Ano(sampleData$obs, clim$clim_obs) acc <- ACC(Mean1Dim(ano_exp, 2), Mean1Dim(ano_obs, 2)) PlotACC(acc$ACC, startDates) } +\references{ +Joliffe and Stephenson (2012). Forecast Verification: A Practitioner's Guide in Atmospheric Science. Wiley-Blackwell. +} \author{ History:\cr -0.1 - 2013-08 (V. Guemas, \email{virginie.guemas at ic3.cat}) - Original code\cr -1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at ic3.cat}) - Formatting to CRAN\cr -1.1 - 2013-09 (C. Prodhomme, \email{chloe.prodhomme at ic3.cat}) - optimization\cr -1.2 - 2014-08 (V. Guemas, \email{virginie.guemas at ic3.cat}) - Bug-fixes : handling of NA & selection of domain + Simplification of code\cr -1.3.0 - 2014-08 (V. Guemas, \email{virginie.guemas at ic3.cat}) - Boostrapping over members\cr -1.3.1 - 2014-09 (C. Prodhomme, chloe.prodhomme at ic3.cat) - Add comments and minor style changes\cr -1.3.2 - 2015-02 (N. Manubens, nicolau.manubens at ic3.cat) - Fixed ACC documentation and examples\cr +0.1 - 2013-08 (V. Guemas, \email{virginie.guemas at bsc.es}) - Original code\cr +1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at bsc.es}) - Formatting to CRAN\cr +1.1 - 2013-09 (C. Prodhomme, \email{chloe.prodhomme at bsc.es}) - optimization\cr +1.2 - 2014-08 (V. Guemas, \email{virginie.guemas at bsc.es}) - Bug-fixes : handling of NA & selection of domain + Simplification of code\cr +1.3.0 - 2014-08 (V. Guemas, \email{virginie.guemas at bsc.es}) - Boostrapping over members\cr +1.3.1 - 2014-09 (C. Prodhomme, chloe.prodhomme at bsc.es) - Add comments and minor style changes\cr +1.3.2 - 2015-02 (N. Manubens, nicolau.manubens at bsc.es) - Fixed ACC documentation and examples\cr } \keyword{datagen} diff --git a/man/Ano_CrossValid.Rd b/man/Ano_CrossValid.Rd index 0261e863cb1746463503df9e20951741e173eb63..7e488ca2e8a7ce923dadadb8957e5f70f18f754c 100644 --- a/man/Ano_CrossValid.Rd +++ b/man/Ano_CrossValid.Rd @@ -2,7 +2,7 @@ \alias{Ano_CrossValid} \title{Computes Anomalies In Cross-Validation Mode} \description{ -This function computes anomalies from experimental and observational matrices output from \code{load()} by subtracting the climatologies computed in a cross-validation mode and with a per-pair method. +Computes the anomalies from the arrays of the experimental and observational data output from \code{load()} by subtracting the climatologies computed with a cross-validation technique and a per-pair method. } \usage{ Ano_CrossValid(var_exp, var_obs, memb = TRUE) diff --git a/man/BrierScore.Rd b/man/BrierScore.Rd index 43843523c2d45a29c2f9e479cd02aef0a1446509..e126729aadedd98349e231190b5aa5a9aac53c99 100644 --- a/man/BrierScore.Rd +++ b/man/BrierScore.Rd @@ -1,16 +1,18 @@ \name{BrierScore} \alias{BrierScore} +\alias{.BrierScore} \title{ Compute Brier Score And Its Decomposition And Brier Skill Score } \description{ -Returns the values of the BS and its standard decomposition as well as the addition of the two winthin-bin extra components (Stephenson et al., 2008). It also solves the bias-corrected decomposition of the BS (Ferro and Fricker, 2012). BSS having the climatology as the reference forecast. -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.\cr +Computes the Brier score (BS) and the components of its standard decomposition as well 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 having the climatology as the reference forecast.\cr +\cr +.BrierScore provides the same functionality, but taking a matrix of ensemble members (exp) as input. } \usage{ BrierScore(obs, pred, thresholds = seq(0, 1, 0.1)) + +.BrierScore(exp, obs, thresholds = seq(0, 1, 0.1)) } \arguments{ \item{obs}{ @@ -22,38 +24,74 @@ Vector of probablistic predictions with values in the range [0,1] \item{thresholds}{ Values used to bin the forecasts. By default the bins are {[0,0.1), [0.1, 0.2), ... [0.9, 1]} } + \item{exp}{ +Matrix of predictions with values in the range [0,1] for the .BrierScore function + } } \value{ -$rel: standard reliability\cr -$res: standard resolution\cr -$unc: standard uncertainty\cr -$bs: Brier score\cr -$bs_check_res: rel-res+unc\cr -$bss_res: res-rel/unc\cr -$gres: generalized resolution\cr -$bs_check_gres: rel-gres+unc\cr -$bss_gres: gres-rel/unc\cr -$rel_bias_corrected: bias-corrected rel\cr -$gres_bias_corrected: bias-corrected gres\cr -$unc_bias_corrected: bias-corrected unc\cr -$bss_bias_corrected: gres_bias_corrected-rel_bias_corrected/unc_bias_corrected\cr -$nk: number of forecast in each bin\cr -$fkbar: average probability of each bin\cr -$okbar: relative frequency that the observed event occurred\cr -$bins: bins used\cr -$pred: values with which the forecasts are verified\cr -$obs: probability forecasts of the event\cr +Both BrierScore and .Brier score provide the same outputs: + \itemize{ + \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} + \item{$bins}{bins used} + \item{$pred}{values with which the forecasts are verified} + \item{$obs}{probability forecasts of the event} + } } \examples{ +# Minimalist examples with BrierScore a <- runif(10) b <- round(a) x <- BrierScore(b, a) x$bs - x$bs_check_res x$bs - x$bs_check_gres x$rel_bias_corrected - x$gres_bias_corrected + x$unc_bias_corrected + \dontrun{ +a <- runif(10) +b <- cbind(round(a),round(a)) # matrix containing 2 identical ensemble members... +x2 <- BrierScore(a, b) + } + +# Example of BrierScore using UltimateBrier +# See ?UltimateBrier for more information +example(Load) +clim <- Clim(sampleData$mod, sampleData$obs) +ano_exp <- Ano(sampleData$mod, clim$clim_exp) +ano_obs <- Ano(sampleData$obs, clim$clim_obs) +bs <- UltimateBrier(ano_exp, ano_obs, thr = c(1/3, 2/3)) + + \dontrun{ +# Example of .BrierScore with veriApply +require(easyVerification) +BrierScore2 <- s2dverification:::.BrierScore +bins_ano_exp <- ProbBins(ano_exp, thr = c(1/3, 2/3), posdates = 3, posdim = 2) +bins_ano_obs <- ProbBins(ano_obs, thr = c(1/3, 2/3), posdates = 3, posdim = 2) +bs2 <- veriApply("BrierScore2", bins_ano_exp, Mean1Dim(bins_ano_obs, 3), + tdim = 2, ensdim = 3) + } +} +\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.\cr } \author{ History:\cr -0.1 - 2012-04 (L. Rodrigues, \email{lrodrigues@ic3.cat}) - Original code\cr +0.1 - 2012-04 (L. Rodrigues, \email{lrodrigues at ic3.cat}) - Original code\cr +0.2 - 2017-02 (A. Hunter, \email{alasdair.hunter at bsc.es}) - Adapted to veriApply() } \keyword{datagen} diff --git a/man/Clim.Rd b/man/Clim.Rd index 03a3982bc41b67ae33946369e3b50b5e3dee6354..a4225653bd63a7ab384c5f95d24beb03607bcb20 100644 --- a/man/Clim.Rd +++ b/man/Clim.Rd @@ -1,6 +1,6 @@ \name{Clim} \alias{Clim} -\title{Computes Per-pair/Kharin/Fuckar Climatologies} +\title{Computes Bias Corrected Climatologies} \description{ This function computes only per-pair climatologies from the experimental and observational matrices output from \code{Load()}. To compute plain climatologies from only experimental or observational data from \code{Load()}, the following code can be used: diff --git a/man/Cluster.Rd b/man/Cluster.Rd index 1f798635ae96d1fbcd59752cf9ee0031ca62df7c..164477312a677c2bc20d9117da3cff557d4e0aea 100644 --- a/man/Cluster.Rd +++ b/man/Cluster.Rd @@ -1,7 +1,7 @@ \name{Cluster} \alias{Cluster} \title{ -The K-means cluster analysis. +K-means Clustering. } \description{ This function computes cluster centers and their time series of occurrences, @@ -12,9 +12,7 @@ to time. Specifically, it partitions the array along time axis in K groups or clusters in which each space vector/array belongs to (i.e., is a member of) the cluster with the nearest center or centroid. This function relies on the NbClust -package (Charrad et al., 2014 JSS). For more information about the K-means see -Chapter 15 Cluster Analysis in Wilks, 2011, Statistical Methods in the -Atmospheric Sciences, 3rd ed., Elsevire, pp 676. +package (Charrad et al., 2014 JSS). } \usage{ Cluster(var, weights, nclusters = NULL, index = 'sdindex', posdates = 1) @@ -115,6 +113,9 @@ res2 <- Cluster(var = a1, weights = array(1, dim = dim(a1)[2])) print(res2$cluster) print(res2$centers) } +\references{ +Wilks, 2011, Statistical Methods in the Atmospheric Sciences, 3rd ed., Elsevire, pp 676. +} \author{ History: 1.0 # 2014-10 (N.S. Fuckar, neven.fuckar@bsc.es) # Original code diff --git a/man/Consist_Trend.Rd b/man/Consist_Trend.Rd index 66deedf85c2ea03f25455af38c7cb1a9264ecc23..6f964bf83bf3350c9ecb4f0952371ac55867eb28 100644 --- a/man/Consist_Trend.Rd +++ b/man/Consist_Trend.Rd @@ -4,9 +4,9 @@ Computes Trends Using Only Model Data For Which Observations Are Available } \description{ -Computes trends by least square fitting together with the associated error interval for both the observational and model data.\cr +Computes the trend coefficients for a time series by least square fitting, together with the associated error interval for both the observational and model data.\cr Provides also the detrended observational and modeled data.\cr -The trend is computed along the second dimension, expected to be the start date dimension (the user is supposed to perform an ensemble averaging operation with \code{Mean1Dim()} prior to using \code{Consist_trend()}). +By default, the trend is computed along the second dimension of the input array, which is expected to be the start date dimension. For arrays containing multiple model members, the user will first have to calculate the ensemble average using \code{Mean1Dim()} or elsewhise (see the example). } \usage{ Consist_Trend(var_exp, var_obs, interval = 1) @@ -29,7 +29,7 @@ Number of months/years between 2 start dates. Default = 1. The trends will be pr } \value{ \item{$trend}{ -Trends of model and observational data with dimensions:\cr +Trend coefficients of model and observational data with dimensions:\cr c(nmod/nexp + nobs, 3, nltime) up to\cr c(nmod/nexp + nobs, 3, nltime, nlevel, nlat, nlon)\cr The length 3 dimension corresponds to the lower limit of the 95\% confidence interval, the slope of the trends and the upper limit of the 95\% confidence interval. diff --git a/man/Corr.Rd b/man/Corr.Rd index 1c5641e0568e4fddcae239172f8ae9222a303ff8..ae4b9f0005733d329afb75d1482ae67b5c1b15e5 100644 --- a/man/Corr.Rd +++ b/man/Corr.Rd @@ -1,88 +1,121 @@ \name{Corr} \alias{Corr} +\alias{.Corr} \title{ -Computes Correlation Skill Measure (Temporal Correlation Along Start Dates) +Computes the correlation coefficient between an array of forecasts and their corresponding observations. } \description{ -Matrix var_exp & var_obs should have the same dimensions except along posloop dimension where the length can be different, with the number of experiments/models for var_exp (nexp) and the number of obserational datasets for var_obs (nobs).\cr -Corr computes the correlation skill of each jexp in 1:nexp against each jobs in 1:nobs which gives nexp x nobs correlation skill measures for each other grid point of the matrix (each latitude/longitude/level/leadtime). \cr -The correlations are computed along the poscor dimension which should correspond to the startdate dimension. If compROW is given, the correlations are computed only if rows along the compROW dimension are complete between limits[1] and limits[2], that mean with no NA between limits[1] and limits[2]. This option can be activated if the user wishes to account only for the forecasts for which observations are available at all leadtimes. \cr -Default: limits[1] = 1 and limits[2] = length(compROW dimension).\cr -The confidence interval is computed by a Fisher transformation.\cr -The significance level relies on a one-sided student-T distribution.\cr -We can modifiy the treshold of the test modifying siglev (default value=0.95). +Calculates the correlation coefficient (Pearson, Kendall or Spearman) for an array of forecasts and observations. The input should be an array with dimensions c(no. of datasets, no. of start dates, no. of forecast times, no. of lons, no. of lats.), where the longitude and latitude dimensions are optional. The correlations are computed along the poscor dimension which should correspond to the startdate dimension. If compROW is given, the correlations are computed only if rows along the compROW dimension are complete between limits[1] and limits[2], i.e. there are no NAs between limits[1] and limits[2]. This option can be activated if the user wishes to account only for the forecasts for which observations are available at all leadtimes. \cr Default: limits[1] = 1 and limits[2] = length(compROW dimension).\cr The confidence interval is computed by a Fisher transformation.\cr The significance level relies on a one-sided student-T distribution.\cr We can modifiy the treshold of the test modifying siglev (default value=0.95). \cr +\cr +.Corr calculates the correlation between the ensemble mean and the observations, using an N by M matrix (exp) of forecasts and a vector of observations (obs) as input. } \usage{ -Corr(var_exp, var_obs, posloop = 1, poscor = 2, compROW = NULL, - limits = NULL, siglev = 0.95, method = 'pearson', - conf = TRUE, pval = TRUE) +Corr(var_exp, var_obs, posloop = 1, poscor = 2, compROW = NULL, + limits = NULL, siglev = 0.95, method = 'pearson', + conf = TRUE, pval = TRUE) + +.Corr(exp, obs, siglev = 0.95, method = 'pearson', + conf = TRUE, pval = TRUE) } \arguments{ \item{var_exp}{ -Matrix of experimental data. +Array of experimental data. } \item{var_obs}{ -Matrix of observational data, same dimensions as var_exp except along posloop dimension, where the length can be nobs instead of nexp. +Array of observational data, same dimensions as var_exp except along posloop dimension, where the length can be nobs instead of nexp. } \item{posloop}{ -Dimension nobs and nexp. +Dimension nobs and nexp. } \item{poscor}{ -Dimension along which correlation are to be computed (the dimension of the start dates). +Dimension along which correlation are to be computed (the dimension of the start dates). } \item{compROW}{ -Data taken into account only if (compROW)th row is complete.\cr -Default = NULL. +Data taken into account only if (compROW)th row is complete.\cr Default = NULL. } \item{limits}{ -Complete between limits[1] & limits[2]. Default = NULL. +Complete between limits[1] & limits[2]. Default = NULL. } \item{siglev}{ -Significance level according. Default = 0.95. +Significance level. Default = 0.95. } \item{method}{ Type of correlation: 'pearson', 'spearman' or 'kendall'. Default='pearson' } \item{conf}{ -Whether to compute confidence intervals (TRUE; default) or not (FALSE). +Whether to compute confidence intervals (default = 'TRUE') or not (FALSE). } \item{pval}{ -Whether to compute statistical significance p-value (TRUE; default) or not (FALSE). +Whether to compute statistical significance p-value (default = 'TRUE') or not (FALSE). + } + \item{exp}{ +N by M matrix of N forecasts from M ensemble members. + } + \item{obs}{ +Vector of the corresponding observations of length N. } } \value{ -Matrix with dimensions :\cr -c(# of datasets alogn posloop in var_exp, # of datasets along posloop in var_obs, 4, all other dimensions of var_exp & var_obs except poscor).\cr -The third dimension, of length 4 maximum, contains to the lower limit of the 95\% confidence interval, the correlation, the upper limit of the 95\% confidence interval and the 95\% significance level given by a one-sided T-test. If the p-value is disabled via \code{pval = FALSE}, this dimension will be of length 3. If the confidence intervals are disabled via \code{conf = FALSE}, this dimension will be of length 2. If both are disabled, this will be of length 2. +Corr: Array with dimensions :\cr c(# of datasets along posloop in var_exp, # of datasets along posloop in var_obs, 4, all other dimensions of var_exp & var_obs except poscor).\cr +The third dimension, of length 4 maximum, contains to the lower limit of the 95\% confidence interval, the correlation, the upper limit of the 95\% confidence interval and the 95\% significance level given by a one-sided T-test. If the p-value is disabled via \code{pval = FALSE}, this dimension will be of length 3. If the confidence intervals are disabled via \code{conf = FALSE}, this dimension will be of length 2. If both are disabled, this will be of length 2. \cr +\cr +.Corr: + \itemize{ + \item{$corr}{ +The correlation statistic. + } + \item{$p_val}{ +Corresponds to the p values for the \code{siglev}\% (only present if \code{pval = TRUE}) for the correlation. + } + \item{$conf_low}{ +Corresponds to the upper limit of the \code{siglev}\% (only present if \code{conf = TRUE}) for the correlation. + } + \item{$conf_high}{ +Corresponds to the lower limit of the \code{siglev}\% (only present if \code{conf = TRUE}) for the correlation. + } + } } \examples{ -# Load sample data as in Load() example: -example(Load) -clim <- Clim(sampleData$mod, sampleData$obs) -ano_exp <- Ano(sampleData$mod, clim$clim_exp) -ano_obs <- Ano(sampleData$obs, clim$clim_obs) -runmean_months <- 12 -dim_to_smooth <- 4 # Smooth along lead-times -smooth_ano_exp <- Smoothing(ano_exp, runmean_months, dim_to_smooth) -smooth_ano_obs <- Smoothing(ano_obs, runmean_months, dim_to_smooth) -dim_to_mean <- 2 # Mean along members -required_complete_row <- 3 # Discard start dates which contain any NA lead-times -leadtimes_per_startdate <- 60 -corr <- Corr(Mean1Dim(smooth_ano_exp, dim_to_mean), - Mean1Dim(smooth_ano_obs, dim_to_mean), - compROW = required_complete_row, - limits = c(ceiling((runmean_months + 1) / 2), - leadtimes_per_startdate - floor(runmean_months / 2))) -PlotVsLTime(corr, toptitle = "correlations", ytitle = "correlation", - monini = 11, limits = c(-1, 2), listexp = c('CMIP5 IC3'), - listobs = c('ERSST'), biglab = FALSE, hlines = c(-1, 0, 1), +# Load sample data as in Load() example: +example(Load) +clim <- Clim(sampleData$mod, sampleData$obs) +ano_exp <- Ano(sampleData$mod, clim$clim_exp) +ano_obs <- Ano(sampleData$obs, clim$clim_obs) +runmean_months <- 12 +dim_to_smooth <- 4 +# Smooth along lead-times +smooth_ano_exp <- Smoothing(ano_exp, runmean_months, dim_to_smooth) +smooth_ano_obs <- Smoothing(ano_obs, runmean_months, dim_to_smooth) +dim_to_mean <- 2 # Mean along members +required_complete_row <- 3 # Discard start dates which contain any NA lead-times +leadtimes_per_startdate <- 60 +corr <- Corr(Mean1Dim(smooth_ano_exp, dim_to_mean), + Mean1Dim(smooth_ano_obs, dim_to_mean), + compROW = required_complete_row, + limits = c(ceiling((runmean_months + 1) / 2), + leadtimes_per_startdate - floor(runmean_months / 2))) +PlotVsLTime(corr, toptitle = "correlations", ytitle = "correlation", + monini = 11, limits = c(-1, 2), listexp = c('CMIP5 IC3'), + listobs = c('ERSST'), biglab = FALSE, hlines = c(-1, 0, 1), fileout = 'tos_cor.eps') + +# The following example uses veriApply combined with .Corr instead of Corr + \dontrun{ +require(easyVerification) +Corr2 <- s2dverification:::.Corr +corr2 <- veriApply("Corr2", + smooth_ano_exp, + # see ?veriApply for how to use the 'parallel' option + Mean1Dim(smooth_ano_obs, dim_to_mean), + tdim = 3, ensdim = 2) + } } \author{ History:\cr -0.1 - 2011-04 (V. Guemas, \email{vguemas at ic3.cat}) - Original code\cr -1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at ic3.cat}) - Formatting to R CRAN\cr -1.1 - 2014-10 (M. Menegoz, \email{martin.menegoz at ic3.cat}) - Adding siglev argument\cr -1.2 - 2015-03 (L.P. Caron, \email{louis-philippe.caron at ic3.cat}) - Adding method argument +0.1 - 2011-04 (V. Guemas, \email{vguemas at bsc.es}) - Original code\cr +1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at bsc.es}) - Formatting to R CRAN\cr +1.1 - 2014-10 (M. Menegoz, \email{martin.menegoz at bsc.es}) - Adding siglev argument\cr +1.2 - 2015-03 (L.P. Caron, \email{louis-philippe.caron at bsc.es}) - Adding method argument\cr +1.3 - 2017-02 (A. Hunter, \email{alasdair.hunter at bsc.es}) - Adapted to veriApply() } \keyword{datagen} diff --git a/man/EOF.Rd b/man/EOF.Rd index 6b3fc9a476746ee74435a88775f75b7dd14cc952..7c914e1570c2c415ae41d570448d0385290aa461 100644 --- a/man/EOF.Rd +++ b/man/EOF.Rd @@ -5,7 +5,7 @@ Area-Weighted Empirical Orthogonal Function Analysis Using SVD } \description{ Performs an area-weighted EOF analysis using SVD based on a covariance matrix -by default, based on a correlation matrix if \code{corr} argument is set to +by default, based on the correlation matrix if \code{corr} argument is set to \code{TRUE}. } \usage{ diff --git a/man/Eno.Rd b/man/Eno.Rd index bd93f9588997c147e09eef4144a6b25fb5ed5995..8a5e0b5f9a40582fbe573e0a9f2f74965b2e6207 100644 --- a/man/Eno.Rd +++ b/man/Eno.Rd @@ -4,8 +4,8 @@ Computes Effective Sample Size With Classical Method } \description{ -Computes the effective number of independant data along the posdim dimension of a matrix.\cr -This effective number of independant date may be required to perform statistical/inference tests.\cr +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. } \usage{ diff --git a/man/EnoNew.Rd b/man/EnoNew.Rd index a1dab6681d61df44198afced3b35783c78ea50f2..0cec23b452c87389481ddca2ee422d6097313262 100644 --- a/man/EnoNew.Rd +++ b/man/EnoNew.Rd @@ -2,14 +2,14 @@ \alias{EnoNew} \title{Computes Effective Sample Size Following Guemas et al, BAMS, 2013b} \description{ -This function computes the equivalent number of independent data in the xdata array following the method described in Guemas V., Auger L., Doblas-Reyes F., JAMC, 2013. The method relies on the Trenberth (1984) formula combined with a reduced uncertainty of the estimated autocorrelation function compared to the original approach.} +This function computes the effective number of independent values in the xdata array following the method described in Guemas V., Auger L., Doblas-Reyes F., JAMC, 2013. \code{EnoNew} provides similar functionality to \code{Eno} but with the added options to remove the linear trend or filter the frequency.} \usage{ EnoNew(xdata, detrend = FALSE, filter = FALSE) } \arguments{ - \item{xdata}{Timeseries from which the equivalent number of independent data is requested} - \item{detrend}{TRUE applies a linear detrending to xdata prior to the estimation of the equivalent number of independant data.} - \item{filter}{TRUE applies a filtering of any frequency peak prior to the estimation of the equivalent number of independant data.} + \item{xdata}{a numeric vector} + \item{detrend}{should the linear trend be removed from the data prior to the estimation of the equivalent number of independent values.} + \item{filter}{should a filtering of the frequency peaks be applied prior to the estimation of the equivalent number of independant data.} } \examples{ # See examples on Load() to understand the first lines in this example @@ -50,6 +50,9 @@ sampleData <- s2dverification:::.LoadSampleData('tos', c('experiment'), eno <- EnoNew(sampleData$mod[1, 1, , 1, 2, 3]) print(eno) } +\references{ +Guemas V, Auger L, Doblas-Reyes FJ, Rust H, Ribes A, 2014, Dependencies in Statistical Hypothesis Tests for Climate Time Series. Bulletin of the American Meteorological Society, 95 (11), 1666-1667. +} \author{ History:\cr 0.1 - 2012-06 (V. Guemas, \email{virginie.guemas at ic3.cat}) - Original code\cr diff --git a/man/Filter.Rd b/man/Filter.Rd index 6029692d2f29f442923c28917d837ba5f2146765..bb1bcf44f79f867407eeae6d0854e8158c8aa15a 100644 --- a/man/Filter.Rd +++ b/man/Filter.Rd @@ -2,9 +2,9 @@ \alias{Filter} \title{Filter Frequency Peaks From An Array} \description{ -This function filters from the xdata array, the signal of frequency freq.\cr -The filtering is performed by dichotomy, seeking for the frequency around freq and the phase that maximizes the signal to subtract to xdata.\cr -The maximization of the signal to subtract relies on a minimization of the mean square differences between xdata and a cosine of given frequency and phase. +This function filters out the selected frequency from a time series.\cr +The filtering is performed by dichotomy, seeking for a frequency around the parameter \code{freq} and the phase that maximizes the signal to subtract from the time series.\cr +The maximization of the signal to subtract relies on a minimization of the mean square differences between the time series (xdata) and the cosine of the specified frequency and phase. } \usage{ Filter(xdata, freq) diff --git a/man/FitAcfCoef.Rd b/man/FitAcfCoef.Rd index 1be3af50ccbd7cbbf0765f0ebe99690418372aba..73db6dfc96fe820d76c583082d6e2db476d50f5c 100644 --- a/man/FitAcfCoef.Rd +++ b/man/FitAcfCoef.Rd @@ -3,8 +3,8 @@ \title{Fits an AR1 AutoCorrelation Function Using the Cardano Formula} \description{ This function finds the minimum point of the fourth order polynom (a - x)2 + 0.25(b - x2)2 written to fit the two autoregression coefficients a and b.\cr -Thanks to the Cardano formula, provided a and b in [0 1], the problem is well posed, delta > 0 and there is only one solution to the minimum.\cr\cr -This function is called in Alpha() to minimize the mean square differences between the theoretical autocorrelation function of an AR1 and the first guess of estimated autocorrelation function estacf, using only the first two lags.} +A consequence of the Cardano formula is that, provided a and b are in [0 1], the problem is well posed, delta > 0 and there is only one minimum.\cr\cr +This function is called in Alpha() to minimize the mean square differences between the theoretical autocorrelation function of an AR1 and the first guess of the estimated autocorrelation function estacf, using only the first two lags.} \usage{FitAcfCoef(a, b)} \arguments{ \item{a}{Coefficient a : first estimate of the autocorrelation at lag 1} diff --git a/man/FitAutocor.Rd b/man/FitAutocor.Rd index 2474453fab60f04fb4a557665078a25712227fe5..d5031fdb74167cc7e3fcf770e835d2c548c965f9 100644 --- a/man/FitAutocor.Rd +++ b/man/FitAutocor.Rd @@ -1,12 +1,12 @@ \name{FitAutocor} \alias{FitAutocor} \title{Fits an AR1 Autocorrelation Function Using Dichotomy} -\description{This function fits the theoretical autocorrelation function of an AR1 to the first guess of estimated autocorrelation function estacf containing any number of lags. The fitting relies on a dichotomial minimisation of the mean square differences between both autocorrelation functions. It returns the autocorrelation at lag 1 of the fitted AR1 process.} +\description{This function fits the theoretical autocorrelation function of an AR1 to the first guess of the estimated autocorrelation function estacf containing any number of lags. The fitting relies on a dichotomial minimisation of the mean square differences between both autocorrelation functions. It returns the autocorrelation at lag 1 of the fitted AR1 process.} \usage{ FitAutocor(estacf, window = c(-1, 1), prec = 0.01) } \arguments{ - \item{estacf}{First guess of the autocorrelation function} + \item{estacf}{First guess for the autocorrelation function} \item{window}{Interval in which the autocorrelation at lag 1 should be found.} \item{prec}{Precision to which the autocorrelation function at lag 1 is to be estimated.} } diff --git a/man/GenSeries.Rd b/man/GenSeries.Rd index 2162e92782a14e06ea2cba92dd57cb9f21162394..e616d6a13c528366040c48764317a937bdec797d 100644 --- a/man/GenSeries.Rd +++ b/man/GenSeries.Rd @@ -1,7 +1,7 @@ \name{GenSeries} \alias{GenSeries} \title{Generates An AR1 Time Series} -\description{This functions generates AR1 processes containing n data, with alpha as autocorrelation at lag 1, and mean and standard deviation provided by the mean and std arguments.} +\description{This function generates AR1 processes containing n data points, where alpha is the autocorrelation at lag 1, and the mean and standard deviation are specified by the mean and std arguments.} \usage{ GenSeries(n, alpha, mean, std) } diff --git a/man/Histo2Hindcast.Rd b/man/Histo2Hindcast.Rd index cf9e2aa8b9b3fb79f67541e60b081852b4c20a65..3764f66134d5ac463840108b903c50001a56fec1 100644 --- a/man/Histo2Hindcast.Rd +++ b/man/Histo2Hindcast.Rd @@ -11,7 +11,7 @@ Histo2Hindcast(varin, sdatesin, sdatesout, nleadtimesout) } \arguments{ \item{varin}{ -Input model or observational data:\cr +Array of model or observational data with dimensions:\cr c(nmod/nexp/nobs, nmemb/nparam, nsdates, nltimes) up to\cr c(nmod/nexp/nobs, nmemb/nparam, nsdates, nltimes, nlevel, nlat, nlon) } @@ -26,7 +26,7 @@ Number of leadtimes in the output matrix. } } \value{ -A matrix with the same number of dimensions as the input one, the same dimensions 1 and 2 and potentially the same dimensions 5 to 7. Dimensions 3 and 4 are set by the arguments sdatesout and nleadtimesout. +An array with the same number of dimensions as varin, the same dimensions 1 and 2 and potentially the same dimensions 5 to 7. Dimensions 3 and 4 are set by the arguments sdatesout and nleadtimesout. } \examples{ # See examples on Load() to understand the first lines in this example diff --git a/man/IniListDims.Rd b/man/IniListDims.Rd index a171dd968f0ad63aac268c6e4319d761cebbab23..0af0253c16911fa06ee6c8ee364383460e45d50f 100644 --- a/man/IniListDims.Rd +++ b/man/IniListDims.Rd @@ -4,7 +4,7 @@ Creates A List Of Integer Ranges } \description{ -This function generates a list of arrays where those arrays contain integers from 1 to various numbers. This list of arrays is used in the other functions as a list of indices of the elements of the matrices. +This function generates a list of arrays containing integers greater than or equal to 1. This list of arrays is used in other functions as a list of indices of the elements of the matrices. } \usage{ IniListDims(dims, lenlist) diff --git a/man/InsertDim.Rd b/man/InsertDim.Rd index d3708b899688703ce3c3689ed0ca5ab34a5b783a..6b08a5830d50125f3a11633fba3de67dd51f1d74 100644 --- a/man/InsertDim.Rd +++ b/man/InsertDim.Rd @@ -1,10 +1,10 @@ \name{InsertDim} \alias{InsertDim} \title{ -Adds A Dimension To A Matrix +Adds A Dimension To An Array } \description{ -Add one dimension to the matrix 'var' in position 'posdim' with length 'lendim' and which correspond to 'lendim' repetitions of the 'var' matrix. +Inserts an extra dimension into an array at position 'posdim' with length 'lendim' and which correspond to 'lendim' repetitions of the 'var' array. } \usage{ InsertDim(var, posdim, lendim) diff --git a/man/LeapYear.Rd b/man/LeapYear.Rd index b09d21573ba91cfeedc1f2762abe0473c1513e22..faf414567b68a50863a5ebe56f51a95cd141e26c 100644 --- a/man/LeapYear.Rd +++ b/man/LeapYear.Rd @@ -2,14 +2,14 @@ \alias{LeapYear} \title{Checks Whether A Year Is Leap Year} \description{ -This function tells whether a year is leap year or not. +This function tells whether a year is a leap year or not. } \usage{ LeapYear(year) } \arguments{ \item{year}{ -The year to tell whether is leap year or not. +A numeric value indicating the year in the Gregorian calendar. } } \value{ diff --git a/man/Mean1Dim.Rd b/man/Mean1Dim.Rd index e1df907c7e7fe3cb79ab3a9e9fd29f6cf6fbf723..0d56d083327c9db4d8f6b80de3157878d9fbdd9c 100644 --- a/man/Mean1Dim.Rd +++ b/man/Mean1Dim.Rd @@ -1,10 +1,10 @@ \name{Mean1Dim} \alias{Mean1Dim} \title{ -Averages A Matrix Along A Dimension +Averages An Array Along A Dimension } \description{ -Averages the matrix var along the posdim dimension between limits [1] and limits [2] if limits argument is provided by the user. +Averages the array along the posdim dimension along the user specified dimension. The user can specify a subset of the dimension to take the mean along. } \usage{ Mean1Dim(var, posdim, narm = TRUE, limits = NULL) @@ -20,11 +20,11 @@ Dimension to average along. Ignore NA (TRUE) values or not (FALSE). } \item{limits}{ -Limits to average between. +Limits to average between. Default is to take the mean along the entire dimension. } } \value{ -Matrix with one dimension less than the input one containing the average along posdim dimension. +Array with one dimension less than the input array, containing the average along the posdim dimension. } \examples{ a <- array(rnorm(24), dim = c(2, 3, 4)) diff --git a/man/MeanListDim.Rd b/man/MeanListDim.Rd index 547d30db7a24fb4748534485899e6f5ba64c4f26..f3667d7543e394ee1c1131e533054ba62e985556 100644 --- a/man/MeanListDim.Rd +++ b/man/MeanListDim.Rd @@ -1,17 +1,17 @@ \name{MeanListDim} \alias{MeanListDim} \title{ -Averages A Matrix Along Various Dimensions +Averages An Array Along Multiple Dimensions } \description{ -Averages the matrix var along a set of dimensions given by the argument dims. +Averages an array along a set of dimensions given by the argument dims. } \usage{ MeanListDim(var, dims, narm = TRUE) } \arguments{ \item{var}{ -Matrix to average. +Input array. } \item{dims}{ List of dimensions to average along. @@ -21,7 +21,7 @@ Ignore NA (TRUE) values or not (FALSE). } } \value{ -Matrix with as many dimensions less than the input matrix as provided by the list dims and containing the average along this list of dimensions. +The averaged array, with the dimensions specified in \code{dims} removed. } \examples{ a <- array(rnorm(24), dim = c(2, 3, 4)) diff --git a/man/NAO.Rd b/man/NAO.Rd index 1f3360aeba26b8a1354167b59cb7c2b838b21b36..c0963652762d4fade328e86cb17bca4927424260 100644 --- a/man/NAO.Rd +++ b/man/NAO.Rd @@ -1,11 +1,11 @@ \name{NAO} \alias{NAO} \title{ -Compute the North Atlantic Oscillation (NAO) Index +Computes the North Atlantic Oscillation (NAO) Index } \description{ Compute the North Atlantic Oscillation (NAO) index based on -the leading EOF of sea level pressure (SLP) anomalies over the +the leading EOF of the sea level pressure (SLP) anomalies over the north Atlantic region (20N-80N, 80W-40E). The PCs are obtained by projecting the forecast and observed anomalies onto the observed EOF pattern (Pobs) or the forecast anomalies onto the EOF pattern of the other years of the forecast @@ -13,9 +13,6 @@ forecast anomalies onto the EOF pattern of the other years of the forecast 1-month lead seasonal forecasts that can be plotted with BoxPlot(). Returns cross-validated PCs of the NAO index for forecast (ano_exp) and observations (ano_obs) based on the leading EOF pattern.\cr -See Doblas-Reyes, F.J., Pavan, V. and Stephenson, D. (2003). The skill of -multi-model seasonal forecasts of the wintertime North Atlantic Oscillation. -Climate Dynamics, 21, 501-514. DOI: 10.1007/s00382-003-0350-4 } \usage{ NAO(ano_exp = NULL, ano_obs = NULL, lon, lat, ftime_average = 2:4, @@ -115,6 +112,11 @@ nao <- NAO(ano$ano_exp, ano$ano_obs, sampleData$lon, sampleData$lat) PlotBoxWhisker(nao$NAO_exp, nao$NAO_obs, "NAO index, DJF", "NAO index (PC1) TOS", monini = 12, yearini = 1985, freq = 1, "Exp. A", "Obs. X") } +\references{ +Doblas-Reyes, F.J., Pavan, V. and Stephenson, D. (2003). The skill of +multi-model seasonal forecasts of the wintertime North Atlantic Oscillation. +Climate Dynamics, 21, 501-514. DOI: 10.1007/s00382-003-0350-4 +} \author{ History:\cr 0.1 - 2013-08 (F. Lienert, \email{flienert at ic3.cat}) - Original code\cr diff --git a/man/PlotACC.Rd b/man/PlotACC.Rd index 21cfd80a8fae8f5888976656e2956f4ed1b3fb2a..dc4b5c8759e70c29eed21e081c219ca0368f6c8d 100644 --- a/man/PlotACC.Rd +++ b/man/PlotACC.Rd @@ -4,9 +4,9 @@ Plot Plumes/Timeseries Of Anomaly Correlation Coefficients } \description{ -Plots plumes/timeseries of ACC from a matrix with dimensions (output from \code{ACC()}): \cr +Plots plumes/timeseries of ACC from an array with dimensions (output from \code{ACC()}): \cr c(nexp, nobs, nsdates, nltime, 4)\cr -with the fourth dimension of length 4 containing the lower limit of the 95\% confidence interval, the ACC, the upper limit of the 95\% confidence interval and the 95\% significance level given by a one-sided T-test. +where the fourth dimension is of length 4 and contains the lower limit of the 95\% confidence interval, the ACC, the upper limit of the 95\% confidence interval and the 95\% significance level given by a one-sided T-test. } \usage{ PlotACC(ACC, sdates, toptitle = "", sizetit = 1, ytitle = "", limits = NULL, diff --git a/man/PlotAno.Rd b/man/PlotAno.Rd index ec3d6de74935c5d5e458846449647f6a3de05d7a..c0185e41b0b9c4fec168adf252041cb6d883c5bb 100644 --- a/man/PlotAno.Rd +++ b/man/PlotAno.Rd @@ -4,9 +4,7 @@ Plot Raw Or Smoothed Anomalies } \description{ -Plots timeseries of raw or smoothed anomalies of any index output from \code{Load()} or \code{Ano()} or or \code{Ano_CrossValid()} or \code{Smoothing()} and organized in matrices with dimensions:\cr - c(nmod/nexp, nmemb/nparam, nsdates, nltime) for the model data\cr - c(nobs, nmemb, nsdates, nltime) for the observational data +Plots timeseries of raw or smoothed anomalies of any variable output from \code{Load()} or \code{Ano()} or or \code{Ano_CrossValid()} or \code{Smoothing()}. } \usage{ PlotAno(exp_ano, obs_ano = NULL, sdates, toptitle = c("", "", "", "", "", "", diff --git a/man/PlotEquiMap.Rd b/man/PlotEquiMap.Rd index 9ba2ca4909ffdcb5da86232d6ffa99437ac14366..cae7ec3f1f5709d2d9b9664ab70639b7d4caeda0 100644 --- a/man/PlotEquiMap.Rd +++ b/man/PlotEquiMap.Rd @@ -4,7 +4,7 @@ Maps A Two-Dimensional Variable On A Cylindrical Equidistant Projection } \description{ -Map longitude-latitude array (on a regular rectangular or gaussian grid) on a cylindrical equidistant latitude and longitude world projection with coloured grid cells. Only the region for which data has been provided is displayed. A colour bar (legend) can be plotted and adjusted. It is possible to draw superimposed arrows, dots, symbols, contour lines and boxes. A number of options is provided to adjust the position, size and colour of the components. This plot function is compatible with figure layouts if colour bar is disabled. +Map longitude-latitude array (on a regular rectangular or gaussian grid) on a cylindrical equidistant latitude and longitude projection with coloured grid cells. Only the region for which data has been provided is displayed. A colour bar (legend) can be plotted and adjusted. It is possible to draw superimposed arrows, dots, symbols, contour lines and boxes. A number of options is provided to adjust the position, size and colour of the components. This plot function is compatible with figure layouts if colour bar is disabled. } \usage{ PlotEquiMap(var, lon, lat, varu = NULL, varv = NULL, diff --git a/man/ProbBins.Rd b/man/ProbBins.Rd index 96ec05767efa099e6678380311beff046b6b704b..8bb763301aa539d45def771aa53ed272f5c8fc3e 100644 --- a/man/ProbBins.Rd +++ b/man/ProbBins.Rd @@ -1,13 +1,13 @@ \name{ProbBins} \alias{ProbBins} \title{ -Computes probabilistic information of a forecast relative to a threshold or a quantile. +Computes Probabilistic Information of a Forecast Relative to a Threshold or a Quantile. } \description{ 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 +Array with probabilistic information and dimensions:\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} diff --git a/man/ProjectField.Rd b/man/ProjectField.Rd index 6d49bf2b19d9f20ecb4b6647e01c66ceaaf08853..91d33039da5d432fc3df07f93b95d0900f35478c 100644 --- a/man/ProjectField.Rd +++ b/man/ProjectField.Rd @@ -4,7 +4,7 @@ Project Anomalies onto Modes of Variability } \description{ -Project anomalies onto modes to get temporal evolution of the EOF mode\cr +Project anomalies onto modes of variability to get the temporal evolution of the EOF mode\cr selected. Returns principal components (PCs) by area-weighted projection onto EOF pattern (from \code{EOF()}). Able to handle NAs. } diff --git a/man/RMS.Rd b/man/RMS.Rd index 118400f407167bcd5d3d8b88953669bc882475bc..34980b658275f642a6ac7576a15175992e4c3034 100644 --- a/man/RMS.Rd +++ b/man/RMS.Rd @@ -1,18 +1,22 @@ \name{RMS} \alias{RMS} +\alias{.RMS} \title{ -Computes Root Mean Square Error Skill Measure +Computes Root Mean Square Error } \description{ -Matrix var_exp & var_obs should have the same dimensions except along posloop dimension where the length can be different, with the number of experiments/models for var_exp (nexp) and the number of obserational datasets for var_obs (nobs).\cr -RMS computes the Root Mean Square Error skill of each jexp in 1:nexp against each jobs in 1:nobs which gives nexp x nobs RMSE skill measures for each other grid point of the matrix (each latitude/longitude/level/leadtime).\cr -The RMSE are computed along the posRMS dimension which should correspond to the startdate dimension.\cr -If compROW is given, the RMSE are computed only if rows along the compROW dimension are complete between limits[1] and limits[2], that mean with no NA between limits[1] and limits[2]. This option can be activated if the user wishes to account only for the forecasts for which observations are available at all leadtimes.\cr +Computes the root mean square error for an array of forecasts, var_exp and an array of observations, var_obs, which should have the same dimensions except along the posloop dimension where the lengths can be different, with the number of experiments/models for var_exp (nexp) and the number of obserational datasets for var_obs (nobs).\cr +The RMSE is computed along the posRMS dimension which should correspond to the startdate dimension.\cr +If compROW is given, the RMSE is computed only if rows along the compROW dimension are complete between limits[1] and limits[2], i.e. there are no NAs between limits[1] and limits[2]. This option can be activated if the user wishes to account only for the forecasts for which observations are available at all leadtimes.\cr Default: limits[1] = 1 and limits[2] = length(compROW dimension).\cr -The confidence interval relies on a chi2 distribution. +The confidence interval relies on a chi2 distribution.\cr +\cr +.RMS provides the same functionality but taking a matrix of ensemble members as input (exp). } \usage{ RMS(var_exp, var_obs, posloop = 1, posRMS = 2, compROW = NULL, limits = NULL, siglev = 0.95, conf = TRUE) + +.RMS(exp, obs, siglev = 0.95, conf = TRUE) } \arguments{ \item{var_exp}{ @@ -40,11 +44,28 @@ Confidence level of the computed confidence interval. 0.95 by default. \item{conf}{ Whether to compute confidence interval or not. TRUE by default. } + \item{exp}{ +N by M matrix of N forecasts from M ensemble members. + } + \item{obs}{ +Vector of the corresponding observations of length N. + } } \value{ -Matrix with dimensions:\cr +RMS: Array with dimensions:\cr c(length(posloop) in var_exp, length(posloop) in var_obs, 1 or 3, all other dimensions of var_exp & var_obs except posRMS).\cr -The 3rd dimension corresponds to the lower limit of the 95\% confidence interval (only present if \code{conf = TRUE}), the RMSE, and the upper limit of the 95\% confidence interval (only present if \code{conf = TRUE}). +The 3rd dimension corresponds to the lower limit of the 95\% confidence interval (only present if \code{conf = TRUE}), the RMSE, and the upper limit of the 95\% confidence interval (only present if \code{conf = TRUE}). \cr +\cr +.RMS: + \item{$rms}{ +The root mean square error, + } + \item{$conf_low}{ +Corresponding to the lower limit of the \code{siglev}\% confidence interval (only present if \code{conf = TRUE}) for the rms. + } + \item{$conf_high}{ +Corresponding to the upper limit of the \code{siglev}\% confidence interval (only present if \code{conf = TRUE}) for the rms. + } } \examples{ # Load sample data as in Load() example: @@ -68,10 +89,21 @@ PlotVsLTime(rms, toptitle = "Root Mean Square Error", ytitle = "K", monini = 11, limits = NULL, listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, hlines = c(0), fileout = 'tos_rms.eps') +# The following example uses veriApply combined with .RMS instead of RMS + \dontrun{ +require(easyVerification) +RMS2 <- s2dverification:::.RMS +rms2 <- veriApply("RMS2", + smooth_ano_exp, + # see ?veriApply for how to use the 'parallel' option + Mean1Dim(smooth_ano_obs, dim_to_mean), + tdim = 3, ensdim = 2) + } } \author{ History:\cr 0.1 - 2011-05 (V. Guemas, \email{vguemas at ic3.cat}) - Original code\cr -1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at ic3.cat}) - Formatting to R CRAN +1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at ic3.cat}) - Formatting to R CRAN\cr +1.1 - 2017-02 (A. Hunter, \email{alasdair.hunter at bsc.es}) - Adapted to veriApply() } \keyword{datagen} diff --git a/man/RMSSS.Rd b/man/RMSSS.Rd index 3d77230dbf7156bc3f96646126c5b7b42c637cf9..b7a994a7765784cfcb6c036209261a340606839b 100644 --- a/man/RMSSS.Rd +++ b/man/RMSSS.Rd @@ -1,16 +1,21 @@ \name{RMSSS} \alias{RMSSS} +\alias{.RMSSS} \title{ Computes Root Mean Square Skill Score } \description{ -Arrays var_exp & var_obs should have the same dimensions except along posloop where the length can be different, with the number of experiments/models for var_exp (nexp) and the number of obserational datasets for var_obs (nobs).\cr +Computes the root mean square error skill score between an array of forecasts, var_exp and an array of observations, var_obs, which should have the same dimensions except along posloop where the lengths can be different, with the number of experiments/models for var_exp (nexp) and the number of obserational datasets for var_obs (nobs).\cr RMSSS computes the Root Mean Square Skill Score of each jexp in 1:nexp against each jobs in 1:nobs which gives nexp x nobs RMSSS for each other grid point of the matrix (each latitude/longitude/level/leadtime).\cr The RMSSS are computed along the posRMS dimension which should correspond to the startdate dimension.\cr -The p-value is optionally provided by a one-sided Fisher test. +The p-value is optionally provided by a one-sided Fisher test.\cr +\cr +.RMSSS provides the same functionality but taking a matrix of ensemble members as input (exp). } \usage{ RMSSS(var_exp, var_obs, posloop = 1, posRMS = 2, pval = TRUE) + +.RMSSS(exp, obs, pval = TRUE) } \arguments{ \item{var_exp}{ @@ -28,11 +33,27 @@ Dimension along which the RMSE are to be computed (the dimension of the start da \item{pval}{ Whether to compute or not the p-value of the test Ho : RMSSS = 0. TRUE by default. } + \item{exp}{ +N by M matrix of N forecasts from M ensemble members. + } + \item{obs}{ +Vector of the corresponding observations of length N. + } } \value{ -Array with dimensions:\cr +RMSSS: Array with dimensions:\cr c(length(posloop) in var_exp, length(posloop) in var_obs, 1 or 2, all other dimensions of var_exp & var_obs except posRMS).\cr -The 3rd dimension corresponds to the RMSSS and, if \code{pval = TRUE}, the p-value of the one-sided Fisher test with Ho: RMSSS = 0. +The 3rd dimension corresponds to the RMSSS and, if \code{pval = TRUE}, the p-value of the one-sided Fisher test with Ho: RMSSS = 0.\cr +\cr +.RMSSS: + \itemize{ + \item{$rmsss}{ +The RMSSS. + } + \item{$p_val}{ +Corresponds to the p values (only present if \code{pval = TRUE}) for the RMSSS. + } + } } \examples{ # Load sample data as in Load() example: @@ -41,17 +62,27 @@ clim <- Clim(sampleData$mod, sampleData$obs) ano_exp <- Ano(sampleData$mod, clim$clim_exp) ano_obs <- Ano(sampleData$obs, clim$clim_obs) rmsss <- RMSSS(Mean1Dim(ano_exp, 2), Mean1Dim(ano_obs, 2)) -rmsss2 <- array(dim = c(dim(rmsss)[1:2], 4, dim(rmsss)[4])) -rmsss2[, , 2, ] <- rmsss[, , 1, ] -rmsss2[, , 4, ] <- rmsss[, , 2, ] -PlotVsLTime(rmsss, toptitle = "Root Mean Square Skill Score", ytitle = "", +rmsss_plot <- array(dim = c(dim(rmsss)[1:2], 4, dim(rmsss)[4])) +rmsss_plot[, , 2, ] <- rmsss[, , 1, ] +rmsss_plot[, , 4, ] <- rmsss[, , 2, ] +PlotVsLTime(rmsss_plot, toptitle = "Root Mean Square Skill Score", ytitle = "", monini = 11, limits = c(-1, 1.3), listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, hlines = c(-1, 0, 1), fileout = 'tos_rmsss.eps') +# The following example uses veriApply combined with .RMSSS instead of RMSSS + \dontrun{ +require(easyVerification) +RMSSS2 <- s2dverification:::.RMSSS +rmsss2 <- veriApply("RMSSS2", ano_exp, + # see ?veriApply for how to use the 'parallel' option + Mean1Dim(ano_obs, 2), + tdim = 3, ensdim = 2) + } } \author{ History:\cr -0.1 - 2012-04 (V. Guemas, \email{vguemas at ic3.cat}) - Original code\cr -1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at ic3.cat}) - Formatting to R CRAN +0.1 - 2012-04 (V. Guemas, \email{vguemas at bsc.es}) - Original code\cr +1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at bsc.es}) - Formatting to R CRAN\cr +1.1 - 2017-02 (A. Hunter, \email{alasdair.hunter at bsc.es}) - Adapted to veriApply() } \keyword{datagen} diff --git a/man/RatioRMS.Rd b/man/RatioRMS.Rd index ee116561f7662d6dfa0ffc15bc3ddf5ceff06b22..aa179b769e2b9603fcf1dd056c113feed22c84f6 100644 --- a/man/RatioRMS.Rd +++ b/man/RatioRMS.Rd @@ -1,25 +1,30 @@ \name{RatioRMS} \alias{RatioRMS} +\alias{.RatioRMS} \title{ -Computes The Ratio Between The RMSE Scores of 2 Experiments. +Computes the Ratio Between The RMSE of Two Experiments } \description{ -Matrix var_exp1 / var_exp2 / var_obs should have the same dimensions.\cr -The ratio RMSE(var_exp1, var_obs) / RMSE(var_exp2, var_obs) is output.\cr -The p-value is provided by a two-sided Fischer test. +Calculates the ratio of the RMSE for two forecasts of the same observations.\cr +The ratio RMSE(ens, obs) / RMSE(ens.ref, obs) is output.\cr +The p-value is provided by a two-sided Fischer test.\cr +\cr +.RatioRMS provides the same functionality but taking two matrices of ensemble members (ens and ens.ref) as input. } \usage{ RatioRMS(var_exp1, var_exp2, var_obs, posRMS = 1, pval = TRUE) + +.RatioRMS(exp, exp_ref, obs, pval = TRUE) } \arguments{ \item{var_exp1}{ -Matrix of experimental data 1. +Array of experimental data 1. } \item{var_exp2}{ -Matrix of experimental data 2, same dimensions as var_exp1. +Array of experimental data 2. } \item{var_obs}{ -Matrix of observational data, same dimensions as var_exp1. +Array of observations. } \item{posRMS}{ Dimension along which the RMSE are to be computed = the position of the start dates. @@ -27,10 +32,25 @@ Dimension along which the RMSE are to be computed = the position of the start da \item{pval}{ Whether to compute the p-value of Ho : RMSE1/RMSE2 = 1 or not. TRUE by default. } + \item{exp}{ +Matrix of experimental data 1. + } + \item{exp_ref}{ +Matrix of experimental data 2. + } + \item{obs}{ +Vector of observations. + } } \value{ -Matrix with the same dimensions as var_exp1/var_exp2/var_obs except along posRMS where the dimension has length 2 if \code{pval = TRUE}, or 1 otherwise.\cr -The dimension corresponds to the ratio between the RMSE (RMSE1/RMSE2) and the p.value of the two-sided Fisher test with Ho: RMSE1/RMSE2 = 1. +RatioRMS:\cr +Matrix with the same dimensions as var_exp1/var_exp2/var_obs except along posRMS where the dimension has length 2 if 'pval = TRUE', or 1 otherwise. The dimension of length 2 corresponds to the ratio between the RMSE (RMSE1/RMSE2) and the p-value of the two-sided Fisher test with Ho: RMSE1/RMSE2 = 1.\cr +\cr +.RatioRMS:\cr + \itemize{ + \item{ratiorms}{The ratio of the RMSE of the two experimental datasets} + \item{p_val}{The p-value} + } } \examples{ # See examples on Load() to understand the first lines in this example @@ -57,6 +77,7 @@ sampleData <- s2dverification:::.LoadSampleData('tos', c('experiment'), latmin = 27, latmax = 48, lonmin = -12, lonmax = 40) } +# Compute DJF seasonal means and anomalies. leadtimes_dimension <- 4 initial_month <- 11 mean_start_month <- 12 @@ -68,15 +89,39 @@ sampleData$obs <- Season(sampleData$obs, leadtimes_dimension, initial_month, clim <- Clim(sampleData$mod, sampleData$obs) ano_exp <- Ano(sampleData$mod, clim$clim_exp) ano_obs <- Ano(sampleData$obs, clim$clim_obs) -rrms <- RatioRMS(Mean1Dim(ano_exp[ , 1:2, , , , ], 1)[, 1, , ], - ano_exp[ , 3, , , , ][, 1, , ], - Mean1Dim(ano_obs, 2)[1, , 1, , ], 1) +# Generate two experiments with 2 and 1 members from the only experiment +# available in the sample data. Take only data values for a single forecast +# time step. +ano_exp_1 <- Subset(ano_exp, 'member', c(1, 2)) +ano_exp_2 <- Subset(ano_exp, 'member', c(3)) +ano_exp_1 <- Subset(ano_exp_1, c('dataset', 'ftime'), list(1, 1), drop = 'selected') +ano_exp_2 <- Subset(ano_exp_2, c('dataset', 'ftime'), list(1, 1), drop = 'selected') +ano_obs <- Subset(ano_obs, c('dataset', 'ftime'), list(1, 1), drop = 'selected') +# Compute ensemble mean and provide as inputs to RatioRMS. +rrms <- RatioRMS(Mean1Dim(ano_exp_1, 1), + Mean1Dim(ano_exp_2, 1), + Mean1Dim(ano_obs, 1)) +# Plot the RatioRMS for the first forecast time step. PlotEquiMap(rrms[1, , ], sampleData$lon, sampleData$lat, toptitle = 'Ratio RMSE') + +# The following example uses veriApply combined with .RatioRMS instead of RatioRMS + \dontrun{ +require(easyVerification) +# The name of the function has to end in 'ss' in order for veriApply() to +# detect it as a skill score. +RatioRMSss <- s2dverification:::.RatioRMS +rrms2 <- veriApply("RatioRMSss", ano_exp_1, + # see ?veriApply for how to use the 'parallel' option + Mean1Dim(ano_obs, 1), + ano_exp_2, + tdim = 2, ensdim = 1) + } } \author{ History:\cr -0.1 - 2011-11 (V. Guemas, \email{vguemas at ic3.cat}) - Original code\cr -1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at ic3.cat}) - Formatting to R CRAN +0.1 - 2011-11 (V. Guemas, \email{vguemas at bsc.es}) - Original code\cr +1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at bsc.es}) - Formatting to R CRAN\cr +1.1 - 2017-02 (A. Hunter, \email{alasdair.hunter at bsc.es}) - Adapted to veriApply() } \keyword{datagen} diff --git a/man/RatioSDRMS.Rd b/man/RatioSDRMS.Rd index 3dcfcd96aeefbaf498e179e1859f918409e2eb58..1e6a9f9da7adeb4908582cf71d144d653c797cc4 100644 --- a/man/RatioSDRMS.Rd +++ b/man/RatioSDRMS.Rd @@ -1,21 +1,26 @@ \name{RatioSDRMS} \alias{RatioSDRMS} +\alias{.RatioSDRMS} \title{ -Computes The Ratio Between the Ensemble Spread and the RMSE of the Ensemble Mean +Computes the ratio between the ensemble spread and RMSE } \description{ -Matrices var_exp & var_obs should have dimensions between\cr +Arrays var_exp & var_obs should have dimensions between\cr c(nmod/nexp, nmemb/nparam, nsdates, nltime)\cr and\cr c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon)\cr The ratio between the standard deviation of the members around the ensemble mean in var_exp and the RMSE between var_exp and var_obs is output for each experiment and each observational dataset.\cr -The p-value is provided by a one-sided Fischer test. +The p-value is provided by a one-sided Fischer test.\cr +\cr +.RatioSDRMS provides the same functionality but taking a matrix of ensemble members as input (exp). } \usage{ RatioSDRMS(var_exp, var_obs, pval = TRUE) + +.RatioSDRMS(exp, obs, pval = TRUE) } \arguments{ - \item{var_exp}{ + \item{var_exp}{ Model data:\cr c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to\cr c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon) @@ -28,27 +33,56 @@ Observational data:\cr \item{pval}{ Whether to compute the p-value of Ho : SD/RMSE = 1 or not. } + \item{exp}{ +N by M matrix of N forecasts from M ensemble members. + } + \item{obs}{ +Vector of the corresponding observations of length N. + } } \value{ -Matrix with dimensions c(nexp/nmod, nobs, 1 or 2, nltime) up to - c(nexp/nmod, nobs, 1 or 2, nltime, nlevel, nlat, nlon). -The 3rd dimension corresponds to the ratio (SD/RMSE) and the p.value (only present if \code{pval = TRUE}) of the one-sided Fisher test with Ho: SD/RMSE = 1. +RatioSDRMS: Array with dimensions c(nexp/nmod, nobs, 1 or 2, nltime) up to + c(nexp/nmod, nobs, 1 or 2, nltime, nlevel, nlat, nlon). \cr +The 3rd dimension corresponds to the ratio (SD/RMSE) and the p.value (only present if \code{pval = TRUE}) of the one-sided Fisher test with Ho: SD/RMSE = 1.\cr +\cr +.RatioSDRMS: + \itemize{ + \item{$ratio}{ +The ratio of the ensemble spread and RMSE, + } + \item{$p_val}{ +Corresponds to the p values of the ratio (only present if \code{pval = TRUE}). + } + } } \examples{ # Load sample data as in Load() example: example(Load) rsdrms <- RatioSDRMS(sampleData$mod, sampleData$obs) -rsdrms2 <- array(dim = c(dim(rsdrms)[1:2], 4, dim(rsdrms)[4])) -rsdrms2[, , 2, ] <- rsdrms[, , 1, ] -rsdrms2[, , 4, ] <- rsdrms[, , 2, ] -PlotVsLTime(rsdrms2, toptitle = "Ratio ensemble spread / RMSE", ytitle = "", +# Reorder the data in order to plot it with PlotVsLTime +rsdrms_plot <- array(dim = c(dim(rsdrms)[1:2], 4, dim(rsdrms)[4])) +rsdrms_plot[, , 2, ] <- rsdrms[, , 1, ] +rsdrms_plot[, , 4, ] <- rsdrms[, , 2, ] +PlotVsLTime(rsdrms_plot, toptitle = "Ratio ensemble spread / RMSE", ytitle = "", monini = 11, limits = c(-1, 1.3), listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, siglev = TRUE, fileout = 'tos_rsdrms.eps') + +# The following example uses veriApply combined with .RatioSDRMS instead of RatioSDRMS + \dontrun{ +require(easyVerification) +RatioSDRMS2 <- s2dverification:::.RatioSDRMS +rsdrms2 <- veriApply("RatioSDRMS2", + sampleData$mod, + # see ?veriApply for how to use the 'parallel' option + Mean1Dim(sampleData$obs, 2), + tdim = 3, ensdim = 2) + } } \author{ History:\cr 0.1 - 2011-12 (V. Guemas, \email{virginie.guemas at ic3.cat}) - Original code\cr -1.0 - 2013-09 (N. Manubens, \email{nicolau-manubens at ic3.cat}) - Formatting to CRAN +1.0 - 2013-09 (N. Manubens, \email{nicolau-manubens at ic3.cat}) - Formatting to CRAN\cr +1.1 - 2017-02 (A. Hunter, \email{alasdair.hunter at bsc.es}) - Adapted to veriApply() } \keyword{datagen} diff --git a/man/Regression.Rd b/man/Regression.Rd index 0024299334dbf307e32d60dac9fe7c4d815b83f0..da4dca37fa59bae7ab33fc737abdb1700d8f9ffb 100644 --- a/man/Regression.Rd +++ b/man/Regression.Rd @@ -1,7 +1,7 @@ \name{Regression} \alias{Regression} \title{ -Computes The Regression Of A Matrix On Another Along A Dimension +Computes The Regression Of An Array On Another Along A Dimension } \description{ Computes the regression of the input matrice vary on the input matrice varx along the posREG dimension by least square fitting. Provides the slope of the regression, the associated confidence interval, and the intercept.\cr @@ -13,10 +13,10 @@ Regression(vary, varx, posREG = 2) } \arguments{ \item{vary}{ -Matrix of any number of dimensions up to 10. +Array of any number of dimensions up to 10. } \item{varx}{ -Matrix of any number of dimensions up to 10. Same dimensions as vary. +Array of any number of dimensions up to 10. Same dimensions as vary. } \item{posREG}{ Position along which to compute the regression. @@ -24,7 +24,7 @@ Position along which to compute the regression. } \value{ \item{$regression}{ -Matrix with same dimensions as varx and vary except along posREG dimension which is replaced by a length 4 dimension, corresponding to the lower limit of the 95\% confidence interval, the slope, the upper limit of the 95\% confidence interval and the intercept. +Array with same dimensions as varx and vary except along posREG dimension which is replaced by a length 4 dimension, corresponding to the lower limit of the 95\% confidence interval, the slope, the upper limit of the 95\% confidence interval and the intercept. } \item{$filtered}{ Same dimensions as vary filtered out from the regression onto varx along the posREG dimension. diff --git a/man/Season.Rd b/man/Season.Rd index 0e19b49054b316e25354a7921733ca62d2dab7d4..aa143e27366238a89d12610ad4d300ddd05d1bc7 100644 --- a/man/Season.Rd +++ b/man/Season.Rd @@ -4,20 +4,20 @@ Computes Seasonal Means } \description{ -Computes seasonal means on timeseries organized in a matrix of any number of dimensions up to 10 dimensions where the time dimension is one of those 10 dimensions. +Computes seasonal means on timeseries organized in a array of any number of dimensions up to 10 dimensions where the time dimension is one of those 10 dimensions. } \usage{ Season(var, posdim = 4, monini, moninf, monsup) } \arguments{ \item{var}{ -Matrix containing the timeseries along one of its dimensions. +Array containing the timeseries along one of its dimensions. } \item{posdim}{ Dimension along which to compute seasonal means = Time dimension } \item{monini}{ -First month of the time-series: 1 to 12. +First month of the time series: 1 to 12. } \item{moninf}{ Month when to start the seasonal means: 1 to 12. @@ -27,7 +27,7 @@ Month when to stop the seasonal means: 1 to 12. } } \value{ -Matrix with the same dimensions as var except along the posdim dimension which length corresponds to the number of seasons. Partial seasons are not accounted for. +Array with the same dimensions as var except along the posdim dimension whose length corresponds to the number of seasons. Partial seasons are not accounted for. } \examples{ # Load sample data as in Load() example: diff --git a/man/SelIndices.Rd b/man/SelIndices.Rd index 11c5fdf0b880323720e5c460b1ff7f3f37487126..b6b1795f988172037ca423a5e0e97a9870464bfa 100644 --- a/man/SelIndices.Rd +++ b/man/SelIndices.Rd @@ -4,24 +4,24 @@ Slices A Matrix Along A Dimension } \description{ -This function allows to select a subensemble from a matrix of any dimensions, providing the dimension along which the user aims at cutting the input matrix and between which indices. +This function selects a subset of ensemble members from an array containing any number of dimensions. } \usage{ SelIndices(var, posdim, limits) } \arguments{ \item{var}{ -A matrix of any dimensions. +An array with any number of dimensions. } \item{posdim}{ -The dimension along which a submatrix should be selected. +The dimension along which the ensemble subset should be selected. } \item{limits}{ -The lower and upper indice of the selection along the posdim dimension. +The lower and upper limits for the selection of ensemble members along the posdim dimension. } } \value{ -The sliced matrix. +The subsetted array. } \examples{ a <- array(rnorm(24), dim = c(2, 3, 4, 1)) diff --git a/man/Smoothing.Rd b/man/Smoothing.Rd index 77ecee4001e96c43e38d8e9254d7fa01613dac11..841fb0526f786d56fd96e87d9f9c8d796e0e7b3a 100644 --- a/man/Smoothing.Rd +++ b/man/Smoothing.Rd @@ -1,10 +1,10 @@ \name{Smoothing} \alias{Smoothing} \title{ -Smoothes A Matrix Along A Dimension +Smoothes An Array Along A Dimension } \description{ -Smoothes a matrix of any number of dimensions along one of its dimensions +Smoothes an array of any number of dimensions along one of its dimensions } \usage{ Smoothing(var, runmeanlen = 12, numdimt = 4) diff --git a/man/Spectrum.Rd b/man/Spectrum.Rd index f99567ad362ce9a43c78b42d3eaa3d5c87bce958..eda82779bb128bc1777ba7e9444a91fdf3ed0e16 100644 --- a/man/Spectrum.Rd +++ b/man/Spectrum.Rd @@ -2,7 +2,7 @@ \alias{Spectrum} \title{Estimates Frequency Spectrum} \description{ -This function estimates the frequency spectrum of the xdata array together with its 95\% and 99\% significance level. The output is provided as a matrix with dimensions c(number of frequencies, 4). The column contains the frequency values, the power, the 95\% significance level and the 99\% one.\cr +This function estimates the frequency spectrum of the xdata array together with its 95\% and 99\% significance level. The output is provided as an array with dimensions c(number of frequencies, 4). The column contains the frequency values, the power, the 95\% significance level and the 99\% one.\cr The spectrum estimation relies on a R built-in function and the significance levels are estimated by a Monte-Carlo method. } \usage{Spectrum(xdata)} diff --git a/man/StatSeasAtlHurr.Rd b/man/StatSeasAtlHurr.Rd index 5ee2fc2e71c15e8ab8708ee1888788228ce4b0ab..3acb0edb02e82bc817a42eda236909e10299df97 100644 --- a/man/StatSeasAtlHurr.Rd +++ b/man/StatSeasAtlHurr.Rd @@ -5,11 +5,6 @@ Compute one of G. Villarini's statistically downscaled measure of mean Atlantic hurricane activity and its variance. The hurricane activity is estimated using seasonal averages of sea surface temperature anomalies over the tropical Atlantic (bounded by 10N-25N and 80W-20W) and the tropics at large (bounded by 30N-30S). The anomalies are for the JJASON season.\cr The estimated seasonal average is either 1) number of hurricanes, 2) number of tropical cyclones with lifetime >=48h or 3) power dissipation index (PDI; in 10^11 m^3 s^{-2}).\cr The statistical models used in this function are described in\cr -Villarini et al. (2010) Mon Wea Rev, 138, 2681-2705.\cr -Villarini et al. (2012) Mon Wea Rev, 140, 44-65.\cr -Villarini et al. (2012) J Clim, 25, 625-637.\cr -An example of how the function can be used in hurricane forecast studies is given in\cr -Caron, L.-P. et al. (2014) Multi-year prediction skill of Atlantic hurricane activity in CMIP5 decadal hindcasts. Climate Dynamics, 42, 2675-2690. doi:10.1007/s00382-013-1773-1. } \usage{ StatSeasAtlHurr(atlano = NULL, tropano = NULL, hrvar = 'HR') @@ -62,6 +57,13 @@ TropAno <- matrix(c(-0.22, -.13, 0.07, -0.16, -0.15, hr_count <- StatSeasAtlHurr(atlano = AtlAno, tropano = TropAno, hrvar = 'HR')$mean print(hr_count) } +\references{ +Villarini et al. (2010) Mon Wea Rev, 138, 2681-2705.\cr +Villarini et al. (2012) Mon Wea Rev, 140, 44-65.\cr +Villarini et al. (2012) J Clim, 25, 625-637.\cr +An example of how the function can be used in hurricane forecast studies is given in\cr +Caron, L.-P. et al. (2014) Multi-year prediction skill of Atlantic hurricane activity in CMIP5 decadal hindcasts. Climate Dynamics, 42, 2675-2690. doi:10.1007/s00382-013-1773-1. +} \author{ History:\cr 0.1 - 2015-11 (Louis-Philippe Caron, \email{louis-philippe.caron@bsc.es}) - Original code diff --git a/man/ToyModel.Rd b/man/ToyModel.Rd index 54c3cf79b0dbb45481c2e2f9c2ecbe0d7b3c7990..b009ef0fdb6d4c772d01775c912d32b41b26c988 100644 --- a/man/ToyModel.Rd +++ b/man/ToyModel.Rd @@ -1,7 +1,7 @@ \name{ToyModel} \alias{ToyModel} \title{ -Synthetic forecast generator imitating seasonal to decadal forecasts. Tg +Synthetic forecast generator imitating seasonal to decadal forecasts. The components of a forecast: (1) predictabiltiy (2) forecast error (3) non-stationarity and (4) ensemble generation. The forecast can be computed for real observations or observations generated artifically. diff --git a/man/Trend.Rd b/man/Trend.Rd index 959e32a99f9f269db9f80dc203d756a3a940d190..716e05439d435289b2d9fbef241894eec5485dfb 100644 --- a/man/Trend.Rd +++ b/man/Trend.Rd @@ -1,23 +1,28 @@ \name{Trend} \alias{Trend} +\alias{.Trend} \title{ -Computes Trends +Computes the Trend of the Ensemble Mean } \description{ -Computes the trend along the posTR dimension of the matrix var by least square fitting, and the associated an error interval.\cr -Provide also the detrended data.\cr +Computes the trend along the forecast time of the ensemble mean by least square fitting, and the associated error interval.\cr +Trend() also provides the time series of the detrended ensemble mean forecasts.\cr The confidence interval relies on a student-T distribution. + +.Trend provides the same functionality but taking a matrix ensemble members as input (exp). } \usage{ Trend(var, posTR = 2, interval = 1, siglev = 0.95, conf = TRUE) + +.Trend(exp, interval = 1, siglev = 0.95, conf = TRUE) } \arguments{ - \item{var}{ -Matrix of any number of dimensions up to 10. - } - \item{posTR}{ -Position along which to compute the trend. + \item{var}{ +Array of any number of dimensions up to 10. } + \item{exp}{ +M by N matrix of M forecasts from N ensemble members. + } \item{interval}{ Number of months/years between 2 points along posTR dimension.\cr Default = 1.\cr @@ -29,11 +34,17 @@ Confidence level for the computation of confidence interval. 0.95 by default. \item{conf}{ Whether to compute the confidence levels or not. TRUE by default. } + \item{posTR}{ +Position along which to compute the trend. + } } \value{ \item{$trend}{ -Same dimensions as var except along the posTR dimension which is replaced by a dimension of length 2 or 4, corresponding to the lower limit of the \code{siglev}\% confidence interval (only present if \code{conf = TRUE}, trends, the upper limit of the \code{siglev}\% confidence interval (only present if \code{conf = TRUE}), and intercept of the trend for each point of the matrix along all the other dimensions. +The intercept and slope coefficients for the least squares fitting of the trend, } + \item{$conf.int}{ +Corresponding to the limits of the \code{siglev}\% confidence interval (only present if \code{conf = TRUE}) for the slope coefficient. + } \item{$detrended}{ Same dimensions as var with linearly detrended var along the posTR dimension. } @@ -54,6 +65,7 @@ PlotAno(trend$detrended, NULL, startDates, \author{ History:\cr 0.1 - 2011-05 (V. Guemas, \email{virginie.guemas at ic3.cat}) - Original code\cr -1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at ic3.cat}) - Formatting to CRAN +1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at ic3.cat}) - Formatting to CRAN\cr +2.0 - 2017-02 (A. Hunter, \email{alasdair.hunter at bsc.es}) - Adapt to veriApply() } \keyword{datagen}