diff --git a/R/BrierScore.R b/R/BrierScore.R index c64434ed5ea97810f72b03b90ef18535943462a4..cce753770921b116b73b5d2cdb56af39d8d8c9a6 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,98 @@ 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(ens, obs, thresholds = seq(0, 1, 0.1)) { + if (max(ens) > 1 | min(ens) < 0) { + stop("Predictions outside [0,1] range. Are you certain this is a probability forecast? \n") + } else if (max(obs) != 1 & min(obs) != 0) { + .message("Binary events must be either 0 or 1. Are you certain this is a binary event? ") + } else { + nbins <- length(thresholds) - 1 # Number of bins + n <- dim(ens)[1] # Number of observations + ens.mean <- rowMeans(ens, na.rm = TRUE) + n.ens <- seq(1,dim(ens)[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(ens, 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..ce576b84b805a99302a765dba6730719e94ae38a 100644 --- a/R/Corr.R +++ b/R/Corr.R @@ -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,57 @@ Corr <- function(var_exp, var_obs, posloop = 1, poscor = 2, compROW = NULL, # CORR } + +.Corr <- function(ens, obs, siglev = 0.95, method = 'pearson', + conf = TRUE, pval = TRUE) { + + if (method != "kendall" && method != "spearman" && method != "pearson") { + stop("Wrong correlation method") + + + # Check the siglev arguments: + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + if (siglev > 1 || siglev < 0) { + stop("siglev need to be higher than O and lower than 1") + } + } + + + p <- c() + conflow <- c() + confhigh <- c() + + ens.mean <- rowMeans(ens) + 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 <- 1 - pt(t, eno-2) + p.val <- p + + } + 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] + } + + invisible(result <- list(corr = CORR, p.val = p, conf_low = conflow, conf_high = confhigh)) + # Output + # ~~~~~~~~ + # + +} diff --git a/R/RMS.R b/R/RMS.R index 1a63b48fdbf4eac911db3c8f47d2ade99bf8a5a9..04bf52b24e72df8ffee853be1b27cd8d6d491454 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,35 @@ RMS <- function(var_exp, var_obs, posloop = 1, posRMS = 2, compROW = NULL, # enlrms } + +.RMS <- function(ens, obs, limits = NULL, siglev = 0.95, conf = TRUE) { + + + # + # RMS & its confidence interval computation + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # + if (conf) { + + conf_low <- (1 - siglev) / 2 + conf_high <- 1 - conf_low + } + + dif <- rowMeans(ens) - 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..74a2a3877c699b5e2377c7b43f2ff3001213f789 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 diff --git a/R/RatioRMS.R b/R/RatioRMS.R index a97ec02f8fb59be249cb5b8cef54f9a8a4aa9cba..a5bc5b4889602bda5cd3cb579bdb95c6ca979936 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,37 @@ RatioRMS <- function(var_exp1, var_exp2, var_obs, posRMS = 1, pval = TRUE) { # enlratiorms } + +.RatioRMS <- function(ens, ens.ref, obs, pval = TRUE) { + + # + # RMS ratio and its pvalue computation + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # + + + dif1 <- rowMeans(ens, na.rm = TRUE) - obs + dif2 <- rowMeans(ens.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..543d6e05950ad7bbbaa0e66767656db2e6f787de 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,46 @@ RatioSDRMS <- function(var_exp, var_obs, pval = TRUE) { # enlratiormssd } + +.RatioSDRMS <- function(ens, obs, pval = TRUE) { + + # + # Ratio RMSE / SD and its significance level + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # + ensmean <- rowMeans(ens) + + std <- sd(ensmean) + enosd <- Eno(ensmean,1) + + dif <- ensmean - obs + rms <- mean(dif ** 2, 1, narm = 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..b8ef4d0cffaa855210c4fc7d6c97ac00210ce7d6 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,26 @@ Trend <- function(var, posTR = 2, interval = 1, siglev = 0.95, conf = TRUE) { # invisible(list(trend = enltrend, detrended = enldetrend)) } + +.Trend <- function(ens, interval = 1, siglev = 0.95, conf = TRUE) { + + ensmean <- rowMeans(ens, 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..2c92fe3ef7222490eb415ea50c53c47ae1eef845 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 (Centred) } \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 jexp in 1:nexp and each jobs in 1:nobs which gives nexp x nobs ACsC 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. @@ -57,6 +56,7 @@ 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,6 +92,9 @@ 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 diff --git a/man/BrierScore.Rd b/man/BrierScore.Rd index 43843523c2d45a29c2f9e479cd02aef0a1446509..d31998c70cc40d5a157c34189fe4bb480c7af689 100644 --- a/man/BrierScore.Rd +++ b/man/BrierScore.Rd @@ -1,5 +1,6 @@ \name{BrierScore} \alias{BrierScore} +\alias{.BrierScore} \title{ Compute Brier Score And Its Decomposition And Brier Skill Score } @@ -8,9 +9,13 @@ Returns the values of the BS and its standard decomposition as well as the addit 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 + +.BrierScore provides the same functionality, but taking a matrix of ensemble members (ens) as input. } \usage{ BrierScore(obs, pred, thresholds = seq(0, 1, 0.1)) + +.BrierScore(ens, obs, thresholds = seq(0, 1, 0.1)) } \arguments{ \item{obs}{ @@ -18,6 +23,9 @@ Vector of binary observations (1 or 0) } \item{pred}{ Vector of probablistic predictions with values in the range [0,1] + } + \item{ens}{ +Matrix of predictions with values in the range [0,1] for the .BrierScore function } \item{thresholds}{ Values used to bin the forecasts. By default the bins are {[0,0.1), [0.1, 0.2), ... [0.9, 1]} diff --git a/man/Corr.Rd b/man/Corr.Rd index 1c5641e0568e4fddcae239172f8ae9222a303ff8..3b2cae755db55ac5f317751d0a1694f0533774ed 100644 --- a/man/Corr.Rd +++ b/man/Corr.Rd @@ -1,83 +1,103 @@ \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 (ens) 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(ens, obs, siglev = 0.95, + method = 'pearson', conf = TRUE, pval = TRUE) } \arguments{ - \item{var_exp}{ -Matrix of experimental data. + \item{var_exp}{ +Array of experimental data. + } + + \item{var_obs}{ +Array of observational data, same dimensions as var_exp except along posloop dimension, where the length can be nobs instead of nexp. + } + \item{ens}{ +N by M matrix of N forecasts from M ensemble members. } - \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. + \item{obs}{ +Vector of the corresponding observations of length N. } - \item{posloop}{ -Dimension nobs and nexp. - } - \item{poscor}{ -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. - } - \item{limits}{ -Complete between limits[1] & limits[2]. Default = NULL. + \item{posloop}{ +Dimension nobs and nexp. + } + \item{poscor}{ +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. + } + \item{limits}{ +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). } } -\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. -} -\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), - fileout = 'tos_cor.eps') -} + + \value{ 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 + +.Corr: + \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), + fileout = 'tos_cor.eps') } + \author{ History:\cr 0.1 - 2011-04 (V. Guemas, \email{vguemas at ic3.cat}) - Original code\cr diff --git a/man/RMS.Rd b/man/RMS.Rd index 118400f407167bcd5d3d8b88953669bc882475bc..663694d53c6334850ed15b9cece7184fba0557d6 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 + +.RMS provides the same functionality but taking a matrix of ensemble members as input (ens). } \usage{ RMS(var_exp, var_obs, posloop = 1, posRMS = 2, compROW = NULL, limits = NULL, siglev = 0.95, conf = TRUE) + +.RMS(ens, obs, limits = NULL, siglev = 0.95, conf = TRUE) } \arguments{ \item{var_exp}{ @@ -20,6 +24,12 @@ Matrix 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. + } + \item{ens}{ +N by M matrix of N forecasts from M ensemble members. + } + \item{obs}{ +Vector of the corresponding observations of length N. } \item{posloop}{ Dimension nobs and nexp. @@ -42,9 +52,19 @@ Whether to compute confidence interval or not. TRUE by default. } } \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 +.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: diff --git a/man/RMSSS.Rd b/man/RMSSS.Rd index 3d77230dbf7156bc3f96646126c5b7b42c637cf9..7e6188364c80c5c20ca7271c6b0ec1ac1783fb3a 100644 --- a/man/RMSSS.Rd +++ b/man/RMSSS.Rd @@ -4,13 +4,14 @@ 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. } \usage{ RMSSS(var_exp, var_obs, posloop = 1, posRMS = 2, pval = TRUE) + } \arguments{ \item{var_exp}{ diff --git a/man/RatioRMS.Rd b/man/RatioRMS.Rd index ee116561f7662d6dfa0ffc15bc3ddf5ceff06b22..b8a7486dbd2e48d35e6b99cfb38a2bb99365a133 100644 --- a/man/RatioRMS.Rd +++ b/man/RatioRMS.Rd @@ -1,27 +1,41 @@ \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 +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. + +.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(ens, ens.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{ens}{ +Matrix of experimental data 1. } - \item{posRMS}{ + \item{ens.ref}{ +Matrix of experimental data 2. + } + \item{obs}{ +Vector of observations. + } + \item{posRMS}{ Dimension along which the RMSE are to be computed = the position of the start dates. } \item{pval}{ @@ -29,8 +43,18 @@ Whether to compute the p-value of Ho : RMSE1/RMSE2 = 1 or not. TRUE by default. } } \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: A list containing the following components :\cr +ratiorms - the ratio of the rms of the two ensembles \cr +p.val - the p value \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.\cr +.RatioRMS: +\item{$ratiorms}{ +The ratio of the root mean square errors, + } + \item{$conf.int}{ +Corresponding to the p values of the ratio of the rmse statistics (only present if \code{pval = TRUE}). + } + } \examples{ # See examples on Load() to understand the first lines in this example @@ -77,6 +101,7 @@ PlotEquiMap(rrms[1, , ], sampleData$lon, sampleData$lat, \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 +1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens at ic3.cat}) - Formatting to R CRAN\cr +2.0 - 2016-08 (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..e38bd2c757b8cf1fb76a40e85bd333a0204dd72c 100644 --- a/man/RatioSDRMS.Rd +++ b/man/RatioSDRMS.Rd @@ -1,21 +1,25 @@ \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 +.RatioSDRMS provides the same functionality but taking a matrix of ensemble members as input (ens). } \usage{ RatioSDRMS(var_exp, var_obs, pval = TRUE) + +.RatioSDRMS(ens, 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) @@ -24,15 +28,29 @@ Model data:\cr Observational data:\cr c(nobs, nmemb, nsdates, nltime) up to\cr c(nobs, nmemb, nsdates, nltime, nlevel, nlat, nlon) + } + \item{ens}{ +N by M matrix of N forecasts from M ensemble members. + } + \item{obs}{ +Vector of the corresponding observations of length N. } \item{pval}{ Whether to compute the p-value of Ho : SD/RMSE = 1 or not. } } \value{ -Matrix with dimensions c(nexp/nmod, nobs, 1 or 2, nltime) up to +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). -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. +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 +.RatioSDRMS: +\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: diff --git a/man/Trend.Rd b/man/Trend.Rd index 959e32a99f9f269db9f80dc203d756a3a940d190..c2981e0f3c8c8cb20cae4d06974ab1b35704cb83 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 (ens). } \usage{ Trend(var, posTR = 2, interval = 1, siglev = 0.95, conf = TRUE) + +.Trend(ens, 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{ens}{ +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 - 2016-08 (A. Hunter, \email{alasdair.hunter at bsc.es}) - Adapt to veriApply() } \keyword{datagen}